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 event info

load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');

parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
for ind_param = 1:numel(parameter)
    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
end

%% load erp

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

time = erpgroup.scene_category.manmade.time;
elec = erpgroup.scene_category.manmade.elec;
channels = erpgroup.scene_category.manmade.label;

%idx_chans = find(ismember(channels, {'O1', 'Oz', 'O2'}));
%idx_chans = find(ismember(channels, {'PO7', 'PO8'}));
%idx_chans = find(ismember(channels, {'Pz', 'CPz', 'P1'}));

mergeddata = cat(4, erpgroup.scene_category.manmade.avg, ...
    erpgroup.scene_category.natural.avg);

%% plot topography of N1

% 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.08 & time <0.12,1:2),3),1),4))'; 
[~, sortidx] = sort(plotData.powspctrm, 'ascend');
idx_chans = sortidx(1:6);

cfg.marker = 'off'; 
cfg.highlight = 'yes';
cfg.highlightchannel = plotData.label(idx_chans);
cfg.highlightcolor = [1 0 0];
cfg.highlightsymbol = '.';
cfg.highlightsize = 18;
cfg.zlim = [-5 5]*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 N1 over negative maximum

% 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;
% 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')
xlim([-1 2]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',12)

h = figure('units','centimeters','position',[0 0 15 10]);
cla; hold on;
% 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';
xlabel('Time (s) from stim onset')
xlim([0 .3]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',12)

% across all channels
condAvg = squeeze(nanmean(nanmean(mergeddata(:,:,:,1:2),2),4));

h = figure('units','centimeters','position',[0 0 15 10]);
cla; hold on;
% new value = old value ? subject average + grand average
curData = squeeze(nanmean(mergeddata(:,:,:,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(:,:,:,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';
xlabel('Time (s) from stim onset')
xlim([0 .3]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',12)

%% 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';
plotData.powspctrm = squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>0.08 & time <0.12,2),3),1),4))'...
    -squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>0.08 & time <0.12,1),3),1),4))'; 
[~, sortidx] = sort(plotData.powspctrm, 'ascend');
idx_chans = sortidx(1:6);

cfg.marker = 'off'; 
cfg.highlight = 'yes';
cfg.highlightchannel = plotData.label(idx_chans);
cfg.highlightcolor = [1 0 0];
cfg.highlightsymbol = '.';
cfg.highlightsize = 18;
%cfg.zlim = [-5 5]*10^-4; 
cfg.figure = h;
ft_topoplotER(cfg,plotData);
cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');


%% visualize for negative pls channels

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

h = figure('units','centimeters','position',[0 0 10 8]);
cla; hold on;
% 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')
xlim([-1 2]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',12)

%% visualize for positive pls channels

idx_chans = [21:23, 58:62];
condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));

h = figure('units','centimeters','position',[0 0 10 8]);
cla; hold on;
% 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')
xlim([-1 2]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',12)