

% for 600 DPI
set(0, 'defaultlinelinewidth', 1.5);
set(0, 'defaulterrorbarlinewidth', 1.5);
set(0, 'defaultaxeslinewidth', 1.2);



global defaultfontsize
defaultfontsize = 14;

set(0,'defaultaxesfontsize',defaultfontsize);



global epoch_fmri_1_rest epoch_fmri_1_pf epoch_fmri_1_uf epoch_fmri_2_rest epoch_fmri_2_pf epoch_fmri_2_uf
global epoch_fmrs_1_rest epoch_fmrs_1_pf epoch_fmrs_1_uf epoch_fmrs_2_rest epoch_fmrs_2_pf epoch_fmrs_2_uf
epoch_fmri_1_rest = [1, 30; 61, 90; 121, 150; 181, 210];
epoch_fmri_1_pf   = [31, 60; 151, 180];
epoch_fmri_1_uf   = [91, 120; 211, 240];
epoch_fmri_2_rest = epoch_fmri_1_rest;
epoch_fmri_2_pf   = epoch_fmri_1_uf;
epoch_fmri_2_uf   = epoch_fmri_1_pf;
epoch_fmrs_1_rest = [1, 240; 481, 720];
epoch_fmrs_1_pf   = [241, 480];
epoch_fmrs_1_uf   = [721, 960];
epoch_fmrs_2_rest = epoch_fmrs_1_rest;
epoch_fmrs_2_pf   = epoch_fmrs_1_uf;
epoch_fmrs_2_uf   = epoch_fmrs_1_pf;






global epoch_title
epoch_title = {'{\bf{fMRI}} Run 1', '{\bf{fMRI}} Run 2', '{\bf{1H-fMRS}} Run 1', '{\bf{1H-fMRS}} Run 2'};


load('./fig_1_data/eyetracker.mat');

load('./fig_1_data/data_performance_trials.mat', 'data_performance_trials'); 

load('./fig_1_data/data_performance.mat', 'data_performance'); 

%%

threshold = 4;  %  just a threshold for visualiaztion purposes (heatmap)

hm_epochs = squeeze(nanmean(hm)); % 3D mean across subjects for epochs(1,2,3)=(local,mrs1,mrs2)

hm_total = squeeze(nanmean(hm_epochs)); % 2D mean across subjects and epochs



%%

global duration_fmri duration_fmrs

duration_fmri = 240;
duration_fmrs = 960;




%image = imread('stimulation_screenshot.png');
%image = imresize(image, [image_size, image_size]);

load('./fig_1_data/data_gaze.mat', 'data_gaze'); 

figure;
set(gcf, 'color', 'white');

subplot(3, 3, 1); 
[pos_y1, pos_h1] = plot_hm(hm_epochs(1, :, :), threshold);

subplot(3, 3, 2);
plot_gaze(1, pos_y1, pos_h1);

subplot(3, 3, 3);
plot_gaze_barplots(squeeze(data_gaze(1, :, :))); 


figure;
set(gcf, 'color', 'white');

subplot(3, 3, 1); 
[pos_y2, pos_h2] = plot_hm(hm_epochs(2, :, :), threshold);

subplot(3, 3, 2);
plot_gaze(2, pos_y2, pos_h2);

subplot(3, 3, 3);
plot_gaze_barplots(squeeze(data_gaze(2, :, :))); 



figure;
set(gcf, 'color', 'white');

subplot(3, 3, 1); 
[pos_y3, pos_h3] = plot_hm(hm_epochs(3, :, :), threshold);

subplot(3, 3, 2);
plot_gaze(4, pos_y3, pos_h3);

subplot(3, 3, 3);
plot_gaze_barplots(squeeze(data_gaze(3, :, :))); 


figure;
set(gcf, 'color', 'white');

subplot(3, 3, 1); 
[pos_y4, pos_h4] = plot_hm(hm_epochs(4, :, :), threshold);

subplot(3, 3, 2);
plot_gaze(4, pos_y4, pos_h4);

subplot(3, 3, 3);
plot_gaze_barplots(squeeze(data_gaze(4, :, :))); 



%%

load('./fig_1_data/data_pupil.mat', 'data_pupil'); 


figure;
set(gcf, 'color', 'white');

subplot(3, 2, 1); 
plot_pupil(1, pos_y1, pos_h1);

subplot(3, 2, 2);
plot_pupil_barplots(squeeze(data_pupil(1, :, :))); 

figure;
set(gcf, 'color', 'white');

subplot(3, 2, 1); 
plot_pupil(2, pos_y2, pos_h2);

subplot(3, 2, 2);
plot_pupil_barplots(squeeze(data_pupil(2, :, :))); 


figure;
set(gcf, 'color', 'white');

subplot(3, 2, 1); 
plot_pupil(3, pos_y3, pos_h3);

subplot(3, 2, 2);
plot_pupil_barplots(squeeze(data_pupil(3, :, :))); 


figure;
set(gcf, 'color', 'white');

subplot(3, 2, 1); 
plot_pupil(4, pos_y4, pos_h4);

subplot(3, 2, 2);
plot_pupil_barplots(squeeze(data_pupil(4, :, :))); 



% -------------------

lim_performance = [0, 800];
lim_pupil = [2, 15];
lim_gaze = [0, 1.5];

label_performance = 'Response Delay (ms)';
label_pupil = {'Pupil';'Diameter';'(mm)'};
label_gaze = {'Gaze';'Displacement';'(deg)'};


aux_2 = 0.91;

figure;
set(gcf, 'color', 'white');

subplot(3, 3, 1);
plot_performance_barplots(1, data_performance);
subplot(3, 3, 3);
plot_scatter(squeeze(data_performance(1, :, :)), squeeze(data_pupil(1, :, :)), ...
    lim_performance, lim_pupil, label_performance, label_pupil, aux_2);
subplot(3, 3, 2);
plot_scatter(squeeze(data_performance(1, :, :)), squeeze(data_gaze(1, :, :)), ...
    lim_performance, lim_gaze, label_performance, label_gaze, 1);

figure;
set(gcf, 'color', 'white');

subplot(3, 3, 1);
plot_performance_barplots(2, data_performance);
subplot(3, 3, 3);
plot_scatter(squeeze(data_performance(2, :, :)), squeeze(data_pupil(2, :, :)), ...
    lim_performance, lim_pupil, label_performance, label_pupil, aux_2);
subplot(3, 3, 2);
plot_scatter(squeeze(data_performance(2, :, :)), squeeze(data_gaze(2, :, :)), ...
    lim_performance, lim_gaze, label_performance, label_gaze, 1);

figure;
set(gcf, 'color', 'white');

subplot(3, 3, 1);
plot_performance_barplots(3, data_performance);
subplot(3, 3, 3);
plot_scatter(squeeze(data_performance(3, :, :)), squeeze(data_pupil(3, :, :)), ...
    lim_performance, lim_pupil, label_performance, label_pupil, aux_2);
subplot(3, 3, 2);
plot_scatter(squeeze(data_performance(3, :, :)), squeeze(data_gaze(3, :, :)), ...
    lim_performance, lim_gaze, label_performance, label_gaze, 1);

figure;
set(gcf, 'color', 'white');

subplot(3, 3, 1);
plot_performance_barplots(4, data_performance);
subplot(3, 3, 3);
plot_scatter(squeeze(data_performance(4, :, :)), squeeze(data_pupil(4, :, :)), ...
    lim_performance, lim_pupil, label_performance, label_pupil, aux_2);
subplot(3, 3, 2);
plot_scatter(squeeze(data_performance(4, :, :)), squeeze(data_gaze(4, :, :)), ...
    lim_performance, lim_gaze, label_performance, label_gaze, 1);


% ----------------------------------------------------------------


load('./fig_1_data/data_BOLD_in_voxel.mat', 'data_BOLD_in_voxel'); 
load('./fig_1_data/data_BOLD_in_atlas.mat', 'data_BOLD_in_atlas'); 


ticklabels = {
    ['100%',10,'UF (30 Hz)'], ...
    ['0%',10,'PF (7.5 Hz)'], ...
    ['52%',10,'PF (7.5 Hz)'], ...
    ['62%',10,'PF (7.5 Hz)'], ...
    ['70%',10,'PF (7.5 Hz)'], ...
    ['76%',10,'PF (7.5 Hz)'], ...
    ['82%',10,'PF (7.5 Hz)'], ...
    ['100%',10,'PF (7.5 Hz)'], ...
};

load('./fig_1_data/data_BOLD_var.mat', 'bold_UF', 'bold_PF', 'contrast');

bold_UF_mean = mean(bold_UF);
bold_UF_std = std(bold_UF);

bold_PF_mean = mean(bold_PF, 2);
bold_PF_std = std(bold_PF, 0, 2);


plot_col_PF = [0.2, 0.2, 0.2];
plot_col_UF = [0, 0.7, 0.8];


fit = polyfit(contrast', squeeze(bold_PF_mean), 1);


figure;
set(gcf, 'color', 'white');

subplot(2, 4, 1);
hold on;

range = 0:100;

myshadedplot(range, bold_UF_mean.*ones(1, numel(range)), bold_UF_std.*ones(1, numel(range)), plot_col_UF); 

contrast_rel_UF = (bold_UF_mean-fit(2))./fit(1);

%line([contrast_rel_UF, 100], [bold_UF_mean, bold_UF_mean], 'color', plot_col_UF, 'linestyle', '--');
line([contrast_rel_UF, contrast_rel_UF], [0,8], 'color', plot_col_UF, 'linestyle', '--','HandleVisibility','off');


errorbar(contrast, bold_PF_mean, bold_PF_std, 'color', plot_col_PF, 'marker', 's', 'markeredgecolor', 'k', 'markerfacecolor', 'w', 'linestyle', 'none','markersize',10);
%errorbar(100, bold_UF_mean, bold_UF_std, 'color', plot_col_UF, 'marker', 'o', 'markeredgecolor', plot_col_UF, 'markerfacecolor', plot_col_UF,'markersize',10);

ylim([0,8]);
xlim([-5, 105]);

xlabel('Contrast (%)');
ylabel('BOLD (%)');



plot(range, polyval(fit, range), 'color', [0.5, 0.5, 0.5], 'linestyle', '--');



%text(contrast_rel_UF-5, 0.85, sprintf('BOLD-matching at'), 'color', 'k', 'horizontalalignment', 'right');
%text(contrast_rel_UF-5, 0.45, ['\sim{\bf{', sprintf('%.0f%%', contrast_rel_UF), '}} contrast'], 'color', plot_col_UF, 'horizontalalignment', 'right', 'fontsize', 12);

legend({'UF @100%', 'PF', 'Linear Fit'}, 'location', 'northwest', 'fontsize', 12);



%%


subplot(2, 4, [2,3]);
plot_extracted_fmri_response;
pos = get(gca, 'position');
set(gca,'position',[1.05.*pos(1), pos(2), pos(3), pos(4)]);
ylim([-0.5, 1.5]);
text(15, 1.43, '{\bf{^1H-fMRS Voxel}}','color','k','horizontalalignment','center','fontsize',12);

subplot(3, 4, 11);
i1 = imread('./fig_1_data/mni_voxel.png');
imshow(i1);
pos = get(gca,'position');
set(gca,'position',[0.72.*pos(1), 7.*pos(2), 0.55.*pos(3), pos(4)]);


%%

fdark_col = 0.95; 

subplot(2, 5, 5);

plot_col_T = fdark_col.*[0.85, 0.9, 0.9];

nsubjects = 19;
data_BOLD_in_atlas(:, 1, 3, 1) = data_BOLD_in_atlas(:, 1, 3, 1)+0.1;
data_BOLD_in_atlas(:, 2, 3, 1) = data_BOLD_in_atlas(:, 2, 3, 1)-0.1;

[h, p] = ttest2(squeeze(data_BOLD_in_voxel(:, 2, 1)), squeeze(data_BOLD_in_voxel(:, 3, 1)));


for s=1:nsubjects
   
    hold on;
    plot(1:2, data_BOLD_in_voxel(s, 2:3, 1), ...
        'color', plot_col_T, 'marker', 's', 'markeredgecolor', plot_col_T, 'markerfacecolor', plot_col_T, ...
        'linestyle', '-');
    
    
end

data_BOLD_in_voxel_mean = squeeze(mean(data_BOLD_in_voxel, 1));
data_BOLD_in_voxel_std = squeeze(std(data_BOLD_in_voxel, 0, 1));

plot(1:2, data_BOLD_in_voxel_mean(2:3, 1), 'marker', 's', 'markersize', 10, 'markeredgecolor', plot_col_UF, 'markerfacecolor', plot_col_UF);
errorbar(1:2, data_BOLD_in_voxel_mean(2:3, 1), data_BOLD_in_voxel_std(2:3, 1), 'color', plot_col_UF);

fprintf('IN VOXEL');
data_BOLD_in_voxel_mean(2:3, 1)
data_BOLD_in_voxel_std(2:3, 1)

xlim([0.5, 2.5]);
ylim([-0.5, 2.5]);
set(gca,'xtick',1:2);
set(gca,'xticklabel',{'PF', 'UF'});
%xtickangle(30);
xlabel('Condition');
ylabel('BOLD (%)');

line([1, 2], [1.8, 1.8], 'color', 'k');
text(1.5, 1.8, sprintf('n.s.\n(p=%.2f)', p),'color','k','horizontalalignment','center','fontsize',12);

text(1.5, 2.4, '{\bf{^1H-fMRS Voxel}}','color','k','horizontalalignment','center','fontsize',12);




load('./fig_1_data/data_BOLD_in_atlas.mat', 'data_BOLD_in_atlas'); 
figure;

plot_col_PF = [0.2, 0.2, 0.2];
plot_col_UF = [0, 0.7, 0.8];
fdark_col = 0.95;

subplot(3, 6, 16);

plot_col_T = fdark_col.*[0.85, 0.9, 0.9];

nsubjects = 19;

[h, p] = ttest2(squeeze(data_BOLD_in_atlas(:, 1, 2, 1)), squeeze(data_BOLD_in_atlas(:, 1, 3, 1)));

for s=1:nsubjects
   
    hold on;
    plot(1:2, squeeze(data_BOLD_in_atlas(s, 1, 2:3, 1)), ...
        'color', plot_col_T, 'marker', 's', 'markeredgecolor', plot_col_T, 'markerfacecolor', plot_col_T, ...
        'linestyle', '-');
    
    
end

data_BOLD_in_atlas_mean = squeeze(mean(data_BOLD_in_atlas, 1));
data_BOLD_in_atlas_std = squeeze(std(data_BOLD_in_atlas, 0, 1));

plot(1:2, squeeze(data_BOLD_in_atlas_mean(1, 2:3, 1)), 'marker', 's', 'markersize', 10, 'markeredgecolor', plot_col_UF, 'markerfacecolor', plot_col_UF);
errorbar(1:2, squeeze(data_BOLD_in_atlas_mean(1, 2:3, 1)), squeeze(data_BOLD_in_atlas_std(1, 2:3, 1)), 'color', plot_col_UF);

fprintf('IN ATLAS BA17');
squeeze(data_BOLD_in_atlas_mean(1, 2:3, 1))
squeeze(data_BOLD_in_atlas_std(1, 2:3, 1))

xlim([0.5, 2.5]);
ylim([-1, 3]);
set(gca,'xtick',1:2);
set(gca,'xticklabel',{'PF', 'UF'});
%xtickangle(30);
xlabel('Condition');
ylabel('BOLD (%)');

line([1, 2], [2.3, 2.3], 'color', 'k');
text(1.5, 2.3, sprintf('n.s.\n(p=%.2f)', p),'color','k','horizontalalignment','center','fontsize',12);

text(1.5, 3.2, '{\bf{BA17 (V1)}}','color','k','horizontalalignment','center','fontsize',12);




subplot(3, 6, 18);

plot_col_T = fdark_col.*[0.85, 0.9, 0.9];

nsubjects = 19;

[h, p] = ttest2(squeeze(data_BOLD_in_atlas(:, 2, 2, 1)), squeeze(data_BOLD_in_atlas(:, 2, 3, 1)));


for s=1:nsubjects
   
    hold on;
    plot(1:2, squeeze(data_BOLD_in_atlas(s, 2, 2:3, 1)), ...
        'color', plot_col_T, 'marker', 's', 'markeredgecolor', plot_col_T, 'markerfacecolor', plot_col_T, ...
        'linestyle', '-');
    
    
end

data_BOLD_in_atlas_mean = squeeze(mean(data_BOLD_in_atlas, 1));
data_BOLD_in_atlas_std = squeeze(std(data_BOLD_in_atlas, 0, 1));

plot(1:2, squeeze(data_BOLD_in_atlas_mean(2, 2:3, 1)), 'marker', 's', 'markersize', 10, 'markeredgecolor', plot_col_UF, 'markerfacecolor', plot_col_UF);
errorbar(1:2, squeeze(data_BOLD_in_atlas_mean(2, 2:3, 1)), squeeze(data_BOLD_in_atlas_std(2, 2:3, 1)), 'color', plot_col_UF);

fprintf('IN ATLAS BA18/19');
squeeze(data_BOLD_in_atlas_mean(2, 2:3, 1))
squeeze(data_BOLD_in_atlas_std(2, 2:3, 1))

xlim([0.5, 2.5]);
ylim([-1, 3]);
set(gca,'xtick',1:2);
set(gca,'xticklabel',{'PF', 'UF'});
%xtickangle(30);
xlabel('Condition');
ylabel('BOLD (%)');

line([1, 2], [2.3, 2.3], 'color', 'k');
text(1.5, 2.3, ['{\ast}',sprintf('\n(p=%.3f)', p)],'color','k','horizontalalignment','center','fontsize',12);

text(1.5, 3.2, '{\bf{BA18/19}}','color','k','horizontalalignment','center','fontsize',12);

%%

i1 = imread('./fig_1_data/mni_ba17_d12D_91_39_70.png');
i2 = imread('./fig_1_data/mni_ba1819_d12D_91_39_70.png');

subplot(3, 6, 10);
imshow(i1);
pos = get(gca,'position');
set(gca,'position',[pos(1)+0.9.*pos(3), pos(2)-pos(4), 0.75.*pos(3), pos(4)]);

subplot(3, 6, 12);
imshow(i2);

pos = get(gca,'position');
set(gca,'position',[pos(1)+0.9.*pos(3), pos(2)-pos(4), 0.75.*pos(3), pos(4)]);


%--------------------


load('./fig_1_data/extracted_BOLD_invoxel.mat');

plot_col_PF = [0.2, 0.2, 0.2];
plot_col_UF = [0.6, 0.6, 0.6];
plot_col_X = [0, 0.7, 0.8];

% [pvalue, tbl, stats] = anova1([mean(Low); mean(High)]')


% remove 100%
High = High - 100;
Low = Low - 100;

pf_uf = {'PF', 'UF'};

[npoints, nsubjects] = size(Low); % same for High

% time
time_step = 1.5;
time_offset = -15;
time_duration = 60;
time = time_offset:time_step:(time_duration-time_offset);
time(end) = [];



% remove points before stimulus onset
remove = find(time<0);
Low(remove, :) = [];
High(remove, :) = [];

% ----------------

% find peak
pLow = max(Low);
pHigh = max(High);

% find ttp as fraction of peak
fraction = 0.05:0.05:0.45;
nfractions = numel(fraction);
onsetLow = nan(nfractions, nsubjects);
onsetHigh = nan(nfractions, nsubjects);
for f=1:nfractions
    for s=1:nsubjects
        aux = find(Low(:, s)>=fraction(f)*pLow(s));
        onsetLow(f, s) = aux(1);
        aux = find(High(:, s)>=fraction(f)*pHigh(s));
        onsetHigh(f, s) = aux(1);
    end
end
% mean and convert in secs
onsetLow_mean = mean(onsetLow, 2).*time_step;
onsetHigh_mean = mean(onsetHigh, 2).*time_step;
onsetLow_std = std(onsetLow, 0, 2).*time_step;
onsetHigh_std = std(onsetHigh, 0, 2).*time_step;


% [pvalue, tbl, stats] = anova1([onsetLow_mean, onsetHigh_mean])


% ------------

% find ttp
ttpLow = nan(1, nsubjects);
ttpHigh = nan(1, nsubjects);
for s=1:nsubjects
    aux = find(Low(:, s)>=pLow(s));
    ttpLow(s) = aux(1);
    aux = find(High(:, s)>=pHigh(s));
    ttpHigh(s) = aux(1);
end
% mean and convert in secs
ttpLow_mean = mean(ttpLow).*time_step;
ttpHigh_mean = mean(ttpHigh).*time_step;
ttpLow_std = std(ttpLow).*time_step;
ttpHigh_std = std(ttpHigh).*time_step;

[h, pvalue] = ttest(ttpLow, ttpHigh)

% --------------------------------------------
% plot

figure;

percent_fraction = fraction*100;

subplot(3, 5, 1);
hold on;
errorbar(percent_fraction, onsetLow_mean, onsetLow_std, nan([1, numel(onsetLow_std)]), 'o-', 'color', plot_col_PF, 'markerfacecolor', 'w', 'markeredgecolor', plot_col_PF, 'markersize', 8);
errorbar(percent_fraction, onsetHigh_mean, nan([1, numel(onsetHigh_std)]), onsetHigh_std, 'o-', 'color', plot_col_UF, 'markerfacecolor', 'w', 'markeredgecolor', plot_col_UF, 'markersize', 8);
set(gca, 'Layer', 'top');
hold off;
xlabel('Peak Fraction (%)', 'fontsize', defaultfontsize);
ylabel('Onset Time (s)', 'fontsize', defaultfontsize);
legend(pf_uf, 'location', 'northwest');
ylim([0, 20]);
xlim([0, 50]);
pos1 = get(gca, 'Position');
ax1 = gca;

set(gca, 'Layer', 'top');



subplot(3, 10, 3);
hold on;
%bar(1, ttpLow_mean, 'facecolor', 'none', 'edgecolor', plot_col_PF);
%bar(2, ttpHigh_mean, 'facecolor', 'none', 'edgecolor', plot_col_UF);
bar(1, ttpLow_mean, 'edgecolor', 'none', 'facecolor', plot_col_PF);
bar(2, ttpHigh_mean, 'edgecolor', 'none', 'facecolor', plot_col_UF);
errorbar(1, ttpLow_mean, ttpLow_std, 'color', plot_col_PF);
errorbar(2, ttpHigh_mean, ttpHigh_std, 'color', plot_col_UF);
hold off;
ylim([0, 40]);
xlim([0, 3]);
set(gca, 'XTick', 1:2);
set(gca, 'XTickLabel', pf_uf);
ylabel('Time-to-Peak (s)', 'fontsize', defaultfontsize);
xlabel('Condition', 'fontsize', defaultfontsize);
line([1, 2], [37, 37], 'color', 'k');
text(1.5, 40, 'n.s.', 'horizontalalignment', 'center', 'fontsize', 12);
pos2 = get(gca, 'Position');
set(ax1, 'Position', [pos1(1), pos2(2), 0.9.*pos1(3), pos2(4)]);

set(gca, 'Layer', 'top');



return



% ---- HELPER FUNCTIONS
%%

function plot_extracted_fmri_response


load('./fig_1_data/extracted_BOLD_invoxel.mat');

%figure; 
High = High - 100;
Low = Low -100;
% 
% High(:,6) = [];
% Low(:,6) = [];

mH = mean(High,2);
sH = std(High,[],2)./sqrt(size(High,2));

mL = mean(Low,2);
sL = std(Low,[],2)./sqrt(size(Low,2));

x = -15:1.5:45;
x(end) = [];

%SL = shadedErrorBar(x,mL,sL,'lineProps',{'color',[0.2,0.2,0.2],'markerfacecolor',[0.2,0.2,0.2]},'transparent',1,'patchSaturation',0.15);
%hold on;
%SH = shadedErrorBar(x,mH,sH,'lineProps',{'color',[0.6,0.6,0.6],'markerfacecolor',[0.6,0.6,0.6]},'transparent',1,'patchSaturation',0.10);


cond = (mL>mH);

sL1 = sL; sL2 = sL;
sL1(cond) = nan;
sL2(not(cond)) = nan;

sH1 = sH; sH2 = sH;
sH1(not(cond)) = nan;
sH2(cond) = nan;

%h=area(x, mL-mH, 'facecolor', [0.85, 0.9, 0.9], 'facealpha', 0.5,'edgecolor', 'none');
hold on;

errorbar(x, mL, sL1, sL2, 'color',[0.2,0.2,0.2], 'linestyle', '-', 'marker', 'o', 'markerfacecolor', 'w', 'markersize', 8);

errorbar(x, mH, sH1, sH2, 'color',[0.6,0.6,0.6], 'linestyle', '-', 'marker', 'o', 'markerfacecolor', 'w', 'markersize', 8);


set(gca, 'Layer', 'top');

xlim([-15,43.5]);
ylim([-0.5 1]);
plot_task([0.4 0.4 0.4]);

box on;
xlabel('Time (s)');
ylabel('BOLD (%)');


set(gca,'Layer','top');
legend({'PF','UF'}, 'location', 'northeast');

box off;

return
end

function plot_task(task_color)
stim_number = 1;
stim_start = [0];
width = 30; %28.5;
y_lim = ylim;
wei_y = 0.2; %y_lim(2) - y_lim(1);
for n = 1:stim_number
        pos = [stim_start(n) y_lim(1) width (wei_y)];
        r = rectangle('Position',pos,'FaceColor',task_color ,'LineStyle','none');
        uistack(r,'bottom')    
end
text(stim_start+width./2, y_lim(1)+wei_y./2, '{\bf{PF}} (7.5 Hz) or {\bf{UF}} (30 Hz)', 'color', 'w', 'fontsize', 14,'horizontalalignment','center');

return
end



% x, y, e(error), c(color)
% shaded area will be from y-e to y+e

function myshadedplot(x, y, e, c)

hold on;

h = area(x, [y-e; 2.*e]','HandleVisibility','off');

plot(x, y, 'Color', c);

set(h(1), 'FaceColor', 'none');
set(h, 'LineStyle', 'none');

ca = c+0.85*(1-c);
ca(ca>0.9) = 0.9;

set(h(2), 'FaceColor', ca);

set(gca, 'Layer', 'top');

return
end


%%

function plot_scatter(a, b, alim, blim, alabel, blabel, aux)

plot_col = [0, 0.7, 0.8; 0.2, 0.2, 0.2; 0.6, 0.6, 0.6];

hold on;
for i=1:3
    da = squeeze(a(:, i));
    db = squeeze(b(:, i));
    plot(da, db, 'o', 'color', plot_col(i, :));
    %{
    [~, tfa] = rmoutliers(da);
    [~, tfb] = rmoutliers(db);
    tf = not(tfa | tfb);
    da = da(tf);
    db = db(tf);
    %}
    [rho, pval] = corrcoef(da, db, 'Rows', 'pairwise');
end

xlim(alim);
ylim(blim);

xlabel(alabel);%, 'fontsize', 12);
ylabel(blabel);%, 'fontsize', 9);

legend({'Rest', 'PF', 'UF'}, 'orientation', 'horizontal');
%legend boxoff;

pos = get(gca, 'position');
set(gca, 'position', [aux.*pos(1), pos(2), 0.7.*pos(3), pos(4)]);


return
end


%%
function plot_performance_barplots(idx, data)

data = squeeze(data(idx, :, :));

[nsubjects, nsessions] = size(data);

[p,~,stats] = anova1(data, {'rest','pf','uf'}, 'off');
%comp = multcompare(stats, 'display', 'off')

plot_col_p = [0.85, 0.9, 0.9];
plot_col_m = [0.2, 0.2, 0.2];
plot_col_z = [0, 0.7, 0.8];


for e=1:nsessions

    d = data(:, e);
    
    dm = nanmean(squeeze(d));
    ds = nanstd(squeeze(d));   
    
    bar(e, dm, 0.65, 'edgecolor', plot_col_m, 'facecolor', plot_col_p); hold on;
    plot(e.*ones(1, nsubjects), d, 'o', 'color', plot_col_z);    
    errorbar(e, dm, ds, 'color', plot_col_m);
           
end

xlim([0.25, nsessions+0.75]);
set(gca,'xtick',1:nsessions);
set(gca,'xticklabel',{'Rest', 'PF', 'UF'}); xlabel('Epoch');
ylim([0, 1500]);
ylabel({'Response';'Delay';'(ms)'});
box off;

text(1, 1300, sprintf('One-Way ANOVA'), 'color', 'k', 'fontsize', 8);
text(1, 1150, sprintf('(df=%d, p=%.4f)', stats.df, p), 'color', 'k', 'fontsize', 10);

pos = get(gca, 'position');
set(gca, 'position', [2.*pos(1), pos(2), 0.4.*pos(3), pos(4)]);


return
end


%%

function plot_gaze_barplots(data)

[nsubjects, nsessions] = size(data);

[p,~,stats] = anova1(data, {'rest','pf','uf'}, 'off');


plot_col_p = [0.85, 0.9, 0.9];
plot_col_m = [0.2, 0.2, 0.2];
plot_col_z = [0, 0.7, 0.8];


for e=1:nsessions

    d = data(:, e);
    
    dm = nanmean(squeeze(d));
    ds = nanstd(squeeze(d));   
    
    bar(e, dm, 0.65, 'edgecolor', plot_col_m, 'facecolor', plot_col_p); hold on;
    plot(e.*ones(1, nsubjects), d, 'o', 'color', plot_col_z);    
    errorbar(e, dm, ds, 'color', plot_col_m);
           
end

global defaultfontsize

xlim([0.25, nsessions+0.75]);
set(gca,'xtick',1:nsessions);
set(gca,'xticklabel',{'Rest', 'PF', 'UF'}); 
xlabel('Epoch','fontsize',defaultfontsize);
ylim([0, 4]);
ylabel({'Gaze';'Displacement';'(deg)'},'fontsize',defaultfontsize);
box off;

text(1, 2.5, sprintf('One-Way ANOVA'), 'color', 'k', 'fontsize', 8);
text(1, 2, sprintf('(df=%d, p=%.4f)', stats.df, p), 'color', 'k', 'fontsize', 10);

pos = get(gca, 'position');
set(gca, 'position', [pos(1), pos(2), 0.4.*pos(3), pos(4)]);

return
end


%%

function plot_pupil_barplots(data)


[p,~,stats] = anova1(data, {'rest','pf','uf'}, 'off');


[nsubjects, nsessions] = size(data);

plot_col_p = [0.85, 0.9, 0.9];
plot_col_m = [0.2, 0.2, 0.2];
plot_col_z = [0, 0.7, 0.8];


for e=1:nsessions

    d = data(:, e);
    
    dm = nanmean(squeeze(d));
    ds = nanstd(squeeze(d));   
    
    bar(e, dm, 0.65, 'edgecolor', plot_col_m, 'facecolor', plot_col_p); hold on;
    plot(e.*ones(1, nsubjects), d, 'o', 'color', plot_col_z);    
    errorbar(e, dm, ds, 'color', plot_col_m);
    
    
end

xlim([0.25, nsessions+0.75]);
set(gca,'xtick',1:nsessions);
set(gca,'xticklabel',{'Rest', 'PF', 'UF'}); xlabel('Epoch');
ylim([0, 15]);
ylabel({'Pupil';'Diameter';'(mm)'});
box off;

text(1, 13.5, sprintf('One-Way ANOVA'), 'color', 'k', 'fontsize', 8);
text(1, 12, sprintf('(df=%d, p=%.4f)', stats.df, p), 'color', 'k', 'fontsize', 10);

pos = get(gca, 'position');
set(gca, 'position', [0.93.*pos(1), pos(2), 0.25.*pos(3), pos(4)]);

return
end

%%

function [pos_y, pos_h] = plot_hm(hm, threshold)
global field_deg 
image_size = 800; % pixel X pixel
tick_position = [1, 400, 800];
tick_label = [-field_deg/2, 0, field_deg/2];
% set(gca,'ColorScale','log')
alpha = 1;
isigma = 1.25;
c_location = 'EastOutside';
c_string = 'Time (%)';
%
hm = imgaussfilt(hm, isigma);
hm = squeeze(hm./nansum(hm(:)).*100);
hm = imresize(hm, [image_size, image_size]);
colormap('jet');
h = imagesc(hm);
set(gca, 'CLim', [0, threshold]);
daspect([1 1 1]);
set(gca,'xtick',tick_position); set(gca,'xticklabel',tick_label);
set(gca,'ytick',tick_position); set(gca,'yticklabel',-tick_label);
h.AlphaData = alpha.*ones(size(hm));
c = colorbar('location', c_location);
c.Label.String = c_string;
c.Label.FontSize = 12;
xlabel('Position (deg)');
ylabel('Position (deg)');
grid on;
set(gca, 'gridcolor', [1,1,1]);
pos = get(gca, 'position');
set(gca,'position', [1.4.*pos(1), pos(2), pos(3), pos(4)]);
pos_y = pos(2);
pos_h = pos(4);
return
end

%%

function plot_pupil(idx, pos_y, pos_h)
global epoch_title
%
name = {'fmri_run1', 'fmri_run2', '1hfmrs_run1', '1hfmrs_run2'};
load(sprintf('./fig_1_data/pupil_mean_sd_%s.mat', name{idx}), 'time', 'pupil', 'pupil_sd'); 
%plot
plotcol = [0, 0.7, 0.8];
myshadedplot(time, pupil, pupil_sd, plotcol); 
legend(epoch_title{idx});
legend boxoff;
ylim([0, 15]);
xlim([0, max(time)]);
ylabel({'Pupil';'Diameter';'(mm)'});
xlabel('Time (s)');
plot_epochs(idx, [0, 2], 8);
pos = get(gca,'position');
if not(isempty(pos_y))
    set(gca,'position',[pos(1), pos_y, pos(3), pos_h]);
end
return
end

%%

function plot_gaze(idx, pos_y, pos_h)
global epoch_title
%
name = {'fmri_run1', 'fmri_run2', '1hfmrs_run1', '1hfmrs_run2'};
load(sprintf('./fig_1_data/gaze_mean_sd_%s.mat', name{idx}), 'time', 'gaze', 'gaze_sd'); 
% plot
plotcol = [0, 0.7, 0.8];
myshadedplot(time, gaze, gaze_sd, plotcol); 
legend(epoch_title{idx});
legend boxoff;
ylim([0, 4]);
xlim([0, max(time)]);
ylabel({'Gaze';'Displacement';'(deg)'});
xlabel('Time (s)');
if (idx<3)
    plot_epochs(idx, [2.5, 3], 6);
else
    plot_epochs(idx, [2.5, 3], 8);
end
pos = get(gca,'position');
if not(isempty(pos_y))
    set(gca,'position',[pos(1), pos_y, pos(3), pos_h]);
end
return
end

%%



function plot_epochs(idx, yc, fs)

%epoch_rest = idx_to_epoch_rest(idx); % unused here
epoch_pf = idx_to_epoch_pf(idx);
epoch_uf = idx_to_epoch_uf(idx);

col_pf = [0.2, 0.2, 0.2];
col_uf = [0.6, 0.6, 0.6];

hold on;
plot_epochs_aux(epoch_pf, col_pf, yc, '7.5 Hz (PF)', fs);
plot_epochs_aux(epoch_uf, col_uf, yc, '30 Hz (UF)', fs);

return
end

function plot_epochs_aux(epoch, fcol, yc, t, fs)
ymin = yc(1);
ymax = yc(2);
[howmany,~] = size(epoch);
for i=1:howmany
    
    pos = epoch(i, :);
    rectangle('position', [pos(1), ymin, pos(2)-pos(1), ymax-ymin], 'edgecolor', 'none', 'facecolor', fcol);
    
    text((pos(1)+pos(2))/2, (ymin+ymax)/2, t, 'color', [1,1,1], ...
        'HorizontalAlignment', 'center', 'verticalalignment', 'middle', ...
        'fontsize', fs); 
    
end
return
end

function epoch_rest = idx_to_epoch_rest(idx)
%
global epoch_fmri_1_rest epoch_fmri_2_rest 
global epoch_fmrs_1_rest epoch_fmrs_2_rest 
%
switch idx
    case 1
        epoch_rest = epoch_fmri_1_rest;
    case 2
        epoch_rest = epoch_fmri_2_rest;  
    case 3
        epoch_rest = epoch_fmrs_1_rest;     
    case 4
        epoch_rest = epoch_fmrs_2_rest;
    otherwise
        error('undefined experiment index');
end
return
end

function epoch_pf = idx_to_epoch_pf(idx)
%
global epoch_fmri_1_pf epoch_fmri_2_pf 
global epoch_fmrs_1_pf epoch_fmrs_2_pf 
%
switch idx
    case 1
        epoch_pf   = epoch_fmri_1_pf;
    case 2
        epoch_pf   = epoch_fmri_2_pf;   
    case 3
        epoch_pf   = epoch_fmrs_1_pf;   
    case 4
        epoch_pf   = epoch_fmrs_2_pf;
    otherwise
        error('undefined experiment index');
end
return
end

function epoch_uf = idx_to_epoch_uf(idx)
%
global epoch_fmri_1_uf epoch_fmri_2_uf
global epoch_fmrs_1_uf epoch_fmrs_2_uf
%
switch idx
    case 1
        epoch_uf   = epoch_fmri_1_uf;
    case 2
        epoch_uf   = epoch_fmri_2_uf;     
    case 3
        epoch_uf   = epoch_fmrs_1_uf;     
    case 4
        epoch_uf   = epoch_fmrs_2_uf; 
    otherwise
        error('undefined experiment index');
end
return
end

