palm_takeargs.m

function [opts,plm] = palm_takeargs(varargin)
% Handle the inputs for PALM.
%
% _____________________________________
% Anderson M. Winkler
% FMRIB / University of Oxford
% Oct/2014
% http://brainder.org

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% PALM -- Permutation Analysis of Linear Models
% Copyright (C) 2015 Anderson M. Winkler
%
% 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
% 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

% Load the defaults
opts = palm_defaults;

% As varargin is actually from another function, fix it.
if nargin == 1
    if exist(varargin{1},'file')
        vararginx = palm_configrw(varargin{1});
    else
        error('Unknown option or file not found: %s',varargin{1});
    end
else
    vararginx = varargin;
    idxa = find(strcmpi(vararginx,'-o'));
    if isempty(idxa)
        otmp = opts.o;
    else
        otmp = vararginx{idxa+1};
    end
    if ~ strcmp(otmp(end),'_')
        otmp = horzcat(otmp,'_');
    end
    cfgname = horzcat(otmp,'palmconfig.txt');
    [opth,~,~] = fileparts(cfgname);
    if ~isempty(opth) && ~exist(opth,'dir')
        mkdir(opth);
    end
    palm_configrw(vararginx,cfgname);
end

% Number of input images/masks/surfaces
% These are NOT meant to be edited.
Ni     = sum(strcmp(vararginx,'-i'));        % number of data inputs
Nm     = sum(strcmp(vararginx,'-m'));        % number of masks
Ns     = sum(strcmp(vararginx,'-s'));        % number of surfaces
Nd     = sum(strcmp(vararginx,'-d'));        % number of design files
Nimiss = sum(strcmp(vararginx,'-imiss'));    % number of missing indicators for inputs
Ndmiss = sum(strcmp(vararginx,'-dmiss'));    % number of missing indicators for designs
Nt     = sum(strcmp(vararginx,'-t'));        % number of t-contrast files
Nf     = sum(strcmp(vararginx,'-f'));        % number of F-test files
Ncon   = sum(strcmp(vararginx,'-con'));      % number of contrast files (t or F, mset format)
Nevd   = sum(strcmp(vararginx,'-evperdat')); % number of EV per datum inputs
opts.i     = cell(Ni,1);   % Input files (to constitute Y later)
opts.m     = cell(Nm,1);   % Mask file(s)
opts.s     = cell(Ns,1);   % Surface file(s)
opts.sa    = cell(Ns,1);   % Area file(s) or weight(s)
opts.d     = cell(Nd,1);   % Design file(s)
opts.imiss = cell(Nd,1);   % Design file(s)
opts.dmiss = cell(Nd,1);   % Design file(s)
opts.t     = cell(Nt,1);   % t contrast file(s)
opts.f     = opts.t;       % F contrast file(s)
opts.Ccon  = cell(Ncon,1); % Contrast file(s) (t or F, mset format)
opts.Dcon  = cell(Ncon,1); % Contrast file(s) (multivariate, mset format)
opts.eb       = [];       % File with definition of exchangeability blocks
opts.vg       = [];       % File with definition of variance groups
opts.EE       = false;    % To be filled below (don't edit this!)
opts.ISE      = false;    % To be filled below (don't edit this!)
opts.within   = false;    % To be filled below (don't edit this!)
opts.whole    = false;    % To be filled below (don't edit this!)
opts.conskipcount = 0;    % When saving the contrasts, skip how many from 1?
opts.singlevg = true;     % Make sure that sigle VG will be used if nothing is supplied (this is NOT a "default" setting, and it's not a setting at all, but hard coded. Don't edit it!)
opts.subjidx  = [];       % Filename of the indices of subjects to keep
plm.subjidx   = [];       % Indices of subjects to keep

% These are to be incremented below
i = 1; m = 1; d = 1;
t = 1; s = 1;
con = 1; ev = 1;
imiss = 1; dmiss = 1;

% Remove trailing empty arguments. This is useful for some Octave versions.
while numel(vararginx) > 0 && isempty(vararginx{1})
    vararginx(1) = [];
end
narginx = numel(vararginx);

% Take the input arguments
a = 1;
while a <= narginx
    switch vararginx{a}
        case {'-help','-?','-basic','-advanced'}
            
            % Do nothing, as these options are parsed separately,
            % and should anyway be given without any other argument.
            a = a + 1;
            
        case '-i' % basic
            
            % Get the filenames for the data.
            opts.i{i} = vararginx{a+1};
            i = i + 1;
            a = a + 2;
            
        case '-m' % basic
            
            % Get the filenames for the masks, if any.
            opts.m{m} = vararginx{a+1};
            m = m + 1;
            a = a + 2;

        case {'-s','-surf'} % basic
            
            % Get the filenames for the surfaces, if any.
            opts.s{s}  = vararginx{a+1};
            if nargin == a+1 || (narginx>a+1 && strcmp(vararginx{a+2}(1),'-'))
                opts.sa{s} = [];
                a = a + 2;
            else
                opts.sa{s} = vararginx{a+2};
                a = a + 3;
            end
            s = s + 1;
            
        case '-d' % basic
            
            % Get the design matrix file.
            opts.d{d} = vararginx{a+1};
            d = d + 1;
            a = a + 2;
            
        case '-evperdat' % advanced
            
            % Use one EV per datum?
            opts.evperdat = true;
            opts.evdatfile{ev} = vararginx{a+1};
            if nargin == a + 1 || ...
                    ischar(vararginx{a+2}) && ...
                    strcmpi(vararginx{a+2}(1),'-')
                opts.evpos{ev}(1) = 1;  % EV position
                opts.evpos{ev}(2) = 1;  % Design number
                a = a + 2;
            elseif nargin == a + 2 || ...
                    ischar(vararginx{a+3}) && ...
                    strcmpi(vararginx{a+3}(1),'-')
                if ischar(vararginx{a+2}) % EV position
                    opts.evpos{ev}(1) = eval(vararginx{a+2});
                else
                    opts.evpos{ev}(1) = vararginx{a+2};
                end
                opts.evpos{ev}(2) = 1;  % Design number
                a = a + 3;
            else
                if ischar(vararginx{a+2}) % EV position
                    opts.evpos{ev}(1) = eval(vararginx{a+2});
                else
                    opts.evpos{ev}(1) = vararginx{a+2};
                end
                if ischar(vararginx{a+3}) % Design number
                    opts.evpos{ev}(2) = eval(vararginx{a+3});
                else
                    opts.evpos{ev}(2) = vararginx{a+3};
                end
                a = a + 4;
            end
            ev = ev + 1;
            
        case '-imiss' % basic
            
            % Get the filenames for the missing data indicators (inputs).
            opts.imiss{imiss} = vararginx{a+1};
            imiss = imiss + 1;
            a = a + 2;
            
        case '-dmiss' % basic
            
            % Get the filenames for the missing data indicators (designs).
            opts.dmiss{dmiss} = vararginx{a+1};
            dmiss = dmiss + 1;
            a = a + 2;
            
        case '-mcar'
            
            % For the missing data, treat as missing completely at random.
            opts.mcar = true;
            a = a + 1;
            
        case '-t' % basic
            
            % Get the t contrast files.
            opts.t{t} = vararginx{a+1};
            t = t + 1;
            a = a + 2;
            
        case '-f' % basic
            
            % Get the F contrast files.
            if t == 1
                error('The option "-f" cannot be specified before its respective "-t".');
            end
            opts.f{t-1} = vararginx{a+1};
            a = a + 2;
            
        case '-con' % advanced
            
            % Get the contrast files from an .mset file or
            % pair of files. If a pair, the 1st is for Cset
            % and the second for Dset.
            opts.Ccon{con} = vararginx{a+1};
            if nargin == a + 1 || ...
                    ischar(vararginx{a+2}) && ...
                    strcmpi(vararginx{a+2}(1),'-')
                opts.Dcon{con} = [];
                a = a + 2;
            else
                opts.Dcon{con} = vararginx{a+2};
                a = a + 3;
            end
            con = con + 1;
            
        case '-conskipcount' % advanced
            
            % Numbers to skip when saving the contrasts
            opts.conskipcount = vararginx{a+1};
            if ischar(opts.conskipcount)
                opts.conskipcount = str2double(opts.conskipcount);
            end
            a = a + 2;
            
        case '-tonly' % advanced
            
            % Run only the t-contrasts
            opts.tonly = true;
            a = a + 1;
            
        case '-fonly' % basic
            
            % Run only the F-contrasts
            opts.fonly = true;
            a = a + 1;
            
        case '-eb' % basic
            
            % Get the exchangeability blocks file.
            opts.eb = vararginx{a+1};
            a = a + 2;
            
        case '-vg' % basic
            
            % Get the variance groups file.
            opts.vg = vararginx{a+1};
            if     ischar(opts.vg) && ...
                    any(strcmpi(opts.vg,{'single'}))
                opts.vg = 'single';
                opts.singlevg = true;
            elseif ischar(opts.vg) && ...
                    any(strcmpi(opts.vg,{'auto','automatic'}))
                opts.vg = 'auto';
                opts.singlevg = false;
            else
                opts.singlevg = false;
            end
            a = a + 2;
            
        case '-swe' % advanced
            
            % Compute one (of various possible) sandwich estimators
            opts.SwE = true;
            a = a + 1;
            
        case '-o' % basic
            
            % Output prefix for the files to be saved.
            opts.o = vararginx{a+1};
            a = a + 2;
            
        case '-n' % basic
            
            % Number of permutations
            opts.nP0 = vararginx{a+1};
            if ischar(opts.nP0)
                opts.nP0 = str2double(opts.nP0);
            end
            a = a + 2;
            
        case '-C' % basic
            
            % Threshold for cluster extent, univariate, NPC and MV
            opts.cluster.uni.do = true;
            opts.cluster.uni.thr = vararginx{a+1};
            if ischar(opts.cluster.uni.thr)
                opts.cluster.uni.thr = str2double(opts.cluster.uni.thr);
            end
            opts.cluster.npc.do = true;
            opts.cluster.npc.thr = vararginx{a+1};
            if ischar(opts.cluster.npc.thr)
                opts.cluster.npc.thr = str2double(opts.cluster.npc.thr);
            end
            opts.cluster.mv.do  = true;
            opts.cluster.mv.thr = vararginx{a+1};
            if ischar(opts.cluster.mv.thr)
                opts.cluster.mv.thr = str2double(opts.cluster.mv.thr);
            end
            a = a + 2;
            
        case '-Cuni' % advanced
            
            % Threshold for cluster statistic, univariate
            opts.cluster.uni.do = true;
            opts.cluster.uni.thr = vararginx{a+1};
            if ischar(opts.cluster.uni.thr)
                opts.cluster.uni.thr = str2double(opts.cluster.uni.thr);
            end
            a = a + 2;
            
        case '-Cnpc' % advanced
            
            % Threshold for cluster statistic, NPC
            opts.NPC = true;
            opts.cluster.npc.do = true;
            opts.cluster.npc.thr = vararginx{a+1};
            if ischar(opts.cluster.npc.thr)
                opts.cluster.npc.thr = str2double(opts.cluster.npc.thr);
            end
            a = a + 2;
            
        case '-Cmv' % advanced
            
            % Threshold for cluster statistic, MV
            opts.MV = true;
            opts.cluster.mv.do = true;
            opts.cluster.mv.thr = vararginx{a+1};
            if ischar(opts.cluster.mv.thr)
                opts.cluster.mv.thr = str2double(opts.cluster.mv.thr);
            end
            a = a + 2;
            
        case '-Cstat' % advanced
            
            % Type of cluster statistic
            opts.cluster.stat = vararginx{a+1};
            if ~ any(strcmp(opts.cluster.stat,{'extent','mass','density','tippett','pivotal'}))
                error('Cluster statistic "%s" unknown.',opts.cluster.stat);
            end
            a = a + 2;
            
        case '-T' % basic
            
            % Do TFCE for univariate, NPC and MV?
            opts.tfce.uni.do = true;
            opts.tfce.npc.do = true;
            opts.tfce.mv.do = true;
            opts.tfce.stat = 'tfce';
            a = a + 1;
            
        case '-Tstat' % not in the help
            
            % Type of cluster statistic
            opts.tfce.stat = vararginx{a+1};
            if ~ any(strcmp(opts.tfce.stat,{'tfce','density','tippett'}))
                error('TFCE statistic "%s" unknown.',opts.tfce.stat);
            end
            a = a + 2;
            
        case '-Tuni' % advanced
            
            % Do TFCE for uni?
            opts.tfce.uni.do = true;
            a = a + 1;
            
        case '-Tnpc' % advanced
            
            % Do TFCE for NPC?
            opts.NPC = true;
            opts.tfce.npc.do = true;
            a = a + 1;
            
        case '-Tmv' % advanced
            
            % Do TFCE for MV?
            opts.MV = true;
            opts.tfce.mv.do = true;
            a = a + 1;
            
        case {'-tfce1D','-tfce1d'} % basic
            
            % Shortcut for -tfce_H 2 -tfce_E 2 -tfce_C 6,
            % i.e., parameters for TFCE in 2D mode
            opts.tfce.H      = 2;
            opts.tfce.E      = 2;
            opts.tfce.conn   = 6;
            a = a + 1;
            
        case {'-tfce2D','-tfce2d'} % basic
            
            % Shortcut for -tfce_H 2 -tfce_E 1 -tfce_C 26,
            % i.e., parameters for TFCE in 2D mode
            opts.tfce.H      = 2;
            opts.tfce.E      = 1;
            opts.tfce.conn   = 26;
            a = a + 1;
            
        case {'-tfce_H','-tfce_h'} % advanced
            
            % TFCE H parameter
            opts.tfce.H = vararginx{a+1};
            if ischar(opts.tfce.H)
                opts.tfce.H = str2double(opts.tfce.H);
            end
            a = a + 2;
            
        case {'-tfce_E','-tfce_e'} % advanced
            
            % TFCE E parameter
            opts.tfce.E = vararginx{a+1};
            if ischar(opts.tfce.E)
                opts.tfce.E = str2double(opts.tfce.E);
            end
            a = a + 2;
            
        case {'-tfce_C','-tfce_c'} % advanced
            
            % TFCE connectivity
            opts.tfce.conn = vararginx{a+1};
            if ischar(opts.tfce.conn)
                opts.tfce.conn = str2double(opts.tfce.conn);
            end
            a = a + 2;
            
        case '-tfce_dh' % advanced
            
            % TFCE delta-h parameter
            opts.tfce.deltah = vararginx{a+1};
            if ischar(opts.tfce.deltah)
                if strcmpi(opts.tfce.deltah,'auto')
                    opts.tfce.deltah = 0;
                else
                    opts.tfce.deltah = str2double(opts.tfce.deltah);
                end
            end
            a = a + 2;
            
        case '-tableasvolume' % basic
            
            % Treat tables (e.g., CSV inputs) as volume, such that TFCE can
            % be calculated. This is useful for TFCE over timeseries.
            opts.tableasvolume = true;
            a = a + 1;
            
        case '-within' % basic
            
            % Define whether should permute blocks as a whole or not
            opts.within = true;
            a = a + 1;
            
        case '-whole' % basic
            
            % Define whether should permute blocks as a whole or not
            opts.whole = true;
            a = a + 1;
            
        case '-ee' % basic
            
            % Exchangeable errors (EE)?
            % If yes, this means permutations.
            opts.EE = true;
            a = a + 1;
            
        case '-ise' % basic
            
            % Independent and symmetric errors (ISE)?
            % If yes, this means sign-flippings.
            opts.ISE = true;
            a = a + 1;
            
        case '-cmcp' % advanced
            
            % Define whether Conditional Monte Carlo should be used or not,
            % that is, ignoring repeated elements in the permutation set.
            opts.cmcp = true;
            a = a + 1;
            
        case '-cmcx' % advanced
            
            % Define whether repeated rows in X should be ignored or not
            % when defining the permutations, which constitutes another
            % form of CMC
            opts.cmcx = true;
            a = a + 1;
            
        case '-twotail' % basic
            
            % Do a two-tailed test for all t-contrasts?
            opts.twotail = true;
            a = a + 1;
            
        case '-concordant' % basic
            
            % For the NPC, favour alternatives with the same sign?
            opts.concordant = true;
            a = a + 1;
            
        case '-reversemasks' % basic
            
            % Reverse masks.
            opts.reversemasks = true;
            a = a + 1;
            
        case '-corrmod' % basic
            
            % Correct over modalities.
            opts.corrmod = true;
            a = a + 1;
            
        case '-corrcon' % basic
            
            % Correct over contrasts.
            opts.corrcon = true;
            a = a + 1;
            
        case '-saveparametric' % advanced
            
            % If the user wants to have also the parametric p-values.
            opts.savepara = true;
            a = a + 1;
            
        case '-saveglm' % advanced
            
            % If the user wants, save COPE and VARCOPEs in the 1st
            % permutation.
            opts.saveglm = true;
            a = a + 1;
            
        case {'-savecdf','-save1-p'} % basic
            
            % Save 1-p values (CDF) instead of the P-values
            opts.savecdf = true;
            a = a + 1;
            
        case '-logp'  % basic
            
            % Convert the P-values or (1-P)-values to -log10 before saving.
            opts.savelogp = true;
            a = a + 1;
            
        case '-savemask' % advanced
            
            % If the user wants to have also the masks used for each.
            % modality
            opts.savemask = true;
            a = a + 1;
            
        case '-rmethod' % advanced
            
            % Which method to use for the regression/permutation?
            if narginx > a
                methlist = {           ...
                    'Draper-Stoneman', ...
                    'Still-White',     ...
                    'Freedman-Lane',   ...
                    'terBraak',        ...
                    'Kennedy',         ... % Kennedy won't be in the help
                    'Manly',           ...
                    'Huh-Jhun',        ...
                    'Dekker'};
                methidx = strcmpi(vararginx{a+1},methlist);
                if ~any(methidx)
                    error('Regression/Permutation method "%s" unknown.',vararginx{a+1});
                else
                    a = a + 2;
                end
                opts.rmethod = methlist{methidx};
            else
                error([...
                    'The option -rmethod requires a method to be specified.\n'...
                    'Consult the documentation.%s'],'');
            end
            
        case '-npc' % basic
            
            % This is a shortcut to enable NPC with the default settings.
            opts.NPC = true;
            opts.npcmod = true;
            a = a + 1;
            
        case '-npcmethod' % basic
            
            % Do the non-parametric combination?
            if nargin == a || (nargin > a && strcmp(vararginx{a+1}(1),'-'))
                error('The option "-npcmethod" requires a combining method to be indicated.');
                
            elseif nargin > a
                
                % Which combining function to use for the combination?
                methlist = {               ...
                    'Tippett',             ...
                    'Fisher',              ...
                    'Pearson-David',       ...
                    'Stouffer',            ...
                    'Wilkinson',           ...
                    'Winer',               ...
                    'Edgington',           ...
                    'Mudholkar-George',    ...
                    'Friston',             ...
                    'Darlington-Hayes',    ...
                    'Zaykin',              ...
                    'Dudbridge-Koeleman',  ...
                    'Dudbridge-Koeleman2', ...
                    'Taylor-Tibshirani',   ...
                    'Jiang'};
                methidx = strcmpi(vararginx{a+1},methlist);
                
                % Check if method exists, and load extra parameters if needed
                if ~any(methidx)
                    error('Combining method "%s" unknown.',vararginx{a+1});
                elseif any(strcmpi(vararginx{a+1},{...
                        'Wilkinson',       ...
                        'Zaykin',          ...
                        'Jiang'}))
                    if ischar(vararginx{a+2}) && ...
                            strcmpi(vararginx{a+2}(1),'-')
                        plm.npcparm = 0.05;
                        a = a + 2;
                    elseif ischar(vararginx{a+2})
                        a = a + 3;
                        plm.npcparm = eval(vararginx{a+2});
                    else
                        plm.npcparm = vararginx{a+2};
                        a = a + 3;
                    end
                elseif any(strcmpi(vararginx{a+1},{...
                        'Darlington-Hayes',   ...
                        'Dudbridge-Koeleman', ...
                        'Jiang'}))
                    if nargin == a + 1 || ...
                            ischar(vararginx{a+2}) && ...
                            strcmpi(vararginx{a+2}(1),'-')
                        plm.npcparm = 1;
                        a = a + 2;
                    elseif ischar(vararginx{a+2})
                        plm.npcparm = eval(vararginx{a+2});
                        a = a + 3;
                    else
                        plm.npcparm = vararginx{a+2};
                        a = a + 3;
                    end
                elseif strcmpi(vararginx{a+1},'Friston')
                    if nargin == a + 1 || ...
                            ischar(vararginx{a+2}) && ...
                            strcmpi(vararginx{a+2}(1),'-')
                        plm.npcparm = 1;
                        a = a + 2;
                    elseif ischar(vararginx{a+2})
                        plm.npcparm = eval(vararginx{a+2});
                        a = a + 3;
                    else
                        plm.npcparm = vararginx{a+2};
                        a = a + 3;
                    end
                elseif strcmpi(vararginx{a+1},'Dudbridge-Koeleman2')
                    if nargin == a + 1 || ...
                            ischar(vararginx{a+2}) && ...
                            strcmpi(vararginx{a+2}(1),'-')
                        plm.npcparm  = 1;
                        plm.npcparm2 = 0.05;
                        a = a + 2;
                    else
                        if ischar(vararginx{a+2})
                            plm.npcparm = eval(vararginx{a+2});
                        else
                            plm.npcparm = vararginx{a+2};
                        end
                        if nargin == a + 2 || ...
                                ischar(vararginx{a+3}) && ...
                                strcmpi(vararginx{a+3}(1),'-')
                            plm.npcparm2 = 0.05;
                        elseif ischar(vararginx{a+3})
                            plm.npcparm2 = eval(vararginx{a+3});
                        else
                            plm.npcparm2 = vararginx{a+3};
                        end
                        a = a + 4;
                    end
                else
                    a = a + 2;
                end
                opts.npcmethod = methlist{methidx};
            end
            
        case '-npcmod' % basic
            
            % NPC over modalities.
            opts.NPC    = true;
            opts.npcmod = true;
            a = a + 1;
            
        case '-npccon' % basic
            
            % NPC over contrasts -- that is, all contrasts, even contrasts
            % in different designs (if more than one -d is supplied).
            opts.NPC       = true;
            opts.npccon    = true;
            opts.syncperms = true;
            a = a + 1;
            
        case '-mv' % basic
            
            % Compute classic multivariate statistics
            if nargin == a
                opts.MV = true;
                a = a + 1;
                
            elseif nargin > a && strcmp(vararginx{a+1}(1),'-')
                opts.MV = true;
                a = a + 1;
                
            elseif nargin > a
                
                % Which multivariate statistic to use?
                methlist = {            ...
                    'auto',             ...
                    'HotellingTsq',     ...
                    'Wilks',            ...
                    'Lawley',           ...
                    'Lawley-Hotelling', ...
                    'Pillai',           ...
                    'Roy',              ...
                    'Roy-ii',           ...
                    'Roy-iii',          ...
                    'CCA',              ...
                    'PLS'};
                methidx = strcmpi(vararginx{a+1},methlist);
                
                % Check if method exists, and load extra parameters if needed
                if ~any(methidx)
                    error('Multivariate statistic "%s" unknown.',vararginx{a+1});
                elseif strcmpi(vararginx{a+1},'CCA')
                    opts.MV  = false;
                    opts.CCA = true;
                    opts.PLS = false;
                    if nargin == a + 1 || ...
                            ischar(vararginx{a+2}) && ...
                            strcmpi(vararginx{a+2}(1),'-')
                        opts.ccaorplsparm = 1;
                        a = a + 2;
                    elseif ischar(vararginx{a+2})
                        opts.ccaorplsparm = eval(vararginx{a+2});
                        a = a + 3;
                    else
                        opts.ccaorplsparm = vararginx{a+2};
                        a = a + 3;
                    end
                elseif strcmpi(vararginx{a+1},'PLS')
                    opts.MV  = false;
                    opts.CCA = false;
                    opts.PLS = true;
                    if nargin == a + 1 || ...
                            ischar(vararginx{a+2}) && ...
                            strcmpi(vararginx{a+2}(1),'-')
                        opts.ccaorplsparm = 1;
                        a = a + 2;
                    elseif ischar(vararginx{a+2})
                        opts.ccaorplsparm = eval(vararginx{a+2});
                        a = a + 3;
                    else
                        opts.ccaorplsparm = vararginx{a+2};
                        a = a + 3;
                    end
                else
                    opts.MV  = true;
                    opts.CCA = false;
                    opts.PLS = false;
                    a = a + 2;
                end
                opts.mvstat = methlist{methidx};
            end
            
        case '-fdr' % basic
            
            % Compute FDR-adjusted p-values
            opts.FDR = true;
            a = a + 1;
            
        case {'-accel','-approx'} % advanced
            
            % Choose a method to do the approximation of p-values
            if narginx > a && ~strcmpi(vararginx{a+1}(1),'-')
                methlist = {   ...
                    'negbin',  ...
                    'tail',    ...
                    'noperm',  ...
                    'gamma',   ...
                    'lowrank'};
                methidx = strcmpi(vararginx{a+1},methlist);
                if ~ any(methidx)
                    error('Approximation method "%s" unknown.',vararginx{a+1});
                end
                for mm = 1:numel(methlist)
                    opts.accel.(methlist{mm}) = methidx(mm);
                end
                
                % Extra parameters
                if opts.accel.negbin
                    
                    % Number of exceedances:
                    if narginx > a+1 && ~strcmpi(vararginx{a+2}(1),'-')
                        if ischar(vararginx{a+2})
                            opts.accel.negbin = str2double(vararginx{a+2});
                        else
                            opts.accel.negbin = vararginx{a+2};
                        end
                        a = a + 3;
                    else
                        opts.accel.negbin = opts.accel.negbin_nexced;
                        a = a + 2;
                    end
                    
                elseif opts.accel.tail
                    
                    % Define whether include or not the unpermuted stat:
                    if narginx > a+1 && ~strcmpi(vararginx{a+2}(1),'-')
                        if ischar(vararginx{a+2})
                            if     any(strcmpi(vararginx{a+2},{'out','G1out','T1out','true', '1'}))
                                opts.accel.G1out = true;
                                opts.saveuncorrected = false; % defensive, as the uncorrected will be invalid here.
                            elseif any(strcmpi(vararginx{a+2},{'in', 'G1in', 'T1in', 'false','0'}))
                                opts.accel.G1out = false;
                            end
                        else
                            if vararginx{a+2}
                                opts.accel.G1out = true;
                                opts.saveuncorrected = false; % defensive, as the uncorrected will be invalid here.
                            else
                                opts.accel.G1out = false;
                            end
                        end
                        a = a + 3;
                    else
                        a = a + 2;
                    end
                    
                elseif opts.accel.gamma
                    
                    % Define whether include or not the unpermuted stat:
                    if narginx > a+1 && ~strcmpi(vararginx{a+2}(1),'-')
                        if ischar(vararginx{a+2})
                            if     any(strcmpi(vararginx{a+2},{'out','G1out','T1out','true', '1'}))
                                opts.accel.G1out = true;
                                opts.saveuncorrected = false;
                            elseif any(strcmpi(vararginx{a+2},{'in', 'G1in', 'T1in', 'false','0'}))
                                opts.accel.G1out = false; % defensive, as the uncorrected will be invalid here.
                            end
                        else
                            if vararginx{a+2}
                                opts.accel.G1out = true;
                                opts.saveuncorrected = false; % defensive, as the uncorrected will be invalid here.
                            else
                                opts.accel.G1out = false;
                            end
                        end
                        a = a + 3;
                    else
                        a = a + 2;
                    end
                    
                elseif opts.accel.lowrank
                    
                    % Fraction of voxels to be sampled (if < 1) or actual
                    % number of voxels to be sampled.
                    if narginx > a+1 && ~strcmpi(vararginx{a+2}(1),'-')
                        if ischar(vararginx{a+2})
                            opts.accel.lowrank_val = str2double(vararginx{a+2});
                        else
                            opts.accel.lowrank_val = vararginx{a+2};
                        end
                        a = a + 3;
                    else
                        a = a + 2;
                    end
                    
                else
                    a = a + 2;
                end
            else
                error([...
                    'The options "-accel" and "-approx" require a method to.\n' ...
                    'be specified. Consult the documentation.%s'],'');
            end
            
        case {'-noniiclass','-nonifticlass'} % advanced
            
            % Disable using the NIFTI class
            opts.useniiclass = false;
            a = a + 1;
            
        case '-precision' % advanced
            
            % Precision to use?
            if narginx > a && ~strcmpi(vararginx{a+1}(1),'-')
                methlist = {'single','double'};
                methidx = strcmpi(vararginx{a+1},methlist);
                if ~any(methidx)
                    error('Precision "%s" unknown. Use "single" or "double".',vararginx{a+1});
                else
                    a = a + 2;
                end
                opts.precision = methlist{methidx};
            else
                error([...
                    'The option "-precision" requires a method to be specified.\n'...
                    'Use "-precision double" or "-precision single".']);
            end
            
        case '-saveperms' % advanced
            
            % Save the permutations
            opts.saveperms = true;
            a = a + 1;
            
        case '-savemax' % advanced
            
            % Save the permutations
            opts.savemax = true;
            a = a + 1;
            
        case '-savemetrics' % advanced
            
            % Save a file with the number of permutations, average
            % Hamming distance, etc.
            opts.savemetrics = true;
            a = a + 1;
            
        case '-inormal' % advanced
            
            % Inverse-normal transformation?
            opts.inormal = true;
            
            % Take the parameters given to -inormal
            parms = {};
            if narginx - a >= 1
                if ~strcmp(vararginx{a+1}(1),'-')
                    parms{1} = vararginx{a+1};
                end
            end
            if narginx - a >= 2
                if ~strcmp(vararginx{a+2}(1),'-')
                    parms{2} = vararginx{a+2};
                end
            end
            a = a + 1 + numel(parms);
            
            % Then modify the variables accordingly
            methlist = {   ...
                'Blom',    ...
                'Tukey',   ...
                'Bliss',   ...
                'Waerden', ...
                'SOLAR'};
            for p = 1:numel(parms)
                methidx = strcmpi(parms{p},methlist);
                if any(methidx)
                    opts.inormal_meth = parms{p};
                elseif any(parms{p},{'quali','qualitative','discrete'})
                    opts.inormal_quanti = false;
                elseif any(parms{p},{'quanti','quantitative','continuous'})
                    opts.inormal_quanti = true;
                else
                    error('Parameter "%s" unknown for the "-inormal" option.',parms{p});
                end
            end
            
        case '-probit' % advanced
            
            % Probit transformation?
            opts.probit = true;
            a = a + 1;
            
        case '-seed' % advanced
            
            % Seed for the random number generator
            opts.seed = vararginx{a+1};
            if ischar(opts.seed) && ...
                    ~any(strcmpi(opts.seed,{'shuffle','twist','reset'}))
                opts.seed = str2double(opts.seed);
            end
            a = a + 2;
            
        case '-demean' % basic
            
            % Demean data and design. Additionally, remove
            % a global intercept, if any, from the design.
            opts.demean = true;
            a = a + 1;
            
        case '-vgdemean' % advanced
            
            % Demean data and design within VG. Additionally, remove
            % a global intercept, if any, from the design.
            opts.vgdemean = true;
            a = a + 1;
            
        case '-ev4vg' % advanced
            
            % Add to the design matrix one EV for each variance group.
            opts.ev4vg = true;
            a = a + 1;
            
        case '-removevgbysize' % advanced
            
            % Remove from the analysis observations that are the only
            % in their variance group.
            opts.removevgbysize = vararginx{a+1};
            if ischar(opts.removevgbysize)
                opts.removevgbysize = str2double(opts.removevgbysize);
            end
            a = a + 2;
            
        case '-zstat' % advanced
            
            % Convert the statistic for each test to a z-score
            opts.zstat = true;
            a = a + 1;
            
        case '-pearson' % basic
            
            % Compute the Pearson's correlation coefficient (R^2 if rank(C)>1).
            opts.pearson = true;
            a = a + 1;
            
        case '-noranktest' % advanced
            
            % Don't test the rank(Y) before doing MANOVA/MANCOVA.
            opts.noranktest = true;
            a = a + 1;
            
        case '-transposedata' % advanced
            
            % Transpose the data if it's 2D?
            opts.transposedata = true;
            a = a + 1;
            
        case '-inputmv' % advanced
            
            % Treat the (sole) input as multivariate, that is,
            % each column is a variable in a multivariate model,
            % as opposed to independent univariate tests.
            % Useful with non-imaging data.
            opts.inputmv = true;
            a = a + 1;
            
        case '-verbosefilenames' % advanced
            
            % Don't simplify filenames when otherwise it would be possible
            opts.verbosefilenames = true;
            a = a + 1;
            
        case '-syncperms' % advanced
            
            % Synchronize permutations regardless of other options?
            opts.syncperms = true;
            a = a + 1;
            
        case {'-designperinput','-erie'} % advanced
            
            % Use one design matrix for each input dataset?
            % This enables
            % syncP regardless.
            opts.designperinput = true;
            opts.syncperms      = true;
            a = a + 1;
            
        case '-subjidx' % advanced
            
            % Indices of the subjects to keep in the design
            opts.subjidx = vararginx{a+1};
            a = a + 2;
            
        case '-quiet' % basic
            
            % Don't print lines indicating progress
            opts.showprogress = false;
            a = a + 1;
            
        case '-nounivariate' % advanced
            
            % Save or not univariate tests? Useful with MV/NPC/CCA
            opts.saveunivariate = false;
            a = a + 1;
            
        case '-nouncorrected'  % advanced
            
            % Save or not uncorrected p-vals
            opts.saveuncorrected = false;
            a = a + 1;
            
        case '-saveuncorrected'  % advanced
            
            % Save or not uncorrected p-vals
            opts.saveuncorrected = true;
            a = a + 1;
            
        case '-savedof' % advanced
            
            % Save a file with the degrees of freedom?
            opts.savedof = true;
            a = a + 1;
            
        case '-pmethodp' % advanced
            
            % Which method to use to partition the model when defining
            % the permutations?
            if narginx > a && ~strcmpi(vararginx{a+1}(1),'-')
                methlist = {    ...
                    'Guttman',  ...
                    'Beckmann', ...
                    'Ridgway',  ...
                    'none'};
                methidx = strcmpi(vararginx{a+1},methlist);
                if ~any(methidx)
                    error('Partition method "%s" unknown.',vararginx{a+1});
                else
                    a = a + 2;
                end
                opts.pmethodp = methlist{methidx};
            else
                error([...
                    'The option "-pmethodp" requires a method to be specified.\n'...
                    'Consult the documentation.']);
            end
            
        case '-pmethodr' % advanced
            
            % Which method to use to partition the model when defining
            % doing the actual regression?
            if narginx > a && ~strcmpi(vararginx{a+1}(1),'-')
                methlist = {    ...
                    'Guttman',  ...
                    'Beckmann', ...
                    'Ridgway',  ...
                    'none'};
                methidx = strcmpi(vararginx{a+1},methlist);
                if ~any(methidx)
                    error('Partition method "%s" unknown.',vararginx{a+1});
                else
                    a = a + 2;
                end
                opts.pmethodr = methlist{methidx};
            else
                error([...
                    'The option "-pmethodr" requires a method to be specified.\n' ...
                    'Consult the documentation.']);
            end
            
        case '-forceintersectionmask' % currently not in the help
            
            % Force the generation of an intersection mask, even if there
            % is no MV or NPC (useful for mediation analysis).
            opts.forcemaskinter = true;
            a = a + 1;
            
        otherwise
            error('Unknown option: "%s"',vararginx{a});
    end
end

% Check for the existence of other programs for input/output
palm_checkprogs;

if Ni == 0
    error('Missing input data (missing "-i").');
end

% Make sure the NPC options make sense
% - if -npc only, it's -npcmod
% - if -npcmod only, it's -npcmod
% - if -npccon only, it's -npccon
% - if -npcmod and -npccon, it's both
if opts.NPC && ~ opts.npccon
    opts.npcmod = true;
end

% A quick check for the case of 1 EV per column of Y.
if opts.evperdat
    if any([...
            opts.ev4vg
            opts.pearson])
        error([...
            'The option "-evperdat" is incompatible with the options listed below:\n' ...
            '"-ev4vg"\n' ...
            '"-pearson"%s'],'');
    end
    if strcmpi(opts.rmethod,'terBraak')
        error('The option "-evperdat" cannot be used with the ter Braak method (not implemented)');
    end
end

% No spatial statistics if no univariate results will be saved anyway
if ~ opts.saveunivariate
    opts.cluster.uni.do = false;
    opts.tfce.uni.do    = false;
end

% This simplifies some tests later
opts.spatial.do  = false;
opts.spatial.uni = false;
opts.spatial.npc = false;
opts.spatial.mv  = false;
if any([ ...
        opts.cluster.uni.do  ...
        opts.cluster.npc.do  ...
        opts.cluster.mv.do   ...
        opts.tfce.uni.do     ...
        opts.tfce.npc.do     ...
        opts.tfce.mv.do]')
    opts.spatial.do = true;
    if any([ ...
            opts.cluster.uni.do  ...
            opts.tfce.uni.do]')
        opts.spatial.uni = true;
    end
    if any([ ...
            opts.cluster.npc.do  ...
            opts.tfce.npc.do]')
        opts.spatial.npc = true;
    end
    if any([ ...
            opts.cluster.mv.do   ...
            opts.tfce.mv.do]')
        opts.spatial.mv = true;
    end
end

% Adjust the delta-h.
if any([ ...
        opts.tfce.uni.do
        opts.tfce.npc.do
        opts.tfce.mv.do ]) && ...
        strcmpi(opts.tfce.deltah,'auto')
    opts.tfce.deltah = 0;
end

% Sanity checks for the acceleration modes.
if sum(logical([ ...
        opts.accel.negbin, ...
        opts.accel.tail,   ...
        opts.accel.noperm, ...
        opts.accel.gamma,  ...
        opts.accel.lowrank])) > 1
    error('Only one approximation method can be used for a given run.');
end
if opts.accel.negbin
    if opts.accel.negbin < 2
        error('The parameter r given to "-accel negbin <r>" must be >= 2.')
    end
    if opts.nP0 < 3
        error('The option "-accel negbin <r>" needs at least 3 permutations.')
    end
    if ~ opts.saveuncorrected
        error('The option "-nouncorrected" cannot be used with "-accel negbin".');
    end
    if (opts.corrmod || opts.corrcon) && ~ opts.FDR
        error('The option "-accel negbin" cannot be used with FWER-correction, only FDR.');
    end
    if opts.NPC
        error('The option "-accel negbin" cannot be used with NPC.');
    end
    if opts.spatial.do
        error('The option "-accel negbin" cannot be used with spatial statistics.');
    end
    if opts.saveperms
        error('The option "-saveperms" cannot be used together with "-accel negbin".');
    end
end
if opts.accel.noperm
    if ~ opts.saveuncorrected
        error('The option "-nouncorrected" cannot be used with "-accel noperm".');
    end
    if (opts.corrmod || opts.corrcon) && ~ opts.FDR
        error('The option "-accel noperm" cannot be used with FWER-correction, only FDR.');
    end
    if opts.NPC
        error('The option "-accel noperm" cannot be used with NPC.');
    end
    if opts.spatial.do
        error('The option "-accel noperm" cannot be used with spatial statistics.');
    end
    if opts.MV
        if ~ any(strcmpi(opts.mvstat,{'auto','pillai'}))
            warning([...
                'With multivariate tests, the option "-accel noperm" can\n' ...
                '         only be used with the Pillai'' trace statistic.\n' ...
                '         Changing automatically to Pillai.%s'],'');
        end
        opts.mvstat = 'pillai';
    end
    if Ni > 1 && ~ opts.MV
        error([...
            'The option "-accel noperm" needs to be used with a single modality\n'...
            '       modality or with "-mv".%s'],'');
    end
end
if opts.accel.lowrank
    if opts.pearson
        error('The option "-accel lowrank" cannot be used with "-pearson".');
    end
    if opts.MV
        error('The option "-accel lowrank" cannot be used with MV.');
    end
    if opts.NPC
        error('The option "-accel lowrank" cannot be used with NPC.');
    end
    if opts.spatial.do
        error('The option "-accel lowrank" cannot be used with spatial statistics.');
    end
    if opts.evperdat
        error('The option "-accel lowrank" cannot be used with "-evperdat".');
    end
    if opts.missingdata
        error('The option "-accel lowrank" cannot be used with missing data.');
    end
    if opts.saveglm
        error('The option "-accel lowrank" cannot be used with "-saveglm".');
    end
end

% Some options can't be used with missing data
if Nimiss || Ndmiss
    opts.missingdata = true;
    if opts.MV
        error('The option "-mv" cannot be used with missing data.');
    end
    if ~ opts.zstat && ~ opts.mcar
        warning([...
            'With missing data MAR/MNAR, the option "-zstat" is mandatory.\n' ...
            '         Adding it automatically.%s'],'');
        opts.zstat = true;
    end
    if ~ opts.cmcx && ~ opts.mcar
        warning([...
            'With missing data, the option "-cmcx" is mandatory.\n' ...
            '         Adding it automatically.%s'],'');
        opts.cmcx = true;
    end
    if opts.demean && ~ opts.mcar
        warning([...
            'With missing data, the option "-demean" must not be used.\n' ...
            '         Removing it automatically.%s'],'');
        opts.demean = false;
    end
    if ~ strcmpi(opts.pmethodp,'guttman') || ~ strcmpi(opts.pmethodr,'guttman')
        warning([...
            'With missing data, the partitioning must use the "Guttman".\n' ...
            '         method. Adding automatically the options\n' ...
            '         "-pmethodp Guttman" and "-pmethodr Guttman".%s'],'');
        opts.pmethodp = 'guttman';
        opts.pmethodr = 'guttman'; 
    end
    if opts.ev4vg || opts.removevgbysize > 0
        error('Missing data cannot be used with "-ev4vg" or "-removevgbysize".')
    end
end

% Some NPC methods don't have an analytical form for the parametric p-value
if opts.NPC && any(strcmpi(opts.npcmethod,{'Darlington-Hayes','Jiang'}))
    plm.nonpcppara = true;
    if opts.savepara
        warning([...
            'No parametric combination p-value will be saved for the\n', ...
            '         Darlington-Hayes or Jiang methods%s'],'');
    end
    if opts.spatial.npc
        warning([ ...
            'No NPC spatial statistic will be produced for the\n', ...
            '         Darlington-Hayes or Jiang methods%s'],'');
        opts.cluster.npc.do = false;
        opts.tfce.npc.do    = false;
        opts.spatial.npc    = false;
    end
else
    plm.nonpcppara = false;
end

% Likewise, some MV methods don't have an analytical form for the parametric p-value
if opts.MV && strcmpi(opts.mvstat,'Roy_iii')
    plm.nomvppara = true;
    if opts.savepara
        warning('No parametric p-value will be saved for the Roy_iii method.%s','');
    end
    if opts.spatial.mv
        warning([ ...
            'No multivariate cluster-level or TFCE statistic will be produced\n', ...
            '         for the Roy_iii statistic%s'],'');
        opts.cluster.mv.do = false;
        opts.tfce.mv.do    = false;
        opts.spatial.mv    = false;
    end
elseif (opts.CCA || opts.PLS) && opts.savepara
    plm.nomvppara = true;
    warning([...
        'No parametric p-value will be saved for CCA or PLS.\n', ...
        '         Use Wilks'' lambda instead.%s'],'');
else
    plm.nomvppara = false;
end

% Some more warnings and sanity checks
if opts.savecdf && opts.savelogp
    error('Should not use "-save1-p" together with "-logp"');
end
if ~ opts.inputmv && opts.designperinput && Ni ~= Nd
    error([
        'To use the option "-designperinput", the number of design files must\n' ...
        'match the number of inputs.\n%s'],'');
end
if (Nt || Nf) && Ncon
    error('Cannot mix options "-t" or "-f" with "-con".');
end
if Nt || Nf
    if Nt > Nd
        error('More t-contrast files (%d) than valid design files (%d) were supplied.',Nt,Nd);
    end
    if Nt ~= 1 && Nt ~= Nd
        error(['The number of supplied t-contrast files (option "-t") must be 1 or\n',...
            'the same as the number of design files (option "-d") (%d).'],Nd);
    end
    if Nf > Nt
        error('More F-contrast files (%d) than t-contrast files (%d) were supplied.',Nf,Nt);
    end
elseif Ncon
    if Ncon > Nd
        error('More contrast files (%d) than design files (%d) were supplied.',Nt,Nd);
    end
    if Ncon ~= 1 && Ncon ~= Nd
        error(['The number of supplied contrast files (option "-con") must be 1 or\n',...
            'the same as the number of design files (option "-d") (%d).'],Nd);
    end
end
if opts.pearson && (opts.NPC || opts.MV)
    error([ ...
        'It isn''t possible to compute the Pearson''s r or R^2 together with NPC or\n', ...
        '         multivariate methods.%s'],'');
end
if opts.pearson && ~ opts.demean
    warning([ ...
        'To compute Pearson''s "r" or the "R^2", the data and the design\n' ...
        '         must be mean centered. Adding option "-demean".%s'],'');
    opts.demean = true;
end
if opts.pearson && opts.savepara
    error([ ...
        'The option "-saveparametric" cannot be used with "-pearson".\n' ...
        '         For a parametric p-value, drop "-pearson" and use the\n' ...
        '         respective results for the t or F-statistics.%s'],'');
end
if opts.pearson && ~ any(strcmpi(opts.pmethodr,{'beckmann','ridgway'}))
    warning([ ...
        'To compute Pearson''s "r" or the "R^2", the design must be\n' ...
        '         partitioned using the Beckmann or Ridgway schemes.'...
        '         Adding the option "-pmethodr Beckmann".%s'],'');
    opts.pmethodr = 'beckmann';
end
if (opts.CCA || opts.PLS) && ~ opts.demean
    warning([ ...
        'To perform CCA or PLS, the data and the design\n' ...
        '         must be mean centered. Adding option "-demean".%s'],'');
    opts.demean = true;
end
if (opts.CCA || opts.PLS) && ~ any(strcmpi(opts.pmethodr,{'beckmann','ridgway'}))
    warning([ ...
        'To perform CCA or PLS, the design must be\n' ...
        '         partitioned using the Beckmann or Ridgway schemes.'...
        '         Adding the option "-pmethodr Beckmann".%s'],'');
    opts.pmethodr = 'beckmann';
end
if opts.demean && opts.vgdemean && ~opts.pearson && ~opts.CCA && ~opts.PLS
    warning([...
        'Cannot use the option "-demean" together with "-vgdemean"\n'...
        '         Ignoring the option "-vgdemean".%s'],'');
    opts.vgdemean = false;
end
if opts.ev4vg && opts.vgdemean
    error('Cannot use the option "-ev4vg" together with "-vgdemean"');
end
if opts.MV && ~ opts.singlevg
    error('Currently MV is only allowed with a single VG, that is, assuming homoscedasticity.');
end
if opts.designperinput && opts.MV
    error('It''s not possible to use the option "-designperinput" together with the option "-mv".');
end
if ~opts.saveunivariate && ~opts.MV && ~opts.NPC && ~opts.CCA && ~opts.PLS
    error('The option "-nounivariate" can only be used with "-mv" or "-npc".');
end
if opts.MV && (opts.CCA || opts.PLS)
    error('Cannot do classical MANOVA at the same time as CCA or PLS.');
end
if Ni > 1 && opts.inputmv
    error('Option "-inputmv" cannot be used with more than one "-i".')
end
if opts.concordant && ~ opts.NPC
    error('The option "-concordant" is for use with NPC only.');
end
if opts.concordant && opts.twotail
    error(['Cannot use "-concordant" together with "-twotail" (inadmissible).\n'...
        'Use either of these, but not both together.%s'],'');
end
if opts.tonly && opts.fonly
    error('Cannot use "-tonly" together with "-fonly".');
end
if opts.probit && opts.inormal
    error('Cannot use "-probit" together with "-inormal".');
end
if opts.FDR && ~ opts.saveuncorrected
    error(['Option "-fdr" cannot be used together with "-nouncorrected".\n'...
        'Use either of these, but not both together. If you are using tail or\n'...
        'gamma approximations, consider keeping "-nouncorrected", and leave\n'...
        '"-fdr" for a separate call in which tail or gamma are not used.%s'],'');
end

% Initialize the random number generator (if nP = 0, no need for that)
if opts.nP0
    if palm_isoctave
        if any(strcmpi(opts.seed,{'reset','shuffle','twist'}))
            opts.seed = 'reset';
        end
        rand('state',opts.seed); %#ok
    else
        if any(strcmpi(opts.seed,{'reset','shuffle','twist'}))
            opts.seed = 'shuffle';
        end
        rng('default');
        rng(opts.seed);
    end
end

% Read and organise the surfaces. If no surfaces have been loaded, but the
% user wants cluster extent, cluster mass, or TFCE, an error will be
% printed later down in the code.
% At this stage the average area from native geometry is also loaded if it
% was provided. If a weight was given, such as 1, use this weight, so that
% all faces are treated as if having the same size. If nothing was
% supplied, use the actual area of the surface.
if opts.spatial.do && Ns > 0
    plm.srf     = cell(Ns,1);
    plm.srfarea = cell(Ns,1);
    plm.srfadj  = cell(Ns,1);
    for s = 1:Ns
        
        % Load surface
        fprintf('Loading surface %d/%d: %s\n',s,Ns,opts.s{s});
        plm.srf{s} = palm_miscread(opts.s{s},[],[],[],true);
        
        % Load areas
        if isempty(opts.sa{s})
            % No area means area of the actual surface file
            plm.srfarea{s}.data = [];
        elseif exist(opts.sa{s},'file')
            % A file with the average areas from native geometry
            plm.srfarea{s} = palm_miscread(opts.sa{s},opts.useniiclass,opts.o,opts.precision,false);
        elseif ~ isnan(str2double(opts.sa{s}))
            % A weight (such as 1)
            plm.srfarea{s}.data = str2double(opts.sa{s});
        else
            % Otherwise gives a helpful error message:
            error('Unknown option given to "-s" or file doesn''t exist:\n%s',opts.sa{s});
        end
    end
end

% There should be no more masks than modalities, and the number of
% masks needs to be either 1 or the same number of modalities.
if Nm > Ni
    error([...
        'There are more masks supplied with "-m" (%d masks) than\n'...
        'modalities supplied with "-i" (%d modalities)'],Nm,Ni);
elseif Nm > 1 && Nm ~= Ni
    error([...
        'The number of masks supplied with "-m" (%d masks) is larger than 1,\n'...
        'but still not the same as the number of modalities supplied with\n'...
        'the option "-i" (%d modalities).'],Nm,Ni);
end

% Read and organise the masks. If there are no masks specified, one for
% each modality will be created after each modality is loaded.
plm.masks = cell(Ni,1);
for m = 1:Nm
    plm.masks{m} = palm_miscread(opts.m{m},opts.useniiclass,opts.o,opts.precision,false);
    if strcmp(plm.masks{m}.readwith,'nifticlass')
        plm.masks{m}.data = double(plm.masks{m}.data);
    end
    if opts.reversemasks
        plm.masks{m}.data(isnan(plm.masks{m}.data)) = 1;
        plm.masks{m}.data(isinf(plm.masks{m}.data)) = 1;
        plm.masks{m}.data = ~ logical(plm.masks{m}.data);
    else
        plm.masks{m}.data(isnan(plm.masks{m}.data)) = 0;
        plm.masks{m}.data(isinf(plm.masks{m}.data)) = 0;
        plm.masks{m}.data = logical(plm.masks{m}.data);
    end
end
if Nm == 1
    for i = 2:Ni
        plm.masks{i} = plm.masks{1};
    end
end

% Indices of the subjects that will be kept
if ~ isempty(opts.subjidx)
    plm.subjidx = palm_miscread(opts.subjidx);
    plm.subjidx = round(plm.subjidx.data);
    u = unique(plm.subjidx);
    if numel(u) == 1
        error('Subject indices are all identical in the -subjidx file.');
    end
end

% Read and organise the data.
plm.Yset     = cell(Ni,1);  % Regressands (Y)
plm.Yisvol   = false(Ni,1); % Is Y a volume image?
plm.Yissrf   = false(Ni,1); % Is Y a surface-based image (DPX)?
plm.Yisvtx   = false(Ns,1); % Is vertexwise?
plm.Yisfac   = false(Ns,1); % is facewise? (this is currently dichotomous with Yisvtx, but later there may be edges/connectivity too)
plm.Yarea    = cell(Ns,1);  % To store area per face or per vertex (used for cluster-level & TFCE).
plm.Ykindstr = cell(Ni,1);  % string to save the files later
for i = 1:Ni
    
    % Read input file
    fprintf('Reading input %d/%d: %s\n',i,Ni,opts.i{i});
    if Nm == 0
        maskstruct = [];
    else
        maskstruct = plm.masks{i};
    end
    [plm.Yset{i},plm.masks{i},plm.Yisvol(i),plm.Yissrf(i),plm.Ykindstr{i}] = ...
        palm_ready(opts.i{i},maskstruct,opts,true);
    
    % Select subjects
    if ~ isempty(plm.subjidx)
        plm.Yset{i} = plm.Yset{i}(plm.subjidx,:);
    end
    
    % For the first input data, keep the size to
    % compare with the others, then check the size
    if i == 1
        plm.N = size(plm.Yset{i},1);
    end
    if size(plm.Yset{i},1) ~= plm.N
        error([
            'At least two of the input data files do not have\n' ...
            'compatible sizes:\n' ...
            '- File %d (%s) has %d observations\n' ...
            '- File %d (%s) has %d observations'], ...
            1,opts.i{1},plm.N, ...
            i,opts.i{i},size(plm.Yset{i},1));
    end
    
    % If this is a DPX/curvature file, and if one of the spatial
    % statistics has been invoked, check if surfaces are available
    % and with compatible size, then compute the area (dpv or dpf).
    % Also take this opportunity to compute the adjacency matrix.
    if opts.tableasvolume
        plm.Yisvol(i) = true;
    else
        if opts.spatial.do && Ns > 0
            
            
            % String defining the types, for the filenames and other tasks.
            if Ns == 1, s = 1; else, s = i; end
            %Changed John Williams JCW 04/03/2023
            thisYisVol = plm.Yisvol(i);
            if(~thisYisVol)
                if any(size(plm.srf{s}.data.vtx,1) == ...
                        size(plm.masks{i}.data))
                    plm.Yisvol(i)   = false;
                    plm.Yissrf(i)   = true;
                    plm.Yisvtx(i)   = true;
                    plm.Yisfac(i)   = false;
                    plm.Ykindstr{i} = '_dpv';
                elseif any(size(plm.srf{s}.data.fac,1) == ...
                        size(plm.masks{i}.data))
                    plm.Yisvol(i)   = false;
                    plm.Yissrf(i)   = true;
                    plm.Yisvtx(i)   = false;
                    plm.Yisfac(i)   = true;
                    plm.Ykindstr{i} = '_dpf';
                else
                    error([...
                        'Surface file does not match the input data:\n' ...
                        '- Surface file %d has %d vertices and %d faces (%s)\n' ...
                        '- Input data file %d has %d points (%s)'],...
                        s,size(plm.srf{s}.data.vtx,1),size(plm.srf{s}.data.fac,1),opts.s{s},...
                        i,max(size(plm.masks{i}.data)),opts.i{i});
                end

                % Surface area, to be used by the spatial statistics
                if isempty(plm.srfarea{s}.data)

                    % No area file given, use the actual surface area
                    plm.Yarea{i} = palm_calcarea(plm.srf{s}.data,plm.Yisvtx(i));

                elseif numel(plm.srfarea{s}.data) == 1

                    % A weight given (such as 1): use that for each vertex or face,
                    % treating as if all had the same area.
                    if plm.Yisvtx(i)
                        plm.Yarea{i} = plm.srfarea{s}.data .* ...
                            ones(size(plm.srf{s}.data.vtx,1),1);
                    elseif plm.Yisfac(i)
                        plm.Yarea{i} = plm.srfarea{s}.data .* ...
                            ones(size(plm.srf{s}.data.fac,1),1);
                    end

                else

                    % Otherwise, just use the data from the file (already loaded).
                    plm.Yarea{i} = plm.srfarea{s}.data;
                end

                % Compute the adjacency matrix
                plm.Yadjacency{i} = palm_adjacency(plm.srf{s}.data.fac,plm.Yisvtx(i));
            end %End JCW Edit 04/03/2023
            
        elseif opts.spatial.do && Ns == 0 && any(strcmpi(plm.masks{i}.readwith,{'load','fs_load_mgh','gifti'}))
            error([ ...
                'To use spatial statistics with vertexwise or facewise data it is\n'...
                'necessary to provide the surface files (with the option "-s").%s'],'');
        end
    end
end
plm.nmasks = numel(plm.masks);
plm.nY     = numel(plm.Yset); % this is redefined below if opts.inputmv is set.

% Some extra packages for Octave
if palm_isoctave
    if opts.spatial.do && any(plm.Yisvol)
        pkg load image
    end
    if opts.accel.lowrank || opts.zstat || opts.corrcon
        pkg load statistics
    end
end

% Read the EV per datum:
if opts.evperdat
    
    % Some sanity check:
    opts.evpos = cat(1,opts.evpos{:});
    if size(unique(opts.evpos,'rows'),1) ~= size(opts.evpos,1)
        error([
            'Some EV per datum have been defined for the same\n'...
            'position in the same design matrices.%s'],'');
    end
    plm.EVset  = cell(Nevd,1);
    plm.nEVdat = Nevd;
    
    % If there's one design per input, use the same masks as
    % those of the input files. Otherwise, create them.
    if ~ opts.designperinput && Nm == 1
        for ev = 1:plm.nEVdat
            plm.masksEV{ev} = plm.masks{1};
        end
    elseif opts.designperinput || (plm.nY == 1 && Nd == 1)
        for ev = 1:plm.nEVdat
            plm.masksEV{ev} = plm.masks{opts.evpos(ev,2)};
        end
    else
        plm.masksEV = cell(plm.nEVdat,1);
    end
    
    % Read input file & select subjects
    for ev = 1:plm.nEVdat
        fprintf('Reading EV per datum %d/%d: %s\n',ev,plm.nEVdat,opts.evdatfile{ev});
        [plm.EVset{ev},plm.masksEV{ev}] = palm_ready(opts.evdatfile{ev},plm.masksEV{ev},opts,opts.demean);
        if ~ isempty(plm.subjidx)
            plm.EVset{ev} = plm.EVset{ev}(plm.subjidx,:);
        end
    end
    plm.nmasksEV = numel(plm.masksEV);
    
    % Make the intersection of the EVs that will go all in the same design
    Dlist = unique(opts.evpos(:,2),'rows');
    for d = 1:numel(Dlist,1)
        evidx = find(opts.evpos(:,2) == Dlist(d));
        if numel(evidx) > 1
            newmask = true(size(plm.masksEV{1}.data));
            for ev = evidx'
                newmask = newmask & plm.masksEV{ev}.data;
            end
            for ev = evidx'
                plm.EVset{ev} = plm.EVset{ev}(:,newmask(plm.masksEV{ev}.data(:)));
                plm.masksEV{ev}.data = newmask;
            end
        end
    end
    
    % Then make the intersections with the respective masks of the input files
    if (opts.designperinput || (plm.nY == 1 && Nd == 1)) && ...
            ~ opts.npcmod && ~ opts.npccon && ~ opts.MV
        for d = unique(opts.evpos(:,2))'
            EV = find(opts.evpos(:,2) == d);
            newmask = zeros(numel(plm.masks{d}.data),numel(EV)+1);
            c = 1;
            for ev = EV'
                newmask(:,c) = plm.masksEV{ev}.data(:);
                c = c + 1;
            end
            newmask(:,c) = plm.masks{d}.data(:);
            newmask = reshape(all(newmask,2),size(plm.masks{d}.data));
            for ev = EV'
                plm.EVset{ev} = plm.EVset{ev}(:,newmask(plm.masksEV{ev}.data(:)));
                plm.masksEV{ev}.data = newmask;
            end
            plm.Yset{d} = plm.Yset{d}(:,newmask(plm.masks{d}.data(:)));
            plm.masks{d}.data    = newmask;
        end
    else
        sizy = zeros(plm.nmasks,1);
        for y = 1:plm.nmasks
            sizy(y) = numel(plm.masks{y}.data);
        end
        sizev = zeros(plm.nmasksEV,1);
        for ev = 1:plm.nmasksEV
            sizev(ev) = numel(plm.masksEV{ev}.data);
        end
        if  numel(unique(sizy)) > 1 || ...
                numel(unique(sizev)) > 1 || ...
                sizy(1) ~= sizev(1)
            error([...
                'For multiple "-i" and/or "-evperdat", with "-npccon",\n',...
                'and without the option "-designperinput", the inputs, EVs \n'...
                'and masks need to be all of the same sizes.%s'],'');
        end
        newmask = true(size(plm.masksEV{1}.data));
        for y = 1:plm.nmasks
            newmask = newmask & plm.masks{y}.data;
        end
        for ev = 1:plm.nmasksEV
            newmask = newmask & plm.masksEV{ev}.data;
        end
        for y = 1:plm.nmasks
            plm.Yset{y} = plm.Yset{y}(:,newmask(plm.masks{y}.data(:)));
            plm.masks{y}.data = newmask;
        end
        for ev = 1:plm.nmasksEV
            plm.EVset{ev} = plm.EVset{ev}(:,newmask(plm.masksEV{ev}.data(:)));
            plm.masksEV{ev}.data = newmask;
        end
    end
end
clear('newmask');

% Create an intersection mask if NPC or MV is to be done, and further apply
% to the data that was previously masked above, as needed.
if opts.npcmod || opts.MV || opts.CCA || opts.PLS || opts.forcemaskinter
    if plm.nmasks > 1
        
        % If there is one mask per modality, make an instersection mask.
        maskinter = true(size(plm.masks{1}.data));
        for m = 1:plm.nmasks
            maskinter = maskinter & plm.masks{m}.data;
        end
        if opts.evperdat
            for ev = 1:plm.nEVdat
                maskinter = maskinter & plm.masksEV{ev}.data;
            end
        end
        
        % Note that this line below uses Ytmp, which is from the previous loop.
        % This can be used here because with NPC all data has the same size.
        plm.maskinter = palm_maskstruct(maskinter(:)',plm.masks{1}.readwith,plm.masks{1}.extra);
        
        % Apply it to further subselect data points
        for y = 1:plm.nY
            plm.Yset{y} = plm.Yset{y}(:,plm.maskinter.data(plm.masks{y}.data));
        end
        if opts.evperdat
            for ev = 1:plm.nEVdat
                plm.EVset{ev} = plm.EVset{ev}(:,plm.maskinter.data(plm.masksEV{ev}.data));
            end
        end
    else
        
        % If only one mask was given.
        plm.maskinter = plm.masks{1};
        for y = 1:plm.nY
            plm.Yset{y} = plm.Yset{y}(:,plm.maskinter.data(plm.masks{1}.data));
        end
    end
end

% Sizes for later
if opts.evperdat
    plm.EVsiz = zeros(plm.nEVdat,1);
    for ev = 1:plm.nEVdat
        plm.EVsiz(ev) = size(plm.EVset{ev},2);
    end
end

% Make sure that all data have the same size if NPC or MV are to be done
if opts.npcmod || opts.MV || opts.forcemaskinter
    % The plm.Ysiz is redefined below.
    for y = 1:plm.nY
        plm.Ysiz(y) = size(plm.Yset{y},2);
    end
    [usiz,uidx] = unique(plm.Ysiz);
    if numel(usiz) > 2 || (numel(usiz) == 2 && min(usiz) ~= 1)
        error('The sizes of some of the imaging modalities don''t match');
    elseif numel(usiz) == 2 && usiz(1) == 1
        for y = 1:plm.nY
            if plm.Ysiz(y) == 1
                fprintf('Expanding modality #%d to match the size of the others.\n',y);
                plm.Yset{y} = repmat(plm.Yset{y},[1 usiz(2)]);
                if plm.nmasks > 1
                    if numel(plm.masks{y}.data) == 1
                        plm.masks{y}.data = plm.masks{uidx(2)}.data;
                    else
                        error('Modality expansion is only allowed for single input variables.')
                    end
                end
            end
        end
    end
    clear('usiz');
end

% Make sure none of the modalities is empty
for y = 1:plm.nY
    if any(size(plm.Yset{y}) == 0)
        error('Modality %d has no data.\n',y);
    end
end

% If the multiple columns of the (sole) input are to be treated
% in a multivariate fashion
if opts.inputmv
    tmp1 = plm.Yset{1};
    tmp2 = plm.Ykindstr{1};
    nY   = size(tmp1,2);
    plm.Yset     = cell(nY,1);
    plm.Ykindstr = cell(nY,1);
    for y = 1:nY
        plm.Yset{y}     = tmp1(:,y);
        plm.Ykindstr{y} = tmp2;
    end
    clear tmp1 tmp2 nY;
    plm.nY = numel(plm.Yset);
end

% A variable with the cumulative sizes of all modalities will be handy later
plm.Ysiz = zeros(plm.nY,1);
for y = 1:plm.nY
    plm.Ysiz(y) = size(plm.Yset{y},2);
end
plm.Ycumsiz = vertcat(0,cumsum(plm.Ysiz));

% Take this opportunity to save the masks if the user requested.
if opts.savemask
    for y = 1:plm.nmasks
        M = plm.masks{y};
        if plm.nY == 1 || plm.nmasks == 1
            M.filename = sprintf('%smask',opts.o);
        else
            M.filename = sprintf('%smask_m%d',opts.o,y);
        end
        M.data = double(M.data);
        palm_miscwrite(M,true);
    end
    if plm.nY > 1 && (opts.npcmod || opts.MV || opts.forcemaskinter)
        M          = plm.maskinter;
        M.filename = sprintf('%sintersection_mask',opts.o);
        M.data     = double(M.data);
        palm_miscwrite(M,true);
    end
end

% If MV was selected, make sure that Y is full rank.
if opts.MV && ~ opts.noranktest
    fprintf('Testing rank of the data for the MV tests. To skip, use -noranktest.\n')
    Y = cat(3,plm.Yset{:});
    Y = permute(Y,[1 3 2]);
    failed = false(1,size(Y,3));
    for v = 1:size(Y,3)
        if rank(Y(:,:,v)) ~= plm.nY
            failed(v) = true;
        end
    end
    if any(failed)
        fname = sprintf('%smv_illconditioned',opts.o);
        palm_quicksave(double(failed),0,opts,plm,[],[],[],fname);
        error([
            'One or more datapoints have ill-conditioned data. It is\n' ...
            'not possible to run multivariate analyses as MANOVA/MANCOVA.\n' ...
            'Please, see these datapoints marked as 1 in the file:\n' ...
            '%s.*\n'],fname); %#ok
    end
end; clear Y;

% Applies an inverse-normal transformation to the modalities if the user requested
if opts.inormal
    for y = 1:plm.nY
        plm.Yset{y} = palm_inormal( ...
            plm.Yset{y},            ...
            opts.inormal_meth,      ...
            opts.inormal_quanti);
    end
end

% Applies a probit transformation to the modalities if the user requested
if opts.probit
    for y = 1:plm.nY
        if min(plm.Yset{y}(:)) < 0 || max(plm.Yset{y}(:)) > 1
            error([
                'Probit transformation can only be used with data in the interval [0 1].\n' ...
                'This fails for at least modality #%d.'],y);
        end
        plm.Yset{y} = erfinv(2*(plm.Yset{y}*0.999999999999999 + .5e-15)-1)*sqrt(2);
    end
end

% Make the adjustments for the EE and ISE options.
% - if the user gives nothing, its EE by default.
% - if the user gives ISE only, it's ISE only
% - if the user gives EE only, it's EE only
% - if the user gives both, it's both
if ~opts.EE && ~opts.ISE
    opts.EE  = true;
end

% Read and assemble the design matrices.
fprintf('Reading design matrix and contrasts.\n');
if opts.evperdat
    plm.Mset = cell(max(Nd,max(opts.evpos(:,2))),1);
    for m = 1:numel(plm.Mset)
        plm.Mset{m} = zeros(plm.N,1);
    end
else
    plm.Mset = cell(max(Nd,1),1);
end
plm.nM = numel(plm.Mset);
if Nd == 0 && ~ opts.evperdat
    plm.Mset{1} = ones(plm.N,1);
    opts.EE     = false;
    opts.ISE    = true;
elseif Nd > 0
    for m = 1:Nd
        Mtmp = palm_miscread(opts.d{m},[],[],opts.precision);
        plm.Mset{m} = Mtmp.data;
        if ~ isempty(plm.subjidx) && size(plm.Mset{m},1) ~= plm.N
            plm.Mset{m} = plm.Mset{m}(plm.subjidx,:);
        end
        if size(plm.Mset{m},1) ~= plm.N
            error([
                'The number of rows in the design matrix does\n' ...
                'not match the number of observations in the data.\n' ...
                '- Rows in the matrix: %d\n' ...
                '- Observations in the data: %d\n' ...
                'In file %s\n'], ...
                size(plm.Mset{m},1),plm.N,opts.d{m});
        end
        if any(isnan(plm.Mset{m}(:))) || any(isinf(plm.Mset{m}(:)))
            error([
                'The design matrix cannot contain NaN or Inf.\n' ...
                'In file %s\n'],opts.d{m});
        end
    end
end

% Include the EV per datum
if opts.evperdat
    for ev = 1:plm.nEVdat
        if ndims(plm.Mset{opts.evpos(ev,2)}) == 2 %#ok<ISMAT>
            plm.Mset{opts.evpos(ev,2)} = ...
                repmat(plm.Mset{opts.evpos(ev,2)},[1 1 plm.EVsiz(ev)]);
        end
        plm.Mset{opts.evpos(ev,2)}(:,opts.evpos(ev,1),:) = ...
            permute(plm.EVset{ev},[1 3 2]);
    end
    plm = rmfield(plm,{'EVset','EVsiz'});
end

% Some related sanity checks
if opts.evperdat
    for m = 1:plm.nM
        if opts.designperinput, loopY = m; else, loopY = 1:plm.nY; end
        for y = loopY
            if size(plm.Yset{y},2) ~= size(plm.Mset{m},3)
                error([
                    'The size of the data and the size of the EV per datum\n' ...
                    'don''t match.%s'],'')
            end
        end
    end
end

% Read and organise the contrasts for each design.
plm.Cset = cell(plm.nM,1);
plm.Dset = cell(plm.nM,1);
plm.nC   = zeros(plm.nM,1);
plm.nD   = zeros(plm.nM,1);
if Nt || Nf
    
    % Load FSL style t contrasts
    tcon = cell(Nt,1);
    for t = 1:Nt
        tmp = palm_miscread(opts.t{t},[],[],opts.precision);
        if any(strcmp(tmp.readwith,{'vestread','csvread','load'}))
            tcon{t} = tmp.data;
        else
            error('Invalid t contrast file: %s',opts.t{t});
        end
    end
    
    % Load FSL style F contrasts
    fcon = cell(Nt,1);
    for t = 1:Nt
        if ~ isempty(opts.f{t})
            tmp = palm_miscread(opts.f{t},[],[],opts.precision);
            if any(strcmp(tmp.readwith,{'vestread','csvread','load'}))
                fcon{t} = tmp.data;
            else
                error('Invalid F contrast file: %s',opts.f{t});
            end
        end
    end
    
    % For each valid design, assemble the contrasts.
    for m = 1:plm.nM
        if Nt == 1
            t = 1;
        else
            t = m;
        end
        c = 1;
        for j = 1:size(tcon{t},1)
            plm.Cset{m}{c} = tcon{t}(j,:)';
            c = c + 1;
        end
        for j = 1:size(fcon{t},1)
            if ~ isempty(fcon{t})
                plm.Cset{m}{c} = tcon{t}(logical(fcon{t}(j,:)),:)';
                c = c + 1;
            end
        end
        plm.nC(m) = numel(plm.Cset{m});
        for c = 1:plm.nC(m)
            plm.Dset{m}{c} = eye(plm.nY);
        end
        plm.nD(m) = numel(plm.Dset{m});
    end
    
elseif Ncon
    
    % Load MSET style contrasts (all contrast pairs)
    Ccon = cell(Ncon,1);
    Dcon = cell(Ncon,1);
    for con = 1:Ncon
        tmp = palm_miscread(opts.Ccon{con},[],[],opts.precision);
        if strcmpi(tmp.readwith,'mset')
            Ccon{con} = tmp.data;
        else
            error(['Files given to the option "-con" must be in .mset format.\n' ...
                'For .csv/.con/.fts files, use "-t" or "-f".%s'],'');
        end
        if isempty(opts.Dcon{con})
            for c = 1:numel(Ccon{con})
                Dcon{con}{c} = eye(plm.nY);
            end
        else
            tmp = palm_miscread(opts.Dcon{con},[],[],opts.precision);
            if strcmpi(tmp.readwith,'mset')
                Dcon{con} = tmp.data;
            else
                error(['Files given to the option "-con" must be in .mset format.\n' ...
                    'For .csv/.con/.fts files, use "-t" or "-f".%s'],'');
            end
        end
    end
    
    % Assign the contrast sets to the design matrices
    for m = 1:plm.nM
        if Ncon == 1
            con = 1;
        else
            con = m;
        end
        plm.Cset{m} = Ccon{con};
        plm.Dset{m} = Dcon{con};
        plm.nC(m) = numel(plm.Cset{m});
        plm.nD(m) = numel(plm.Dset{m});
        for c = 1:plm.nC(m)
            plm.Cset{m}{c} = plm.Cset{m}{c}';
        end
        for d = 1:plm.nD(m)
            plm.Dset{m}{d} = plm.Dset{m}{d}';
        end
    end
else
    % If no constrasts were at all specified:
    for m = 1:plm.nM
        if size(plm.Mset{m},2) == 1
            
            % If there is only 1 regressor, test its effect both
            % positive and negative.
            % The statistic will be t or v, depending on the number of VGs.
            plm.Cset{m}{1} = 1;
            plm.Cset{m}{2} = -1;
            plm.Dset{m}{1} = eye(plm.nY);
            plm.Dset{m}{2} = eye(plm.nY);
        else
            % Otherwise, run an F-test over all regressors in the design matrix.
            % The statistic will be F or G, depending on the number of VGs.
            plm.Cset{m}{1} = eye(size(plm.Mset{m},2));
            plm.Dset{m}{1} = eye(plm.nY);
        end
        plm.nC(m) = numel(plm.Cset{m});
        plm.nD(m) = numel(plm.Dset{m});
    end
end

% Ranks of the contrasts
plm.rC = cell(plm.nM,1);
plm.rD = plm.rC;
for m = 1:plm.nM
    plm.rC{m} = zeros(plm.nC(m),1);
    plm.rD{m} = plm.rC{m};
    for c = plm.nC(m):-1:1
        rC = rank(plm.Cset{m}{c});
        rD = rank(plm.Dset{m}{c});
        if rC == 0 || rD == 0
            plm.Cset{m}(c) = [];
            plm.Dset{m}(c) = [];
            plm.rC{m}(c)   = [];
            plm.rD{m}(c)   = [];
            warning('For design %d, contrast %d has rank 0 and will be removed.',m,c);
        else
            plm.rC{m}(c) = rC;
            plm.rD{m}(c) = rD;
        end
    end
    plm.nC(m) = numel(plm.Cset{m});
    plm.nD(m) = numel(plm.Dset{m});
end
plm.rC0 = plm.rC; % the rC can be changed for z-scores, but not rC0.

% If only the t or F tests are to be performed
if opts.tonly
    for m = 1:plm.nM
        for c = plm.nC(m):-1:1
            if plm.rC{m}(c) > 1
                plm.Cset{m}(c) = [];
                plm.Dset{m}(c) = [];
                plm.rC{m}(c)   = [];
                plm.rC0{m}(c)  = [];
            end
        end
        plm.nC(m) = numel(plm.Cset{m});
        plm.nD(m) = numel(plm.Dset{m});
    end
elseif opts.fonly
    for m = 1:plm.nM
        for c = plm.nC(m):-1:1
            if plm.rC{m}(c) == 1
                plm.Cset{m}(c) = [];
                plm.Dset{m}(c) = [];
                plm.rC{m}(c)   = [];
                plm.rC0{m}(c)  = [];
            end
        end
        plm.nC(m) = numel(plm.Cset{m});
        plm.nD(m) = numel(plm.Dset{m});
    end
end

% Some more sanity checks
for m = 1:plm.nM
    for c = 1:plm.nC(m)
        if any(isnan(plm.Cset{m}{c}(:))) || any(isinf(plm.Cset{m}{c}(:)))
            error('The constrasts cannot contain NaN or Inf.');
        end
        if size(plm.Cset{m}{c},1) ~= size(plm.Mset{m},2)
            error('The size of one or more contrasts don''t match the size of the respective design matrix.')
        end
    end
    for c = 1:plm.nD(m)
        if any(isnan(plm.Dset{m}{c}(:))) || any(isinf(plm.Dset{m}{c}(:)))
            error('The constrasts cannot contain NaN or Inf.');
        end
    end
end
if opts.MV
    if any(plm.nC ~= plm.nD)
        error('The number of C and D contrasts must be the same');
    end
    for m = 1:plm.nM
        for d = 1:plm.nD(m)
            if size(plm.Dset{m}{d},1) ~= plm.nY
                error('The size of one or more MV contrasts don''t match the number of modalities.');
            end
        end
    end
end
for m = 1:plm.nM
    for c = 1:plm.nC(m)
        if opts.concordant && plm.rC{m}(c) > 1
            error(['Cannot use the "-concordant" option with F-tests (inadmissible).\n'...
                'Use "-tonly" to run just the t-tests, or remove "-concordant".%s'],'');
        end
    end
end

% Make sure EVs of interest aren't represented also in the nuisance
% Note that some lines below the same is done for the cases in which 
% data and design are mean-centered
if ~ opts.demean && ~ opts.vgdemean && ~ opts.noranktest
    testrank(plm)
end

% Check if the contrasts have all the same rank for correction over
% contrasts. If not, convert to zstat.
if opts.corrcon && ~ opts.zstat
    rC1 = plm.rC{1}(1);
    rD1 = plm.rD{1}(1);
    for m = 1:plm.nM
        brflag = false;
        for c = 1:plm.nC(m)
            if rC1 ~= plm.rC{m}(c) || rD1 ~= plm.rD{m}(c)
                warning([...
                    'Not all contrasts have the same rank, and the option "-corrcon" was used.\n' ...
                    '         Adding the option "-zstat" automatically.%s'],'');
                opts.zstat = true;
                brflag = true;
                break;
            end
        end
        if brflag
            break;
        end
    end
end

% Partition the model according to the contrasts and design matrix.
% The partitioning needs to be done now, because of the need for
% synchronised permutations/sign-flips
if opts.corrcon || opts.npccon
    opts.syncperms = true;
end
if ~ opts.cmcx
    seqtmp = zeros(plm.N,sum(plm.nC));
    j = 1;
    plm.seq = cell(plm.nM,1);
    for m = 1:plm.nM
        plm.seq{m} = cell(plm.nC(m),1);
        for c = 1:plm.nC(m)
            Xtmp = palm_partition(plm.Mset{m}(:,:,1),plm.Cset{m}{c},opts.pmethodp);
            [~,~,plm.seq{m}{c}] = unique(Xtmp,'rows');
            seqtmp(:,j) = plm.seq{m}{c};
            j = j + 1;
        end
    end
    tmp = sum(diff(seqtmp,1,2).^2,2);
    if (opts.corrcon || opts.npccon || opts.syncperms) && any(tmp(:) ~= 0)
        warning([ ...
            'You chose to correct over contrasts, or run NPC\n'             ...
            '         between contrasts, but with the design(s) and,\n'     ...
            '         contrasts given it is not possible to run\n'          ...
            '         synchronised permutations without ignoring repeated\n'...
            '         elements in the design matrix (or matrices). To\n'    ...
            '         solve this, adding the option "-cmcx" automatically.%s\n'],'');
        opts.cmcx = true;
    end
end
if opts.cmcx % note can't use 'else' here because opts.cmcx is modified in the 'if' above
    plm.seq = cell(plm.nM,1);
    for m = 1:plm.nM
        plm.seq{m} = cell(plm.nC(m),1);
        for c = 1:plm.nC(m)
            plm.seq{m}{c} = (1:plm.N)';
        end
    end
end

% Make sure not too many components are asked if CCA or PLS are used
if opts.CCA || opts.PLS
    if opts.ccaorplsparm > plm.nY
        error(['Cannot ask more canonical correlations (CCA) or \n', ...
            'score vectors (PLS) (k=%d) than the number of modalities (#(-i)=%d).\n'],...
            opts.ccaorplsparm,plm.nY);
    end
    for m = 1:plm.nM
        for c = 1:plm.nC(m)
            if opts.ccaorplsparm > plm.rC{m}(c)
                error(['Cannot ask more canonical correlations (for CCA) or score \n', ...
                    'vectors (for PLS) than the rank of the contrast (k=%d > rank=%d).\n', ...
                    'Please check design %d, contrast %d.'],opts.ccaorplsparm,plm.rC{m}(c),m,c);
            end
        end
    end
end

% Read the exchangeability blocks. If none is specified, all observations
% are assumed to be in the same large block. Also treat the legacy format of
% a single column for the EBs.
if isempty(opts.eb)
    plm.EB = [];
    if opts.within || opts.whole
        error([ ...
            'Options -within and/or -whole require a file defining\n' ...
            '         the exchangeability blocks (option -eb).\n%s'],'');
    end
else
    plm.EB = palm_miscread(opts.eb);
    plm.EB = plm.EB.data;
    if ~ isempty(plm.subjidx)
        plm.EB = plm.EB(plm.subjidx,:);
    end
    if isvector(plm.EB)
        if opts.within && opts.whole % within + whole block shuffling
            plm.EB = [+ones(plm.N,1) +plm.EB(:) (1:plm.N)'];
        elseif opts.whole             % whole-block shuffling
            plm.EB = [+ones(plm.N,1) -plm.EB(:) (1:plm.N)'];
        else                          % within-block shuffling (this is the default, and not meant to be changed)
            plm.EB = [-ones(plm.N,1) +plm.EB(:) (1:plm.N)'];
        end
    elseif opts.within || opts.whole
        warning([ ...
            'Options -within and/or -whole ignored, as the file defining\n' ...
            '         the exchangeability blocks (option -eb) already \n' ...
            '         defines how the data should be shuffled.%s'],'');
    elseif ~ isempty(plm.subjidx)
        error('Cannot use "-subjidx" with multi-level blocks.');
    end
    plm.EB = palm_reindex(plm.EB,'fixleaves');
end

% Load/define the variance groups.
if opts.singlevg
    % If single VG, it's all ones
    plm.VG = ones(plm.N,1);
elseif strcmpi(opts.vg,'auto')
    if isempty(plm.EB)
        % If auto, but there are no exchangeability blocks, it's all ones too
        plm.VG = ones(plm.N,1);
    else
        % Generate an initial dependence tree, to be used to define variance groups.
        % The tree used for the permutations later require the design matrix, and
        % varies for each contrast -- all to be taken care of later.
        Ptree  = palm_tree(plm.EB,(1:plm.N)');
        plm.VG = palm_ptree2vg(Ptree);
    end
else
    % The automatic variance groups can be overriden if the user specified
    % a file with the custom definitions.
    plm.VG = palm_miscread(opts.vg);
    plm.VG = plm.VG.data;
end
if ~ isempty(plm.subjidx) && size(plm.VG,1) ~= plm.N
    plm.VG = plm.VG(plm.subjidx,:);
end
[tmp,~,plm.VG] = unique(plm.VG);
plm.nVG = numel(tmp);
if plm.nVG == 1, opts.singlevg = true; end

% MV can't yet be used if nVG>1, although NPC remains an option
if opts.MV && plm.nVG > 1
    error('There are more than one variance group. MV cannot be used (but NPC can).');
end

% There should be no more missing indicators than modalities, or designs.
% These need to be either 1 or the same as the corresponding numbers of
% modalities/designs.
if Nimiss > Ni
    error([...
        'There are more missing indicators supplied with "-imiss" (%d) than\n'...
        'modalities supplied with "-i" (%d)'],Nimiss,Ni);
elseif Nimiss > 1 && Nimiss ~= Ni
    error([...
        'The number of missing indicators supplied with "-imiss" (%d) is larger,\n'...
        'than 1, but still not the same as the number of modalities supplied with\n'...
        'the option "-i" (%d).'],Nimiss,Ni);
end
if Ndmiss > Nd
    error([...
        'There are more missing indicators supplied with "-dmiss" (%d) than\n'...
        'designs supplied with "-d" (%d)'],Nimiss,Ni);
elseif Ndmiss > 1 && Ndmiss ~= Nd
    error([...
        'The number of missing indicators supplied with "-dmiss" (%d) is larger,\n'...
        'than 1, but still not the same as the number of modalities supplied with\n'...
        'the option "-d" (%d).'],Ndmiss,Nd);
end

% Load the missing indicators for the data ("imiss"):
for i = 1:Nimiss
    if strcmpi(opts.imiss{i},'none')
        if isempty(plm.subjidx)
            tmp = zeros(plm.N,1);
        else
            tmp = zeros(size(plm.subjidx,1),1);
        end
    else
        tmp = palm_miscread(opts.imiss{i});
        tmpfname = tmp.filename;
        tmp = tmp.data;
        if ~ isempty(plm.subjidx) && size(tmp,1) ~= plm.N
            tmp = tmp(plm.subjidx,:);
        end
        checkmiss(tmp,tmpfname,plm.N);
    end
    plm.Ymiss{i} = tmp;
end
if Nimiss == 1
    for i = 2:Ni
        plm.Ymiss{i} = plm.Ymiss{1};
    end
end

% Load the missing indicators for the design ("dmiss"):
for d = 1:Ndmiss
    if strcmpi(opts.dmiss{d},'none')
        if isempty(plm.subjidx)
            tmp = zeros(plm.N,1);
        else
            tmp = zeros(size(plm.subjidx,1),1);
        end
    else
        tmp = palm_miscread(opts.dmiss{d});
        tmpfname = tmp.filename;
        tmp = tmp.data;
        if ~ isempty(plm.subjidx) && size(tmp,1) ~= plm.N
            tmp = tmp(plm.subjidx,:);
        end
        checkmiss(tmp,tmpfname,plm.N);
    end
    plm.Mmiss{d} = tmp;
end
if Ndmiss == 1
    for d = 2:Nd
        plm.Mmiss{d} = plm.Mmiss{1};
    end
end
for d = 1:Ndmiss
    if any(size(plm.Mmiss{d}) ~= size(plm.Mset{d}))
        if strcmpi(opts.dmiss{d},'none')
            plm.Mmiss{d} = repmat(plm.Mmiss{d},[1 size(plm.Mset{d},2)]);
        else
            error([ ...
                'The missing data indicator ("-dmiss") must have\n', ...
                'the same size as the respective design.%s'],'');
        end
    end
end

% If only data or design missing indicators are missing, fill the other
% with all-false indicators.
if     Nimiss && ~ Ndmiss
    for m = 1:plm.nM
        plm.Mmiss{m} = false(size(plm.Mset{m}));
    end
elseif Ndmiss && ~ Nimiss
    for y = 1:plm.nY
        plm.Ymiss{y} = false(size(plm.Yset{y}));
    end
end

% Remove the variance groups with tiny sample sizes?
if plm.nVG > 1 && ~ opts.removevgbysize && (opts.vgdemean || opts.ev4vg) && ...
        any(sum(bsxfun(@eq,plm.VG,unique(plm.VG)'),1) == 1)
    warning([...
        'The options "-vgdemean" and "-ev4vg" require that observations\n' ...
        '         in variance groups of size 1 are removed.\n' ...
        '         Enabling the option "-removevgbysize 1"%s.'],'');
    opts.removevgbysize = 1;
end
if ~ opts.removevgbysize
    tmpwarned = false;
    for u = 1:plm.nVG
        if sum((plm.VG == u),1) == 1
            if ~ tmpwarned
                warning([...
                    'There are variance groups with just one observation.\n' ...
                    '         Consider using the option "-removevgbysize 1" to improve the\n' ...
                    '         variance estimates (at the cost of reducing sample size).%s'],'');
                tmpwarned = true;
            end
        end
    end
end
if opts.removevgbysize > 0
    
    % Indices of the observations to keep
    uVG = unique(plm.VG)';
    idxvg = sum(bsxfun(@eq,plm.VG,uVG),1) <= opts.removevgbysize;
    idx   = any(bsxfun(@eq,plm.VG,uVG(~idxvg)),2);
    
    % Modify all data as needed
    for y = 1:plm.nY
        plm.Yset{y} = plm.Yset{y}(idx,:);
    end
    if ~ isempty(plm.EB)
        plm.EB = plm.EB(idx,:);
    end
    for m = 1:plm.nM
        plm.Mset{m} = plm.Mset{m}(idx,:);
    end
    for m = 1:plm.nM
        for c = 1:plm.nC(m)
            plm.seq{m}{c} = plm.seq{m}{c}(idx,:);
        end
    end
    plm.N = sum(idx);
    [tmp,~,plm.VG] = unique(plm.VG(idx));
    plm.nVG = numel(tmp);
end

% Add one regressor for each variance group if requested
if opts.ev4vg
    for m = 1:plm.nM
        Mvg = zeros(plm.N,plm.nVG);
        V = unique(plm.VG);
        for v = 1:plm.nVG
            Mvg(plm.VG == V(v),v) = 1;
        end
        rM   = round(sum(diag(plm.Mset{m}*pinv(plm.Mset{m}))));
        Mnew = horzcat(plm.Mset{m},Mvg);
        if round(sum(diag(Mnew*pinv(Mnew)))) == (rM + plm.nVG)
            plm.Mset{m} = Mnew;
            nadded      = plm.nVG;
        else
            Mnew = Mnew(:,1:end-1);
            if round(sum(diag(Mnew*pinv(Mnew)))) == (rM + plm.nVG - 1)
                plm.Mset{m} = Mnew;
                nadded      = plm.nVG - 1;
            else
                error([ ...
                    'It was not possible to add one regressor for each variance group\n' ...
                    'perhaps because they already exist in the design. Check your design\n' ...
                    'matrix and maybe consider including these regressors manually.%s'],'');
            end
        end
        for c = 1:plm.nC(m)
            plm.Cset{m}{c} = vertcat(plm.Cset{m}{c},...
                zeros(nadded,size(plm.Cset{m}{c},2)));
        end
    end
end

% Make sure that mean-centering won't damage a contrast that tests the intercept.
if opts.demean || opts.vgdemean
    for m = 1:plm.nM
        for c = 1:plm.nC(m)
            % The design may be 3D here if -evperdat was used, hence can't
            % do matrix multiplication of design by contrast
            for t = 1:size(plm.Cset{m}{c},2)
                Xtmp = sum(bsxfun(@times,plm.Mset{m},plm.Cset{m}{c}(:,t)'),2);
                siz  = size(Xtmp);
                isintercp = all(bsxfun(@eq,reshape(Xtmp(1,:),[1 siz(2:end)]),Xtmp),1);
                if isintercp
                    error([ ...
                        'Contrast %d for design %d (and perhaps others) tests the intercept.\n' ...
                        'This means that the options "-demean" and "-vgdemean" cannot be used.\n' ...
                        'If "-demean" was added to calculate Pearson''s "r" or the "R^2"\n' ...
                        'note that these statistics cannot be computed for constant variables.%s'],c,m,'');
                end
            end
        end
    end
end

% Mean center data and design (-demean)
if opts.demean
    for m = 1:plm.nM
        plm.Mset{m} = bsxfun(@minus,plm.Mset{m},mean(plm.Mset{m},1));
    end
    for y = 1:plm.nY
        plm.Yset{y} = bsxfun(@minus,plm.Yset{y},mean(plm.Yset{y},1));
    end
end

% Mean center data and design, within VG
if opts.vgdemean
    
    % For each VG
    V = unique(plm.VG);
    for v = 1:plm.nVG
        vidx = plm.VG == V(v);
        
        % Demean design within VG
        for m = 1:plm.nM
            plm.Mset{m}(vidx,:) = bsxfun(@minus,...
                plm.Mset{m}(vidx,:),mean(plm.Mset{m}(vidx,:),1));
        end
        
        % Demean data within VG
        for y = 1:plm.nY
            plm.Yset{y}(vidx,:) = bsxfun(@minus,...
                plm.Yset{y}(vidx,:),mean(plm.Yset{y}(vidx,:),1));
        end
    end
end

% Make sure EVs of interest aren't represented also in the nuisance
% Note that some lines above the same is done for the cases in which there
% is no mean-centering
if (opts.demean || opts.vgdemean) && ~ opts.noranktest
    testrank(plm)
end

% Number of tests to be selected for the low rank approximation
if opts.accel.lowrank
    if plm.nVG > 1
        error('The option "-accel lowrank" cannot be used with more than one variance group.');
    end
    if opts.nP0 == 0
        error('With lowrank approximation you must indicate a larger-than-zero number of permutations.');
    end
    if opts.nP0 < plm.N*(plm.N+1)/2
        error([ ...
            'Too few permutations selected to use with lowrank approximation.\n' ...
            'Use at least N*(N+1)/2 = %d to note a speed difference and have reasonably accurate results.\n'...
            'Otherwise, don''t bother using lowrank approximation.\n'],plm.N*(plm.N+1)/2);
    end
    if opts.spatial.do
        warning([ ...
            'There isn''t much benefit in using lowrank approximation with spatial statistics\n' ...
            '         like TFCE and cluster extent and/or mass. These cannot be accelerated with this\n' ...
            '         method, and the overall gain will be minimal. Consider other approximation methods,\n' ...
            '         or run the full permutation test, or just drop spatial statistics.%s'],'');
    end
    plm.nsel = zeros(plm.nY,1);
    if opts.accel.lowrank_val <= 1
        for y = 1:plm.nY
            plm.nsel(y) = ceil(opts.accel.lowrank_val*plm.Ysiz(y));
        end
    elseif opts.accel.lowrank_val > 1
        plm.nsel(1:end) = ceil(opts.accel.lowrank_val);
    elseif isnan(opts.accel.lowrank_val)
        plm.nsel(1:end) = plm.N*(plm.N+1)/2;
    end
end

% ==============================================================
function testrank(plm)
% Test the complicated case in which design is rank deficient, with
% redundant dimensions being equally represented in EVs of interest and in
% nuisance.
badcon = cell(plm.nM,1);
for m = 1:plm.nM
    badcon{m} = zeros(plm.nC(m),1);
    for c = 1:plm.nC(m)
        [Xgutt,Zgutt] = palm_partition(plm.Mset{m}(:,:,1),plm.Cset{m}{c},'Guttman');
        if size(Zgutt,2) > 0
            cc = palm_cca(Xgutt,Zgutt,plm.N);
            if cc(1) == 1
                badcon{m}(c) = sum(cc == 1);
            end
        end
    end
end
idxbad = vertcat(badcon{:});
if any(idxbad)
    badlist = zeros(2,sum(idxbad > 0));
    j = 1;
    for m = 1:plm.nM
        for c = 1:plm.nC(m)
            if badcon{m}(c)
                badlist(1,j) = m;
                badlist(2,j) = c;
                badlist(3,j) = badcon{m}(c);
                j = j + 1;
            end
        end
    end
    badmsg = sprintf('- Design %d, Contrast %d, at least %d regressors\n',badlist);
    error([...
        'The following contrasts try to test regressor(s) also fully represented by\n'...
        'by nuisance variable(s), but such tests are not possible (rank deficiency):\n%s'],badmsg); %#ok<SPERR>
end

% ==============================================================
function checkmiss(A,Afname,N)
% Check if the missing data indicators are sane.
for a = 1:size(A,2)
    U = unique(A(:,a));
    if size(A,1) ~= N
        error([ ...
            'The missing data indicators ("-imiss" and "-dmiss") must have\n', ...
            'the same number of observations as the data and design.'],'');
    elseif ...
            (numel(U) >  2) || ...
            (numel(U) == 2 && ~ any(U == 0)) || ...
            (numel(U) == 2 && ~ any(U(U~=0) == [-1 1 2])) || ...
            (numel(U) == 1 && U ~= 0)
        error([ ...
            'The missing data indicators ("-imiss" and "-dmiss") must have\n', ...
            'no more than two unique values per column, one being 0, the\n', ...
            'the other being either -1, 1, or 2.\n', ...
            'Consult the documentation for details.\n'...
            '- Filename: %s\n',...
            '- Column: %d (possibly also others)'],Afname,a);
    end
end