%% Second level statistics - across subjects
%specify: ROI name
% testType : 'memory' | not for category, memorability etc atm!

%loads timextime correlation scores per subject
% for Encoding_Recognition correlation of
%1 ) all 2)rem 3)forg 4)rem-forg stimuli

function second_level_stats(subjects,testType,nrep,rsa_tag,roi_tag)
    
    fs = '\';
    baseDir = '..\Data';
    
    %load folder names of subjects
    load('..\subject_list.mat'); %subNames for easier folder accessing
    
    
    %% load ROI first-level stats for each subject and store empirical and permutation t-values in one struct
    %valType = {'all','rem','forg','mem'};

    for phase = 1:2
        groupstats = {};
        groupstats.subnum = [];
        groupstats.channum = [];
        if phase == 1
            phaseTag = 'enc';
            disp(['encoding'])
        elseif phase == 2
            phaseTag = 'rec';
            disp(['recognition'])
        end

        for subject = 1:length(subjects)
            first_lvl_stats_dir = [baseDir fs subNames{subjects(subject)} fs 'stats' fs 'BtwItem' fs rsa_tag fs roi_tag fs 'Btw_within_1D'];

            load([first_lvl_stats_dir fs phaseTag '_btw_within_emp.mat']); emp = emp_stats;
            load([first_lvl_stats_dir fs phaseTag '_btw_within_perm.mat']); perm = perm_stats;
            load([baseDir fs subNames{subjects(subject)} fs 'stats' fs 'BtwItem' fs rsa_tag fs roi_tag fs 'BtwInfo.mat'])

            %add subnum and channels to groupstats file
            groupstats.subnum(end+1,1) = subjects(subject);
            groupstats.channum(end+1,1) = length(info.channels);

            groupstats.sub.emp(subject,:) = emp.t; 
            groupstats.sub.perm(subject,:,:) = squeeze(perm.t(:,:,:)); 

        end

        %% perform second-level ttests across subjects  
        for win = 1:size(groupstats.sub.emp,2) %encoding time

            [~,P,~,STATS] = ttest(squeeze(groupstats.sub.emp(:,win)));

            groupstats.emp.t(win) = STATS.tstat;
            groupstats.emp.p(win) = P;

            for rep = 1:nrep
                [~,P,~,STATS] = ttest(squeeze(groupstats.sub.perm(:,rep,win)));

                groupstats.perm.t(rep,win) = STATS.tstat;
                groupstats.perm.p(rep,win) = P;
            end
        end

        %% Cluster-statistics
        %empirical clusters
        tmap = groupstats.emp.t;
        pmap = groupstats.emp.p;
        groupstats.emp_cluster = extract_clusterstats(tmap,pmap);

        %permutation clusters
        for rep = 1:nrep
            tmap = squeeze(groupstats.perm.t(rep,:,:));
            pmap = squeeze(groupstats.perm.p(rep,:,:));

            groupstats.perm_cluster(rep) = extract_clusterstats(tmap,pmap);

            perm_pos_cluster(rep) = groupstats.perm_cluster(rep).pos.sum;
            perm_neg_cluster(rep) = groupstats.perm_cluster(rep).neg.sum;
        end

        % permutation cluster size distribution and ranking of empirical in
        % surrogate data
        perm_pos_dist = sort(perm_pos_cluster,'ascend');
        perm_neg_dist = sort(perm_neg_cluster,'descend');

        groupstats.emp_cluster.pos.rank = nearest(perm_pos_dist, groupstats.emp_cluster.pos.sum);
        groupstats.emp_cluster.pos.p = 1 - (groupstats.emp_cluster.pos.rank / nrep);

        groupstats.emp_cluster.neg.rank = nearest(perm_neg_dist, groupstats.emp_cluster.neg.sum);
        groupstats.emp_cluster.neg.p = 1 - (groupstats.emp_cluster.neg.rank / nrep);
    %     %% plot clusterstats results
    %     f = plot_clusterstats(groupstats,rsawidth,rsastep,maxT,ROI,testType);

        %% save stats 
        second_level_stats_dir = ['..\STATS' fs 'Btw_within' fs rsa_tag fs '1D' fs testType];
        status = mkdir(second_level_stats_dir);
        
        save([second_level_stats_dir fs roi_tag '_' phaseTag '_stats.mat'],'groupstats');
    end
end

    
    
            
            
            
            
        
        