%% 02 ENCODNG-RETRIEVAL SIMILARITY ANALYSIS PIPELINE (030722)
% 02_Encoding-Retrieval Similarity (ERS)
    %set parameters of RSA features
    %get all subjects who contribute to specified ROI
    %compute ERS for each subject (EncodingRecognitionSimilarity_030722)
    %compute first-level analysis (ERS_first_level_stats_030722)
    %compute second-lvl analysis (ERS_second_level_stats_030822)
    %plot clusterstats and post-hoc condition differences (plot_clusterstats_030822, plot_cluster_barscatter, plot_cluster_raincloud, posthoc_ersitem_avg, posthoc_ers_param)

% to reproduce ERS findings, call the function for each region of interest
% atlasType: 'BA' indicates electrode selection based on Brodmann Areas
% atlasIdx: [17 18 19],[37],[20],[39 40] %BA codes for Occ, pTC, aTC, LPC
% ERS_ANALYSIS_PIPELINE('BA',[17 18 19])
    
% Contact: Elias Rau, elias.rau@rub.de
    
function ERS_memoryreinstatement(atlasType,atlasIdx)
    tic
    fs = '\';

    %% ANALYSIS parameters
    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']); %load subj-specific information

    %% RSA Sliding window settings
    rsawidth = 300; %window size
    rsastep = 50;   %window slide step
    maxT = 1500;    %max time
    %set directory names for storing (ers settings)
    rsa_tag = [num2str(rsawidth) 'width_' num2str(rsastep) 'step_' num2str(maxT) 'ms'];

    %% ROI (Define ROI for which to compute ERS analysis)
    roi_tag = [atlasType char(replace(join(string(atlasIdx)),' ','_'))];

    %% find all subjects that contribute channels to current ROI 
    [subjects,numChans,channels] = FindSubjOfROI_012423(all_subjects,atlasType,atlasIdx,'both','all');
    disp([num2str(length(subjects)) ' subjects contribute ' num2str(sum(numChans)) ' channels in ' roi_tag])

    %% Compute ERS
    disp('COMPUTE ENCODING RECOGNITION SIMILARITY MATRICES')
    parfor num = 1:length(subjects)
        EncodingRecognitionSimilarity(subjects(num),channels{num},atlasType,atlasIdx,rsawidth,rsastep,maxT,rsa_tag,roi_tag)
    end

    %% First level stats
    %executes first-level stats for specified ROIs and each subject

    %first-level stats 
    control = 'memcat'; %Item-specific as SAME - avg( all / category / memory / memcat)
    diff_type = 'row';  %consider DIFF of same: row, col, cross
    nrep = 1000;        %n permutation iterations

    %set directory names for storing (1st level stats settings)
    stats_tag = ['item-specific' fs 'diff_' control '_' diff_type '_' num2str(nrep) 'rep'];
    disp(['RUN FIRST LEVEL STATS'])
    disp(stats_tag)

    parfor num = 1:length(subjects)
        tic
        ERS_first_level_stats(subjects(num),control,diff_type,nrep,rsa_tag,roi_tag,stats_tag)
    end

    %% second-level-stats
    testType = 'memory';
    dir_tag = 'ERS'; %folder name under which 2nd lvl stats are stored (.. /STATS/ERS/ .. )
    disp(['second lvl stats for ' roi_tag])
    ERS_second_level_stats(subjects,testType,nrep,rsa_tag,roi_tag,stats_tag,dir_tag)
    %disp([' second level stats took ' num2str(toc/60) ' minutes'])

    %% clusterplots - plots timextime cluster statistics across subjects
    disp(['plot clusterstats for ' roi_tag])
    plot_clusterstats(rsawidth,rsastep,maxT,testType,nrep,rsa_tag,roi_tag,stats_tag,dir_tag)
    disp(['ERS ANALYSIS FOR ' roi_tag ' took ' num2str(toc/60) ' minutes'])
    
    %% post-hoc cluster - plot subj-avg correlation scores for SAME (rem forg) and DIFF (rem forg)
    if strcmp(roi_tag,'BA20')
        direction = 'pos'; %aTC Reinstatement 
    elseif strcmp(roi_tag,'BA39_40')
        direction = 'neg'; %LPC Transformation
    else
        direction = 'pos'; %look for reinstatement effect in other ROIs
    end
    
    posthoc_cluster_barscatter(subjects,rsa_tag,roi_tag,stats_tag,testType,direction)
    posthoc_cluster_raincloud(subjects,rsa_tag,roi_tag,stats_tag,testType,direction)
    
    %% plot avg ERSitem scores across subjects for all timeXtime points
    posthoc_ersitem_avg(subjects,atlasType,atlasIdx,rsa_tag,roi_tag,stats_tag) 
    
    %% plot post-hoc cluster t-stat by param (age, performance)
    posthoc_ers_param(roi_tag,rsa_tag,stats_tag,'age','spearman',direction)
    posthoc_ers_param(roi_tag,rsa_tag,stats_tag,'performance','spearman',direction)
end

% % % %%% post-hoc median split of age for BA40 transformation cluster
% %get subjects from BA40 roi
% roi = 'BA39_40';
% dir = 'neg';
% 
% all_subjects = [1:22 24:28 30:37 39:49];
% [subjects,numChans] = FindSubjOfROI_030822(all_subjects,'BA',[39 40]);
% load(['F:\iEEG_PFC\ECoG_MemDev\subject_list.mat']);
% 
% 
% %%%% MD SPLIT AGE
% %get their age and perform a median split
% age_subjects = age(subjects);
% md = median(age_subjects);
% disp(['cut off age of median split is: ' num2str(md)])
% old_subjects = subjects(age_subjects > md);
% young_subjects = subjects(age_subjects <= md);
% 
% %perform post-hoc barscatter with young and old
% %old
% posthoc_cluster_barscatter(old_subjects,'300width_50step_1500ms',roi,'item-specific\diff_memcat_row_1000rep','memory',dir)
% %young
% posthoc_cluster_barscatter(young_subjects,'300width_50step_1500ms',roi,'item-specific\diff_memcat_row_1000rep','memory',dir)
% 
% % %est cluster-avg ttest against 0 for young and old 
% % %load cluster-avg subject-specific t-value
% load(['F:\iEEG_PFC\ECoG_MemDev\STATS\ERS\300width_50step_1500ms\item-specific\diff_memcat_row_1000rep\memory\BA39_40_stats.mat']);
% clustidx = find(groupstats.emp_cluster.neg.map == groupstats.emp_cluster.neg.idx);
% cluster_t = squeeze(nanmean(groupstats.sub.emp(:,clustidx),2));
% %old
% [h,p,ci,stats]= ttest(cluster_t(find(age_subjects > md)))
% disp(['old subs, cluster t: t = ' num2str(stats.tstat) ' , p = ' num2str(p) ' , df = ' num2str(stats.df)])
% %young
% [h,p,ci,stats]= ttest(cluster_t(find(age_subjects <= md)))
% disp(['young subs, cluster t: t = ' num2str(stats.tstat) ' , p = ' num2str(p) ' , df = ' num2str(stats.df)])
% 
% % %%%% MD SPLIT PERFORMANCE
% %get performance and perform a median split
% perf_subjects = recognition(subjects);
% md = median(perf_subjects);
% disp(['cut off performance of median split is: ' num2str(md)])
% good_subjects = subjects(perf_subjects > md);
% bad_subjects = subjects(perf_subjects <= md);
% 
% %perform post-hoc barscatter with young and old
% %good
% posthoc_cluster_barscatter(good_subjects,'300width_50step_1500ms',roi,'item-specific\diff_memcat_row_1000rep','memory',dir)
% posthoc_cluster_raincloud(good_subjects,'300width_50step_1500ms',roi,'item-specific\diff_memcat_row_1000rep','memory',dir)
% %bad
% posthoc_cluster_barscatter(bad_subjects,'300width_50step_1500ms',roi,'item-specific\diff_memcat_row_1000rep','memory',dir)
% posthoc_cluster_raincloud(bad_subjects,'300width_50step_1500ms',roi,'item-specific\diff_memcat_row_1000rep','memory',dir)
% 
% %%% show LPC ERS effect (rem-forg) separately for low and high performers
% bad_diff = [0.0048,-0.0723,-0.0697,-0.0124,-0.0196,-0.0270,0.0010,-0.0193,-0.0803,-0.0395,-0.1070,-0.0228,-0.0090,-0.0435,-0.0101,-0.0281,-0.0377,-0.0627,-0.0294]
% good_diff = [-0.0109,0.0342,0.0324,-0.0294,0.0392,0.0403,-0.0378,0.0573,-0.0054,0.0962,-0.0189,-0.0743,-0.0612,0.011,-0.06180,0.0025,0.0566,-0.0469,-0.0049];
% data = {bad_diff, good_diff};
% figure('Position',[700 550 350 450]);
% cols = cbrewer('qual','Paired',12);
% cols(cols>1) = 1; cols(cols<0) = 0;
% bad_perf_c = cols(9,:);
% good_perf_c = cols(10,:);
% h1 = raincloud_plot(data{1}, 'box_on', 1, 'color', bad_perf_c, 'alpha', 0.8,'box_dodge',1,'box_dodge_amount',0.2,'dot_dodge_amount', 0.2, 'box_col_match', 1);
% h2 = raincloud_plot(data{2}, 'box_on', 1, 'color', good_perf_c, 'alpha', 0.8,'box_dodge',1,'box_dodge_amount',0.45,'dot_dodge_amount', 0.45, 'box_col_match', 1);
% % 
% %adjust raindrops
% h1{2}.SizeData = 25; h1{2}.MarkerEdgeColor=[0.8 0.8 0.8]; h1{2}.MarkerFaceColor = [0.8 0.8 0.8];
% h2{2}.SizeData = 25; h2{2}.MarkerEdgeColor=[0.8 0.8 0.8]; h2{2}.MarkerFaceColor = [0.8 0.8 0.8];
% 
% %remove line of density dist
% h1{1}.EdgeColor = 'none';
% h2{1}.EdgeColor = 'none';
% 
% %add legend 
% legend([h1{1} h2{1}], {'low performance','high performance'})
% legend('AutoUpdate','off');
% %add labels and title
% title('LPC Transformation');
% yticklabels({});
% xlabel('ERS (rem-forg)');
% %view mode
% set(gca,'XLim', [-0.2 0.2], 'YLim', [-10 15]);
% camroll(90);
% set(gcf,'color','white')
% xline(0,'--')
% 
% 
% %% make a final plot of low and high performers. For each show (rem-forg) diff for SAME and DIFF items
% % low will show sig. deviation from 0 for same and diff, high will show
% % distributions around zero
% good_same = [-0.0115,0.0288,-0.056,0.0082,0.0405,0.0272,0.0238,0.0280,-0.02390,0.0068,0.0056,-0.0023,-0.0532,0.0129,-0.0254,0.0058,0.0416,-0.0491,-0.0147]
% good_diff = [-0.0006,-0.0054,-0.0884,0.0376,0.0013,-0.0131,0.055,-0.0293,-0.0095,-0.0723,0.0246,0.0719,0.008,0.0019,0.0364,0.0033,-0.015,-0.0022,-0.0098]
% bad_same = [0.0108,-0.0650,-0.0396,0.0128,-0.0128,-0.0278,-0.0017,-0.0232,-0.0634,-0.0268,-0.0135,-0.0022,-0.0063,-0.0148,-0.0053,-0.0270,-0.0535,-0.0536,-0.026]
% bad_diff = [0.0060,0.0073,0.0301,0.0252,0.0067,-0.0008,-0.0027,-0.0039,0.01690,0.0126,0.0936,0.0206,0.0027,0.0072,0.0048,0.001,-0.0158,0.0091,0.0033]
% 
% data = {bad_same, bad_diff, good_same, good_diff};
% 
% figure('Position',[700 550 350 450]);
% 
% cols = cbrewer('qual','Paired',12);
% cols(cols>1) = 1; cols(cols<0) = 0;
% same_c = cols(2,:);
% diff_c = cols(8,:);
% 
% h1 = raincloud_plot(data{1}, 'box_on', 1, 'color', same_c, 'alpha', 0.8,'box_dodge',1,'box_dodge_amount',0.2,'dot_dodge_amount', 0.2, 'box_col_match', 1);
% h2 = raincloud_plot(data{2}, 'box_on', 1, 'color', diff_c, 'alpha', 0.8,'box_dodge',1,'box_dodge_amount',0.45,'dot_dodge_amount', 0.45, 'box_col_match', 1);
% h3 = raincloud_plot(data{3}, 'box_on', 1, 'color', same_c, 'alpha', 0.8,'box_dodge',1,'box_dodge_amount',0.8,'dot_dodge_amount', 0.8, 'box_col_match', 1);
% h4 = raincloud_plot(data{4}, 'box_on', 1, 'color', diff_c, 'alpha', 0.8,'box_dodge',1,'box_dodge_amount',1.05,'dot_dodge_amount', 1.05, 'box_col_match', 1);
% 
% %adjust raindrops
% h1{2}.SizeData = 25; h1{2}.MarkerEdgeColor=[0.8 0.8 0.8]; h1{2}.MarkerFaceColor = [0.8 0.8 0.8];
% h2{2}.SizeData = 25; h2{2}.MarkerEdgeColor=[0.8 0.8 0.8]; h2{2}.MarkerFaceColor = [0.8 0.8 0.8];
% h3{2}.SizeData = 25; h3{2}.MarkerEdgeColor=[0.8 0.8 0.8]; h3{2}.MarkerFaceColor = [0.8 0.8 0.8];
% h4{2}.SizeData = 25; h4{2}.MarkerEdgeColor=[0.8 0.8 0.8]; h4{2}.MarkerFaceColor = [0.8 0.8 0.8];
% 
% %remove line of density dist
% h1{1}.EdgeColor = 'none';
% h2{1}.EdgeColor = 'none';
% h3{1}.EdgeColor = 'none';
% h4{1}.EdgeColor = 'none';
% 
% %add legend 
% legend([h1{1} h2{1} h3{1} h4{1}], {'SAME','DIFF','SAME','DIFF'})
% legend('AutoUpdate','off');
% %add labels and title
% title('Low vs High Performers');
% yticklabels({});
% xlabel('ERS (rem-forg)');
% %view mode
% set(gca,'XLim', [-0.2 0.2], 'YLim', [-20 30]);
% camroll(90);
% set(gcf,'color','white')
% xline(0,'--')
    

% % %est cluster-avg ttest against 0 for young and old 
% % %load cluster-avg subject-specific t-value
% load(['F:\iEEG_PFC\ECoG_MemDev\STATS\ERS\300width_50step_1500ms\item-specific\diff_memcat_row_1000rep\memory\' roi '_stats.mat']);
% clustidx = find(groupstats.emp_cluster.neg.map == groupstats.emp_cluster.neg.idx);
% cluster_t = squeeze(nanmean(groupstats.sub.emp(:,clustidx),2));
% %good
% [h,p,ci,stats]= ttest(cluster_t(find(perf_subjects > md)))
% disp(['old subs, cluster t: t = ' num2str(stats.tstat) ' , p = ' num2str(p) ' , df = ' num2str(stats.df)])
% %bad
% [h,p,ci,stats]= ttest(cluster_t(find(perf_subjects <= md)))
% disp(['young subs, cluster t: t = ' num2str(stats.tstat) ' , p = ' num2str(p) ' , df = ' num2str(stats.df)])