functions/popfunc/pop_runica.m

% POP_RUNICA - Run an ICA decomposition of an EEG dataset using RUNICA, 
%                BINICA, or another ICA or other linear decomposition. 
% Usage:
%   >> OUT_EEG = pop_runica( EEG ); % pops-up a data entry window
%   >> OUT_EEG = pop_runica( EEG, 'key', 'val' ); % no pop_up
%
% Graphic interface:
%   "ICA algorithm to use" - [edit box] The ICA algorithm to use for 
%                 ICA decomposition. Command line equivalent: 'icatype'
%   "Commandline options" - [edit box] Command line options to forward
%                 to the ICA algorithm. Command line equivalent: 'options' 
% Inputs:
%   EEG         - input EEG dataset or array of datasets
%
% Optional inputs:
%   'icatype'   - ['runica'|'binica'|'jader'|'fastica'] ICA algorithm 
%                 to use for the ICA decomposition. The nature of any 
%                 differences in the results of these algorithms have 
%                 not been well characterized. {default: BINICA, if
%                 found, else RUNICA}
%   'dataset'   - [integer array] dataset index or indices.
%   'chanind'   - [integer array or cell array] subset of channel indices 
%                 for running the ICA decomposition. Alternatively, you may
%                 also enter channel types here in a cell array.
%   'concatenate' - ['on'|'off'] 'on' concatenate all input datasets 
%                 (assuming there are several). 'off' run ICA independently
%                 on each dataset. Default is 'off'.
%   'concatcond'  - ['on'|'off'] 'on' concatenate conditions for input datasets 
%                 of the same sessions and the same subject. Default is 'off'.
%   'reorder'   - ['on'|'off'] re-order components by variance if that's not
%                 already the case. Default is 'on'.
%   'key','val' - ICA algorithm options (see ICA routine help messages).
% 
% Adding a new algorithm:
% Add the algorithm to the list of algorithms line 366 to 466, for example
%
%    case 'myalgo', [EEG.icaweights] = myalgo( tmpdata, g.options{:} );
%
% where "myalgo" is the name of your algorithm (and Matlab function). 
% tmpdata is the 2-D array containing the EEG data (channels x points) and
% g.options{} contains custom options for your algorithm (there is no
% predetermined format for these options). The output EEG.icaweights is the
% mixing matrix (or inverse of the unmixing matrix).
%
% Note:
% 1) Infomax (runica, binica) is the ICA algorithm we use most. It is based 
%    on Tony Bell's infomax algorithm as implemented for automated use by 
%    Scott Makeig et al. using the natural gradient of Amari et al. It can 
%    also extract sub-Gaussian sources using the (recommended) 'extended' option 
%    of Lee and Girolami. Function RUNICA is the all-Matlab version; function 
%    BINICA calls the (1.5x faster) binary version (a separate download) 
%    translated into C from RUNICA by Sigurd Enghoff.
% 2) JADER calls the JADE algorithm of Jean-Francois Cardoso. This is 
%    included in the EEGLAB toolbox by his permission. See >> help jader
% 3) To run FASTICA, download the fastICA toolbox from its website,
%    http://www.cis.hut.fi/projects/ica/fastica/, and make it available 
%    in your Matlab path. According to its authors, default parameters
%    are not optimal: Try args 'approach', 'sym' to estimate components 
%    in parallel.
%
% Outputs:
%   OUT_EEG = The input EEGLAB dataset with new fields icaweights, icasphere 
%             and icachansind (channel indices). 
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: RUNICA, BINICA, JADER, SOBI, FASTICA.

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad 
% 03-07-02 add the eeglab options -ad
% 03-18-02 add other decomposition options -ad
% 03-19-02 text edition -sm

function [ALLEEG, com] = pop_runica( ALLEEG, varargin )

com = '';
if nargin < 1   
    help pop_runica;
    return;
end

% find available algorithms
% -------------------------
allalgs(1).name = 'runica';
allalgs(1).description = 'Extended Infomax (runica.m; default)';
allalgs(1).options     = '''extended'', 1, ''rndreset'', ''yes''';
allalgs(end+1).name = 'runica';
allalgs(end).description = 'Robust Extended Infomax (runica.m; slow)';
allalgs(end).options     = '''extended'', 1, ''lrate'', 1e-5, ''maxstep'', 2000';
allalgs(end).help        = 'See this <a href="https://sccn.ucsd.edu/wiki/Makoto%27s_useful_EEGLAB_code#How_to_obtain_practically_reproducible_ICA_results_.2809.2F26.2F2022_added.29">reference</a> for ICA conservative parameters ';
allalgs(end+1).name = 'amica';
allalgs(end).description = 'AMICA (slowest; best)';
allalgs(end).options     = '''maxiter'', 2000';
allalgs(end).help        = 'See this <a href="https://github.com/sccn/amica/wiki/AMICA">reference</a> for AMICA';
allalgs(end+1).name = 'picard';
allalgs(end).description = 'Infomax picard.m';
allalgs(end).options     = '''maxiter'', 500, ''mode'', ''standard''';
allalgs(end).help        = 'See this <a href="https://github.com/pierreablin/picard">page</a> for Picard parameters ';
allalgs(end+1).name = 'picard';
allalgs(end).description = 'FastICA picard.m (fastest)';
allalgs(end).options     = '''maxiter'', 500';
allalgs(end).help        = 'See this <a href="https://github.com/pierreablin/picard">page</a> for Picard parameters ';

allalgs(end+1).name = 'sobi'; allalgs(end).description = 'SOBI (sobi.m function)';
allalgs(end+1).name = 'acsobiro'; allalgs(end).description = 'SOBI (acsobiro.m function)';
allalgs(end+1).name = 'pearson_ica'; allalgs(end).description = 'Pearson (pearson_ica.m function)';
allalgs(end+1).name = 'jader';     allalgs(end).description = 'Jade (jader.m function)';
allalgs(end+1).name = 'jadeop';    allalgs(end).description = 'Jade (jadeop.m function)';
allalgs(end+1).name = 'jade_td_p'; allalgs(end).description = 'Jade (jade_td_p.m function)';
allalgs(end+1).name = 'MatlabshibbsR'; allalgs(end).description = 'MatlabshibbsR';
allalgs(end+1).name = 'fastica'; allalgs(end).description = 'FastICA (fastica.m function)';
allalgs(end+1).name = 'tica'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'simbec'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'unica'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'amuse'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'fobi'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'evd'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'evd24'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'sons'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'ng_ol'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'acrsobibpf'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'egld_ica'; allalgs(end).description = allalgs(end).name; % do not use egld_ica => too slow
allalgs(end+1).name = 'eeA'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'tfbss'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'icaML'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'icaMS'; allalgs(end).description = allalgs(end).name;
allalgs(end+1).name = 'binica'; allalgs(end).description = 'Executable runica';
defaultopts = allalgs(1).options;

for index = length(allalgs):-1:1
    if ~exist(allalgs(index).name, 'file')
        allalgs(index) = [];
    end
end

% special AMICA (deprecated use the pop_runamica function)
% -------------
selectamica = 0;
if nargin > 1
    if ischar(varargin{1})
        if strcmpi(varargin{1}, 'selectamica')
            selectamica = 1;
            allalgs = { 'amica' allalgs{:} };
            defaultopts = sprintf('''outdir'', ''%s''', fullfile(pwd, 'amicaout'));
        elseif strcmpi(varargin{1}, 'selectamicaloc')
            selectamica = 1;
            allalgs = { 'amica' allalgs{:} };
            defaultopts = sprintf('''outdir'', ''%s'', ''qsub'', ''off''', fullfile(pwd, 'amicaout'));
        end
    end
end

% popup window parameters
% -----------------------
fig = [];
if nargin < 2 || selectamica
    commandchans = [ 'tmpchans = get(gcbf, ''userdata'');' ...
                     'tmpchans = tmpchans{1};' ...
                     'set(findobj(gcbf, ''tag'', ''chantype''), ''string'', ' ...
                     '       int2str(pop_chansel( tmpchans )));' ...
                     'clear tmpchans;' ];
    commandtype = ['tmptype = get(gcbf, ''userdata'');' ...
                   'tmptype = tmptype{2};' ...
                   'if ~isempty(tmptype),' ...
                   '    [tmps,tmpv, tmpstr] = listdlg2(''PromptString'',''Select type(s)'', ''ListString'', tmptype);' ...
				   '    if tmpv' ...
				   '        set(findobj(''parent'', gcbf, ''tag'', ''chantype''), ''string'', tmpstr);' ...
				   '    end;' ...
                   'else,' ...
                   '    warndlg2(''No channel type'', ''No channel type'');' ...
                   'end;' ...
				   'clear tmps tmpv tmpstr tmptype tmpchans;' ];
    cb_ica = ['tmpalgos = get(gcbf, ''userdata'');' ...
              'tmpalgos = tmpalgos{3};' ...
              'set(findobj(gcbf, ''tag'', ''params''), ''string'',  tmpalgos( get(gcbo, ''value'') ).options );' ... 
              'disp(char(tmpalgos( get(gcbo, ''value'') ).help));' ...
              'clear tmpalgos;' ];
%     [ 'if get(gcbo, ''value'') < 3, ' ...
%                '     set(findobj(gcbf, ''tag'', ''params''), ''string'', ''''''extended'''', 1'');' ...
%                'else,' ...
%                '    tmpStr =  get(gcbo, ''string'');' ....
%                '    tmpAlgo = tmpStr(get(gcbo, ''value''));' ...
%                '    if strcmpi(tmpAlgo, ''picard''),' ...
%                '        set(findobj(gcbf, ''tag'', ''params''), ''string'', ''''''maxiter'''', 500'');' ...
%                '    else,' ...
%                '        set(findobj(gcbf, ''tag'', ''params''), ''string'', '''');' ...
%                '    end;' ...
%                'end;' ];
               
    promptstr    = { { 'style' 'text'       'string' [ 'ICA algorithm to use (click to select)' 10 ] } ...
                     { 'style' 'listbox'    'string' { allalgs.description } 'callback', cb_ica } ...
                     { 'style' 'text'       'string' 'Commandline options (See help messages)' } ...
                     { 'style' 'edit'       'string' defaultopts 'tag' 'params' } ...
                     { 'style'  'checkbox'  'string' 'Reorder components by variance (if that''s not already the case)' 'value' 1 } ...
                     { 'style' 'text'       'string' 'Use only channel type(s) or indices' } ...
                     { 'style' 'edit'       'string' '' 'tag' 'chantype' }  ...
                     { 'style' 'pushbutton' 'string' '... types' 'callback' commandtype } ...
                     { 'style' 'pushbutton' 'string' '... channels' 'callback' commandchans } };
    geometry = { [2 1.5] [2 1.5] [1] [2 1 1 1] };
    geomvert = [ 3 1 1 1];
    if length(ALLEEG) > 1
        cb1 = 'set(findobj(''parent'', gcbf, ''tag'', ''concat2''), ''value'', 0);';
        cb2 = 'set(findobj(''parent'', gcbf, ''tag'', ''concat1''), ''value'', 0);';
        promptstr = { promptstr{:}, ...
                     { 'style' 'text'       'string' 'Concatenate all datasets (check=yes; uncheck=run ICA on each dataset)?' }, ...
                     { 'style' 'checkbox'   'string' '' 'value' 0 'tag' 'concat1' 'callback' cb1 }, ...
                     { 'style' 'text'       'string' 'Concatenate datasets for the same subject and session (check=yes)?' }, ...
                     { 'style' 'checkbox'   'string' '' 'value' 1 'tag' 'concat2' 'callback' cb2 } };
        geometry = { geometry{:} [ 2 0.2 ] [ 2 0.2 ]};
        geomvert = [ geomvert 1 1];
    end
    
    % channel types
    % -------------
    if isfield(ALLEEG(1).chanlocs, 'type')
        tmpchanlocs = ALLEEG(1).chanlocs;
        alltypes = { tmpchanlocs.type };
        indempty = cellfun('isempty', alltypes);
        alltypes(indempty) = '';
        try
            alltypes = unique_bc(alltypes);
        catch
            alltypes = '';
        end
    else
        alltypes = '';
    end
    
    % channel labels
    % --------------
    if ~isempty(ALLEEG(1).chanlocs)
        tmpchanlocs = ALLEEG(1).chanlocs;        
        alllabels = { tmpchanlocs.labels };
    else
        for index = 1:ALLEEG(1).nbchan
            alllabels{index} = int2str(index);
        end
    end
    
    % gui
    % ---
    result       = inputgui( 'geometry', geometry, 'geomvert', geomvert, 'uilist', promptstr, ...
                             'helpcom', 'pophelp(''pop_runica'')', ...
                             'title', 'Run ICA decomposition -- pop_runica()', 'userdata', { alllabels alltypes allalgs } );
    if isempty(result)
        return;
    end
    options = { 'icatype' allalgs(result{1}).name 'dataset' [1:length(ALLEEG)] 'options' eval( [ '{' result{2} '}' ]) 'reorder' fastif(result{3}, 'on', 'off') };
    if ~isempty(result{4})
        if ~isempty(str2num(result{4})), options = { options{:} 'chanind' str2num(result{4}) };
        else                             options = { options{:} 'chanind' parsetxt(result{4}) }; 
        end
    end
    if length(result) > 4
        options = { options{:} 'concatenate' fastif(result{5}, 'on', 'off') };
        options = { options{:} 'concatcond'  fastif(result{6}, 'on', 'off') };
    end
else 
    if mod(length(varargin),2) == 1
        options = { 'icatype' varargin{1:end} };
    else
        options = varargin;
    end
end

% decode input arguments
% ----------------------
[ g, addoptions ] = finputcheck( options, { 'icatype'        'string'  unique({ allalgs.name })   'runica'; ...
                            'dataset'        'integer' []        [1:length(ALLEEG)];
                            'options'        'cell'    []        {};
                            'concatenate'    'string'  { 'on','off' }   'off';
                            'concatcond'     'string'  { 'on','off' }   'off';
                            'reorder'        'string'  { 'on','off' }   'on';
                            'chanind'        { 'cell','integer' } { [] [] }        [];}, ...
                            'pop_runica', 'ignore');
if ischar(g), error(g); end
if ~isempty(addoptions), g.options = { g.options{:} addoptions{:}}; end

% select datasets, create new big dataset if necessary
% ----------------------------------------------------
if length(g.dataset) == 1
    EEG = ALLEEG(g.dataset);
    EEG = eeg_checkset(EEG, 'loaddata');
    if isfield(EEG.etc, 'ic_classification')
        EEG.etc = rmfield(EEG.etc, 'ic_classification');
    end
elseif length(ALLEEG) > 1 && ~strcmpi(g.concatenate, 'on') && ~strcmpi(g.concatcond, 'on')
    [ ALLEEG, com ] = eeg_eval( 'pop_runica', ALLEEG, 'warning', 'off', 'params', ...
           { 'icatype' g.icatype 'options' g.options 'chanind' g.chanind } );
    return;
elseif length(ALLEEG) > 1 && strcmpi(g.concatcond, 'on')
    allsubjects = { ALLEEG.subject };
    allsessions = { ALLEEG.session };
    alltags     = zeros(1,length(allsubjects));
    if any(cellfun('isempty', allsubjects))
        errordlg2( [ 'Aborting: Subject names missing from at least one dataset file.' 10 ...
                     'Subject names must be stored within the datasets. To do so,' 10 ...
                     'use the STUDY > Edit STUDY Info menu and check the box' 10 ...
                     '"Dataset info (condition, group, ...) differs from study info..."' ]);
        return
    end
    dats = {};
    for index = 1:length(allsubjects)
        if ~alltags(index)
            allinds = strmatch(allsubjects{index}, allsubjects, 'exact');
            rmind = [];
            % if we have different sessions they will not be concatenated
            for tmpi = setdiff_bc(allinds,index)'
                if ~isequal(allsessions(index), allsessions(tmpi))
                    rmind = [rmind tmpi];
                end
            end
            allinds = setdiff_bc(allinds, rmind);
            fprintf('Found %d datasets for subject ''%s'' session %d\n', length(allinds), allsubjects{index}, allsessions{index});
            dats = { dats{:} allinds };
            alltags(allinds) = 1;
        end
    end
    fprintf('**************************\nNOW RUNNING ALL DECOMPOSITIONS\n****************************\n');
    eeglab_options;
    if option_parallel
        for index = 1:length(dats)
            TMPEEG{index} = ALLEEG(dats{index});
        end
        parfor index = 1:length(dats)
            TMPEEG{index} = pop_runica(TMPEEG{index}, 'icatype', g.icatype, ...
                'options', g.options, 'chanind', g.chanind, 'concatenate', 'on');
        end
        for index = 1:length(dats)
            ALLEEG(dats{index}) = TMPEEG{index};
            for idat = 1:length(dats{index})
                ALLEEG(dats{index}(idat)).saved = 'no';
                pop_saveset(ALLEEG(dats{index}(idat)), 'savemode', 'resave');
                ALLEEG(dats{index}(idat)).saved = 'yes';
            end
        end
    else
        for index = 1:length(dats)
            ALLEEG(dats{index}) = pop_runica(ALLEEG(dats{index}), 'icatype', g.icatype, ...
                'options', g.options, 'chanind', g.chanind, 'concatenate', 'on');
            for idat = 1:length(dats{index})
                ALLEEG(dats{index}(idat)).saved = 'no';
                pop_saveset(ALLEEG(dats{index}(idat)), 'savemode', 'resave');
                ALLEEG(dats{index}(idat)).saved = 'yes';
                ALLEEG(dats{index}(idat)) = update_datafield(ALLEEG(dats{index}(idat)));
            end
        end
    end
    com = sprintf('EEG = pop_runica(EEG, %s);', ...
              vararg2str({ 'icatype' g.icatype 'concatcond' 'on' 'options' g.options }) );
    return;
else
    % copy data
    % ---------
    if all([ALLEEG(g.dataset).trials] == 1)
        EEG = eeg_checkset(ALLEEG, 'loaddata');
        EEG = pop_mergeset(EEG, 1:length(EEG));
    else
        disp('Concatenating datasets...');
        EEG = ALLEEG(g.dataset(1));
        
        cpnts = 1;
        % compute total data size
        % -----------------------
        totalpnts = 0;
        for i = g.dataset
            totalpnts = totalpnts+ALLEEG(g.dataset(i)).pnts*ALLEEG(g.dataset(i)).trials;
        end
        EEG.data = zeros(EEG.nbchan, totalpnts);
            
        for i = g.dataset
            tmplen = ALLEEG(g.dataset(i)).pnts*ALLEEG(g.dataset(i)).trials;
            TMP = eeg_checkset(ALLEEG(g.dataset(i)), 'loaddata');
            EEG.data(:,cpnts:cpnts+tmplen-1) = reshape(TMP.data, size(TMP.data,1), size(TMP.data,2)*size(TMP.data,3));
            cpnts = cpnts+tmplen;
        end
    end
    EEG.icaweights = [];
    EEG.trials = 1;
    EEG.pnts   = size(EEG.data,2);
    EEG.saved  = 'no';
    if isfield(EEG.etc, 'ic_classification')
        EEG.etc = rmfield(EEG.etc, 'ic_classification');
    end
end

% Store and then remove current EEG ICA weights and sphere
% ---------------------------------------------------
fprintf('\n');
if ~isempty(EEG.icaweights)
    fprintf('Saving current ICA decomposition in "EEG.etc.oldicaweights" (etc.).\n');
    if ~isfield(EEG,'etc'), EEG.etc = []; end
    if ~isfield(EEG.etc,'oldicaweights')
        EEG.etc.oldicaweights = {};
        EEG.etc.oldicasphere = {};
        EEG.etc.oldicachansind = {};
    end
    tmpoldicaweights  = EEG.etc.oldicaweights;
    tmpoldicasphere   = EEG.etc.oldicasphere;
    tmpoldicachansind = EEG.etc.oldicachansind;
    EEG.etc.oldicaweights = { EEG.icaweights    tmpoldicaweights{:} };
    EEG.etc.oldicasphere  = { EEG.icasphere     tmpoldicasphere{:}  };
    EEG.etc.oldicachansind  = { EEG.icachansind tmpoldicachansind{:}  };
    fprintf('               Decomposition saved as entry %d.\n',length(EEG.etc.oldicaweights));
end
EEG.icaweights = [];
EEG.icasphere  = [];
EEG.icawinv    = [];
EEG.icaact     = [];

% select sub_channels
% -------------------
if isempty(g.chanind)
    g.chanind = 1:EEG.nbchan;
end
if iscell(g.chanind)
    datatype = {EEG.chanlocs.type};
    tmpChanInd = [];
    for iChan = 1:length(datatype)
        if ~isempty(datatype{iChan}) && ~isempty(strmatch(datatype{iChan}, g.chanind))
            tmpChanInd = [ tmpChanInd iChan ];
        end
    end
    g.chanind = tmpChanInd;
end
EEG.icachansind = g.chanind;

% is pca already an option?
% -------------------------
pca_opt = 0;
for i = 1:length(g.options)
    if ischar(g.options{i})
        if strcmpi(g.options{i}, 'pca')
            pca_opt = 1;
            pca_ind = i;
        end
    end
end
if pca_opt
    fprintf([ 'Warning: you have used PCA to reduce dimensionality so ICA\n' ...
        '         is not modeling the entire data, only the PCA-reduced data.\n' ]);
end

%------------------------------
% compute ICA on a definite set
% -----------------------------
tmpdata = reshape( EEG.data(g.chanind,:,:), length(g.chanind), EEG.pnts*EEG.trials);
tmprank = getrank(double(tmpdata(:,1:min(3000, size(tmpdata,2)))), pca_opt);
tmpdata = tmpdata - repmat(mean(tmpdata,2), [1 size(tmpdata,2)]); % zero mean 
if ~strcmpi(g.icatype, 'binica')
    try
        disp('Attempting to convert data matrix to double precision for more accurate ICA results.')
        tmpdata = double(tmpdata);
        tmpdata = tmpdata - repmat(mean(tmpdata,2), [1 size(tmpdata,2)]); % zero mean (more precise than single precision)
    catch
        disp('*************************************************************')
        disp('Not enough memory to convert data matrix to double precision.')
        disp('All computations will be done in single precision. Matlab 7.x')
        disp('under 64-bit Linux and others is imprecise in this mode.')
        disp('We advise use of "binica" instead of "runica."')
        disp('*************************************************************')
    end
end
switch lower(g.icatype)
    case 'runica' 
        % make sure we are using the correct ICA function
        runicaCurrentLoc = fileparts(which('runica'));
        runicaDesiredLoc = fullfile(fileparts(fileparts(which('pop_runica'))), 'sigprocfunc');
        if ~isequal(runicaCurrentLoc, runicaDesiredLoc)
            addpath(runicaDesiredLoc); % put back to the beginning of the path
        end
        try 
            if ismatlab && nargin < 2, g.options = {  g.options{:}, 'interrupt', 'on' }; end
        catch, end
        if tmprank == size(tmpdata,1) || pca_opt
            [EEG.icaweights,EEG.icasphere] = runica( tmpdata, 'lrate', 0.001,  g.options{:} );
        else 
            if nargin < 2
                uilist = { { 'style' 'text' 'string' [ 'EEGLAB has detected that the rank of your data matrix' 10 ...
                                                       'is lower the number of input data channels. This might' 10 ...
                                                       'be because you are including a reference channel or' 10 ...
                                                       'because you are running a second ICA decomposition.' 10 ...
                                                       sprintf('The proposed dimension for ICA is %d (out of %d channels).', tmprank, size(tmpdata,1)) 10 ...
                                                       'Rank computation may be inaccurate so you may edit this' 10 ...
                                                       'number below. If you do not understand, simply press OK.' ] } { } ...
                           { 'style' 'text' 'string' 'Proposed rank:' } ...
                           { 'style' 'edit' 'string' num2str(tmprank) } };
                res = inputgui('uilist', uilist, 'geometry', { [1] [1] [1 1] }, 'geomvert', [6 1 1]);
                if isempty(res), return; end
                tmprank = str2num(res{1});
                g.options = [g.options  { 'pca' tmprank }];
            else
                g.options = [g.options  {'pca' tmprank }]; % automatic for STUDY (batch processing)
            end
            disp(['Data rank (' int2str(tmprank) ') is smaller than the number of channels (' int2str(size(tmpdata,1)) ').']);
            [EEG.icaweights,EEG.icasphere] = runica( tmpdata, 'lrate', 0.001, g.options{:} );
        end
     case 'binica'
        icadefs;
        fprintf(['Warning: If the binary ICA function does not work, check that you have added the\n' ...
                 'binary file location (in the EEGLAB directory) to your Unix /bin directory (.cshrc file)\n']);
        if exist(ICABINARY) ~= 2
            error('Pop_runica(): binary ICA executable not found. Edit icadefs.m file to specify the ICABINARY location');
        end
        tmprank = getrank(tmpdata(:,1:min(3000, size(tmpdata,2))));
        if tmprank == size(tmpdata,1) || pca_opt
            [EEG.icaweights,EEG.icasphere] = binica( tmpdata, 'lrate', 0.001, g.options{:} );
        else 
            disp(['Data rank (' int2str(tmprank) ') is smaller than the number of channels (' int2str(size(tmpdata,1)) ').']);
            [EEG.icaweights,EEG.icasphere] = binica( tmpdata, 'lrate', 0.001, 'pca', tmprank, g.options{:} );
        end
    case 'amica' 
         if ~exist('pop_runamica')
             if nargin < 2
                 errordlg2('You must install the AMICA plugin first to use AMICA');
             else
                error('You must install the AMICA plugin first to use AMICA');
             end
         end
         EEG = pop_runamica(EEG, g.options{:});
     case 'picard' 
         if ~exist('picard')
             if nargin < 2
                 errordlg2('You must install the picard plugin first to use Picard');
             else
                error('You must install the picard plugin first to use Picard');
             end
         end
         options2 = g.options;
         if pca_opt
             if g.options{pca_ind+1} < 0
                 fprintf('Decreasing data dimension from %d to %d using PCA...\n', size(tmpdata,1), size(tmpdata,1)+g.options{pca_ind+1});
                 [tmpdata,eigvec] = runpca(tmpdata, size(tmpdata,1)+g.options{pca_ind+1});
             else
                 fprintf('Decreasing data dimension from %d to %d using PCA...\n', size(tmpdata,1), g.options{pca_ind+1});
                 [tmpdata,eigvec] = runpca(tmpdata, g.options{pca_ind+1});
             end
             options2(pca_ind:pca_ind+1) = [];
         end
        [~, EEG.icaweights] = picard( tmpdata, 'verbose', true, options2{:});
         if pca_opt
            EEG.icaweights = EEG.icaweights*pinv(eigvec);
         end
     case 'pearson_ica' 
        if isempty(g.options)
            disp('Warning: EEGLAB default for pearson ICA is 1000 iterations and epsilon=0.0005');
            [~, EEG.icaweights] = pearson_ica( tmpdata, 'maxNumIterations', 1000,'epsilon',0.0005);
        else    
            [~, EEG.icaweights] = pearson_ica( tmpdata, g.options{:});
        end
     case 'egld_ica', disp('Warning: This algorithm is very slow!!!');
        [~, EEG.icaweights] = egld_ica( tmpdata, g.options{:} );
     case 'tfbss' 
        if  isempty(g.options)
             [~, EEG.icaweights] = tfbss( tmpdata, size(tmpdata,1), 8, 512 );
        else    
             [~, EEG.icaweights] = tfbss( tmpdata, g.options{:} );
        end
     case 'jader',         [EEG.icaweights] = jader( tmpdata, g.options{:} );
     case 'matlabshibbsr', [EEG.icaweights] = MatlabshibbsR( tmpdata, g.options{:} );
     case 'eea',           [EEG.icaweights] = eeA( tmpdata, g.options{:} );
     case 'icaml',         [~, EEG.icawinv] = icaML( tmpdata, g.options{:} );
     case 'icams',         [~, EEG.icawinv] = icaMS( tmpdata, g.options{:} );
     case 'fastica',       [~, EEG.icawinv, EEG.icaweights] = fastica( tmpdata, 'displayMode', 'off', g.options{:} );
     case { 'tica' 'erica' 'simbec' 'unica' 'amuse' 'fobi' 'evd' 'sons' ...
            'jadeop' 'jade_td_p' 'evd24' 'sobi' 'ng_ol' 'acsobiro' 'acrsobibpf' } 
        fig = figure('tag', 'alg_is_run', 'visible', 'off');
        
        if isempty(g.options), g.options = { size(tmpdata,1) }; end
        switch lower(g.icatype)
         case 'tica',     EEG.icaweights = tica( tmpdata, g.options{:} );
         case 'erica',    EEG.icaweights = erica( tmpdata, g.options{:} );
         case 'simbec',   EEG.icaweights = simbec( tmpdata, g.options{:} );
         case 'unica',    EEG.icaweights = unica( tmpdata, g.options{:} );
         case 'amuse',    EEG.icaweights = amuse( tmpdata );
         case 'fobi',     [~, EEG.icaweights] = fobi( tmpdata, g.options{:} );
         case 'evd',      EEG.icaweights = evd( tmpdata, g.options{:} );
         case 'sons',     EEG.icaweights = sons( tmpdata, g.options{:} );
         case 'jadeop',   EEG.icaweights = jadeop( tmpdata, g.options{:} );
         case 'jade_td_p',EEG.icaweights = jade_td_p( tmpdata, g.options{:} );
         case 'evd24',    EEG.icaweights = evd24( tmpdata, g.options{:} );
         case 'sobi',     EEG.icawinv    = sobi( tmpdata, g.options{:} );
         case 'ng_ol',    [~, EEG.icaweights] = ng_ol( tmpdata, g.options{:} );
         case 'acsobiro', EEG.icawinv   = acsobiro( tmpdata, g.options{:} );
         case 'acrsobibpf', EEG.icawinv = acrsobibpf( tmpdata, g.options{:} );
        end
        clear tmp;
        close(fig);
     otherwise, error('Pop_runica: unrecognized algorithm');
end

% update weight and inverse matrices etc...
% -----------------------------------------
if ~isempty(fig), try close(fig); catch, end; end
if isempty(EEG.icaweights)
    EEG.icaweights = pinv(EEG.icawinv);
end
if isempty(EEG.icasphere)
    EEG.icasphere  = eye(size(EEG.icaweights,2));
end
if isempty(EEG.icawinv)
    EEG.icawinv    = pinv(EEG.icaweights*EEG.icasphere); % a priori same result as inv
end

% Reorder components by variance
% ------------------------------
meanvar = sum(EEG.icawinv.^2).*sum(transpose((EEG.icaweights *  EEG.icasphere)*EEG.data(EEG.icachansind,:)).^2)/((length(EEG.icachansind)*EEG.pnts)-1);
[~, windex] = sort(meanvar);
windex = windex(end:-1:1); % order large to small
meanvar = meanvar(windex);
EEG.icaweights = EEG.icaweights(windex,:);
EEG.icawinv    = pinv( EEG.icaweights *  EEG.icasphere );
if ~isempty(EEG.icaact)
    EEG.icaact = EEG.icaact(windex,:,:);
end

% copy back data to datasets if necessary
% ---------------------------------------
if length(g.dataset) > 1
    for i = g.dataset
        ALLEEG(i).icaweights = EEG.icaweights;
        ALLEEG(i).icasphere  = EEG.icasphere;
        ALLEEG(i).icawinv    = EEG.icawinv;
        ALLEEG(i).icachansind = g.chanind;
    end          
    ALLEEG = eeg_checkset(ALLEEG);
else
    EEG = eeg_checkset(EEG);
    ALLEEG = eeg_store(ALLEEG, EEG, g.dataset);
end

if nargin < 2 || selectamica
    if ~isempty(g.options)
        com = sprintf('EEG = pop_runica(EEG, ''icatype'', ''%s'', %s);', g.icatype, vararg2str(g.options) ); %vararg2str({ 'icatype' g.icatype 'dataset' g.dataset 'options' g.options }) );
    else
        com = sprintf('EEG = pop_runica(EEG, ''icatype'', ''%s'');',g.icatype );
    end
end

return;

function tmprank2 = getrank(tmpdata, pca_opt)
    
    tmprank = rank(tmpdata);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Here: alternate computation of the rank by Sven Hoffman
    %tmprank = rank(tmpdata(:,1:min(3000, size(tmpdata,2)))); old code
    covarianceMatrix = cov(tmpdata', 1);
    [~, D] = eig (covarianceMatrix);
    rankTolerance = 1e-7;
    tmprank2=sum (diag (D) > rankTolerance);
    if tmprank ~= tmprank2
        if nargin >= 2 && pca_opt ~= 0
            fprintf('Warning: fixing rank computation inconsistency (%d vs %d) most likely because running under Linux 64-bit Matlab, taking the minimum\n', tmprank, tmprank2);
        end
        tmprank2 = min(tmprank, tmprank2);
    end
            
function EEG = update_datafield(EEG)
    if ~isfield(EEG, 'datfile'), EEG.datfile = ''; end
    if ~isempty(EEG.datfile)
        EEG.data = EEG.datfile;
    else 
        EEG.data = 'in set file';
    end
    EEG.icaact = [];            