%% Script to reproduce all plots in Figure 5 of the main FIT text
clear all; close all;
% Notice: in the script FIT is often referred to as rDFI according to an 
% previous nomenclature we used for the measure
% In this version of the script we take maximum over delays (instead of
% mean)
% Added max between null hypotheses and bootstrap

% Load results 
load('metadata.mat');

% Computed in computeMI
load('miEyeVisXEEG2bins.mat');

% Computed in analyzeTransferBothEyesContralateral
BothEyesTransfer = 'FIT_contralateral_conditioned';
BothEyesTransferSimpleNull = 'FIT_contralateral_simple';

% Computed in analyzeTransferBothEyesIpsilateral
BothEyesNoTransfer = 'FIT_ipsilateral';
load(BothEyesNoTransfer, 'leftEyeInf', 'rightEyeInf');
ipsiLeftEye = leftEyeInf;
ipsiRightEye = rightEyeInf;
clear leftEyeInf rightEyeInf
load(BothEyesTransfer, 'leftEyeInf', 'rightEyeInf');
simpleNull = load(BothEyesTransferSimpleNull, 'leftEyeInf', 'rightEyeInf');
contraLeftEye = leftEyeInf;
contraRightEye = rightEyeInf;

%% Fig. 4C left and right

leftEyeInf.FIT = squeeze(leftEyeInf.FIT);
nSubjects = size(leftEyeInf.FIT, 1);

plotAllSubjectsAverage(leftEyeInf.FIT, leftEyeRightSensorMI, leftEyeLeftSensorMI, time, 'Left Eye - ROT to LOT');

% figure;
rightEyeInf.FIT = squeeze(rightEyeInf.FIT);

plotAllSubjectsAverage(rightEyeInf.FIT, rightEyeLeftSensorMI, rightEyeRightSensorMI, time, 'Right Eye - LOT to ROT');

%% Fig. 4D (Left Eye)
rng(1)
clusteringParams = [0.01, 0.99]; % params for the cluster permutation test

% Compute the cluster stat with all surrogates put together
delaymax = 60;
nShuff = 10;
nBoot = 500;

sel_shuff = randperm(size(leftEyeInf.FITSh,5));
FIT = squeeze(leftEyeInf.FIT(:,:,1:delaymax));
FITcondSh = squeeze(leftEyeInf.FITSh(1,:,:,1:delaymax,sel_shuff(1:nShuff)));
FITsimpleSh = squeeze(simpleNull.leftEyeInf.FITSh(1,:,:,1:delaymax,1:nShuff));
FITcondShBtsp = btstrp_shuff(FITcondSh,nBoot);
FITsimpleShBtsp = btstrp_shuff(FITsimpleSh,nBoot);

infQuant = squeeze(mean(FIT, 1));
infQuantSh = max(cat(4,FITcondShBtsp,FITsimpleShBtsp),[],4);
%infQuantSh = squeeze(mean(FITsimpleSh,1));

% Use averages over subjects surrogate data as null hypotheses
mask = clusterStat_v2(infQuant, infQuantSh, clusteringParams(1), clusteringParams(2));

% Plot 

figure('units','normalized','outerposition',[0 0 1 1]);

i = imagesc(squeeze(infQuant(151:351,1:delaymax))');
i.AlphaData = mask(151:351, 1:delaymax)';
set(gca,'XTick',[0 50 100 150]);
set(gca,'XTickLabel',[time(151) time(201) time(251) time(301)],'fontsize',14);
xlabel('Time (ms)','fontsize', 16);

yTickLabels = get(gca, 'YTickLabel');
for i = 1:length(yTickLabels)
    yTickLabels{i} = 2 * str2num(yTickLabels{i});
end
 ylabel('delay (ms)','fontsize', 16);
yTicks = delaymax - get(gca, 'YTick');
yTicks = [delaymax yTicks(1:(length(yTicks) - 1))];
yTicks = sort(yTicks);
yTicks = [ yTicks(1), yTicks(3), yTicks(5) ];
yTickLabels = [ yTickLabels(1), yTickLabels(3), yTickLabels(5),'fontsize',14];
set(gca, 'YTick', yTicks);
set(gca,'YTickLabel',yTickLabels);
c = colorbar;
ylabel(c, '(bits)','fontsize',14);
ca = caxis;
caxis([0 ca(2)]);
title('Significant LE FIT','fontsize', 14);

leftMask = mask;

%% Fig. 4D (Right Eye)
% The same for Right Eye

sel_shuff = randperm(size(rightEyeInf.FITSh,5));
FIT = squeeze(rightEyeInf.FIT(:,:,1:delaymax));
FITcondSh = squeeze(rightEyeInf.FITSh(1,:,:,1:delaymax,sel_shuff(1:nShuff)));
FITsimpleSh = squeeze(simpleNull.rightEyeInf.FITSh(1,:,:,1:delaymax,1:nShuff));
FITcondShBtsp = btstrp_shuff(FITcondSh,nBoot);
FITsimpleShBtsp = btstrp_shuff(FITsimpleSh,nBoot);

infQuant = squeeze(mean(FIT, 1));
infQantSh = max(cat(4,FITcondShBtsp,FITsimpleShBtsp),[],4);

% Use averages over subjects surrogate data as null hypotheses
mask = clusterStat_v2(infQuant, infQuantSh, clusteringParams(1), clusteringParams(2));

% Plot 

delaymax = 60;
figure('units','normalized','outerposition',[0 0 1 1]);

i = imagesc(squeeze(infQuant(151:351,1:delaymax))');
i.AlphaData = mask(151:351, 1:delaymax)';
set(gca,'XTick',[0 50 100 150]);
set(gca,'XTickLabel',[time(151) time(201) time(251) time(301)]);
xlabel('Time (ms)');

yTickLabels = get(gca, 'YTickLabel');
for i = 1:length(yTickLabels)
    yTickLabels{i} = 2 * str2num(yTickLabels{i});
end
 ylabel('delay (ms)','fontsize', 16);
yTicks = delaymax - get(gca, 'YTick');
yTicks = [delaymax yTicks(1:(length(yTicks) - 1))];
yTicks = sort(yTicks);
yTicks = [ yTicks(1), yTicks(3), yTicks(5) ];
yTickLabels = [ yTickLabels(1), yTickLabels(3), yTickLabels(5) ];

set(gca, 'YTick', yTicks);
set(gca,'YTickLabel',yTickLabels);
c = colorbar;
ylabel(c, '(bits)');
ca = caxis;
caxis([0 ca(2)]);
title('Significant RE FIT','fontsize',14);

rightMask = mask;

%% Fig. 4B
% Plot FIT for all directions of transfer averaged across delays
clear leftEyeInf rightEyeInf

load(BothEyesTransfer, 'leftEyeInf', 'rightEyeInf');
leftEyeRightToLeft.FIT = squeeze(leftEyeInf.FIT);
rightEyeLeftToRight.FIT = squeeze(rightEyeInf.FIT);
leftEyeRightToLeft.FITSh = squeeze(leftEyeInf.FITSh);
rightEyeLeftToRight.FITSh = squeeze(rightEyeInf.FITSh);

clear leftEyeInf rightEyeInf

load(BothEyesNoTransfer, 'leftEyeInf', 'rightEyeInf');
leftEyeLeftToRight.FIT = squeeze(leftEyeInf.FIT);
rightEyeRightToLeft.FIT = squeeze(rightEyeInf.FIT);
leftEyeLeftToRight.FITSh = squeeze(leftEyeInf.FITSh);
rightEyeRightToLeft.FITSh = squeeze(rightEyeInf.FITSh);


%%
params.tMin = 151; % Window of interest after the stimulus presentation
params.tMax = 296; % 
params.time = (params.tMin-151:2:2*(params.tMax-params.tMin)+params.tMin-151);
params.adjSigTime = 10; % Consecutive time points that have to be significant to be plotted (i.e. at least 10ms)
params.sign_threshold = 0.025;

colors = [0, 114/255, 189/255; 217/255, 83/255, 25/255; 237/255, 177/255, 32/255; 126/255, 47/255, 142/255];
numShuff = size(leftEyeLeftToRight.FITSh,4);
maxTime = size(leftEyeRightToLeft.FIT, 2);
avgFITLEr2l = zeros(maxTime, 1);
avgFITREl2r = zeros(maxTime, 1);
avgFITLEl2r = zeros(maxTime, 1);
avgFITREr2l = zeros(maxTime, 1);
stdFITLEr2l = zeros(maxTime, 1);
stdFITREl2r = zeros(maxTime, 1);
stdFITLEl2r = zeros(maxTime, 1);
stdFITREr2l = zeros(maxTime, 1);

% Compute average FIT temporal profiles
for t = 2:maxTime
    if t > delaymax
        avgFITLEr2l(t) = squeeze(mean(mean(leftEyeRightToLeft.FIT(:,t,1:delaymax), 3), 1));
        avgFITLEl2r(t) = squeeze(mean(mean(leftEyeLeftToRight.FIT(:,t,1:delaymax), 3), 1));
        avgFITREr2l(t) = squeeze(mean(mean(rightEyeRightToLeft.FIT(:,t,1:delaymax), 3), 1));
        avgFITREl2r(t) = squeeze(mean(mean(rightEyeLeftToRight.FIT(:,t,1:delaymax), 3), 1));
        
        stdFITLEr2l(t) = squeeze(std(mean(leftEyeRightToLeft.FIT(:,t,1:delaymax), 3), [], 1))/sqrt(nSubjects);
        stdFITLEl2r(t) = squeeze(std(mean(leftEyeLeftToRight.FIT(:,t,1:delaymax), 3), [], 1))/sqrt(nSubjects);
        stdFITREr2l(t) = squeeze(std(mean(rightEyeRightToLeft.FIT(:,t,1:delaymax), 3), [], 1))/sqrt(nSubjects);
        stdFITREl2r(t) = squeeze(std(mean(rightEyeLeftToRight.FIT(:,t,1:delaymax), 3), [], 1))/sqrt(nSubjects);
    else
        avgFITLEr2l(t) = squeeze(mean(mean(leftEyeRightToLeft.FIT(:,t,1:(t-1)), 3), 1));
        avgFITLEl2r(t) = squeeze(mean(mean(leftEyeLeftToRight.FIT(:,t,1:(t-1)), 3), 1));
        avgFITREr2l(t) = squeeze(mean(mean(rightEyeRightToLeft.FIT(:,t,1:(t-1)), 3), 1));
        avgFITREl2r(t) = squeeze(mean(mean(rightEyeLeftToRight.FIT(:,t,1:(t-1)), 3), 1));  
        
        stdFITLEr2l(t) = squeeze(std(mean(leftEyeRightToLeft.FIT(:,t,1:(t-1)), 3), [], 1))/sqrt(nSubjects);
        stdFITLEl2r(t) = squeeze(std(mean(leftEyeLeftToRight.FIT(:,t,1:(t-1)), 3), [], 1))/sqrt(nSubjects);
        stdFITREr2l(t) = squeeze(std(mean(rightEyeRightToLeft.FIT(:,t,1:(t-1)), 3), [], 1))/sqrt(nSubjects);
        stdFITREl2r(t) = squeeze(std(mean(rightEyeLeftToRight.FIT(:,t,1:(t-1)), 3), [], 1))/sqrt(nSubjects);
    end
end

% sigFITLEr2l = time_significance_mask(avgFITLEr2l(151:296)', avgFITLEr2lSh(151:296,:), params);
% sigFITLEl2r = time_significance_mask(avgFITLEl2r(151:296)', avgFITLEl2rSh(151:296,:), params);
% sigFITREr2l = time_significance_mask(avgFITREr2l(151:296)', avgFITREr2lSh(151:296,:), params);
% sigFITREl2r = time_significance_mask(avgFITREl2r(151:296)', avgFITREl2rSh(151:296,:), params);
options.error = 'sem';
options.handle = figure(1);
options.alpha = 0.1;
options.line_width = 1;
options.x_axis = 1:146;

width = 2.15;
% Plot
figure(1);
hold on;
sd = std(nonzeros(avgFITLEr2l));
mn = mean(nonzeros(avgFITLEr2l));
h(1)=plot((avgFITLEr2l(151:296)), 'linewidth', width, 'color', colors(1,:));
%sigTimes = find(sigFITLEr2l==1);
%plot(sigTimes+1,(avgFITLEr2l(sigTimes+151)),'o','color',colors(1,:), 'MarkerFaceColor',colors(1,:));
%plot(prctile(avgFITLEr2lSh(151:296,:),97.5,2),'--','linewidth',1.5,'color',colors(1,:))
%shadedErrorBar([],avgFITLEr2l(151:296),stdFITLEr2l(151:296),'lineProps',{'markerfacecolor','k'},'patchSaturation',0.1);
data=mean(leftEyeRightToLeft.FIT(:,151:296,1:delaymax), 3);
options.color_area = colors(1,:);
options.color_line = colors(1,:);

plot_areaerrorbar(data,options);

hold on
sd = std(nonzeros(avgFITLEl2r));
mn = mean(nonzeros(avgFITLEl2r));
h(2)=plot((avgFITLEl2r(151:296)), 'linewidth', width, 'color', colors(2,:));
data=mean(leftEyeLeftToRight.FIT(:,151:296,1:delaymax), 3);
options.color_area = colors(2,:);
options.color_line = colors(2,:);
plot_areaerrorbar(data,options);

hold on
sd = std(nonzeros(avgFITREr2l));
mn = mean(nonzeros(avgFITREr2l));
h(3)=plot((avgFITREr2l(151:296)), 'linewidth', width, 'color', colors(3,:));
data=mean(rightEyeRightToLeft.FIT(:,151:296,1:delaymax), 3);
options.color_area = colors(3,:);
options.color_line = colors(3,:);
plot_areaerrorbar(data,options);

hold on
sd = std(nonzeros(avgFITREl2r));
mn = mean(nonzeros(avgFITREl2r));
h(4)=plot((avgFITREl2r(151:296)), 'linewidth', width, 'color', colors(4,:));
data=mean(rightEyeLeftToRight.FIT(:,151:296,1:delaymax), 3);
options.color_area = colors(4,:);
options.color_line = colors(4,:);
plot_areaerrorbar(data,options);

%axis([0 200 -Inf Inf]);
set(gca,'XTick',[1 26 51 76 101 126]);
set(gca,'XTickLabel',[time(151) time(176) time(201) time(226) time(251) time(276)]);
ax = gca;
ax.FontSize = 18;
ylabel('(bits)', 'fontsize', 16);
xlabel('Time (ms)', 'fontsize', 16);
title('FIT over time', 'fontsize', 22);
lgd = legend([h(1) h(2) h(3) h(4)],'LE R\rightarrowL', 'LE L\rightarrowR', 'RE R\rightarrowL', 'RE L\rightarrowR');
lgd.FontSize = 16;
xlim([0,145])
    

%% Panel E
 
% Boxes used to pick information values (both DI and FIT) for the left and
% the right eye. Boxes have been chosen to include the significant clusters
% for both eyes.

sizes.TDbox = [221, 271; 10, 46]; % TimeBins of 2ms --> boxes = 

% Import di, dfi and FIT in the different conditions
clear leftEyeInf rightEyeInf

load(BothEyesTransfer, 'leftEyeInf', 'rightEyeInf');
FIT1 = squeeze(leftEyeInf.FIT);
dfi1 = squeeze(leftEyeInf.dfi);
diR2L = squeeze(leftEyeInf.di);

FIT3 = squeeze(rightEyeInf.FIT);
dfi3 = squeeze(rightEyeInf.dfi);
diL2R = squeeze(rightEyeInf.di);

clear leftEyeInf rightEyeInf
load(BothEyesNoTransfer, 'leftEyeInf', 'rightEyeInf');
FIT2 = squeeze(leftEyeInf.FIT);
dfi2 = squeeze(leftEyeInf.dfi);

FIT4 = squeeze(rightEyeInf.FIT);
dfi4 = squeeze(rightEyeInf.dfi);

% Initialize masks
selBox = sizes.TDbox;
    
% Pick an information value for each subject by averaging over delays and
% taking the maximum over time
for idx = 1:nSubjects
    leR2LMean(idx) = max(mean(squeeze(FIT1(idx, selBox(1,1):selBox(1,2),selBox(2,1):selBox(2,2))),2));
    leL2RMean(idx) = max(mean(squeeze(FIT2(idx, selBox(1,1):selBox(1,2),selBox(2,1):selBox(2,2))),2));
    reR2LMean(idx) = max(mean(squeeze(FIT4(idx, selBox(1,1):selBox(1,2),selBox(2,1):selBox(2,2))),2));
    reL2RMean(idx) = max(mean(squeeze(FIT3(idx, selBox(1,1):selBox(1,2),selBox(2,1):selBox(2,2))),2));
end

%% FIT panel E
select_subjects = ones(1,15);
select_subjects = logical(select_subjects);

figure()
subplot(1,2,1);

% MEAN +- SEM
MleR2L = mean(leR2LMean(select_subjects));
SleR2L = std(leR2LMean(select_subjects))/sqrt(numel(leR2LMean(select_subjects)));
leR2Lplot = [MleR2L-SleR2L, MleR2L, MleR2L+SleR2L];
 
MleL2R = mean(leL2RMean(select_subjects));
SleL2R = std(leL2RMean(select_subjects))/sqrt(numel(leL2RMean(select_subjects)));
leL2Rplot = [MleL2R-SleL2R, MleL2R, MleL2R+SleL2R];
 
MreR2L = mean(reR2LMean(select_subjects));
SreR2L = std(reR2LMean(select_subjects))/sqrt(numel(reR2LMean(select_subjects)));
reR2Lplot = [MreR2L-SreR2L, MreR2L, MreR2L+SreR2L];
 
%h=boxplot([leR2Lplot', leL2Rplot', reR2Lplot'], [zeros(1, numel(leR2Lplot)) ones(1, numel(leL2Rplot)) (ones(1, numel(reR2Lplot)) + 1)], 'Colors', 'rkb', 'Symbol', 'g+');
h = errorbar([1, 2, 3], [MleR2L, MleL2R, MreR2L],[SleR2L, SleL2R, SreR2L],'o','MarkerFaceColor','blue');
set(h,{'linew'},{2})
[h,p] = ttest(leR2LMean(select_subjects),leL2RMean(select_subjects));
diff1x2P = p;
[h,p] = ttest(leR2LMean(select_subjects),reR2LMean(select_subjects));
diff1x3P = p;


ps = [diff1x2P, diff1x3P];

hold on
maxY = max([leR2Lplot, leL2Rplot]);
pvalues_plot(diff1x2P,1.5,maxY,maxY/7,1,22,0.45,'k',0)
maxY = max([leR2Lplot, reR2Lplot]);
pvalues_plot(diff1x3P,2,maxY,maxY/7,2,22,0.9,'k',0)

title('FIT LE','fontsize',16);
set(gca, 'XTickLabel', {'LE R->L', 'LE L->R', 'RE R->L'});
xlim([0.5,3.5])
ylim([0, 0.0035])

subplot(1,2,2);

% MEAN +- SEM
MreL2R = mean(reL2RMean(select_subjects));
SreL2R = std(reL2RMean(select_subjects))/sqrt(numel(reL2RMean(select_subjects)));
reL2Rplot = [MreL2R-SreL2R, MreL2R, MreL2R+SreL2R];
 
MleL2R = mean(leL2RMean(select_subjects));
SleL2R = std(leL2RMean(select_subjects))/sqrt(numel(leL2RMean(select_subjects)));
leL2Rplot = [MleL2R-SleL2R, MleL2R, MleL2R+SleL2R];
 
MreR2L = mean(reR2LMean(select_subjects));
SreR2L = std(reR2LMean(select_subjects))/sqrt(numel(reR2LMean(select_subjects)));
reR2Lplot = [MreR2L-SreR2L, MreR2L, MreR2L+SreR2L];
 
%h=boxplot([reL2Rplot', leL2Rplot', reR2Lplot'], [zeros(1, numel(reL2Rplot)) ones(1, numel(leL2Rplot)) (ones(1, numel(reR2Lplot)) + 1)], 'Colors', 'rkb', 'Symbol', 'g+');
h = errorbar([1, 2, 3], [MreL2R, MleL2R, MreR2L],[SreL2R, SleL2R, SreR2L],'o','MarkerFaceColor','blue');
set(h,{'linew'},{2})

[h,p] = ttest(reL2RMean(select_subjects),leL2RMean(select_subjects));
diff1x2P = p;
[h,p] = ttest(reL2RMean(select_subjects),reR2LMean(select_subjects));
diff1x3P = p;

ps = [diff1x2P, diff1x3P];

hold on
maxY = max([reL2Rplot, leL2Rplot]);
pvalues_plot(diff1x2P,1.5,maxY,maxY/5,1,22,0.45,'k',0)
maxY = max([reL2Rplot, reR2Lplot]);
pvalues_plot(diff1x3P,2,maxY,maxY/5,2,22,0.9,'k',0)

title('FIT RE','fontsize',16);
set(gca, 'XTickLabel', {'RE L->R', 'LE L->R', 'RE R->L'});
xlim([0.5,3.5])
ylim([0, 0.0035])

%% DI panel
for idx = 1:nSubjects
    R2LMean(idx) = max(mean(squeeze(diR2L(idx, selBox(1,1):selBox(1,2),selBox(2,1):selBox(2,2))),2));
    L2RMean(idx) = max(mean(squeeze(diL2R(idx, selBox(1,1):selBox(1,2),selBox(2,1):selBox(2,2))),2));
end

figure('Position',[360,198,291,420])
% MEAN +- SEM
MR2L = mean(R2LMean);
SR2L = std(R2LMean)/sqrt(numel(R2LMean));
R2Lplot = [MR2L-SR2L, MR2L, MR2L+SR2L];
 
ML2R = mean(L2RMean);
SL2R = std(L2RMean)/sqrt(numel(L2RMean));
L2Rplot = [ML2R-SL2R, ML2R, ML2R+SL2R];
 
%h=boxplot([leR2Lplot', leL2Rplot', reR2Lplot'], [zeros(1, numel(leR2Lplot)) ones(1, numel(leL2Rplot)) (ones(1, numel(reR2Lplot)) + 1)], 'Colors', 'rkb', 'Symbol', 'g+');
h = errorbar([1, 2], [MR2L, ML2R],[SR2L, SL2R],'o','MarkerFaceColor','blue');
set(h,{'linew'},{2})

[h,p] = ttest(R2LMean,L2RMean);
diff1x2P = p;

ps = [diff1x2P];

hold on
maxY = max([R2Lplot, L2Rplot]);
pvalues_plot(diff1x2P,1.5,maxY,maxY/12,1,22,0.45,'k',0)

title('DI','fontsize',16);
set(gca, 'XTick', [1 2], 'XTickLabel', {'R->L', 'L->R'});
xlim([0.5,2.5])
ylim([0,0.068])

%% DFI panel (Fig. S15C)

for idx = 1:nSubjects
    leR2LMean(idx) = max(mean(squeeze(dfi1(idx, selBox(1,1):selBox(1,2),selBox(2,1):selBox(2,2))),2));
    leL2RMean(idx) = max(mean(squeeze(dfi2(idx, selBox(1,1):selBox(1,2),selBox(2,1):selBox(2,2))),2));
    reR2LMean(idx) = max(mean(squeeze(dfi3(idx, selBox(1,1):selBox(1,2),selBox(2,1):selBox(2,2))),2));
    reL2RMean(idx) = max(mean(squeeze(dfi4(idx, selBox(1,1):selBox(1,2),selBox(2,1):selBox(2,2))),2));
end

select_subjects = ones(1,15);
select_subjects = logical(select_subjects);

figure()
subplot(1,2,1);

% MEAN +- SEM
MleR2L = mean(leR2LMean(select_subjects));
SleR2L = std(leR2LMean(select_subjects))/sqrt(numel(leR2LMean(select_subjects)));
leR2Lplot = [MleR2L-SleR2L, MleR2L, MleR2L+SleR2L];
 
MleL2R = mean(leL2RMean(select_subjects));
SleL2R = std(leL2RMean(select_subjects))/sqrt(numel(leL2RMean(select_subjects)));
leL2Rplot = [MleL2R-SleL2R, MleL2R, MleL2R+SleL2R];
 
MreR2L = mean(reR2LMean(select_subjects));
SreR2L = std(reR2LMean(select_subjects))/sqrt(numel(reR2LMean(select_subjects)));
reR2Lplot = [MreR2L-SreR2L, MreR2L, MreR2L+SreR2L];
 
%h=boxplot([leR2Lplot', leL2Rplot', reR2Lplot'], [zeros(1, numel(leR2Lplot)) ones(1, numel(leL2Rplot)) (ones(1, numel(reR2Lplot)) + 1)], 'Colors', 'rkb', 'Symbol', 'g+');
h = errorbar([1, 2, 3], [MleR2L, MleL2R, MreR2L],[SleR2L, SleL2R, SreR2L],'o','MarkerFaceColor','blue');
set(h,{'linew'},{2})
[h,p] = ttest(leR2LMean(select_subjects),leL2RMean(select_subjects));
diff1x2P = p;
[h,p] = ttest(leR2LMean(select_subjects),reR2LMean(select_subjects));
diff1x3P = p;


ps = [diff1x2P, diff1x3P];

hold on
maxY = max([leR2Lplot, leL2Rplot]);
pvalues_plot(diff1x2P,1.5,maxY,maxY/7,1,22,0.45,'k',0)
maxY = max([leR2Lplot, reR2Lplot]);
pvalues_plot(diff1x3P,2,maxY,maxY/7,2,22,0.9,'k',0)

title('DFI LE','fontsize',16);
set(gca, 'XTickLabel', {'LE R->L', 'LE L->R', 'RE R->L'});
xlim([0.5,3.5])
ylim([-0.03,0])

subplot(1,2,2);

% MEAN +- SEM
MreL2R = mean(reL2RMean(select_subjects));
SreL2R = std(reL2RMean(select_subjects))/sqrt(numel(reL2RMean(select_subjects)));
reL2Rplot = [MreL2R-SreL2R, MreL2R, MreL2R+SreL2R];
 
MleL2R = mean(leL2RMean(select_subjects));
SleL2R = std(leL2RMean(select_subjects))/sqrt(numel(leL2RMean(select_subjects)));
leL2Rplot = [MleL2R-SleL2R, MleL2R, MleL2R+SleL2R];
 
MreR2L = mean(reR2LMean(select_subjects));
SreR2L = std(reR2LMean(select_subjects))/sqrt(numel(reR2LMean(select_subjects)));
reR2Lplot = [MreR2L-SreR2L, MreR2L, MreR2L+SreR2L];
 
%h=boxplot([reL2Rplot', leL2Rplot', reR2Lplot'], [zeros(1, numel(reL2Rplot)) ones(1, numel(leL2Rplot)) (ones(1, numel(reR2Lplot)) + 1)], 'Colors', 'rkb', 'Symbol', 'g+');
h = errorbar([1, 2, 3], [MreL2R, MleL2R, MreR2L],[SreL2R, SleL2R, SreR2L],'o','MarkerFaceColor','blue');
set(h,{'linew'},{2})

[h,p] = ttest(reL2RMean(select_subjects),leL2RMean(select_subjects));
diff1x2P = p;
[h,p] = ttest(reL2RMean(select_subjects),reR2LMean(select_subjects));
diff1x3P = p;

ps = [diff1x2P, diff1x3P];

hold on
maxY = max([reL2Rplot, leL2Rplot]);
pvalues_plot(diff1x2P,1.5,maxY,maxY/5,1,22,0.45,'k',0)
maxY = max([reL2Rplot, reR2Lplot]);
pvalues_plot(diff1x3P,2,maxY,maxY/5,2,22,0.9,'k',0)

title('DFI RE','fontsize',16);
set(gca, 'XTickLabel', {'RE L->R', 'LE L->R', 'RE R->L'});
xlim([0.5,3.5])
ylim([-0.03,0])
