%% Mutual information analysis of visual ROIs 
% Script used to select ROIs in the visual cortical network used for the 
% MEG information transmission analysis

% We first compute stimulus and choice information for all visual ROIs of
% the Glasser et al. atlas, then we select ROIs based on high visual
% stimulus information and on spatial location
% We plot temporal profiles of stimulus information (used for Fig.3C) and
% compare information quantities in the visual network (Fig. S6A)

clear all; clc;
params.randSeed = 1; % for reproducibility

mainPath = pwd;
resultsPath = [mainPath '\Results'];
if ~exist(resultsPath,'dir')
    mkdir(resultsPath)
end

scriptsPath = [mainPath '\Scripts'];
load([mainPath,'\data\rois_info.mat'])
load([mainPath, '\data\times'])

dataPath = [mainPath,'data']; % path of neural data
behavPath = [mainPath,'data']; % path of behavioral data
addpath(genpath(scriptsPath))

% info params
opts.verbose = false;
opts.method = "dr";
opts.bias = 'naive';
opts.btsp = 0;
opts.bin_method = 'eqpop';
opts.n_bins = 2;
%opts.Sbins = 2;
opts.doShuff = 1;
opts.nShuff = 1;
opts.binning_pool = 'trials';
opts.doNiliNull = 0;

% General settings
nSubj = 15;
nWorkers = 1;
hemis_labels = {'left','right','sum','diff'};
info_type = {'S','C'};
correrrLabels = {'all','corr','err'};
opts.clusterParams = [0.05,0.9];
% load([dataPath, '\roinames_pair']) % ROIs in neural data are ordered as left/right/left/right/...

sel_subgroups = [1,2,3,4]; % cortical groups corresponding to primary visual, dorsolateral visual stream, ventral visual stream and MT+ complex, respectively
selROIs = ismember(subgroups_all,sel_subgroups);
selROIs = find(selROIs);
selROIs_labels = ordered_rois(selROIs);
nROIs = numel(selROIs);

%% Compute stim and choice info

if (nWorkers > 1)
    if (isempty(gcp('nocreate')))
        parpool(nWorkers);
    else
        poolTemp = gcp;
        if(poolTemp.NumWorkers < nWorkers)
            delete(gcp)
            parpool(nWorkers);
        end
    end
end

% pre-allocate info maps 
for ceLab = correrrLabels
    for iIdx = 1:numel(info_type)
        iLab = info_type{iIdx};
        for hIdx = 1:numel(hemis_labels)
            hLab = hemis_labels{hIdx};
            MI_maps.(ceLab{1}).(iLab).(hLab) = cell(1,nSubj);
            if opts.doShuff
                MIsh_maps.(ceLab{1}).(iLab).(hLab) = cell(1,nSubj);
            end
        end
    end
end

% Compute stimulus and choice info
for subIdx = 1:nSubj
    behav = load([behavPath, '\S', num2str(subIdx)]);
    neural = load([dataPath, '\S', num2str(subIdx)]);
    
    opts.timeSteps = size(neural.data_subj{1},2);
    % Initialize info structures
    for ceLab = correrrLabels
        for iIdx = 1:numel(info_type)
            iLab = info_type{iIdx};
            for hIdx = 1:numel(hemis_labels)
                hLab = hemis_labels{hIdx};
                MI_maps.(ceLab{1}).(iLab).(hLab){subIdx} = zeros(4,nROIs,opts.timeSteps);
                if opts.doShuff
                    MIsh_maps.(ceLab{1}).(iLab).(hLab){subIdx} = zeros(4,nROIs,opts.timeSteps,opts.nShuff);
                end
            end
        end
    
        for sessIdx = 1:4 
            stim = double(mean(behav.v_cons{sessIdx},2)>0.5)';
            choice = (behav.v_choice{sessIdx}+1)/2;
            %stim = eqpop(stim,opts.Sbins)';

            sel_trials = [];
            if strcmp(ceLab{1},'all')
                sel_trials = logical(ones(1,numel(stim)));
            elseif strcmp(ceLab{1},'err')
                sel_trials = stim~=choice;
            elseif strcmp(ceLab{1},'corr')
                % Random subsamp of correct trials
                nErr = sum(stim~=choice);
                nTrials = numel(stim);
                sel_trials = zeros(1,nTrials);
                nCorr = nTrials - nErr;
                
                correct_idxs = find(stim==choice);
                rng(params.randSeed)
                rIdx = randperm(nCorr);
                rIdx = rIdx(1:nErr);
                sel_trials(correct_idxs(rIdx)) = 1;
                sel_trials = logical(sel_trials);
            end

            opts.nTrials = numel(stim);
            opts.timeSteps = size(neural.data_subj{sessIdx},2);

            % Coherent shuffling across ROIs and time points
            shIdxs = zeros(opts.nTrials,opts.nShuff);
            for shIdx = 1:opts.nShuff
                shIdxs(:,shIdx) = randperm(opts.nTrials);
            end

            % Binning neural signals
            for hIdx = 1:numel(hemis_labels)
                hLab = hemis_labels{hIdx};
                bNeural{sessIdx}.(hLab) = zeros(sum(sel_trials),nROIs,opts.timeSteps);
            end
            for roiIdx = 1:nROIs
                tmpROI = selROIs(roiIdx); % ROI idx according to subgroups ordering
                roiLab = ordered_rois(tmpROI); % ROI label
                roi = find(strcmp(neural.clusters,['/Pair/',roiLab{1},'_RH'])); % search idx of right hemisphere as reference
                for t = 1:opts.timeSteps
                    
                    Rleft = eqpop(neural.data_subj{sessIdx}(sel_trials,t,(roi-1)),opts.n_bins);
                    bNeural{sessIdx}.left(:,t,roiIdx) = equi_binning2(Rleft, opts.n_bins, 0);
                    Rright = eqpop(neural.data_subj{sessIdx}(sel_trials,t,(roi)),opts.n_bins);
                    bNeural{sessIdx}.right(:,t,roiIdx) = equi_binning2(Rright, opts.n_bins, 0);
                    
                end
            end
            
            % Compute mutual information

            for hIdx = 1:numel(hemis_labels)
                hLab = hemis_labels{hIdx};
                disp(['Computing subj ',num2str(subIdx),' sess ', num2str(sessIdx),' ', hLab,' hemispheres - ',ceLab{1}, ' trials'])

                for featLab = fieldnames(bNeural{sessIdx})'
                    for roiIdx = 1:nROIs
                        tmpROI = selROIs(roiIdx); % ROI idx according to subgroups ordering
                        roiLab = ordered_rois(tmpROI); % ROI label

                        tmpbNeural = [];
                        tmpbNeural = squeeze(bNeural{sessIdx}.(hLab)(:,:,roiIdx));
                        tmpbNeural = squeeze(bNeural{sessIdx}.(hLab)(:,:,roiIdx));

    %                     disp(['Subj ',num2str(subIdx),' sess ', num2str(sessIdx), ' roi ', roiLab{1}]);
                        for t = 1:opts.timeSteps
                            % Compute info theoretic quantities

                            tmp_neuSig = squeeze(tmpbNeural(:,t))';
                            
                            % stimulus info
                            [M_rs, nt] = buildr(stim,tmp_neuSig);
                            opts.nt = nt;
                            tmpInfo = information(M_rs, opts, 'I');
                            MI_maps.S.(hLab){subIdx}(sessIdx,roiIdx,t) = tmpInfo;
                            
                            % choice info
                            [M_rc, nt] = buildr(choice,tmp_neuSig);
                            opts.nt = nt;
                            tmpInfo = information(M_rc, opts, 'I');
                            MI_maps.C.(hLab){subIdx}(sessIdx,roiIdx,t) = tmpInfo;
                            
                        end
                    end
                end
            end
        end % end loop over sessions
    end % end loop over corr/err trials
end
        
save([resultsPath,'\allSubj_visual_MI_S_C_hemisCombined_',date])

%% Average results across subjects and sessions
% Initialize grand averaged/subj. average structurs
load([resultsPath,'\allSubj_visual_MI_S_C_hemisCombined_',date,'.mat'])

%% Time parameters to compute mean/max
opts.tMin.early = 0; % s
opts.tMax.early = 0.5; % s
minFrame.early = find(times == opts.tMin.early);
maxFrame.early = find(times == opts.tMax.early);
opts.tMin.late = 0.5; % s
opts.tMax.late = 1; % s
minFrame.late = find(times == opts.tMin.late);
maxFrame.late = find(times == opts.tMax.late);
opts.tMin.plot = -0.1; % s
opts.tMax.plot = 1.1; % s
minFrame.plot = find(times == opts.tMin.plot);
maxFrame.plot = find(times == opts.tMax.plot);

selSubgroups = subgroups_all(selROIs);

dls_rois = find(selSubgroups == 2); % Dorsolateral stream regions indxs
mt_rois = find(selSubgroups == 4);  % MT+ complex stream regions indxs
dls_mt_rois = [dls_rois;mt_rois];
DLS_MT_labels = selROIs_labels(dls_mt_rois);
selSubgroups = selSubgroups(dls_mt_rois);

time_label={'early','late'};

for ceLab = correrrLabels
    for iIdx = 1:numel(info_type)
        iLab = info_type{iIdx};
        for hIdx = 1:numel(hemis_labels)
            hLab = hemis_labels{hIdx};
            % MI structure for dorsolateral and MT+ complex regions
            avgMI_dls_mt_maps.(ceLab{1}).(iLab).(hLab) = zeros(15,4,numel(dls_mt_rois),opts.timeSteps);
            % MI structure for all visual regions
            avgMI_maps.(ceLab{1}).(iLab).(hLab) = zeros(15,4,nROIs,opts.timeSteps);
        end
    end

    % Compute averages across subj
    for subIdx = 1:nSubj
        for iIdx = 1:numel(info_type)
            iLab = info_type{iIdx};
            for hIdx = 1:numel(hemis_labels)
                hLab = hemis_labels{hIdx};
                avgMI_dls_mt_maps.(ceLab{1}).(iLab).(hLab)(subIdx,:,:,:) = MI_maps.(ceLab{1}).(iLab).(hLab){subIdx}(:,dls_mt_rois,:);
                avgMI_maps.(ceLab{1}).(iLab).(hLab)(subIdx,:,:,:) = MI_maps.(ceLab{1}).(iLab).(hLab){subIdx}(:,:,:);
            end
        end
    end

    avgMI_dls_mt_maps.(ceLab{1}).S.hemisAvg = 0.5*(avgMI_dls_mt_maps.(ceLab{1}).S.right+avgMI_dls_mt_maps.(ceLab{1}).S.left);
    avgMI_dls_mt_maps.(ceLab{1}).C.hemisAvg = 0.5*(avgMI_dls_mt_maps.(ceLab{1}).C.right+avgMI_dls_mt_maps.(ceLab{1}).C.left);
    avgMI_maps.(ceLab{1}).S.hemisAvg = 0.5*(avgMI_maps.(ceLab{1}).S.right+avgMI_maps.(ceLab{1}).S.left);
    avgMI_maps.(ceLab{1}).C.hemisAvg = 0.5*(avgMI_maps.(ceLab{1}).C.right+avgMI_maps.(ceLab{1}).C.left);
end

nerual_features = fieldnames(avgMI_dls_mt_maps.all.S)';
for tLab = time_label
    for hIdx = 1:numel(nerual_features)
        hLab = nerual_features{hIdx};
        % average across subj and sess in the selected time window
        sessInfoProfiles = squeeze(avgMI_dls_mt_maps.all.S.(hLab)(:,:,:,minFrame.(tLab{1}):maxFrame.(tLab{1})));
        sessInfoProfiles = reshape(sessInfoProfiles,[],numel(dls_mt_rois),numel(minFrame.(tLab{1}):maxFrame.(tLab{1})));
        
        % Find peak of mutual information over time
        subjAvjMI.(tLab{1}).(hLab) = squeeze(mean(sessInfoProfiles,1));
        [maxMI.(tLab{1}).(hLab)] = max(subjAvjMI.(tLab{1}).(hLab),[],2);
        [vals,rankIdx.(tLab{1}).max.(hLab)]=sort(maxMI.(tLab{1}).(hLab),'descend');

        sessMaxMI.(tLab{1}).(hLab) = max(sessInfoProfiles,[],3);
        sessMeanMI.(tLab{1}).(hLab) = mean(sessInfoProfiles,3);
    end
end
% Compute overall max
for hIdx = 1:numel(nerual_features)
    hLab = nerual_features{hIdx};
    % average across subj and sess in the selected time window
    sessInfoProfiles = squeeze(avgMI_dls_mt_maps.all.S.(hLab));
    sessInfoProfiles = reshape(sessInfoProfiles,[],numel(dls_mt_rois),opts.timeSteps);

    [overallMax.(hLab),overallMaxIdxs.(hLab)] = max(sessInfoProfiles,[],3);
end
%% Plot ranked ROIs in dorsolateral and MT+ cortex 
% (we discard primary areas for vicinity to V1 and vental stram because
% they are known not to be sensitive to contrast)

% Colors for the 4 visual subgroups
cols = [11, 13, 186; 0, 119, 182; 0, 180, 216; 144, 224, 239]./255;

cmap = jet(256);
% Make lowest one black
cmap(129,:) = 1;

% selSubgroups = subgroups_all(selROIs);
% selROIs_labels = selROIs_labels(sel_rois);

rankFeat = 'max'; % either 'max' or 'mean'
rankTiming = 'early';
%% Plot ranked Dorsolateral stream and MT+ complex regions

ranked_selROIsLab = DLS_MT_labels(rankIdx.(rankTiming).max.hemisAvg);
ranked_subgroups = selSubgroups(rankIdx.(rankTiming).max.hemisAvg)';

figure()
hold on
imagesc(zeros(size(sum(ranked_subgroups == mode(ranked_subgroups)),3)))

groupsVals = unique(ranked_subgroups);
count = zeros(1,numel(groupsVals));

for roiIdx = 1:numel(dls_mt_rois)
    roiLab = ranked_selROIsLab(roiIdx);
    tmpIdx = find(groupsVals == ranked_subgroups(roiIdx));
    count(tmpIdx) = count(tmpIdx)+1;
    text(tmpIdx,count(tmpIdx),roiLab{1},'color',cols(ranked_subgroups(roiIdx),:))
end

colormap(cmap);
xline([1.5,2.5])
xlim([0.5,2.5])
yline([.5:1:sum(ranked_subgroups == mode(ranked_subgroups))])
yline(30.5,'linewidth',2)
xlabel('Hemispheres feature')
ylabel('ROIs rank by stim. info peak')
ylim([.5,sum(ranked_subgroups == mode(ranked_subgroups))])
set(gca, 'XTick',1:2, 'XTickLabel', {'Dorsolateral stream','MT+ complex'},'FontSize',14)
set(gca,'YDir','normal')
%     title(['ROIs ranked by average TF II in ',selFeat],'FontSize',16)


%% Compute info timing metrics (plot info, compute rois carrying sig. info, compute info max and ramping)
% l_idxs = 1:2:opts.nROIs;
% r_idxs = 2:2:opts.nROIs;
subselROIs = {'V1','V3A','LO3'}; % subselection of top visual ROIs (V1 is the main visual info broadcaster, V3A is the region carrying maximum stim. info in DLS vis., LO3 is the 2nd regions carrying max stim. info in the MT+ complex (V3CD too close to V3A)

% Set params and paths for saving
params.save_folder = []; % name of folder to save figures
paths.figuresPath = []; % path to save figures
rois_save_label = subselROIs{1};
for i = 1:numel(subselROIs)-1
    rois_save_label = [rois_save_label,'_',subselROIs{i+1}];
end
paths.figuresPath = [paths.figuresPath,params.save_folder,'\',rois_save_label];

if ~exist(paths.figuresPath,'dir')
   mkdir(paths.figuresPath) 
end

params.save_figures = 1;

% Plot info time profiles
sel_info = 'S';
sel_hem = 'hemisAvg';

figPairs=figure('Position',[360,198,277,420]);
for roiIdx = 1:numel(subselROIs)
    subselROIsIdx = find(ismember(selROIs_labels,subselROIs(roiIdx)));
    
    subplot(numel(subselROIs),1,roiIdx)
    hold on
    roiProfiles = squeeze(avgMI_maps.all.(sel_info).(sel_hem)(:,:,subselROIsIdx,:));
    roiProfiles = reshape(roiProfiles, [], 135);
    meanProfile = mean(roiProfiles(:,minFrame.plot:maxFrame.plot),1);
    semProfile = std(roiProfiles(:,minFrame.plot:maxFrame.plot),[],1)/sqrt(size(roiProfiles,1));
    plot(times(:,minFrame.plot:maxFrame.plot),meanProfile)
    shadedErrorBar(times(:,minFrame.plot:maxFrame.plot),meanProfile,semProfile)
    title(subselROIs(roiIdx))
    xlabel('Time [s]')
    xlim([opts.tMin.plot,opts.tMax.plot])
end
sgtitle('Stim info profiles (sem sess)')

if params.save_figures
    fname = [paths.figuresPath,'\selROIs_time_profiles_',date,'.png'];
    saveas(figPairs,fname)
    fname = [paths.figuresPath,'\selROIs_time_profiles_',date,'.fig'];
    saveas(figPairs,fname)
    fname = [paths.figuresPath,'\selROIs_time_profiles_',date,'.svg'];
    saveas(figPairs,fname)
end

figPairs=figure('Position',[360,198,277,420]);
for roiIdx = 1:numel(subselROIs)
    subselROIsIdx = find(ismember(selROIs_labels,subselROIs(roiIdx)));
    
    subplot(numel(subselROIs),1,roiIdx)
    hold on
    roiProfiles = squeeze(mean(avgMI_maps.all.(sel_info).(sel_hem)(:,:,subselROIsIdx,:),2));
    roiProfiles = reshape(roiProfiles, [], 135);
    meanProfile = mean(roiProfiles(:,minFrame.plot:maxFrame.plot),1);
    semProfile = std(roiProfiles(:,minFrame.plot:maxFrame.plot),[],1)/sqrt(size(roiProfiles,1));
    plot(times(minFrame.plot:maxFrame.plot),meanProfile)
    shadedErrorBar(times(minFrame.plot:maxFrame.plot),meanProfile,semProfile)
    title(subselROIs(roiIdx))
    xlabel('Time [s]')
    xlim([opts.tMin.plot,opts.tMax.plot])
end
sgtitle('Stimulus information profiles (sem across subj)')

% Plot superimposed time profiles
figPairs=figure('Position',[281,370,560,218]);
cols = [0 0 0; 1 0 0; 0 1 0; 0 0 1];
for roiIdx = 1:numel(subselROIs)
    subselROIsIdx = find(ismember(selROIs_labels,subselROIs(roiIdx)));
    
    hold on
    roiProfiles = squeeze(avgMI_maps.all.(sel_info).(sel_hem)(:,:,subselROIsIdx,:));
    roiProfiles = reshape(roiProfiles, [], 135);
    meanProfile = mean(roiProfiles(:,minFrame.plot:maxFrame.plot),1);
    semProfile = std(roiProfiles(:,minFrame.plot:maxFrame.plot),[],1)/sqrt(size(roiProfiles,1));
    h(roiIdx) = plot(times(:,minFrame.plot:maxFrame.plot),meanProfile,'color',cols(roiIdx,:),'LineWidth',1.5);
    shadedErrorBar(times(:,minFrame.plot:maxFrame.plot),meanProfile,semProfile,'LineProps',{'color',cols(roiIdx,:)},'patchSaturation',0.05)
    xlabel('Time [s]')
    xlim([opts.tMin.plot,opts.tMax.plot])
end
legend([h(1), h(2), h(3), h(4)], subselROIs,'FontSize',8,'Location','NorthEast')

title('Stim info profiles (sem sess)')

if params.save_figures
    fname = [paths.figuresPath,'\selROIs_superimposed_profiles_',date,'.png'];
    saveas(figPairs,fname)
    fname = [paths.figuresPath,'\selROIs_superimposed_profiles_',date,'.fig'];
    saveas(figPairs,fname)
    fname = [paths.figuresPath,'\selROIs_superimposed_profiles_',date,'.svg'];
    saveas(figPairs,fname)
end
%% Correct vs error comparison

totCorr = []; totErr = [];
for roiIdx = 1:numel(subselROIs)
    subselROIsIdx = find(ismember(selROIs_labels,subselROIs(roiIdx)));

    corrInfo = squeeze(max(avgMI_maps.corr.(sel_info).(sel_hem)(:,:,subselROIsIdx,minFrame.(rankTiming):maxFrame.(rankTiming)),[],4));
    errInfo = squeeze(max(avgMI_maps.err.(sel_info).(sel_hem)(:,:,subselROIsIdx,minFrame.(rankTiming):maxFrame.(rankTiming)),[],4));
    
    corrInfo = corrInfo(:);
    totCorr = [totCorr; corrInfo];
    errInfo = errInfo(:);
    totErr = [totErr; errInfo];
    
end

figOverall=figure('Position',[137,233,284,376]);
hold on
% Plot overall bar
meanCorr = mean(totCorr);
semCorr = std(totCorr)/sqrt(numel(totCorr));

meanErr = mean(totErr);
semErr = std(totErr)/sqrt(numel(totErr));

scatter(1,meanCorr,12,'filled','MarkerFaceColor',[0 0.7 0])
scatter(2,meanErr,12,'filled','MarkerFaceColor',[0.7 0.7 0.7])
    
legend('Stimulus','Choice','Autoupdate','off','Location','SouthEast')

errorbar(1,meanCorr,semCorr,'Color',[0 0.7 0])
errorbar(2,meanErr,semErr,'Color',[0.7 0.7 0.7])
maxY = max([meanCorr+semCorr,meanErr+semErr]);

[a,p] = ttest(totCorr,totErr);
pvalues_plot(p,1.5,(5/4)*maxY,maxY/10,maxY/10,14,0.5,'k',0)
    
ylim([0,(3/2)*maxY])
xlim([0.5,2.5])
title(['Stimulus info correct vs error (max)'])

if params.save_figures  
    fname = [paths.figuresPath,'\stim_corr_vs_err_overall_max_',date,'.png'];
    saveas(figOverall,fname)
    fname = [paths.figuresPath,'\stim_corr_vs_err_overall_max_',date,'.fig'];
    saveas(figOverall,fname)
    fname = [paths.figuresPath,'\stim_corr_vs_err_overall_max_',date,'.svg'];
    figOverall.Renderer = 'painters';
    saveas(figOverall,fname)
end
%% Stim vs choice

totStim = []; totChoice = [];
for roiIdx = 1:numel(subselROIs)
    subselROIsIdx = find(ismember(selROIs_labels,subselROIs(roiIdx)));

    stimInfo = squeeze(max(avgMI_maps.all.S.(sel_hem)(:,:,subselROIsIdx,minFrame.(rankTiming):maxFrame.(rankTiming)),[],4));
    choiceInfo = squeeze(max(avgMI_maps.all.C.(sel_hem)(:,:,subselROIsIdx,minFrame.(rankTiming):maxFrame.(rankTiming)),[],4));
    
    stimInfo = stimInfo(:);
    totStim = [totStim; stimInfo];
    
    choiceInfo = choiceInfo(:);
    totChoice = [totChoice; choiceInfo];
    
end

figOverall=figure('Position',[137,233,284,376]);
hold on
% Plot overall bar
meanStim = mean(totStim);
semStim = std(totStim)/sqrt(numel(totStim));

meanChoice = mean(totChoice);
semChoice = std(totChoice)/sqrt(numel(totChoice));

scatter(1,meanStim,12,'filled','MarkerFaceColor',[0.2 0.6 1])
scatter(2,meanChoice,12,'filled','MarkerFaceColor',[0.7 0 0.7])
% bar(numel(subselROIs)+1+x_shift,meanStim,'FaceColor',[0.2 0.6 1])
% bar(numel(subselROIs)+1+x_shift+1,meanChoice,'FaceColor',[0.7 0 0.7])

legend('Stimulus','Choice','Autoupdate','off','Location','SouthEast')

errorbar(1,meanStim,semStim,'color',[0.2 0.6 1])
errorbar(2,meanChoice,semChoice,'color',[0.7 0 0.7])
maxY = max([meanStim+semStim,meanChoice+semChoice]);

[a,p] = ttest(totStim,totChoice);
pvalues_plot(p,1.5,(5/4)*maxY,maxY/10,maxY/10,14,0.5,'k',0)
    
ylim([0,(3/2)*maxY])
xlim([0.5,2.5])
title(['Stimulus vs choice info (max)'])

if params.save_figures   
    fname = [paths.figuresPath,'\stim_vs_choice_overall_max_',date,'.png'];
    saveas(figOverall,fname)
    fname = [paths.figuresPath,'\stim_vs_choice_overall_max_',date,'.fig'];
    saveas(figOverall,fname)
    fname = [paths.figuresPath,'\stim_vs_choice_overall_max_',date,'.svg'];
    figOverall.Renderer = 'painters';
    saveas(figOverall,fname)
end

%% Early vs late

totEarly = []; totLate = [];
for roiIdx = 1:numel(subselROIs)
    subselROIsIdx = find(ismember(selROIs_labels,subselROIs(roiIdx)));

    earlyInfo = squeeze(max(avgMI_maps.all.(sel_info).(sel_hem)(:,:,subselROIsIdx,minFrame.early:maxFrame.early),[],4));
    lateInfo = squeeze(max(avgMI_maps.all.(sel_info).(sel_hem)(:,:,subselROIsIdx,minFrame.late:maxFrame.late),[],4));

    earlyInfo = earlyInfo(:);
    totEarly = [totEarly; earlyInfo];
    
    lateInfo = lateInfo(:);
    totLate = [totLate; lateInfo];
    
end

figOverall=figure('Position',[137,233,284,376]);
hold on
% Plot overall bar
meanEarly = mean(totEarly);
semEarly = std(totEarly)/sqrt(numel(totEarly));

meanLate = mean(totLate);
semLate = std(totLate)/sqrt(numel(totLate));
scatter(1,meanEarly,12,'filled','MarkerFaceColor',[185 140 100]/256)
scatter(2,meanLate,12,'filled','MarkerFaceColor',[40 21 6]/256)
% bar(numel(subselROIs)+1+x_shift,meanEarly,'FaceColor',[0.2 0.6 1])
% bar(numel(subselROIs)+1+x_shift+1,meanLate,'FaceColor',[0.7 0 0.7])

legend('Early','Late','Autoupdate','off','Location','SouthEast')

errorbar(1,meanEarly,semEarly,'color',[185 140 100]/256)
errorbar(2,meanLate,semLate,'color',[40 21 6]/256)
maxY = max([meanEarly+semEarly,meanLate+semLate]);

[a,p] = ttest(totEarly,totLate);
pvalues_plot(p,1.5,(5/4)*maxY,maxY/10,maxY/10,14,0.5,'k',0)
    
ylim([0,(3/2)*maxY])
xlim([0.5,2.5])
title(['Early vs late info (max)'])

if params.save_figures   
    fname = [paths.figuresPath,'\early_vs_late_overall_max_',date,'.png'];
    saveas(figOverall,fname)
    fname = [paths.figuresPath,'\early_vs_late_overall_max_',date,'.fig'];
    saveas(figOverall,fname)
    fname = [paths.figuresPath,'\early_vs_late_overall_max_',date,'.svg'];
    figPairs.Renderer='painters';
    saveas(figOverall,fname)
end
