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.tools = fullfile(rootpath, 'tools');
    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
    addpath(fullfile(pn.tools, 'BrewerMap'));
    addpath(fullfile(pn.tools, 'shadedErrorBar'));
pn.figures = fullfile(rootpath, 'figures');

%% load erp

for ind_id = 1:33
    id = sprintf('sub-%03d', ind_id);
    load(fullfile(pn.data_erp, [id,'_erp_bl.mat']));
    for ind_option = 1:numel(conds.old)
        if ind_id == 1
             erpgroup.old.(conds.old{ind_option}) = erp_bl.old{ind_option};
             erpgroup.old.(conds.old{ind_option}) = ...
                 rmfield(erpgroup.old.(conds.old{ind_option}), {'avg', 'var', 'dof'});
             erpgroup.old.(conds.old{ind_option}).dimord = 'sub_chan_time';
        end
        erpgroup.old.(conds.old{ind_option}).avg(ind_id,:,:) = erp_bl.old{ind_option}.avg;
    end
end

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

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

%% plot ERP topography
% 
% % set custom colormap
% cBrew = brewermap(500,'RdBu');
% cBrew = flipud(cBrew);
% colormap(cBrew)
% 
% h = figure('units','centimeters','position',[0 0 10 10]);
% set(gcf,'renderer','Painters')
% 
% cfg = [];
% cfg.layout = 'biosemi64.lay';
% cfg.parameter = 'powspctrm';
% cfg.comment = 'no';
% cfg.colormap = cBrew;
% cfg.colorbar = 'EastOutside';
% 
% plotData = [];
% plotData.label = elec.label; % {1 x N}
% plotData.dimord = 'chan';
% plotData.powspctrm = squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>=0.3& time <=0.5,1:2),3),1),4))'; 
% 
% cfg.zlim = [-14 14]*10^-4; 
% cfg.figure = h;
% ft_topoplotER(cfg,plotData);
% cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
% % figureName = ['xxx'];
% % saveas(h, fullfile(pn.figures, figureName), 'epsc');
% % saveas(h, fullfile(pn.figures, figureName), 'png');

%% visualize frontal ERP

idx_chans = [38];

% avg across channels and conditions
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 = [-8 2]*10^-4;
    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);
% ax = gca; ax.YDir = 'reverse';
legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
    'location', 'southwest'); legend('boxoff')
xlabel('Time (s) from stim onset')
xlim([-.05 .6]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',14)

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

%% visualize central ERP

idx_chans = [32, 19, 56];

% avg across channels and conditions
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 = [-8 4]*10^-4;
    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);
% ax = gca; ax.YDir = 'reverse';
legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
    'location', 'southwest'); legend('boxoff')
xlabel('Time (s) from stim onset')
xlim([-.05 .6]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',14)

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

%% plot by group

idx_chans = [4, 38];

% avg across channels and conditions
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 = [-8 2]*10^-4;
        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);
    % ax = gca; ax.YDir = 'reverse';
    legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
        'location', 'southwest'); legend('boxoff')
    xlabel('Time (s) from stim onset')
    xlim([-.05 .6]); %ylim([-.03 .18])
    ylabel({'ERP';'(microVolts)'});
    xlabel({'Time (s)'});
    title('Initial 17')
    set(findall(gcf,'-property','FontSize'),'FontSize',14)
subplot(1,2,2); cla; hold on;
    % avg across channels and conditions
    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 = [-8 2]*10^-4;
        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);
    % ax = gca; ax.YDir = 'reverse';
    legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
        'location', 'southwest'); legend('boxoff')
    xlabel('Time (s) from stim onset')
    xlim([-.05 .6]); %ylim([-.03 .18])
    ylabel({'ERP';'(microVolts)'});
    xlabel({'Time (s)'});
    title('Final 16')
    set(findall(gcf,'-property','FontSize'),'FontSize',14)

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

%% plot difference

idx_chans = [38];

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 = [-10 5]*10^-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(mergeddata(:,idx_chans,:,2)-mergeddata(:,idx_chans,:,1));
standError = nanstd(curData,1)./sqrt(size(curData,1));
l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
legend([l1.mainLine],{'new-old'}, ...
    'location', 'southwest'); legend('boxoff')
xlabel('Time (s) from stim onset')
xlim([-.05 .6]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',14)


%% plot difference between conditions

h = figure('units','centimeters','position',[0 0 10 10]);
set(gcf,'renderer','Painters')

cfg = [];
cfg.layout = 'biosemi64.lay';
cfg.parameter = 'powspctrm';
cfg.comment = 'no';
cfg.colormap = cBrew;
cfg.colorbar = 'EastOutside';

plotData = [];
plotData.label = elec.label; % {1 x N}
plotData.dimord = 'chan';
oldminnew = mergeddata(:,:,:,1)-mergeddata(:,:,:,2);
plotData.powspctrm = squeeze(nanmean(nanmean(oldminnew(:,:,time>=.3 & time <=.5),3),1))'; 

cfg.figure = h;
ft_topoplotER(cfg,plotData);
cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');