function dccn_univar2( varargin)

% WHAT:     Simple univariate analyses (averaged beta maps per condition)
% INPUTS:   [shape] or [size], [f]orward or [b]ackward train<->test runs
% FILE:     uni_[trainmodel]_[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;                      % decide whether to plot results at all or not.   
show.rois   =   [];                     % 'LOC','EVC','PPA','TOS','RSC','FFA'
show.hemi   =   {'LOC', 'EVC'}; 
show.wplots =   [13, 14];               %[1:4, 13, 14] % listed on line 101
show.fixedy =   1;
show.saveit =   1;

% select which metric to report
show.null   =   'perm';             % 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)

% subjects
subjects    =   1:24;

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

% test name
mytest      =   'uni';                  

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

% 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)

% 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
labels_train    =   {'mel_near', 'box_near', 'mel_far',  'box_far' };
labels_test     =   labels_train;

models          =   {'train', 'test'};

conds           =   {'mel_near',  [1 0 0 0]; ...
                     'box_near',  [0 1 0 0]; ...
                     'mel_far',   [0 0 1 0]; ...
                     'box_far',   [0 0 0 1]; ...
                     'near',      [1 1 0 0]; ...
                     'far',       [0 0 1 1]; ...
                     'mel',       [1 0 1 0]; ...
                     'box',       [0 1 0 1]} ;
                 
difconds        =   {'shape_near', [2 1 0 0 0 0 0 0]; ...      % mel_near - box_near
                     'shape_far',  [0 0 2 1 0 0 0 0]; ...      % mel_far  - box_far
                     'size_mel',   [1 0 2 0 0 0 0 0]; ...      % mel_near - mel_far
                     'size_box',   [0 1 0 2 0 0 0 0]; ...      % box_near - box-far
                     'shape',      [0 0 0 0 0 0 2 1]; ...
                     'size',       [0 0 0 0 1 2 0 0]} ;

intconds        =   {'shape_x',    [1 2 0 0 0 0]; ...
                     'size_x',     [0 0 1 2 0 0]; ...
                     'main_x',     [0 0 0 0 1 2]};
                 
                 
%% 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 =   100;
   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);
            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{1,2});
nrois       =   length( rois);
myconds     =   {conds{:,1}};

% determine or create output directory (group level)
paths.res  =   fullfile( paths.local, 'univariate_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('%s_hemi_%s_%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('%s_hemi_%s_%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
sn = 'a':'z';
paths.temp = sn(randi(numel( sn), 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 conditions
        f2      =   dir( paths.sim);
        if isempty( f2)
            def.runit = 1; % extract betas nonetheless if no matching file is found
        else
            paths.f2  =   fullfile( paths.res, f2(1).name); 
            load( paths.f2, 'data');
        end
    end
end        


%% Extract betas

% 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( genpath( fullfile( tdt_ver.folder, tdt_ver(end).name))); 
    addpath( '/home/common/matlab/spm12');
    
    % check whether cosmo mvpa path is set (and add path if not)
    if ~exist('cosmo_wtf.m','file')
        cd( paths.cosmo);
        cosmo_set_path();
        savepath();
        cd( paths.proj)
    end
    
    % 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};
            
            for nroi = 1:length( rois)
                check_it(nroi,nsub)  =   contains( sprintf( '%s.nii', rois{nroi}), name_list);
            end
        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
        
    
    % Start iteration over beta extraction jobs %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fprintf( '\nExtracting betas ...\n');
    tempaths = [];
    
    % set memory/cpu requirement for cluster job
    mem_ext   =   1.5;      % cluster job memory requirement per participant, in gb
    time_ext  =   3;        % cluster job walltime requirement per participant, in h
    
    % create empty struct for data storage
    data.train  =   NaN( nSubjects, nrois, length( labels_train));
    data.test   =   NaN( nSubjects, nrois, length( labels_test));
    data.rois   =   rois;
    
    % iterate over subjects
    for n_sub = 1 : nSubjects
        
        % set current subject info
        subj            =   sprintf( '%02g', subjects(n_sub));
        subject_dir     =   fullfile( paths.subj,  subj);
        sub_roi_dir     =   fullfile( subject_dir, roi_dir);
                    
        % iterate over models
        for n_mod = 1: length( models)
            
            % set current model
            model           =   eval( sprintf( 'model_%s',           models{n_mod}));
            labels          =   eval( sprintf( 'labels_%s',          models{n_mod}));
            beta_dir        =   fullfile( subject_dir, 'univariate', model);
            
            % determine regressor names (run-based or trial-based models)
            if strcmp( 'st_', model(1:3))
                
                load( fullfile( beta_dir, 'betas', sprintf( '%s_beta_info.mat', subj)));
                
                regressor_names =   trialinfo(2:end, [6, 2, 8, 4])'; %#ok<*FNCOLND>
                beta_dir        =   cellfun( @(x) fullfile( beta_dir, 'betas', x), regressor_names(3,:), 'UniformOutput', 0)';
            else
                regressor_names =   design_from_spm( beta_dir);
            end
            
            % change label names to match model
            if any( contains( beta_dir, 'test')) && all( contains( labels{n_lab}, '_'))
                labels =  cellfun( @(x) [x(strfind(x,'_')+1:end), '_', x(1:strfind(x,'_')-1)], labels, 'UniformOutput', 0);
            end
            
            % load SPM file, list beta maps
            load( fullfile( beta_dir, 'SPM.mat'), 'SPM');
                            
            % iterate over sets
            for n_lab = 1 : length( labels)
                
                % list beta maps
                maps        =   find( contains( regressor_names(1,:), labels{n_lab}));
                nmaps       =   length( maps);
                
                % set tempfile
                tempath     =  fullfile( paths.res, 'temp', sprintf( '%s-%s-%s-%g.mat', paths.temp, subj, models{n_mod}, n_lab));
                tempaths    =  [tempaths; {tempath, n_sub, models{n_mod}, n_lab}];
                
                % run job
                qsubfeval( @get_uni, beta_dir, sub_roi_dir, rois, maps, tempath, 'memreq', mem_ext*1024^3, 'timreq', time_ext*3600);
                
            end % end of label loop
        end % end of model loop
    end % end of subject loop
    
    
    % retrieve univariate results
    fprintf( '\nBetas exctracted. Retrieving from cluster: ... ');
    
    % repeatedly check for new files
    waitfordata = 5;
    while ~isempty( tempaths)
        
        % iterate over (missing) datafiles
        for nf = size( tempaths, 1):-1:1   %NB: iterate backwards to circumvent position shifts due to deletions
            
            % check wether temp file was created for current roi/metric
            if exist( tempaths{nf,1}, 'file')
                
                % load and stack data
                while 1
                    try
                        load( tempaths{nf,1}, 'meanvals');
                        break
                    catch
                        disp( 'Mat-file is still being written, trying again shortly.')
                        pause( 1);
                    end
                end
                
                % store all extracted betas in struct
                data.(tempaths{nf,3})(tempaths{nf,2}, :, tempaths{nf,4}) = meanvals;
                
                % delete temporary data
                delete( tempaths{nf,1});
                tempaths(nf,:) = [];
                
                % reduce waiting time and report progress
                waitfordata = 0.5;
                fprintf( '%g... ', size( tempaths, 1));
                if mod( size( tempaths, 1), 10)==0
                    fprintf( '\n\t');
                end
            end
        end
        
        % wait N seconds before checking for new data again
        pause(waitfordata);
    end
    
    % save univariate results
    fprintf( '\Betas retrieved, and saved to results folder.');
    save( paths.f, 'data');
end

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


%% Create conditions and organize data per ROI

if def.runit || ~exist( 'RS', 'var') 
    
    % preallocate data matrices/structs
    ts  =   struct();
    RS  =   struct();

    % compute derived metrics
    for nm = 1 : length( models)
        for nc = 1 : size( conds, 1)
            ts.(models{nm}).(conds{nc,1})       =   mean( data.(models{nm})(:, :, logical( conds{nc, 2})), 3);
        end
        for nd = 1 : size( difconds, 1)
            ts.(models{nm}).(difconds{nd,1})    =   ts.(models{nm}).(conds{difconds{nd, 2}==1}) - ts.(models{nm}).(conds{difconds{nd, 2}==2});
        end
        for ne = 1 : size( intconds, 1)
            ts.(models{nm}).(intconds{ne,1})    =   ts.(models{nm}).(difconds{intconds{ne, 2}==1}) - ts.(models{nm}).(difconds{intconds{ne, 2}==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.labels.test       =   labels_test;
        RS.(gen_roi).info.labels.train      =   labels_train;
        RS.(gen_roi).info.labels.conditions =   conds;
        RS.(gen_roi).info.labels.diffscores =   difconds;
        
        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.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
        for nm = 1 : length( models)
            conditions     =   fields( ts.(models{nm}));
            for nc = 1 : length( conditions)
                RS.(gen_roi).data.vox.(models{nm}).(conditions{nc}) =  ts.(models{nm}).(conditions{nc})(:,voxroi,:);
                RS.(gen_roi).data.all.(models{nm}).(conditions{nc}) =  ts.(models{nm}).(conditions{nc})(:,allroi,:);
            end
        end
        
    end % end of loop over main ROI regions

    % 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);
        
        % iterate over voxel groupings, ..
        for nd = 1 : length( nvoxgr)
            
            models = fields( RS.(gr_roi).data.(nvoxgr{nd}));
            
            % models
            for nm = 1 : length( models)
                
                conditions = fields( RS.(gr_roi).data.(nvoxgr{nd}).(models{nm}));
                
                % and conditions
                for nc = 1 : length( conditions)
                    
                    % iterate over ROI constituents
                    for nmat = 2: length( cs_rois)
                        
                        RS.(gr_roi).data.(nvoxgr{nd}).(models{nm}).(conditions{nc})  =   RS.(gr_roi).data.(nvoxgr{nd}).(models{nm}).(conditions{nc}) + RS.(cs_rois{nmat}).data.(nvoxgr{nd}).(models{nm}).(conditions{nc});
                        
                    end % end of constituing ROIs
                    
                    % divide by number of constituents to create average
                    RS.(gr_roi).data.(nvoxgr{nd}).(models{nm}).(conditions{nc})  =   RS.(gr_roi).data.(nvoxgr{nd}).(models{nm}).(conditions{nc}) / length(cs_rois);
                    
                end % end of condition loop
            end % end of models 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');
    
    % 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                            =   {'perm'};               % 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
    plotuni2;
    
end % end of plotting

end % end of main function
