% PURPOSE  : 	Averages bin-based ERPsets (grand average)
%
% FORMAT   :
%
% ERP = pop_gaverager( ALLERP , Parameters);
%
% EXAMPLE  :
%
% ERP = pop_gaverager( ALLERP , 'Criterion',100, 'ERPindex',1:4, 'Stdev',...
% 'on', 'Warning', 'on', 'Weighted', 'on' );
%
% Or
%
% ERP = pop_gaverager( ALLERP , 'Criterion',100, 'Loadlist',...
% 'C:\Users\Work\Documents\MATLAB\test.txt', 'Stdev', 'on', 'Warning', 'on',
%
% INPUTS   :
%
% ALLERP            - structure array of ERP structures (ERPsets)
%
% PARAMETERS     :
%
% 'ERPindex'        - index(es) pointing to ERP structures within ALLERP
% 'Loadlist'        - name of a text file containing the list of ERPset filenames
%                    (with path)
% 'Weighted'        - 'on' means apply weighted-average, 'off' means classic average.
% 'Stdev'           - Get standard deviation. 'on' or 'off'
% 'Warning'         - Warning 'on' or 'off'
% 'Criterion'       - Max allowed mean artifact detection proportion
%
% OUTPUTS  :
%
% ERP               - data structure containing the average of all specified datasets.
%
%
% GUI: grandaveragerGUI.m ; SUBROUTINE: --
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon & Steven Luck
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2009

%b8d3721ed219e65100184c6b95db209bb8d3721ed219e65100184c6b95db209b
%
% ERPLAB Toolbox
% Copyright © 2007 The Regents of the University of California
% Created by Javier Lopez-Calderon and Steven Luck
% Center for Mind and Brain, University of California, Davis,
% javlopez@ucdavis.edu, sjluck@ucdavis.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 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/>.

function [ERP erpcom] = pop_gaverager(ALLERP, varargin)

erpcom = '';
ERP = preloadERP;
if nargin<1
        help pop_gaverager
        return
end
if nargin==1  % GUI
        if isstruct(ALLERP)
                iserp = iserpstruct(ALLERP(1));
                if ~iserp
                        ALLERP = [];
                end
                actualnset = length(ALLERP);
        else
                ALLERP = [];
                actualnset = 0;
        end
        
        def  = erpworkingmemory('pop_gaverager');
        if isempty(def)
                def = { actualnset 0 '' 100 0 1};
        else
                def{1}=actualnset;
        end
        if isnumeric(def{3}) && ~isempty(ALLERP)
                if max(def{3})>length(ALLERP)
                        def{3} = def{3}(def{3}<=length(ALLERP));
                end
                if isempty(def{3})
                        def{3} = 1;
                end
        end
        
        %
        % Open Grand Average GUI
        %
        answer = grandaveragerGUI(def);
        if isempty(answer)
                disp('User selected Cancel')
                return
        end
        optioni  = answer{1}; %1 means from a filelist, 0 means from erpsets menu
        erpset   = answer{2};
        artcrite = answer{3}; % max percentage of rej to be included in the gaverage
        wavg     = answer{4}; % 0;1
        stdev    = answer{5}; % 0;1
        jk       = answer{6}; % 0;1
        jkerpname  = answer{7}; % erpname for JK grand averages
        jkfilename = answer{8}; % filename for JK grand averages
        
        if optioni==1 % from files
                filelist = erpset;
                erpindx = [];
        else % from erpsets
                filelist = '';
                erpindx  = erpset;
        end
        
        def = {actualnset, optioni, erpset, artcrite, wavg, stdev};
        erpworkingmemory('pop_gaverager', def);
        if stdev==1
                stdsstr = 'on';
        else
                stdsstr = 'off';
        end
        if wavg==1
                wavgstr = 'on';
        else
                wavgstr = 'off';
        end
        if jk==1
                %jkstr = 'on';
                [ALLERP erpcom]  = pop_jkgaverager(ALLERP, 'ERPindex', erpindx, 'Loadlist', filelist,'Criterion', artcrite,...
                        'Stdev', stdsstr, 'Weighted', wavgstr, 'Erpname', jkerpname, 'Filename', jkfilename);
                
                assignin('base','ALLERP', ALLERP);  % save to workspace
                updatemenuerp(ALLERP); % add a new erpset to the erpset menu
        else
                %jkstr = 'off';
                [ERP erpcom]  = pop_gaverager(ALLERP, 'ERPindex', erpindx, 'Loadlist', filelist,'Criterion', artcrite,...
                        'Stdev', stdsstr, 'Weighted', wavgstr, 'Saveas', 'on');
        end
        pause(0.1)
        return
else  % Scripting
        p = inputParser;
        p.FunctionName  = mfilename;
        p.CaseSensitive = false;
        p.addRequired('ALLERP');
        
        p.addParamValue('ERPindex', [], @isnumeric);
        p.addParamValue('Loadlist', '');
        p.addParamValue('Criterion', 100, @isnumeric);  %  0-100; threshold of artifacts, in %, for delivering a warning
        p.addParamValue('Weighted', 'off', @ischar);  % 'on', 'off'
        p.addParamValue('Stdev', 'off', @ischar);     % 'on', 'off'
        p.addParamValue('Saveas', 'off', @ischar);    % 'on', 'off'
        %p.addParamValue('Jackknife','off', @ischar); % 'on', 'off'
        p.addParamValue('Warning', 'on', @ischar);    % 'on', 'off'
        
        p.parse(ALLERP, varargin{:});
        
        erpsetArray = p.Results.ERPindex;
        filelist    = p.Results.Loadlist;
        artcrite    = p.Results.Criterion;
        
        if ismember({p.Results.Stdev}, {'on','yes'})
                stdev    = 1;
        else
                stdev    = 0;
        end
        if ismember({p.Results.Weighted}, {'on','yes'})
                wavg    = 1;
        else
                wavg    = 0;
        end
        if ismember({p.Results.Saveas}, {'on','yes'})
                saveas    = 1;
        else
                saveas    = 0;
        end
        if ismember({p.Results.Warning}, {'on','yes'})
                warn    = 1;
        else
                warn    = 0;
        end
        if iserpstruct(ALLERP);
                if isempty(erpsetArray) && isempty(filelist)
                        error('ERPLAB:pop_gaverager', '\nERPLAB says: ERPindex is empty')
                end
        else
                if isempty(filelist)
                        error('ERPLAB:pop_gaverager', '\nERPLAB says: Loadlist is empty (ALLERP is empty as well...)\n')
                end
        end
end
if isempty(erpsetArray) && ~isempty(filelist)
        if iscell(filelist)
                inputfname = strtrim(filelist);
                nfile = length(inputfname);
        else
                disp(['pop_gaverager(): For file-List, user selected ', filelist])
                try
                        fid_list   = fopen( filelist );
                        inputfname   = textscan(fid_list, '%[^\n]','CommentStyle','#');
                        fclose(fid_list);
                        
                        inputfname    = strtrim(cellstr(inputfname{:}));
                        nfile = length(inputfname);
                catch
                        error('ERPLAB:pop_gaverager', '\nERPLAB says: "%s" couldn''t be loaded', filelist)
                end
        end
        optioni=1;
elseif  ~isempty(erpsetArray) && isempty(filelist)
        nfile  = length(erpsetArray);
        optioni=0;
elseif  isempty(erpsetArray) && isempty(filelist)
        error('ERPLAB:pop_gaverager', '\nERPLAB says: Both ERPindex and Loadlist are empty')
else
        error('ERPLAB:pop_gaverager', '\nERPLAB says: ERPindex and Loadlist were specified.ERPLAB got confused by this...')
end

% ERPaux = ERP; % original ERP
naccepted = [];
nrejected = [];
ninvalid  = [];
narflags  = [];
ds=1;

%
% Reads and sums ERPsets
%
for j=1:nfile
        if optioni==1
                ERPTX = load(inputfname{j}, '-mat');
                ERPT  = ERPTX.ERP;
        else
                ERPT = ALLERP(erpsetArray(j));
        end
        
        [ERPT conti serror] = olderpscan(ERPT, warn);
        
        if conti==0
                break
        end
        if serror
                msgboxText =  ['Your erpset %s is not compatible at all with the current ERPLAB version.\n'...
                        'Please, try upgrading your ERP structure.'];
                title = 'ERPLAB: erp_loaderp() Error';
                errorfound(sprintf(msgboxText, ERP.filename), title);
                break
        end
        if j>1
                % basic test for number of channels (for now...)
                if  pre_nchan ~= ERPT.nchan
                        msgboxText =  sprintf('Erpsets #%g and #%g have different number of channels!', j-1, j);
                        title = 'ERPLAB: pop_gaverager() Error';
                        errorfound(msgboxText, title);
                        return
                end
                
                % basic test for number of points (for now...)
                if pre_pnts  ~= ERPT.pnts
                        msgboxText =  sprintf('Erpsets #%g and #%g have different number of points!', j-1, j);
                        title = 'ERPLAB: pop_gaverager() Error';
                        errorfound(msgboxText, title);
                        return
                end
                
                % basic test for number of bins (for now...)
                if pre_nbin  ~= ERPT.nbin
                        msgboxText =  sprintf('Erpsets #%g and #%g have different number of bins!', j-1, j);
                        title = 'ERPLAB: pop_gaverager() Error';
                        errorfound(msgboxText, title);
                        return
                end
        end
        
        pre_nchan = ERPT.nchan;
        pre_pnts  = ERPT.pnts;
        pre_nbin  = ERPT.nbin;
        
        %
        % Sum ERPsets
        %
        if isempty(ERPT.pexcluded)
                ERPT.pexcluded = 0;
        end
        
        %
        % Artifact detection threshold test
        %        
        if ERPT.pexcluded> artcrite  % test the mean artifact detection prop per ERPset
                callwarning(sprintf('%d) ', j),'');
                fprintf(' ERPset # %d "%s" has %.1f %% of mean artifact detection.\n', j, ERPT.erpname, ERPT.pexcluded);
                fprintf(' ERP #%g is being included anyway...\n', j);
        else
                fprintf(' %d) Including ERP #%g...\n', j, j);
        end
        
        %
        % Preparing ERP struct (first ERP)
        %
        if ds==1
                workfileArray    = ERPT.workfiles;
                [nch npnts nbin] = size(ERPT.bindata);
                sumERP           = zeros(nch,npnts,nbin);
                if stdev
                        sumERP2  = zeros(nch,npnts,nbin);
                end
                naccepted    = zeros(1,nbin);
                nrejected    = zeros(1,nbin);
                ninvalid     = zeros(1,nbin);
                narflags     = zeros(nbin,8);
                ERP.erpname  = [];
                ERP.filename = [];
                ERP.filepath = [];
                ERP.subject  = ERPT.subject;
                ERP.nchan    = ERPT.nchan;
                ERP.nbin     = ERPT.nbin;
                ERP.pnts     = ERPT.pnts;
                ERP.srate    = ERPT.srate;
                ERP.xmin     = ERPT.xmin;
                ERP.xmax     = ERPT.xmax;
                ERP.times    = ERPT.times;
                ERP.bindata  = [];
                ERP.binerror = [];
                ERP.chanlocs = ERPT.chanlocs;
                ERP.ref      = ERPT.ref;
                ERP.bindescr = ERPT.bindescr;
                ERP.history  = ERPT.history;
                ERP.saved    = ERPT.saved;
                ERP.isfilt   = ERPT.isfilt;
                ERP.version  = ERPT.version;
                EVEL         = ERPT.EVENTLIST;
                ALLEVENTLIST(1:length(EVEL)) = EVEL;
        else
                workfileArray = [workfileArray ERPT.workfiles];
                EVEL =  ERPT.EVENTLIST;
                ALLEVENTLIST(end+1:end+length(EVEL)) = EVEL;
        end
        
        countbinOK = [ERPT.ntrials.accepted]; % These work as weights
        
        if wavg
                for bb=1:pre_nbin
                        sumERP(:,:,bb)  = sumERP(:,:,bb)  + ERPT.bindata(:,:,bb).*countbinOK(bb);      % weighted sum: Sum(wi*xi)
                        if stdev
                                sumERP2(:,:,bb) = sumERP2(:,:,bb)  + (ERPT.bindata(:,:,bb).^2).*countbinOK(bb); % weighted sum of squares: Sum(wi*xi^2)
                        end
                end
        else
                sumERP    = sumERP + ERPT.bindata;              % general sum: Sum(xi)
                if stdev
                        sumERP2   = sumERP2 + ERPT.bindata.^2;    % general sum of squares: Sum(xi^2)
                end
        end
        
        naccepted = naccepted + countbinOK;                   % sum of (current) weights: Sum(wi)
        nrejected = nrejected + ERPT.ntrials.rejected;
        ninvalid  = ninvalid  + ERPT.ntrials.invalid;
        
        %
        % Sum flags counter per bin
        %
        if isfield(ERPT.ntrials, 'arflags')
                narflags  = narflags  + ERPT.ntrials.arflags;
        end
        
        ds = ds+1; % counter for clean ERPsets (keep it that way until Steve might eventually decide to reject subjects automatically...)
        %       else
        %             callwarning(sprintf('%d) ', j),'');
        %             fprintf(' ERPset # %d "%s" was skipped because it had %.1f %% of mean artifact detection.\n', j, ERPT.erpname, ERPT.pexcluded);
        %       end
end
if conti==0
        return
end
if serror==1
        fprintf('pop_gaverager() cancelled.\n');
        return
end
if ~exist('sumERP', 'var')
        error('ERPLAB:emptyERP','\nERPLAB says: Your grand average is empty which means none of your ERPset met the artifact proportion criterion you specified.\n\n')
end

ERP.erpname   = [];
ERP.workfiles = workfileArray;
ERP.EVENTLIST = ALLEVENTLIST;

%
% Gets the grand averaged ERPset
%
if wavg
        for bb=1:pre_nbin
                if naccepted(bb)>0
                        
                        ERP.bindata(:,:,bb) = sumERP(:,:,bb)./naccepted(bb); % get ERP!  --> weighted sum is divided by the sum of the weights
                        
                        if stdev
                                if bb==1
                                        fprintf('\nEstimating weighted standard deviation of data...\n');
                                end
                                insqrt = sumERP2(:,:,bb).*naccepted(bb) - sumERP(:,:,bb).^2;
                                if nnz(insqrt<0)>0
                                        ERP.binerror(:,:,bb)= zeros(nch,npnts,1);
                                else
                                        ERP.binerror(:,:,bb)= (1/naccepted(bb))*sqrt(insqrt);
                                end
                        end
                else
                        ERP.bindata(:,:,bb) = zeros(nch,npnts,1);            % makes a zero erp
                        if stdev
                                ERP.binerror(:,:,bb)= zeros(nch,npnts,1);
                        end
                end
        end
else
        ERP.bindata   = sumERP./nfile; % get ERP!  --> general sum is divided by the number of files (erpsets)
        if stdev
                fprintf('\nEstimating standard deviation of data...\n');
                ERP.binerror  = sqrt(sumERP2.*(1/nfile) - ERP.bindata.^2) ; % ERP stdev
        end
end
if ~isempty(ERP.binerror)
        if nnz(ERP.binerror>0)==0
                fprintf('\n*******************************************************\n');
                fprintf('WARNING: There is not variance in the data! So, ERP.binerror = [] \n');
                fprintf('*******************************************************\n\n');
                ERP.binerror = [];
        end
end

ERP.ntrials.accepted = naccepted;
ERP.ntrials.rejected = nrejected;
ERP.ntrials.invalid  = ninvalid;
ERP.ntrials.arflags  = narflags;
vali = round(1000*(sum(nrejected)/(sum(naccepted)+sum(nrejected))))/10;
ERP.pexcluded = vali;

[ERP serror] = sorterpstruct(ERP);
if serror
        error('ERPLAB:pop_gaverager','\nERPLAB says: Your ERPsets are not compatibles')
end
if wavg
        fprintf('\n %g ERPs were (weight) averaged.\n', ds-1);
else
        fprintf('\n %g ERPs were (classic) averaged.\n', ds-1);
end

%
% Completion statement
%
msg2end

if saveas
        [ERP issave erpcom_save] = pop_savemyerp(ERP,'gui','erplab');
else
        issave = 1; % for getting command history from scripting job
end
if issave>0
        %                 if iscell(artcrite)
        skipfields = {'ALLERP', 'Saveas'};
        fn     = fieldnames(p.Results);
        if ischar(ALLERP)
                inputvari = 'loadlist'; % Thanks to Felix Bacigalupo for this suggestion. Dic 12, 2011
        else
                inputvari = inputname(1);
        end
        erpcom = sprintf( 'ERP = pop_gaverager( %s ', inputvari);
        for q=1:length(fn)
                fn2com = fn{q};
                if ~ismember(fn2com, skipfields)
                        fn2res = p.Results.(fn2com);
                        if ~isempty(fn2res)                                
                                if ischar(fn2res)
                                        if ~strcmpi(fn2res,'off')
                                                erpcom = sprintf( '%s, ''%s'', ''%s''', erpcom, fn2com, fn2res);
                                        end
                                else
                                        if iscell(fn2res)
                                                if ischar([fn2res{:}])
                                                        fn2resstr = sprintf('''%s'' ', fn2res{:});
                                                else
                                                        fn2resstr = vect2colon(cell2mat(fn2res), 'Sort','on');
                                                end
                                                fnformat = '{%s}';
                                        else
                                                fn2resstr = vect2colon(fn2res, 'Sort','on');
                                                fnformat = '%s';
                                        end                                        
                                        if strcmpi(fn2com,'Criterion')
                                                if p.Results.Criterion<100
                                                        erpcom = sprintf( ['%s, ''%s'', ' fnformat], erpcom, fn2com, fn2resstr);
                                                end
                                        else
                                                erpcom = sprintf( ['%s, ''%s'', ' fnformat], erpcom, fn2com, fn2resstr);
                                        end
                                end
                        end
                end
        end
        erpcom = sprintf( '%s );', erpcom);
        
        %artcritestr = vect2colon(cell2mat(artcrite), 'Sort','yes','Delimiter','off');
        %erpcom = sprintf( 'ERP = pop_averager( %s, %s, {%s}, %s );', inputname(1), setindexstr,...
        %        artcritestr, num2str(stdev));
        %         else
        %
        %                         erpcom = sprintf( 'ERP = pop_averager( %s, %s, %s, %s );', inputname(1), setindexstr,...
        %                                 num2str(artcrite), num2str(stdev));
        %                 end
        if issave==2
                erpcom = sprintf('%s\n%s', erpcom, erpcom_save);
                msgwrng = '*** Your ERPset was saved on your hard drive.***';
                fprintf('\n%s\n\n', msgwrng)
        end
else
        msgwrng = 'ERPLAB Warning: Your changes were not saved';
        try cprintf([1 0.52 0.2], '%s\n\n', msgwrng); catch,fprintf('%s\n\n', msgwrng);end ;
        return
end
return
