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'));
    
%% load erp

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

time = erpgroup.behavior.hit.time;
elec = erpgroup.behavior.hit.elec;
channels = erpgroup.behavior.hit.label;

mergeddata = cat(4, erpgroup.behavior.hit.avg, ...
    erpgroup.behavior.miss.avg);

smoothdur = 10; % 5 = 10 ms

%% get max. channels

load(fullfile(pn.data, 'd1_taskpls_erp.mat'),...
    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

stat.mask = stat.mask;
stat.prob = stat.prob.*-1;
result.vsc = result.vsc.*-1;
result.usc = result.usc.*-1;

plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,stat.time>.35 & stat.time<.55).*...
    stat.prob(:,:,stat.time>.35 & stat.time<.55),3)); 
[~, topo1] = sort(plotData.powspctrm, 'descend');

plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,stat.time>.6 & stat.time<.8).*...
    stat.prob(:,:,stat.time>0.6 & stat.time<.8),3)); 
[~, topo2] = sort(plotData.powspctrm, 'descend');

plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,stat.time>1.2 & stat.time<1.7).*...
    stat.prob(:,:,stat.time>1.2 & stat.time<1.7),3)); 
[~, topo3] = sort(plotData.powspctrm, 'descend');

%% visualize intial changes

idx_chans = topo1(1);
% 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.35 0.55];
patches.colorVec = [.8 .95 1];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-7 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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    '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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
% ax = gca; ax.YDir = 'reverse';
legend([l1.mainLine, l2.mainLine],{'hit', 'miss'}, ...
    'location', 'southwest'); legend('boxoff')
xlabel('Time (s) from stim onset')
xlim([-.25 1.9]); ylim(YLim)
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',14)

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

idx_chans = topo1(end-1:end);
% 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.35 0.55];
patches.colorVec = [.8 .95 1];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-3 3]*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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    '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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
% ax = gca; ax.YDir = 'reverse';
legend([l1.mainLine, l2.mainLine],{'hit', 'miss'}, ...
    'location', 'southwest'); legend('boxoff')
xlabel('Time (s) from stim onset')
xlim([-.25 1.9]); ylim(YLim)
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',14)

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

%% visualize second topo

idx_chans = topo2(1);
% 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.6 1.0];
patches.colorVec = [1 .95 .8];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-5 5]*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

% % highlight relevant phase in background
% patches.timeVec = [1.2 1.7];
% patches.colorVec = [1 .8 .7];
% for indP = 1:size(patches.timeVec,2)-1
%     YLim = [-5 5]*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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    '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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
% ax = gca; ax.YDir = 'reverse';
legend([l1.mainLine, l2.mainLine],{'hit', 'miss'}, ...
    'location', 'southwest'); legend('boxoff')
xlabel('Time (s) from stim onset')
xlim([-.25 1.9]); ylim(YLim)
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',14)

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

%% topo 3: negative

idx_chans = [63, 26, 21, 58];
% 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.6 1.0];
patches.colorVec = [1 .95 .8];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-2 12]*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

% highlight relevant phase in background
patches.timeVec = [1.2 1.7];
patches.colorVec = [1 .8 .7];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-2 12]*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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    '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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
% ax = gca; ax.YDir = 'reverse';
legend([l1.mainLine, l2.mainLine],{'hit', 'miss'}, ...
    'location', 'southwest'); legend('boxoff')
xlabel('Time (s) from stim onset')
xlim([-.25 1.9]); ylim(YLim)
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',14)

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

%% topo 3: positive lateral

idx_chans = [16, 53, 15, 54, 52, 51];
% 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.6 1.0];
% patches.colorVec = [1 .95 .8];
% for indP = 1:size(patches.timeVec,2)-1
%     YLim = [-4 2.5]*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

% highlight relevant phase in background
patches.timeVec = [1.2 1.7];
patches.colorVec = [1 .8 .7];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-4 2.5]*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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    '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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
% ax = gca; ax.YDir = 'reverse';
legend([l1.mainLine, l2.mainLine],{'hit', 'miss'}, ...
    'location', 'southwest'); legend('boxoff')
xlabel('Time (s) from stim onset')
xlim([-.25 1.9]); ylim(YLim)
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',14)

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

%% topo 3: positive posterior

idx_chans = [28];
% 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.6 1.0];
patches.colorVec = [1 .95 .8];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-6 2.5]*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

% highlight relevant phase in background
patches.timeVec = [1.2 1.7];
patches.colorVec = [1 .8 .7];
for indP = 1:size(patches.timeVec,2)-1
    YLim = [-6 2.5]*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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    '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,smoothts(nanmean(curData,1),'b',smoothdur),smoothts(standError,'b',smoothdur),...
    'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
% ax = gca; ax.YDir = 'reverse';
legend([l1.mainLine, l2.mainLine],{'hit', 'miss'}, ...
    'location', 'southwest'); legend('boxoff')
xlabel('Time (s) from stim onset')
xlim([-.25 1.9]); ylim(YLim)
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',14)

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