%%% Test whether Between Item Similarities of remembered items differ from
%%% forgotten items, separately for encoding and recognition

% compute unpaired t-test of pairwise similarities in tempGen (enctime X
% enctime)

function btwitem_memoryreinstatement(atlasType, atlasIdx)
    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']);

    %% 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)
    atlasType = atlasType; % 'BA' for Brodmann, or 'HO' for HarvardOxford
    atlasIdx = atlasIdx;

    %set directory names for storing (roi selection)
    if strcmp(atlasType,'BA')
        roi_tag = [atlasType char(replace(join(string(atlasIdx)),' ','_'))];
    end

    %% find all subjects that contribute to current ROI analysis
    [subjects,numChans,channels] = FindSubjOfROI_012423(all_subjects,atlasType,atlasIdx,'both','all');
    disp([num2str(length(subjects)) ' subjects contribute ' num2str(sum(numChans)) ' channels in ' roi_tag])
    %for all further analyses run from this script, only subjects who
    %contribute channels are included!
    
    %% Between Item Similarity

    %computes stimXstim similarity matrices across channels and freqs
    %for each timestep of sliding window of rsa settings
    %separately for encoding and recognition
    %trial order is sorted by memory and category

    for num = 1:length(subjects)
        %check whether BtwItem Similarity of rsa features and for roi
        %already exists!
        file_dir = ['..\Data' fs subNames{subjects(num)} fs 'stats' fs 'BtwItem' fs rsa_tag fs roi_tag];
        BetweenItemSimilarity(subjects(num),channels{num},atlasType,atlasIdx,rsawidth,rsastep,maxT,rsa_tag,roi_tag);
    end
    
    %% First-level Between Item Similarity of Neural Geometries
    disp([' ----- FIRST LEVEL STATS ----- '])
    tic
    nrep = 1000;
    
    parfor num = 1:length(subjects)
        Btw_within_first_level_stats_1D(subjects(num),rsa_tag,roi_tag,nrep)
    end
    disp(['first level stats took ' num2str(toc/60) ' minutes'])
    
    %% second-level-stats
    testType = 'memory';
    dir_tag = 'BtwWithin'; %folder name under which 2nd lvl stats are stored (.. /STATS/ERS_btwitem/ .. )
    disp(['second lvl stats for ' roi_tag])
    Btw_within_second_level_stats_1D(subjects,testType,nrep,rsa_tag,roi_tag)
    %disp([' second level stats took ' num2str(toc/60) ' minutes'])

    %% clusterplots
    %1d plot
    %plot_1d_clusterstats({'BA17_18_19'},{'Occ'})

end

%%% plot results
% plot_1d_clusterstats({'BA17_18_19','BA37','BA20','BA39_40'},{'Occ','pTC','aTC','LPC'})


%% post-hoc test of correlation with age/performance per timestep
% %%% correlate age with t-stat across time
% cmap = cbrewer('div','RdBu',40); cmap = flipud(cmap); cmap(cmap>1) = 1; cmap(cmap<0) = 0;
% subage = age(groupstats.subnum);
% 
% param = recognition(groupstats.subnum);
% param = age(groupstats.subnum)';
% 
% rs = zeros(25,1); ps = zeros(25,1);
% for t1 = 1:25
%     [rho,p] = corr(param,squeeze(groupstats.sub.emp(:,t1)),'rows','complete','type','spearman');
%     rs(t1) = rho; ps(t1) = p;
% end
% figure; 
% subplot(1,3,1); plot(groupstats.emp.t); colorbar(); set(gca,'Ydir','normal'); title('BtwERS'); colormap(cmap); hold on;
% subplot(1,3,2); plot(rs); colorbar(); set(gca,'Ydir','normal'); title('age correlation (rho)'); colormap(cmap); hold on;
% subplot(1,3,3); plot(ps); colorbar(); set(gca,'Ydir','normal'); caxis([0 0.05]); title('correlation p-value'); colormap(cmap);
% 
% %% age-related analysis
% clustidx = find(groupstats.emp_cluster.neg.map == groupstats.emp_cluster.neg.idx);
% temp = zeros(25); temp(clustidx) = 1; %imagesc(temp)
% 
% subt = squeeze(nanmean(groupstats.sub.emp(:,clustidx),2));
% 
% plot_correlation(param,subt,'age','t-stat','spearman')
% plot_correlation_raw(param,subt,'performance','RRS','spearman')
% 
% 
%% post-hoc plot barscatter of avg correlation across subjects
% rem = [];
% forg = [];
% 
% for num = 1:length(subjects)
%     sub = subNames{subjects(num)};
%     load(['F:\iEEG_PFC\ECoG_MemDev\DATA' fs sub fs 'stats' fs 'BtwItem' fs rsa_tag fs roi_tag fs 'Btw_within_1D' fs 'enc_btw_within_emp.mat'])
%     
%     rem(num) = nanmean(emp_stats.avg.rem(8:25));
%     forg(num) = nanmean(emp_stats.avg.forg(8:25));
% end
% 
% %% plot barscatter for remembered vs forgotten of ITEM
% %figure('Renderer','painters','Position',[300 600 300 350]);
% X1 = ones(length(rem),1)*1;
% X2 = ones(length(forg),1)*2;
% %array values per subject
% Y1 = rem;
% Y2 = forg;
% %colors
% rem_c = [0.7412 0.2745 0.4431]; %prev [0 0.5 0.2]
% forg_c = [0.3294 0.4471 0.7686]; %prev [0.7 0.2 0]
% cls = [rem_c; forg_c];
% grey = [0.8 0.8 0.8];
% 
% % make it rain REM VS FORG (ITEM = SAME - DIFF)   
% f = figure('Position', [650 550 420 350]);
% d = {Y1, Y2};   
% %plot rainclouds
% h1 = raincloud_plot(d{1}, 'box_on', 1, 'color', rem_c, 'alpha', 0.5,'box_dodge', 1, 'box_dodge_amount', 0.15, 'dot_dodge_amount', .15, 'box_col_match', 1);
% h2 = raincloud_plot(d{2}, 'box_on', 1, 'color', forg_c, 'alpha', 0.5,'box_dodge', 1, 'box_dodge_amount', 0.45, 'dot_dodge_amount', .45, 'box_col_match', 1);
% %add mean point
% scatter([mean(d{1}) mean(d{2})],[0 0],40,[rem_c; forg_c],'o','filled','MarkerEdgeColor','black','LineWidth',2)
% 
% %adjust raindrops
% h1{2}.SizeData = 30; h1{2}.MarkerEdgeColor=grey; h1{2}.MarkerFaceColor = grey;
% h2{2}.SizeData = 30; h2{2}.MarkerEdgeColor=grey; h2{2}.MarkerFaceColor = grey; 
% %make scatter align on one line and draw rep connecting lines
% h1{2}.YData = repmat(h1{5}.YData(1), length(h1{2}.YData), 1);
% h2{2}.YData = repmat(h2{5}.YData(1), length(h2{2}.YData), 1);
% 
% %remove line of density dist
% h1{1}.EdgeColor = 'none';
% h2{1}.EdgeColor = 'none';
% 
% %add rep mes lines
% lines = plot([h1{2}.XData; h2{2}.XData],[h1{2}.YData; h2{2}.YData],'Color',grey);
% uistack(lines,'bottom'); 
% %add legend 
% legend([h1{1} h2{1}], {'Remembered', 'Forgotten'})
% legend('AutoUpdate','off');
% %add labels and title
% title('Rem vs Forg');
% yticklabels({});
% xlabel('EES (spearman rho)');
% %view mode
% set(gca,'XLim', [-0.15 0.35], 'YLim', [-10 11]);
% if strcmp(roi_tag,'BA20')
%     set(gca,'XLim', [-0.2 0.15], 'YLim', [-10 18]);
% end
% camroll(90);
% set(gcf,'color','white')
% xline(0,'--');
% 
% %post-hoc ttest of rem vs forg
% [~,p,~,stats]= ttest(Y1,Y2);
% disp(['rem vs forg t = ' num2str(stats.tstat) ' , p = ' num2str(p) ' , df = ' num2str(stats.df)]);
% 
% title(['t = ' num2str(stats.tstat) ' , p = ' num2str(p) ' , df = ' num2str(stats.df)]);
% 
% %test remembered and forgotten against 0!
% [h,p,ci,stats] = ttest(Y1);
% disp(['rem vs 0: t = ' num2str(stats.tstat) ' , p = ' num2str(p) ' , df = ' num2str(stats.df)]);
% [h,p,ci,stats] = ttest(Y2);
% disp(['forg vs 0: t = ' num2str(stats.tstat) ' , p = ' num2str(p) ' , df = ' num2str(stats.df)]);
% disp(' ');
    
