clear all; cla; clc;

currentFile = mfilename('fullpath');
[pathstr,~,~] = fileparts(currentFile);
cd(fullfile(pathstr,'..'))
rootpath = pwd;

pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
pn.data_erp = fullfile(rootpath, 'data', 'erp');
pn.data_erf = fullfile(rootpath, 'data', 'erf');
pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
pn.tools = fullfile(rootpath, 'tools');
    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
    addpath(fullfile(pn.tools, 'BrewerMap'));
    addpath(fullfile(pn.tools, 'shadedErrorBar'));
pn.figures = fullfile(rootpath, 'figures');

%% add seed for reproducibility

rng(0, 'twister');

%% load erp

for ind_id = 1:33
    id = sprintf('sub-%03d', ind_id); disp(id)
    load(fullfile(pn.data_erf, [id,'_erf.mat']));
    for ind_option = 1:numel(conds.old)
        if ind_id == 1
             erpgroup.old.(conds.old{ind_option}) = erf.old{ind_option};
             erpgroup.old.(conds.old{ind_option}) = ...
                 rmfield(erpgroup.old.(conds.old{ind_option}), {'powspctrm'});
             erpgroup.old.(conds.old{ind_option}).dimord = 'sub_chan_freq_time';
             freq = erpgroup.old.(conds.old{ind_option}).freq;
        end
        idx_freq_t = freq >2 & freq <8;
        thetagroup.old.(conds.old{ind_option}).avg(ind_id,:,:,:) = squeeze(nanmean(erf.old{ind_option}.powspctrm(:,idx_freq_t,:),2));
        idx_freq_a = freq >7 & freq <13;
        alphagroup.old.(conds.old{ind_option}).avg(ind_id,:,:,:) = squeeze(nanmean(erf.old{ind_option}.powspctrm(:,idx_freq_a,:),2));
    end
end

time = erpgroup.old.old.time;
channels = erpgroup.old.old.label;

% identify frontal channels
%idx_chans_t = startsWith(channels, 'F') | startsWith(channels, 'A') | startsWith(channels, 'C') & ~startsWith(channels, 'CP');

% identify parieto-occipital channels
idx_chans_a = startsWith(channels, 'P') | startsWith(channels, 'O');

%% visualize for frontal theta

mergeddata = cat(4, thetagroup.old.old.avg, thetagroup.old.new.avg);

%figure; imagesc(time, [],squeeze(nanmean(thetagroup.old.new.avg-thetagroup.old.old.avg,1)))

idx_chans = [33,4,38]; 
condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));

h = figure('units','centimeters','position',[0 0 10 8]);
cla; hold on;
% highlight relevant phase in background
patches.timeVec = [0.3 0.5];
patches.colorVec = [1 .95 .8];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-4.95 -4.8];
    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
    p.EdgeColor = 'none';
end
% new value = old value ? subject average + grand average
curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
xlabel('Time (s) from stim onset')
legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
    'location', 'southeast'); legend('boxoff')
xlim([-.5 1.5]); ylim(YLim)
%xlim([-.05 .3]);
ylabel({'theta power';'(log10)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',12)

figureName = ['c_novelty_frontaltheta'];
saveas(h, fullfile(pn.figures, figureName), 'epsc');
saveas(h, fullfile(pn.figures, figureName), 'png');

%% visualize for two groups

idx_chans = [33,4,38]; 
condAvg = squeeze(nanmean(nanmean(mergeddata(1:17,idx_chans,:,1:2),2),4));

h = figure('units','centimeters','position',[0 0 20 8]);
subplot(1,2,1); cla; hold on;
% highlight relevant phase in background
patches.timeVec = [0.3 0.5];
patches.colorVec = [1 .95 .8];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-4.95 -4.75];
    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
    p.EdgeColor = 'none';
end
% new value = old value ? subject average + grand average
curData = squeeze(nanmean(mergeddata(1:17,idx_chans,:,1),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
curData = squeeze(nanmean(mergeddata(1:17,idx_chans,:,2),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
xlabel('Time (s) from stim onset')
legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
    'location', 'NorthEast'); legend('boxoff')
xlim([-.5 1.5]); ylim(YLim)
%xlim([-.05 .3]);
ylabel({'theta power';'(log10)'});
xlabel({'Time (s)'});
title('Initial 17');
set(findall(gcf,'-property','FontSize'),'FontSize',12)

subplot(1,2,2); cla; hold on;
condAvg = squeeze(nanmean(nanmean(mergeddata(18:end,idx_chans,:,1:2),2),4));

% highlight relevant phase in background
patches.timeVec = [0.3 0.5];
patches.colorVec = [1 .95 .8];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-4.95 -4.75];
    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
    p.EdgeColor = 'none';
end
% new value = old value ? subject average + grand average
curData = squeeze(nanmean(mergeddata(18:end,idx_chans,:,1),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
curData = squeeze(nanmean(mergeddata(18:end,idx_chans,:,2),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
xlabel('Time (s) from stim onset')
legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
    'location', 'NorthEast'); legend('boxoff')
xlim([-.5 1.5]); ylim(YLim)
%xlim([-.05 .3]);
ylabel({'theta power';'(log10)'});
xlabel({'Time (s)'});
title('Final 16');
set(findall(gcf,'-property','FontSize'),'FontSize',12)

figureName = ['c_novelty_frontaltheta_bygroup'];
saveas(h, fullfile(pn.figures, figureName), 'epsc');
saveas(h, fullfile(pn.figures, figureName), 'png');

%% visualize for posterior alpha

mergeddata = cat(4, alphagroup.old.old.avg, alphagroup.old.new.avg);

idx_chans = idx_chans_a;
condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));

h = figure('units','centimeters','position',[0 0 10 8]);
cla; hold on;
% highlight relevant phase in background
patches.timeVec = [0.3 0.5];
patches.colorVec = [1 .95 .8];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-5.2 -4.5];
    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
    p.EdgeColor = 'none';
end
% new value = old value ? subject average + grand average
curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'k','linewidth', 3}, 'patchSaturation', .1);
curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'r','linewidth', 3}, 'patchSaturation', .1);
xlabel('Time (s) from stim onset')
legend([l1.mainLine, l2.mainLine],{'old', 'new'}, 'location', 'northeast'); legend('boxoff')
xlim([-.5 1.5]); ylim(YLim)
%xlim([-.05 .3]);
ylabel({'alpha power';'(log10)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',15)

figureName = ['c_novelty_alpha'];
saveas(h, fullfile(pn.figures, figureName), 'epsc');
saveas(h, fullfile(pn.figures, figureName), 'png');

%% visualize for two groups

idx_chans = idx_chans_a;
condAvg = squeeze(nanmean(nanmean(mergeddata(1:17,idx_chans,:,1:2),2),4));

h = figure('units','centimeters','position',[0 0 20 8]);
subplot(1,2,1); cla; hold on;
% highlight relevant phase in background
patches.timeVec = [0.3 0.5];
patches.colorVec = [1 .95 .8];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-5.2 -4.5];
    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
    p.EdgeColor = 'none';
end
% new value = old value ? subject average + grand average
curData = squeeze(nanmean(mergeddata(1:17,idx_chans,:,1),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
curData = squeeze(nanmean(mergeddata(1:17,idx_chans,:,2),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
xlabel('Time (s) from stim onset')
legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
    'location', 'NorthEast'); legend('boxoff')
xlim([-.5 1.5]); ylim(YLim)
%xlim([-.05 .3]);
ylabel({'alpha power';'(log10)'});
xlabel({'Time (s)'});
title('Initial 17');
set(findall(gcf,'-property','FontSize'),'FontSize',12)

subplot(1,2,2); cla; hold on;
condAvg = squeeze(nanmean(nanmean(mergeddata(18:end,idx_chans,:,1:2),2),4));

% highlight relevant phase in background
patches.timeVec = [0.3 0.5];
patches.colorVec = [1 .95 .8];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-5.2 -4.5];
    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
    p.EdgeColor = 'none';
end
% new value = old value ? subject average + grand average
curData = squeeze(nanmean(mergeddata(18:end,idx_chans,:,1),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
curData = squeeze(nanmean(mergeddata(18:end,idx_chans,:,2),2));
curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
xlabel('Time (s) from stim onset')
legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
    'location', 'NorthEast'); legend('boxoff')
xlim([-.5 1.5]); ylim(YLim)
%xlim([-.05 .3]);
ylabel({'alpha power';'(log10)'});
xlabel({'Time (s)'});
title('Final 16');
set(findall(gcf,'-property','FontSize'),'FontSize',12)

figureName = ['c_novelty_alpha_bygroup'];
saveas(h, fullfile(pn.figures, figureName), 'epsc');
saveas(h, fullfile(pn.figures, figureName), 'png');

%% plot superimposed

mergeddata = cat(4, thetagroup.old.old.avg, thetagroup.old.new.avg);

idx_chans = [33,4,47]; 
h = figure('units','centimeters','position',[0 0 10 8]);
cla; hold on;
% new value = old value ? subject average + grand average
curData = squeeze(nanmean(nanmean(mergeddata(1:17,idx_chans,:,1:2),4),2));
standError = nanstd(curData,1)./sqrt(size(curData,1));
l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
curData = squeeze(nanmean(nanmean(mergeddata(18:end,idx_chans,:,1:2),4),2));
standError = nanstd(curData,1)./sqrt(size(curData,1));
l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
xlabel('Time (s) from stim onset')
legend([l1.mainLine, l2.mainLine],{'initial', 'final'}, ...
    'location', 'southeast'); legend('boxoff')
xlim([-.5 1.5]); ylim([-5.1 -4.75])
ylabel({'theta power';'(log10)'});
set(findall(gcf,'-property','FontSize'),'FontSize',12)

mergeddata = cat(4, alphagroup.old.old.avg, alphagroup.old.new.avg);

idx_chans = idx_chans_a;

h = figure('units','centimeters','position',[0 0 10 8]);
cla; hold on;
% new value = old value ? subject average + grand average
curData = squeeze(nanmean(nanmean(mergeddata(1:17,idx_chans,:,1:2),4),2));
standError = nanstd(curData,1)./sqrt(size(curData,1));
l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
curData = squeeze(nanmean(nanmean(mergeddata(18:end,idx_chans,:,1:2),4),2));
standError = nanstd(curData,1)./sqrt(size(curData,1));
l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
    {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
xlabel('Time (s) from stim onset')
legend([l1.mainLine, l2.mainLine],{'initial', 'final'}, ...
    'location', 'southeast'); legend('boxoff')
xlim([-.5 1.5]); ylim([-5.2 -4.5])
ylabel({'alpha power';'(log10)'});
set(findall(gcf,'-property','FontSize'),'FontSize',12)