text:
function [varargout] = ft_selectdata(cfg, varargin)

% FT_SELECTDATA makes a selection in the input data along specific data
% dimensions, such as channels, time, frequency, trials, etc. It can also
% be used to average the data along each of the specific dimensions.
%
% Use as
%  [data] = ft_selectdata(cfg, data, ...)
%
% The cfg argument is a configuration structure which can contain
%   cfg.tolerance   = scalar, tolerance value to determine equality of time/frequency bins (default = 1e-5)
%
% For data with trials or subjects as repetitions, you can specify
%   cfg.trials      = 1xN, trial indices to keep, can be 'all'. You can use logical indexing, where false(1,N) removes all the trials
%   cfg.avgoverrpt  = string, can be 'yes' or 'no' (default = 'no')
%
% For data with a channel dimension you can specify
%   cfg.channel     = Nx1 cell-array with selection of channels (default = 'all'), see FT_CHANNELSELECTION
%   cfg.avgoverchan = string, can be 'yes' or 'no' (default = 'no')
%   cfg.nanmean     = string, can be 'yes' or 'no' (default = 'no')
%
% For data with channel combinations you can specify
%   cfg.channelcmb     = Nx2 cell-array with selection of channels (default = 'all'), see FT_CHANNELCOMBINATION
%   cfg.avgoverchancmb = string, can be 'yes' or 'no' (default = 'no')
%
% For data with a time dimension you can specify
%   cfg.latency     = scalar or string, can be 'all', 'minperiod', 'maxperiod', 'prestim', 'poststim', or [beg end], specify time range in seconds
%   cfg.avgovertime = string, can be 'yes' or 'no' (default = 'no')
%   cfg.nanmean     = string, can be 'yes' or 'no' (default = 'no')
%
% For data with a frequency dimension you can specify
%   cfg.frequency   = scalar or string, can be 'all', or [beg end], specify frequency range in Hz
%   cfg.avgoverfreq = string, can be 'yes' or 'no' (default = 'no')
%   cfg.nanmean     = string, can be 'yes' or 'no' (default = 'no')
%
% When you average over a dimension, you can choose whether to keep that dimension in
% the data representation or remove it alltogether.
%   cfg.keeprptdim     = 'yes' or 'no' (default is automatic)
%   cfg.keepchandim    = 'yes' or 'no' (default = 'yes')
%   cfg.keepchancmbdim = 'yes' or 'no' (default = 'yes')
%   cfg.keeptimedim    = 'yes' or 'no' (default = 'yes')
%   cfg.keepfreqdim    = 'yes' or 'no' (default = 'yes')
%
% If multiple input arguments are provided, FT_SELECTDATA will adjust the individual
% inputs such that either the INTERSECTION across inputs is retained (i.e. only the
% channel, time, and frequency points that are shared across all input arguments), or
% that the UNION across inputs is retained (replacing missing data with nans). In
% either case, the order of the channels is made consistent across inputs. The
% behavior can be specified with
%   cfg.select      = string, can be 'intersect' or 'union' (default = 'intersect')
% For raw data structures you cannot make the union.
%
% See also FT_DATATYPE, FT_CHANNELSELECTION, FT_CHANNELCOMBINATION

% Undocumented options
%   cfg.avgoverpos
%   cfg.keepposdim     = 'yes' or 'no' (default = 'yes')

% Copyright (C) 2012-2022, Robert Oostenveld & Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

ft_defaults
ft_preamble init
ft_preamble debug

ft_preamble loadvar varargin
ft_preamble provenance varargin

% determine the characteristics of the input data
dtype = ft_datatype(varargin{1});
for i=2:length(varargin)
  % ensure that all subsequent inputs are of the same type
  ok = ft_datatype(varargin{i}, dtype);
  if ~ok, ft_error('input data should be of the same datatype'); end
end

% this only works with certain data types, it is not meant for descriptive fields such as elec, grad, opto, layout, etc.
assert(~ismember(dtype, {'elec', 'grad', 'opto', 'layout'}), 'invalid input data type "%s"', dtype);

% ensure that the user does not give invalid selection options
cfg = ft_checkconfig(cfg, 'forbidden', {'foi', 'toi'});

cfg = ft_checkconfig(cfg, 'renamed',    {'selmode',    'select'});
cfg = ft_checkconfig(cfg, 'renamed',    {'toilim',     'latency'});
cfg = ft_checkconfig(cfg, 'renamed',    {'foilim',     'frequency'});
cfg = ft_checkconfig(cfg, 'renamed',    {'avgoverroi', 'avgoverpos'});
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'avg.pow', 'pow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'avg.mom', 'mom'});
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'avg.nai', 'nai'});
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'trial.pow', 'pow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'trial.mom', 'mom'});
cfg = ft_checkconfig(cfg, 'renamedval', {'parameter', 'trial.nai', 'nai'});

cfg.tolerance = ft_getopt(cfg, 'tolerance', 1e-5);        % default tolerance for checking equality of time/freq axes
cfg.select    = ft_getopt(cfg, 'select',   'intersect');  % default is to take intersection, alternative 'union'
if isequal(dtype, 'raw') && isequal(cfg.select, 'union')
  ft_error('using cfg.select=''union'' in combination with ''raw'' datatype is not supported');
end

if strcmp(dtype, 'volume') || strcmp(dtype, 'segmentation')
  % it must be a source representation, not a volume representation
  for i=1:length(varargin)
    varargin{i} = ft_checkdata(varargin{i}, 'datatype', 'source');
  end
  dtype = 'source';
else
  % check that the data is according to the latest FieldTrip representation
  for i=1:length(varargin)
    varargin{i} = ft_checkdata(varargin{i});
  end
end

% this function only works for the upcoming (not yet standard) source representation without sub-structures
% update the old-style beamformer source reconstruction to the upcoming representation
if strcmp(dtype, 'source')
  if isfield(varargin{1}, 'avg') && isstruct(varargin{1}.avg)
    restoreavg = fieldnames(varargin{1}.avg);
  else
    restoreavg = {};
  end
  for i=1:length(varargin)
    varargin{i} = ft_datatype_source(varargin{i}, 'version', 'upcoming');
  end
end

cfg.latency = ft_getopt(cfg, 'latency', 'all', 1);
if isnumeric(cfg.latency) && numel(cfg.latency)==2 && cfg.latency(1)==cfg.latency(2)
  % this is better specified by a single number
  cfg.latency = cfg.latency(1);
end

cfg.channel = ft_getopt(cfg, 'channel', 'all', 1);
cfg.trials  = ft_getopt(cfg, 'trials',  'all', 1);

if length(varargin)>1 && ~isequal(cfg.trials, 'all')
  ft_error('it is ambiguous to make a subselection of trials while at the same time concatenating multiple data structures')
end

cfg.frequency = ft_getopt(cfg, 'frequency', 'all', 1);
if isnumeric(cfg.frequency) && numel(cfg.frequency)==2 && cfg.frequency(1)==cfg.frequency(2)
  % this is better specified by a single number
  cfg.frequency = cfg.frequency(1);
end

datfield  = fieldnames(varargin{1});
for i=2:length(varargin)
  % only consider fields that are present in all inputs
  datfield = intersect(datfield, fieldnames(varargin{i}));
end
datfield  = setdiff(datfield, {'label' 'labelcmb'});  % these fields will be used for selection, but are not treated as data fields
datfield  = setdiff(datfield, {'dim'});               % not used for selection, also not treated as data field
datfield  = setdiff(datfield, ignorefields('selectdata'));
orgdim1   = datfield(~cellfun(@isempty, regexp(datfield, 'label$')) & cellfun(@isempty, regexp(datfield, '^csd'))); % xxxlabel, with the exception of csdlabel
datfield  = setdiff(datfield, orgdim1);
datfield  = datfield(:)';
orgdim1   = datfield(~cellfun(@isempty, regexp(datfield, 'dimord$'))); % xxxdimord
datfield  = setdiff(datfield, orgdim1);
datfield  = datfield(:)';

sel = strcmp(datfield, 'cumtapcnt');
if any(sel)
  % move this field to the end, as it is needed to make the selections in the other fields
  datfield(sel) = [];
  datfield = [datfield {'cumtapcnt'}];
end

orgdim2 = cell(size(orgdim1));
for i=1:length(orgdim1)
  orgdim2{i} = varargin{1}.(orgdim1{i});
end

dimord = cell(size(datfield));
for i=1:length(datfield)
  dimord{i} = getdimord(varargin{1}, datfield{i});
end

% do not consider fields of which the dimensions are unknown
% sel = cellfun(@isempty, regexp(dimord, 'unknown'));
% for i=find(~sel)
%   fprintf('not including "%s" in selection\n', datfield{i});
% end
% datfield = datfield(sel);
% dimord   = dimord(sel);

% determine all dimensions that are present in all data fields
dimtok = {};
for i=1:length(datfield)
  dimtok = cat(2, dimtok, tokenize(dimord{i}, '_'));
end
dimtok = unique(dimtok);

hasspike   = any(ismember(dimtok, 'spike'));
haspos     = any(ismember(dimtok, {'pos' '{pos}'}));
haschan    = any(ismember(dimtok, {'chan' '{chan}'}));
haschancmb = any(ismember(dimtok, 'chancmb'));
hasfreq    = any(ismember(dimtok, 'freq'));
hastime    = any(ismember(dimtok, 'time'));
hasrpt     = any(ismember(dimtok, {'rpt' 'subj' '{rpt}'}));
hasrpttap  = any(ismember(dimtok, 'rpttap'));

if hasspike
  % cfg.latency is used to select individual spikes, not to select from a continuously sampled time axis
  hastime = false;
end

clear dimtok

haspos     = haspos     && isfield(varargin{1}, 'pos');
haschan    = haschan    && isfield(varargin{1}, 'label');
haschancmb = haschancmb && isfield(varargin{1}, 'labelcmb');
hasfreq    = hasfreq    && isfield(varargin{1}, 'freq');
hastime    = hastime    && isfield(varargin{1}, 'time');

% do a sanity check on all input arguments
if haspos,     assert(all(cellfun(@isfield, varargin, repmat({'pos'},      size(varargin)))), 'not all input arguments have a "pos" field'); end
if haschan,    assert(all(cellfun(@isfield, varargin, repmat({'label'},    size(varargin)))), 'not all input arguments have a "label" field'); end
if haschancmb, assert(all(cellfun(@isfield, varargin, repmat({'labelcmb'}, size(varargin)))), 'not all input arguments have a "labelcmb" field'); end
if hasfreq,    assert(all(cellfun(@isfield, varargin, repmat({'freq'},     size(varargin)))), 'not all input arguments have a "freq" field'); end
if hastime,    assert(all(cellfun(@isfield, varargin, repmat({'time'},     size(varargin)))), 'not all input arguments have a "time" field'); end

avgoverpos     = istrue(ft_getopt(cfg, 'avgoverpos',  false)); % at some places it is also referred to as roi (region-of-interest)
avgoverchan    = istrue(ft_getopt(cfg, 'avgoverchan', false));
avgoverchancmb = istrue(ft_getopt(cfg, 'avgoverchancmb', false));
avgoverfreq    = istrue(ft_getopt(cfg, 'avgoverfreq', false));
avgovertime    = istrue(ft_getopt(cfg, 'avgovertime', false));
avgoverrpt     = istrue(ft_getopt(cfg, 'avgoverrpt',  false));

% do a sanity check for the averaging options
if avgoverpos,     assert(haspos,     'there are no source positions, so averaging is not possible'); end
if avgoverchan,    assert(haschan,    'there is no channel dimension, so averaging is not possible'); end
if avgoverchancmb, assert(haschancmb, 'there are no channel combinations, so averaging is not possible'); end
if avgoverfreq,    assert(hasfreq,    'there is no frequency dimension, so averaging is not possible'); end
if avgovertime,    assert(hastime,    'there is no time dimension, so averaging over time is not possible'); end
if avgoverrpt,     assert(hasrpt||hasrpttap, 'there are no repetitions, so averaging is not possible'); end

% set averaging function
cfg.nanmean = ft_getopt(cfg, 'nanmean', 'no');
if strcmp(cfg.nanmean, 'yes')
  average = @nanmean;
else
  average = @mean;
end

% by default we keep most of the dimensions in the data structure when averaging over them
keepposdim     = istrue(ft_getopt(cfg, 'keepposdim',  true));
keepchandim    = istrue(ft_getopt(cfg, 'keepchandim', true));
keepchancmbdim = istrue(ft_getopt(cfg, 'keepchancmbdim', true));
keepfreqdim    = istrue(ft_getopt(cfg, 'keepfreqdim', true));
keeptimedim    = istrue(ft_getopt(cfg, 'keeptimedim', true));
keeprptdim     = istrue(ft_getopt(cfg, 'keeprptdim', ~avgoverrpt));

if ~keepposdim,     assert(avgoverpos,     'removing a dimension is only possible when averaging'); end
if ~keepchandim,    assert(avgoverchan,    'removing a dimension is only possible when averaging'); end
if ~keepchancmbdim, assert(avgoverchancmb, 'removing a dimension is only possible when averaging'); end
if ~keepfreqdim,    assert(avgoverfreq,    'removing a dimension is only possible when averaging'); end
if ~keeptimedim,    assert(avgovertime,    'removing a dimension is only possible when averaging'); end
if ~keeprptdim,     assert(avgoverrpt,     'removing a dimension is only possible when averaging'); end

% trim the selection to all inputs, rpt and rpttap are dealt with later
if hasspike,   [selspike,   cfg] = getselection_spike  (cfg, varargin{:}); end
if haspos,     [selpos,     cfg] = getselection_pos    (cfg, varargin{:}, cfg.tolerance, cfg.select); end
if haschan,    [selchan,    cfg] = getselection_chan   (cfg, varargin{:}, cfg.select); end
if haschancmb, [selchancmb, cfg] = getselection_chancmb(cfg, varargin{:}, cfg.select); end
if hasfreq,    [selfreq,    cfg] = getselection_freq   (cfg, varargin{:}, cfg.tolerance, cfg.select); end
if hastime,    [seltime,    cfg] = getselection_time   (cfg, varargin{:}, cfg.tolerance, cfg.select); end

% this is to keep track of all fields that should be retained in the output
keepfield = datfield;

for i=1:numel(varargin)

  for j=1:numel(datfield)
    dimtok = tokenize(dimord{j}, '_');

    % the rpt selection should only work with a single data argument
    % in case tapers were kept, selrpt~=selrpttap, otherwise selrpt==selrpttap
    [selrpt{i}, dum, rptdim{i}, selrpttap{i}] = getselection_rpt(cfg, varargin{i}, dimord{j});

    % check for the presence of each dimension in each datafield
    fieldhasspike   = ismember('spike',   dimtok);
    fieldhaspos     = ismember('pos',     dimtok) || ismember('{pos}', dimtok);
    fieldhaschan    = (ismember('chan',    dimtok) || ismember('{chan}', dimtok)) && isfield(varargin{1}, 'label');
    fieldhaschancmb = ismember('chancmb', dimtok);
    fieldhastime    = ismember('time',    dimtok) && ~hasspike;
    fieldhasfreq    = ismember('freq',    dimtok);
    fieldhasrpt     = ismember('rpt',     dimtok) | ismember('subj', dimtok) | ismember('{rpt}', dimtok);
    fieldhasrpttap  = ismember('rpttap',  dimtok);

    % cfg.latency is used to select individual spikes, not to select from a continuously sampled time axis

    if fieldhasspike,   varargin{i} = makeselection(varargin{i}, datfield{j}, dimtok, find(strcmp(dimtok,'spike')),             selspike{i},   false,           'intersect', average); end
    if fieldhaspos,     varargin{i} = makeselection(varargin{i}, datfield{j}, dimtok, find(ismember(dimtok, {'pos', '{pos}'})), selpos{i},     avgoverpos,      cfg.select, average);  end
    if fieldhaschan,    varargin{i} = makeselection(varargin{i}, datfield{j}, dimtok, find(ismember(dimtok,{'chan' '{chan}'})), selchan{i},    avgoverchan,     cfg.select, average);  end
    if fieldhaschancmb, varargin{i} = makeselection(varargin{i}, datfield{j}, dimtok, find(strcmp(dimtok,'chancmb')),           selchancmb{i}, avgoverchancmb,  cfg.select, average);  end
    if fieldhastime,    varargin{i} = makeselection(varargin{i}, datfield{j}, dimtok, find(strcmp(dimtok,'time')),              seltime{i},    avgovertime,     cfg.select, average);  end
    if fieldhasfreq,    varargin{i} = makeselection(varargin{i}, datfield{j}, dimtok, find(strcmp(dimtok,'freq')),              selfreq{i},    avgoverfreq,     cfg.select, average);  end
    if fieldhasrpt,     varargin{i} = makeselection(varargin{i}, datfield{j}, dimtok, rptdim{i},                                selrpt{i},     avgoverrpt,      'intersect', average); end
    if fieldhasrpttap,  varargin{i} = makeselection(varargin{i}, datfield{j}, dimtok, rptdim{i},                                selrpttap{i},  avgoverrpt,      'intersect', average); end

    % update the fields that should be kept in the structure as a whole
    % and update the dimord for this specific datfield
    keepdim = true(size(dimtok));

    if avgoverchan && ~keepchandim
      keepdim(strcmp(dimtok, 'chan')) = false;
      keepfield = setdiff(keepfield, 'label');
    else
      keepfield = [keepfield 'label'];
    end

    if avgoverchancmb && ~keepchancmbdim
      keepdim(strcmp(dimtok, 'chancmb')) = false;
      keepfield = setdiff(keepfield, 'labelcmb');
    else
      keepfield = [keepfield 'labelcmb'];
    end

    if avgoverfreq && ~keepfreqdim
      keepdim(strcmp(dimtok, 'freq')) = false;
      keepfield = setdiff(keepfield, 'freq');
    else
      keepfield = [keepfield 'freq'];
    end

    if avgovertime && ~keeptimedim
      keepdim(strcmp(dimtok, 'time')) = false;
      keepfield = setdiff(keepfield, 'time');
    else
      keepfield = [keepfield 'time'];
    end

    if avgoverpos && ~keepposdim
      keepdim(strcmp(dimtok, 'pos'))   = false;
      keepdim(strcmp(dimtok, '{pos}')) = false;
      keepdim(strcmp(dimtok, 'dim'))   = false;
      keepfield = setdiff(keepfield, {'pos' '{pos}' 'dim'});
    elseif avgoverpos && keepposdim
      keepfield = setdiff(keepfield, {'dim'}); % this should be removed anyway
    else
      keepfield = [keepfield {'pos' '{pos}' 'dim'}];
    end

    if avgoverrpt && ~keeprptdim
      keepdim(strcmp(dimtok, 'rpt'))    = false;
      keepdim(strcmp(dimtok, 'rpttap')) = false;
      keepdim(strcmp(dimtok, 'subj'))   = false;
    end

    % update the sampleinfo, if possible, and needed
    if strcmp(datfield{j}, 'sampleinfo') && ~isequal(cfg.latency, 'all')
      if iscell(seltime{i}) && numel(seltime{i})==size(varargin{i}.sampleinfo,1)
        for k = 1:numel(seltime{i})
          varargin{i}.sampleinfo(k,:) = varargin{i}.sampleinfo(k,1) - 1 + seltime{i}{k}([1 end]);
        end
      elseif ~iscell(seltime{i}) && ~isempty(seltime{i}) && ~all(isnan(seltime{i}))
        nrpt       = size(varargin{i}.sampleinfo,1);
        seltime{i} = seltime{i}(:)';
        varargin{i}.sampleinfo = varargin{i}.sampleinfo(:,[1 1]) - 1 + repmat(seltime{i}([1 end]),nrpt,1);
      end
    end

    varargin{i}.(datfield{j}) = squeezedim(varargin{i}.(datfield{j}), ~keepdim);

  end % for datfield

  % also update the fields that describe the dimensions, time/freq/pos have been dealt with as data
  if haschan,    varargin{i} = makeselection_chan   (varargin{i}, selchan{i}, avgoverchan); end % update the label field
  if haschancmb, varargin{i} = makeselection_chancmb(varargin{i}, selchancmb{i}, avgoverchancmb); end % update the labelcmb field

end % for varargin

if strcmp(cfg.select, 'union')
  % create the union of the descriptive axes
  if haspos,      varargin = makeunion(varargin, 'pos'); end
  if haschan,     varargin = makeunion(varargin, 'label'); end
  if haschancmb,  varargin = makeunion(varargin, 'labelcmb'); end
  if hastime,     varargin = makeunion(varargin, 'time'); end
  if hasfreq,     varargin = makeunion(varargin, 'freq'); end
end

% remove all fields from the data structure that do not pertain to the selection
sel = strcmp(keepfield, '{pos}'); if any(sel), keepfield(sel) = {'pos'}; end
sel = strcmp(keepfield, 'chan');  if any(sel), keepfield(sel) = {'label'}; end
sel = strcmp(keepfield, 'chancmb');  if any(sel), keepfield(sel) = {'labelcmb'}; end

if avgoverrpt
  % these are invalid after averaging
  keepfield = setdiff(keepfield, {'cumsumcnt' 'cumtapcnt' 'trialinfo' 'sampleinfo'});
end

if avgovertime
  % these are invalid after averaging or making a latency selection
  keepfield = setdiff(keepfield, {'sampleinfo'});
end

for i=1:numel(varargin)
  varargin{i} = keepfields(varargin{i}, [keepfield ignorefields('selectdata')']);
end

% restore the original dimord fields in the data
for i=1:length(orgdim1)
  dimtok = tokenize(orgdim2{i}, '_');

  % using a setdiff may result in double occurrences of chan and pos to
  % disappear, so this causes problems as per bug 2962
  % if ~keeprptdim, dimtok = setdiff(dimtok, {'rpt' 'rpttap' 'subj'}); end
  % if ~keepposdim, dimtok = setdiff(dimtok, {'pos' '{pos}'}); end
  % if ~keepchandim, dimtok = setdiff(dimtok, {'chan'}); end
  % if ~keepfreqdim, dimtok = setdiff(dimtok, {'freq'}); end
  % if ~keeptimedim, dimtok = setdiff(dimtok, {'time'}); end
  if ~keeprptdim, dimtok = dimtok(~ismember(dimtok, {'rpt' 'rpttap' 'subj'})); end
  if ~keepposdim, dimtok = dimtok(~ismember(dimtok, {'pos' '{pos}'}));         end
  if ~keepchandim, dimtok = dimtok(~ismember(dimtok, {'chan'})); end
  if ~keepfreqdim, dimtok = dimtok(~ismember(dimtok, {'freq'})); end
  if ~keeptimedim, dimtok = dimtok(~ismember(dimtok, {'time'})); end

  dimord = sprintf('%s_', dimtok{:});
  dimord = dimord(1:end-1); % remove the trailing _
  for j=1:length(varargin)
    varargin{j}.(orgdim1{i}) = dimord;
  end
end

% restore the source.avg field, this keeps the output reasonably consistent with the
% old-style source representation of the input
if strcmp(dtype, 'source') && ~isempty(restoreavg)
  for i=1:length(varargin)
    varargin{i}.avg = keepfields(varargin{i}, restoreavg);
    varargin{i}     = removefields(varargin{i}, restoreavg);
  end
end

varargout = varargin;

ft_postamble debug

ft_postamble previous varargin
ft_postamble provenance varargout
ft_postamble history varargout
ft_postamble savevar varargout

% the varargout variable can be cleared when written to outputfile
if exist('varargout', 'var') && ft_nargout>numel(varargout)
  % also return the input cfg with the combined selection over all input data structures
  varargout{end+1} = cfg;
end

end % function ft_selectdata

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = makeselection(data, datfield, dimtok, seldim, selindx, avgoverdim, selmode, average)

if numel(seldim) > 1
  for k = 1:numel(seldim)
    data = makeselection(data, datfield, dimtok, seldim(k), selindx, avgoverdim, selmode, average);
  end
  return;
end

if isnumeric(data.(datfield))
  if isrow(data.(datfield)) && seldim==1
    if length(dimtok)==1
      seldim = 2; % switch row and column
    end
  elseif iscolumn(data.(datfield)) && seldim==2
    if length(dimtok)==1
      seldim = 1; % switch row and column
    end
  end
elseif iscell(data.(datfield))
  if isrow(data.(datfield){1}) && seldim==2
    if length(dimtok)==2
      seldim = 3; % switch row and column
    end
  elseif iscolumn(data.(datfield){1}) && seldim==3
    if length(dimtok)==2
      seldim = 2; % switch row and column
    end
  end
end

% an empty selindx means that nothing(!) should be selected and hence everything should be removed, which is different than keeping everything
% the selindx value of NaN indicates that it is not needed to make a selection

switch selmode
  case 'intersect'
    if iscell(selindx)
      % there are multiple selections in multipe vectors, the selection is in the matrices contained within the cell-array
      for j=1:numel(selindx)
        if ~isempty(selindx{j}) && all(isnan(selindx{j}))
          % no selection needs to be made
        else
          data.(datfield){j} = cellmatselect(data.(datfield){j}, seldim-1, selindx{j}, numel(dimtok)==1);
        end
      end

    else
      % there is a single selection in a single vector
      if ~isempty(selindx) && all(isnan(selindx))
        % no selection needs to be made
      else
        data.(datfield) = cellmatselect(data.(datfield), seldim, selindx, numel(dimtok)==1);
      end
    end

    if avgoverdim
      data.(datfield) = cellmatmean(data.(datfield), seldim, average);
    end

  case 'union'
    if ~isempty(selindx) && all(isnan(selindx))
      % no selection needs to be made
    elseif isequal(seldim,1) && any(strcmp({'time' 'freq'}, datfield))
      % treat this as an exception, because these fields should only be
      % unionized along the second dimension, so here also no selection
      % needs to be made
    else
      tmp = data.(datfield);
      siz = size(tmp);
      siz(seldim) = numel(selindx);
      data.(datfield) = nan(siz);
      sel = isfinite(selindx);
      switch seldim
        case 1
          data.(datfield)(sel,:,:,:,:,:) = tmp(selindx(sel),:,:,:,:,:);
        case 2
          data.(datfield)(:,sel,:,:,:,:) = tmp(:,selindx(sel),:,:,:,:);
        case 3
          data.(datfield)(:,:,sel,:,:,:) = tmp(:,:,selindx(sel),:,:,:);
        case 4
          data.(datfield)(:,:,:,sel,:,:) = tmp(:,:,:,selindx(sel),:,:);
        case 5
          data.(datfield)(:,:,:,:,sel,:) = tmp(:,:,:,:,selindx(sel),:);
        case 6
          data.(datfield)(:,:,:,:,:,sel) = tmp(:,:,:,:,:,selindx(sel));
        otherwise
          ft_error('unsupported dimension (%d) for making a selection for %s', seldim, datfield);
      end
    end

    if avgoverdim
      data.(datfield) = average(data.(datfield), seldim);
    end
end % switch

end % function makeselection

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = makeselection_chan(data, selchan, avgoverchan)
if isempty(selchan)
  %error('no channels were selected');
  data.label = {};
elseif avgoverchan && all(isnan(selchan))
  str = sprintf('%s, ', data.label{:});
  str = str(1:end-2);
  str = sprintf('mean(%s)', str);
  data.label = {str};
elseif avgoverchan && ~any(isnan(selchan))
  str = sprintf('%s, ', data.label{selchan});
  str = str(1:end-2);
  str = sprintf('mean(%s)', str);
  data.label = {str};                 % remove the last '+'
elseif all(isfinite(selchan))
  data.label = data.label(selchan);
  data.label = data.label(:);
elseif numel(selchan)==1 && any(~isfinite(selchan))
  % do nothing
elseif numel(selchan)>1  && any(~isfinite(selchan))
  tmp = cell(numel(selchan),1);
  for k = 1:numel(tmp)
    if isfinite(selchan(k))
      tmp{k} = data.label{selchan(k)};
    end
  end
  data.label = tmp;
else
  % this should never happen
  ft_error('cannot figure out how to select channels');
end
end % function makeselection_chan

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data = makeselection_chancmb(data, selchancmb, avgoverchancmb)
if isempty(selchancmb)
  ft_error('no channel combinations were selected');
elseif avgoverchancmb && all(isnan(selchancmb))
  % naming the channel combinations becomes ambiguous, but should not
  % suggest that the mean was computed prior to combining
  str1 = sprintf('%s, ', data.labelcmb{:,1});
  str1 = str1(1:end-2);
  % str1 = sprintf('mean(%s)', str1);
  str2 = sprintf('%s, ', data.labelcmb{:,2});
  str2 = str2(1:end-2);
  % str2 = sprintf('mean(%s)', str2);
  data.label = {str1, str2};
elseif avgoverchancmb && ~any(isnan(selchancmb))
  % naming the channel combinations becomes ambiguous, but should not
  % suggest that the mean was computed prior to combining
  str1 = sprintf('%s, ', data.labelcmb{selchancmb,1});
  str1 = str1(1:end-2);
  % str1 = sprintf('mean(%s)', str1);
  str2 = sprintf('%s, ', data.labelcmb{selchancmb,2});
  str2 = str2(1:end-2);
  % str2 = sprintf('mean(%s)', str2);
  data.label = {str1, str2};
elseif all(isfinite(selchancmb))
  data.labelcmb = data.labelcmb(selchancmb,:);
elseif numel(selchancmb)==1 && any(~isfinite(selchancmb))
  % do nothing
elseif numel(selchancmb)>1  && any(~isfinite(selchancmb))
  tmp = cell(numel(selchancmb),2);
  for k = 1:size(tmp,1)
    if isfinite(selchan(k))
      tmp{k,1} = data.labelcmb{selchan(k),1};
      tmp{k,2} = data.labelcmb{selchan(k),2};
    end
  end
  data.labelcmb = tmp;
else
  % this should never happen
  ft_error('cannot figure out how to select channelcombinations');
end
end % function makeselection_chancmb

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [chanindx, cfg] = getselection_chan(cfg, varargin)

selmode  = varargin{end};
ndata    = numel(varargin)-1;
varargin = varargin(1:ndata);

% loop over data once to initialize
chanindx = cell(ndata,1);
label    = cell(1,0);

if ndata==1 && (isequal(cfg.channel, 'all') || isequal(cfg.channel, varargin{1}.label))
  % the loop across data arguments, as well as the expensive calls to
  % FT_CHANNELSELECTION can be avoided if there is only a single data
  % argument and if 'all' channels are to be returned in the output
  label = varargin{1}.label(:);

else
  for k = 1:ndata
    selchannel      = cell(0,1);
    selgrad         = [];
    selelec         = [];
    selopto         = [];
    if isfield(varargin{k}, 'grad') && isfield(varargin{k}.grad, 'type')
      % this makes channel selection more robust, e.g. when using wildcards in cfg.channel
      [selgrad, dum] = match_str(varargin{k}.label, varargin{k}.grad.label);
      selchannel     = cat(1, selchannel, ft_channelselection(cfg.channel, varargin{k}.label(selgrad), varargin{k}.grad.type));
    end
    if isfield(varargin{k}, 'elec') && isfield(varargin{k}.elec, 'type')
      % this makes channel selection more robust, e.g. when using wildcards in cfg.channel
      [selelec, dum] = match_str(varargin{k}.label, varargin{k}.elec.label);
      selchannel     = cat(1, selchannel, ft_channelselection(cfg.channel, varargin{k}.label(selelec), varargin{k}.elec.type));
    end
    if isfield(varargin{k}, 'opto') && isfield(varargin{k}.opto, 'type')
      % this makes channel selection more robust, e.g. when using wildcards in cfg.channel
      [selopto, dum] = match_str(varargin{k}.label, varargin{k}.opto.label);
      selchannel     = cat(1, selchannel, ft_channelselection(cfg.channel, varargin{k}.label(selopto), varargin{k}.opto.type));
    end
    selrest    = setdiff((1:numel(varargin{k}.label))', [selgrad; selelec; selopto]);
    selchannel = cat(1, selchannel, ft_channelselection(cfg.channel, varargin{k}.label(selrest)));
    label      = union(label, selchannel);
  end
  label = label(:);   % ensure that this is a column array

  % this call to match_str ensures that that labels are always in the
  % order of the first input argument see bug_2917, but also temporarily keep
  % the labels from the other data structures not present in the first one
  % (in case selmode = 'union')
  [ix, iy] = match_str(varargin{1}.label, label);
  label1   = varargin{1}.label(:); % ensure column array
  label    = [label1(ix); label(setdiff(1:numel(label),iy))];

end % if ndata==1 and all channels are to be returned

indx = nan+zeros(numel(label), ndata);
for k = 1:ndata
  [ix, iy] = match_str(label, varargin{k}.label);
  indx(ix,k) = iy;
end

switch selmode
  case 'intersect'
    sel      = sum(isfinite(indx),2)==ndata;
    indx     = indx(sel,:);
    label    = varargin{1}.label(indx(:,1));
  case 'union'
    % don't do a subselection
  otherwise
    ft_error('invalid value for cfg.select');
end % switch

ok = false(size(indx,1),1);
for k = 1:ndata
  % loop through the columns to preserve the order of the channels, where
  % the order of the input arguments determines the final order
  ix = find(~ok);
  [srt,srtix] = sort(indx(ix,k));
  indx(ix,:)  = indx(ix(srtix),:);
  ok = ok | isfinite(indx(:,k));
end

for k = 1:ndata
  % do a sanity check on double occurrences
  if numel(unique(indx(isfinite(indx(:,k)),k)))<sum(isfinite(indx(:,k)))
    ft_error('the selection of channels across input arguments leads to double occurrences');
  end
  chanindx{k} = indx(:,k);
end

for k = 1:ndata
  if isequal(chanindx{k}, (1:numel(varargin{k}.label))')
    % no actual selection is needed for this data structure
    chanindx{k} = nan;
  end
end

cfg.channel = label;

end % function getselection_chan

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [chancmbindx, cfg] = getselection_chancmb(cfg, varargin)

selmode  = varargin{end};
ndata    = numel(varargin)-1;
varargin = varargin(1:ndata);

chancmbindx = cell(ndata,1);

if ~isfield(cfg, 'channelcmb')
  for k=1:ndata
    % the nan return value specifies that no selection was specified
    chancmbindx{k} = nan;
  end

else

  switch selmode
    case 'intersect'
      haslabel = false(ndata,1);
      for k=1:ndata
        haslabel = isfield(varargin{k}, 'label');
      end
      if all(haslabel)
        for k=1:ndata
          cfg.channelcmb = ft_channelcombination(cfg.channelcmb, varargin{k}.label);
        end
        cfgcmb = cellfun(@sprintf,repmat({'%s_%s'},size(cfg.channelcmb,1),1),cfg.channelcmb(:,1),cfg.channelcmb(:,2),'UniformOutput',false);
      elseif all(~haslabel)
        % the data already has labelcmb, and thus needs a slightly different way to
        % preset the cfg.channelcmb
        chancmb = cellfun(@sprintf,repmat({'%s_%s'},size(varargin{1}.labelcmb,1),1),varargin{1}.labelcmb(:,1),varargin{1}.labelcmb(:,2),'UniformOutput',false);
        for k=2:ndata
          tmp = cellfun(@sprintf,repmat({'%s_%s'},size(varargin{k}.labelcmb,1),1),varargin{k}.labelcmb(:,1),varargin{k}.labelcmb(:,2),'UniformOutput',false);
          chancmb = intersect(chancmb, tmp);
        end
        cfgcmb = unique(chancmb);

        if isequal(cfg.channelcmb, {'all' 'all'})
          % nothing needed here
        else
          cfg.channelcmb = cellfun(@sprintf,repmat({'%s_%s'},size(cfg.channelcmb,1),1),cfg.channelcmb(:,1),cfg.channelcmb(:,2),'UniformOutput',false);

          cfgcmb = intersect(cfg.channelcmb, cfgcmb);
        end

      else
        ft_error('a combination of data with and without label field is not possible');
      end

      for k=1:ndata
        datcmb = cellfun(@sprintf,repmat({'%s_%s'},size(varargin{k}.labelcmb,1),1),varargin{k}.labelcmb(:,1),varargin{k}.labelcmb(:,2),'UniformOutput',false);

        % return the order according to the (joint) configuration, not according to the (individual) data
        % FIXME this should adhere to the general code guidelines, where
        % the order returned will be according to the first data argument!
        [dum, chancmbindx{k}] = match_str(cfgcmb, datcmb);
      end

    case 'union'
      % FIXME this is not yet implemented
      ft_error('union of channel combination is not yet supported');

    otherwise
      ft_error('invalid value for cfg.select');
  end % switch

end

end % function getselection_chancmb

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [spikeindx, cfg] = getselection_spike(cfg, varargin)
% possible specifications are
% cfg.latency = string -> 'all'
% cfg.latency = [beg end]
% cfg.trials  = string -> 'all'
% cfg.trials  = vector with indices

ndata    = numel(varargin);
varargin = varargin(1:ndata);

if isequal(cfg.latency, 'all') && isequal(cfg.trials, 'all')
  spikeindx = cell(1,ndata);
  for i=1:ndata
    spikeindx{i} = num2cell(nan(1, length(varargin{i}.time)));
  end
  return
end

trialbeg = varargin{1}.trialtime(:,1);
trialend = varargin{1}.trialtime(:,2);
for i=2:ndata
  trialbeg = cat(1, trialbeg, varargin{1}.trialtime(:,1));
  trialend = cat(1, trialend, varargin{1}.trialtime(:,2));
end

% convert string into a numeric selection
if ischar(cfg.latency)
  switch cfg.latency
    case 'all'
      cfg.latency = [-inf inf];
    case 'maxperiod'
      cfg.latency = [min(trialbeg) max(trialend)];
    case 'minperiod'
      cfg.latency = [max(trialbeg) min(trialend)];
    case 'prestim'
      cfg.latency = [min(trialbeg) 0];
    case 'poststim'
      cfg.latency = [0 max(trialend)];
    otherwise
      ft_error('incorrect specification of cfg.latency');
  end % switch
end

spikeindx = cell(1,ndata);
for i=1:ndata
  nchan = length(varargin{i}.time);
  spikeindx{i} = cell(1,nchan);
  for j=1:nchan
    selbegtime = varargin{i}.time{j}>=cfg.latency(1);
    selendtime = varargin{i}.time{j}<=cfg.latency(2);
    if isequal(cfg.trials, 'all')
      seltrial = true(size(varargin{i}.trial{j}));
    else
      seltrial = ismember(varargin{i}.trial{j}, cfg.trials);
    end
    spikeindx{i}{j} = find(selbegtime & selendtime & seltrial);
  end
end

end % function getselection_spiketime

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [timeindx, cfg] = getselection_time(cfg, varargin)
% possible specifications are
% cfg.latency = value     -> can be 'all'
% cfg.latency = [beg end]

if ft_datatype(varargin{1}, 'spike')
  ft_error('latency selection in spike data is not supported')
end

selmode  = varargin{end};
tol      = varargin{end-1};
ndata    = numel(varargin)-2;
varargin = varargin(1:ndata);

if isequal(cfg.latency, 'all') && iscell(varargin{1}.time)
  % for raw data this means that all trials should be selected as they are
  % for timelock/freq data it is still needed to make the intersection between data arguments
  timeindx = cell(1,ndata);
  for i=1:ndata
    % the nan return value specifies that no selection was specified
    timeindx{i} = num2cell(nan(1, length(varargin{i}.time)));
  end
  return
end

% if there is a single timelock/freq input, there is one time vector
% if there are multiple timelock/freq inputs, there are multiple time vectors
% if there is a single raw input, there are multiple time vectors
% if there are multiple raw inputs, there are multiple time vectors

% collect all time axes in one large cell-array
alltimecell = {};
if iscell(varargin{1}.time)
  for k = 1:ndata
    alltimecell = [alltimecell varargin{k}.time{:}];
  end
else
  for k = 1:ndata
    alltimecell = [alltimecell {varargin{k}.time}];
  end
end

% the nan return value specifies that no selection was specified
timeindx = repmat({nan}, size(alltimecell));

% loop over data once to determine the union of all time axes
alltimevec = zeros(1,0);
for k = 1:length(alltimecell)
  alltimevec = union(alltimevec, round(alltimecell{k}/tol)*tol);
end

indx = nan(numel(alltimevec), numel(alltimecell));
for k = 1:numel(alltimecell)
  [dum, ix, iy] = intersect(alltimevec, round(alltimecell{k}/tol)*tol);
  indx(ix,k) = iy;
end

if iscell(varargin{1}.time) && ~isequal(cfg.latency, 'minperiod')
  % if the input data arguments are of type 'raw', temporarily set the
  % selmode to union, otherwise the potentially different length trials
  % will be truncated to the shortest epoch, prior to latency selection.
  selmode = 'union';
elseif ischar(cfg.latency) && strcmp(cfg.latency, 'minperiod')
  % enforce intersect
  selmode = 'intersect';
end
switch selmode
  case 'intersect'
    sel        = sum(isfinite(indx),2)==numel(alltimecell);
    indx       = indx(sel,:);
    alltimevec = alltimevec(sel);
  case 'union'
    % don't do a subselection
  otherwise
    ft_error('invalid value for cfg.select');
end

% Note that cfg.toilim handling has been removed, as it was renamed to cfg.latency

% convert a string selection into a numeric selection
if ischar(cfg.latency)
  switch cfg.latency
    case {'all' 'maxperlen' 'maxperiod'}
      cfg.latency = [min(alltimevec) max(alltimevec)];
    case 'prestim'
      cfg.latency = [min(alltimevec) 0];
    case 'poststim'
      cfg.latency = [0 max(alltimevec)];
    case 'minperiod'
      % the time vector has been pruned above
      cfg.latency = [min(alltimevec) max(alltimevec)];
    otherwise
      ft_error('incorrect specification of cfg.latency');
  end % switch
end

% deal with numeric selection
if isempty(cfg.latency)
  for k = 1:numel(alltimecell)
    % FIXME I do not understand this
    % this signifies that all time bins are deselected and should be removed
    timeindx{k} = [];
  end

elseif numel(cfg.latency)==1
  % this single value should be within the time axis of each input data structure
  if numel(alltimevec)>1
    tbin = nearest(alltimevec, cfg.latency, true, true); % determine the numerical tolerance
  else
    tbin = nearest(alltimevec, cfg.latency, true, false); % don't consider tolerance
  end
  cfg.latency = alltimevec(tbin);

  for k = 1:ndata
    timeindx{k} = indx(tbin, k);
  end

elseif numel(cfg.latency)==2
  % the [min max] range can be specifed with +inf or -inf, but should
  % at least partially overlap with the time axis of the input data
  mintime = min(alltimevec);
  maxtime = max(alltimevec);
  if all(cfg.latency<mintime) || all(cfg.latency>maxtime)
    ft_error('the selected time range falls outside the time axis in the data');
  end
  tbeg = nearest(alltimevec, cfg.latency(1), false, false);
  tend = nearest(alltimevec, cfg.latency(2), false, false);
  cfg.latency = alltimevec([tbeg tend]);

  for k = 1:numel(alltimecell)
    timeindx{k} = indx(tbeg:tend, k);
    % if the input data arguments are of type 'raw', the non-finite values
    % need to be removed from the individual cells to ensure correct
    % behavior
    if iscell(varargin{1}.time)
      timeindx{k} = timeindx{k}(isfinite(timeindx{k}));
    end
  end

elseif size(cfg.latency,2)==2
  % this may be used for specification of the computation, not for data selection

else
  ft_error('incorrect specification of cfg.latency');
end

for k = 1:numel(alltimecell)
  if ~iscell(varargin{1}.time)
    if isequal(timeindx{k}(:)', 1:length(alltimecell{k}))
      % no actual selection is needed for this data structure
      timeindx{k} = nan;
    end
  else
    % if the input data arguments are of type 'raw', they need to be
    % handled differently, because the individual trials can be of
    % different length
  end
end

if iscell(varargin{1}.time)
  % split all time axes again over the different input raw data structures
  dum = cell(1,ndata);
  for k = 1:ndata
    sel = 1:length(varargin{k}.time);
    dum{k} = timeindx(sel); % get the first selection
    timeindx(sel) = []; % remove the first selection
  end
  timeindx = dum;
else
  % no splitting is needed, each input data structure has one selection
end

end % function getselection_time

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [freqindx, cfg] = getselection_freq(cfg, varargin)
% possible specifications are
% cfg.frequency = value     -> can be 'all'
% cfg.frequency = [beg end]

selmode  = varargin{end};
tol      = varargin{end-1};
ndata    = numel(varargin)-2;
varargin = varargin(1:ndata);

% loop over data once to initialize
freqindx = cell(ndata,1);
freqaxis = zeros(1,0);
for k = 1:ndata
  % the nan return value specifies that no selection was specified
  freqindx{k} = nan;

  % update the axis along which the frequencies are defined
  freqaxis = union(freqaxis, round(varargin{k}.freq(:)/tol)*tol);
end

indx = nan+zeros(numel(freqaxis), ndata);
for k = 1:ndata
  [dum, ix, iy] = intersect(freqaxis, round(varargin{k}.freq(:)/tol)*tol);
  indx(ix,k) = iy;
end

switch selmode
  case 'intersect'
    sel      = sum(isfinite(indx),2)==ndata;
    indx     = indx(sel,:);
    freqaxis = varargin{1}.freq(indx(:,1));
  case 'union'
    % don't do a subselection
  otherwise
    ft_error('invalid value for cfg.select');
end

if isfield(cfg, 'frequency')
  % deal with string selection
  % some of these do not make sense, but are here for consistency with ft_multiplotER
  if ischar(cfg.frequency)
    if strcmp(cfg.frequency, 'all')
      cfg.frequency = [min(freqaxis) max(freqaxis)];
    elseif strcmp(cfg.frequency, 'maxmin')
      cfg.frequency = [min(freqaxis) max(freqaxis)]; % the same as 'all'
    elseif strcmp(cfg.frequency, 'minzero')
      cfg.frequency = [min(freqaxis) 0];
    elseif strcmp(cfg.frequency, 'maxabs')
      cfg.frequency = [-max(abs(freqaxis)) max(abs(freqaxis))];
    elseif strcmp(cfg.frequency, 'zeromax')
      cfg.frequency = [0 max(freqaxis)];
    elseif strcmp(cfg.frequency, 'zeromax')
      cfg.frequency = [0 max(freqaxis)];
    else
      ft_error('incorrect specification of cfg.frequency');
    end
  end

  % deal with numeric selection
  if isempty(cfg.frequency)
    for k = 1:ndata
      % FIXME I do not understand this
      % this signifies that all frequency bins are deselected and should be removed
      freqindx{k} = [];
    end

  elseif numel(cfg.frequency)==1
    % this single value should be within the frequency axis of each input data structure
    if numel(freqaxis)>1
      fbin = nearest(freqaxis, cfg.frequency, true, true); % determine the numerical tolerance
    else
      fbin = nearest(freqaxis, cfg.frequency, true, false); % don't consider tolerance
    end
    cfg.frequency = freqaxis(fbin);

    for k = 1:ndata
      freqindx{k} = indx(fbin,k);
    end

  elseif numel(cfg.frequency)==2
    % the [min max] range can be specifed with +inf or -inf, but should
    % at least partially overlap with the freq axis of the input data
    minfreq = min(freqaxis);
    maxfreq = max(freqaxis);
    if all(cfg.frequency<minfreq) || all(cfg.frequency>maxfreq)
      ft_error('the selected range falls outside the frequency axis in the data');
    end
    fbeg = nearest(freqaxis, cfg.frequency(1), false, false);
    fend = nearest(freqaxis, cfg.frequency(2), false, false);
    cfg.frequency = freqaxis([fbeg fend]);

    for k = 1:ndata
      freqindx{k} = indx(fbeg:fend,k);
    end

  elseif size(cfg.frequency,2)==2
    % this may be used for specification of the computation, not for data selection

  else
    ft_error('incorrect specification of cfg.frequency');
  end
end % if cfg.frequency

for k = 1:ndata
  if isequal(freqindx{k}, 1:length(varargin{k}.freq))
    % the cfg was updated, but no selection is needed for the data
    freqindx{k} = nan;
  end
end

end % function getselection_freq

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rptindx, cfg, rptdim, rpttapindx] = getselection_rpt(cfg, varargin)
% this should deal with cfg.trials

dimord   = varargin{end};
ndata    = numel(varargin)-1;
data     = varargin{1:ndata}; % this syntax ensures that it will only work on a single data input

dimtok = tokenize(dimord, '_');
rptdim = find(strcmp(dimtok, '{rpt}') | strcmp(dimtok, 'rpt') | strcmp(dimtok, 'rpttap') | strcmp(dimtok, 'subj'));

if isequal(cfg.trials, 'all')
  rptindx    = nan; % the nan return value specifies that no selection was specified
  rpttapindx = nan; % the nan return value specifies that no selection was specified

elseif isempty(rptdim)
  % FIXME should [] not mean that none of the trials is to be selected?
  rptindx    = nan; % the nan return value specifies that no selection was specified
  rpttapindx = nan; % the nan return value specifies that no selection was specified

else
  rptindx = ft_getopt(cfg, 'trials');

  if islogical(rptindx)
    % convert from booleans to indices
    rptindx = find(rptindx);
  end

  rptindx = unique(sort(rptindx));

  if strcmp(dimtok{rptdim}, 'rpttap') && isfield(data, 'cumtapcnt')
    % there are tapers in the data

    % determine for each taper to which trial it belongs
    nrpt = size(data.cumtapcnt, 1);
    taper = zeros(nrpt, 1);
    sumtapcnt = cumsum([0; data.cumtapcnt(:)]);
    begtapcnt = sumtapcnt(1:end-1)+1;
    endtapcnt = sumtapcnt(2:end);
    for i=1:nrpt
      taper(begtapcnt(i):endtapcnt(i)) = i;
    end
    rpttapindx = find(ismember(taper, rptindx));

  else
    % there are no tapers in the data
    rpttapindx = rptindx;
  end
end

end % function getselection_rpt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [posindx, cfg] = getselection_pos(cfg, varargin)
% possible specifications are <none>

ndata   = numel(varargin)-2;
tol     = varargin{end-1}; % FIXME this is still ignored
selmode = varargin{end};   % FIXME this is still ignored
data    = varargin(1:ndata);

for i=1:ndata
  if ~isequal(varargin{i}.pos, varargin{1}.pos)
    % FIXME it would be possible here to make a selection based on intersect or union
    ft_error('not yet implemented');
  end
end

if strcmp(cfg.select, 'union')
  % FIXME it would be possible here to make a selection based on intersect or union
  ft_error('not yet implemented');
end

for i=1:ndata
  posindx{i} = nan;    % the nan return value specifies that no selection was specified
end
end % function getselection_pos

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = squeezedim(x, dim)
siz = size(x);
for i=(numel(siz)+1):numel(dim)
  % all trailing singleton dimensions have length 1
  siz(i) = 1;
end
if isvector(x) && ~(isrow(x) && dim(1) && numel(x)>1)
  % there is no harm to keep it as it is, unless the data matrix is 1xNx1x1
elseif istable(x)
  % there is no harm to keep it as it is
else
  x = reshape(x, [siz(~dim) 1]);
end
end % function squeezedim

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = makeunion(x, field)
old = cellfun(@getfield, x, repmat({field}, size(x)), 'uniformoutput', false);
if iscell(old{1})
  % empty is indicated to represent missing value for a cell-array (label, labelcmb)
  new = old{1};
  for i=2:length(old)
    sel = ~cellfun(@isempty, old{i});
    new(sel) = old{i}(sel);
  end
else
  % nan is indicated to represent missing value for a numeric array (time, freq, pos)
  new = old{1};
  for i=2:length(old)
    sel = ~isnan(old{i});
    new(sel) = old{i}(sel);
  end
end
x = cellfun(@setfield, x, repmat({field}, size(x)), repmat({new}, size(x)), 'uniformoutput', false);
end % function makeunion

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION to make a selextion in data representations like {pos}_ori_time
% FIXME this will fail for {xxx_yyy}_zzz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = cellmatselect(x, seldim, selindx, maybevector)
if nargin<4
  % some fields are a vector with an unspecified singleton dimension, these can be transposed
  % if the singleton dimension represents something explicit, they should not be transposed
  % they might for example represent a single trial, or a single channel
  maybevector = true;
end
if iscell(x)
  if seldim==1
    x = x(selindx);
  else
    for i=1:numel(x)
      if isempty(x{i})
        continue
      end
      switch seldim
        case 2
          if maybevector && isvector(x{i})
            % sometimes the data is 1xN, whereas the dimord describes only the first dimension
            % in this case a row and column vector can be interpreted as equivalent
            x{i} = x{i}(selindx);
          elseif istable(x)
            % multidimensional indexing is not supported
            x{i} = x{i}(selindx,:);
          else
            x{i} = x{i}(selindx,:,:,:,:);
          end
        case 3
          x{i} = x{i}(:,selindx,:,:,:);
        case 4
          x{i} = x{i}(:,:,selindx,:,:);
        case 5
          x{i} = x{i}(:,:,:,selindx,:);
        case 6
          x{i} = x{i}(:,:,:,:,selindx);
        otherwise
          ft_error('unsupported dimension (%d) for making a selection', seldim);
      end % switch
    end % for
  end
else
  switch seldim
    case 1
      if maybevector && isvector(x)
        % sometimes the data is 1xN, whereas the dimord describes only the first dimension
        % in this case a row and column vector can be interpreted as equivalent
        x = x(selindx);
      elseif istable(x)
        % multidimensional indexing is not supported
        x = x(selindx,:);
      else
        x = x(selindx,:,:,:,:,:);
      end
    case 2
      x = x(:,selindx,:,:,:,:);
    case 3
      x = x(:,:,selindx,:,:,:);
    case 4
      x = x(:,:,:,selindx,:,:);
    case 5
      x = x(:,:,:,:,selindx,:);
    case 6
      x = x(:,:,:,:,:,selindx);
    otherwise
      ft_error('unsupported dimension (%d) for making a selection', seldim);
  end
end
end % function cellmatselect

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION to take an average in data representations like {pos}_ori_time
% FIXME this will fail for {xxx_yyy}_zzz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = cellmatmean(x, seldim, average)
if iscell(x)
  if seldim==1
    for i=2:numel(x)
      x{1} = x{1} + x{i};
    end
    x = {x{1}/numel(x)};
  else
    for i=1:numel(x)
      x{i} = average(x{i}, seldim-1);
    end % for
  end
elseif istable(x)
  try
    % try to convert to an array, depending on the table content this might fail
    x = average(table2array(x), seldim);
  catch
    % construct an appropriately sized array with NaN values
    s = size(x);
    s(seldim) = 1;
    x = nan(s);
  end
  % convert back to table
  x = array2table(x);
else
  x = average(x, seldim);
end
end % function cellmatmean

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%DEMO_BAYESIANOPTIMIZATION1  A demonstration program for Bayesian
%                            optimization in 1 dimension
%
% The set of BO demos
%  Part 1:  this file
%  One dimensional example
%
%  Part 2:  see demo_bayesoptimization2
%  Two dimensional example
%
%  Part 3:  see demo_bayesoptimization3
%  Two dimensional example with constraints
%  * The implementation of constraints follows Gelbart et al. (2014)
%
%  References:
%    Jones, D., Schonlau, M., & Welch, W. (1998). Efficient global
%    optimization of expensive black-box functions. Journal of Global
%    Optimization, 13(4), 455-492. doi:10.1023/a:1008306431147
%
%    Michael A. Gelbart, Jasper Snoek, and Ryan P. Adams
%    (2014). Bayesian Optimization with Unknown Constraints.
%    http://arxiv.org/pdf/1403.5607v1.pdf
%
%    Snoek, J., Larochelle, H, Adams, R. P. (2012). Practical Bayesian
%    Optimization of Machine Learning Algorithms. NIPS 25
%
%  Copyright (c) 2015-2017 Jarno Vanhatalo

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

%%  Part 1:
%  One dimensional example
% For testing purposes:
stack = dbstack;
if (~isempty(stack) && (strcmp(stack(end).name, 'runtestset') || strcmp(stack(end).name, 'runtests'))) test = 1; else test = 0; end;

% Construct a function to be optimized
xl = linspace(0,10,100)';
fx = @(x) 0.6*x -0.1*x.^2 + sin(2*x);

% construct GP
cfse = gpcf_sexp('lengthScale',1,'magnSigma2',1,'magnSigma2_prior',prior_sqrtt('s2',10^2));
lik = lik_gaussian('sigma2', 0.001, 'sigma2_prior', prior_fixed);
gp = gp_set('cf', {cfse}, 'lik', lik);

% ----- conduct Bayesian optimization -----
% draw initial point

% Set the options for optimizer of the acquisition function
optimf = @fmincon;
optdefault=struct('GradObj','on','LargeScale','off','Algorithm','SQP','TolFun',1e-6,'TolX',1e-3);
opt=optimset(optdefault);
lb=0;     % lower bound of the input space
ub=10;    % upper bound of the input space

% draw initial point
rng(3)
x = 10*rand;
y = fx(x);

figure, % figure for visualization
i1 = 1;
maxiter = 15;
improv = inf;   % improvement between two successive query points
while i1 < maxiter && improv>1e-6
%while i1 < maxiter

    % Train the GP model for objective function and calculate variables
    % that are needed when calculating the Expected improvement
    % (Acquisition function)
    if i1>1
        gp = gp_optim(gp,x,y);
    end
    [K, C] = gp_trcov(gp,x);
    invC = inv(C);
    a = C\y;
    fmin = min( fx(x) );

    % Calculate EI and posterior of the function for visualization purposes
    EI = expectedimprovement_eg(xl, gp, x, a, invC, fmin);
    [Ef,Varf] = gp_pred(gp, x, y, xl);

    % optimize acquisition function
    %    Note! Opposite to the standard notation we minimize negative Expected
    %    Improvement since Matlab optimizers seek for functions minimum
    % Here we use multiple starting points for the optimization so that we
    % don't crash into suboptimal mode
    fh_eg = @(x_new) expectedimprovement_eg(x_new, gp, x, a, invC, fmin); % The function handle to the Expected Improvement function
    indbest = find(y == fmin);
    xstart = [linspace(0.5,9.5,5) x(indbest)+0.1*randn(1,2)];
    for s1=1:length(xstart)
        x_new(s1) = optimf(fh_eg, xstart(s1), [], [], [], [], lb, ub, [], opt);
    end
    EIs = expectedimprovement_eg(x_new(:), gp, x, a, invC, fmin);
    x_new = x_new( find(EIs==min(EIs),1) ); % pick up the point where Expected Improvement is maximized

    % put new sample point to the list of evaluation points
    x(end+1) = x_new;
    y(end+1) = fx(x(end));  % calculate the function value at query point
    x=x(:);y=y(:);

    % visualize
    clf
    subplot(2,1,1),hold on, title('function to be optimized and GP fit')
    %plot(xl,fx(xl))
    box on
    plot(xl,fx(xl),'r')
    % The function evaluations so far
    plot(x(1:end-1),y(1:end-1), 'ko')
    % The new sample location
    plot(x(end),y(end), 'ro')
    % the posterior of the function
    plot(xl,Ef, 'k')
    plot(xl,Ef + 2*sqrt(Varf), 'k--')
    plot(xl,Ef - 2*sqrt(Varf), 'k--')
    legend('objective function', 'function evaluations', 'next query point', 'GP mean', 'GP 95% interval','location','southwest')
    % The expected information
    subplot(2,1,2)
    plot(xl,EI, 'r'), hold on
    plot(x(end),0, 'r*')
    plot(x(end)*[1 1],ylim, 'r--')
    title('acquisition function')

    improv = abs(y(end) - y(end-1));
    i1=i1+1;

    if test == 0
        pause
    end
end
%subplot(2,1,1)
%plot(xl,fx(xl),'r')

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function dir_names=get_sub_dir_names(root_dir)

file_infos=dir(root_dir);
dir_names=[];
dir_dates=[];

for f_idx=1:length(file_infos)

    one_f_info=file_infos(f_idx);
    if one_f_info.isdir
        dir_name=one_f_info.name;

        if ~strcmp(dir_name(1), '.')
            dir_names=cat(1, dir_names, {dir_name});
            dir_dates=cat(1, dir_dates, datenum(one_f_info.date));
        end
    end

end

% sort by create date:
[~, tmp_sort_idxes]=sort(dir_dates);
dir_names=dir_names(tmp_sort_idxes);

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
clear all

% Build mex function
mex -R2018a Wigner_d_fast.c

%% Compare different Wigner-d methods
N = 512
beta =rand*pi;

% past method
tic
a = Wigner_D(N,beta);
toc

% new faster method
tic
b = Wigner_d_fast(N,beta);
toc

% create cell array of all Wigner-ds with harmonic degree smaller then L
tic
c = Wigner_d_recursion(N,beta);
toc

% get new Wigner-d from previous two
tic
d = Wigner_d_recursion(c{end-1},c{end-2},beta);
toc

% get new Wigner-d from previous two (faster)
tic
e = Wigner_d_fast(c{end-1},c{end-2},beta);
toc

max(max(abs(a-b)))

max(max(abs(d-e)))

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function value = meixner ( n, beta, c, x )

%*****************************************************************************80
%
%% MEIXNER evaluates Meixner polynomials at a point.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Walter Gautschi,
%    Orthogonal Polynomials: Computation and Approximation,
%    Oxford, 2004,
%    ISBN: 0-19-850672-4,
%    LC: QA404.5 G3555.
%
%  Parameters:
%
%    Input, integer N, the maximum order of the polynomial.
%    N must be at least 0.
%
%    Input, real BETA, the Beta parameter.  0 < BETA.
%
%    Input, real C, the C parameter.  0 < C < 1.
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE(N+1), the value of the polynomials at X.
%
  if ( beta <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MEIXNER - Fatal error!\n' );
    fprintf ( 1, '  Parameter BETA must be positive.\n' );
    error ( 'MEIXNER - Fatal error!\n');
  end

  if ( c <= 0.0 || 1.0 <= c )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MEIXNER - Fatal error!\n' );
    fprintf ( 1, '  Parameter C must be strictly between 0 and 1.\n' );
    error ( 'MEIXNER - Fatal error!\n');
  end

  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MEIXNER - Fatal error!\n' );
    fprintf ( 1, '  Parameter N must be nonnegative.\n' );
    error ( 'MEIXNER - Fatal error!\n');
  end

  OFFSET = 1;

  value(0+OFFSET) = 1.0;

  if ( n == 0 )
    return
  end

  value(1+OFFSET) = ( c - 1 ) * x / beta / c + 1.0;

  if ( n == 1 )
    return
  end

  for i = 1 : n - 1
    value(i+1+OFFSET) = ( ...
      ( ( c - 1.0 ) * x + ( 1.0 + c ) * i + beta * c ) * value(i+OFFSET) ...
      - i * value(i-1+OFFSET) ...
      ) / ( i + beta );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function z = spadaptdemo
% SPADAPTDEMO   Simple example of dimension-adaptive interpolation
%    A 2D-example for dimension-adaptive multi-linear sparse grid
%    interpolation using the Chebyshev grid and vectorized processing
%    of the model function.
%
%    See also SPINTERP, SPVALS.

% Author : Andreas Klimke, Universitaet Stuttgart
% Version: 1.1
% Date   : September 29, 2003

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

% The Branin function is used here
f = inline(['(5/pi*x-5.1/(4*pi^2)*x.^2+y-6).^2 + 10*(1-1/(8*pi))*' ...
	    ' cos(x)+10']);

% Define objective box
range = [-5 10; 0 15];

% Define problem dimension
d = 2;

% Create full grid for plotting
gs = 33;
[X,Y] = meshgrid(linspace(range(1,1),range(1,2),gs),...
								 linspace(range(2,1),range(2,2),gs));

% Set options: Switch vectorized processing on.
options = spset('Vectorized', 'on','RelTol', 1e-2, ...
								'GridType', 'Chebyshev', ...
								'DimensionAdaptive', 'on', ...
								'DimadaptDegree', 1, ...
                'MinPoints', 10);

% Compute sparse grid weights over range
z = spvals(f, d, range, options);

% Compute inpterpolated values at full grid
ip = spinterp(z, X, Y);

% Plot original function, interpolation, and error
subplot(2,2,1);
mesh(X,Y,f(X,Y));
title('original');
axis tight;

subplot(2,2,2);
mesh(X,Y,ip);
title('interpolated');
axis tight;

subplot(2,2,3);
mesh(X,Y,abs(f(X,Y)-ip));
title('absolute error');
axis tight;

subplot(2,2,4);
plotindices(z);
title('resulting index sets');


text:
function [ po, pc, pe ] = lagrange_complete ( d, n, r, nd, xd )

%*****************************************************************************80
%
%% LAGRANGE_COMPLETE: Complete Lagrange polynomial basis from data.
%
%  Discussion:
%
%    This function represents algorithm 4.1 in the reference.
%
%    This function is given XD, a set of ND distinct data points in a
%    D dimensional space, and returns information defining a set of
%    ND Lagrange polynomials L(i)(X) with the property that:
%
%      L(i)(XD(j)) = delta(i,j)
%
%    This function does a very weak form of pivoting, by selecting
%    the very next polynomial that is "nonzero" at the current point.
%    An improved version of this function chooses instead the polynomial
%    with maximum function value at the current point.
%
%    In order for this computation to be carried out, it is necessary that
%    ND, the number of data points, is equal to R, the dimension of the
%    space of polynomials in D dimensions and total degree N or less, that is:
%
%      ND = R = Choose ( N + D, N )
%
%    There will be ND polynomials returned.  Each polynomial can have
%    as many as R coefficients.
%
%    Each polynomial is given as a vector, with each entry corresponding
%    to a nonzero coefficient.  In particular, for polynomial L(i)(X):
%
%      PO(i) is the order, that is, the number of nonzero coefficients;
%      PC(i,j), for 1 <= j <= PO(i), is the coefficient of the J-th term.
%      PE(i,j), for 1 <= j <= PO(i), encodes the exponents of the J-th term.
%
%    The exponent codes are a compact way of recording the exponent vector
%    associated with each monomial.  If PE(i,j) = k, then the corresponding
%    vector of D exponents can be determined by:
%
%      E = mono_unrank_grlex ( D, k );
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Tomas Sauer, Yuan Xu,
%    On multivariate Lagrange interpolation,
%    Mathematics of Computation,
%    Volume 64, Number 211, July 1995, pages 1147-1170.
%
%  Parameters:
%
%    Input, integer D, the spatial dimension.
%
%    Input, integer N, the maximum total degree.
%
%    Input, integer R, the number of monomials in D dimensions
%    of total degree N or less.
%
%    Input, integer ND, the number of data points.
%    This function requires that the ND is equal to R.
%
%    Input, real XD(D,ND), the data points, which must be distinct.
%
%    Output, integer PO(ND), the order (number of nonzero coefficients) for the
%    Lagrange basis polynomials.
%
%    Output, real PC(ND,R), the coefficients for the
%    Lagrange basis polynomials.
%
%    Output, integer PE(ND,R), the exponent indices for the
%    Lagrange basis polynomials.
%

%
%  Verify that R is correct.
%
  if ( r ~= mono_upto_enum ( d, n ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LAGRANGE_COMPLETE - Fatal error!\n' );
    fprintf ( 1, '  The value R is not correct.\n' );
    error ( 'LAGRANGE_COMPLETE - Fatal error!' );
  end

  if ( r ~= nd )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LAGRANGE_COMPLETE - Fatal error!\n' );
    fprintf ( 1, '  The value R = %d does not equal ND = %d.\n', r, nd );
    error ( 'LAGRANGE_COMPLETE - Fatal error!' );
  end
%
%  Verify that the points are sufficiently distinct.
%
  [ d_min, d_max ] = r8col_separation ( d, nd, xd );
  d_tol = sqrt ( eps );

  if ( d_min < d_tol )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LAGRANGE_COMPLETE - Fatal error!\n' );
    fprintf ( 1, '  Some points are too close!\n' );
    fprintf ( 1, '  Minimum data point separation is = %g.\n', d_min );
    error ( 'LAGRANGE_COMPLETE - Fatal error!' );
  end
%
%  Initialize the polynomials Q to be all monomials of degree N or less.
%
  qo = zeros ( r, 1 );
  qc = zeros ( r, r );
  qe = zeros ( r, r );

  for k = 1 : r
    qo(k,1) = 1;
    qc(k,1) = 1.0;
    qe(k,1) = k;
  end
%
%  Now set up the P polynomials.
%
  po = zeros ( r, 1 );
  pc = zeros ( r, r );
  pe = zeros ( r, r );

  for k = 1 : nd
%
%  Find the first polynomial Q(K:R)(X) which is nonzero at X(K).
%
    i = r + 1;

    for j = k : r
      o = qo(j);
      value = polynomial_value ( d, o, qc(j,1:o), qe(j,1:o), 1, xd(1:d,k) );
      if ( value ~= 0.0 )
        i = j;
        break
      end
    end

    if ( i == r + 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'LAGRANGE_COMPLETE - Fatal error!\n' );
      fprintf ( 1, '  I = R+1.\n' );
      error ( 'LAGRANGE_COMPLETE - Fatal error!' );
    end
%
%  Define P(K)(X) = Q(I)(X) / Q(I)(X(k)
%
    o = qo(i);
    po(k) = qo(i);
    pc(k,1:o) = qc(i,1:o) / value;
    pe(k,1:o) = qe(i,1:o);
%
%  Modify P(1:k-1)(X).
%
    for j = 1 : k - 1

      oj = po(j);
      ok = po(k);

      value = polynomial_value ( d, oj, pc(j,1:oj), pe(j,1:oj), 1, xd(1:d,k) );

      [ o, c, e ] = polynomial_axpy ( - value, ok, pc(k,1:ok), pe(k,1:ok), ...
        oj, pc(j,1:oj), pe(j,1:oj) );

      po(j) = o;
      pc(j,1:o) = c(1:o)';
      pe(j,1:o) = e(1:o)';

    end
%
%  Modify Q(I:downto:K+1)
%
    for j = i : -1 : k + 1

      oj = qo(j-1);
      ok = po(k);

      value = polynomial_value ( d, oj, qc(j-1,1:oj), qe(j-1,1:oj), ...
        1, xd(1:d,k) );

      [ o, c, e ] = polynomial_axpy ( - value, ok, pc(k,1:ok), pe(k,1:ok), ...
        oj, qc(j-1,1:oj), qe(j-1,1:oj) );

      qo(j) = o;
      qc(j,1:o) = c(1:o)';
      qe(j,1:o) = e(1:o)';

    end
%
%  Modify Q(I+1:R)
%
    for j = i + 1 : r

      oj = qo(j);
      ok = po(k);

      value = polynomial_value ( d, oj, qc(j,1:oj), qe(j,1:oj), ...
        1, xd(1:d,k) );

      [ o, c, e ] = polynomial_axpy ( - value, ok, pc(k,1:ok), pe(k,1:ok), ...
        oj, qc(j,1:oj), qe(j,1:oj) );

      qo(j) = o;
      qc(j,1:o) = c(1:o)';
      qe(j,1:o) = e(1:o)';

    end

  end
%
%  Get rid of tiny coefficients.
%
  for i = 1 : nd
    oi = po(i);
    [ o, c, e ] = polynomial_compress ( ...
      po(i), pc(i,1:oi), pe(i,1:oi) );
    po(i) = o;
    pc(i,1:o) = c(1:o)';
    pe(i,1:o) = e(1:o)';
  end

  return
end


text:
function jsubs = cfg2jobsubs(item, subs)

% function jsubs = cfg2jobsubs(item, subs)
% Return the subscript into the job tree for a given subscript vector into
% the val part of the cfg tree. In a cfg_choice, this is a struct reference
% to a field with the name of the tag of the corresponding child node.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg2jobsubs.m 1862 2008-06-30 14:12:49Z volkmar $

rev = '$Rev: 1862 $'; %#ok

% Only de-reference subscripts into .val{1} field
if isequal(subs(1:2), substruct('.','val','{}',{1}))
    % return possible indices into .val{1} field
    citem = subsref(item, subs(1:2));
    if numel(subs) > 2
        jsubs1 = cfg2jobsubs(citem, subs(3:end));
    else
        jsubs1 = [];
    end;
    jsubs = [substruct('.', gettag(citem)) jsubs1];
else
    cfg_message('matlabbatch:cfg2jobsubs:wrongsubs', 'Inappropriate subscript reference in item ''%s''.', item.tag);
    jsubs = struct('type',{},'subs',{});
end;

author: spm


text:
function methods = get_method_configs()
% Get properties of all proposal methods. For examples, this includes paths to
% candidate files, intermediate files, and also information about how to sort
% candidates.

% If you want to add your own method, add it at the bottom.

colormap = [ ...
0, 0, 0; ...
228, 229, 97 ; ...
163, 163, 163 ; ...
218, 71, 56 ; ...
219, 135, 45 ; ...
145, 92, 146 ; ...
83, 136, 173 ; ...
135, 130, 174 ; ...
225, 119, 174 ; ...
142, 195, 129 ; ...
138, 180, 66 ; ...
223, 200, 51 ; ...
92, 172, 158 ; ...
177,89,40;
0, 255, 255;
188, 128, 189;
255, 255, 0;
] ./ 256;

 % config = get_config();
 % precomputed_prefix = config.precomputed_candidates;

methods = [];
  %{
   i = numel(methods) + 1;
  methods(i).name = 'LabelTransfer';
  methods(i).short_name = 'LT';
  prefix = ['/home/gneelima/work/data/intermediate/objectProposals/labelTransfer/nonflann'];
  methods(i).candidate_dir = [prefix ];
  methods(i).best_voc07_candidates_file = [prefix 'best_candidates.mat'];
  methods(i).best_imagenet_candidates_file = [prefix 'best_candidates_imagenet.mat'];
  methods(i).order = 'none';
  methods(i).abo_voc07_candidates_file=[prefix 'abo_candidates.mat'];
% methods(i).extract = @run_prop;
  methods(i).num_candidates = false;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'Sel.Search';
  methods(i).short_name = 'SS';
  prefix = ['/home/gneelima/Selective_search_imagenet/val/'];
  methods(i).candidate_dir = [prefix ];
  methods(i).best_voc07_candidates_file = [prefix 'best_candidates.mat'];
  methods(i).abo_voc07_candidates_file=[prefix 'abo_candidates.mat'];
  methods(i).order = 'ascend';
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;
%}
  i = numel(methods) + 1;
  methods(i).name = 'EdgeBoxes';
  methods(i).short_name = 'EB';
  prefix = ['/home/gneelima/work/data/output/objectProposals/edgeBoxes/imagenetBoxesIOU7/'];
  methods(i).candidate_dir = [prefix ];
  methods(i).best_voc07_candidates_file = [prefix 'best_candidates.mat'];
  methods(i).abo_voc07_candidates_file=[prefix 'abo_candidates.mat'];
  methods(i).order = 'descend';
  methods(i).extract = @run_edge_boxes;
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;
%{
  i = numel(methods) + 1;
  methods(i).name = 'Edge Boxes 50';
  methods(i).short_name = 'EB50';
  prefix = ['/home/gneelima/work/data/output/objectProposals/edgeBoxes/imagenetBoxesIOU5/'];
  methods(i).candidate_dir = [prefix ];
  methods(i).best_voc07_candidates_file = [prefix 'best_candidates.mat'];
  methods(i).abo_voc07_candidates_file=[prefix 'abo_candidates.mat'];
  methods(i).order = 'descend';
  methods(i).num_candidates = true;
  methods(i).color = colormap(i,:);
  methods(i).is_baseline = false;

  i = numel(methods) + 1;
  methods(i).name = 'Edge Boxes 90';
  methods(i).short_name = 'EB90';
  prefix = ['/home/gneelima/work/data/output/objectProposals/edgeBoxes/imagenetBoxesIOU9/'];
  methods(i).candidate_dir = [prefix ];
  methods(i).abo_voc07_candidates_file=[prefix 'abo_candidates.mat'];
  methods(i).best_imagenet_candidates_file = [prefix 'best_candidates_imagenet.mat'];
  methods(i).order = 'descend';
  methods(i).is_baseline = false;
 %}
  % add your own method here:
  if false
  i = numel(methods) + 1;
  methods(i).name = 'The full name of your method';
  methods(i).short_name = 'a very short version of the name';
  prefix = [precomputed_prefix 'ours-wip/'];
  methods(i).candidate_dir = [prefix ];
  methods(i).repeatability_candidate_dir = [prefix 'repeatability_mat'];
  methods(i).best_voc07_candidates_file = [prefix 'best_candidates.mat'];
  methods(i).abo_voc07_candidates_file=[prefix 'abo_candidates.mat'];
  methods(i).best_imagenet_candidates_file = [prefix 'best_candidates_imagenet.mat'];
  methods(i).repeatability_matching_file = [prefix 'repeatability_matching.mat'];
  % This specifies how to order candidates so that the first n, are the best n
  % candidates. For example we run a method for 10000 candidates and then take
  % the first 10, instead of rerunning for 10 candidates. Valid orderings are:
  %   none: candidates are already sorted, do nothing
  %   ascend/descend: sort by score descending or ascending
  %   random: random order
  %   biggest/smallest: sort by size of the bounding boxes
  methods(i).order = 'descend';
  % A function pointer to a method that runs your proposal detector.
  methods(i).extract = @run_edge_boxes90;
  % If your method supports sorting this should be empty. If your method has to
  % be rerun for every number of candidates we want, specify the number of
  % candidates here:
  methods(i).rerun_num_candidates = []; % ceil(10 .^ (2:0.5:4));
  % Specifies whether or not your method takes the desired number of candidates
  % as an input.
  % TODO(hosang): Is this actually used anywhere?
  methods(i).num_candidates = true;
  % color for drawing
  methods(i).color = colormap(i,:);
  % This should be false. Is used for drawing baselines dashed.
  methods(i).is_baseline = false;
  end

  % do the sorting dance
  sort_keys = [num2cell([methods.is_baseline])', {methods.name}'];
  for i = 1:numel(methods)
    sort_keys{i,1} = sprintf('%d', sort_keys{i,1});
  end
  [~,idx] = sortrows(sort_keys);
  for i = 1:numel(methods)
    methods(idx(i)).sort_key = i;
  end
end

author: Cloud-CV

repo: object-proposals


text:
function chi = chi_measure ( dim_num, n, z, ns, sample_routine, seed_init )

%*****************************************************************************80
%
%% CHI_MEASURE determines the pointset quality measure CHI.
%
%  Discussion:
%
%    The CHI measure of point distribution quality for a set Z of
%    N points in an DIM_NUM-dimensional region is defined as follows:
%
%    Assign every point X in the region to the nearest element
%    Z(I) of the point set.  For each Z(I), let H(I) be the maximum
%    distance between Z(I) and any point assigned to it by this process.
%
%    For each point Z(I), we determine the nearest distinct element of
%    the pointset by
%
%      GAMMA(I) = minimum ( 1 <= J <= N, I /= J ) distance ( Z(I), Z(J) )
%
%    Then
%
%      CHI(I) = 2 * H(I) / GAMMA(I)
%
%    and
%
%      CHI = maximum ( 1 <= I <= N ) CHI(I)
%
%    This quantity can be estimated by using sampling to pick a large
%    number of points in the region, rather than all of them.
%
%    For an ideally regular mesh, all the CHI(I)'s will be equal.
%    Any deviation from regularity increases the value of some entries
%    of CHI; thus, given two meshes, the one with a lower value of
%    CHI is to be recommended.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2004
%
%  Author:
%
%    Max Gunzburger
%    John Burkardt
%
%  Reference:
%
%    Max Gunzburger and John Burkardt,
%    Uniformity Measures for Point Samples in Hypercubes.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real ( kind = 8 ) Z(DIM_NUM,N), the points.
%
%    Input, integer NS, the number of sample points.
%
%    Input, character SAMPLE_ROUTINE, the name of the function to be used to
%    return sample points in the region.  The function must have the form
%      [ x, seed ] = sample_routine ( dim_num, n, seed )
%
%    Input, integer SEED_INIT, the initial value of the random number seed.
%
%    Output, real CHI, the CHI quality measure.
%
  seed = seed_init;
  h(1:n) = 0.0;

  for k = 1 : ns

    [ x, seed ] = feval ( sample_routine, dim_num, 1, seed );

    closest = find_closest ( dim_num, n, 1, x, z );

    dist = sum ( ( x(1:dim_num) - z(1:dim_num,closest(1)) ).^2 );

    h(closest(1)) = max ( h(closest(1)), dist );

  end

  gamma = pointset_spacing ( dim_num, n, z );

  chi_vec(1:n) = 2.0 * sqrt ( h(1:n) ) ./ gamma(1:n);

  chi = max ( chi_vec(1:n) );

  return
end


text:
function [vclassifier, hclassifier] = mcmcTrainSegmentRegressor(features, labelprobs, weights)

ntrees = 20;
nnodes = 8;

maxdata = 30000;

for f = 1:numel(labelprobs)
    vprobs{f} = [labelprobs{f}(:, 1) sum(labelprobs{f}(:, 2:6), 2) labelprobs{f}(:, 7)];
    ind = find(vprobs{f}(:, 2)~=0);
    hprobs{f} = zeros(size(vprobs{f}, 1), 5);
    hprobs{f}(ind, :) = labelprobs{f}(ind, 2:6) ./ repmat(vprobs{f}(ind, 2), [1 5]);
end

[vdata, vlab, vw] = formatData(features, vprobs, weights, maxdata);
disp(numel(vw))
[hdata, hlab, hw] = formatData(features, hprobs, weights, maxdata);
disp(numel(hw))

vnames = {'000', '090', 'sky'};
hnames = {'045', '090', '135', 'por', 'sol'};

catids = 75;

vclassifier =  train_boosted_dt_mc(vdata, catids, vnames(vlab)', ntrees, nnodes, vw, vnames);

hclassifier =  train_boosted_dt_mc(hdata, catids, hnames(hlab)', ntrees, nnodes, hw, hnames);

function [data, lab, w] = formatData(features, probs, weights, maxdata)
% concatenate data

nimages = numel(features);
nlabels = size(probs{1}, 2);

[tmp, nvars] = size(features{1});

% count segments
nseg = 0;
for f = 1:nimages
    nseg = nseg + size(features{f}, 1);
end
disp(num2str(nseg))

data = zeros(nseg, nvars);
labprob = zeros(nseg, nlabels);
w = zeros(nseg, 1);

% concatenate data
vc = 0;
for f = 1:nimages
    ind = [1:size(features{f}, 1)];
    data(vc+1:vc+numel(ind), :) = features{f}(ind, :);
    labprob(vc+1:vc+numel(ind), :) = probs{f}(ind, :);
    w(vc+1:vc+numel(ind)) = weights{f}(ind);
    vc = vc + numel(ind);
%     for y = 1:nlabels
%         data(vc+1:vc+numel(ind), :) = features{f}(ind, :);
%         lab(vc+1:vc+numel(ind)) = y;
%         % weight according to (image area) * (label percentage)
%         w(vc+1:vc+numel(ind)) = weights{f}(ind).*probs{f}(ind, y);
%         vc = vc + numel(ind);
%     end
end

if nseg > maxdata
    rind = randperm(nseg);
    data = data(rind(1:maxdata), :);
    labprob = labprob(rind(1:maxdata), :);
    w = w(rind(1:maxdata));
    nseg = maxdata;
end
tmpw = w;

w = zeros(nseg*nlabels, 1);
lab = zeros(nseg*nlabels, 1);

data = repmat(data, [nlabels 1]);
for y = 1:nlabels
    lab((y-1)*nseg+1:y*nseg) = y;
    w((y-1)*nseg+1:y*nseg) = tmpw .* labprob(:, y);
end

% remove zero weight data
ind = find(w==0);
data(ind, :) = [];
lab(ind) = [];
w(ind) = [];

w = w / sum(w);

author: zouchuhang

repo: LayoutNet


text:
function [W] = MW2W(MW)
% Convert power from megawatts to watts.
% Chad A. Greene 2012
W = MW*1000000 ;

author: Sable

repo: mcbench-benchmarks


text:
function varargout=lpcrf2lo(varargin)
%
% For calling details please see v_lpcrf2lo.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcrf2lo.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcrf2lo.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcrf2lo(varargin{:});
else
    v_lpcrf2lo(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function center = graphCenter(v, e, l)
%GRAPHCENTER Center of a graph
%
%   CENTER = graphCenter(V, E)
%   Computes the center of the graph given by V and E. The center of the
%   graph is the set of vertices whose eccentricity is minimal. The
%   function returns indices of center vertices.
%
%   CENTER = graphCenter(V, E, L)
%   Specifies the weight of each edge for computing the distances. Default
%   is to consider a weight of 1 for each edge.
%
%   Example
%     nodes = [20 20;20 50;20 80;50 50;80 20;80 50;80 80];
%     edges = [1 2;2 3;2 4;4 6;5 6;6 7];
%     figure; drawGraph(nodes, edges);
%     axis([0 100 0 100]); axis equal; hold on
%     C = graphCenter(nodes, edges)
%     C =
%         4
%
%   See Also
%   grPropagateDistance, grVertexEccentricity
%   graphRadius, graphDiameter, graphPeripheralVertices
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-09-07,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% ensure there is a valid length array
if nargin<3
    l = ones(size(e,1), 1);
end

g = grVertexEccentricity(v, e, l);

center = find(g==min(g));

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function fact_test ( )

%*****************************************************************************80
%
%% FACT_TEST tests the use of the MEX file FACT.C
%
%  Discussion:
%
%    The file fact.c is a C function which computes the factorial.
%
%    This M file "compiles" fact.c, and then shows how it can be called.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2006
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FACT_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate a simple use of the MEX compiler,\n' );
  fprintf ( 1, '  which allows MATLAB to call C functions.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Get a directory listing.  The file FACT.C should be there.\n' );

  ls

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compile the file FACT.C.\n' );

  mex fact.c

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Get a directory listing.  A new file should show up,\n' );
  fprintf ( 1, '  containing the compiled information.\n' );

  ls

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now use FACT as though it were a MATLAB M-file function.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   N  (N Factorial)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    j = fact ( i );

    fprintf ( 1, '  %2d  %10d\n', i, j );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FACT_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( )

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end