% EEG_GETDATACT - get EEG data from a specified dataset or

%                  component activity

%

% Usage:

%       >> signal = eeg_getdatact( EEG );

%       >> signal = eeg_getdatact( EEG, 'key', 'val');

%

% Inputs:

%   EEG       - Input dataset

%

% Optional input:

%   'channel'   - [integer array] read only specific channels.

%                 Default is to read all data channels.

%   'component' - [integer array] read only specific components

%   'projchan'  - [integer or cell array] channel(s) onto which the component

%                 should be projected.

%   'rmcomps'   - [integer array] remove selected components from data

%                 channels. This is only to be used with channel data not

%                 when selecting components.

%   'trialindices' - [integer array] only read specific trials. Default is

%                 to read all trials.

%   'samples'   - [integer array] only read specific samples. Default is

%                 to read all samples.

%   'reshape'   - ['2d'|'3d'] reshape data. Default is '3d' when possible.

%   'verbose'   - ['on'|'off'] verbose mode. Default is 'on'.

%

% Outputs:

%   signal      - EEG data or component activity

%

% Author: Arnaud Delorme, SCCN & CERCO, CNRS, 2008-

%

% See also: EEG_CHECKSET

% Copyright (C) 15 Feb 2002 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.

function [data, boundaries] = eeg_getdatact( EEG, varargin)

data = [];

if nargin < 1

    help eeg_getdatact;

    return;

end

% reading data from several datasets and concatenating it

% -------------------------------------------------------

if iscell(EEG) || (~ischar(EEG) && length(EEG) > 1)

    % decode some arguments

    % ---------------------

    trials  = cell(1,length(EEG));

    rmcomps = cell(1,length(EEG));

    for iArg = length(varargin)-1:-2:1

        if strcmpi(varargin{iArg}, 'trialindices')

            trials = varargin{iArg+1};

            varargin(iArg:iArg+1) = [];

        elseif strcmpi(varargin{iArg}, 'rmcomps')

            rmcomps = varargin{iArg+1};

            varargin(iArg:iArg+1) = [];

        end

    end

    if isnumeric(rmcomps), rmtmp = rmcomps; rmcomps = cell(1,length(EEG)); rmcomps(:) = { rmtmp }; end

    % concatenate datasets

    % --------------------

    data = [];

    boundaries = [];

    for dat = 1:length(EEG)

        if iscell(EEG)

             [tmpdata, datboundaries] = eeg_getdatact(EEG{dat}, 'trialindices', trials{dat}, 'rmcomps', rmcomps{dat}, varargin{:} );

        else [tmpdata, datboundaries] = eeg_getdatact(EEG(dat), 'trialindices', trials{dat}, 'rmcomps', rmcomps{dat}, varargin{:} );

        end

        if isempty(data)

            data       = tmpdata;

            boundaries = datboundaries;

        else

            if (isstruct(EEG) && all([EEG.trials] == 1)) || size(tmpdata,3) == 1 || size(data,2) ~= size(tmpdata,2) % continuous data (if same number of data points, consider 1 trial dataset)

                if size(tmpdata,3) ~= 1 || size(data,3) ~= 1, error('Trying to concatenate continuous and epoched datasets'); end

                if size(data,1) ~= size(tmpdata,1), error('Datasets to be concatenated do not have the same number of channels'); end

                % adding boundaries

                if ~isempty(datboundaries)

                    boundaries = [boundaries datboundaries size(data,2)];

                else

                    boundaries = [boundaries size(data,2)];

                end

                data(:,end+1:end+size(tmpdata,2)) = tmpdata; % concatenating trials

            else

                if size(data,1) ~= size(tmpdata,1), error('Datasets to be concatenated do not have the same number of channels'); end

                if size(data,2) ~= size(tmpdata,2), error('Datasets to be concatenated do not have the same number of time points'); end

                data(:,:,end+1:end+size(tmpdata,3)) = tmpdata; % concatenating trials

            end

        end

    end

    return;

end

% if string load dataset

% ----------------------

if ischar(EEG)

    EEG = pop_loadset('filename', EEG, 'loadmode', 'info');

end

opt = finputcheck(varargin, { ...

    'channel'   'integer' {} [];

    'verbose'   'string'  { 'on','off' } 'on';

    'reshape'   'string'  { '2d','3d' }  '3d';

    'projchan'  {'integer','cell' } { {} {} } [];

    'component' 'integer' {} [];

    'samples'   'integer' {} [];

    'interp'    'struct'  { }        struct([]);

    'trialindices' {'integer','cell'} { {} {} } [];

    'rmcomps'      {'integer','cell'} { {} {} } [] }, 'eeg_getdatact');

if ischar(opt), error(opt); end

channelNotDefined = 0;

if isempty(opt.channel)

    opt.channel = [1:EEG.nbchan]; 

    channelNotDefined = 1;

elseif isequal(opt.channel, [1:EEG.nbchan]) && ~isempty(opt.interp)

    channelNotDefined = 1;

end

if iscell( opt.trialindices), opt.trialindices = opt.trialindices{1}; end

if isempty(opt.trialindices), opt.trialindices = [1:EEG.trials]; end

if iscell(opt.rmcomps     ), opt.rmcomps      = opt.rmcomps{1};      end

if (~isempty(opt.rmcomps) || ~isempty(opt.component)) && isempty(EEG.icaweights)

    error('No ICA weight in dataset');

end

if strcmpi(EEG.data, 'in set file')

    EEGTMP = pop_loadset('filename', EEG.filename, 'filepath', EEG.filepath, 'verbose', 'off');

    EEG.data = EEGTMP.data;

end

% get data boundaries if continuous data

% --------------------------------------

boundaries = [];

if nargout > 1 && EEG.trials == 1 && ~isempty(EEG.event) && isfield(EEG.event, 'type') && ischar(EEG.event(1).type)

    if ~isempty(opt.samples)

        disp('WARNING: eeg_getdatact.m, boundaries are not accurate when selecting data samples');

    end

    tmpevent = EEG.event;

    tmpbound = strmatch('boundary', lower({ tmpevent.type }));

    if ~isempty(tmpbound)

        boundaries = [tmpevent(tmpbound).latency ]-0.5;

    end

end

% getting channel or component activation

% ---------------------------------------

filename = fullfile(EEG.filepath, [ EEG.filename(1:end-4) '.icaact' ] );

if ~isempty(opt.component) && ~isempty(EEG.icaact)

    data = EEG.icaact(opt.component,:,:);

elseif ~isempty(opt.component) && exist(filename)

    % reading ICA file

    % ----------------

    data = repmat(single(0), [ length(opt.component) EEG.pnts EEG.trials ]);

    fid = fopen( filename, 'r', 'ieee-le'); %little endian (see also pop_saveset)

    if fid == -1, error( ['file ' filename ' could not be open' ]); end

    for ind = 1:length(opt.component)

        fseek(fid, (opt.component(ind)-1)*EEG.pnts*EEG.trials*4, -1);

        data(ind,:) = fread(fid, [EEG.trials*EEG.pnts 1], 'float32')';

    end

    fclose(fid);

elseif ~isempty(opt.component)

    if isempty(EEG.icaact)

        data = eeg_getdatact( EEG );

        data = (EEG.icaweights(opt.component,:)*EEG.icasphere)*data(EEG.icachansind,:);

    else

        data = EEG.icaact(opt.component,:,:);

    end

else

    if isnumeric(EEG.data) % channel

        data = EEG.data;

    else % channel but no data loaded string

        [~,EEG.data,ext] = fileparts(EEG.data); % remove path if present

        EEG.data = [ EEG.data ext];

        filename = fullfile(EEG.filepath, EEG.data);

        fid = fopen( filename, 'r', 'ieee-le'); %little endian (see also pop_saveset)

        if fid == -1

            filename = fullfile(EEG.data);

            fid = fopen( filename, 'r', 'ieee-le'); %little endian (see also pop_saveset)

            if fid == -1

                filename = fullfile(EEG.filepath, [EEG.filename(1:end-3) 'fdt' ]);

                fid = fopen( filename, 'r', 'ieee-le'); %little endian (see also pop_saveset)

                if fid == -1

                    error( ['file ' filename ' not found. If you have renamed/moved' 10 ...

                        'the .set file, you must also rename/move the associated data file.' ]);

                end

            end

        end

        if strcmpi(opt.verbose, 'on')

            fprintf('Reading float file ''%s''...\n', filename);

        end

        % old format = .fdt; new format = .dat (transposed)

        % -------------------------------------------------

        datformat = 0;

        if length(filename) > 3

            if strcmpi(filename(end-2:end), 'dat')

                datformat = 1;

            end

        end

        EEG.datfile = EEG.data;

        % reading data file

        % -----------------

        eeglab_options;

        if length(opt.channel) == EEG.nbchan && option_memmapdata

            fclose(fid);

            data = mmo(filename, [EEG.nbchan EEG.pnts EEG.trials], false);

            %data = memmapdata(filename, [EEG.nbchan EEG.pnts EEG.trials]);

        else

            if datformat

                if length(opt.channel) == EEG.nbchan || ~isempty(opt.interp)

                    data = fread(fid, [EEG.trials*EEG.pnts EEG.nbchan], 'float32')';

                else

                    data = repmat(single(0), [ length(opt.channel) EEG.pnts EEG.trials ]);

                    for ind = 1:length(opt.channel)

                        fseek(fid, (opt.channel(ind)-1)*EEG.pnts*EEG.trials*4, -1);

                        data(ind,:) = fread(fid, [EEG.trials*EEG.pnts 1], 'float32')';

                    end

                    opt.channel = [1:size(data,1)];

                end

            else

                data = fread(fid, [EEG.nbchan Inf], 'float32');

            end

            fclose(fid);

        end

    end

    % subtracting components from data

    % --------------------------------

    if ~isempty(opt.rmcomps)

        if strcmpi(opt.verbose, 'on')

            fprintf('Removing %d artifactual components\n', length(opt.rmcomps));

        end

        rmcomps = eeg_getdatact( EEG, 'component', opt.rmcomps); % loaded from file

        rmchan    = [];

        rmchanica = [];

        for index = 1:length(opt.channel)

            tmpicaind = find(opt.channel(index) == EEG.icachansind);

            if ~isempty(tmpicaind)

                rmchan    = [ rmchan    index ];

                rmchanica = [ rmchanica tmpicaind ];

            end

        end

        data(rmchan,:) = data(rmchan,:) - EEG.icawinv(rmchanica,opt.rmcomps)*rmcomps(:,:);

        %EEG = eeg_checkset(EEG, 'loaddata');

        %EEG = pop_subcomp(EEG, opt.rmcomps);

        %data = EEG.data(opt.channel,:,:);

        %if strcmpi(EEG.subject, 'julien') & strcmpi(EEG.condition, 'oddball') & strcmpi(EEG.group, 'after')

        %    jjjjf

        %end

    end

    if ~isempty(opt.interp)

        EEG.data   = data;

        EEG.event  = [];

        EEG.epoch  = [];

        EEG = eeg_interp(EEG, opt.interp, 'spherical');

        data = EEG.data;

        if channelNotDefined, opt.channel = [1:EEG.nbchan]; end

    end

    if ~isequal(opt.channel, [1:EEG.nbchan])

        data = data(intersect(opt.channel,[1:EEG.nbchan]),:,:);

    end

end

% projecting components on data channels

% --------------------------------------

if ~isempty(opt.projchan)

    if iscell(opt.projchan)

        opt.projchan = std_chaninds(EEG, opt.projchan);

    end

    finalChanInds = [];

    for iChan = 1:length(opt.projchan)

        tmpInd = find(EEG.icachansind == opt.projchan(iChan));

        if isempty(tmpInd)

            error(sprintf('Warning: can not backproject component on channel %d (not used for ICA)\n', opt.projchan(iChan)));

        end

        finalChanInds = [ finalChanInds tmpInd ];

    end

    data = EEG.icawinv(finalChanInds, opt.component)*data(:,:);

end

if size(data,2)*size(data,3) ~= EEG.pnts*EEG.trials

    disp('WARNING: The file size on disk does not correspond to the dataset, file has been truncated');

end

try

    if EEG.trials == 1, EEG.pnts = size(data,2); end

    if  strcmpi(opt.reshape, '3d')

         data = reshape(data, size(data,1), EEG.pnts, EEG.trials);

    else data = reshape(data, size(data,1), EEG.pnts*EEG.trials);

    end

catch

    error('The file size on disk does not correspond to the dataset information.');

end

% select trials

% -------------

if length(opt.trialindices) ~= EEG.trials

    data = data(:,:,opt.trialindices);

end

if ~isempty(opt.samples)

    data = data(:,opt.samples,:);

end