ft_componentanalysis.m

function [comp] = ft_componentanalysis(cfg, data)

% FT_COMPONENTANALYSIS performs independent component analysis or other
% spatio-temporal decompositions of EEG or MEG data. This function computes
% the topography and timecourses of the components. The output of this
% function can be further analyzed with FT_TIMELOCKANALYSIS or
% FT_FREQANALYSIS.
%
% Use as
%   [comp] = ft_componentanalysis(cfg, data)
% where cfg is a configuration structure and the input data is obtained from
% FT_PREPROCESSING or from FT_TIMELOCKANALYSIS.
%
% The configuration should contain
%   cfg.method       = 'runica', 'fastica', 'binica', 'pca', 'svd', 'jader',
%                      'varimax', 'dss', 'cca', 'sobi', 'white' or 'csp'
%                      (default = 'runica')
%   cfg.channel      = cell-array with channel selection (default = 'all'),
%                      see FT_CHANNELSELECTION for details
%   cfg.split        = cell-array of channel types between which covariance
%                      is split, it can also be 'all' or 'no' (default = 'no')
%   cfg.trials       = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.numcomponent = 'all' or number (default = 'all')
%   cfg.demean       = 'no' or 'yes', whether to demean the input data (default = 'yes')
%   cfg.updatesens   = 'no' or 'yes' (default = 'yes')
%   cfg.feedback     = 'no', 'text', 'textbar', 'gui' (default = 'text')
%
% The runica method supports the following method-specific options. The
% values that these options can take can be found with HELP RUNICA.
%   cfg.runica.extended
%   cfg.runica.pca
%   cfg.runica.sphering
%   cfg.runica.weights
%   cfg.runica.lrate
%   cfg.runica.block
%   cfg.runica.anneal
%   cfg.runica.annealdeg
%   cfg.runica.stop
%   cfg.runica.maxsteps
%   cfg.runica.bias
%   cfg.runica.momentum
%   cfg.runica.specgram
%   cfg.runica.posact
%   cfg.runica.verbose
%   cfg.runica.logfile
%   cfg.runica.interput
%
% The fastica method supports the following method-specific options. The
% values that these options can take can be found with HELP FASTICA.
%   cfg.fastica.approach
%   cfg.fastica.numOfIC
%   cfg.fastica.g
%   cfg.fastica.finetune
%   cfg.fastica.a1
%   cfg.fastica.a2
%   cfg.fastica.mu
%   cfg.fastica.stabilization
%   cfg.fastica.epsilon
%   cfg.fastica.maxNumIterations
%   cfg.fastica.maxFinetune
%   cfg.fastica.sampleSize
%   cfg.fastica.initGuess
%   cfg.fastica.verbose
%   cfg.fastica.displayMode
%   cfg.fastica.displayInterval
%   cfg.fastica.firstEig
%   cfg.fastica.lastEig
%   cfg.fastica.interactivePCA
%   cfg.fastica.pcaE
%   cfg.fastica.pcaD
%   cfg.fastica.whiteSig
%   cfg.fastica.whiteMat
%   cfg.fastica.dewhiteMat
%   cfg.fastica.only
%
% The binica method supports the following method-specific options. The
% values that these options can take can be found with HELP BINICA.
%   cfg.binica.extended
%   cfg.binica.pca
%   cfg.binica.sphering
%   cfg.binica.lrate
%   cfg.binica.blocksize
%   cfg.binica.maxsteps
%   cfg.binica.stop
%   cfg.binica.weightsin
%   cfg.binica.verbose
%   cfg.binica.filenum
%   cfg.binica.posact
%   cfg.binica.annealstep
%   cfg.binica.annealdeg
%   cfg.binica.bias
%   cfg.binica.momentum
%
% The dss method requires the following method-specific option and supports
% a whole lot of other options. The values that these options can take can
% be found with HELP DSS_CREATE_STATE.
%   cfg.dss.denf.function
%   cfg.dss.denf.params
%
% The sobi method supports the following method-specific options. The
% values that these options can take can be found with HELP SOBI.
%   cfg.sobi.n_sources
%   cfg.sobi.p_correlations
%
% The csp method implements the common-spatial patterns method. For CSP, the
% following specific options can be defined:
%   cfg.csp.classlabels = vector that assigns a trial to class 1 or 2.
%   cfg.csp.numfilters  = the number of spatial filters to use (default: 6).
%
% The icasso method implements icasso. It runs fastica a specified number of
% times, and provides information about the stability of the components found
% The following specific options can be defined, see ICASSOEST:
%   cfg.icasso.mode
%   cfg.icasso.Niter
%
% Instead of specifying a component analysis method, you can also specify
% a previously computed unmixing matrix, which will be used to estimate the
% component timecourses in this data. This requires
%   cfg.unmixing     = NxN unmixing matrix
%   cfg.topolabel    = Nx1 cell-array with the channel labels
%
% You may specify a particular seed for random numbers called by
% rand/randn/randi, or the random state used by a previous call to this
% function to replicate results. For example:
%   cfg.randomseed   = integer seed value of user's choice
%   cfg.randomseed   = comp.cfg.callinfo.randomseed (from previous call)
%
% 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_TOPOPLOTIC, FT_REJECTCOMPONENT, FASTICA, RUNICA, BINICA, SVD,
% JADER, VARIMAX, DSS, CCA, SOBI, ICASSO

% Copyright (C) 2003-2012, 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$

% undocumented cfg options:
%   cfg.cellmode = string, 'no' or 'yes', allows to run in cell-mode, i.e.
%     no concatenation across trials is needed. This is based on experimental
%     code and only supported for 'dss', 'fastica' and 'bsscca' as methods.

% 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
ft_preamble provenance data
ft_preamble randomseed

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

% check if the input data is valid for this function
istimelock = ft_datatype(data, 'timelock');
data = ft_checkdata(data, 'datatype', 'raw', 'feedback', 'yes');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels', 'trial'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'forbidden',  {'detrend'});
cfg = ft_checkconfig(cfg, 'renamed',    {'blc', 'demean'});
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'predetermined mixing matrix', 'predetermined unmixing matrix'});
cfg = ft_checkconfig(cfg, 'deprecated', {'topo'});

% set the defaults
cfg.method          = ft_getopt(cfg, 'method',       'runica');
cfg.demean          = ft_getopt(cfg, 'demean',       'yes');
cfg.trials          = ft_getopt(cfg, 'trials',       'all', 1);
cfg.channel         = ft_getopt(cfg, 'channel',      'all');
cfg.split           = ft_getopt(cfg, 'split',        'no');
cfg.numcomponent    = ft_getopt(cfg, 'numcomponent', 'all');
cfg.normalisesphere = ft_getopt(cfg, 'normalisesphere', 'yes');
cfg.cellmode        = ft_getopt(cfg, 'cellmode',     'no');
cfg.doscale         = ft_getopt(cfg, 'doscale',      'yes');
cfg.updatesens      = ft_getopt(cfg, 'updatesens',   'yes');
cfg.feedback        = ft_getopt(cfg, 'feedback',     'text');

% select channels, has to be done prior to handling of previous (un)mixing matrix
cfg.channel = ft_channelselection(cfg.channel, data.label);

if istrue(cfg.cellmode)
  ft_hastoolbox('cellfunction', 1);
end

if isfield(cfg, 'topo') && isfield(cfg, 'topolabel')
  ft_warning(['Specifying cfg.topo (= mixing matrix) to determine component '...
    'timecourses in specified data is deprecated; please specify an '...
    'unmixing matrix instead with cfg.unmixing. '...
    'Using cfg.unmixing=pinv(cfg.topo) for now to reproduce old behavior.']);
  
  cfg.unmixing = pinv(cfg.topo);
  cfg = rmfield(cfg, 'topo');
end

if isfield(cfg, 'unmixing') && isfield(cfg, 'topolabel')
  % use the previously determined unmixing matrix on this dataset
  
  % test whether all required channels are present in the data
  [datsel, toposel] = match_str(cfg.channel, cfg.topolabel);
  if length(toposel)~=length(cfg.topolabel)
    ft_error('not all channels that are required for the unmixing are present in the data');
  end
  
  % ensure that all data channels not used in the unmixing should be removed from the channel selection
  tmpchan = match_str(cfg.channel, cfg.topolabel);
  cfg.channel = cfg.channel(tmpchan);
  
  % update some settings where there is no further choice to be made by the user
  cfg.numcomponent = 'all';
  cfg.method       = 'predetermined unmixing matrix';
end

% add the options for the specified methods to the configuration, only if needed
switch cfg.method
  case 'icasso'
    cfg.icasso        = ft_getopt(cfg,        'icasso', []);
    cfg.icasso.mode   = ft_getopt(cfg.icasso, 'mode',   'both');
    cfg.icasso.Niter  = ft_getopt(cfg.icasso, 'Niter',  15);
    cfg.icasso.method = ft_getopt(cfg.icasso, 'method', 'fastica');
    
    cfg.fastica       = ft_getopt(cfg, 'fastica', []);
  case 'fastica'
    % additional options, see FASTICA for details
    cfg.fastica = ft_getopt(cfg, 'fastica', []);
  case 'runica'
    % additional options, see RUNICA for details
    cfg.runica       = ft_getopt(cfg,        'runica',  []);
    cfg.runica.lrate = ft_getopt(cfg.runica, 'lrate',   0.001);
  case 'binica'
    % additional options, see BINICA for details
    cfg.binica       = ft_getopt(cfg,        'binica',  []);
    cfg.binica.lrate = ft_getopt(cfg.binica, 'lrate',   0.001);
  case 'dss'
    % additional options, see DSS for details
    cfg.dss               = ft_getopt(cfg,          'dss',      []);
    cfg.dss.denf          = ft_getopt(cfg.dss,      'denf',     []);
    cfg.dss.denf.function = ft_getopt(cfg.dss.denf, 'function', 'denoise_fica_tanh');
    cfg.dss.denf.params   = ft_getopt(cfg.dss.denf, 'params',   []);
    cfg.dss.preprocf      = ft_getopt(cfg.dss,      'preprocf', []);
    cfg.dss.preprocf.function = ft_getopt(cfg.dss.preprocf, 'function', 'pre_sphere');
    cfg.dss.preprocf.params   = ft_getopt(cfg.dss.preprocf, 'params', []);
  case 'csp'
    % additional options, see CSP for details
    cfg.csp = ft_getopt(cfg, 'csp', []);
    cfg.csp.numfilters = ft_getopt(cfg.csp, 'numfilters', 6);
    cfg.csp.classlabels = ft_getopt(cfg.csp, 'classlabels');
  case 'bsscca'
    % additional options, see BSSCCA for details
    cfg.bsscca           = ft_getopt(cfg,        'bsscca', []);
    cfg.bsscca.refdelay  = ft_getopt(cfg.bsscca, 'refdelay', 1);
    cfg.bsscca.chandelay = ft_getopt(cfg.bsscca, 'chandelay', 0);
    if strcmp(cfg.cellmode, 'no')
      ft_error('cfg.mehod = ''bsscca'' requires cfg.cellmode = ''yes''');
    end
  otherwise
    % do nothing
end

% select trials of interest
tmpcfg = keepfields(cfg, {'trials', 'channel', 'tolerance', 'showcallinfo', 'trackcallinfo', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo'});
data   = ft_selectdata(tmpcfg, data);
% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

% deal with different chantypes if requested
if isequal(cfg.split, 'no')
  chantype = {};
elseif isequal(cfg.split, 'all')
  chantype = unique(ft_chantype(data.label));
else
  chantype = cfg.split;
end

if numel(chantype)>0
  % recurse per specified chantype
  tmpdata = cell(1, numel(chantype));
  for k = 1:numel(chantype)
    tmpcfg         = cfg;
    tmpcfg.channel = data.label(ft_chantype(data.label, lower(chantype{k})));
    tmpcfg.split   = 'no';
    tmpcfg.chantype = lower(chantype{k}); % makes the output labels unique, to allow appending later on
    tmpdata{1,k}   = ft_componentanalysis(tmpcfg, data);
  end
  comp = ft_appenddata([], tmpdata{:});
  return;
else
  %
end

Ntrials  = length(data.trial);
Nchans   = length(data.label);
if Nchans==0
  ft_error('no channels were selected');
end

% default is to compute just as many components as there are channels in the data
if strcmp(cfg.numcomponent, 'all')
  defaultNumCompsUsed = true(1);
  cfg.numcomponent = length(data.label);
else
  defaultNumCompsUsed = false(1);
end

% determine the size of each trial, they can be variable length
Nsamples = zeros(1,Ntrials);
for trial=1:Ntrials
  Nsamples(trial) = size(data.trial{trial},2);
end

if strcmp(cfg.demean, 'yes')
  % optionally perform baseline correction on each trial
  ft_info('baseline correcting data \n');
  for trial=1:Ntrials
    data.trial{trial} = ft_preproc_baselinecorrect(data.trial{trial});
  end
end

if strcmp(cfg.doscale, 'yes')
  % determine the scaling of the data, scale it to approximately unity
  % this will improve the performance of some methods, esp. fastica
  trlidx = 1;
  tmp    = data.trial{trlidx};
  while all(isnan(tmp(:))) % if all data in this trial is NaN
      trlidx  = trlidx + 1; % try next trial
      tmp     = data.trial{trlidx}; % overwrite tmp with next trial
  end
  tmp(~isfinite(tmp)) = 0; % ensure that the scaling is a finite value
  scale = norm((tmp*tmp')./size(tmp,2)); clear tmp;
  scale = sqrt(scale);
  if scale ~= 0
    ft_info('scaling data with 1 over %f\n', scale);
    for trial=1:Ntrials
      data.trial{trial} = data.trial{trial} ./ scale;
    end
  else
    ft_info('no scaling applied, since factor is 0\n');
  end
else
  ft_info('no scaling applied to the data\n');
end

if strcmp(cfg.method, 'sobi')
  
  % concatenate all the data into a 3D matrix respectively 2D (sobi)
  ft_info('concatenating data');
  Nsamples = Nsamples(1);
  dat = zeros(Ntrials, Nchans, Nsamples);
  % all trials should have an equal number of samples
  % and it is assumed that the time axes of all trials are aligned
  for trial=1:Ntrials
    ft_info('.');
    dat(trial,:,:) = data.trial{trial};
  end
  ft_info('\n');
  ft_info('concatenated data matrix size %dx%dx%d\n', size(dat,1), size(dat,2), size(dat,3));
  if Ntrials == 1
    dummy = 0;
    [dat, dummy] = shiftdim(dat);
  else
    dat = shiftdim(dat,1);
  end
  
elseif strcmp(cfg.method, 'csp')
  
  % concatenate the trials into two data matrices, one for each class
  sel1 = find(cfg.csp.classlabels==1);
  sel2 = find(cfg.csp.classlabels==2);
  if min(length(sel1), length(sel2)) == 0
    ft_error('CSP requires class labels!');
  end
  if length(sel1)+length(sel2)~=length(cfg.csp.classlabels)
    ft_warning('not all trials belong to class 1 or 2');
  end
  dat1 = cat(2, data.trial{sel1});
  dat2 = cat(2, data.trial{sel2});
  ft_info('concatenated data matrix size for class 1 is %dx%d\n', size(dat1,1), size(dat1,2));
  ft_info('concatenated data matrix size for class 2 is %dx%d\n', size(dat2,1), size(dat2,2));
  
elseif ~strcmp(cfg.method, 'predetermined unmixing matrix') && strcmp(cfg.cellmode, 'no')
  % concatenate all the data into a 2D matrix unless we already have an
  % unmixing matrix or unless the user request it otherwise
  ft_info('concatenating data');
  
  dat = zeros(Nchans, sum(Nsamples));
  ft_progress('init', cfg.feedback, 'concatenating trials...');
  for trial=1:Ntrials
    ft_progress(trial/Ntrials, 'Concatenating trial %d from %d', trial, Ntrials);
    begsample = sum(Nsamples(1:(trial-1))) + 1;
    endsample = sum(Nsamples(1:trial));
    dat(:,begsample:endsample) = data.trial{trial};
  end
  ft_progress('close')
  ft_info('concatenated data matrix size %dx%d\n', size(dat,1), size(dat,2));
  
  hasdatanans = any(~isfinite(dat(:)));
  if hasdatanans && strcmp(cfg.method, 'dss')
    ft_error('DSS does not work with nans or inf in the data');
  elseif hasdatanans
    ft_info('data contains nan or inf, only using the samples without nan or inf\n');
    finitevals = sum(~isfinite(dat))==0;
    if ~any(finitevals)
      ft_error('no samples remaining');
    else
      dat = dat(:,finitevals);
    end
  end
else
  ft_info('not concatenating data\n');
  dat = data.trial;
  % FIXME cellmode processing is not nan-transparent yet
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform the component analysis
ft_info('starting decomposition using %s\n', cfg.method);
switch cfg.method
  
  case 'icasso'
    % check whether the required low-level toolboxes are installed
    ft_hastoolbox('icasso',  1);
    
    if strcmp(cfg.icasso.method, 'fastica')
      ft_hastoolbox('fastica', 1);
      cfg.fastica.numOfIC = cfg.numcomponent;
      
      optarg = ft_cfg2keyval(cfg.(cfg.icasso.method));
      sR     = icassoEst(cfg.icasso.mode, dat, cfg.icasso.Niter, optarg{:});
    elseif strcmp(cfg.icasso.method, 'dss')
      % recurse into ft_componentanalysis and do some post processing
      tmpcfg        = rmfield(cfg, 'icasso');
      tmpcfg.method = cfg.icasso.method;
      tmpdata       = data;
      
      % initialize the variables to hold the output
      sR.W     = cell(cfg.icasso.Niter, 1);
      sR.A     = cell(cfg.icasso.Niter, 1);
      sR.index = zeros(0,2);
      for k = 1:cfg.icasso.Niter
        tmp = ft_componentanalysis(tmpcfg, tmpdata);
        sR.W{k}  = tmp.unmixing;
        sR.A{k}  = tmp.topo;
        sR.index = cat(1, sR.index, [k*ones(size(tmp.topo,2),1) (1:size(tmp.topo,2))']);
        sR.whiteningMatrix   = tmp.cfg.dss.V;
        sR.dewhiteningMatrix = tmp.cfg.dss.dV;
      end
      sR.signal = dat;
      sR.mode   = cfg.icasso.mode;
      sR.rdim   = size(tmp.topo,2);
    else
      ft_error('only ''fastica'' or ''dss'' is supported as method for icasso');
    end
    
    % do the rest of the icasso related processing
    sR = icassoCluster(sR, 'strategy', 'AL', 'simfcn', 'abscorr', 's2d', 'sim2dis', 'L',cfg.numcomponent);
    sR = icassoProjection(sR, 'cca', 's2d', 'sqrtsim2dis', 'epochs', 75);
    [Iq, mixing, unmixing, dum, index2centrotypes] = icassoResult(sR,cfg.numcomponent);
    
    % this step is done, because in icassoResult mixing is determined to be
    % pinv(unmixing), which yields strange results. Better take it from the
    % individual iterations. NOTE: as a consequence unmixing*mixing is not
    % necessarily identity anymore !!!
    for k = 1:size(mixing,2)
      ix = sR.index(index2centrotypes(k),:);
      mixing(:,k) = sR.A{ix(1)}(:,ix(2));
    end
    
    %[Iq, mixing, unmixing, dat] = icassoShow(sR, 'estimate', 'off', 'L', cfg.numcomponent);
    
    % sort the output according to Iq
    [srt, ix] = sort(-Iq); % account for NaNs
    mixing    = mixing(:, ix);
    unmixing  = unmixing(ix, :);
    
    cfg.icasso.Iq = Iq(ix);
    cfg.icasso.sR = rmfield(sR, 'signal'); % keep the rest of the information
    
  case 'fastica'
    % check whether the required low-level toolboxes are installed
    ft_hastoolbox('fastica', 1);       % see http://www.cis.hut.fi/projects/ica/fastica
    
    if ~defaultNumCompsUsed &&...
        (~isfield(cfg, 'fastica') || ~isfield(cfg.fastica, 'numOfIC'))
      % user has specified cfg.numcomponent and not specified
      % cfg.fastica.numOfIC, so copy cfg.numcomponent over
      cfg.fastica.numOfIC = cfg.numcomponent;
    elseif ~defaultNumCompsUsed &&...
        isfield(cfg, 'fastica') && isfield(cfg.fastica, 'numOfIC')
      % user specified both cfg.numcomponent and cfg.fastica.numOfIC,
      % unsure which one to use
      ft_error('you can specify either cfg.fastica.numOfIC or cfg.numcomponent (they will have the same effect), but not both');
    end
    
    try
      % construct key-value pairs for the optional arguments
      optarg = ft_cfg2keyval(cfg.fastica);
      [mixing, unmixing] = fastica(dat, optarg{:});
    catch
      % the "catch me" syntax is broken on MATLAB74, this fixes it
      me = lasterror;
      % give a hopefully instructive error message
      ft_info(['If you get an out-of-memory in fastica here, and you use fastica 2.5, change fastica.m, line 482: \n' ...
        'from\n' ...
        '  if ~isempty(W)                  %% ORIGINAL VERSION\n' ...
        'to\n' ...
        '  if ~isempty(W) && nargout ~= 2  %% if nargout == 2, we return [A, W], and NOT ICASIG\n']);
      % forward original error
      rethrow(me);
    end
    
  case 'runica'
    % check whether the required low-level toolboxes are installed
    % see http://www.sccn.ucsd.edu/eeglab
    ft_hastoolbox('eeglab', 1);
    
    if ~defaultNumCompsUsed &&...
        (~isfield(cfg, 'runica') || ~isfield(cfg.runica, 'pca'))
      % user has specified cfg.numcomponent and not specified
      % cfg.runica.pca, so copy cfg.numcomponent over
      cfg.runica.pca = cfg.numcomponent;
    elseif ~defaultNumCompsUsed &&...
        isfield(cfg, 'runica') && isfield(cfg.runica, 'pca')
      % user specified both cfg.numcomponent and cfg.runica.pca,
      % unsure which one to use
      ft_error('you can specify either cfg.runica.pca or cfg.numcomponent (they will have the same effect), but not both');
    end
    
    % construct key-value pairs for the optional arguments
    optarg = [ft_cfg2keyval(cfg.runica) {'reset_randomseed' 0}]; % let FieldTrip deal with the random seed handling
    [weights, sphere] = runica(dat, optarg{:});
    
    % scale the sphering matrix to unit norm
    if strcmp(cfg.normalisesphere, 'yes')
      sphere = sphere./norm(sphere);
    end
    
    unmixing = weights*sphere;
    mixing = [];
    
  case 'binica'
    % check whether the required low-level toolboxes are installed
    % see http://www.sccn.ucsd.edu/eeglab
    ft_hastoolbox('eeglab', 1);
    
    if ~defaultNumCompsUsed &&...
        (~isfield(cfg, 'binica') || ~isfield(cfg.binica, 'pca'))
      % user has specified cfg.numcomponent and not specified
      % cfg.binica.pca, so copy cfg.numcomponent over
      cfg.binica.pca = cfg.numcomponent;
    elseif ~defaultNumCompsUsed &&...
        isfield(cfg, 'binica') && isfield(cfg.binica, 'pca')
      % user specified both cfg.numcomponent and cfg.binica.pca,
      % unsure which one to use
      ft_error('you can specify either cfg.binica.pca or cfg.numcomponent (they will have the same effect), but not both');
    end
    
    % construct key-value pairs for the optional arguments
    optarg = ft_cfg2keyval(cfg.binica);
    [weights, sphere] = binica(dat, optarg{:});
    
    % scale the sphering matrix to unit norm
    if strcmp(cfg.normalisesphere, 'yes')
      sphere = sphere./norm(sphere);
    end
    
    unmixing = weights*sphere;
    mixing = [];
    
  case 'jader'
    % check whether the required low-level toolboxes are installed
    % see http://www.sccn.ucsd.edu/eeglab
    ft_hastoolbox('eeglab', 1);
    
    unmixing = jader(dat, cfg.numcomponent);
    mixing = [];
    
  case 'varimax'
    % check whether the required low-level toolboxes are installed
    % see http://www.sccn.ucsd.edu/eeglab
    ft_hastoolbox('eeglab', 1);
    
    unmixing = varimax(dat);
    mixing = [];
    
  case 'cca'
    % check whether the required low-level toolboxes are installed
    % see http://www.sccn.ucsd.edu/eeglab
    ft_hastoolbox('cca', 1);
    
    [y, w] = ccabss(dat);
    unmixing = w';
    mixing = [];
    
  case 'pca'
    % compute data cross-covariance matrix
    if iscell(dat)
      C  = zeros(size(dat{1},1));
      nC = 0;
      for k = 1:numel(dat)
        C  = C + (dat{k}*dat{k}');
        nC = nC + size(dat{k},2);
      end
      C = C./(nC-1);
    else
      C = (dat*dat')./(size(dat,2)-1);
    end

    % eigenvalue decomposition (EVD)
    [E,D] = eig(C);
    
    % sort eigenvectors in descending order of eigenvalues
    d = cat(2,(1:1:Nchans)',diag(D));
    d = sortrows(d, -2);
    
    % return the desired number of principal components
    unmixing = E(:,d(1:cfg.numcomponent,1))';
    mixing = [];
    
    clear C D E d
    
  case 'kpca'
    
    % linear kernel (same as normal covariance)
    %kern = @(X,y) (sum(bsxfun(@times, X, y),2));
    
    % polynomial kernel degree 2
    %kern = @(X,y) (sum(bsxfun(@times, X, y),2).^2);
    
    % RBF kernel
    kern = @(X,y) (exp(-0.5* sqrt(sum(bsxfun(@minus, X, y).^2, 2))));
    
    % compute kernel matrix
    C = zeros(Nchans,Nchans);
    ft_progress('init', cfg.feedback, 'computing kernel matrix...');
    for k = 1:Nchans
      ft_progress(k/Nchans, 'computing kernel matrix %d from %d', k, Nchans);
      C(k,:) = kern(dat, dat(k,:));
    end
    ft_progress('close');
    
    % eigenvalue decomposition (EVD)
    [E,D] = eig(C);
    
    % sort eigenvectors in descending order of eigenvalues
    d = cat(2,(1:1:Nchans)',diag(D));
    d = sortrows(d, -2);
    
    % return the desired number of principal components
    unmixing = E(:,d(1:cfg.numcomponent,1))';
    mixing = [];
    
    clear C D E d
    
  case 'svd'
    % it is more memory efficient to use the (non-scaled) covariance
    if cfg.numcomponent<Nchans
      % compute only the first components
      [u, s, v] = svds(dat*dat', cfg.numcomponent);
    else
      % compute all components
      [u, s, v] = svd(dat*dat', 0);
    end
    clear s v % not needed
    
    unmixing = u';
    mixing = [];
    
  case 'dss'
    % check whether the required low-level toolboxes are installed
    % see http://www.cis.hut.fi/projects/dss
    ft_hastoolbox('dss', 1);
    
    params         = removefields(struct(cfg.dss), {'V' 'dV' 'W' 'indx'});
    params.denf.h  = str2func(cfg.dss.denf.function);
    params.preprocf.h = str2func(cfg.dss.preprocf.function);
    if isfield(cfg.dss, 'wdim') && ~isempty(cfg.dss.wdim)
      params.wdim = cfg.dss.wdim;
    end
    if ~ischar(cfg.numcomponent)
      params.sdim = cfg.numcomponent;
      if isfield(params, 'wdim')
        params.sdim = min(params.sdim, params.wdim);
      end
    end
    

    if isfield(params.denf, 'params') && isfield(params.denf.params, 'artifact')
      % this may require the sampleinfo in the params structure, to keep the sampling bookkeeping correct
      params.denf.params.sampleinfo = data.sampleinfo;
    end
    
    % create the state
    state   = dss_create_state(dat, params);
    if isfield(cfg.dss, 'V') && ~isempty(cfg.dss.V)
      state.V = cfg.dss.V;
      state.Y = cfg.dss.V*dat;
    end
    if isfield(cfg.dss, 'dV') && ~isempty(cfg.dss.dV)
      state.dV = cfg.dss.dV;
    end
    if isfield(cfg.dss, 'W') && ~isempty(cfg.dss.W)
      state.W = cfg.dss.W;
    end
    if isfield(cfg.dss, 'indx') && ~isempty(cfg.dss.indx)
      state.indx = cfg.dss.indx; %may be needed for dss_core_mim
    end
    
    % increase the amount of information that is displayed on screen
    % state.verbose = 3;
    % start the decomposition
    state   = denss(state);  % this is for the DSS toolbox version 1.0
    
    mixing   = state.A;
    unmixing = state.B;
  
    % remember the updated configuration details
    cfg.dss.denf      = state.denf;
    cfg.dss.orthof    = state.orthof;
    cfg.dss.preprocf  = state.preprocf;
    cfg.dss.stopf     = state.stopf;
    cfg.dss.W         = state.W;
    cfg.dss.V         = state.V;
    cfg.dss.dV        = state.dV;
    if isfield(state, 'D'), cfg.dss.D = state.D(1:min([state.sdim size(state.dV)])); end
    cfg.numcomponent  = min([state.sdim size(state.dV)]);
    
  case 'sobi'
    % check whether the required low-level toolboxes are installed
    % see http://www.sccn.ucsd.edu/eeglab
    ft_hastoolbox('eeglab', 1);
    
    % check for additional options, see SOBI for details
    if ~isfield(cfg, 'sobi')
      mixing = sobi(dat, cfg.numcomponent);
    elseif isfield(cfg.sobi, 'n_sources') && isfield(cfg.sobi, 'p_correlations')
      mixing = sobi(dat, cfg.sobi.n_sources, cfg.sobi.p_correlations);
    elseif isfield(cfg.sobi, 'n_sources')
      mixing = sobi(dat,cfg.sobi.n_sources);
    else
      ft_error('unknown options for SOBI component analysis');
    end
    
    unmixing = [];
    
  case 'predetermined unmixing matrix'
    % check which labels from the cfg are identical to those of the data
    % this gives us the rows of cfg.topo (the channels) and of
    % data.trial (also channels) that we are going to use later
    [datsel, chansel] = match_str(data.label, cfg.topolabel);
    
    % ensure 1:1 corresponcence between cfg.topolabel & data.label
    % otherwise we cannot compute the components (if source channels are
    % missing) or will have a problem when projecting it back (because we
    % dont have a marker to say that there are channels in data.label
    % which we did not use and thus can't recover from source-space)
    
    if length(cfg.topolabel)<length(chansel)
      ft_error('cfg.topolabels do not uniquely correspond to data.label, please check')
    end
    if length(data.label)<length(datsel)
      ft_error('cfg.topolabels do not uniquely correspond to data.label, please check')
    end
    
    % reorder the mixing matrix so that the channel order matches the order in the data
    cfg.unmixing  = cfg.unmixing(:,chansel);
    cfg.topolabel = cfg.topolabel(chansel);
    
    unmixing = cfg.unmixing;
    mixing   = [];
    
  case 'white'
    % compute the covariance matrix and an unmixing matrix that makes the data white
    c = dat*dat';
    c = c./(size(dat,2)-1);
    [u, s] = svd(c);
    % split the singular values into half
    for i=1:size(s)
      if (s(i,i)/s(1,1))>(100*eps)
        s(i,i) = 1./sqrt(s(i,i));
      else
        s(i,i) = 0;
      end
    end
    unmixing = s * u';
    mixing   = [];
    
  case 'csp'
    C1 = cov(dat1');
    C2 = cov(dat2');
    unmixing = csp(C1, C2, cfg.csp.numfilters);
    mixing   = [];  % will be computed below
    
  case 'bsscca'
    % this method relies on time shifting of the original data, in much the
    % same way as ft_denoise_tsr. as such it is more natural to represent
    % the data in the cell-array, because the trial-boundaries are clear.
    % if represented in a concatenated array one has to keep track of the
    % trial boundaries
    
    optarg          = ft_cfg2keyval(cfg.bsscca);
    optarg          = cat(2,optarg, {'time', data.time});
    [unmixing, mixing, rho, compdata, time] = bsscca(dat, optarg{:});
    data.trial = mixing*compdata;
    data.time  = time;
    data       = removefields(data, 'sampleinfo');
 
    if size(mixing,1)>numel(data.label)
      for m = 1:(size(mixing,1)-numel(data.label))
        data.label{end+1} = sprintf('refchan%03d',m);
      end
    end
    
    % remember the canonical correlations
    cfg.bsscca.rho = rho;
    
  case 'parafac'
    ft_error('parafac is not supported anymore in ft_componentanalysis');
    
  otherwise
    ft_error('unknown method for component analysis');
end % switch method

% make sure we have both mixing and unmixing matrices
% if not, compute (pseudo-)inverse to go from one to the other
if isempty(unmixing) && ~isempty(mixing)
  if (size(mixing,1)==size(mixing,2))
    unmixing = inv(mixing);
  else
    unmixing = pinv(mixing);
  end
elseif isempty(mixing) && ~isempty(unmixing)
  if (size(unmixing,1)==size(unmixing,2)) && rank(unmixing)==size(unmixing,1)
    mixing = inv(unmixing);
  else
    mixing = pinv(unmixing);
  end
elseif isempty(mixing) && isempty(unmixing)
  % this sanity check is needed to catch convergence problems in fastica
  % see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1519
  ft_error('the component unmixing failed');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% collect the results and construct data structure

comp = keepfields(data, {'time' 'fsample'});

% make sure we don't return more components than were requested
% (some methods respect the maxcomponent parameters, others just always
% return a fixed (i.e., numchans) number of components)
if size(unmixing,1) > cfg.numcomponent
  unmixing(cfg.numcomponent+1:end,:) = [];
end
if size(mixing,2) > cfg.numcomponent
  mixing(:,cfg.numcomponent+1:end) = [];
end

% compute the activations in each trial
if strcmp(cfg.doscale, 'yes')
  for trial=1:Ntrials
    comp.trial{trial} = scale * unmixing * data.trial{trial};
  end
else
  for trial=1:Ntrials
    comp.trial{trial} = unmixing * data.trial{trial};
  end
end

% store mixing/unmixing matrices in structure
comp.topo = mixing;
comp.unmixing = unmixing;

% get the labels
if strcmp(cfg.method, 'predetermined unmixing matrix')
  prefix = 'component';
else
  prefix = cfg.method;
end

st = dbstack;
if numel(st)>1 && isequal(st(2).name, 'ft_componentanalysis')
  % this is a recursive call, as per the cfg.split option, add something
  % extra to the prefix
  chantype = ft_getopt(cfg, 'chantype', '');
  prefix   = [prefix chantype];
end

for k = 1:size(comp.topo,2)
  comp.label{k,1} = sprintf('%s%03d', prefix, k);
end
comp.topolabel = data.label(:);

sensfield = cell(0,1);
if isfield(data, 'grad')
  sensfield{end+1} = 'grad';
end
if isfield(data, 'elec')
  sensfield{end+1} = 'elec';
end
if isfield(data, 'opto')
  sensfield{end+1} = 'opto';
end

% apply the linear projection also to the sensor description
if ~isempty(sensfield)
  if  strcmp(cfg.updatesens, 'yes')
    % construct a montage and apply it to the sensor description
    montage          = [];
    montage.labelold = data.label;
    montage.labelnew = comp.label;
    montage.tra      = unmixing;
    
    for m = 1:numel(sensfield)
      ft_info('also applying the unmixing matrix to the %s structure\n', sensfield{m});
      comp.(sensfield{m}) = ft_apply_montage(data.(sensfield{m}), montage, 'balancename', 'comp', 'keepunused', 'yes');
      
      % The output sensor array cannot simply be interpreted as the input
      % sensor array, hence the type should be removed to allow autodetection
      % See also http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1806
      if isfield(comp.(sensfield{m}), 'type')
        comp.(sensfield{m}) = rmfield(comp.(sensfield{m}), 'type');
      end
    end
  else
    for m = 1:numel(sensfield)
      ft_info('not applying the unmixing matrix to the %s structure\n', sensfield{m});
      % simply copy it over
      comp.(sensfield{m}) = data.(sensfield{m});
    end
  end
end % if sensfield

% copy the sampleinfo into the output
if isfield(data, 'sampleinfo')
  comp.sampleinfo = data.sampleinfo;
end

% copy the trialinfo into the output
if isfield(data, 'trialinfo')
  comp.trialinfo = data.trialinfo;
end

% convert back to input type if necessary
if istimelock
  comp = ft_checkdata(comp, 'datatype', 'timelock+comp');
end

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