

function dnn_second_level_btw(subjects,dnnType,corr_type,rsa_tag,roi_tag,nrep)
    dataDir = '..\Data';
    fs = '\';
        
    %load folder names of subjects
    load('..\subject_list.mat'); %subNames for easier folder accessing

    %% DNN specifics
    switch dnnType
        case 'AlexNet'
            layerNum = 8;
            layerNames = {'conv1','conv2','conv3','conv4','conv5','fc6','fc7','fc8'}; %%% ADD FOR OTHERS TOO!
        case 'AlexNet_avg'
            layerNum = 8;
            layerNames = {'conv1','conv2','conv3','conv4','conv5','fc6','fc7','fc8'}; %%% ADD FOR OTHERS TOO!
        case 'VGG16'
            layerNum = 16;
            layerNames = {'conv1_1','conv1_2','conv2_1','conv2_2','conv3_1','conv3_2','conv3_3','conv4_1','conv4_2','conv4_3','conv5_1','conv5_2','conv5_3','fc6', 'fc7', 'fc8a'};
        case 'VGG16_avg'
            layerNum = 16;
            layerNames = {'conv1_1','conv1_2','conv2_1','conv2_2','conv3_1','conv3_2','conv3_3','conv4_1','conv4_2','conv4_3','conv5_1','conv5_2','conv5_3','fc6', 'fc7', 'fc8a'};
        case 'Google'
            layerNum = 1;
            layerNames = {'GoogleEncoder'};
        case 'GoogleSemantic'
            layerNum = 1;
            layerNames = {'GoogleSemantic'};
    end

    %set colormap
    cmap = flipud(cbrewer('div','RdBu',51)); cmap(cmap>1) = 1; cmap(cmap<0) = 0;
    statsDir = ['stats' fs 'BtwItemDNN' fs rsa_tag fs roi_tag];
    
    %% COMPUTE THE EMPIRICAL STATS
    %iterate over subj
    for phase = 1:2
        if phase == 1
            ph = 'enc';
            phTag = 'Encoding';
            disp(phTag)
        else
            ph = 'rec';
            phTag = 'Recognition';
            disp(phTag)
        end
        
        %init groupstats info file
        groupstats = {};
        groupstats.subs = subjects;
        groupstats.chanNum = [];
        
        groupstats.sub.emp.all.t = [];
        groupstats.sub.emp.rem.t = [];
        groupstats.sub.emp.forg.t = [];
        groupstats.sub.emp.mem.t = [];
        
        groupstats.sub.perm.all.t = [];
        groupstats.sub.perm.rem.t = [];
        groupstats.sub.perm.forg.t = [];
        groupstats.sub.perm.mem.t = [];

        %load subjects empirical and permutation data
        for num = 1:length(subjects)
            subject = subjects(num);
            %set folderName of current subject
            fileName = subNames{subject}; 

            %load raw correlation data for empirical and surrogate
            %
            load([dataDir fs fileName fs statsDir fs 'Item' dnnType '_' corr_type '_' ph '.mat'])
            subInfo = stats.info;
            
            %add num of chans of current subject
            groupstats.chanNum(num,1) = length(subInfo.channels);
                        
            %empirical correlations
            groupstats.sub.emp.all.t(num,:,:) = stats.emp.all;      %all trials against 0
            groupstats.sub.emp.rem.t(num,:,:) = stats.emp.rem;      %remembered trials against 0
            groupstats.sub.emp.forg.t(num,:,:) = stats.emp.forg;    %forgotten against 0
            groupstats.sub.emp.mem.t(num,:,:) = stats.emp.mem;      %remembered vs forgotten 
                    
            %permutation correlations
            groupstats.sub.perm.all.t(num,:,:,:) = stats.perm.all;      %all trials against 0
            groupstats.sub.perm.rem.t(num,:,:,:) = stats.perm.rem;      %remembered trials against 0
            groupstats.sub.perm.forg.t(num,:,:,:) = stats.perm.forg;    %forgotten against 0
            groupstats.sub.perm.mem.t(num,:,:,:) = stats.perm.mem;      %remembered vs forgotten 
        end

        %init struct fields for t-test of correlations across subjects
        groupstats.emp.all.t = []; groupstats.emp.all.p = [];
        groupstats.emp.rem.t = []; groupstats.emp.rem.p = [];
        groupstats.emp.forg.t = [];groupstats.emp.forg.p = [];
        groupstats.emp.mem.t = []; groupstats.emp.mem.p = [];
        
        groupstats.perm.all.t = []; groupstats.perm.all.p = [];
        groupstats.perm.rem.t = []; groupstats.perm.rem.p = [];
        groupstats.perm.forg.t = []; groupstats.perm.forg.p = [];
        groupstats.perm.mem.t = []; groupstats.perm.mem.p = [];
        
        %% perform t-tests across subjects for empirical and permutation data at each TimeXLayer point
        disp('perform t-tests across subjects')
        for win = 1:size(stats.emp.all,1)
            for layer = 1:layerNum
                %all
                subs_data = squeeze(groupstats.sub.emp.all.t(:,win,layer));
                [~,P,~,cluster_stats] = ttest(subs_data);
                groupstats.emp.all.t(win,layer) = cluster_stats.tstat;
                groupstats.emp.all.p(win,layer) = P;

                %rem
                subs_data = squeeze(groupstats.sub.emp.rem.t(:,win,layer));
                [~,P,~,cluster_stats] = ttest(subs_data);
                groupstats.emp.rem.t(win,layer) = cluster_stats.tstat;
                groupstats.emp.rem.p(win,layer) = P;

                %forg
                subs_data = squeeze(groupstats.sub.emp.forg.t(:,win,layer));
                [~,P,~,cluster_stats] = ttest(subs_data);
                groupstats.emp.forg.t(win,layer) = cluster_stats.tstat;
                groupstats.emp.forg.p(win,layer) = P;

                %memory (remembered - forgotten)
                subs_data = squeeze(groupstats.sub.emp.mem.t(:,win,layer));
                [~,P,~,cluster_stats] = ttest(subs_data);
                groupstats.emp.mem.t(win,layer) = cluster_stats.tstat;
                groupstats.emp.mem.p(win,layer) = P;
                
                %repeat for each permutation iteration
                for rep = 1:nrep
                    %all
                    subs_data = squeeze(groupstats.sub.perm.all.t(:,rep,win,layer));
                    [~,P,~,cluster_stats] = ttest(subs_data);
                    groupstats.perm.all.t(rep,win,layer) = cluster_stats.tstat;
                    groupstats.perm.all.p(rep,win,layer) = P;
                    %rem
                    subs_data = squeeze(groupstats.sub.perm.rem.t(:,rep,win,layer));
                    [~,P,~,cluster_stats] = ttest(subs_data);
                    groupstats.perm.rem.t(rep,win,layer) = cluster_stats.tstat;
                    groupstats.perm.rem.p(rep,win,layer) = P;
                    %forg
                    subs_data = squeeze(groupstats.sub.perm.forg.t(:,rep,win,layer));
                    [~,P,~,cluster_stats] = ttest(subs_data);
                    groupstats.perm.forg.t(rep,win,layer) = cluster_stats.tstat;
                    groupstats.perm.forg.p(rep,win,layer) = P;
                    %memory (remembered - forgotten)
                    subs_data = squeeze(groupstats.sub.perm.mem.t(:,rep,win,layer));
                    [~,P,~,cluster_stats] = ttest(subs_data);
                    groupstats.perm.mem.t(rep,win,layer) = cluster_stats.tstat;
                    groupstats.perm.mem.p(rep,win,layer) = P;
                end
            end
        end

        %% save results
        groupstatsDir = ['..\STATS\ItemDNN' fs dnnType '_' corr_type fs rsa_tag fs roi_tag];
        status = mkdir(groupstatsDir);
        save([groupstatsDir fs 'groupstats_' ph '.mat'],'groupstats') 
        
        %% perform cluster statistics separately for each layer and trialType
        types = {'all','rem','forg','mem'};
        for trialtype = 1:4 %all, rem, forg, mem
            field =types{trialtype};
            clusterstats = {};
            for layer = 1:layerNum
                tmap = groupstats.emp.(field).t(:,layer);
                pmap = groupstats.emp.(field).p(:,layer);
                emp_cluster = extract_clusterstats(tmap,pmap);

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

                    perm_cluster = extract_clusterstats(tmap,pmap);

                    perm_pos_cluster(rep) = perm_cluster.pos.sum;
                    perm_neg_cluster(rep) = perm_cluster.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');
                
                if emp_cluster.pos.sum ~= 0
                    emp_cluster.pos.rank = nearest(perm_pos_dist,emp_cluster.pos.sum);
                    emp_cluster.pos.p = 1 - (emp_cluster.pos.rank / nrep);
                else
                    emp_cluster.pos.rank = NaN;
                    emp_cluster.pos.p = NaN;
                end
                if emp_cluster.neg.sum ~= 0
                    emp_cluster.neg.rank = nearest(perm_neg_dist, emp_cluster.neg.sum);
                    emp_cluster.neg.p = 1 - (emp_cluster.neg.rank / nrep);
                else
                    emp_cluster.neg.rank = NaN;
                    emp_cluster.neg.p = NaN;
                end

                clusterstats.([layerNames{layer}]) = emp_cluster;
            end

            %create a map indicating all sig clusters!
            tempmat = zeros(size(groupstats.sub.emp.(field).t,2),layerNum);

            for layer = 1:layerNum
                %look for pos and neg cluster
                cs = clusterstats.(layerNames{layer});
                if cs.pos.sum ~= 0 && cs.pos.p < 0.05
                    tempmat(find(cs.pos.map == cs.pos.idx),layer) = 1;
                end
                if cs.neg.sum ~= 0 && cs.neg.p < 0.05
                    tempmat(find(cs.neg.map == cs.neg.idx),layer) = 1;
                end
            end

            %% plot results in line plot style and mark times of sig correlations
            figure('Renderer','painters','Position',[100 100 1200 800]);
            cmap = cbrewer('seq','YlGnBu',51); cmap(cmap>1) = 1; cmap(cmap<0) = 0;
            
            %prepare x axis labels
            x = 1:size(groupstats.sub.emp.all.t,2);
            maxT = str2num(extractAfter(extractBefore(rsa_tag,'ms'),'step_')); %get max T from roi_tag
            rsawidth = str2double(extractBefore(rsa_tag,'width'));
            rsastep = str2double(extractBefore(extractAfter(rsa_tag,'width_'),'step'));
            timepoints = rsawidth/2:rsastep:maxT;
            
            %plot across-subject average for each layer
            subavg = squeeze(nanmean(groupstats.sub.emp.(field).t,1));
            for layer = 1:layerNum
                if strcmp(dnnType,'AlexNet') || strcmp(dnnType,'AlexNet_avg')
                    subavg_lay = squeeze(subavg(:,layer));
                else
                    subavg_lay = subavg;
                end

                plot(x,subavg_lay,'color',[0.8 0.8 0.8],'LineWidth',3); hold on;
            end
            ylim([min(min(subavg))-0.01 max(max(subavg))+0.01]);
            
            % add legend and mark layers with sig timepoints with *
            legLabels = layerNames;
            %find sig layers;
            [r,c] = find(tempmat); sigLays = unique(c);
            for lay = 1:length(sigLays)
                legLabels{sigLays(lay)} = char([legLabels{sigLays(lay)} ' * ']);
            end
            leg = legend(legLabels,'AutoUpdate','off');
            %plot yline at 0
            yline(0,'black'); hold on;

            %also plot shaded area of sig timepoints
            %first find sig points and layers
%             [time,lay] = find(tempmat);
%             for un = 1:length(unique(lay))
%                 uniques = unique(lay);
%                 siglay = uniques(un);
%                 times = time(find(lay==siglay));
% 
%                 ys = subavg(times,siglay);
% 
%                 a = area(times,ys,'FaceColor',[0.9 0.9 0.9]); a.FaceAlpha = 0.5; 
%                 a.EdgeColor = 'none';
%                 hold on;
%             end
            %set title and axis labels
            set(gca, 'Children', flipud(get(gca, 'Children')) )
            ylabel('subj-avg t-value'); xlabel('time');
            xticks(x); xticklabels(timepoints); xtickangle(90);
            sgtitle([phTag ' | ' field ' trials | ' replace(roi_tag,'_',' ') ' | ' num2str(length(groupstats.subs)) ' subjects | ' num2str(sum(groupstats.chanNum)) ' electrodes '])
            
            %% save stats and figs
            save([groupstatsDir fs 'clusterstats_' ph '_' field '.mat'],'clusterstats')
            saveas(gcf,[groupstatsDir fs phTag '_DNNsim_' field '.fig'])
        end
    end
end  