%% after computation of aTC (BA20) and LPC (BA39+BA40), extract trial-wise ERS
% for Mixed Model Analysis

%% ANALYSIS SETTINGS
%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

% DEFINE  REGIONS OF INTEREST
ROI1 = [20]; %V1 V2
roi1_tag = ['BA' char(replace(join(string(ROI1)),' ','_'))];
%roi1_tag = ['BA40']
roi1_direction = 'pos';

ROI2 = [39 40]; %[8 44 45 46 47]; %8 dlPFC 44-47 IFG
roi2_tag = ['BA' char(replace(join(string(ROI2)),' ','_'))];
roi2_direction = 'neg';

%regions_of_interest = {ROI1; ROI2};
%regions_of_interest_labels = {ROI1_label; ROI2_label};

% 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'); %POW ROI
[roi2_subs,roi2_channum,roi2_chanix] = FindSubjOfROI_012423(all_subjects,'BA',ROI2,'both','all'); %ERS 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])

%% 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);
[roi2_idx,roi2_enct,roi2_rect] = extract_ersclusteridx_032122(rsa_tag,stats_tag,testType,roi2_tag,roi2_direction);

%plot for safety
t = zeros(numWins); 
t1_map = t;
if strcmp(roi1_direction,'pos')
    t1_map(roi1_idx) = 1;
else
    t1_map(roi1_idx) = -1;
end
t2_map = t;
if strcmp(roi2_direction,'pos')
    t2_map(roi2_idx) = 1;
else
    t2_map(roi2_idx) = -1;
end

%overlap
t3_map = t;
t3_map(intersect(roi1_idx,roi2_idx)) = 1;

% figure; 
subplot(1,3,1); imagesc(t1_map); title([roi1_tag ' ' roi1_direction]); set(gca,'Ydir','normal'); colorbar(); colormap(cmap);
subplot(1,3,2); imagesc(t2_map); title([roi2_tag ' ' roi2_direction]); set(gca,'Ydir','normal'); colorbar(); colormap(cmap);
subplot(1,3,3); imagesc(t3_map); title(['overlap of ' roi1_tag ' and ' roi2_tag]); set(gca,'Ydir','normal'); colorbar(); colormap(cmap);

%%% uncomment to use only timepoints where clusters overlap
roi1_idx = intersect(roi1_idx,roi2_idx);
roi2_idx = intersect(roi1_idx,roi2_idx);

%% load subj-specific data and correlate across trials
subs_rho = zeros(length(subs_in_both),4); %all, rem, forg, rem-forg

%figure('renderer','painters','position',[500 200 1300 800]);
hold on;
ns = numSubplots(length(subs_in_both));
all_trials = {};
rem_trials = {};
forg_trials = {};

for num = 1:length(subs_in_both)
    fileName = subNames{subs_in_both(num)};
    
    %load ERS info
    ERSinfo = load([dataDir fs fileName fs 'stats' fs 'ERS' fs rsa_tag fs roi1_tag fs 'ERSinfo.mat']); ERSinfo = ERSinfo.info; clear info;
    rem_ix = find(ERSinfo.trialinfo(:,2)==1);
    forg_ix = find(ERSinfo.trialinfo(:,2) == 2);

    
    %add info about indoor/outdoor category! (R2.24)
    in_ix = find(ERSinfo.trialinfo(:,5) == 1);
    out_ix = find(ERSinfo.trialinfo(:,5) == 2);
    cat_bin = ERSinfo.trialinfo(:,5);
    
    %add info about memory category as a predictor (R3.12)
    mem_bin = ERSinfo.trialinfo(:,2);
    mem_bin(find(mem_bin==2)) = 0;
    
    %load item-specific stats for roi1
    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)
    %avg over timepoints of sig cluster
    ers_roi1 = squeeze(nanmean(itemrnst1(:,roi1_idx),2));
    
    %load item-specific stats for roi2
    load([dataDir fs fileName fs 'stats' fs 'ERS' fs rsa_tag fs roi2_tag fs stats_tag fs 'itemspecific_reinstatement.mat'])
    %reshape for proper indexing of timepoints (to stim X time X time)
    itemrnst2 = permute(Reinstatement.ITEM,[3,1,2]); %permute dimensions (last (stim) first, then first(time1) and second(time2)
    %avg over timepoints of sig cluster
    ers_roi2 = squeeze(nanmean(itemrnst2(:,roi2_idx),2));
    
    %%%% add trials to trial-overall structure
    all_trials{num,1} = [ers_roi1,ers_roi2,mem_bin,cat_bin];
    rem_trials{num,1} = [ers_roi1(rem_ix),ers_roi2(rem_ix)];
    forg_trials{num,1} = [ers_roi1(forg_ix),ers_roi2(forg_ix)];
    

    %% correlate across all trials
    %all trials
    [rho,~] = corr(ers_roi1,ers_roi2,'rows','complete','type','spearman');
    subs_rho(num,1) = rho;
    %rem trials
    [rho2,~] = corr(ers_roi1(rem_ix),ers_roi2(rem_ix),'rows','complete','type','spearman');
    subs_rho(num,2) = rho2;
    %forg trials
    [rho3,~] = corr(ers_roi1(forg_ix),ers_roi2(forg_ix),'rows','complete','type','spearman');
    subs_rho(num,3) = rho3;
    %memory diff (rem - forg)
    subs_rho(num,4) = rho2-rho3;
end

%% CREATE ONE MODEL WITH MEMORY AND CATEGORY AS BINARY PREDICTOR 
subnums = subs_in_both;
all_trials_mat = cell2mat(all_trials);

%create according subnum
subnum_vec = [];
memory_bin = [];

for k = 1:length(subnums)
    n_rem = length(rem_trials{k});
    n_forg = length(forg_trials{k});
    bin = [repmat(1,1,n_rem) repmat(0,1,n_forg)]';
    
    for trial = 1:length(all_trials{k})
        subnum_vec(end+1,1) = subnums(k);
        %subage_vec(end+1,1) = age(subnums(k));
        %subperf_vec(end+1,1) = recognition(subnums(k));
    end       
    memory_bin(end+1:end+length(bin),1) = bin;
end

data_all = table(subnum_vec, memory_bin, all_trials_mat(:,1), all_trials_mat(:,2), all_trials_mat(:,4),'VariableNames',{'subject','memory',roi1_tag,roi2_tag,'category'});

%do mixed model
lme1 = fitlme(data_all, 'BA20 ~ BA39_40 * memory + category + (1|subject)'); %full model info
anova_lme1 = anova(lme1); % stats and p-values of interest

%rem trials
data_hit = data_all(data_all.memory==1,:);
lme2 = fitlme(data_hit, 'BA20 ~ BA39_40 + category + (1|subject)'); %full model info
anova_lme2 = anova(lme2); % stats and p-values of interest

%forg trials
data_miss = data_all(data_all.memory == 0,:);
lme3 = fitlme(data_miss, 'BA20 ~ BA39_40 + category + (1|subject)'); %full model info
anova_lme3 = anova(lme3); % stats and p-values of interest

%%% save data as csv
file_name = ['ERS_' roi1_tag '-' roi2_tag '_all_trials_SAME_memcat.csv'];
writetable(data_all,['..\STATS\EES_ERS_Connectivity\' file_name]);

%% Correlate across all trials (not within subjects!)
all_trials_mat = cell2mat(all_trials);
plot_correlation_raw(all_trials_mat(:,1), all_trials_mat(:,2),[roi1_tag ' reinstatement'],[roi2_tag ' reinstatement'],'spearman');

rem_trials_mat = cell2mat(rem_trials);
plot_correlation_raw(rem_trials_mat(:,1), rem_trials_mat(:,2),[roi1_tag ' reinstatement'],[roi2_tag ' reinstatement'],'spearman');

forg_trials_mat = cell2mat(forg_trials);
plot_correlation_raw(forg_trials_mat(:,1), forg_trials_mat(:,2),[roi1_tag ' reinstatement'],[roi2_tag ' reinstatement'],'spearman');

