text:
function [source] = ft_sourceanalysis(cfg, data, baseline)

% FT_SOURCEANALYSIS performs beamformer dipole analysis on EEG or MEG data
% after preprocessing and a timelocked or frequency analysis
%
% Use as
%   [source] = ft_sourceanalysis(cfg, freq)
% or
%   [source] = ft_sourceanalysis(cfg, timelock)
%
% where the second input argument with the data should be organised in a structure
% as obtained from the FT_FREQANALYSIS or FT_TIMELOCKANALYSIS function. The
% configuration "cfg" is a structure containing information about source positions
% and other options.
%
% The different source reconstruction algorithms that are implemented are
%   cfg.method     = 'lcmv'    linear constrained minimum variance beamformer
%                    'sam'     synthetic aperture magnetometry
%                    'dics'    dynamic imaging of coherent sources
%                    'pcc'     partial canonical correlation/coherence
%                    'mne'     minimum norm estimation
%                    'rv'      scan residual variance with single dipole
%                    'music'   multiple signal classification
%                    'sloreta' standardized low-resolution electromagnetic tomography
%                    'eloreta' exact low-resolution electromagnetic tomography
% The DICS and PCC methods are for frequency or time-frequency domain data, all other
% methods are for time domain data. ELORETA can be used both for time, frequency and
% time-frequency domain data.
%
% The complete grid with dipole positions and optionally precomputed leadfields is
% constructed using FT_PREPARE_SOURCEMODEL. It can be specified as as a regular 3-D
% grid that is aligned with the axes of the head coordinate system using
%   cfg.xgrid               = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.ygrid               = vector (e.g. -20:1:20) or 'auto' (default = 'auto')
%   cfg.zgrid               = vector (e.g.   0:1:20) or 'auto' (default = 'auto')
%   cfg.resolution          = number (e.g. 1 cm) for automatic grid generation
% If the source model destribes a triangulated cortical sheet, it is described as
%   cfg.sourcemodel.pos     = N*3 matrix with the vertex positions of the cortical sheet
%   cfg.sourcemodel.tri     = M*3 matrix that describes the triangles connecting the vertices
% Alternatively the position of a few dipoles at locations of interest can be
% user-specified, for example obtained from an anatomical or functional MRI
%   cfg.sourcemodel.pos     = N*3 matrix with position of each source
%   cfg.sourcemodel.inside  = N*1 vector with boolean value whether grid point is inside brain (optional)
%   cfg.sourcemodel.dim     = [Nx Ny Nz] vector with dimensions in case of 3-D grid (optional)
%
% Besides the source positions, you may also include previously computed
% spatial filters and/or leadfields using
%   cfg.sourcemodel.filter
%   cfg.sourcemodel.leadfield
%
% The following strategies are supported to obtain statistics for the source parameters using
% multiple trials in the data, either directly or through a resampling-based approach
%   cfg.rawtrial      = 'no' or 'yes'   construct filter from single trials, apply to single trials. Note that you also may want to set cfg.keeptrials='yes' to keep all trial information, especially if using in combination with sourcemodel.filter
%   cfg.jackknife     = 'no' or 'yes'   jackknife resampling of trials
%   cfg.pseudovalue   = 'no' or 'yes'   pseudovalue resampling of trials
%   cfg.bootstrap     = 'no' or 'yes'   bootstrap resampling of trials
%   cfg.numbootstrap  = number of bootstrap replications (e.g. number of original trials)
% If none of these options is specified, the average over the trials will
% be computed prior to computing the source reconstruction.
%
% To obtain statistics over the source parameters between two conditions, you
% can also use a resampling procedure that reshuffles the trials over both
% conditions. In that case, you should call the function with two datasets
% containing single trial data like
%   [source] = ft_sourceanalysis(cfg, freqA, freqB)
%   [source] = ft_sourceanalysis(cfg, timelockA, timelockB)
% and you should specify
%   cfg.randomization      = 'no' or 'yes'
%   cfg.permutation        = 'no' or 'yes'
%   cfg.numrandomization   = number, e.g. 500
%   cfg.numpermutation     = number, e.g. 500 or 'all'
%
% If you have not specified a sourcemodel with pre-computed leadfields, the leadfield
% for each source position will be computed on the fly. In that case you can modify
% the leadfields by reducing the rank (i.e. remove the weakest orientation), or by
% normalizing each column.
%   cfg.reducerank      = 'no', or number (default = 3 for EEG, 2 for MEG)
%   cfg.backproject     = 'yes' or 'no',  determines when reducerank is applied whether the
%                         lower rank leadfield is projected back onto the original linear
%                         subspace, or not (default = 'yes')
%   cfg.normalize       = 'yes' or 'no' (default = 'no')
%   cfg.normalizeparam  = depth normalization parameter (default = 0.5)
%   cfg.weight          = number or Nx1 vector, weight for each dipole position to compensate
%                         for the size of the corresponding patch (default = 1)
%
% Other configuration options are
%   cfg.channel       = Nx1 cell-array with selection of channels (default = 'all'), see FT_CHANNELSELECTION for details
%   cfg.frequency     = single number (in Hz)
%   cfg.latency       = single number in seconds, for time-frequency analysis
%   cfg.lambda        = number or empty for automatic default
%   cfg.kappa         = number or empty for automatic default
%   cfg.tol           = number or empty for automatic default
%   cfg.refchan       = reference channel label (for coherence)
%   cfg.refdip        = reference dipole location (for coherence)
%   cfg.supchan       = suppressed channel label(s)
%   cfg.supdip        = suppressed dipole location(s)
%   cfg.keeptrials    = 'no' or 'yes'
%   cfg.keepleadfield = 'no' or 'yes'
%   cfg.projectnoise  = 'no' or 'yes'
%   cfg.keepfilter    = 'no' or 'yes'
%   cfg.keepcsd       = 'no' or 'yes'
%   cfg.keepmom       = 'no' or 'yes'
%   cfg.feedback      = 'no', 'text', 'textbar', 'gui' (default = 'text')
%
% The volume conduction model of the head should be specified as
%   cfg.headmodel     = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%
% The EEG or MEG sensor positions can be present in the data or can be specified as
%   cfg.elec          = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.grad          = structure with gradiometer definition or filename, see FT_READ_SENS
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_SOURCEDESCRIPTIVES, FT_SOURCESTATISTICS, FT_PREPARE_LEADFIELD,
% FT_PREPARE_HEADMODEL, FT_PREPARE_SOURCEMODEL

% Undocumented local options:
% cfg.numcomponents
% cfg.refchannel
% cfg.trialweight   = 'equal' or 'proportional'
% cfg.powmethod     = 'lambda1' or 'trace'

% Copyright (c) 2003-2008, F.C. Donders Centre, Robert Oostenveld
%
% 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;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data baseline
ft_preamble provenance data baseline

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% the baseline data can be passed as input argument or can be read from disk
hasbaseline = exist('baseline', 'var');

% check if the input data is valid for this function
data = ft_checkdata(data, 'datatype', {'timelock', 'freq', 'comp'}, 'feedback', 'yes');

if hasbaseline
  baseline = ft_checkdata(baseline, 'datatype', {'timelock', 'freq', 'comp'}, 'feedback', 'yes');
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'forbidden',  {'parallel', 'trials'});
cfg = ft_checkconfig(cfg, 'forbidden',  {'foi', 'toi'});
cfg = ft_checkconfig(cfg, 'renamed',    {'toilim', 'latency'});
cfg = ft_checkconfig(cfg, 'renamed',    {'foilim', 'frequency'});
cfg = ft_checkconfig(cfg, 'renamed',    {'jacknife', 'jackknife'});
cfg = ft_checkconfig(cfg, 'renamed',    {'refchannel', 'refchan'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'power',           'dics'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'coh_refchan',     'dics'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'coh_refdip',      'dics'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'dics_cohrefchan', 'dics'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'dics_cohrefdip',  'dics'});
cfg = ft_checkconfig(cfg, 'renamed',    {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'vol',     'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'grid',    'sourcemodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'elecfile', 'elec'});
cfg = ft_checkconfig(cfg, 'renamed',    {'gradfile', 'grad'});
cfg = ft_checkconfig(cfg, 'renamed',    {'optofile', 'opto'});

% determine the type of input data
istimelock = ft_datatype(data, 'timelock');
isfreq     = ft_datatype(data, 'freq');
iscomp     = ft_datatype(data, 'comp');
if ~any([isfreq iscomp istimelock])
  ft_error('input data is not recognized');
end

% set the defaults
if istimelock
  cfg.method = ft_getopt(cfg, 'method', 'lcmv');
elseif isfreq
  cfg.method = ft_getopt(cfg, 'method', 'dics');
end

% ensure that the method is specified
ft_checkopt(cfg, 'method', 'char');

if isequal(cfg.method, 'harmony')
  ft_error('The harmony implementation does not work at present. Please contact the main developer of this method directly');
end

% put the low-level options pertaining to the source reconstruction method in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg',  cfg.method);
% move some fields from cfg.method back to the top-level configuration
cfg = ft_checkconfig(cfg, 'createtopcfg', cfg.method);

% put the low-level options pertaining to the dipole grid in their own field
cfg = ft_checkconfig(cfg, 'renamed', {'tightgrid', 'tight'});  % this is moved to cfg.sourcemodel.tight by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'renamed', {'sourceunits', 'unit'}); % this is moved to cfg.sourcemodel.unit  by the subsequent createsubcfg

% put the low-level options pertaining to the sourcemodel in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg', {'sourcemodel'});
% move some fields from cfg.sourcemodel back to the top-level configuration
cfg = ft_checkconfig(cfg, 'createtopcfg', {'sourcemodel'});

% get the low-level options for the inverse estimation method, these are method specific
cfg.(cfg.method)               = ft_getopt(cfg, cfg.method);
cfg.(cfg.method).keepfilter    = ft_getopt(cfg.(cfg.method), 'keepfilter',    'no');
cfg.(cfg.method).keepcsd       = ft_getopt(cfg.(cfg.method), 'keepcsd',       'no');
cfg.(cfg.method).keepmom       = ft_getopt(cfg.(cfg.method), 'keepmom',       'yes');
cfg.(cfg.method).projectnoise  = ft_getopt(cfg.(cfg.method), 'projectnoise',  'no');
cfg.(cfg.method).feedback      = ft_getopt(cfg.(cfg.method), 'feedback',      'text');
cfg.(cfg.method).lambda        = ft_getopt(cfg.(cfg.method), 'lambda',        []);
cfg.(cfg.method).kappa         = ft_getopt(cfg.(cfg.method), 'kappa',         []);
cfg.(cfg.method).tol           = ft_getopt(cfg.(cfg.method), 'tol',           []);
cfg.(cfg.method).invmethod     = ft_getopt(cfg.(cfg.method), 'invmethod',     []);
cfg.(cfg.method).powmethod     = ft_getopt(cfg.(cfg.method), 'powmethod',     []);

% get any further options
cfg.keepleadfield    = ft_getopt(cfg, 'keepleadfield',  'no');
cfg.keeptrials       = ft_getopt(cfg, 'keeptrials',     'no');
cfg.trialweight      = ft_getopt(cfg, 'trialweight',    'equal');
cfg.jackknife        = ft_getopt(cfg, 'jackknife',      'no');
cfg.pseudovalue      = ft_getopt(cfg, 'pseudovalue',    'no');
cfg.bootstrap        = ft_getopt(cfg, 'bootstrap',      'no');
cfg.singletrial      = ft_getopt(cfg, 'singletrial',    'no');
cfg.rawtrial         = ft_getopt(cfg, 'rawtrial',       'no');
cfg.randomization    = ft_getopt(cfg, 'randomization',  'no');
cfg.numrandomization = ft_getopt(cfg, 'numrandomization', 100);
cfg.permutation      = ft_getopt(cfg, 'permutation',    'no');
cfg.numpermutation   = ft_getopt(cfg, 'numpermutation', 100);
cfg.wakewulf         = ft_getopt(cfg, 'wakewulf',       'yes');
cfg.killwulf         = ft_getopt(cfg, 'killwulf',       'yes');
cfg.channel          = ft_getopt(cfg, 'channel',        'all');
cfg.latency          = ft_getopt(cfg, 'latency',        'all');
cfg.frequency        = ft_getopt(cfg, 'frequency',      'all');
% these only apply to DICS and PCC
cfg.refdip           = ft_getopt(cfg, 'refdip', []);
cfg.supdip           = ft_getopt(cfg, 'supdip', []);
cfg.refchan          = ft_getopt(cfg, 'refchan', []);
cfg.supchan          = ft_getopt(cfg, 'supchan', []);

if hasbaseline && (strcmp(cfg.randomization, 'no') && strcmp(cfg.permutation, 'no'))
  ft_error('input of two conditions only makes sense if you want to randomize or permute');
elseif ~hasbaseline && (strcmp(cfg.randomization, 'yes') || strcmp(cfg.permutation, 'yes'))
  ft_error('randomization or permutation requires that you give two conditions as input');
end

if sum([strcmp(cfg.jackknife, 'yes'), strcmp(cfg.bootstrap, 'yes'), strcmp(cfg.pseudovalue, 'yes'), strcmp(cfg.singletrial, 'yes'), strcmp(cfg.rawtrial, 'yes'), strcmp(cfg.randomization, 'yes'), strcmp(cfg.permutation, 'yes')])>1
  ft_error('jackknife, bootstrap, pseudovalue, singletrial, rawtrial, randomization and permutation are mutually exclusive');
end

if strcmp(cfg.rawtrial, 'yes') && isfield(cfg, 'sourcemodel') && ~isfield(cfg.sourcemodel, 'filter')
  ft_warning('Using each trial to compute its own filter is not currently recommended. Use this option only with precomputed filters in cfg.sourcemodel.filter');
end

if ~isempty(cfg.refchan)
  cfg.refchan = ft_channelselection(cfg.refchan, data.label);
  assert(numel(cfg.refchan)>0, 'cfg.refchan is not present in the data');
end
if ~isempty(cfg.supchan)
  cfg.supchan = ft_channelselection(cfg.supchan, data.label);
  assert(numel(cfg.supchan)>0, 'cfg.supchan is not present in the data');
end

% spectrally decomposed data can have label and/or labelcmb
if ~isfield(data, 'label') && isfield(data, 'labelcmb')
  % the code further down assumes that data.label is present
  % we can construct it from all channel combinations
  data.label = unique(data.labelcmb(:));
end

% make the selection of channels consistent with the data
cfg.channel = ft_channelselection(cfg.channel, data.label);
% keep the refchan and supchan
cfg.channel = ft_channelselection([cfg.channel(:); cfg.refchan(:)], data.label);
cfg.channel = ft_channelselection([cfg.channel(:); cfg.supchan(:)], data.label);

% start with an empty output structure
source = [];

if istimelock
  tmpcfg = keepfields(cfg, {'channel', 'latency', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  % keep the time axis in the output
  tmpcfg.avgovertime = 'no';
  data = ft_selectdata(tmpcfg, data);
  % restore the provenance information
  [cfg, data] = rollback_provenance(cfg, data);

  % copy the descriptive fields to the output
  source = copyfields(data, source, {'time'});

elseif isfreq
  tmpcfg = keepfields(cfg, {'channel', 'latency', 'frequency', 'nanmean', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});

  if ismember(cfg.method, {'pcc' 'dics'})
    tmpcfg.avgoverfreq = 'yes';
    if isfield(data, 'time')
      tmpcfg.avgovertime = 'yes';
    end
  end
  % include the refchan and supchan if specified
  tmpcfg.channel = ft_channelselection([cfg.channel(:); cfg.refchan(:); cfg.supchan(:)], data.label);
  data = ft_selectdata(tmpcfg, data);
  % restore the provenance information
  [cfg, data] = rollback_provenance(cfg, data);

  % copy the descriptive fields to the output
  source = copyfields(data, source, {'time', 'freq', 'cumtapcnt'});

  if ismember(cfg.method, {'pcc' 'dics'})
    cfg.frequency = data.freq; % should be a single number here
    if isfield(data, 'time')
      cfg.latency   = data.time; % should be a single number here
    end

  end
end

if isfreq && isfield(data, 'labelcmb')
  % ensure that the cross-spectral densities are chan_chan_therest,
  % otherwise the latency and frequency selection could fail, so we don't
  % need to worry about linearly indexed cross-spectral densities below
  % this point, this step may take some time, if multiple trials are
  % present in the data
  fprintf('converting the linearly indexed channelcombinations into a square CSD-matrix\n');
  data = ft_checkdata(data, 'cmbstyle', 'full');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prepare the sourcemodel, headmodel, sensors and/or leadfields
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isempty(cfg.refdip) || ~isempty(cfg.supdip)
  ft_notice('computing the leadfields on the fly');

  % the leadfields for refdip and supdip have to be computed on the fly; to ensure
  % that the leadfields for the grid are consistent, these will also be computed on
  % the fly

  if isfield(cfg.sourcemodel, 'leadfield')
    ft_warning('ignoring the precomputed leadfield that are provided');
    cfg.sourcemodel = rmfield(cfg.sourcemodel, 'leadfield');
  end
  if isfield(cfg.sourcemodel, 'filter')
    ft_warning('ignoring the precomputed filters that are provided');
    cfg.sourcemodel = rmfield(cfg.sourcemodel, 'filter');
  end

  % collect and preprocess the electrodes/gradiometer and head model
  [headmodel, sens, cfg] = prepare_headmodel(cfg, data);

  % construct the dipole positions on which the source reconstruction will be done
  tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid', 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'reducerank', 'backproject', 'normalize', 'normalizeparam', 'weight', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.headmodel = headmodel;
  if ft_senstype(sens, 'eeg')
    tmpcfg.elec = sens;
  elseif ft_senstype(sens, 'meg')
    tmpcfg.grad = sens;
  end
  sourcemodel = ft_prepare_sourcemodel(tmpcfg);

elseif isfield(cfg.sourcemodel, 'filter')
  ft_notice('using precomputed filters, not computing any leadfields');
  sourcemodel = keepfields(cfg.sourcemodel, {'pos', 'tri', 'dim', 'inside', 'filter', 'filterdimord', 'label', 'cfg'});

  if ~isfield(sourcemodel, 'label')
    ft_warning('the labels are missing for the precomputed filters, assuming that they were computed with the same channel selection');
    sourcemodel.label = cfg.channel;
  end

  % select the channels corresponding to the data and the user configuration
  tmpcfg = keepfields(cfg, 'channel');
  sourcemodel = ft_selectdata(tmpcfg, sourcemodel);

  % sort the channels to be consistent with the data
  [dum, chansel] = match_str(data.label, sourcemodel.label);
  sourcemodel.label = sourcemodel.label(chansel);
  for i=1:numel(sourcemodel.filter)
    if ~isempty(sourcemodel.filter{i})
      sourcemodel.filter{i} = sourcemodel.filter{i}(:, chansel);
    end
  end

  % ensure that the channels are consistent with the data
  if isempty(ft_getopt(cfg, 'refchan')) && isempty(ft_getopt(cfg, 'supchan'))
    assert(isequal(sourcemodel.label(:), cfg.channel(:)), 'cannot match the channels in the sourcemodel to those in the data');
  else
    % the data and cfg also includes the recfchan or supchan
    assert(all(ismember(sourcemodel.label, cfg.channel)), 'cannot match the channels in the sourcemodel to those in the data');
  end

  % no forward computations are needed
  headmodel = [];
  sens = [];
  cfg = removefields(cfg, {'headmodel', 'elec', 'grad'});

elseif isfield(cfg.sourcemodel, 'leadfield')
  ft_notice('using precomputed leadfields');
  sourcemodel = keepfields(cfg.sourcemodel, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label', 'cfg', 'subspace'});

  if ~isfield(sourcemodel, 'label')
    ft_warning('the labels are missing for the precomputed leadfields, assuming that they were computed with the same channel selection');
    sourcemodel.label = cfg.channel;
  end

  % select the channels corresponding to the data and the user configuration
  tmpcfg = keepfields(cfg, 'channel');
  sourcemodel = ft_selectdata(tmpcfg, sourcemodel);

  % sort the channels to be consistent with the data
  [dum, chansel] = match_str(data.label, sourcemodel.label);
  sourcemodel.label = sourcemodel.label(chansel);
  for i=1:numel(sourcemodel.leadfield)
    if ~isempty(sourcemodel.leadfield{i})
      sourcemodel.leadfield{i} = sourcemodel.leadfield{i}(chansel, :);
    end
  end

  % ensure that the channels are consistent with the data
  if isempty(ft_getopt(cfg, 'refchan')) && isempty(ft_getopt(cfg, 'supchan'))
    assert(isequal(sourcemodel.label(:), cfg.channel(:)), 'cannot match the channels in the sourcemodel to those in the data');
  else
    % the data and cfg also includes the recfchan or supchan
    assert(all(ismember(sourcemodel.label, cfg.channel)), 'cannot match the channels in the sourcemodel to those in the data');
  end

  % no forward computations are needed
  headmodel = [];
  sens = [];
  cfg = removefields(cfg, {'headmodel', 'elec', 'grad'});

elseif istrue(cfg.keepleadfield) || istrue(cfg.permutation) || istrue(cfg.randomization) || istrue(cfg.bootstrap) || istrue(cfg.jackknife) || istrue(cfg.pseudovalue) || istrue(cfg.singletrial) || istrue(cfg.rawtrial)
  ft_notice('computing the leadfields in advance');

  % collect and preprocess the electrodes/gradiometer and head model
  [headmodel, sens, cfg] = prepare_headmodel(cfg, data);

  % construct the dipole positions on which the source reconstruction will be done
  tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid', 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'reducerank', 'backproject', 'normalize', 'normalizeparam', 'weight', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.headmodel = headmodel;
  if ft_senstype(sens, 'eeg')
    tmpcfg.elec = sens;
  elseif ft_senstype(sens, 'meg')
    tmpcfg.grad = sens;
  end
  sourcemodel = ft_prepare_leadfield(tmpcfg);

  % no further forward computations are needed, but keep them in the cfg
  needheadmodel = false;
  headmodel = [];
  sens = [];

else
  ft_notice('computing the leadfields on the fly');

  % collect and preprocess the electrodes/gradiometer and head model
  [headmodel, sens, cfg] = prepare_headmodel(cfg, data);

  % construct the dipole positions on which the source reconstruction will be done
  tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid', 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'reducerank', 'backproject', 'normalize', 'normalizeparam', 'weight', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.headmodel = headmodel;
  if ft_senstype(sens, 'eeg')
    tmpcfg.elec = sens;
  elseif ft_senstype(sens, 'meg')
    tmpcfg.grad = sens;
  end
  sourcemodel = ft_prepare_sourcemodel(tmpcfg);

end % if refdip/supdip, precomputed filter, leadfield, keepfilter, keepleadfield, or so

% It might be that the number of channels in the data, the number of
% channels in the electrode/gradiometer definition and the number of
% channels in the localspheres volume conduction model are different.
% Hence a subset of the data channels will be used.
Nchans = length(cfg.channel);
if ~iscomp && contains(data.dimord, 'freq')
  Nfreq = numel(data.freq);
else
  Nfreq = 1;
end
if ~iscomp && contains(data.dimord, 'time')
  Ntime = numel(data.time);
else
  Ntime = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do frequency domain source reconstruction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfreq && any(strcmp(cfg.method, {'dics', 'pcc', 'eloreta', 'mne','harmony', 'rv', 'music'}))

  switch cfg.method
    case 'pcc'

      if hasbaseline
        ft_error('not supported')
      end

      tmpcfg         = keepfields(cfg, {'keeptrials', 'rawtrial', 'refchan', 'channel'});
      tmpcfg.refchan = []; % for PCC prepare_freq_matrices should not know explicitly about the refchan

      % select the data in the channels and the frequency of interest
      [Cf, Cr, Pr, Ntrials, tmpcfg] = prepare_freq_matrices(tmpcfg, data);

      if ~isempty(cfg.refchan)
        [dum, refchanindx] = match_str(cfg.refchan, tmpcfg.channel);
      else
        refchanindx = [];
      end
      if ~isempty(cfg.supchan)
        [dum, supchanindx] = match_str(cfg.supchan, tmpcfg.channel);
      else
        supchanindx = [];
      end
      Nchans = length(tmpcfg.channel); % update the number of channels

      % if the input data has a complete Fourier spectrum, it can be projected through the filters
      if isfield(data, 'fourierspctrm')
        [dum, datchanindx] = match_str(tmpcfg.channel, data.label);
        fbin = nearest(data.freq, cfg.frequency);
        if strcmp(data.dimord, 'chan_freq')
          avg = data.fourierspctrm(datchanindx, fbin);
        elseif strcmp(data.dimord, 'rpt_chan_freq') || strcmp(data.dimord, 'rpttap_chan_freq')
          avg = transpose(data.fourierspctrm(:, datchanindx, fbin));
        elseif strcmp(data.dimord, 'chan_freq_time')
          tbin = nearest(data.time, cfg.latency);
          avg = data.fourierspctrm(datchanindx, fbin, tbin);
        elseif strcmp(data.dimord, 'rpt_chan_freq_time') || strcmp(data.dimord, 'rpttap_chan_freq_time')
          tbin = nearest(data.time, cfg.latency);
          avg  = transpose(data.fourierspctrm(:, datchanindx, fbin, tbin));
        end
      else
        avg = [];
      end

    case {'eloreta' 'mne' 'rv' 'music' 'harmony'}
      % these can handle both a csd matrix and a fourier matrix
      tmpcfg = keepfields(cfg, {'keeptrials', 'rawtrial', 'refchan', 'channel'});
      [Cf, Cr, Pr, Ntrials, tmpcfg] = prepare_freq_matrices(tmpcfg, data);

      % if the input data has a complete fourier spectrum, it can be projected through the filters
      if isfield(data, 'fourierspctrm')
        [dum, datchanindx] = match_str(tmpcfg.channel, data.label);
        fbin = nearest(data.freq, cfg.frequency);
        if numel(fbin)==1, fbin = fbin.*[1 1]; end
        if strcmp(data.dimord, 'chan_freq')
          avg = data.fourierspctrm(datchanindx, fbin);
        elseif strcmp(data.dimord, 'rpt_chan_freq') || strcmp(data.dimord, 'rpttap_chan_freq')
          avg = permute(data.fourierspctrm(:, datchanindx, fbin(1):fbin(2)), [2 1 3]);
        elseif strcmp(data.dimord, 'chan_freq_time')
          tbin = nearest(data.time, cfg.latency);
          if numel(tbin)==1, tbin = tbin.*[1 1]; end
          avg = data.fourierspctrm(datchanindx, fbin(1):fbin(2), tbin(1):tbin(2));
        elseif strcmp(data.dimord, 'rpt_chan_freq_time') || strcmp(data.dimord, 'rpttap_chan_freq_time')
          tbin = nearest(data.time, cfg.latency);
          if numel(tbin)==1, tbin = tbin.*[1 1]; end
          avg  = permute(data.fourierspctrm(:, datchanindx, fbin(1):fbin(2), tbin(1):tbin(2)), [2 1 3 4]);
        end
      else % The input data is a CSD matrix, this is enough for computing source power, coherence and residual power.
        ft_warning('no fourierspctra in the input data, so the frequency domain dipole moments cannot be computed');
        avg = [];
      end

    case 'dics'
      tmpcfg         = keepfields(cfg, {'keeptrials', 'rawtrial', 'refchan', 'channel'});
      tmpcfg.channel = setdiff(cfg.channel, cfg.refchan, 'stable'); % remove the refchan, ensure that the ordering does not change, see https://github.com/fieldtrip/fieldtrip/issues/1587

      % select the data in the channels and the frequency of interest
      [Cf, Cr, Pr, Ntrials, tmpcfg] = prepare_freq_matrices(tmpcfg, data);

      Nchans = length(tmpcfg.channel); % update the number of channels

      % assign a descriptive name to each of the dics sub-methods, the default is power only
      if ~isempty(cfg.refdip)
        submethod = 'dics_refdip';
      elseif ~isempty(cfg.refchan)
        submethod = 'dics_refchan';
      else
        submethod = 'dics_power';
      end

    otherwise
      ft_error('unsupported cfg.method');
  end

  % fill these with NaNs, so that I dont have to treat them separately
  if isempty(Cr), Cr = nan(Ntrials, Nchans, Nfreq, Ntime); end
  if isempty(Pr), Pr = nan(Ntrials, Nfreq,  Ntime); end

  if hasbaseline
    % repeat the conversion for the baseline condition
    [bCf, bCr, bPr, Nbaseline, tmpcfg] = prepare_freq_matrices(tmpcfg, baseline);
    % fill these with NaNs, so that I dont have to treat them separately
    if isempty(bCr), bCr = nan(Nbaseline, Nchans, 1); end
    if isempty(bPr), bPr = nan(Nbaseline, 1, 1); end
    % rename the active condition for convenience
    aCf = Cf;
    aCr = Cr;
    aPr = Pr;
    % this is required for averaging 2 conditions using prepare_resampled_data
    cfg2 = [];
    cfg2.numcondition = 2;
    % this is required for randomizing/permuting 2 conditions using prepare_resampled_data
    cfg.numcondition = 2;
  end

  % prepare the resampling of the trials, or average the data if multiple trials are present and no resampling is necessary
  if (Ntrials<=1) && (strcmp(cfg.jackknife, 'yes') || strcmp(cfg.bootstrap, 'yes') || strcmp(cfg.pseudovalue, 'yes') || strcmp(cfg.singletrial, 'yes') || strcmp(cfg.rawtrial, 'yes') || strcmp(cfg.randomization, 'yes') || strcmp(cfg.permutation, 'yes'))
    ft_error('multiple trials required in the data\n');

  elseif strcmp(cfg.permutation, 'yes')
    % compute the cross-spectral density matrix without resampling
    [dum, avg_aCf, avg_aCr, avg_aPr, avg_bCf, avg_bCr, avg_bPr] = prepare_resampled_data(cfg2 , aCf, aCr, aPr, bCf, bCr, bPr);
    % compute the cross-spectral density matrix with random permutation
    [dum, rnd_aCf, rnd_aCr, rnd_aPr, rnd_bCf, rnd_bCr, rnd_bPr] = prepare_resampled_data(cfg, aCf, aCr, aPr, bCf, bCr, bPr);
    % concatenate the different resamplings
    Cf = cat(1, reshape(avg_aCf, [1 Nchans Nchans]), reshape(avg_bCf, [1 Nchans Nchans]), rnd_aCf, rnd_bCf);
    Cr = cat(1, reshape(avg_aCr, [1 Nchans 1     ]), reshape(avg_bCr, [1 Nchans 1     ]), rnd_aCr, rnd_bCr);
    Pr = cat(1, reshape(avg_aPr, [1 1      1     ]), reshape(avg_bPr, [1 1      1     ]), rnd_aPr, rnd_bPr);
    % clear temporary working copies
    clear avg_aCf avg_aCr avg_aPr avg_bCf avg_bCr avg_bPr
    clear rnd_aCf rnd_aCr rnd_aPr rnd_bCf rnd_bCr rnd_bPr
    % the order of the resamplings should be [avgA avgB rndA rndB rndA rndB rndA rndB ....]
    Nrepetitions = 2*cfg.numpermutation + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    Cf = Cf(order,:,:);
    Cr = Cr(order,:,:);
    Pr = Pr(order,:,:);

  elseif strcmp(cfg.randomization, 'yes')
    % compute the cross-spectral density matrix without resampling
    [dum, avg_aCf, avg_aCr, avg_aPr, avg_bCf, avg_bCr, avg_bPr] = prepare_resampled_data(cfg2 , aCf, aCr, aPr, bCf, bCr, bPr);
    % compute the cross-spectral density matrix with random resampling
    [dum, rnd_aCf, rnd_aCr, rnd_aPr, rnd_bCf, rnd_bCr, rnd_bPr] = prepare_resampled_data(cfg, aCf, aCr, aPr, bCf, bCr, bPr);
    % concatenate the different resamplings
    Cf = cat(1, reshape(avg_aCf, [1 Nchans Nchans]), reshape(avg_bCf, [1 Nchans Nchans]), rnd_aCf, rnd_bCf);
    Cr = cat(1, reshape(avg_aCr, [1 Nchans 1     ]), reshape(avg_bCr, [1 Nchans 1     ]), rnd_aCr, rnd_bCr);
    Pr = cat(1, reshape(avg_aPr, [1 1      1     ]), reshape(avg_bPr, [1 1      1     ]), rnd_aPr, rnd_bPr);
    % clear temporary working copies
    clear avg_aCf avg_aCr avg_aPr avg_bCf avg_bCr avg_bPr
    clear rnd_aCf rnd_aCr rnd_aPr rnd_bCf rnd_bCr rnd_bPr
    % the order of the resamplings should be [avgA avgB rndA rndB rndA rndB rndA rndB ....]
    Nrepetitions = 2*cfg.numrandomization + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    Cf = Cf(order,:,:);
    Cr = Cr(order,:,:);
    Pr = Pr(order,:,:);

  elseif strcmp(cfg.jackknife, 'yes')
    % compute the cross-spectral density matrix with jackknife resampling
    [cfg, Cf, Cr, Pr] = prepare_resampled_data(cfg, Cf, Cr, Pr);
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.bootstrap, 'yes')
    % compute the cross-spectral density matrix with bootstrap resampling
    [cfg, Cf, Cr, Pr] = prepare_resampled_data(cfg, Cf, Cr, Pr);
    Nrepetitions = cfg.numbootstrap;

  elseif strcmp(cfg.pseudovalue, 'yes')
    % compute the cross-spectral density matrix with pseudovalue resampling
    [cfg, Cf, Cr, Pr] = prepare_resampled_data(cfg, Cf, Cr, Pr);
    Nrepetitions = Ntrials+1;

  elseif strcmp(cfg.singletrial, 'yes')
    % The idea is that beamformer uses the average covariance to construct the
    % filter and applies it to the single trial covariance/csd. The problem
    % is that beamformer will use the averaged covariance/csd to estimate the
    % power and not the single trial covariance/csd
    ft_error('this option contains a bug, and is therefore not supported at the moment');
    Cf = Cf; % FIXME, should be averaged and repeated for each trial
    Cr = Cr; % FIXME, should be averaged and repeated for each trial
    Pr = Pr; % FIXME, should be averaged and repeated for each trial
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.rawtrial, 'yes')
    % keep all the individual trials, do not average them
    Cf = Cf;
    Cr = Cr;
    Pr = Pr;
    Nrepetitions = Ntrials;

  elseif Ntrials>1
    % compute the average from the individual trials
    Cf = reshape(sum(Cf, 1) / Ntrials, [Nchans Nchans]);
    Cr = reshape(sum(Cr, 1) / Ntrials, [Nchans 1]);
    Pr = reshape(sum(Pr, 1) / Ntrials, [1 1]);
    Nrepetitions = 1;

  elseif Ntrials==1
    % no rearrangement of trials is neccesary, the data already represents the average
    Cf = Cf;
    Cr = Cr;
    Pr = Pr;
    Nrepetitions = 1;
  end

  % reshape so that it also looks like one trial (out of many)
  if Nrepetitions==1
    Cf  = reshape(Cf , [1 Nchans Nchans Nfreq Ntime]);
    Cr  = reshape(Cr , [1 Nchans Nfreq Ntime]);
    Pr  = reshape(Pr , [1 Nfreq Ntime]);
  end

  % get the relevant low level options from the cfg and convert into key-value pairs
  tmpcfg = cfg.(cfg.method);
  % disable console feedback for the low-level function in case of multiple repetitions
  if Nrepetitions > 1
    tmpcfg.feedback = 'none';
  end
  methodopt = ft_cfg2keyval(tmpcfg);

  % construct the low-level options for the leadfield computation as key-value pairs, these are passed to the inverse function and FT_COMPUTE_LEADFIELD
  leadfieldopt = {};
  leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(cfg, 'reducerank'));
  leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(cfg, 'backproject'));
  leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(cfg, 'normalize'));
  leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(cfg, 'normalizeparam'));
  leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(cfg, 'weight'));

  if Nrepetitions > 1
    ft_progress('init', cfg.(cfg.method).feedback, 'scanning repetition...');
  end

  for i=1:Nrepetitions
    size_Cf    = size(Cf);
    squeeze_Cf = reshape(Cf(i,:,:), size_Cf(2:end));

    if Nrepetitions > 1
      ft_progress(i/Nrepetitions, 'scanning repetition %d from %d', i, Nrepetitions);
    end

    switch cfg.method
      case 'dics'
        if strcmp(submethod, 'dics_power')
          dip(i) = ft_inverse_dics(sourcemodel, sens, headmodel, [],  squeeze_Cf, methodopt{:}, leadfieldopt{:});
        elseif strcmp(submethod, 'dics_refchan')
          dip(i) = ft_inverse_dics(sourcemodel, sens, headmodel, [],  squeeze_Cf, methodopt{:}, leadfieldopt{:}, 'Cr', Cr(i,:), 'Pr', Pr(i));
        elseif strcmp(submethod, 'dics_refdip')
          dip(i) = ft_inverse_dics(sourcemodel, sens, headmodel, [],  squeeze_Cf, methodopt{:}, leadfieldopt{:}, 'refdip', cfg.refdip);
        end
      case 'pcc'
        if ~isempty(avg) && istrue(cfg.rawtrial)
          % FIXME added by jansch because an appropriate subselection of avg
          % should be done first (i.e. select the tapers that belong to this
          % repetition
          ft_error('rawtrial in combination with pcc has been temporarily disabled');
        else
          dip(i) = ft_inverse_pcc(sourcemodel, sens, headmodel, avg, squeeze_Cf, methodopt{:}, leadfieldopt{:}, 'refdip', cfg.refdip, 'refchan', refchanindx, 'supdip', cfg.supdip, 'supchan', supchanindx);
        end
      case 'eloreta'
        dip(i) = ft_inverse_eloreta(sourcemodel, sens, headmodel, avg, squeeze_Cf, methodopt{:}, leadfieldopt{:});
      case 'mne'
        dip(i) = ft_inverse_mne(sourcemodel, sens, headmodel, avg, methodopt{:}, leadfieldopt{:});
      case 'harmony'
        dip(i) = ft_inverse_harmony(sourcemodel, sens, headmodel, avg, methodopt{:}, leadfieldopt{:});
        % ft_error(sprintf('method ''%s'' is unsupported for source reconstruction in the frequency domain', cfg.method));
      case {'rv'}
        dip(i) = ft_inverse_rv(sourcemodel, sens, headmodel, avg, methodopt{:}, leadfieldopt{:});
      case {'music'}
        ft_error('method ''%s'' is currently unsupported for source reconstruction in the frequency domain', cfg.method);
      otherwise
        ft_error('unsupported cfg.method');
    end

  end
  if Nrepetitions > 1
    ft_progress('close');
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % do time domain source reconstruction
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif istimelock && any(strcmp(cfg.method, {'lcmv', 'sam', 'mne', 'harmony', 'rv', 'music', 'pcc', 'mvl', 'sloreta', 'eloreta'}))

  % determine the size of the data
  Nsamples = length(data.time);
  Nchans   = length(data.label);
  if isfield(data, 'cov') && length(size(data.cov))==3
    Ntrials = size(data.cov,1);
  elseif isfield(data, 'trial') && length(size(data.trial))==3
    Ntrials = size(data.trial,1);
  else
    Ntrials = 1;
  end

  if isfield(data, 'cov')
    % use the estimated data covariance matrix
    hascovariance = true;
  else
    % add a identity covariance matrix, this simplifies the handling of the different source reconstruction methods
    % since the covariance is only used by some reconstruction methods and might not always be present in the data
    if Ntrials==1
      data.cov = eye(Nchans);
    else
      data.cov = zeros(Ntrials,Nchans,Nchans);
      for i=1:Ntrials
        data.cov(i,:,:) = eye(Nchans);
      end
    end
    hascovariance = false;
    ft_warning('No covariance matrix found, assuming identity covariance matrix');
  end

  if strcmp(cfg.method, 'pcc')

    if hasbaseline
      ft_error('not supported')
    end

    tmpcfg = [];
    tmpcfg.channel = cfg.channel;
    if ~isempty(cfg.refchan)
      tmpcfg.channel = [tmpcfg.channel(:); cfg.refchan(:)];
    end
    if ~isempty(cfg.supchan)
      tmpcfg.channel = [tmpcfg.channel cfg.supchan(:)'];
    end

    % select the data in the channels of interest
    [dum, datchanindx] = match_str(tmpcfg.channel, data.label);
    if Ntrials==1
      data.avg = data.avg(datchanindx,:);
      data.cov = data.cov(datchanindx,datchanindx);
    else
      data.cov   = data.cov(:,datchanindx,datchanindx);
      data.trial = data.trial(:,datchanindx,:);
    end
    data.label = data.label(datchanindx);

    if ~isempty(cfg.refchan)
      [dum, refchanindx] = match_str(cfg.refchan, data.label);
    else
      refchanindx = [];
    end
    if ~isempty(cfg.supchan)
      [dum, supchanindx] = match_str(cfg.supchan, data.label);
    else
      supchanindx = [];
    end
    Nchans = length(tmpcfg.channel); % update the number of channels

  else
    % HACK: use the default code
    % select the channels of interest
    [dum, datchanindx] = match_str(cfg.channel, data.label);
    if strcmp(data.dimord, 'chan_time')
      % It is in principle possible to have timelockanalysis with
      % keeptrial=yes and only a single trial in the raw data.
      % In that case the covariance should be represented as Nchan*Nchan
      data.avg = data.avg(datchanindx,:);
      %data.cov = reshape(data.cov, length(datchanindx), length(datchanindx));
      data.cov = data.cov(datchanindx,datchanindx);
    elseif strcmp(data.dimord, 'rpt_chan_time')
      data.cov   = data.cov(:,datchanindx,datchanindx);
      data.trial = data.trial(:,datchanindx,:);
    else
      ft_error('unexpected dimord');
    end
    data.label = data.label(datchanindx);
    Nchans     = length(data.label);
  end

  if hasbaseline
    % baseline and active are only available together for resampling purposes,
    % or as a noise covariance for SAM beamforming
    % hence I assume here that there are multiple trials in both
    if isfield(baseline, 'avg')
      baseline.avg   = baseline.avg(datchanindx,:);
      baseline.cov   = baseline.cov(datchanindx,datchanindx);
    else
      baseline.cov   = baseline.cov(:,datchanindx,datchanindx);
      baseline.trial = baseline.trial(:,datchanindx,:);
    end

    % this is required for averaging 2 conditions using prepare_resampled_data
    cfg2 = [];
    cfg2.numcondition = 2;
  end

  % prepare the resampling of the trials, or average the data if multiple trials are present and no resampling is necessary
  if (strcmp(cfg.jackknife, 'yes') || strcmp(cfg.bootstrap, 'yes') || strcmp(cfg.pseudovalue, 'yes') || strcmp(cfg.singletrial, 'yes') || strcmp(cfg.rawtrial, 'yes') || strcmp(cfg.randomization, 'yes')) && ~strcmp(data.dimord, 'rpt_chan_time')
    ft_error('multiple trials required in the data\n');

  elseif strcmp(cfg.permutation, 'yes')
    % compute the average and covariance without resampling
    [dum, avgA, covA, avgB, covB] = prepare_resampled_data(cfg2 , data.trial, data.cov, baseline.trial, baseline.cov);
    % compute the average and covariance with random permutation
    [cfg, avRA, coRA, avRB, coRB] = prepare_resampled_data(cfg, data.trial, data.cov, baseline.trial, baseline.cov);
    % concatenate the different resamplings
    avg = cat(1, reshape(avgA, [1 Nchans Nsamples]), reshape(avgB, [1 Nchans Nsamples]), avRA, avRB);
    Cy  = cat(1, reshape(covA, [1 Nchans Nchans  ]), reshape(covB, [1 Nchans Nchans  ]), coRA, coRB);
    % clear temporary working copies
    clear avgA avgB covA covB
    clear avRA avRB coRA coRB
    % the order of the resamplings should be [avgA avgB randA randB randA randB randA randB ....]
    Nrepetitions = 2*cfg.numpermutation + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    avg = avg(order,:,:);
    Cy  = Cy (order,:,:);

  elseif strcmp(cfg.randomization, 'yes')
    % compute the average and covariance without resampling
    [dum, avgA, covA, avgB, covB] = prepare_resampled_data(cfg2 , data.trial, data.cov, baseline.trial, baseline.cov);
    % compute the average and covariance with random resampling
    [cfg, avRA, coRA, avRB, coRB] = prepare_resampled_data(cfg, data.trial, data.cov, baseline.trial, baseline.cov);
    % concatenate the different resamplings
    avg = cat(1, reshape(avgA, [1 Nchans Nsamples]), reshape(avgB, [1 Nchans Nsamples]), avRA, avRB);
    Cy  = cat(1, reshape(covA, [1 Nchans Nchans  ]), reshape(covB, [1 Nchans Nchans  ]), coRA, coRB);
    % clear temporary working copies
    clear avgA avgB covA covB
    clear avRA avRB coRA coRB
    % the order of the resamplings should be [avgA avgB randA randB randA randB randA randB ....]
    Nrepetitions = 2*cfg.numrandomization + 2;
    order = [1 2 3:2:Nrepetitions 4:2:Nrepetitions];
    avg = avg(order,:,:);
    Cy  = Cy (order,:,:);

  elseif strcmp(cfg.jackknife, 'yes')
    % compute the jackknife repetitions for the average and covariance
    [cfg, avg, Cy] = prepare_resampled_data(cfg, data.trial, data.cov);
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.bootstrap, 'yes')
    % compute the bootstrap repetitions for the average and covariance
    [cfg, avg, Cy] = prepare_resampled_data(cfg, data.trial, data.cov);
    Nrepetitions = cfg.numbootstrap;

  elseif strcmp(cfg.pseudovalue, 'yes')
    % compute the pseudovalue repetitions for the average and covariance
    [cfg, avg, Cy] = prepare_resampled_data(cfg, data.trial, data.cov);
    Nrepetitions = Ntrials+1;

  elseif strcmp(cfg.singletrial, 'yes')
    % The idea is that beamformer uses the average covariance to construct the
    % filter and applies it to the single trial covariance/csd. The problem
    % is that beamformer will use the averaged covariance/csd to estimate the
    % power and not the single trial covariance/csd
    ft_error('this option contains a bug, and is therefore not supported at the moment');
    % average the single-trial covariance matrices
    Cy = mean(data.cov,1);
    % copy the average covariance matrix for every individual trial
    Cy = repmat(Cy, [Ntrials 1 1]);
    % keep the single-trial ERFs, rename them to avg for convenience
    avg = data.trial;
    Nrepetitions = Ntrials;

  elseif strcmp(cfg.rawtrial, 'yes')
    % do not do any resampling, keep the single-trial covariances
    Cy = data.cov;
    % do not do any resampling, keep the single-trial ERFs (rename them to avg for convenience)
    avg = data.trial;
    Nrepetitions = Ntrials;

  elseif Ntrials>1
    % average the single-trial covariance matrices
    Cy  = reshape(mean(data.cov,1), [Nchans Nchans]);
    % compute the average ERF
    avg = shiftdim(mean(data.trial,1),1);
    Nrepetitions = 1;

  elseif Ntrials==1
    % select the average covariance matrix
    Cy  = data.cov;
    % select the average ERF
    avg = data.avg;
    Nrepetitions = 1;
  end

  % reshape so that it also looks like one trial (out of many)
  if Nrepetitions==1
    Cy  = reshape(Cy , [1 Nchans Nchans]);
    avg = reshape(avg, [1 Nchans Nsamples]);
  end

  % get the relevant low level options from the cfg and convert into key-value pairs
  tmpcfg = cfg.(cfg.method);
  % disable console feedback for the low-level function in case of multiple repetitions
  if Nrepetitions > 1
    tmpcfg.feedback = 'none';
  end
  methodopt = ft_cfg2keyval(tmpcfg);

  % construct the low-level options for the leadfield computation as key-value pairs, these are passed to the inverse function and FT_COMPUTE_LEADFIELD
  leadfieldopt = {};
  leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(cfg, 'reducerank'));
  leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(cfg, 'backproject'));
  leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(cfg, 'normalize'));
  leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(cfg, 'normalizeparam'));
  leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(cfg, 'weight'));

  size_avg = [size(avg) 1];
  size_Cy  = [size(Cy) 1];
  if strcmp(cfg.method, 'lcmv')% && ~isfield(sourcemodel, 'filter')
    for i = 1:Nrepetitions
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
      fprintf('scanning repetition %d\n', i);
      dip(i) = ft_inverse_lcmv(sourcemodel, sens, headmodel, squeeze_avg, squeeze_Cy, methodopt{:}, leadfieldopt{:});
    end

    % the following has been disabled since it turns out to be wrong (see
    % bugzilla bug 2395)
    %   elseif 0 && strcmp(cfg.method, 'lcmv')
    %     %don't loop over repetitions (slow), but reshape the input data to obtain single trial timecourses efficiently
    %     %in the presence of filters pre-computed on the average (or whatever)
    %     tmpdat = reshape(permute(avg,[2 3 1]),[size_avg(2) size_avg(3)*size_avg(1)]);
    %     tmpdip = ft_inverse_lcmv(sourcemodel, sens, headmodel, tmpdat, squeeze(mean(Cy,1)), methodopt{:}, leadfieldopt{:});
    %     tmpmom = tmpdip.mom{tmpdip.inside(1)};
    %     sizmom = size(tmpmom);
    %
    %     for i=1:length(tmpdip.inside)
    %       indx = tmpdip.inside(i);
    %       tmpdip.mom{indx} = permute(reshape(tmpdip.mom{indx}, [sizmom(1) size_avg(3) size_avg(1)]), [3 1 2]);
    %     end
    %     try, tmpdip = rmfield(tmpdip, 'pow'); end
    %     try, tmpdip = rmfield(tmpdip, 'cov'); end
    %     try, tmpdip = rmfield(tmpdip, 'noise'); end
    %     for i=1:Nrepetitions
    %       dip(i).pos     = tmpdip.pos;
    %       dip(i).inside  = tmpdip.inside;
    %       dip(i).outside = tmpdip.outside;
    %       dip(i).mom     = cell(1,size(tmpdip.pos,1));
    %       if isfield(tmpdip, 'ori')
    %         dip(i).ori   = cell(1,size(tmpdip.pos,1));
    %       end
    %       dip(i).cov     = cell(1,size(tmpdip.pos,1));
    %       dip(i).pow     = nan(size(tmpdip.pos,1),1);
    %       for ii=1:length(tmpdip.inside)
    %         indx             = tmpdip.inside(ii);
    %         tmpmom           = reshape(tmpdip.mom{indx}(i,:,:),[sizmom(1) size_avg(3)]);
    %         dip(i).mom{indx} = tmpmom;
    %         if isfield(tmpdip, 'ori')
    %           dip(i).ori{indx} = tmpdip.ori{indx};
    %         end
    %
    %         % the following recovers the single trial power and covariance, but
    %         % importantly the latency over which the power is defined is the
    %         % latency of the event-related field in the input and not the
    %         % latency of the covariance window, which can differ from the
    %         % former
    %         dip(i).cov{indx} = (tmpmom*tmpmom')./size_avg(3);
    %         if isempty(cfg.lcmv.powmethod) || strcmp(cfg.lcmv.powmethod, 'trace')
    %           dip(i).pow(indx) = trace(dip(i).cov{indx});
    %         else
    %           [tmpu,tmps,tmpv] = svd(dip(i).cov{indx});
    %           dip(i).pow(indx) = tmps(1);
    %         end
    %       end
    %     end

  elseif strcmp(cfg.method, 'sloreta')
    for i=1:Nrepetitions
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
      fprintf('scanning repetition %d\n', i);
      dip(i) = ft_inverse_sloreta(sourcemodel, sens, headmodel, squeeze_avg, squeeze_Cy, methodopt{:}, leadfieldopt{:});
    end

  elseif strcmp(cfg.method, 'eloreta')
    for i=1:Nrepetitions
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
      fprintf('scanning repetition %d\n', i);
      dip(i) = ft_inverse_eloreta(sourcemodel, sens, headmodel, squeeze_avg, squeeze_Cy, methodopt{:}, leadfieldopt{:});
    end
  elseif strcmp(cfg.method, 'sam')
    % convert time in samples for Event Related SAM
    latency = ft_getopt(methodopt, 'latency_toi');
    toi     = ft_getopt(methodopt, 'toi');
    if isempty(toi) && ~isempty(latency) && ~isequal(latency, 'all')
      toi(1) = nearest(data.time, latency(1));
      toi(2) = nearest(data.time, latency(2));
      methodopt = ft_setopt(methodopt, 'toi', toi);
    end
    for i=1:Nrepetitions
      squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      fprintf('scanning repetition %d\n', i);
      dip(i) = ft_inverse_sam(sourcemodel, sens, headmodel, squeeze_avg, squeeze_Cy, methodopt{:}, leadfieldopt{:});
    end
  elseif strcmp(cfg.method, 'pcc')
    for i=1:Nrepetitions
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
      fprintf('scanning repetition %d\n', i);
      dip(i) = ft_inverse_pcc(sourcemodel, sens, headmodel, squeeze_avg, squeeze_Cy, methodopt{:}, leadfieldopt{:}, 'refdip', cfg.refdip, 'refchan', refchanindx, 'supchan', supchanindx);
    end
  elseif strcmp(cfg.method, 'mne')
    for i=1:Nrepetitions
      fprintf('estimating current density distribution for repetition %d\n', i);
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      if hascovariance
        squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
        dip(i) = ft_inverse_mne(sourcemodel, sens, headmodel, squeeze_avg, methodopt{:}, leadfieldopt{:}, 'noisecov', squeeze_Cy);
      else
        dip(i) = ft_inverse_mne(sourcemodel, sens, headmodel, squeeze_avg, methodopt{:}, leadfieldopt{:});
      end
    end
  elseif strcmp(cfg.method, 'harmony')
    for i=1:Nrepetitions
      fprintf('estimating current density distribution for repetition %d\n', i);
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      if hascovariance
        squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
        dip(i) = ft_inverse_harmony(sourcemodel, sens, headmodel, squeeze_avg, methodopt{:}, leadfieldopt{:}, 'noisecov', squeeze_Cy);
      else
        dip(i) = ft_inverse_harmony(sourcemodel, sens, headmodel, squeeze_avg, methodopt{:}, leadfieldopt{:});
      end
    end
  elseif strcmp(cfg.method, 'rv')
    for i=1:Nrepetitions
      fprintf('estimating residual variance at each source position for repetition %d\n', i);
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      dip(i) = ft_inverse_rv(sourcemodel, sens, headmodel, squeeze_avg,      methodopt{:}, leadfieldopt{:});
    end
  elseif strcmp(cfg.method, 'music')
    for i=1:Nrepetitions
      fprintf('computing multiple signal classification for repetition %d\n', i);
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      if hascovariance
        squeeze_Cy  = reshape(Cy(i,:,:), [size_Cy(2)  size_Cy(3)]);
        dip(i) = ft_inverse_music(sourcemodel, sens, headmodel, squeeze_avg, 'cov', squeeze_Cy, methodopt{:}, leadfieldopt{:});
      else
        dip(i) = ft_inverse_music(sourcemodel, sens, headmodel, squeeze_avg,                    methodopt{:}, leadfieldopt{:});
      end
    end
  elseif strcmp(cfg.method, 'mvl')
    for i=1:Nrepetitions
      fprintf('estimating current density distribution for repetition %d\n', i);
      fns = fieldnames(cfg);
      methodopt = cell(1,length(fns));
      n=1;
      for c=1:length(fns)
        methodopt{n} = fns{c};
        methodopt{n+1} = cfg.(fns{c});
        n=n+2;
      end
      squeeze_avg = reshape(avg(i,:,:),[size_avg(2) size_avg(3)]);
      dip(i) = mvlestimate(sourcemodel, sens, headmodel, squeeze_avg, methodopt{:}, leadfieldopt{:});
    end
  else
    ft_error('method ''%s'' is unsupported for source reconstruction in the time domain', cfg.method);
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % do source reconstruction for component topographies
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif iscomp && any(strcmp(cfg.method, {'rv'}))

  % get the relevant low level options from the cfg and convert into key-value pairs
  tmpcfg = cfg.(cfg.method);
  methodopt = ft_cfg2keyval(tmpcfg);

  % construct the low-level options for the leadfield computation as key-value pairs, these are passed to the inverse function and FT_COMPUTE_LEADFIELD
  leadfieldopt = {};
  leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(cfg, 'reducerank'));
  leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(cfg, 'backproject'));
  leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(cfg, 'normalize'));
  leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(cfg, 'normalizeparam'));
  leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(cfg, 'weight'));

  if strcmp(cfg.method, 'rv')
    for i=cfg.component(:)'
      fprintf('estimating residual variance for component %d\n', i);
      dip(i) = ft_inverse_rv(sourcemodel, sens, headmodel, data.topo(:,i), methodopt{:}, leadfieldopt{:});
    end
  elseif strcmp(cfg.method, 'sloreta')
    for i=cfg.component(:)'
      fprintf('estimating sloreta for component %d\n', i);
      dip(i) = ft_inverse_sloreta(sourcemodel, sens, headmodel, data.topo(:,i), [], methodopt{:}, leadfieldopt{:});
    end
  elseif strcmp(cfg.method, 'eloreta')
    for i=cfg.component(:)'
      fprintf('estimating eloreta for component %d\n', i);
      dip(i) = ft_inverse_eloreta(sourcemodel, sens, headmodel, data.topo(:,i), [], methodopt{:}, leadfieldopt{:});
    end
  else
    ft_error('method ''%s'' is unsupported for source reconstruction in the time domain', cfg.method);
  end

else
  ft_error('the specified method ''%s'' combined with the input data of type ''%s'' is not supported', cfg.method, ft_datatype(data));
end % if freq or timelock or comp data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clean up and collect the results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

source = copyfields(sourcemodel, source, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label', 'cfg'});

if exist('dip', 'var')
  % the fields in the dip structure might be more recent than those in the sourcemodel structure
  source = copyfields(dip, source, {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label', 'cfg'});

  % prevent duplication of these fields when copying the content of dip into source.avg or source.trial
  dip    = removefields(dip,       {'pos', 'tri', 'dim', 'inside', 'leadfield', 'leadfielddimord', 'label', 'cfg'});

  if istrue(cfg.(cfg.method).keepfilter) && isfield(dip(1), 'filter')
    for k=1:numel(dip)
      if isfield(sourcemodel, 'label')
        % pre-computed leadfields or filters were used
        dip(k).label = sourcemodel.label;
      else
        % leadfields and filters were computed on the fly
        dip(k).label = sens.label;
      end
      dip(k).filterdimord = '{pos}_ori_chan';
    end
  end
end

if ~istrue(cfg.keepleadfield)
  % remove the precomputed leadfields from the output source (if present)
  source = removefields(source, {'leadfield' 'leadfielddimord' 'label'});
end

cfg.headmodel = headmodel;

% remove the precomputed leadfields from the cfg, regardless of what keepleadfield is saying
% it should not be kept in cfg, since there it takes up too much space
cfg.sourcemodel = removefields(sourcemodel, {'leadfield' 'leadfielddimord' 'filter' 'filterdimord' 'label'});

if strcmp(cfg.jackknife, 'yes')
  source.method = 'jackknife';
  source.trial = dip;
  source.df = Ntrials;
elseif strcmp(cfg.bootstrap, 'yes')
  source.method = 'bootstrap';
  source.trial = dip;
  source.df = Ntrials;
elseif strcmp(cfg.pseudovalue, 'yes')
  source.method = 'pseudovalue';
  source.trial = dip;
elseif strcmp(cfg.singletrial, 'yes')
  source.method = 'singletrial';
  source.trial = dip;
  source.df = Ntrials;     % is this correct?
elseif strcmp(cfg.rawtrial, 'yes')
  source.method = 'rawtrial';
  source.trial = dip;
  source.df = Ntrials;     % is this correct?
elseif strcmp(cfg.randomization, 'yes')
  % assign the randomized resamplings to the output, keeping the special meaning of trial 1 and 2 in mind
  source.method = 'randomization';
  source.avgA = dip(1);
  source.avgB = dip(2);
  source.trialA = dip(1+2*(1:cfg.numrandomization));
  source.trialB = dip(2+2*(1:cfg.numrandomization));
elseif strcmp(cfg.permutation, 'yes')
  % assign the randomized resamplings to the output, keeping the special meaning of trial 1 and 2 in mind
  source.method = 'permutation';
  source.avgA = dip(1);
  source.avgB = dip(2);
  source.trialA = dip(1+2*(1:cfg.numpermutation));
  source.trialB = dip(2+2*(1:cfg.numpermutation));
elseif (strcmp(cfg.jackknife, 'yes') || strcmp(cfg.bootstrap, 'yes') || strcmp(cfg.pseudovalue, 'yes') || strcmp(cfg.singletrial, 'yes') || strcmp(cfg.rawtrial, 'yes')) && strcmp(cfg.keeptrials, 'yes')
  % keep the source reconstruction for each repeated or resampled trial
  source.trial = dip;
elseif exist('dip', 'var')
  % it looks like beamformer analysis was done on an average input, keep the average source reconstruction
  source.method = 'average';
  source.avg = dip;
else
  % apparently no computations were performed
end

% remember the trialinfo
if (strcmp(cfg.keeptrials, 'yes') || strcmp(cfg.method, 'pcc')) && isfield(data, 'trialinfo')
  source.trialinfo = data.trialinfo;
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   data baseline
ft_postamble provenance source
ft_postamble history    source
ft_postamble savevar    source


text:
function riverplot_set_ribbon_property(ribbons, propertyname, propertyval)
%
% riverplot_set_ribbon_property(ribbons, propertyname, propertyval)
%
% e.g., riverplot_set_ribbon_property(ribbons, 'LineStyle', 'none')

if ~iscell(ribbons) || isempty(ribbons) || all(all(cellfun(@isempty, ribbons)))
    return % nothing to do
end

allribbons = cat(1, ribbons{:});

% all ribbon patches
allpatches = cat(1, allribbons(:).patchh);
%set(allpatches, 'LineStyle', 'none');
set(allpatches, propertyname, propertyval);

%
% for i = 1:length(allribbons)
%     line1h{i} = allribbons(i).line1.h;
%     line2h{i} = allribbons(i).line2.h;
%
% end
%
% line1h = cat(1, line1h{:});
% line2h = cat(1, line2h{:});
%
% set(line1h, 'Visible', 'off');
% set(line2h, 'Visible', 'off');

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function sphere_triangle_quad_test06 ( )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_TEST06 tests SPHERE01_TRIANGLE_QUAD_ICOS2V.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2014
%
%  Author:
%
%    John Burkardt
%
  e_test = [ ...
    0, 0, 0; ...
    1, 0, 0; ...
    0, 1, 0; ...
    0, 0, 1; ...
    2, 0, 0; ...
    0, 2, 2; ...
    2, 2, 2; ...
    0, 2, 4; ...
    0, 0, 6; ...
    1, 2, 4; ...
    2, 4, 2; ...
    6, 2, 0; ...
    0, 0, 8; ...
    6, 0, 4; ...
    4, 6, 2; ...
    2, 4, 8; ...
   16, 0, 0 ]';
  n_mc1 = 1000;
  n_mc2 = 10000;
  n_mc3 = 100000;

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_TRIANGLE_QUAD_TEST06\n' );
  fprintf ( 1, '  SPHERE01_TRIANGLE_QUAD_ICOS2V approximates the\n' );
  fprintf ( 1, '  integral of a function over a spherical triangle on\n' );
  fprintf ( 1, '  the surface of the unit sphere using a vertex rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We do not have an exact result, so we compare each\n' );
  fprintf ( 1, '  estimate to the final one.\n' );
%
%  Choose three points at random to define a spherical triangle.
%
  [ v1, seed ] = sphere01_sample ( 1, seed );
  [ v2, seed ] = sphere01_sample ( 1, seed );
  [ v3, seed ] = sphere01_sample ( 1, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Vertices of random spherical triangle:\n' );
  fprintf ( 1, '\n' );
  r8vec_transpose_print ( 3, v1, '  V1:' );
  r8vec_transpose_print ( 3, v2, '  V2:' );
  r8vec_transpose_print ( 3, v3, '  V3:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  FACTOR   N   RESULT\n' );

  for j = 1 : 17

    e(1:3) = e_test(1:3,j);

    e = polyterm_exponent ( 'SET', e );

    polyterm_exponent ( 'PRINT', e );

%   factor = 2 ^ 11;
    factor = 2 ^ 7;
    [ best, node_num ] = sphere01_triangle_quad_icos2v ( v1, v2, v3, factor, ...
      @polyterm_value_3d );

    factor = 1;
    for factor_log = 0 : 5

      [ result, node_num ] = sphere01_triangle_quad_icos2v ( v1, v2, v3, ...
        factor, @polyterm_value_3d );

      error = abs ( result - best );

      fprintf ( 1, '  %4d  %8d  %16.8g  %10.2e\n', ...
        factor, node_num, result, error );

      factor = factor * 2;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function u=diffBurg(Afunc,u0,nu,dx,T,Nt,bndcond)
%
%    Solves the equation u_t + 0.5(u^2)_x  = \nu A(u)_{xx} by operator
%    splitting.  u0 and a must the size of x.
%
%    Output is a matrix of size lentgth(x)\times (T/dt)+1. length of x
%    must be even!
%
%%
% initial setup
if nargin<7
	bndcond='neumann';
end;
dt=T/Nt;
m=length(u0);
u=zeros(m,Nt+1);
u(:,1)=u0;

%%
% Solving by operator splitting
for i=1:Nt,
	u1=EOBurg(u(:,i),dt,dx,1,bndcond);                % Burgers' equation
	u(:,i+1)=nonlindiff(Afunc,nu,u1,dt,dx,1,bndcond); % Nonlinear diffusion
end;

repo: Aero-matlab


text:
function gamma = angle_shift_deg ( alpha, beta )

%*****************************************************************************80
%
%% ANGLE_SHIFT_DEG shifts angle ALPHA to lie between BETA and BETA+360.
%
%  Discussion:
%
%    The input angle ALPHA is shifted by multiples of 360 to lie
%    between BETA and BETA+360.
%
%    The resulting angle GAMMA has all the same trigonometric function
%    values as ALPHA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the angle to be shifted, in degrees.
%
%    Input, real BETA, defines the lower endpoint of
%    the angle range.
%
%    Output, real GAMMA, the shifted angle.
%
  if ( alpha < beta )
    gamma = beta - mod ( beta - alpha, 360.0 ) + 360.0;
  else
    gamma = beta + mod ( alpha - beta, 360.0 );
  end

  return
end


text:
function CPD = update_CPT(CPD)
% Compute the big CPT for an HHMM F node given internal termprob
% function CPD = update_CPT(CPD)

Qsz = CPD.Qsz;
Qpsz = CPD.Qpsz;

% CPT(Qpsz, Q, Fbelow, Fself)
CPT = zeros(Qpsz, Qsz, 2, 2);
CPT(:,:,1,1) = 1; % if Fbelow=1 (off), then Fself=1 (off)
CPT(:,:,2,:) = CPD.termprob;

CPD = set_fields(CPD, 'CPT', CPT);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [acc, mederr] = compute_vp_acc_mederror(view_filename, img_label_filename)

est_views = importdata(view_filename);
object = importdata(img_label_filename);
data = object.data;
gt_views = data(:,2:4);

% CALCULATE ROTATION MATRIX ANGLES
est_views = est_views/180*pi;
gt_views = gt_views/180*pi;
N = size(gt_views,1);
R_angle_results = [];
for j = 1:N
   R_pred = angle2dcm(est_views(j,1), est_views(j,2), est_views(j,3));
   R_label = angle2dcm(gt_views(j,1), gt_views(j,2), gt_views(j,3));
   R_angle = norm(logm(R_pred'*R_label)) / sqrt(2);
   R_angle_results = [R_angle_results; R_angle];
end

acc = sum(R_angle_results < pi/6) / length(R_angle_results);
mederr = median(R_angle_results);

end

author: ShapeNet

repo: RenderForCNN


text:
function mrk = file_readNIRxMarker(mrkName, varargin)
% file_readNIRxMarker - Read marker in NIRS format from single marker file.
%
% Synopsis:
%   [MRK, FS]= file_readNIRxMarker(MRKNAME,<OPT>)
%
% Arguments:
%   MRKNAME: name of marker file (no extension),
%            relative to BTB.RawDir unless beginning with '/'
%
% Properties:
%   'System': NIRS system used (default 'nirx')
%   'Prefix': string prepended to each marker string (default 'S '; e.g.
%             makes a marker "1" become "S 1")
%   'Path' : directory with raw NIRS data (default BTB.RawDir).
%           If data is file contains an absolute path, path is ignored.

% Returns:
%   MRK: struct array of markers with fields desc (marker descriptor=number
%   from 0 to 15) and time (milliseconds relative to recording onset).
%
% Description:
%   Read all marker information from a Nirx generic data format file.
%   The sampling interval is read from the corresponding head file.
%
% Note: Based on eegfile_readBVmarkers.
% See also: nirs_* nirsfile_*
%
% matthias.treder@tu-berlin.de 2011
% Markus Wenzel 2013 (adapted it to the new toolbox)
% Jan Mehnert February 2014 (ready for pubic BTB toolbox) (jan@mehnert.org)

global BTB
props={ 'Path'           BTB.RawDir   'CHAR'
        'System'         'nirx'       'CHAR'
        'Prefix'         'S'          'CHAR'
        'Verbose'        0            'BOOL'
        'MarkerFormat'   'numeric'    'CHAR(string numeric)'};
if nargin==0,
    mrk= props; return
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);

% Clear opt.Path if mrkName contains absolute path (= starts with '/')
if fileutil_isAbsolutePath(mrkName),
  opt.Path=[];
  isdefault.Path=0;
end

opt_checkProplist(opt, props);
misc_checkType(mrkName, 'CHAR');%|CELL{CHAR}'); % single marker file

%% Scan marker file
fullName= fullfile(opt.Path,mrkName);
mrk = struct();

if strcmp(opt.System,'nirx')
  % Structure of file: (col 1) Timestamp (col 2-9) 8 bits signifying the
  % marker (caution! lowest bit is left [not right], need to flip bits before converting to decimal)
  fid = fopen([fullName '.evt'],'r');
  s= textscan(fid,'%d %s','delimiter','\n');
  pos = s{1}';
  % Remove \t's, flip bits and convert to decimal
  descno = cellfun(@(x)(bin2dec(fliplr(strrep(x,sprintf('\t'),'')))), s{2})';
  desc = str_cprintf([opt.Prefix '%3d'], descno);

  % get sampling frequency
  hdr=file_readNIRxHeader(fullName);
  mrk.time= pos*1000/hdr.fs;

  mrk.event.desc = desc';
end
mrk.time = double(mrk.time);

% New. Avoids errors in mrk_defineClasses. Adapted from file_readBVmarkers.
if strcmp(opt.MarkerFormat, 'numeric'),
  [toe,idx]= bbciutil_markerMappingSposRneg(mrk.event.desc);
  mrk.event.desc= zeros(size(mrk.event.desc));
  mrk.event.desc(idx)= toe;
end

author: bbci


text:
% compute covariance matrix of A on GPU in batches to avoid memory errors
function out = gpuBlockXtX(A)

g = gpuDevice;
bytesGPU = g.AvailableMemory;
infoA = whos('A');
[h, w] = size(A);
bytesPerRow = infoA.bytes/h;
% theoretically, to keep the result and two (A dn A') copies of data
% bytesRequired = bytesPerRow*w + 2*batchSize*bytesPerRow;
batchSize = floor((bytesGPU/bytesPerRow-w)/2);
batchSize = floor(batchSize/8); % just to be on the safe side
% also, interestingly with smaller batches the code runs slightly faster
nBatches = ceil(h/batchSize);

startIdx = 1:batchSize:h;
endIdx = startIdx + batchSize-1;
endIdx = min(endIdx, h);

out = gpuArray(zeros(w, 'like', A));
for iBatch = 1:nBatches
%     fprintf('Batch %d/%d\n', iBatch, nBatches);
    batchA = gpuArray(A(startIdx(iBatch):endIdx(iBatch), :));
    out = out + batchA'*batchA;

    % trying to debug OUT OF MEMORY issues on GPU
%     batchAT = batchA';
%     out = out + batchAT*batchA;
%     clear batchA batchAT;

% MK - this is an annoying line of code, but it helps Matlab to do memory
% management on the GPU properly. Otherwise I get OUT OF MEMORY errors
% maybe it is just my GPU is a bit unstable
    wait(gpuDevice);
end

out = gather(out);

return

%% this is a non-GPU version
batchSize = 2^16;
[h, w] = size(A);
out = zeros(w, 'like', A);
nBatches = ceil(h/batchSize);
startIdx = 1:batchSize:h;
endIdx = startIdx + batchSize-1;
endIdx = min(endIdx, h);

for iBatch = 1:nBatches
    fprintf('Batch %d/%d\n', iBatch, nBatches);
    batchA = A(startIdx(iBatch):endIdx(iBatch), :);
    out = batchA'*batchA + out;
end

author: cortex-lab


text:
addpath matlab/;
vl_compilenn('enableGpu',true);

author: jiangqy


text:
function test_suite = test_strjoin
% tests for cosmo_strjoin
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_strjoin_basics
    aeq=@(a,varargin)assertEqual(cosmo_strjoin(varargin{:}),a);

    aeq('a b cc',{'a','b','cc'});
    aeq('a>#<b>#<cc',{'a','b','cc'}, '>#<');
    aeq(sprintf('a\tb\tcc'),{'a','b','cc'}, '\t');
    aeq('a\b\cc',{'a','b','cc'}, '\\');
    aeq('a+b=cc',{'a','b','cc'}, {'+','='});

function test_strjoin_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                        cosmo_strjoin(varargin{:}),'');

    % firat argument must be cell with strings
    aet('a');
    aet(struct)
    aet(2);
    aet([]);
    aet({2});

    % second argument must be string or cell with strings
    aet({'a','b','c'},2);
    aet({'a','b','c'},struct);
    aet({'a','b','c'},[]);
    aet({'a','b','c'},[]);
    aet({'a','b','c'},{'x'});
    aet({'a','b','c'},{1,2});

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function ImageStack_corrected=ImageUnwrapping(ImageStack,MovedPoint,RawPoint)

%% this function is for correcting the image distortion of the MINI2P images by using 2-D piecewise linear geometric transformation. The distortion is caused by the nonlinearity of the MEMS scanning.

%%
      ImageStack_corrected=ImageStack;
      tform = fitgeotrans(RawPoint,MovedPoint,'polynomial',3); %% generate the affine transformation mask
      Z_size=size(ImageStack,3);
      for j=1:1:Z_size
      ImageStack_corrected(:,:,j) = imwarp(ImageStack(:,:,j),tform,'OutputView',imref2d(size(ImageStack(:,:,j))));
%       ImageStack_corrected_16bit=uint16(ImageStack_corrected(:,:,j));
%       imwrite(ImageStack_corrected_16bit,[Result_folder,'\',FileName,'_wrappiing corrected.tif'],'WriteMode','append');
      disp(['Unwrapping image number ',num2str(j)]);
      end

end

author: kavli-ntnu


text:
folder = 'E:\datasets\cohn-kanade-images\keep-align-kaipeng';
addpath('..');
image_list = get_image_list_in_folder(folder);
list_file = 'E:\datasets\cohn-kanade-images\list.txt';

list_fid = fopen(list_file,'w');
for i=1:length(image_list)
    [a,b,c]=fileparts(image_list{i});
    fprintf(list_fid,'%s %c\r\n', image_list{i}, b(end));
end;
fclose(list_fid);

author: happynear

repo: FaceVerification


text:
function refineGrid(im, success, stats, cnt, id0, id1, searchPoints)

    %Refine with a bigger grid
    xymin = -30;
    xystep = 2;
    xymax = 30;
    %Define a delta distance
    delta = 7; %[pixel]
    [x,y] = meshgrid([xymin:xystep:xymax]);
    x = reshape(x, 1, numel(x));
    y = reshape(y, 1, numel(y));
    xy = [x;y];
    usedPtsIndex = ones(1,length(x));

    %Go through the points
    [H, err] = invpersp(searchPoints(2:3,:), searchPoints(4:5,:));
    xyp = homoTrans(H,[xy;ones(1,length(xy))]);
    plot(xyp(1,:),xyp(2,:),'r+')

author: christianwengert


text:
function i4col_sort2_a_test ( )

%*****************************************************************************80
%
%% I4COL_SORT2_A_TEST tests I4COL_SORT2_A;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  m = 6;
  n = 4;

  b = 0;
  c = 20;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4COL_SORT2_A_TEST\n' );
  fprintf ( 1, '  For a rectangular integer matrix:\n' );
  fprintf ( 1, '  I4COL_SORT2_D sorts the elements of the columns.\n' );

  [ a, seed ] = i4mat_uniform_ab ( m, n, b, c, seed );

  i4mat_print ( m, n, a, '  The matrix:' );

  a = i4col_sort2_a ( m, n, a );

  i4mat_print ( m, n, a, '  The element-sorted column matrix:' );

  return
end


text:
function feats = ecgsqi(signal,qrs,fs)
% This function calculates multiple SQIs for electrocardiogram (ECG) signals.
%
% Input:
%     signal            ECG signal (function only supports vectors)
%     fs                sampling frequency (in Hz)
%
% Output:
%     sqi               Resulting SQI for ABP signal
%    qrs_out            QRS samplestamps
%
%
% --
% ECG classification from single-lead segments using Deep Convolutional Neural
% Networks and Feature-Based Approaches - December 2017
%
% Released under the GNU General Public License
%
% Copyright (C) 2017  Fernando Andreotti, Oliver Carr
% University of Oxford, Insitute of Biomedical Engineering, CIBIM Lab - Oxford 2017
% fernando.andreotti@eng.ox.ac.uk
%
%
% For more information visit: https://github.com/fernandoandreotti/cinc-challenge2017
%
% Referencing this work
%
% Andreotti, F., Carr, O., Pimentel, M.A.F., Mahdi, A., & De Vos, M. (2017).
% Comparing Feature Based Classifiers and Convolutional Neural Networks to Detect
% Arrhythmia from Short Segments of ECG. In Computing in Cardiology. Rennes (France).
%
% Last updated : December 2017
%
% This program 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.
%
% This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.

%% Parameters
WIN_ACCEPT = 0.10;    % acceptance interval for FP in QRS detection (in s)
WIN_QRS = 0.1;         % window used to delimit QRS complexes (in s)

% Temporal SQIs
% feat_flat = flatsqi(signal);  % flatline detection
feat_stdsqi  = stdsqi(signal);   % standard deviation
feat_ksqi  = ksqi(signal);   % standard deviation
feat_ssqi  = ssqi(signal);   % standard deviation

% Frequency bands
feat_psqi = psqi(signal,fs,[15 45],[0 100]);

% Detection-based SQIs
combs = nchoosek(1:length(qrs),2);
combs = num2cell(combs,2);
feat_bsqi = cellfun(@(x) bsqi(qrs{x(1)},qrs{x(2)},WIN_ACCEPT,fs),combs);
feat_rsqi = arrayfun(@(x) rsqi(qrs{x},fs,0.96),1:length(qrs));
feat_csqi = arrayfun(@(x) csqi(signal,qrs{x},fs,WIN_QRS),1:length(qrs));
feat_xsqi = arrayfun(@(x) xsqi(signal,qrs{x},fs,WIN_QRS),1:length(qrs));

feats = [feat_stdsqi, feat_ksqi, feat_ssqi, ...
    feat_psqi, feat_bsqi', feat_rsqi, feat_csqi, feat_xsqi];

feats(isnan(feats)) = 0; % making sure there are no NaNs

author: fernandoandreotti


text:
function [compareImg, rng, ttltxt] = rxCompare(A,B,method);
%
% [compareImg, rng, ttltxt] = rxCompare(A,B,method);
%
% Compare two images (generally an interpolated
% slice and reference slice) and generate an
% image reflecting the similarity between the
% two, using the selected method (an integer).
%
% Methods:
%
%
%            1) Subtract
%            2) Overlay (red / green)
%            3) Checkerboard
%			 4) threshold A and overlay on top of B
%				(like a mean map on an inplane in mrVista)
%
%
%
% 02/05 ras.
compareImg = [];
rng = [];

switch method

    case 1, % color overlay (red/blue)
        compareImg = B;
        compareImg(:,:,2) = A;
        compareImg(:,:,3) = A;
        rng = [0 1];
        ttltxt = 'Blue-Green: Interpolated; Red: Reference';

    case 2, % subtract
        compareImg = A-B;
        rng = [-1 1];
        ttltxt = 'Interpolated Minus Reference';

    case 3, % make checkerboard
        compareImg = regMosaic(A,B,10);
        rng = [min(compareImg(:)) max(compareImg(:))];
        ttltxt = 'Interpolated / Reference Mosaic';

	case 4, % thresholded vol as overlay
% 		compareImg = repmat(B, [1 1 3]);
% 		overlay = ind2rgb(A.*255, hot(256));
% 		thresh = min(A(:)) + 0.1*max(A(:));
% 		for z = 1:3
% 			tmp = compareImg(:,:,z);
% 			tmp(A > thresh) = overlay(A > thresh);
% 			compareImg(:,:,z) = tmp;
% 		end

		[A B] = swap(A, B);

		thresh = min(A(:)) + 0.1*max(A(:));
		C = B .* 127; % map underlay to grayscale range, 0-127
		C(A > thresh) = [128 + A(A > thresh).*128]; % overlay range 128-256
		compareImg = ind2rgb(round(C), [gray(128); hot(128)]);

        rng = [0 1];
        ttltxt = 'Reference Slice, w/ Interp Volume Overlay';

end

% now check if stats are requested,
% and if so display them:
hstats = findobj('Tag','quantifyMenu');

if isempty(hstats)
    showStats = 0;
else
    showStats = isequal(get(hstats(end),'Checked'),'on');
end

if showStats
    % (I only do this now, so the code could
    % conceivably be used outside mrRx):
    cfig = findobj('Tag','rxControlFig');
    rx = get(cfig,'UserData');

    % test for valid text fields
    if isfield(rx.ui,'compareStats') && any(A(:)>0) && any(B(:)>0)
        if ishandle(rx.ui.compareStats.corrcoefVal)
            % compute correlation coefficient
            R = corrcoef(A(:),B(:));
            set(rx.ui.compareStats.corrcoefVal,'String',num2str(R(2)));

            % compute Root Mean-Squared Error
            RMSE = sqrt(mse(A,B));
            set(rx.ui.compareStats.rmseVal,'String',num2str(RMSE));
        else
            rx.ui = rmfield(rx.ui,'compareStats');
        end
    end
end

return

function [B A] = swap(A, B);
% switch variable values -- (shortest function ever!).
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function zCart=polarAzEquidistProj2Cart(pAzPts,latLonRef,a,f)
%%POLARAZEQUIDISTPROJ2CART Given points as polar coordinates sof an
%       azimuthal equidistant projection, convert the points to Cartesian
%       coordinates. If a third height coordinate is omitted, the points
%       are assumed to have zero height.
%
%INPUTS: pAzPts A 2XN set of the [ground distance; heading] points, with the
%              heading given in radians East of North, to convert.
%              Alternatively, if heights are also given, this can be a 3XN
%              set of points with the height being the third dimension.
%    latLonRef A 2X1 [latitude;longitude] reference point in radians about
%              which the projection is taken.
%            a The semi-major axis of the reference ellipsoid (in meters).
%              If this argument is omitted or an empty matrix is passed,
%              the value in Constants.WGS84SemiMajorAxis is used.
%            f The flattening factor of the reference ellipsoid. If this
%              argument is omitted or an empty matrix is passed, the value
%              in Constants.WGS84Flattening is used.
%
%OUTPUTS: zCart The 3XN set of converted points in Cartesian coordinates.
%
%This function just calls polarAzEquidistProj2Ellipse and then ellips2Cart.
%
%January 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(f))
    f=Constants.WGS84Flattening;
end

if(nargin<3||isempty(a))
    a=Constants.WGS84SemiMajorAxis;
end

latLonPts=polarAzEquidistProj2Ellipse(pAzPts,latLonRef,a,f);
zCart=ellips2Cart(latLonPts,a,f);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [S,lam] = cov_shrink(X,w)
% Estimate the covariance matrix for some data using analytical shrinkage.
% Sigma = cov_shrink(X,w)
%
% This gives robust estimates even when the ratio of number of observations to number of variables is very low. [1]
%
% X is a matrix where each row is an observation, and each column is a
% variable. The mean is removed from each column before calculating the
% result. The estimate is unbiased, i.e. cov_shrink normalizes by N-1.
%
% w is an optional [Nx1] vector of nonnegative weights.
%
% See also:
%   cov, mean
%
% Dependencies:
%   bsxfun
%
% References:
%   [1] J. Schaefer and K. Strimmer,
%       "A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics."
%       Statist. Appl. Genet. Mol. Biol. 4:32.
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-03

% Copyright (C) Christian Kothe, SCCN, 2010, christian@sccn.ucsd.edu
%
% This program 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 2 of the
% License, or (at your option) any later version.
%
% This program 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 this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

mem_quota = 0.1;        % use at most this fraction of the free memory for temp data
block_increment = 1.5;  % consider incrementally larger numbers of blocks (incremented by this factor)

[n,p] = size(X);
if nargin < 2
    w = ones(n,1)/n;
else
    w = w/sum(w);
end

S = var(X,w)*n/(n-1);

if p ~= 1
    SD = sqrt(S);
    scale = diag(SD);
    % center & standardize
    X = bsxfun(@minus,X,sum(bsxfun(@times,X,w)));
    X = bsxfun(@times,X,1./SD);
    % the following computation can be quite memory intensive; do it in successively smaller blocks (if it fails)...
    VR = zeros(p);
    R = zeros(p);
    % calc candidate data partitionings
    tryblocks = ceil(block_increment.^(0:p));
    % calc according temporary data sizes in bytes
    trysizes = size(X,1)*(ceil(p./tryblocks)).^2*8;
    % try only those that would fit within the free-memory quota...
    try
        bean = java.lang.management.ManagementFactory.getOperatingSystemMXBean();
        mem_free = bean.getFreePhysicalMemorySize();
    catch
        mem_free = 1e10; % make a fairly conservative assumption
    end
    for blocks = tryblocks(trysizes <= mem_free*mem_quota)
        % 2d block length
        blen = ceil(p/blocks);
        % pad X with zeros
        X(:,end+1:blen*blocks) = 0;
        try
            for bx = 0:blocks-1
                for by = 0:blocks-1
                    ix = 1+bx*blen : blen+bx*blen;
                    iy = 1+by*blen : blen+by*blen;
                    % calc a block of the sample covariance tensor
                    B = bsxfun(@times,permute(X(:,ix),[2 3 1]),permute(X(:,iy),[3 2 1]));
                    % calc a block of the per-element weighted variance of the correlation matrix
                    VR(ix,iy) = var(B,w,3)*n^2/((n-1)^3);
                    % and calc a block of the mean of the correlation matrix
                    R(ix,iy) = sum(bsxfun(@times,permute(w,[3 2 1]),B),3)*n/(n-1);
                end
            end
            % remove zeros
            VR = VR(1:p,1:p);
            R = R(1:p,1:p);
            % success...
            break;
         catch
            % out of memory, start again
         end
    end
    % calc shrinkage parameter
    lam = max(0,min(1,sum(sum(tril(VR,-1)))/sum(sum(tril(R,-1).^2))));
    % calc correlation matrix with shrinkage
    R = (1-lam)*R;
    % Round-off error compensation
    R(logical(eye(p))) = 1;
    % scale back to covariance matrix
    S = scale*R*scale;
else
    lam = 0;
end

repo: REST


text:
% [INPUT]
% am = A binary n-by-n matrix {0;1} representing the adjcency matrix.
%
% [OUTPUT]
% bc = A row vector of floats [0,Inf) of length n representing the betweenness centrality values.
% cc = A row vector of floats [0,Inf) of length n representing the closeness centrality values.
% dc = A row vector of floats [0,Inf) of length n representing the degree centrality values.
% ec = A row vector of floats [0,Inf) of length n representing the eigenvector centrality values.
% kc = A row vector of floats [0,Inf) of length n representing the Katz centrality values.
% clc = A row vector of floats [0,Inf) of length n representing the clustering coefficient.
% deg = A row vector of floats [0,Inf) of length n representing the degrees.
% deg_in = A row vector of floats [0,Inf) of length n representing the in-degrees.
% deg_out = A row vector of floats [0,Inf) of length n representing the out-degrees.

function [bc,cc,dc,ec,kc,clc,deg,deg_in,deg_out] = network_centralities(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('am',@(x)validateattributes(x,{'double'},{'real' 'finite' 'binary' '2d' 'square' 'nonempty'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    am = validate_input(ipr.am);

    nargoutchk(9,9);

    [bc,cc,dc,ec,kc,clc,deg,deg_in,deg_out] = network_centralities_internal(am);

end

function [bc,cc,dc,ec,kc,clc,deg,deg_in,deg_out] = network_centralities_internal(am)

    am_len = size(am,1);

    bc = betweenness_centrality(am,am_len);
    [deg,deg_in,deg_out,dc] = degree_centrality(am,am_len);
    cc = closeness_centrality(am,am_len);
    ec = eigenvector_centrality(am);
    kc = katz_centrality(am,am_len);
    clc = clustering_coefficient(am,am_len,deg);

end

function bc = betweenness_centrality(am,am_len)

    bc = zeros(1,am_len);

    for i = 1:am_len
        depth = 0;
        nsp = accumarray([1 i],1,[1 am_len]);
        bfs = false(250,am_len);
        fringe = am(i,:);

        while ((nnz(fringe) > 0) && (depth <= 250))
            depth = depth + 1;
            nsp = nsp + fringe;
            bfs(depth,:) = logical(fringe);
            fringe = (fringe * am) .* ~nsp;
        end

        [rows,cols,v] = find(nsp);
        v = 1 ./ v;

        nsp_inv = accumarray([rows.' cols.'],v,[1 am_len]);

        bcu = ones(1,am_len);

        for depth = depth:-1:2
            w = (bfs(depth,:) .* nsp_inv) .* bcu;
            bcu = bcu + ((am * w.').' .* bfs(depth-1,:)) .* nsp;
        end

        bc = bc + sum(bcu,1);
    end

    bc = bc - am_len;
    bc = (bc .* 2) ./ ((am_len - 1) * (am_len - 2));

end

function cc = closeness_centrality(am,am_len)

    cc = zeros(1,am_len);

    for i = 1:am_len
        paths = dijkstra_shortest_paths(am,am_len,i);
        paths_sum = sum(paths(~isinf(paths)));

        if (paths_sum ~= 0)
            cc(i) = 1 / paths_sum;
        end
    end

    cc = cc .* (am_len - 1);

end

function clc = clustering_coefficient(am,am_len,deg)

    if (issymmetric(am))
        f = 2;
    else
        f = 1;
    end

    clc = zeros(am_len,1);

    for i = 1:am_len
        degree = deg(i);

        if ((degree == 0) || (degree == 1))
            continue;
        end

        k_neighbors = find(am(i,:) ~= 0);
        k_subgraph = am(k_neighbors,k_neighbors);

        if (issymmetric(k_subgraph))
            k_subgraph_trace = trace(k_subgraph);

            if (k_subgraph_trace == 0)
                edges = sum(sum(k_subgraph)) / 2;
            else
                edges = ((sum(sum(k_subgraph)) - k_subgraph_trace) / 2) + k_subgraph_trace;
            end
        else
            edges = sum(sum(k_subgraph));
        end

        clc(i) = (f * edges) / (degree * (degree - 1));
    end

    clc = clc.';

end

function [deg,deg_in,deg_out,dc] = degree_centrality(am,am_len)

    deg_in = sum(am);
    deg_out = sum(am.');

    if (issymmetric(am))
        deg = deg_in + diag(am).';
    else
        deg = deg_in + deg_out;
    end

    dc = deg ./ (am_len - 1);

end

function ec = eigenvector_centrality(am)

    [eigen_vector,eigen_values] = eig(am);
    [~,indices] = max(diag(eigen_values));

    ec = abs(eigen_vector(:,indices)).';
    ec = ec ./ sum(ec);

end

function kc = katz_centrality(am,am_len)

    kc = linsolve(eye(am_len) - (am .* 0.1),ones(am_len,1));
    kc = kc.' ./ (sign(sum(kc)) * norm(kc,'fro'));

end

function paths = dijkstra_shortest_paths(am,am_len,node)

    paths = Inf(1,am_len);
    paths(node) = 0;

    s = 1:am_len;

    while (~isempty(s))
        [~,idx] = min(paths(s));
        s_min = s(idx);

        for i = 1:length(s)
            s_i = s(i);

            offset = am(s_min,s_i);
            offset_sum = offset + paths(s_min);

            if ((offset > 0) && (paths(s_i) > offset_sum))
                paths(s_i) = offset_sum;
            end
        end

        s = setdiff(s,s_min);
    end

end

function am = validate_input(am)

    amv = am(:);

    if (numel(amv) < 4)
        error('The value of ''am'' is invalid. Expected input to be a square matrix with a minimum size of 2x2.');
    end

end

author: TommasoBelluzzo

repo: SystemicRisk