%% Settings
clearvars;
%save to disk
save_to_disk = 0;

% select subjects
all_subjects = [1:22 24:28 30:37 39:49]; %%%  exclude 23 and 29 for negative d' and 38 for different stimset and 2nd time implant
load('..\subject_list.mat'); %subNames for easier folder accessing

% ERS ROI
ROI1 = [20]; %aTC
roi1_tag = ['BA' char(replace(join(string(ROI1)),' ','_'))];
roi1_direction = 'pos';

% EES/RRS ROI
BtwType = 'encoding';
ROI2 = [20]; %aTC
roi2_tag = ['BA' char(replace(join(string(ROI2)),' ','_'))];

% set basedir
fs = '\';
dataDir = '..\DATA';

% RSA feature settings
testType = 'memory';
rsawidth = 300;
rsastep = 50;
maxT = 1500;
numWins = 25;
rsa_tag = [num2str(rsawidth) 'width_' num2str(rsastep) 'step_' num2str(maxT) 'ms'];

% set settings of ERS stats to analyze here
control = 'memcat'; %Item-specific as SAME - avg( all / category / memory / memcat)
diff_type = 'row';  %consider DIFF of same: row, col, cross
nrep = 1000;        
%set directory names for storing (1st level stats settings)
stats_tag = ['item-specific' fs 'diff_' control '_' diff_type '_' num2str(nrep) 'rep'];

% Colormap
cmap = flipud(cbrewer2('div','RdBu',51)); cmap(cmap>1) = 1; cmap(cmap<0) = 0;
    
%% Fetch subjects and plot channels of both ROIs

%find all subjects with respective channels and visualize them
% find all subjects and channels who contribute to both specified ROIs
[roi1_subs,roi1_channum,roi1_chanix] = FindSubjOfROI_012423(all_subjects,'BA',ROI1,'both','all'); %ERS ROI
[roi2_subs,roi2_channum,roi2_chanix] = FindSubjOfROI_012423(all_subjects,'BA',ROI2,'both','all'); %EES ROI

%find their overlap
subs_in_both = intersect(roi1_subs, roi2_subs);
roi1_chans = roi1_chanix(find(ismember(roi1_subs, subs_in_both)));
roi2_chans = roi2_chanix(find(ismember(roi2_subs, subs_in_both)));
disp([num2str(length(subs_in_both)) ' have ' num2str(sum(cellfun(@length, roi1_chans))) ' in ' roi1_tag ' and ' num2str(sum(cellfun(@length, roi2_chans))) ' in ' roi2_tag])

%% avg ers in cluster, avg ees in same timepoints (during ers)
% get the times of significant clusters from ers analysis
[roi1_idx,roi1_enct,roi1_rect] = extract_ersclusteridx_032122(rsa_tag,stats_tag,testType,roi1_tag,roi1_direction);

ers_timepoints = roi1_idx;  %idx of ERS cluster
ers_encoding_timepoints = roi1_enct; %enc wins of ERS cluster (not from EES analysis!)
%ees_timepoints = [15 16]; %aTC EES effect timepoints
ers_recognition_timepoints = roi1_rect; %rec wins of ERS cluster (not from RRS analysis!)

%% load subj-specific data and correlate across trials

all_trials = {};
rem_trials = {};
forg_trials = {};

for num = 1:length(subs_in_both)
    fileName = subNames{subs_in_both(num)};
    
    %% load item-specific ERS 
    ERSinfo = load([dataDir fs fileName fs 'stats' fs 'ERS' fs rsa_tag fs roi1_tag fs 'ERSinfo.mat']); ERSinfo = ERSinfo.info; clear info;
    %fetch indices of trials by memory + inout category
    rem_ix = find(ERSinfo.trialinfo(:,2)==1);
    forg_ix = find(ERSinfo.trialinfo(:,2) == 2);
    indoor_ix = find(ERSinfo.trialinfo(:,5)==1);
    outdoor_ix = find(ERSinfo.trialinfo(:,5)==2);
    rem_in_ix = intersect(rem_ix, indoor_ix);
    rem_out_ix = intersect(rem_ix, outdoor_ix);
    forg_in_ix = intersect(forg_ix, indoor_ix);
    forg_out_ix = intersect(forg_ix, outdoor_ix);
    
    %%% subject and memory vector
    subnum_vec = repmat(subs_in_both(num),length(ERSinfo.trialinfo),1);
    memory_vec = [repmat(1,length(rem_ix),1); repmat(0,length(forg_ix),1)];
    category_vec = [repmat(1,length(indoor_ix),1); repmat(0,length(outdoor_ix),1)];
    
    %% load item-specific ERS
    load([dataDir fs fileName fs 'stats' fs 'ERS' fs rsa_tag fs roi1_tag fs stats_tag fs 'itemspecific_reinstatement.mat'])
    %reshape for proper indexing of timepoints (to stim X time X time)
    itemrnst1 = permute(Reinstatement.ITEM,[3,1,2]); %permute dimensions (last (stim) first, then first(time1) and second(time2)
    samernst1 = permute(Reinstatement.SAME,[3,1,2]);
    diffrnst1 = permute(Reinstatement.DIFF,[3,1,2]);
    
    %avg over timepoints of sig cluster
    item_ers_roi1 = squeeze(nanmean(itemrnst1(:,roi1_idx),2));
    same_ers_roi1 = squeeze(nanmean(samernst1(:,roi1_idx),2));
    diff_ers_roi1 = squeeze(nanmean(diffrnst1(:,roi1_idx),2));

    %% load item-specific BtwSim
    load([dataDir fs fileName fs 'stats' fs 'BtwItem' fs rsa_tag fs roi2_tag fs 'BtwItem_' BtwType '.mat'])

    if strcmp(BtwType,'encoding')
        data = encoding;
        timepoints = 10:12; %10:12;%ers_encoding_timepoints; %should use EES / RRS cluster timepoints, not whole ERS cluster!
    elseif strcmp(BtwType,'recognition')
        data = recognition;
        timepoints = ers_recognition_timepoints;
    end
    
    %avg over timepoints of sig cluster
    btw_roi2 = squeeze(nanmean(data(timepoints,:,:),1));
    
    %now compute item-wise Between-item similarity -> avg of all row-wise
    %elements
    mirrored = btw_roi2 + btw_roi2'; %flip because ees_roi2 is only half matrix
    
    item_btw_roi2 = zeros(size(mirrored,1),1);
    for scene = 1:size(mirrored,1)            %loop over items and fetch avg sim to all other scenes (row-wise)
        scene_row = mirrored(scene,:);        %complete row
        
        %get avg sim to all others from same memory category
        if ismember(scene,rem_ix)
            scene_row_mem = scene_row(rem_ix); %reduce to all rem trials
            scene_row_mem(scene) = [];         %drop diag entry (0)
        elseif ismember(scene,forg_ix)
            scene_row_mem = scene_row(forg_ix); %reduce to all forg trials
            scene_row_mem(scene - length(rem_ix)) = []; %drop diag entry (0) - for forg trials adapt idx!
        end
        
                
        %avg btw sims and store to array
        item_btw_roi2(scene,1) = nanmean(scene_row_mem);       %item-specific
    end
    
    %%%% add trials to trial-overall structure (subnum vec, memory vec, ers
    %%%% vals, ees vals
    all_trials{num,1} = [repmat(subs_in_both(num),length(item_ers_roi1),1),[repmat(1,length(rem_ix),1)' repmat(0,length(forg_ix),1)']',item_ers_roi1,item_btw_roi2,];
    rem_trials{num,1} = [repmat(subs_in_both(num),length(rem_ix),1),repmat(1,length(rem_ix),1),item_ers_roi1(rem_ix),item_btw_roi2(rem_ix)];
    forg_trials{num,1} = [repmat(subs_in_both(num),length(forg_ix),1),repmat(0,length(forg_ix),1),item_ers_roi1(forg_ix),item_btw_roi2(forg_ix)];

    all_trials{num,1} = [subnum_vec, memory_vec, category_vec, item_ers_roi1, same_ers_roi1, diff_ers_roi1, item_btw_roi2];
end
    
%%
all = cell2mat(all_trials);
rem = cell2mat(rem_trials);
forg = cell2mat(forg_trials);

data_all = table(all(:,1), all(:,2), all(:,3), all(:,4), all(:,5), all(:,6),all(:,7),'VariableNames',{'subject','memory','inout',[roi1_tag '_ITEM'],[roi1_tag '_SAME'],[roi1_tag '_DIFF'],[roi2_tag '_EES']});
%data_rem = table(rem(:,1), rem(:,2), rem(:,3), rem(:,4),'VariableNames',{'subject','memory','inout','ERS','EES'});
%data_forg = table(forg(:,1), forg(:,2), forg(:,3), forg(:,4),'VariableNames',{'subject','memory','inout','ERS','EES'});

data_rem = data_all(data_all.memory == 1,:);
data_forg = data_all(data_all.memory==0,:);

file_name = ['ERS_' roi1_tag '-EES_' roi2_tag 'eeswin_all_trials_memcat.csv'];
writetable(data_all,['..\STATS\EES_ERS_Connectivity\' file_name])

plot_correlation_raw(all(:,5),all(:,7),'ERS','EES','spearman')
plot_correlation_raw(rem(:,3),rem(:,4),'ERS','EES','spearman')
plot_correlation_raw(forg(:,3),forg(:,4),'ERS','EES','spearman')

%% ONE MODEL TO RULE THEM(TRIALS) ALL

%% ATC
%do mixed model across all trials, with memory as binary predictor
lme1 = fitlme(data_all, 'BA20_ITEM ~ BA20_EES + inout + memory + (1|subject)'); %full model info
anova_lme1 = anova(lme1); % stats and p-values of interest

%do mixed model across all trials, with memory as binary predictor
lme1 = fitlme(data_all, 'BA39_40_ITEM ~ BA39_40_EES + inout + memory + (1|subject)'); %full model info
anova_lme1 = anova(lme1); % stats and p-values of interest

%%%% separately for rem and forg
%rem
lme2 = fitlme(data_rem, 'BA20_ITEM ~ BA20_EES + inout + (1|subject)'); %full model info
anova_lme2 = anova(lme2); % stats and p-values of interest

%forg
lme3 = fitlme(data_forg, 'BA20_ITEM ~ BA20_EES + inout + (1|subject)'); %full model info
anova_lme3 = anova(lme3); % stats and p-values of interest

%% LPC
%do mixed model across all trials, with memory as binary predictor
lme1 = fitlme(data_all, 'BA39_40_ITEM ~ BA39_40_EES + inout + memory + (1|subject)'); %full model info
anova_lme1 = anova(lme1); % stats and p-values of interest

%% separately for rem and forg
%rem
lme2 = fitlme(data_rem, 'BA39_40_ITEM ~ BA39_40_EES + inout + (1|subject)'); %full model info
anova_lme2 = anova(lme2); % stats and p-values of interest

%forg
lme3 = fitlme(data_forg, 'BA39_40_ITEM ~ BA39_40_EES + inout + (1|subject)'); %full model info
anova_lme3 = anova(lme3); % stats and p-values of interest

