% retrieve data types
datatypes       =   fields( RS.(roi_regs{1}).data);

% set computation time and memory
mem_tfc         =   50;
time_tfc        =   10;
mem_clu         =   5;
time_clu        =   3;
all_paths       =   [];
all_cpaths      =   [];

% iterate over ROI regions
for nrt = 1:length( roi_regs)
    
    % get current ROI
    gen_roi = roi_regs{nrt};
    
    % retrieve some handy variables
    nsteps      =   length( RS.(gen_roi).info.names);
    nsamples    =   RS.(gen_roi).info.stats.nsamples;
    metrics     =   fields( RS.(gen_roi).data.vox);
    
    % iterate over classification metrics
    for nm = 1 : length( metrics)
        
        % retrieve conditions of interest
        conditions = fields( RS.(gen_roi).data.vox.(metrics{nm}));
        
        % iterate over conditions of interest
        for nc = 1 : length( conditions)
            
            % select data for tfce
            mydata      =   RS.(gen_roi).data.vox.(metrics{nm}).(conditions{nc});
            nSubjects   =   size( mydata, 1);
            
            % set tfce settings (only once)
            if nrt == 1 && ~exist( 'DS', 'var')
                
                % 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
                
                % create cosmo struct
                DS.sa.targets           =   ones( nSubjects, 1);
                DS.sa.chunks            =   (1:nSubjects)';
                
                DS.a.meeg.sample_field  =   'trial';
                DS.a.meeg.samples_type  =   'timelock';
                DS.a.meeg.samples_label =   'rpt';
                
                DS.a.fdim.labels        =   {'time'};
                
                % analysis info (note: method defaults to 'tfce')
                Opt                     =   struct();
                Opt.h0_mean             =   0;
                
                % set variable settings (number of ROI sizes)
                DS.a.fdim.values{1}     =   1:nsteps;
                DS.fa.time              =   DS.a.fdim.values{1};
            end
            
            % load observed data
            DS.samples  =   mydata(:,:,1);
            
            % compute neighborhoods
            if ~exist( 'nbrh', 'var')
                nbrh    =   cosmo_cluster_neighborhood(DS);
            end
            
            % iterate over methods for generating a null distribution
            for nmeth = 1 : length( null_methods)
                
                % determine settings for each method
                switch null_methods{nmeth}
                    case 'perm'
                        
                        % null data from permuted data
                        null_data       =   [];
                        numiter         =   nsamples*10;
                        
                    case 'labs'
                        
                        % null data from shuffled test labels
                        null_data       =   {};
                        data_samples    =   datasample( 2:nsamples+1, nsamples, 'Replace', false);
                        for ni = 1 : length( data_samples)
                            null_data{ni}            =   DS;
                            null_data{ni}.samples    =   mydata(:,:,ni+1);
                        end
                        numiter         =   nsamples;
                        
                end % end of switch statements for different tfce methods
                
                % set null data
                if nmeth > 1
                    Opt(nmeth) = Opt(1);
                end
                Opt(nmeth).null    =   null_data;   % no data, or shuffled data
                Opt(nmeth).niter   =   numiter;     % 10*nsamples or nsamples
                
                % set temp file
                tempath{nmeth}  =   fullfile( paths.res, 'temp', sprintf( '%s_%s_%s_%s_%g.mat', paths.temp, gen_roi, metrics{nm}, conditions{nc}, nmeth));
                all_paths       =   [all_paths; {tempath{nmeth}, gen_roi, metrics{nm}, conditions{nc}, nmeth}];
                
            end % end of tfce methods for generating null distribution
            
            % run tfce analyses (both methods)
            while 1
                try
                    qsubfeval( @compute_tfce2, DS, nbrh, Opt, tempath, 'memreq', mem_tfc*1024^3, 'timreq', time_tfc*3600);
                    break;
                catch
                    warning( 'Job not sent to [qsubfeval], trying again shortly');
                    pause(10);
                end
            end
            
        end % end of conditions loop
    end % end of metrics loop
end % end of ROI loop

% retrieve stats options from RS struct
stats       =   RS.(roi_regs{1}).info.stats;

sided       =   stats.sided;
cutoff      =   stats.cutoff;
ttype       =   stats.ttype;
nsamples    =   stats.nsamples;

stats.null_methods          =   null_methods;
[h0_mean, stats.h0_mean]    =   deal( Opt.h0_mean);

% determine significance level for 1/2-sided tests
thr         =  [cutoff./sided, 1-cutoff./sided];

% iterate through general ROIs (e.g., EVC, LOC, PPA, ...)
for nrt = 1 : length( roi_regs)
    
    % get current ROI name and ROI coint
    gen_roi     =   roi_regs{nrt};
    nvox        =   length( RS.(gen_roi).info.names);
    
    % retrieve classification metrics
    metrics = fields( RS.(gen_roi).data.vox);
    
    % iterate over classification metrics
    for nm = 1 : length( metrics)
        
        % retrieve conditions of interest
        conditions = fields( RS.(gen_roi).data.vox.(metrics{nm}));
        
        % iterate over conditions of interest
        for nc = 1 : length( conditions)
            
            % select data
            voxdata     =   RS.(gen_roi).data.vox.(metrics{nm}).(conditions{nc});
            alldata     =   RS.(gen_roi).data.all.(metrics{nm}).(conditions{nc}); %#ok<NASGU>
            
            % remove permuted data (to free up storage space) once used
            RS.(gen_roi).data.vox.(metrics{nm}).(conditions{nc}) = RS.(gen_roi).data.vox.(metrics{nm}).(conditions{nc})(:,:,1);
            RS.(gen_roi).data.all.(metrics{nm}).(conditions{nc}) = RS.(gen_roi).data.all.(metrics{nm}).(conditions{nc})(:,:,1);
            
            % iterate over 'vox' and 'all' rois
            for nd = 1 : length( datatypes)
                
                % pick data
                obs_data    =   eval( sprintf( '%sdata(:,:,1)',     datatypes{nd}));
                nulldata    =   eval( sprintf( '%sdata(:,:,2:end)', datatypes{nd}));
                
                % choose statistical test
                switch ttype
                    
                    % conduct t-test, report results and store means and CIs
                    case 'ttest'
                        
                        % retrieve statistics
                        [~,pval,ci] =   ttest( obs_data, h0_mean, 'alpha', cutoff*(3-sided));
                        
                        pval        =   pval/(3-sided);
                        mu          =   mean( obs_data);
                        
                        
                        % or conduct bootstrap test, and store means and CIs
                    case 'bootstrap'
                        
                        % randomly sample obs_data
                        mus = NaN( nsamples, size( obs_data, 2));
                        for ns = 1 : nsamples
                            mus(ns, :)     =   mean( datasample( obs_data, size( obs_data, 1)));
                        end
                        
                        % retrieve statistic
                        pval    =   min( [mean( mus<=h0_mean); 1-mean( mus<=h0_mean)], [], 1)*sided; % note: different predictions for LOC and EVC (alpha level is accounted for during permutation test)
                        mu      =   mean( mus);
                        ci      =   quantile( mus, thr);
                        
                        
                        % or conduct signed rank test
                    case 'rank'
                        
                        % retrieve rank stats (per roi)
                        pval = NaN( 1, size( obs_data, 2));
                        for nr = 1:size( obs_data, 2)
                            pval(nr)    =   signrank( obs_data(:,nr), h0_mean);
                        end
                        pval = pval/(3-sided);
                        
                        % retrieve bootstrapped mu and ci
                        mus = NaN( nsamples, size( obs_data, 2));
                        for ns = 1 : nsamples
                            mus(ns, :)     =   mean( datasample( obs_data, size( obs_data, 1)));
                        end
                        
                        ci      =   quantile( mus, thr);
                        mu      =   mean( mus);
                        
                end % end of switch statement (ttype)
                
                % stack stats
                RS.(gen_roi).stats.(datatypes{nd}).(metrics{nm}).(conditions{nc}).pvals  =   pval;
                RS.(gen_roi).stats.(datatypes{nd}).(metrics{nm}).(conditions{nc}).ci     =   ci;
                RS.(gen_roi).stats.(datatypes{nd}).(metrics{nm}).(conditions{nc}).mean   =   mu;
                
                % add null distribution based on shuffled data
                null_mus    =   squeeze( mean( nulldata, 1))';
                
                RS.(gen_roi).null.(datatypes{nd}).(metrics{nm}).(conditions{nc}).mean   =   mean( null_mus);
                RS.(gen_roi).null.(datatypes{nd}).(metrics{nm}).(conditions{nc}).ci     =   quantile( null_mus, thr);
                RS.(gen_roi).null.(datatypes{nd}).(metrics{nm}).(conditions{nc}).pvals  =   min( [mean( null_mus<=mu); 1-mean( null_mus<=mu)], [], 1)*sided;
                
            end % end of datatypes loop (vox/all)
            
            
            % cluster stats %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            % get real data (only for data per vox-count)
            real_p      =   RS.(gen_roi).stats.vox.(metrics{nm}).(conditions{nc}).pvals;
            
            % retrieve sign (to separate clusters by sign)
            it_sign         =       (mean( voxdata(:,:,1), 1) - h0_mean) > 0;
            ss              =       zeros( size( real_p));
            ss(strfind( [0,it_sign(:)'],[0 1])) = 1;
            ss(strfind( [0,it_sign(:)'],[1 0])) = 1;
            
            % find cluster sizes and locations in observed data
            ii              =       zeros( size( real_p));
            jj              =       real_p < cutoff;            % all sub threshold rois
            ii(strfind( [0,jj(:)'],[0 1])) = 1;                 % onset of sub threshold rois
            ii              =        ii | (ss & jj);            % taking into account sign changes
            real_idx        =       cumsum(ii).*jj;             % and number the cluster incrementally
            
            % stack cluster sizes of each cluster
            if any( real_idx)
                real_clus   =       cellfun( @length, accumarray( real_idx(jj)', real_p(jj)', [], @(x){x'}));
            else
                real_clus   =       [];
            end
            
            % store cluster indices and cluster number
            RS.(gen_roi).stats.vox.(metrics{nm}).(conditions{nc}).clus.all_k    =   real_idx;
            
            % iterate over methods for generating a null distribution
            for nmeth = 1 : length( null_methods)
                
                % set temp file
                tempath{nmeth}  =   fullfile( paths.res, 'temp', sprintf( '%s_clus_%s_%s_%s_%g.mat', paths.temp, gen_roi, metrics{nm}, conditions{nc}, nmeth));
                all_cpaths      =   [all_cpaths; {tempath{nmeth}, gen_roi, metrics{nm}, conditions{nc}, nmeth}];
            end
            
            % run cluster permutation analyses (both methods)
            while 1
                try
                    qsubfeval( @clusterstats, tempath, voxdata, real_idx, real_clus, stats, 'memreq', mem_clu*1024^3, 'timreq', time_clu*3600);
                    break;
                catch
                    warning( 'Job not sent to [qsubfeval], trying again shortly');
                    pause(10);
                end
            end
            
        end % end of conditions of interest
    end %end of classification metrics
end % end of ROI loop


% wait for all tfce result files to arrive from cluster
waitfordata     =   5;
fprintf( '\nRetrieving tfce results from cluster ...\n\t> Data sets remaining: \n\t%g... ', size( all_paths, 1))

% repeatedly check for new files
while ~isempty( all_paths)
    
    % iterate over (missing) datafiles
    for nf = size( all_paths, 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( all_paths{nf,1}, 'file')
            
            % load and stack data
            while 1
                try
                    load( all_paths{nf,1}, 'pvals');
                    break
                catch
                    disp( 'Mat-file is still being written, trying again shortly.')
                    pause( 1);
                end
            end
            RS.(all_paths{nf,2}).stats.vox.(all_paths{nf,3}).(all_paths{nf,4}).tfce.(null_methods{all_paths{nf,5}}) = pvals;
            
            % delete temporary data
            delete( all_paths{nf,1});
            all_paths(nf,:) = [];
            
            % reduce waiting time and report progress
            waitfordata = 0.5;
            fprintf( '%g... ', size( all_paths, 1));
            if mod( size( all_paths, 1), 10)==0
                fprintf( '\n\t');
            end
        end
    end
    
    % wait N seconds before checking for new data again
    pause(waitfordata);
end

% save classification results
fprintf( '\n\t> TFCE results collected, and saved in "stats" struct.\n\n');

% wait for all tfce result files to arrive from cluster
waitfordata     =   5;
fprintf( '\nRetrieving cluster results from cluster ...\n\t> Data sets remaining: \n\t%g... ', size( all_cpaths, 1))

% repeatedly check for new files
while ~isempty( all_cpaths)
    
    % iterate over (missing) datafiles
    for nf = size( all_cpaths, 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( all_cpaths{nf,1}, 'file')
            
            % load and stack data
            load( all_cpaths{nf,1}, 'clus_idx', 'cl_prob', 'min_clus', 'k_mean');
            
            % store all cluster-based analysis outcomes in struct
            RS.(all_cpaths{nf,2}).stats.vox.(all_cpaths{nf,3}).(all_cpaths{nf,4}).clus.(null_methods{all_cpaths{nf,5}}).sig_k    =   clus_idx;   % indices of all significant clusters
            RS.(all_cpaths{nf,2}).stats.vox.(all_cpaths{nf,3}).(all_cpaths{nf,4}).clus.(null_methods{all_cpaths{nf,5}}).k_prob   =   cl_prob;    % probabilities of all cluster sizes
            RS.(all_cpaths{nf,2}).stats.vox.(all_cpaths{nf,3}).(all_cpaths{nf,4}).clus.(null_methods{all_cpaths{nf,5}}).k_min    =   min_clus;   % minimal cluster size for significance
            RS.(all_cpaths{nf,2}).stats.vox.(all_cpaths{nf,3}).(all_cpaths{nf,4}).clus.(null_methods{all_cpaths{nf,5}}).k_mean   =   k_mean;     % average classification accuracy in significant cluster
            
            % delete temporary data
            delete( all_cpaths{nf,1});
            all_cpaths(nf,:) = [];
            
            % reduce waiting time and report progress
            waitfordata = 0.5;
            fprintf( '%g... ', size( all_cpaths, 1));
            if mod( size( all_cpaths, 1), 10)==0
                fprintf( '\n\t');
            end
        end
    end
    
    % wait N seconds before checking for new data again
    pause(waitfordata);
end

% save classification results
fprintf( '\n\t> Cluster results collected, and saved in "stats" struct.\n\n');

% save tfce statistics
load( paths.f, 'data');
save( paths.f, 'RS', 'data', '-v7.3');

% clean up
cleanit;