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'));
    addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
pn.figures = fullfile(rootpath, 'figures');

%% load erp_bl

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.scene_category)
        if ind_id == 1
             erpgroup.scene_category.(conds.scene_category{ind_option}) = erp_bl.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_bl.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;

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

%% plot the ERPs (for initial inspection)

% manmade = erpgroup.scene_category.manmade;
% manmade.avg = squeeze(nanmean(manmade.avg,1));
% manmade.dimord = 'chan_time';
% 
% natural = erpgroup.scene_category.natural;
% natural.avg = squeeze(nanmean(natural.avg,1));
% natural.dimord = 'chan_time';
% 
% cfg = [];
% cfg.layout = 'EEG1010.lay';
% cfg.interactive = 'yes';
% cfg.showoutline = 'yes';
% cfg.xlim = [-.2 .3];
% ft_multiplotER(cfg, natural, manmade)

%% plot topography of visual 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 = 'EEG1005.lay';
cfg.parameter = 'powspctrm';
cfg.comment = 'no';
cfg.colormap = cBrew;
cfg.colorbar = 'EastOutside';

plotData = [];
plotData.label = elec.label(1:64); % {1 x N}
plotData.dimord = 'chan';
plotData.powspctrm = squeeze(nanmean(nanmean(nanmin(mergeddata(:,:,time>0.04 & time <0.12,1:2),[],3),1),4))'; 
[~, sortidx] = sort(plotData.powspctrm, 'ascend');
idx_chans = sortidx(1);

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

figureName = ['b_topo_minimum'];
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));
condAvg1 = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1),2),4));
condAvg2 = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,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);
% ax = gca; ax.YDir = 'reverse';
xlabel('Time (s) from stim onset')
xlim([-.025 .16]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',12)

%% ERP components for subjects 1-17 and 18-33 are shifted in time!

h = figure('units','centimeters','position',[0 0 10 8]);
cla; hold on;
curData = squeeze(nanmean(mergeddata(1:17,idx_chans,:,1),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(mergeddata(18:end,idx_chans,:,1),2));
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({'Initial 17', 'Final 16'}, 'location', 'NorthWest'); legend('boxoff')
xlabel('Time (s) from stim onset')
xlim([-.025 .16]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (s)'});    
set(findall(gcf,'-property','FontSize'),'FontSize',12)

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

%% align individual subjects N1 to first negative peak

% find individual minimum (avg. across conditions) between 40 and 120 ms

time2search = find(time>0.04 & time <0.12);

newtime = 0-100*(time(2)-time(1)):(time(2)-time(1)):0+100*(time(2)-time(1));

for indID = 1:size(condAvg,1)
    tmp = find(islocalmin(condAvg(indID,time2search), ...
        'FlatSelection', 'center', ...
        'MinSeparation', 25));
    minVal1(indID) = condAvg1(indID,time2search(tmp(1)));
    curmin = time2search(tmp(1));
    alignedN1_1(indID,:) = condAvg1(indID,curmin-100:curmin+100);

    minVal2(indID) = condAvg2(indID,time2search(tmp(1)));
    curmin = time2search(tmp(1));
    alignedN1_2(indID,:) = condAvg2(indID,curmin-100:curmin+100);
    
    alignedTopo(indID,:,:) = squeeze(nanmean(mergeddata(indID,:,curmin-100:curmin+100,1:2),4));
end

% alternatively: consider global minimum
% for indID = 1:size(condAvg,1)
%     [~, minLoc1(indID)] = min(condAvg(indID,time2search));
%     curmin = time2search(minLoc1(indID));
%     minVal1(indID) = condAvg1(indID,curmin);
%     alignedN1_1(indID,:) = condAvg1(indID,curmin-100:curmin+100);
%     [~, minLoc2(indID)] = min(condAvg(indID,time2search));
%     curmin = time2search(minLoc2(indID));
%     minVal1(indID) = condAvg2(indID,curmin);
%     alignedN1_2(indID,:) = condAvg2(indID,curmin-100:curmin+100);
% end

mergeddata_aligned = cat(3, alignedN1_1, alignedN1_2);

[h, p, ci, stats] = ttest(minVal1, minVal2)

% avg across channels and conditions
condAvg_al = squeeze(nanmean(mergeddata_aligned(:,:,1:2),3));

% check that troughs are aligned
% figure; imagesc(zscore(condAvg_al,[],2))
% figure; imagesc(condAvg_al)

h = figure('units','centimeters','position',[0 0 10 8]);
cla; hold on;
% new value = old value ? subject average + grand average
curData = squeeze(mergeddata_aligned(:,:,1));
curData = curData-condAvg_al+repmat(nanmean(condAvg_al,1),size(condAvg_al,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l1 = shadedErrorBar(newtime*1000,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
curData = squeeze(mergeddata_aligned(:,:,2));
curData = curData-condAvg_al+repmat(nanmean(condAvg_al,1),size(condAvg_al,1),1);
standError = nanstd(curData,1)./sqrt(size(curData,1));
l2 = shadedErrorBar(newtime*1000,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
% ax = gca; ax.YDir = 'reverse';
legend({'manmade', 'natural'}, 'location', 'NorthWest'); legend('boxoff')
xlabel('Time (s) from local minimum')
xlim([-100 100]); %ylim([-.03 .18])
ylabel({'ERP';'(microVolts)'});
xlabel({'Time (ms) from local minimum'});    
set(findall(gcf,'-property','FontSize'),'FontSize',12)

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

%% plot topography around detected trough

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

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

plotData = [];
plotData.label = elec.label(1:64); % {1 x N}
plotData.dimord = 'chan';
plotData.powspctrm = squeeze(nanmean(nanmean(alignedTopo(:, :, newtime>-.01 & newtime < .01),3),1))'; 
[~, sortidx] = sort(plotData.powspctrm, 'ascend');
idx_chans = sortidx(1);
idx_chans_visual = idx_chans;

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');
set(findall(gcf,'-property','FontSize'),'FontSize',12)

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

%% plot raincloudplot of extracted N1 peak amplitudes

% Note that the stats reported here may vary from above in the case of
% outliers.

conds = {'manmade'; 'natural'};
uData{1} = [minVal1; minVal2];

cBrew(1,:) = [.6 .6 .6];

idx_outlier = cell(1); idx_standard = cell(1);
for indGroup = 1
    dataToPlot = uData{indGroup}';
    % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
    X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
    IndividualSlopes = uData{indGroup}(2,:)-uData{indGroup}(1,:);
    outliers = isoutlier(IndividualSlopes, 'median');
    idx_outlier{indGroup} = find(outliers);
    idx_standard{indGroup} = find(outliers==0);
end

h = figure('units','centimeter','position',[0 0 10 8]);
for indGroup = 1
    dataToPlot = uData{indGroup}';
    % read into cell array of the appropriate dimensions
    data = []; data_ws = [];
    for i = 1:2
        for j = 1:1
            data{i, j} = dataToPlot(:,i);
            % individually demean for within-subject visualization
            data_ws{i, j} = dataToPlot(:,i)-...
                nanmean(dataToPlot(:,:),2)+...
                repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
            data_nooutlier{i, j} = data{i, j};
            data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
            data_ws_nooutlier{i, j} = data_ws{i, j};
            data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
            % sort outliers to back in original data for improved plot overlap
            data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
        end
    end

    % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!

    set(gcf,'renderer','Painters')
        cla;
        cl = cBrew(indGroup,:);
        [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
        h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
        view([90 -90]);
        axis ij
    box(gca,'off')
    set(gca, 'YTickLabels', {'natural'; 'manmade'});
    yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);

    minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
    xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
    ylabel('scene category'); xlabel({'ERP'; '[Individually centered]'})

    % test linear effect
    curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
    X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
    [~, p, ci, stats] = ttest(IndividualSlopes);
    title(['M:', num2str(round(mean(IndividualSlopes),6)), '; p=', num2str(round(p,6))])
end
set(findall(gcf,'-property','FontSize'),'FontSize',12)

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


%% plot topography for different timewindows
% 
% timewins = [0.03 0.075; ...
%     0.075, 0.1; ...
%     0.1, 0.15; ...
%     0.495, 0.535];
% 
% % 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')
% 
% for indTime = 1:size(timewins,1)
%     subplot(3,2,indTime)
% 
%     cfg = [];
%     cfg.layout = 'EEG1010.lay';
%     cfg.parameter = 'powspctrm';
%     cfg.comment = 'no';
%     cfg.colormap = cBrew;
%     cfg.colorbar = 'EastOutside';
% 
%     plotData = [];
%     plotData.label = elec.label(1:64); % {1 x N}
%     plotData.dimord = 'chan';
%     plotData.powspctrm = squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>=timewins(indTime,1) & time <=timewins(indTime,2),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 = [-10 10]*10^-4; 
%     cfg.figure = h;
%     ft_topoplotER(cfg,plotData);
%     cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
% 
% end