%%%%%%%%%%%%%%
%This script starts from FieldTrip processed .mat files
%This assumes artifacts have been rejected
%%%%%%%%%%%%%%

%% Calculate ERPs for all subjects
subjectList={'1',
'2',
'3',
'4',
'5',
'6',
'7',
'8',
'9',
'10',
'11',
'12',
'13',
'14',
'15',
'16',
'17',
'18',
'19',
'20',
'21',
'22',
'23',
'24',
'25',
'26',
'27',
'28',
'29',
'30',
'31',
'32',
'33',
'34',
'35',
'36'}


%run a loop
openPath                                = '/Users/jsprouse/Desktop/'
savePath                                = '/Users/jsprouse/Desktop/'

for i = 1:length(subjectList)
	%clear whatever happened in the last loop, keep the subjectList and the variable i

	clearvars -except subjectList i openPath savePath
	
	%choose a subject from the list
	subject                             = char(subjectList(i));
    
    %load the file
    load([openPath subject '.N400_epoch.mat'])

    %split triggers into two conditions
    congruent                           = [55 66 77 88];
    incongruent                         = [155 166 177 188];
    
    cfg                                 = [];
    cfg.trials                          = ismember(N400_clean.trialinfo, congruent);
    data_congruent                      = ft_selectdata(cfg, N400_clean)

    cfg.trials                          = ismember(N400_clean.trialinfo, incongruent);
    data_incongruent                    = ft_selectdata(cfg, N400_clean)
    
    %pre-processing steps that are specific to ERPs
    cfg                                 = [];
    cfg.demean                          = 'yes'; %have to set demean to yes to get baselining
    cfg.baselinewindow                  = [-.2 0];
    cfg.lpfilter                        = 'yes';
    cfg.lpfreq                          = 30;
    
    preprocessed_congruent              = ft_preprocessing(cfg, data_congruent)
    preprocessed_incongruent            = ft_preprocessing(cfg, data_incongruent)
    
    %create ERPs
    cfg                                 = [];
    erp_congruent                       = ft_timelockanalysis(cfg, preprocessed_congruent);
    erp_incongruent                     = ft_timelockanalysis(cfg, preprocessed_incongruent);
    
    %save the ERPs
    save([savePath subject '.N400_erp_congruent.mat'],'erp_congruent');
    save([savePath subject '.N400_erp_incongruent.mat'],'erp_incongruent');

end;

%% Calculate a grandaveraged ERP - This is the set of subjects that had trials > 20, N=33
subjectList={'alanna.cooney',
'amanda.womeldorf',
'anees.rehman',
'ashley.clancy',
'austin.ficara',
'brad.dalelo',
'chibuikem.uche',
%'corinna.kraemer',
'danni.dong',
'erik.pettit',
'gilliam weintraub',
%'hao.xu',
'jason.sittambalam',
'Jayvon Snider',
'john.fable',
'john.han',
'john.martin',
'kenil.pobari',
'kimberly.clemens2',
'kristen.fagan',
'madeline sadlowski',
'martin.marx',
'meghan.nalaboff',
'nicholas.pham',
%'nicolas.engelbert',
'nicole.miller',
'noah.fields',
'rebecca.feuerberg',
'sharon.park',
'shinelle.wilkins',
'talia.ward',
'thomas.binghi',
'valerie.wann',
'Vimratha Teepireddy',
'wenqi.sun',
'Zach.Zwick'}

loadPath                                = '/Users/jsprouse/Desktop/2. N400 ERPs/'
savePath                                = '/Users/jsprouse/Desktop/2. N400 ERPs/'
%create an array to save the individual subject ERPs
erparray=cell(2,length(subjectList));

%loop through the subjects to load their ERPs
for k=1:length(subjectList)
    
    %clear whatever happened in the last loop, keep the subjectList and the variable i
	clearvars -except subjectList k loadPath savePath erparray
    
    %choose a subject from the list
	subject                             = char(subjectList(k));
    
    %load the file and save the ERPs to the array
    load([loadPath subject '.N400_erp_congruent.mat']) 
    erparray{1,k}                       = erp_congruent;
    
    load([loadPath subject '.N400_erp_incongruent.mat']) 
    erparray{2,k}                       = erp_incongruent;
end

%save the erparrays. They are useful sometimes (like some ERP plots)
save ([savePath 'N400_erp_array_33.mat'], 'erparray'); 


%Calculate grandaverage ERP without keeping the individuals (so, just for
%plotting, and things like ft_math)
cfg                                     = [];  
cfg.keepindividual                      = 'no';  

erp_grandaverage_keepNO_congruent       = ft_timelockgrandaverage(cfg, erparray{1,:});

erp_grandaverage_keepNO_incongruent     = ft_timelockgrandaverage(cfg, erparray{2,:});

%Calculate grandaverage ERP keeping the individuals (for both plotting and
%stats)

cfg                                     = [];  
cfg.keepindividual                      = 'yes';  

%calculate each grandaverage ERP and save
erp_grandaverage_keepYES_congruent      = ft_timelockgrandaverage(cfg, erparray{1,:});

erp_grandaverage_keepYES_incongruent    = ft_timelockgrandaverage(cfg, erparray{2,:});

%save the grandaverages
save ([savePath 'N400_erp_grandaverage_congruent_33.mat'], 'erp_grandaverage_keepYES_congruent', 'erp_grandaverage_keepNO_congruent'); 

save ([savePath 'N400_erp_grandaverage_incongruent_33.mat'], 'erp_grandaverage_keepYES_incongruent', 'erp_grandaverage_keepNO_incongruent'); 

%% Calculate a grandaveraged ERP - trials>=20, N=29
subjectList={'alanna.cooney',
'amanda.womeldorf',
'anees.rehman',
'ashley.clancy',
'austin.ficara',
'brad.dalelo',
'chibuikem.uche',
%'corinna.kraemer',
%'danni.dong',
'erik.pettit',
'gilliam weintraub',
%'hao.xu',
'jason.sittambalam',
'Jayvon Snider',
%'john.fable',
'john.han',
'john.martin',
%'kenil.pobari',
'kimberly.clemens2',
'kristen.fagan',
'madeline sadlowski',
'martin.marx',
'meghan.nalaboff',
'nicholas.pham',
%'nicolas.engelbert',
'nicole.miller',
'noah.fields',
'rebecca.feuerberg',
'sharon.park',
'shinelle.wilkins',
'talia.ward',
%'thomas.binghi',
'valerie.wann',
'Vimratha Teepireddy',
'wenqi.sun',
'Zach.Zwick'}

loadPath                                = '/Users/jsprouse/Desktop/2. N400 ERPs/'
savePath                                = '/Users/jsprouse/Desktop/2. N400 ERPs/'
%create an array to save the individual subject ERPs
erparray=cell(2,length(subjectList));

%loop through the subjects to load their ERPs
for k=1:length(subjectList)
    
    %clear whatever happened in the last loop, keep the subjectList and the variable i
	clearvars -except subjectList k loadPath savePath erparray
    
    %choose a subject from the list
	subject                             = char(subjectList(k));
    
    %load the file and save the ERPs to the array
    load([loadPath subject '.N400_erp_congruent.mat']) 
    erparray{1,k}                       = erp_congruent;
    
    load([loadPath subject '.N400_erp_incongruent.mat']) 
    erparray{2,k}                       = erp_incongruent;
end

%save the erparrays. They are useful sometimes (like some ERP plots)
save ([savePath 'N400_erp_array_29.mat'], 'erparray'); 


%Calculate grandaverage ERP without keeping the individuals (so, just for
%plotting, and things like ft_math)
cfg                                     = [];  
cfg.keepindividual                      = 'no';  

erp_grandaverage_keepNO_congruent       = ft_timelockgrandaverage(cfg, erparray{1,:});

erp_grandaverage_keepNO_incongruent     = ft_timelockgrandaverage(cfg, erparray{2,:});

%Calculate grandaverage ERP keeping the individuals (for both plotting and
%stats)

cfg                                     = [];  
cfg.keepindividual                      = 'yes';  

%calculate each grandaverage ERP and save
erp_grandaverage_keepYES_congruent      = ft_timelockgrandaverage(cfg, erparray{1,:});

erp_grandaverage_keepYES_incongruent    = ft_timelockgrandaverage(cfg, erparray{2,:});

%save the grandaverages
save ([savePath 'N400_erp_grandaverage_congruent_29.mat'], 'erp_grandaverage_keepYES_congruent', 'erp_grandaverage_keepNO_congruent'); 

save ([savePath 'N400_erp_grandaverage_incongruent_29.mat'], 'erp_grandaverage_keepYES_incongruent', 'erp_grandaverage_keepNO_incongruent'); 

%% load data for plotting and statistics
loadPath                                = '/Users/jsprouse/Desktop/'

load([loadPath 'N400_erp_grandaverage_congruent_29.mat'])
load([loadPath 'N400_erp_grandaverage_incongruent_29.mat'])
load([loadPath 'N400_erp_array_29.mat'])

%% plot ERP, multiplot 
cfg                                     = [];
cfg.showlabels                          = 'yes'; 
cfg.fontsize                            = 14; 
cfg.layout                              = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20180811/template/layout/acticap-64ch-standard2.mat';
lay                                     = ft_prepare_layout(cfg);
cfg.xlim                                = [-.2 1];
%cfg.baseline                            = [.3 .5];
cfg.interactive                         = 'yes';
cfg.showoutline                         = 'yes';

figure;
legend('congruent','incongruent')
ft_multiplotER(cfg, erp_grandaverage_keepNO_congruent, erp_grandaverage_keepNO_incongruent);

%% plot ERP, one channel, all subjects
%this requires an array that holds all of the subjects (erparray)

time = [.3 .5];
ymax = 10;
ymin = -10;
channel = 24;

figure; 
for isub = 1:size(erparray,2)
    subplot(8,4,isub)
    % use the rectangle to indicate a time range of specific interest
    % this is for one channel, given by the number after the .avg
    rectangle('Position',[time(1) ymin (time(2)-time(1)) abs(ymax-ymin)],'FaceColor',[0.7 0.7 0.7]);
    hold on;
    % plot the lines in front of the rectangle
    plot(erparray{1,isub}.time,erparray{1,isub}.avg(channel,:));
    plot(erparray{2,isub}.time,erparray{2,isub}.avg(channel,:),'r');
    title(strcat('subject ',num2str(isub)))
    ylim([ymin ymax])
    xlim([-.2 1])
end 
subplot(8,4,32); 
text(0.5,0.5,'congruent','color','b') ;text(0.5,0.3,'incongruent','color','r')
axis off

%% plot difference between condition for each subject during a time window of interest
%this uses erparray


time = [.3 .5];
ymax = 10;
ymin = -10;
channel = 24;

% find the time points for the effect of interest using the first subject
% of erparray
timesel_congruent = find(erparray{1,1}.time >= time(1) & erparray{1,1}.time <= time(2));
timesel_incongruent = find(erparray{2,1}.time >= time(1) & erparray{2,1}.time <= time(2));
 
% select the individual subject data from the time points and calculate the mean
for isub = 1:size(erparray,2)
    values_congruent(isub)  = mean(erparray{1,isub}.avg(channel,timesel_congruent));   
    values_incongruent(isub)  = mean(erparray{2,isub}.avg(channel,timesel_incongruent));
end
 

% plot to see the effect in each subject
M = [values_congruent',values_incongruent'];
figure; plot(M','o-'); xlim([0.5 2.5])
legend('location','EastOutside');

%% plot a topoplot using topoplotER?

%calculate difference wave
cfg                                     = [];
cfg.operation                           = 'subtract';
cfg.parameter                           = 'avg';
erp_difference                          = ft_math(cfg, erp_grandaverage_keepNO_incongruent, erp_grandaverage_keepNO_congruent);

%plot
cfg                                     = [];
cfg.xlim = [0:0.05:.75]; %this plots a series of averages over time, where the middle number is the window size
cfg.zlim = [-5 5];            
%cfg.baseline = [-0.2 0];          
%cfg.baselinetype = 'db';
cfg.layout                              = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20180811/template/layout/acticap-64ch-standard2.mat';
cfg.comment = 'xlim';
cfg.commentpos = 'title';
figure; ft_topoplotER(cfg, erp_difference); 

ft_hastoolbox('brewermap', 1);         % ensure this toolbox is on the path
colormap(flipud(brewermap(64,'RdBu'))) % change the colormap

%% statistics for ERPs - cluster based permutation test

%rebaseline the data?
%cfg = [];
%cfg.baseline                            = [.3 .5]; 
%cfg.channel                             = {'all'}

%erp_grandaverage_keepYES_congruent    = ft_timelockbaseline(cfg, erp_grandaverage_keepYES_congruent)
%erp_grandaverage_keepYES_incongruent  = ft_timelockbaseline(cfg, erp_grandaverage_keepYES_incongruent)


%define neighbors
cfg = [];
cfg.channel                             = {'all' '-M1' '-M2'};
cfg_neighb.method                       = 'distance'; 
cfg_neighb.neighbourdist               = 0.125;
cfg_neighb.layout                       = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20180811/template/layout/acticap-64ch-standard2.mat';
cfg_neighb.channel                      = {'all' '-M1' '-M2'};
cfg_neighb.enableedit                   = 'no';
cfg_neighb.feedback                     = 'no';	%this shows you the channel map, and then you have to click X to close it to move on.

cfg.neighbours                          = ft_prepare_neighbours(cfg_neighb);

%define design matrix for the experiment
%the first is the number of subjects, which I have softcoded using the
%grandaverage of one of the conditions
subj                                    = size(erp_grandaverage_keepYES_congruent.individual,1);
design                                  = zeros(2,2*subj);

for i = 1:subj
  design(1,i) = i;
end

for i = 1:subj
  design(1,subj+i) = i;
end

design(2,1:subj)            = 1;
design(2,subj+1:2*subj)     = 2;

%configuration for the test
cfg.design                              = design;
cfg.uvar                                = 1; % subject number (unit variable) on line 1 of the design matrix
cfg.ivar                                = 2; %=2; condition number (independent variable) on line 2 of the design matrix

cfg.channel                             = {'all' '-M1' '-M2'};

cfg.method                              = 'montecarlo';
cfg.statistic                           = 'depsamplesT';
cfg.correctm                            = 'cluster';
cfg.clusteralpha                        = 0.05;
cfg.clusterstatistic                    = 'maxsum';
cfg.minnbchan                           = 2;
cfg.tail                                = 0; %two sided
cfg.clustertail                         = 0;
cfg.alpha                               = 0.025; %alpha level per tail
cfg.numrandomization                    = 500;

%tests on specific time windows
cfg.latency                             = [0 1];
full                                    = ft_timelockstatistics(cfg, erp_grandaverage_keepYES_incongruent, erp_grandaverage_keepYES_congruent)

%% plotting statistical results with ft_clusterplot?
%this plots only the significant clusters, which is nice, but it plots over
%the stats (t-values), which may be less desirable
%this works over each sample (so, 1ms), so we may need to downsample, then
%re-run the statistics to use this (or find a way to downsample the
%statistics output)
cfg                                     = [];
cfg.alpha                               = .025;            
cfg.layout                              = '/Users/jsprouse/Documents/MATLAB/eeglab14_1_2b/external/fieldtrip-20180811/template/layout/acticap-64ch-standard2.mat';
figure; ft_clusterplot(cfg, full);

%change to color brewer?
ft_hastoolbox('brewermap', 1);         % ensure this toolbox is on the path
colormap(flipud(brewermap(64,'RdBu'))) 


%% plotting the statistical results using some code from fieldtrip website
%this works from the grandaverages keepNO
%this only looks at one set of clusters (positive or negative) at a time
%may be able to edit it to do both sets of clusters... in the future.

%need to tell it the name of the test that you are going to look at
stat = subject

%rebaseline the data?
%cfg = [];
%cfg.baseline                           = [.3 .5]; 
%cfg.channel                            = {'all'}

%erp_grandaverage_keepNO_congruent      = ft_timelockbaseline(cfg, erp_grandaverage_keepNO_congruent)
%erp_grandaverage_keepNO_incongruent    = ft_timelockbaseline(cfg, erp_grandaverage_keepNO_incongruent)


%calculate difference
cfg                                     = [];
cfg.operation                           = 'subtract';
cfg.parameter                           = 'avg';
erp_difference                          = ft_math(cfg, erp_grandaverage_keepNO_incongruent, erp_grandaverage_keepNO_congruent);


figure;  
% define parameters for plotting
timestep                                = 0.025;      %(in seconds)
sampling_rate                           = 1000; %this shouldn't be hardcoded, but I don't know where to find it in the grandaverages
sample_count                            = length(stat.time);
j                                       = [0:timestep:1];   % Temporal endpoints (in seconds) of the ERP average computed in each subplot
m                                       = [1:timestep*sampling_rate:sample_count];  % temporal endpoints in EEG samples
% get relevant (significant) values
neg_cluster_pvals                       = [stat.negclusters(:).prob];

% In case you have downloaded and loaded the data, ensure stat.cfg.alpha exists:
if ~isfield(stat.cfg,'alpha'); stat.cfg.alpha = 0.025; end; 

neg_signif_clust                        = find(neg_cluster_pvals < stat.cfg.alpha);
neg                                     = ismember(stat.negclusterslabelmat, neg_signif_clust);
  
% First ensure the channels to have the same order in the average and in the statistical output.
% This might not be the case, because ft_math might shuffle the order  
[i1,i2]                                 = match_str(erp_difference.label, stat.label); 

% plot
for k = 1:20; %this is hardcoded, but should vary with the window size and timestep set above
     subplot(4,5,k);   
     cfg = [];   
     cfg.xlim=[j(k) j(k+1)];   
     cfg.zlim = [-8 8];   %this is hardcoded, but should vary with the data
     neg_int = zeros(numel(erp_difference.label),1);
     neg_int(i1) = all(neg(i2, m(k):m(k+1)), 2);
     cfg.highlight = 'on';
     cfg.highlightchannel = find(neg_int);       
     cfg.comment = 'xlim';   
     cfg.commentpos = 'title';   
     cfg.layout = '/Users/jsprouse/Documents/MATLAB/eeglab13_6_5b/external/fieldtrip-20180510/template/layout/acticap-64ch-standard2.mat';
     ft_topoplotER(cfg, erp_difference);
end  