function dccn_hemi_cross2( varargin)

% WHAT:     Cross-decoding of object shape/size (training vs test runs)
% INPUTS:   [shape] or [size], [f]orward or [b]ackward train<->test runs
% FILE:     Cross_[shape/size]_[trainmodel]_to_[testmodel]_[S]steps_([N]_testspecs)
% ROIs:     step:step:median(roi), and all p<.05 voxels 
% Stats:    resampling vs shuffling-based TFCE, & cluster-based permutation


%% preparation

% check if running file as a function (default) or script (for debugging)
try 
    xx = ~nargin;
catch
    cd('/project/2418046.01/TRANSFER_Aug_2020');
    addpath( fileparts( mfilename('fullpath')));
    clearvars; clc;
end

% pilot or not
pilot = 0;


%% Display settings

% Select ROIs to display
plotit      =   1;                  % set to 1 to plot results   
show.rois   =   {'LOC','EVC'};      % 'LOC','EVC','PPA','TOS','RSC','FFA'
show.hemi   =   show.rois;          % plot data per hemifield with density plots
show.fixedy =   1;                  % set to 1 to use user-defined y-axes (otherwise defined automatically to fit data)
show.saveit =   1;                  % set to 1 to save figures as svg files

% select which conditions to display
show.plot   =   [1 2 3 4];          % plot which conditions             (1 match, 2 mismatch, 3 difference, 4 average)
show.bars   =   [3 4];              % plot error bars for ...           
show.indiv  =   [];                 % plot individual data for ...      
show.signi  =   [show.plot];        % plot significance for ...          

% select which metric to report
show.metric =   'rdb';              % 'acc' (standard classification accuracy), 'rdb' (rank-order averaged distance to boundary)
show.null   =   'labs';             % create null distributions by ('perm') permuting data with random sign, or ('labs') by shuffling test labels (applies to tfce and cluster analyses, not to condition-wise statistics)
show.multi  =   'tfce';             % 'tfce', 'clus', or 'none' (for TFCE-based, cluster-based, or no correction for multiple comparisons across voxel sizes
show.stats  =   'ttype';            % 'prob' or 'ttype', for all (or none) statistical tests (prob is probability based on permuted label shufflings, ttype is test ensuing from 'ttype' variable


%% Classification settings

% define ROIs
roi_regs    =   {'LOL',  'LOR',  'EVL',  'EVR'};%,  'PPL',  'PPR',  'TOL',  'TOR',  'RSL',  'RSR',  'FFL',  'FFR'  };
roitypes    =   {'loca', 'loca', 'trai', 'trai'};%, 'loca', 'loca', 'loca', 'loca', 'loca', 'loca', 'loca', 'loca' };
roinames    =   {'LOC',  'LOC',  'EVC',  'EVC'};%,  'PPA',  'PPA',  'TOS',  'TOS',  'RSC',  'RSC',  'FFA',  'FFA'  };  
nsteps      =   20;                             % number of ROI sizes (within a brain region)
rvox        =   1;                              % minimum step size (in number of voxels)

% metrics
def.metrics     =   {'rdb'};
def.null        =   {'labs'};

% subjects
subjects    =   1:24;

% models
model_train =   'training_s3';           
model_test  =   'testing_4c_s3';        

% test name
mytest      =   'shape';         % cross classification of 'shape' or 'size'

% direction of classification
classdir    =   'f';            % ['f']orward or ['b']ackward classification (train -> test, or test -> train)

% set roots
paths.proj      =   pwd;
paths.subj      =   fullfile( paths.proj, 'SUBJECT level');
paths.local     =   fileparts( mfilename( 'fullpath'));

% classifier type and feature selection
classifier          =   'libsvm';       
feature_selection   =   'none';         

% group level analyses
ttype               =   'bootstrap';    % 'bootstrap', 'ttest', 'rank' 
nsamples            =   2000;           % used for bootstrap tests and cluster permutation test
sided               =   1;              % one or two sided test (note: cluster correction ensures alpha level of 0.05 for 1 sided testing)
cutoff              =   0.05;           % alpha level (1st level significance threshold)

% behavior
paths.behav         =   fullfile( paths.subj, '..', 'Behavioral', 'melbox1_behav.mat');
show.cor            =   1;              % corelate with behavior yes (1), no (0), or selection (2)
show.cdata          =   'all';          % neural data is full ROI ('all'), mean of cluster-based significant ROI-sizes ('clus'), or mean of tfce-based significant ROI-sizes ('tfce')
ctype               =   'Kendall';      % 'Kendall', 'Pearson', or 'Spearman'

% Overwrite analyses ?
def.runit           =   0;              % (0) show results (if they exist), or (1) run new analysis (forcefully)


%% Set conditions and labels

% overwrite defaults with user requests
if exist( 'varargin', 'var')
    for nvar = 1:2:length( varargin)
        if exist( varargin{nvar}, 'var')
            switch class( varargin{nvar+1})
                case 'char'
                    eval( sprintf( "%s = '%s';", varargin{nvar}, varargin{nvar+1}));
                case 'double'
                    eval( sprintf( "%s = %g;", varargin{nvar}, varargin{nvar+1}));
            end
        else
            error( '[%s] is not a valid variable name. Program aborted.', varargin{nvar}); 
        end
    end
end
def.plotit  =   plotit;
def.pilot   =   pilot;

% determine labels
switch lower(mytest)
    
    % shape cross-classification
    case 'shape'
    
        labels_train = {{'mel_near', 'box_near'}, ...
                        {'mel_far',  'box_far' }};
 
        labels_test  = {{{'mel_near', 'box_near'}, {'mel_far',  'box_far' }}, ...
                        {{'mel_far',  'box_far' }, {'mel_near', 'box_near'}}};
    
    % size cross-classification            
    case 'size'
        
        labels_train = {{'mel_near', 'mel_far'}, ...
                        {'box_near', 'box_far'}};

        labels_test  = {{{'mel_near', 'mel_far'}, {'box_near', 'box_far'}}, ...
                        {{'box_near', 'box_far'}, {'mel_near', 'mel_far'}}};
        
end

% conditions
conds       =   {'matching', 'mismatch', ...
                 'matching', 'mismatch'};    

% reverse classification direction: [test -> train] OR [train -> test]                
if strcmpi( classdir, 'b')
    [model_train,  model_test ] =   deal( model_test,  model_train );
end


%% Retrieve and/or define ROIs

% get suffix to retrieve correct file types
mysm = strfind( model_test, '_s');
if isempty( mysm)
    suffix = '';
else
    suffix = model_test(mysm(end):mysm(end)+2);
end

% add some relevant paths (functions, qsub, cosmomvpa, roi directory, etc.)
addpath(        fullfile( paths.proj, 'GROUP level', 'SPM', 'functions'));
addpath(        '/home/common/matlab/fieldtrip/qsub');
paths.cosmo =   fullfile( paths.proj, 'Toolboxes', 'CoSMoMVPA-master', 'mvpa');
roi_dir     =   fullfile( 'roi', 'mni2');        
voxfile     =   fullfile( paths.local, 'roi_vox.mat');  % file with max_voxels (total number of significant voxels per participant per ROI (at p<.05 uncorrected contrast threshold)

% def.pilot settings
if def.pilot
   roi_regs =   roi_regs(1:4);
   roitypes =   roitypes(1:4);
   roinames =   roinames(1:4);
   subjects =   subjects(12:18);
   nsamples =   1000;
   nsteps   =   4;
end

% retrieve roi sizes
while 1
    if exist( voxfile, 'file')
        
        v       =   load( voxfile);
        
        if all( contains( roi_regs, v.roi_regs))
            if def.plotit
                roi_sizes( roi_regs, roitypes, subjects, suffix, paths.subj, paths.local, roi_dir, 0, 1);
                
                % save figure
                if show.saveit
                    figpath     =   fullfile( paths.local, 'images');
                    if ~exist( figpath, 'dir')
                        mkdir( figpath);
                    end
                    
                    set( gcf, 'Renderer', 'painters')
                    figname = sprintf( 'Voxelcount_%g_ROIs', length( roi_regs)/2);
                    print( gcf, fullfile( figpath, figname), '-depsc');
                end
            end
            break
        else
            roi_sizes( roi_regs, roitypes, subjects, suffix, paths.subj, paths.local, roi_dir, 1, 0);
        end
        
    else
        roi_sizes( roi_regs, roitypes, subjects, suffix, paths.subj, paths.local, roi_dir, 1, 0);
    end
end

% roi preparation
nregs       =   length( roi_regs);
rois        =   [];
voxmax      =   NaN( 1, nregs);
voxstep     =   voxmax;

% create cell array with ROI names
for nroi = 1 : nregs
    
    % find correct ROI in file with roi sizes
    irois       =   contains( v.roi_regs, roi_regs( contains( roinames, roinames{nroi})));
    
    % retrieve maximum ROI size
    thismax         =   round( mean( median( v.nvox(:, irois))));   % find median ROI size
    
    % get step size, and ensure endpoint difference < .05 
    if (thismax/nsteps) < rvox*20
        voxstep(nroi)   =   round( (thismax/nsteps));
    else
        voxstep(nroi)   =   round( (thismax/nsteps)/rvox) * rvox;       % compute stepsize
    end
    
    % compute new max, and store all ROI sizes
    newmax          =   voxstep(nroi)*nsteps;
    myvoxs          =   voxstep(nroi) : voxstep(nroi) : newmax;     % define voxel counts
    voxmax(nroi)    =   myvoxs(end);
    
    % stack ROI names
    rois            =   [ rois, ...
                          strcat( compose( '%s_%s_%g', roi_regs{nroi}, roitypes{nroi}, myvoxs), suffix), ...
                          sprintf( '%s_%s_all%s', roi_regs{nroi}, roitypes{nroi}, suffix) ];        
end


%% data storage

% handy variables
nSubjects   =   length( subjects);
nPairs      =   length( conds);
nrois       =   length( rois);
n_trains    =   length( labels_train);
myconds     =   unique( conds, 'stable');
wpairs      =   cellfun( @(x) strcmp( x, conds), myconds, 'UniformOutput', 0);

samesvm     =   NaN( 1, nPairs);
it = 0;
for n_train = 1:n_trains
    for n_test = 1:length( labels_test{n_train})
        it = it+1;
        samesvm(1,it) = n_train;
    end
end

% determine or create output directory (group level)
paths.res  =   fullfile( paths.local, 'classification_results');
if ~exist( paths.res,'dir')
    mkdir( paths.res)
end
if ~exist( fullfile( paths.res, 'temp'),'dir')
    mkdir( fullfile( paths.res, 'temp'))
end

% determine filename
paths.f         =   fullfile( paths.res, sprintf('C_cross_hemi_%s_%s_to_%s_%gsteps_(N=%g_%s%g_n%g).mat', mytest, model_train, model_test, nsteps, nSubjects, ttype(1), sided, nsamples));
paths.sim       =   fullfile( paths.res, sprintf('C_cross_hemi_%s_%s_to_%s_%gsteps_(N=%g*.mat',          mytest, model_train, model_test, nsteps, nSubjects));
paths.trial     =   strrep(   paths.sim, '*', '_trial-data)');

% create randomized name for temporary files
s = 'a':'z';
paths.temp = s(randi(numel(s), 1, 10));


%% Determine whether to run or load classification

% load data from file if not forcing to run new classification analysis
if ~def.runit
    if exist( paths.f, 'file')
        load( paths.f); %#ok<*LOAD>
    else
        
        % list files using similar classification
        f2      =   dir( paths.sim);

        if isempty( f2)
            def.runit = 1; % run classification nonetheless if no matching file is found
        else
            paths.f2  =   fullfile( paths.res, f2(1).name); 
            load( paths.f2);
        end
    end
end        


%% Run classification

% Run classification
if def.runit
    
    % Preparation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % add paths for TDT and SPM
    tdt_ver = dir( fullfile( paths.proj, 'Toolboxes', 'decoding_toolbox_*')); % find TDT versions
    addpath( fullfile( tdt_ver.folder, tdt_ver(end).name)); 
    addpath( '/home/common/matlab/spm12');
    
    % initialize SPM
    spm defaults fmri
    spm_jobman initcfg
    spm_get_defaults('cmdline',true)
    
    % clean up
    cleanit;
        
    % Check ROIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Create ROIs (if not yet done)
    firstit = 1; oldsl = inf; 
    while 1
        
        % check whether all requested ROI's exist in all subjects
        check_it        =   NaN(nrois,nSubjects);
        for nsub = 1:nSubjects
            
            temp_list       =   dir( fullfile( paths.subj,  sprintf( '%02g', subjects(nsub)), roi_dir));
            name_list       =   {temp_list(3:end).name};
            
            % check if any ROI is missing for current subject
            check_it(:,nsub) = cellfun( @(x) any( contains( name_list, sprintf( '%s.nii', x))), rois)';
        end
        runsubs = any( check_it==0);
        
        % if any subjects are missing ROIs
        if any( runsubs)
            
            % make ROIs for these specific subjects
            if firstit
                
                % set memory/cpu requirement for cluster job
                mem_roi   =   1;     % cluster job memory requirement per participant, in gb
                time_roi  =   0.5;   % cluster job walltime requirement per participant, in h
                
                % send message
                fprintf( '\nNot all requested ROIs were found. \nNow building ROIs before starting classification.\n')
                
                % add path and run function to create ROIs
                addpath( fullfile( paths.proj, 'GROUP level', 'roi'));
                
                %ROI_pipeline_subfunc( subjects(runsubs), roi_regs, voxstep, voxmax, '2')
                for ss = subjects(runsubs)
                    qsubfeval( @ROI_pipeline_subfunc, ss, roi_regs, voxstep, voxmax, '2', 'memreq', mem_roi*1024^3, 'timreq', time_roi*3600);
                end
                firstit = 0;
            
            else
                % check every N second if all ROIs are found
                if sum(runsubs) < oldsl 
                    fprintf( '\nStill working on ROIs for %g subjects...\n', sum(runsubs));
                    oldsl = sum(runsubs);
                end
                pause(2);
                
            end
        else
            % wait until all ROIs are found
            fprintf( '\nAll ROIs created.\n\n')
            break
        end
    end
    
    
    % Create classification struct %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % load defaults
    cfg                     =   decoding_defaults;
    cfg.testmode            =   0;
    cfg.results.overwrite   =   1;
    
    % classification settings
    cfg.software            =   'SPM12';
    cfg.decoding.software   =   classifier;
    cfg.decoding.method     =   'classification';
    
    if strcmp( classifier, 'neCTon')
        cfg.decoding.train.classification.model_parameters = 0;
    elseif strcmp( classifier, 'libsvm')
        cfg.decoding.train.classification.model_parameters = '-s 0 -t 0 -c 1 -b 0 -q'; % linear classification
    end
    
    % scaling settings
    cfg.scale.method         =  'none'; % 'z', 'min0max1', 'none'
    
    if strcmp( cfg.scale.method, 'none')
        cfg.scale.estimation =  'none';
    else
        cfg.scale.estimation =  'all';
    end
    
    % feature selection
    if contains( feature_selection, 'filter')
        cfg.feature_selection.method = 'filter';
        cfg.feature_selection.filter = feature_selection( strfind( feature_selection,'_')+1:end); % values: F, F0, F1, U, W, external
    elseif contains( feature_selection, 'rankfile')
        cfg.feature_selection.method = 'filter';
        cfg.feature_selection.filter = 'external';
    else
        cfg.feature_selection.filter = nan;
        cfg.feature_selection.method = feature_selection;
    end
    cfg.feature_selection.n_vox             =   0.005:0.005:1;                  % range or 'automatic'
    cfg.feature_selection.nested_n_vox      =   0.005:0.005:1;
    cfg.feature_selection.decoding.method   =   'classification';
    
    % output
     cfg.results.output     =   {'accuracy_minus_chance', 'decision_values', 'predicted_labels'};   % 'SVM_weights'
    
    % verbosity options
    cfg.plot_selected_voxels                    =   0;
    cfg.plot_design                             =   0;
    cfg.verbose                                 =   0; % NB: if set to zero, only local fprintf commands are printed
    
    % set ROI settings and file name in cfg struct
    cfg.analysis    =   'ROI';
    cfg.savename    =   'temp';

    
    % Start iteration over classification jobs %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % set memory/cpu requirement for cluster job
    mem_svm   =   2;    % cluster job memory requirement per participant, in gb
    time_svm  =   4;    % cluster job walltime requirement per participant, in h
    
    % iterate over subjects
    for s = 1 : nSubjects
        
        % set and report current subject
        subj            =   sprintf( '%02g', subjects(s));
        subject_dir     =   fullfile( paths.subj,  subj);
        d.tempath       =   fullfile( paths.res, 'temp', sprintf( '%s-%s.mat', paths.temp, subj));
        
        fprintf( '\nSubject %g / %g [s=%02s]\n', s, nSubjects, subj)
        
        % locate model and ROI
        d.beta_dir_train    =   fullfile( subject_dir, 'univariate', model_train);
        d.beta_dir_test     =   fullfile( subject_dir, 'univariate', model_test);
        d.sub_roi_dir       =   fullfile( subject_dir, roi_dir);
                    
        % determine regressor names (run-based or trial-based models)
        if strcmp( 'st_', model_train(1:3))
            load( fullfile( d.beta_dir_train, 'betas', sprintf( '%s_beta_info.mat', subj)));
            
            d.regressor_names_train     =   trialinfo(2:end, [6, 2, 8, 4])'; %#ok<*FNCOLND>
            d.beta_dir_train            =   cellfun( @(x) fullfile( d.beta_dir_train, 'betas', x), d.regressor_names_train(3,:), 'UniformOutput', 0)'; 
            
        else
            d.regressor_names_train     =   design_from_spm( d.beta_dir_train);
        end
        
        if strcmp( 'st_', model_test(1:3))
            load( fullfile( d.beta_dir_test, 'betas', sprintf( '%s_beta_info.mat', subj)));
            
            d.regressor_names_test      =   trialinfo(2:end, [6, 2, 8, 4])';
            d.beta_dir_test             =   cellfun( @(x) fullfile( d.beta_dir_test, 'betas', x), d.regressor_names_test(3,:), 'UniformOutput', 0)'; 
            
        else
            d.regressor_names_test      =   design_from_spm( d.beta_dir_test);
        end
        
        % determine (or create) output directory
        cfg.results.dir =   fullfile( subject_dir, 'multivariate', sprintf( 'cross_%s', model_test));
        if ~exist( cfg.results.dir, 'dir')
            mkdir( cfg.results.dir);
        end
        
        % Transfer jobs to cluster %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        % pack used variables into struct
        d.rois                      =   rois;
        d.labels_train              =   labels_train;
        d.labels_test               =   labels_test;
        d.p_root                    =   paths.proj;
        d.subj                      =   subj;
                
        % choose classification method
        if strcmp( model_train, model_test)
            d.cl_method        =   'xclass_cv';
        else
            d.cl_method        =   'xclass';
        end
        
        % run decoding (for current participant)
        qsubfeval( @decode_rois, cfg, d, 'memreq', mem_svm*1024^3, 'timreq', time_svm*3600);
      
    end % end of subject loop
    
    
    % wait for all files to arrive from cluster
    waitforsubs =   subjects;
    waitfordata =   5;
    fprintf( '\nRetrieving classification results from cluster ...\n\t> Data sets remaining: %g... ', length(subjects))
    
    % repeatedly check for new files
    while ~isempty( waitforsubs)
        
        % iterate over (missing) subjects
        for ss = waitforsubs
            
            % check wether temp file was created for current subject
            tempath = fullfile( paths.res, 'temp', sprintf( '%s-%02g.mat', paths.temp, ss));
            
            if exist( tempath, 'file')
                
                % load and stack data
                r = load( tempath, 'data');
                for n = 1:nPairs
                    data(n).dist(subjects==ss, :, :)    =   r.data(n).dist; %#ok<*SAGROW>
                    data(n).labels                      =   r.data(n).labels;
                    data(n).chunks                      =   r.data(n).chunks;   
                end
                                
                % delete temporary data
                waitforsubs( waitforsubs==ss) = [];
                delete( tempath);
                
                % reduce waiting time
                waitfordata = 0.5;
                fprintf( '%g... ', length( waitforsubs));
            end    
        end  
        
        % wait N seconds before checking for new data again
        pause(waitfordata);
    end
    
    % save classification results
    fprintf( '\nClasification results collected, and saved to results folder.');
    save( paths.f, 'data');
    
    % clean up big matrices to free memory
    clearvars cfg cfg2;
end

% run stats if classification data was loaded with different specs
if isfield( paths, 'f2')
    def.runit = 1;
end


%% Create classification metrics, permute and organize test data

% set null models
if isempty( def.null)
    def.null    =   {'perm', 'labs'};
end

% set metrics
if isempty( def.metrics)
    def.metrics =   {'acc', 'pfc', 'mfc', 'rdb', 'zdb'};
end

if def.runit || ~exist( 'RS', 'var') 
    
    % preallocate data matrices/structs
    s   =   struct();
    RS  =   struct();
    
    % check whether event-related design or not
    if contains( model_test, 'testing') && contains( model_test, 'st_')
        save_dpt = 1;
    else
        save_dpt = 0;
    end

    % if so: retrieve condition variables per trial 
    if save_dpt
        
        % pre-allocate data matrices / structs
        con_pt  =   struct();
        t       =   struct();
        CT     =   struct(); 
        
        % iterate over subjects
        for ip = 1:nSubjects
            
            % subject name
            subj = sprintf( '%02g',subjects(ip));
            
            % load data of both sessions
            b1_file = fullfile( paths.proj, 'SUBJECT level', subj, 'behavior', 'testing', sprintf( '%s_testing_08.mat', subj));
            b2_file = fullfile( paths.proj, 'SUBJECT level', subj, 'behavior', 'testing', sprintf( '%s_testing_16.mat', subj));
            
            b1= load( b1_file, 'behav');
            b2= load( b2_file, 'behav');
            
            % concatenate and add trial number (to be sure)
            ppconds                         =   cat( 2, b1.behav.cond, b2.behav.cond);
            ppconds(end+1, :)               =   1 : size(ppconds,2);
            ppconds(:, ppconds(1,:)==1)     =   [];
            
            % restructure based on SPM/TDT file order (regressor names file)
            ppconds = [ ppconds(:, b1.behav.cond( 4, :)==1 & b1.behav.cond( 3, :)==1), ...
                        ppconds(:, b1.behav.cond( 4, :)==1 & b1.behav.cond( 3, :)==2), ...
                        ppconds(:, b1.behav.cond( 4, :)==2 & b1.behav.cond( 3, :)==1), ...
                        ppconds(:, b1.behav.cond( 4, :)==2 & b1.behav.cond( 3, :)==2)  ];
            
            % set some handy variables
            testlabs        =   cat(2, labels_test{:});
            con_p_trial     =   cell(1,length( testlabs));
            
            % iterate over test pairs
            for n_te = 1:length( testlabs)
                for n_lab = 1:2
                    
                    this_lab = testlabs{n_te}{n_lab};
                    
                    dis = find( cellfun( @(x) contains( this_lab, x), {'far','near'}));
                    obj = find( cellfun( @(x) contains( this_lab, x), {'mel','box'}));
                    
                    % select trials to include
                    sel_trials  =   ones( 1, size( ppconds, 2));
                    if ~isempty(dis)
                        sel_trials(1, ppconds( 4, :) ~= dis)  =   0;
                    end
                    if ~isempty(obj)
                        sel_trials(1, ppconds( 3, :) ~= obj)  =   0;
                    end
                    sel_trials = find( sel_trials);
                    
                    % stack data from both labels
                    con_p_trial{n_te} =   [con_p_trial{n_te}; ppconds( [5 6 end], sel_trials)']; %#ok<*FNDSB>
                end
            end
            
            % stack over conditions of non-interest
            for nc = 1:length( myconds)                                           % 1. upright vs inverted scenes
                con_pt(ip).(myconds{nc})  =   cat(1, con_p_trial{wpairs{nc}});  % 2. specific scene 1:8
            end                                                                     % 3. trial number
        end
    end
    
    % set memory/cpu requirement for permutation job
    mem_per     =   55;     % cluster job memory requirement per participant, in gb
    time_per    =   3;      % cluster job walltime requirement per participant, in h
    waitfordata =   3;
    ntestex     =   NaN( 1, nPairs);
    
    % iterate over classification pairs and create shuffled data
    for np = 1:nPairs
        
        ntestex(np) = length( data(np).labels);
        
        fprintf( 'Creating permuted data set for classification pair %g/%g.\n', np, nPairs);
        qsubfeval( @permute_data, np, data(np), nsamples, samesvm, paths, def, 'memreq', mem_per*1024^3, 'timreq', time_per*3600);
        
    end % end of loop over classification pairs
    
    % wait for permuted data to be available
    fprintf( '\nRetrieving permuted data. Pairs left: %g...', nPairs);
    old_l = nPairs;
    while 1
        flist = dir( fullfile( paths.res, 'temp', sprintf( '%s_np*', paths.temp)));
        if length( flist) == nPairs
            cleanit;
            break
        else
            if nPairs - length( flist) < old_l
                old_l = nPairs - length( flist);
                fprintf( '%g...', old_l)
                waitfordata = 1;
            else
                pause(waitfordata);
            end
        end
    end
    fprintf( 'Done: all permutation sets have been retrieved.\n');
    
    % iterate over conditions of interest and compute metrics
    for nc = 1:length( myconds)
        
        % retrieve labels, and adapt size/shape to match data
        fulllabs    =   repmat( reshape( cat( 1, data(wpairs{nc}).labels), 1, 1, sum( ntestex(wpairs{nc}))), nSubjects, nrois, 1, nsamples+1);
        thesepairs  =   find( wpairs{nc});
        
        % recreate accuracy data from (un)shuffled distance to bound
        if isempty( def) || any( contains( def.metrics, {'acc'}))
            
            for ni = 1:length( thesepairs)
                np              =   thesepairs(ni);
                tempath         =   fullfile( paths.res, 'temp', sprintf( '%s_np%g.mat', paths.temp, np));
                
                % give several attempts to load file.. (might be writing still)
                while 1
                    try
                        tempdata        =   load( tempath, 'd_data');
                        break
                    catch
                        disp( 'Mat-file is still being written, trying again shortly.')
                        pause( 10);
                    end
                end
                
                d_data(np).dist =   tempdata.d_data.dist;
                clearvars tempdata
            end
            
            fulldata            =   cat(3, d_data(wpairs{nc}).dist);
            s.acc.(myconds{nc}) =   squeeze( 100* mean( fulldata .* fulllabs < 0, 3)) -50;
            clearvars d_data
        end
        
        % compute combined z-scored distance to bound measure
        if isempty( def) || any( contains( def.metrics, {'zdb','mfc'}))
            
            for ni = 1:length( thesepairs)
                np              =   thesepairs(ni);
                tempath         =   fullfile( paths.res, 'temp', sprintf( '%s_np%g.mat', paths.temp, np));
                
                % give several attempts to load file.. (might be writing still)
                while 1
                    try
                        tempdata        =   load( tempath, 'z_data');
                        break
                    catch
                        disp( 'Mat-file is still being written, trying again shortly.')
                        pause( 10);
                    end
                end
                
                z_data(np).dist =   tempdata.z_data.dist;
                clearvars tempdata
            end
            
            if isempty( def) || any( contains( def.metrics, {'zdb'}))
                fulldata            =   cat( 3, z_data(wpairs{nc}).dist);
                s.zdb.(myconds{nc}) =   squeeze( mean( fulldata .* -fulllabs, 3));
                
                if save_dpt     % trial based data
                    t.zdb.(myconds{nc}) =     fulldata(:,:,:,1) .* -fulllabs(:,:,:,1);
                end
            end
            
            % compute classification score (forced, median-based)
            if isempty( def) || any( contains( def.metrics, {'mfc'}))
                s.mfc.(myconds{nc}) =   squeeze( 100* mean( fulldata .* fulllabs < 0, 3)) -50;
            end
            clearvars z_data
        end
        
        % compute classification score (forced, pairwise)
        if isempty( def) || any( contains( def.metrics, {'pfc'}))
            
            for ni = 1:length( thesepairs)
                np              =   thesepairs(ni);
                tempath         =   fullfile( paths.res, 'temp', sprintf( '%s_np%g.mat', paths.temp, np));
                
                % give several attempts to load file.. (might be writing still)
                while 1
                    try
                        tempdata        =   load( tempath, 'f_data');
                        break
                    catch
                        disp( 'Mat-file is still being written, trying again shortly.')
                        pause( 10);
                    end
                end
                
                f_data(np).dist =   tempdata.f_data.dist;
                clearvars tempdata
            end
            
            fulldata            =   cat(3, f_data(wpairs{nc}).dist);
            s.pfc.(myconds{nc}) =   squeeze( 100* mean( fulldata < 0, 3)) -50;
            clearvars f_data
        end
        
        % compute combined rank-based continuous measure
        if isempty( def) || any( contains( def.metrics, {'rdb'}))
            
            for ni = 1:length( thesepairs)
                np              =   thesepairs(ni);
                tempath         =   fullfile( paths.res, 'temp', sprintf( '%s_np%g.mat', paths.temp, np));
                
                % give several attempts to load file.. (might be writing still)
                while 1
                    try
                        tempdata        =   load( tempath, 'r_data');
                        break
                    catch
                        disp( 'Mat-file is still being written, trying again shortly.')
                        pause( 10);
                    end
                end
                
                r_data(np).dist =   tempdata.r_data.dist;
                clearvars tempdata
            end
            
            fulldata            =   cat( 3, r_data(wpairs{nc}).dist);
            s.rdb.(myconds{nc}) =   squeeze( sum( fulldata .* -fulllabs, 3));
            
            if save_dpt     % trial based data
                t.rdb.(myconds{nc}) =     fulldata(:,:,:,1) .* -fulllabs(:,:,:,1);
            end
            clearvars r_data
        end
        
    end % end of condition iteration
        
    % compute derived metrics (e.g., difference and average)
    metrics         =   fields( s);
    for nm = 1 : length( metrics)
        s.(metrics{nm}).difference  =    s.(metrics{nm}).matching - s.(metrics{nm}).mismatch;
        s.(metrics{nm}).average     =   (s.(metrics{nm}).matching + s.(metrics{nm}).mismatch) ./2;
    end

    % do the same for trialwise data
    if save_dpt
        metrics         =   fields( t);
        for nm = 1 : length( metrics)
        t.(metrics{nm}).difference  =    t.(metrics{nm}).matching - t.(metrics{nm}).mismatch;
        t.(metrics{nm}).average     =   (t.(metrics{nm}).matching + t.(metrics{nm}).mismatch) ./2;
        end
    end
    
    % iterate over main ROI regions (e.g., EVC, LOC, etc.)
    for nrt = 1:length( roi_regs)

        % get current ROI
        gen_roi = roi_regs{nrt};

        % store some general info
        RS.(gen_roi).info.names             =   rois( contains( rois, gen_roi) & ~contains( rois, '_all_'));
        RS.(gen_roi).info.voxels            =   compose( '%g', voxstep(nrt):voxstep(nrt):voxmax(nrt));
        RS.(gen_roi).info.regions           =   gen_roi;
        RS.(gen_roi).info.mytest            =   mytest;
        RS.(gen_roi).info.direction         =   classdir;
        
        RS.(gen_roi).info.labels.test       =   labels_test;
        RS.(gen_roi).info.labels.train      =   labels_train;
        RS.(gen_roi).info.labels.pairs      =   conds;
        RS.(gen_roi).info.labels.conditions =   myconds;
        
        RS.(gen_roi).info.stats.sided       =   sided;
        RS.(gen_roi).info.stats.cutoff      =   cutoff;
        RS.(gen_roi).info.stats.ttype       =   ttype;
        RS.(gen_roi).info.stats.corr        =   ctype;
        RS.(gen_roi).info.stats.nsamples    =   nsamples;
        
        % retrieve indices of current ROIs
        allroi = contains( rois, gen_roi) &  contains( rois, '_all_');
        voxroi = contains( rois, gen_roi) & ~contains( rois, '_all_');

        % assign data to correct ROI in struct
        metrics         =   fields( s);
        for nm = 1 : length( metrics)
            conditions     =   fields( s.(metrics{nm}));
            for nc = 1 : length( conditions)
                RS.(gen_roi).data.vox.(metrics{nm}).(conditions{nc}) =  s.(metrics{nm}).(conditions{nc})(:,voxroi,:);
                RS.(gen_roi).data.all.(metrics{nm}).(conditions{nc}) =  s.(metrics{nm}).(conditions{nc})(:,allroi,:);
            end
        end
        
        % do the same for trial based data
        if save_dpt
            
            % save info
            CT.(gen_roi).info = RS.(gen_roi).info; 
            
            % assign data to correct ROI in struct
            metrics         =   fields( t);
            for nm = 1 : length( metrics)
                conditions     =   fields( t.(metrics{nm}));
                for nc = 1 : length( conditions)
                    CT.(gen_roi).data.vox.(metrics{nm}).(conditions{nc}) =  t.(metrics{nm}).(conditions{nc})(:,voxroi,:);
                    CT.(gen_roi).data.all.(metrics{nm}).(conditions{nc}) =  t.(metrics{nm}).(conditions{nc})(:,allroi,:);
                end
            end
        end
        
    end % end of loop over main ROI regions

    % load behavioral results per participant and add to RS struct
    load( paths.behav);
    RS.behav.labels    =   flabels;
    RS.behav.b_data    =   all_data(subjects,:); %#ok<IDISVAR,NODEF>
    
    % average all data from complementart ROIs
    u_roinames  =   unique( roinames, 'stable');
    
    % iterate over grouped ROI names
    for ngroi = 1 : length( u_roinames)
        
        % retrieve ROI names and consituent ROIs
        gr_roi  =   u_roinames{ngroi};
        cs_rois =   roi_regs( contains( roinames, gr_roi));
        
        % copy entire sub-struct (including info-field, etc.)
        RS.(gr_roi) =   RS.(cs_rois{1});
        nvoxgr      =   fields( RS.(gr_roi).data);
        
        % do the same for trial based data
        if save_dpt
            CT.(gr_roi) =   CT.(cs_rois{1});
        end
        
        % iterate over voxel groupings, ..
        for nd = 1 : length( nvoxgr)
            
            metrics = fields( RS.(gr_roi).data.(nvoxgr{nd}));
            
            % metrics
            for nm = 1 : length( metrics)
                
                conditions = fields( RS.(gr_roi).data.(nvoxgr{nd}).(metrics{nm}));
                
                % and conditions
                for nc = 1 : length( conditions)
                    
                    % iterate over ROI constituents
                    for nmat = 2: length( cs_rois)
                        
                        RS.(gr_roi).data.(nvoxgr{nd}).(metrics{nm}).(conditions{nc})  =   RS.(gr_roi).data.(nvoxgr{nd}).(metrics{nm}).(conditions{nc}) + RS.(cs_rois{nmat}).data.(nvoxgr{nd}).(metrics{nm}).(conditions{nc});
                        
                    end % end of constituing ROIs
                    
                    % divide by number of constituents to create average
                    RS.(gr_roi).data.(nvoxgr{nd}).(metrics{nm}).(conditions{nc})  =   RS.(gr_roi).data.(nvoxgr{nd}).(metrics{nm}).(conditions{nc}) / length(cs_rois);
                    
                end % end of condition loop
            end % end of metrics loop
            
            % do the same for trial-based data
            if save_dpt
                
                metrics = fields( CT.(gr_roi).data.(nvoxgr{nd}));
                
                % metrics
                for nm = 1 : length( metrics)
                    
                    conditions = fields( CT.(gr_roi).data.(nvoxgr{nd}).(metrics{nm}));
                    
                    % and conditions
                    for nc = 1 : length( conditions)
                        
                        % iterate over ROI constituents
                        for nmat = 2: length( cs_rois)
                            
                            CT.(gr_roi).data.(nvoxgr{nd}).(metrics{nm}).(conditions{nc})  =   CT.(gr_roi).data.(nvoxgr{nd}).(metrics{nm}).(conditions{nc}) + CT.(cs_rois{nmat}).data.(nvoxgr{nd}).(metrics{nm}).(conditions{nc});
                            
                        end % end of constituing ROIs
                        
                        % divide by number of constituents to create average
                        CT.(gr_roi).data.(nvoxgr{nd}).(metrics{nm}).(conditions{nc})  =   CT.(gr_roi).data.(nvoxgr{nd}).(metrics{nm}).(conditions{nc}) / length(cs_rois);
                        
                    end % end of condition loop
                end % end of metrics loop
                
            end % end of trial-based data loop
            
        end  % end of voxel groupings (vox/all)
    end % end of LR-combined ROIs
      
    % save data in structs
    save( paths.f,     'RS', 'data', '-v7.3');
    if save_dpt
        save( paths.trial, 'CT', 'con_pt');
    end
    
    % clean up to free up memory
    clearvars -except show paths def RS
    delete( fullfile( paths.res, 'temp', sprintf( '%s_np*', paths.temp)));
    cleanit;
    
end % end of 'organize data per main ROI'


%% Stastistics (roi-wise stats, and tfce & cluster-based stats across ROIs)

% retrieve useful data
roi_regs                                =   fields( RS);
roi_regs( strcmp( roi_regs, 'behav'))   =   [];
null_methods                            =   def.null;               % permutation of data (1) or shuffling of labels (2)

% check whether analyses are already conducted or not
if def.runit || ~isfield( RS.(roi_regs{1}), 'stats')
    
    % conduct all statistical tests
    testall;
    
end % end of statistics
        
% clean up
clearvars -except show paths def RS


%% Plot graphs 

% check if plots are requested
if def.plotit==1
    
    % pot all requested graphs and report main results
    plotall3;
    
end % end of plotting

end % end of main function

