% STD_LIMO - Export and run in LIMO the EEGLAB STUDY design.
%           call limo_batch to create all 1st level LIMO_EEG analysis 
%
% Usage:
%   [STUDY LIMO_files] = std_limo(STUDY,ALLEEG,'key',val)
%
% Inputs:
%  STUDY        - studyset structure containing some or all files in ALLEEG
%  ALLEEG       - vector of loaded EEG datasets
%
% Optional inputs:
%  'measure'      - ['daterp'|'icaerp'|'datspec'|'icaspec'|'datersp'|'icaersp']
%                   measure to compute. Currently, only 'daterp' and
%                   'datspec' are supported. Default is 'daterp'.
%  'method'       - ['OLS'|'WTS'|'IRLS'] Ordinary Least Squares (OLS) or Weighted
%                   Least Squares (WTS) or Iterative Reweighted Least Squares'IRLS'.
%                   WTS should be used as it is more robust. IRLS is much slower
%                   and better across subjects than across trials.
%  'design'       - [integer] design index to process. Default is the current
%                   design stored in STUDY.currentdesign.
%  'contnan'      - ['on'|'off'] NaN for continuous variables. When 'on'
%                   NaNs are allowed (and trials removed). When NaNs are 
%                   replaced with 0s.
%  'erase'        - ['on'|'off'] erase previous files. Default is 'on'.
%  'neighboropt'  - [cell] cell array of options for the function computing
%                   the channel neighbox matrix STD_PREPARE_CHANLOCS. The file
%                   is saved automatically if channel location are present.
%                   This option allows to overwrite the defaults when computing
%                   the channel neighbox matrix.
%   'chanloc'     - Channel location structure. Must be used with 'neighbormat',
%                   or it will be ignored. If this option is used, it will
%                   ignore 'neighboropt' if used.
%   'neighbormat' - Neighborhood matrix of electrodes. Must be used with 'chanloc',
%                   or it will be ignored. If this option is used, it will
%                   ignore 'neighboropt' if used.
%   'freqlim'     - Frequency trimming in Hz
%   'timelim'     - Time trimming in millisecond
%
% Outputs:
%  STUDY     - modified STUDY structure (the STUDY.design now contains a list
%              of the limo files)
%  LIMO_files a structure with the following fields
%     LIMO_files.LIMO the LIMO folder name where the study is analyzed
%     LIMO_files.mat a list of 1st level LIMO.mat (with path)
%     LIMO_files.Beta a list of 1st level Betas.mat (with path)
%     LIMO_files.con a list of 1st level con.mat (with path)
%     LIMO_files.expected_chanlocs expected channel location neighbor file for
%                                  correcting for multiple comparisons
% Example:
%  [STUDY LIMO_files] = std_limo(STUDY,ALLEEG,'measure','daterp')
%
% Author: Arnaud Delorme (SCCN) & Cyril Pernet (LIMO Team)
%         Based on previous version from Ramon Martinez-Cancino and Arnaud Delorme
%
% Copyright (C) 2018 Arnaud Delorme
%
% 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 [STUDY, LIMO_files] = std_limo(STUDY,ALLEEG,varargin)

LIMO_files = [];

if isempty(STUDY.filepath)
    STUDY.filepath = pwd;
end
cd(STUDY.filepath);
if nargin < 2
    help std_limo;
    return;
end

warning('off', 'MATLAB:lang:cannotClearExecutingFunction');
if ischar(varargin{1}) && ( strcmpi(varargin{1}, 'daterp') || ...
        strcmpi(varargin{1}, 'datspec') || ...
        strcmpi(varargin{1}, 'dattimef') || ...
        strcmpi(varargin{1}, 'icaerp')|| ...
        strcmpi(varargin{1}, 'icaspec')|| ...
        strcmpi(varargin{1}, 'icatimef'))
    opt.measure  = varargin{1};
    opt.design   = varargin{2};
    opt.ow_chanlocfile = 'no';  % if chanloc file exist, do not overwrite
    opt.erase          = 'on';  % erase previous folders/file with the same name
    opt.method         = 'WLS'; % weighted least squares by default
    opt.zscore         = 1;     % zscore regressors
else
    opt = finputcheck( varargin, ...
        { 'measure'        'string'  { 'daterp' 'datspec' 'dattimef' 'icaerp' 'icaspec' 'icatimef' } 'daterp'; ...
        'method'         'string'  { 'OLS' 'WLS' 'IRLS' } 'WLS';
        'design'         'integer' [] STUDY.currentdesign;
        'erase'          'string'  { 'on','off' }   'off';
        'splitreg'       'string'  { 'on','off' }   'off';
        'interaction'    'string'  { 'on','off' }   'off';
        'contnan'        'string'  { 'on','off' }   'off';
        'freqlim'        'real'    []               [] ;
        'timelim'        'real'    []               [] ;
        'neighboropt'    'cell'    {}               {} ;
        'chanloc'     'struct'  {}  struct('no', {});
        'neighbormat' 'real'    []  [];
        'zscore'         'real'    [0,1]            1  ;
        'ow_chanlocfile' 'string'  {'yes','no'}     'no'},...
        'std_limo');
    if ischar(opt), error(opt); end
end
opt.measureori = opt.measure;
if strcmpi(opt.measure, 'datersp')
    opt.measure = 'dattimef';
end

Analysis     = opt.measure;
design_index = opt.design;

% Make sure paths are ok for LIMO (Consider to move this to eeglab.m in a future)
% -------------------------------------------------------------------------
root = fileparts(which('limo_eeg'));
pathCell = regexp(path, pathsep, 'split');
onPath = all([sum(strcmp([root filesep 'help'],pathCell))~=0,...
    sum(strcmp([root filesep 'limo_cluster_functions'],pathCell))~=0,...
    sum(strcmp([root filesep 'external' filesep 'psom'],pathCell))~=0,...
    sum(strcmp([root filesep 'deprecated'], pathCell))~=0]);
if onPath == 0
    addpath([root filesep 'limo_cluster_functions'])
    addpath([root filesep 'external'])
    addpath([root filesep 'external' filesep 'psom'])
    addpath([root filesep 'help'])
    addpath([root filesep 'deprecated'])
end

% Checking fieldtrip paths to compute gp channel location
skip_chanlocs   = 0;
chanloc_created = 0;
limoChanlocs    = [];
if exist(fullfile([STUDY.filepath filesep 'derivatives'], 'limo_gp_level_chanlocs.mat'),'file')
    limoChanlocs = fullfile([STUDY.filepath filesep 'derivatives'], 'limo_gp_level_chanlocs.mat');
elseif exist(fullfile(STUDY.filepath, 'limo_gp_level_chanlocs.mat'),'file')
    limoChanlocs = fullfile(STUDY.filepath, 'limo_gp_level_chanlocs.mat');
elseif exist(fullfile([STUDY.filepath filesep 'derivatives'], 'limo_chanlocs.mat'),'file')
    limoChanlocs = fullfile([STUDY.filepath filesep 'derivatives'], 'limo_chanlocs.mat');
elseif exist(fullfile(STUDY.filepath, 'limo_chanlocs.mat'),'file')
    limoChanlocs = fullfile(STUDY.filepath, 'limo_chanlocs.mat');
end

if ~isempty(limoChanlocs)
    if ~strcmpi(opt.ow_chanlocfile,'no') % empty or yes
        opt.ow_chanlocfile = questdlg2('channel location file found, do you want to overwrite','overwrite?','yes','no','no');
    end

    if isempty(opt.ow_chanlocfile) || strcmpi(opt.ow_chanlocfile,'no')
        skip_chanlocs = 1;
    end
end

if skip_chanlocs == 0
    if ~exist('ft_prepare_neighbours','file')
        warndlg('std_limo error: Fieldtrip extension should be installed - chanlocs NOT generated', '', 'non-modal');
        skip_chanlocs = 1;
    else
        if ~exist('eeglab2fieldtrip','file')
            root = fileparts(which('ft_prepare_neighbours'));
            addpath([root filesep 'external' filesep 'eeglab']);
        end
    end
end

% Detecting type of analysis
% -------------------------------------------------------------------------
model.defaults.datatype = opt.measureori(4:end);
if ~isempty(strfind(Analysis,'dat')) %#ok<STREMP>
    model.defaults.type = 'Channels';
elseif  ~isempty(strfind(Analysis,'ica')) %#ok<STREMP>
    [STUDY,flags]=std_checkdatasession(STUDY,ALLEEG);
    if sum(flags)>0
        error('some subjects have data from different sessions - can''t do ICA');
    end
    model.defaults.type = 'Components';
end

% Checking if clusters
% -------------------------------------------------------------------------
if strcmp(model.defaults.type,'Components')
    if isempty(STUDY.cluster(1).child)
        warning('Components have not been clustered, LIMO will not match them across subjects');
        model.defaults.icaclustering = 0;
    else
        model.defaults.icaclustering = 1;
    end
end

% computing channel neighbour matrix
% ---------------------------------
if skip_chanlocs == 0
    chanloc_created = 1;
    if isempty(opt.chanloc) && isempty(opt.neighbormat)
        if isfield(ALLEEG(1).chanlocs, 'theta') &&  ~strcmp(model.defaults.type,'Components')
            if  ~isfield(STUDY.etc,'statistic')
                STUDY = pop_statparams(STUDY, 'default');
            end

            try
                [~,~,limoChanlocs] = std_prepare_neighbors(STUDY, ALLEEG, 'force', 'on', opt.neighboropt{:});
                chanlocname = 'limo_gp_level_chanlocs.mat';
            catch neighbors_error
                limoChanlocs = []; chanloc_created = 0;
                warndlg2(neighbors_error.message,'limo_gp_level_chanlocs.mat not created', 'non-modal')
            end
        else
            limoChanlocs = []; chanloc_created = 0;
            if ~isempty(STUDY.cluster(1).child)
                disp('Warning: cannot compute expected channel distance for correction for multiple comparisons');
            end
        end
    else
        limoChanlocs.expected_chanlocs   = opt.chanloc;
        limoChanlocs.channeighbstructmat = opt.neighbormat;
        chanlocname = 'limo_chanlocs.mat';
    end
end

if chanloc_created
    % contains will not work in Octave
    if isempty(strfind(STUDY.filepath,'derivatives'))
        if ~exist([STUDY.filepath filesep 'derivatives'],'dir')
            mkdir([STUDY.filepath filesep 'derivatives']);
        end
        limoChanlocsFile = fullfile([STUDY.filepath filesep 'derivatives'], chanlocname);
    else
        limoChanlocsFile = fullfile(STUDY.filepath, chanlocname);
    end
    % this sometimes happen to be nested in expected_chanlocs, fixing it here
    if all(arrayfun(@(x) any(strcmp(x,{'expected_chanlocs','channeighbstructmat'})), fieldnames(limoChanlocs.expected_chanlocs)))
        limoChanlocs = limoChanlocs.expected_chanlocs;
    end
    save('-mat', limoChanlocsFile, '-struct', 'limoChanlocs');
    fprintf('Saving channel neighbors for correction for multiple comparisons in \n%s\n', limoChanlocsFile);
end

% 1st level analysis
% -------------------------------------------------------------------------
model.cat_files  = [];
model.cont_files = [];

% Cleaning old files from the current design (Cleaning ALL)
% -------------------------------------------------------------------------
% NOTE: Clean up the .lock files to (to be implemented)
% [STUDY.filepath filesep 'derivatives' filesep 'limo_batch_report']
if strcmp(opt.erase,'on')
    [~,filename] = fileparts(STUDY.filename);
    std_limoerase(STUDY.filepath, filename, STUDY.subject, num2str(STUDY.currentdesign));
    STUDY.limo = [];
end

% Check if the measures has been computed
% also find out if the channels are interpolated
% -------------------------------------------------------------------------
interpolated = zeros(1,length(STUDY.datasetinfo));
for iDat = 1:length(STUDY.datasetinfo)
    fileName = fullfile(STUDY.datasetinfo(iDat).filepath, [ STUDY.datasetinfo(iDat).subject '*.' opt.measure ]);
    % fileName should already match unless user moves / rename, hence using dir
    fileName = dir(fileName);
    if isempty(fileName)
        error('std_limo subject %s: Measures must be computed first',STUDY.datasetinfo(iDat).subject);
    else
        if strcmp(model.defaults.type,'Channels')
            tmpChans = load('-mat', fullfile(fileName(1).folder,fileName(1).name), 'labels');
            if length(tmpChans.labels) > ALLEEG(iDat).nbchan, interpolated(iDat) = 1; end
        end
    end
end

measureflags = struct('daterp','off',...
    'datspec','off',...
    'datersp','off',...
    'dattimef','off',...
    'datitc' ,'off',...
    'icaerp' ,'off',...
    'icaspec','off',...
    'icatimef','off',...
    'icaersp','off',...
    'icaitc','off');
measureflags.(lower(opt.measureori))= 'on';
STUDY.etc.measureflags = measureflags;
mergedChanlocs = eeg_mergelocs(ALLEEG.chanlocs);
fprintf('generating temporary files, pulling relevant trials ... \n')

% generate temporary merged datasets needed by LIMO
% -------------------------------------------------
allSubjects    = { STUDY.datasetinfo.subject };
allSessions    = { STUDY.datasetinfo.session };
uniqueSubjects = unique(allSubjects);
nb_subjects    = length(uniqueSubjects);
allSessions(cellfun(@isempty, allSessions)) = { 1 };
allSessions    = cellfun(@num2str, allSessions, 'uniformoutput', false);
uniqueSessions = unique(allSessions);

% by default we create a design matrix with all condition
factors = pop_listfactors(STUDY.design(opt.design), 'gui', 'off', 'level', 'one', 'constant', 'off');

for iSubj = 1:nb_subjects
    for iSess = 1:length(uniqueSessions)
        inds1 = strmatch( uniqueSubjects{iSubj}, allSubjects, 'exact');
        inds2 = strmatch( uniqueSessions{iSess}, allSessions, 'exact');
        inds  = intersect(inds1, inds2);
        if ~isempty(inds)
            if length(inds) ~= 1
                error([ 'Cannot calculate contrast because more than 1 dataset per session' 10 ...
                    'per subject. Merge datasets for each subject and try again.' ]);
            end
            
            % make file-up
            [~,subname] = fileparts(STUDY.datasetinfo(inds).filename);
            if isfield(ALLEEG,'filename')
                if ~strcmp(subname,ALLEEG(inds).filename(1:end-4))
                    error('STUDY and ALLEEG mismatch, can''t figure out which file to use')
                end
            else
                warning('No filename in ALLEEG, pulling data blindly from STUDY')
            end

            if contains(subname,'sub-') 
                if contains(subname,'ses-')
                    filename = [subname '_design' num2str(design_index) '.set'];
                else
                    filename = [subname '_ses-' num2str(iSess) '_design' num2str(design_index) '.set'];
                end
            else
                if contains(subname,'ses-')
                    filename = ['sub-' subname '_design' num2str(design_index) '.set'];
                else
                    filename = ['sub-' subname '_ses-' num2str(iSess) '_design' num2str(design_index) '.set'];
                end
            end

            % Creating fields for limo
            % ------------------------
            fprintf('pulling trials for %s ... \n',filename)
            EEGTMP                     = std_lm_seteegfields(STUDY,ALLEEG(inds), inds,'datatype',model.defaults.type,'format', 'cell');
            ALLEEG                     = eeg_store(ALLEEG, EEGTMP, inds);
            file_fullpath              = rel2fullpath(STUDY.filepath,ALLEEG(inds).filepath);
            model.set_files{inds}      = fullfile(file_fullpath , filename);

            OUTEEG = [];
            if all([ALLEEG(inds).trials] == 1)
                OUTEEG.trials = 1;
            else
                OUTEEG.trials = sum([ALLEEG(inds).trials]);
            end

            filepath_tmp               = rel2fullpath(STUDY.filepath,ALLEEG(inds).filepath);
            OUTEEG.filepath            = filepath_tmp;
            OUTEEG.filename            = filename;
            OUTEEG.srate               = ALLEEG(inds).srate;
            OUTEEG.icaweights          = ALLEEG(inds).icaweights;
            OUTEEG.icasphere           = ALLEEG(inds).icasphere;
            OUTEEG.icawinv             = ALLEEG(inds).icawinv;
            OUTEEG.icachansind         = ALLEEG(inds).icachansind;
            OUTEEG.etc                 = ALLEEG(inds).etc;
            OUTEEG.times               = ALLEEG(inds).times;
            if any(interpolated)
                OUTEEG.chanlocs        = mergedChanlocs;
                OUTEEG.etc.interpolatedchannels = setdiff(1:length(OUTEEG.chanlocs), std_chaninds(OUTEEG, { ALLEEG(inds).chanlocs.labels }));
            else
                OUTEEG.chanlocs        = ALLEEG(inds).chanlocs;
            end

            % update EEG.etc
            OUTEEG.etc.merged{1}       = ALLEEG(inds).filename;

            % Def fields
            OUTEEG.etc.datafiles.daterp   = [];
            OUTEEG.etc.datafiles.datspec  = [];
            OUTEEG.etc.datafiles.datersp  = [];
            OUTEEG.etc.datafiles.dattimef = [];
            OUTEEG.etc.datafiles.datitc   = [];
            OUTEEG.etc.datafiles.icaerp   = [];
            OUTEEG.etc.datafiles.icaspec  = [];
            OUTEEG.etc.datafiles.icaersp  = [];
            OUTEEG.etc.datafiles.icatimef = [];
            OUTEEG.etc.datafiles.icaitc   = [];

            % Filling fields
            single_trials_filename = EEGTMP.etc.datafiles.(opt.measureori);
            if exist(single_trials_filename,'file')
                if strcmpi(measureflags.daterp,'on')
                    OUTEEG.etc.datafiles.daterp = single_trials_filename;
                elseif strcmpi(measureflags.datspec,'on')
                    OUTEEG.etc.datafiles.datspec = single_trials_filename;
                elseif strcmpi(measureflags.datersp,'on')
                    OUTEEG.etc.datafiles.datersp = single_trials_filename;
                elseif strcmpi(measureflags.datitc,'on')
                    OUTEEG.etc.datafiles.datitc = single_trials_filename;
                elseif strcmpi(measureflags.dattimef,'on')
                    OUTEEG.etc.datafiles.dattimef = single_trials_filename;
                elseif strcmpi(measureflags.icaerp,'on')
                    OUTEEG.etc.datafiles.icaerp = single_trials_filename;
                elseif strcmpi(measureflags.icaspec,'on')
                    OUTEEG.etc.datafiles.icaspec = single_trials_filename;
                elseif strcmpi(measureflags.icaersp,'on')
                    OUTEEG.etc.datafiles.icaersp = single_trials_filename;
                elseif strcmpi(measureflags.icaitc,'on')
                    OUTEEG.etc.datafiles.icaitc = single_trials_filename;
                elseif strcmpi(measureflags.icatimef,'on')
                    OUTEEG.etc.datafiles.icatimef = single_trials_filename;
                end
            end

            % Save info
            EEG = OUTEEG;
            save('-mat', fullfile( filepath_tmp, OUTEEG.filename), 'EEG');
            clear OUTEEG filepath_tmp

            % generate data files
            % -------------------
            fprintf('making up statistical model for %s ... \n',filename)
            % save continuous and categorical data files
            trialinfo = std_combtrialinfo(STUDY.datasetinfo, inds);
            % [catMat,contMat,limodesign] = std_limodesign(factors, trialinfo, 'splitreg', opt.splitreg, 'interaction', opt.interaction);
            [catMat,contMat,limodesign] = std_limodesign(factors, trialinfo, 'splitreg', 'off', 'interaction', opt.interaction);
            if strcmpi(opt.splitreg,'on') && ~isempty(contMat)
                for c=size(contMat,2):-1:1
                    splitreg{c} = limo_split_continuous(catMat,contMat(:,c)); % std_limodesign does something else when splitting regressors
                end
                contMat    = cell2mat(splitreg);
                opt.zscore = 0; % regressors are now zscored
            end

            % replace NaNs if necessary
            if strcmpi(opt.contnan, 'off')
                contMat(isnan(contMat)) = 0;
            end

            % copy results
            model.cat_files{inds}                 = catMat;
            model.cont_files{inds}                = contMat;
            if isfield(limodesign, 'categorical')
                STUDY.limo.categorical = limodesign.categorical;
            else
                STUDY.limo.categorical = {};
            end
            if isfield(limodesign, 'continuous')
                STUDY.limo.continuous = limodesign.continuous;
            else
                STUDY.limo.continuous = {};
            end
            STUDY.limo.subjects(inds).subject     = STUDY.datasetinfo(inds(1)).subject;
            STUDY.limo.subjects(inds).cat_file    = catMat;
            STUDY.limo.subjects(inds).cont_file   = contMat;
        end
    end % exit session
end % exit subject

% then we add contrasts for conditions that were merged during design selection
% i.e. multiple categorical variables (factors) and yet not matching the number
% of variables (contrasts are then a weighted sum of the crossed factors)
if ~isempty(factors) && isfield(factors, 'value') && ...
        sum(arrayfun(@(x) ~strcmpi(x.label,'group'),STUDY.design(opt.design).variable)) == 1 % only one non-continuous variable other than group
    if length(STUDY.design(opt.design).variable(1).value) ~= length(factors) % and this var has more values than the number of factors
        limocontrast = zeros(length(STUDY.design(opt.design).variable(1).value),length(factors)+1); % length(factors)+1 to add the constant
        for n=length(factors):-1:1
            factor_names{n} = factors(n).value;
        end

        index = find(arrayfun(@(x) ~strcmpi(x.label,'group'),STUDY.design(opt.design).variable)); % which one is not group
        for c=1:length(STUDY.design(opt.design).variable(index).value)
            limocontrast(c,1:length(factors)) = single(ismember(factor_names,STUDY.design(opt.design).variable(index).value{c}));
            limocontrast(c,1:length(factors)) = limocontrast(c,1:length(factors)) ./ sum(limocontrast(c,1:length(factors))); % scale by the number of variables
        end
    end
end

% transpose
model.set_files  = model.set_files';
model.cat_files  = model.cat_files';
model.cont_files = model.cont_files';
if all(cellfun(@isempty, model.cat_files )), model.cat_files  = []; end
if all(cellfun(@isempty, model.cont_files)), model.cont_files = []; end


% set model.defaults - all conditions no bootstrap
% -----------------------------------------------------------------
if strcmp(Analysis,'daterp') || strcmp(Analysis,'icaerp')
    model.defaults.analysis = 'Time';
    for s=nb_subjects:-1:1
        vs(s) = ALLEEG(s).xmin*1000;
        ve(s) = ALLEEG(s).xmax*1000;
    end
    model.defaults.start    = max(vs);
    model.defaults.end      = min(ve);

    if length(opt.timelim) == 2 && opt.timelim(1) < opt.timelim(end)
        % start value
        if opt.timelim(1) < model.defaults.start
            fprintf('std_limo: Invalid time lower limit, using default value instead');
        else
            model.defaults.start = opt.timelim(1);
        end
        % end value
        if opt.timelim(end) > model.defaults.end
            fprintf('std_limo: Invalid time upper limit, using default value instead');
        else
            model.defaults.end = opt.timelim(end);
        end
    end
    model.defaults.lowf  = [];
    model.defaults.highf = [];

elseif strcmp(Analysis,'datspec') || strcmp(Analysis,'icaspec')

    model.defaults.analysis= 'Frequency';
    if length(opt.freqlim) == 2
        model.defaults.lowf    = opt.freqlim(1);
        model.defaults.highf   = opt.freqlim(2);
    else
        error('std_limo: Frequency limits need to be specified');
    end
    model.defaults.start    = [];
    model.defaults.end      = [];

elseif strcmp(Analysis,'dattimef') || any(strcmp(Analysis,{'icatimef','icaersp'}))
    
    model.defaults.analysis = 'Time-Frequency';
    for s=nb_subjects:-1:1
        vs(s) = ALLEEG(s).xmin*1000;
        ve(s) = ALLEEG(s).xmax*1000;
    end
    model.defaults.start    = max(vs);
    model.defaults.end      = min(ve);
    model.defaults.lowf     = [];
    model.defaults.highf    = [];

    if length(opt.timelim) == 2
        model.defaults.start    = opt.timelim(1);
        model.defaults.end      = opt.timelim(2);
    end
    if length(opt.freqlim) == 2
        model.defaults.lowf     = opt.freqlim(1);
        model.defaults.highf    = opt.freqlim(2);
    else
        error('std_limo: Frequency limits need to be specified');
    end
end

model.defaults.fullfactorial    = 0;                 % all variables
model.defaults.zscore           = opt.zscore;        % done that already
model.defaults.bootstrap        = 0 ;                % only for single subject analyses - not included for studies
model.defaults.tfce             = 0;                 % only for single subject analyses - not included for studies
model.defaults.method           = opt.method;        % default is WLS 
model.defaults.Level            = 1;                 % 1st level analysis
model.defaults.type_of_analysis = 'Mass-univariate'; % option can be multivariate (work in progress)
model.defaults.labels           = factors;

if ~exist('limocontrast','var')
    [LIMO_files, procstatus] = limo_batch('model specification',model,[],STUDY);
else
    contrast.mat = limocontrast;
    [LIMO_files, procstatus] = limo_batch('both',model,contrast,STUDY);
    if exist(fullfile([STUDY.filepath filesep 'derivatives']),'dir')
        save(fullfile([STUDY.filepath filesep 'derivatives'],[STUDY.design(opt.design).name '_contrast.mat']),'limocontrast');
    else
        save(fullfile(STUDY.filepath,[STUDY.design(opt.design).name '_contrast.mat']),'limocontrast');
    end
end

STUDY.limo.model         = model;
STUDY.limo.datatype      = Analysis;
STUDY.limo.chanloc       = limoChanlocs;
STUDY.limo.betas         = pop_listfactors(STUDY, 'gui', 'off', 'level', 'one', 'splitreg', opt.splitreg, 'interaction', opt.interaction);
if exist('limocontrast','var')
    STUDY.limo.contrast      = limocontrast;
end

% generate between session contrasts
% ----------------------------------
index = 1;
for s = 1:nb_subjects
   sess_index = find(cellfun(@(x) strcmpi(x,uniqueSubjects{s}), allSubjects));
   % matches sess_index = find(contains(LIMO_files.mat,uniqueSubjects{s}))
   if length(sess_index) > 1
        fprintf('std_limo, computing additional between sessions contrasts for subject %s\n',uniqueSubjects{s})
        sess_name = allSessions(sess_index);
        pairs = nchoosek(1:length(sess_index),2); % do all session pairs
        for p=1:size(pairs,1)
            strpair = [cell2mat(sess_name(pairs(p,1))) cell2mat(sess_name(pairs(p,2)))];
            strpair(isspace(strpair)) = []; % remove spaces
            filesout{p} = limo_contrast_sessions(cell2mat(LIMO_files.mat(sess_index(pairs(p,1)))), ...
                cell2mat(LIMO_files.mat(sess_index(pairs(p,2)))),strpair);
        end
        
        for f=1:length(filesout)
            for ff=1:length(filesout{f})
                allcon{index} = filesout{f}{ff};
                index = index +1;
            end
        end
        clear filesout
    end
end

% use same glm_name as limo_batch
design_name = STUDY.design(STUDY.currentdesign).name;
design_name(isspace(design_name)) = [];
if findstr(design_name,'STUDY.')
    design_name = design_name(7:end);
end
glm_name = [STUDY.filename(1:end-6) '_' design_name '_GLM_' model.defaults.type '_' model.defaults.analysis '_' model.defaults.method];

% further split that list per regressor and group
if exist('allcon','var')
    maxcon = max(cellfun(@(x) str2double(x(strfind(x,'con_')+4:strfind(x,'sess_')-1)),allcon));
    for con=1:maxcon
        index = find(cellfun(@(x) ~isempty(x),cellfun(@(x) strfind(x,['con_' num2str(con)]),allcon','UniformOutput',false)));
        cell2csv([LIMO_files.LIMO filesep 'Between_sessions_con_' num2str(con) '_' glm_name '.txt'], allcon(index)');
        if length(STUDY.group) > 1
            for g= 1:length(STUDY.group)
                % find subjects of group g
                subset = find(arrayfun(@(x)(strcmpi(x.group,STUDY.group{g})), STUDY.datasetinfo));
                for s=1:length(subset)
                    sub{s} = STUDY.datasetinfo(subset(s)).subject;
                end
                sub = unique(sub);
                % find subjects of group g and contrast con
                subcon = [];
                for s = 1:length(sub)
                    if strfind(sub{s},'sub-')
                        subindex = find(cellfun(@(x) ~isempty(x),(cellfun(@(x) strfind(x,sub{s}),allcon','UniformOutput',false)))); % subject s group g
                    else
                        subindex = find(cellfun(@(x) ~isempty(x),(cellfun(@(x) strfind(x,['sub-' sub{s} ]),allcon','UniformOutput',false)))); % subject s group g
                    end
                    subcon = [subcon;intersect(index,subindex)];
                end
                % save
                if ~isempty(subcon)
                    cell2csv([LIMO_files.LIMO filesep 'Between_sessions_con_' num2str(con) 'Gp' STUDY.group{g} '_' glm_name '.txt'], allcon(subcon)');
                end
            end
        end
    end
end

% Save STUDY - delete tmp files
% ------------------------------
cd(STUDY.filepath);
STUDY      = pop_savestudy( STUDY, [],'filepath', STUDY.filepath,'savemode','resave');
keep_files = 'no';
if all(procstatus)
    disp('All subjects have been successfully processed.')
else
    if sum(procstatus)==0 % not a WLS issue - limo_batch errors for that and tells the user
        errordlg2('all subjects failed to process, check limo batch report')
    else
        warndlg2('some subjects failed to process, check limo batch report','', 'non-modal')
    end
    % cleanup temp files - except for subjects without errors
    db = dbstack;
    if length(db) <= 2
        keep_files = questdlg('Do you want to keep temp files of unsuccessulfully processed subjects','option for manual debugging','yes','no','no');
    end
end

% delete
if isempty(keep_files) || strcmpi(keep_files,'no')
    for s = 1:nb_subjects
        delete(model.set_files{s});
    end
else
    for s = find(procstatus)
        delete(model.set_files{s});
    end
end

% -------------------------------------------------------------------------
% Return full path if 'filepath' is a relative path. The output format will
% fit the one of 'filepath'. That means that if 'filepath' is a cell array,
% then the output will a cell array too, and the same if is a string.
function file_fullpath = rel2fullpath(studypath,filepath)

nit = 1; if iscell(filepath), nit = length(filepath);end

for i = 1:nit
    if iscell(filepath)
        pathtmp = filepath{i};
    else
        pathtmp = filepath;
    end

    if ~isempty(pathtmp) && contains(pathtmp(end),filesep) %#ok<STRIFCND>
        pathtmp = pathtmp(1:end-1);
    end % Getting rid of filesep at the end

    if ~isempty(pathtmp) && (contains(pathtmp(1:2),['.' filesep]) || ...
            (isunix && pathtmp(1) ~= '/') || (ispc && pathtmp(2) ~= ':')) %#ok<STREMP>
        if iscell(filepath)
            file_fullpath{i} = fullfile(studypath,pathtmp(1:end));
        else
            file_fullpath = fullfile(studypath,pathtmp(1:end));
        end
    else
        if iscell(filepath)
            file_fullpath{i} = pathtmp;
        else
            file_fullpath = pathtmp;
        end
    end
end

% get file base name (from std_precomp)
% ------------------
function filebase = getfilename(filepath, subj, sess, fileSuffix, onlyOneSession)
if onlyOneSession
    filebase = fullfile(filepath, [ subj fileSuffix ] );
else
    sesStr   = [ '0' sess ];
    filebase = fullfile(filepath, [ subj '_ses-' sesStr(end-1:end) fileSuffix ] );
end
