%% first level stats for Btw RSA, test within phase
% 1D across temporal domain (25 timewins per phase, not temporal
% generalization

function Btw_first_level(subject,rsa_tag,roi_tag,nrep)
    fs = '\';
    dataDir = '..\Data';
    ftpath = 'C:\Program Files\MATLAB\R2020b\toolbox\fieldtrip';
        
    %load folder names of subjects
    load('..\subject_list.mat'); %subNames for easier folder accessing
    
    fileName = subNames{subject}; %set folderName of current subject
    disp(fileName)    
    
    
    %% load Between Item Similarity Matrices + info
    load([dataDir fs fileName fs 'stats' fs 'BtwItem' fs rsa_tag fs roi_tag fs 'BtwItem_encoding.mat']); %encoding
    load([dataDir fs fileName fs 'stats' fs 'BtwItem' fs rsa_tag fs roi_tag fs 'BtwItem_recognition.mat']); %recognition
    
    BtwInfo = load([dataDir fs fileName fs 'stats' fs  'BtwItem' fs rsa_tag fs roi_tag fs 'BtwInfo.mat']); %info
    BtwInfo = BtwInfo.info; clear info;

    numTrials = size(BtwInfo.trialinfo,1);
    
    %% prepare indices for proper indexing (indices of lower matrix elements)
    exampleMat = ones(numTrials) - triu(ones(numTrials),1) - eye(numTrials); %half matrix of ones and subj-specific size
    lower_idx = find(exampleMat); %indices for all lower matrix elements
    
    %% prep indices of all pairwise cells according to memory, 
    remembered_ix = find(BtwInfo.trialinfo(:,2) == 1);
    forgotten_ix = find(BtwInfo.trialinfo(:,2) == 2);
    memMat = exampleMat;
    memMat(remembered_ix,remembered_ix) = 2;
    memMat(forgotten_ix,forgotten_ix) = 3;
    rem_ix = intersect(lower_idx,find(memMat==2));
    forg_ix = intersect(lower_idx,find(memMat==3));
        
    %% compute contrast of btw item sims of remembered vs forgotten items

    for phase = 1:2
        if phase == 1
            data = encoding;
        elseif phase == 2
            data = recognition;
        end

        emp_stats = {};
        emp_stats.t = zeros(1,size(data,1),1); %empty time matrix
        emp_stats.p = zeros(size(emp_stats.t));
        emp_stats.avg.rem = zeros(size(emp_stats.t));
        emp_stats.avg.forg = zeros(size(emp_stats.t));
        perm_stats = {};
        perm_stats.t = zeros(nrep,size(data,1)); %nrep * time 
        perm_stats.p = zeros(size(perm_stats.t));
        perm_stats.avg.rem = zeros(size(perm_stats.t));
        perm_stats.avg.forg = zeros(size(perm_stats.t));

        for timebin = 1:size(data,1) %time1
            vals = squeeze(data(timebin,:,:));
            rem_vals = vals(rem_ix);
            forg_vals = vals(forg_ix);
            [h,p,ci,stat] = ttest2(rem_vals,forg_vals);
            emp_stats.t(timebin) = stat.tstat;
            emp_stats.p(timebin) = p;
            emp_stats.avg.rem(timebin) = nanmean(rem_vals);
            emp_stats.avg.forg(timebin) = nanmean(forg_vals);

            for rep = 1:nrep
                all_vals = [rem_vals; forg_vals];
                shuff_vals = all_vals(randperm(length(all_vals)));
                shuff_rem = shuff_vals(1:length(rem_vals));
                shuff_forg = shuff_vals(length(rem_vals)+1:end);
                [h,p,ci,stat] = ttest2(shuff_rem, shuff_forg);
                perm_stats.t(rep,timebin) = stat.tstat;
                perm_stats.p(rep,timebin) = p;
                perm_stats.avg.rem = nanmean(shuff_rem);
                perm_stats.avg.forg = nanmean(shuff_forg);
            end
        end

        %% save stats to disk
        saveDir = [dataDir fs fileName fs 'stats' fs 'BtwItem' fs rsa_tag fs roi_tag fs 'Btw_within_1D'];
        status = mkdir(saveDir);

        if phase == 1
            fileName1 = ['enc_btw_within_emp.mat'];
            fileName2 = ['enc_btw_within_perm.mat'];
        elseif phase == 2
            fileName1 = ['rec_btw_within_emp.mat'];
            fileName2 = ['rec_btw_within_perm.mat'];
        end

        save([saveDir fs fileName1],'emp_stats')
        save([saveDir fs fileName2],'perm_stats')
    end
end
    

        