WildLab/code/collectUnfold_plotERPs.m

%%
% This script loads all the deconvolved data and plots them as ERPs


cd '/net/store/nbp/users/agert/itw'
init_itw

projectFolder=['/net/store/nbp/projects/IntoTheWild'];

codingscheme = 'mashup';
phase='WLFO';
filtType='causal';

timewin=[-0.5 1.0];


% LOADING
% needed for chanlocs
if strcmp(phase,'Lab')
    Data=load(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/' phase '/collectedERPs_-5001000_causal_noCar.mat']);
elseif strcmp(phase,'WLFO')
    Data=load(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/' phase '/collectedERPs_-5001000_causal.mat']);
    
end


fns=fieldnames(Data);
collData = Data.(fns{1});

paramNames=[];
%

d2nd = [];

for subj = [8 9 29:36 38:41 44 45 48 51:53];
    fPath = fullfile(projectFolder,['Daten/EEG/' phase '/unfold_' phase '/' filtType '/']);
    
    fPath = fullfile(fPath,['ufresult_subj' num2str(subj) ,'_allElec_' codingscheme '_-500-1000_regularized0.mat']);
    try
        % load subject results
        if ~exist(fPath,'file')
            continue
        end
        fprintf('\n \n loading subject %i \n',subj)
        
        d = load(fPath);
        
        ufresult = d.ufresult;
        
        tmpBeta=[];
        tmpBeta_nodc=[];
        if strcmp(phase,'Lab') && subj==51 % sub 51 does not have predictore sacAmp -> adjust structure accordingly
            tmpBeta(:,:,1:9)=ufresult.beta(:,:,1:9);
            tmpBeta(:,:,10:13)=nan(length(d2nd.chanlocs),length(ufresult.times),4);
            tmpBeta(:,:,14:20)=ufresult.beta(:,:,10:16);
            tmpBeta_nodc(:,:,1:9)=ufresult.beta_nodc(:,:,1:9);
            tmpBeta_nodc(:,:,10:13)=nan(length(d2nd.chanlocs),length(ufresult.times),4);
            tmpBeta_nodc(:,:,14:20)=ufresult.beta_nodc(:,:,10:16);
            
            ufresult.beta=tmpBeta;
            ufresult.beta_nodc=tmpBeta_nodc;
            if ~isempty(d2nd)
                ufresult.param=d2nd.param;
            end
        end
        
        %         first subject: initialize d2nd (2nd-level analysis)
        if isempty(d2nd)
            d2nd = ufresult;
            d2nd.unfold(1) = d2nd.unfold;
            d2nd.subject = subj;
        else
            d2nd.unfold(end+1) = ufresult.unfold;
            d2nd.beta(:,:,:,end+1) = ufresult.beta;
            d2nd.beta_nodc(:,:,:,end+1) = ufresult.beta_nodc;
            d2nd.subject(end+1) = subj;
        end
        
        
    catch e
        display(['problem with subject ',num2str(subj), ':   ' ,e.message])
    end
    
end


%% Calculate ERPs

elecfind = @(x)find(strcmp(x,{d2nd.chanlocs.labels}));
chan = 1:length(d2nd.chanlocs);%[elecfind('PO8')];

if strcmp(codingscheme,'effects')
    [x2n,x2h,n2x,h2x,x2n_nodc,x2h_nodc,n2x_nodc,h2x_nodc,SB]=ITW_unfold_calcERPs(d2nd,codingscheme,chan,phase);
else
    [n2n,n2h,h2n,h2h,n2n_nodc,n2h_nodc,h2n_nodc,h2h_nodc,SB]=ITW_unfold_calcERPs(d2nd,codingscheme,chan,phase);
end


%% Save in ERP structure

ufERP.times=d2nd.times;
ufERP.chanlocs=d2nd.chanlocs;
ufERP.subject=d2nd.subject;
ufERP.data=[];
ufERP.data(:,:,:,1)=n2n;
ufERP.data(:,:,:,2)=n2h;
ufERP.data(:,:,:,3)=h2n;
ufERP.data(:,:,:,4)=h2h;
ufERP.cond={'n2n','n2h','h2n','h2h'};


save(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/' phase '/unfold_' phase '/' filtType '/unfoldERPs_-5001000.mat'],'ufERP')



%% PLOT ERPs
%  load(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/WLFO/unfold_WLFO/' filtType '/unfoldERPs_-5001000.mat'],'ufERP')

elecfind = @(x)find(strcmp(x,{d2nd.chanlocs.labels}));
chan =[elecfind('PO8'),elecfind('P8'),elecfind('PO7'),elecfind('P7')];

% plot 2x2 results of deconvolution
    figure('units','normalized','outerposition',[0 0 1 1])
    hold on,  plot(d2nd.times,squeeze(mean(mean(n2n(:,chan,:)),2)),'b--','LineWidth',2);
    hold on, plot(d2nd.times,squeeze(mean(mean(n2h(:,chan,:)),2)),'r--','LineWidth',2);
    hold on, plot(d2nd.times,squeeze(mean(mean(h2n(:,chan,:)),2)),'b','LineWidth',2);
    hold on, plot(d2nd.times,squeeze(mean(mean(h2h(:,chan,:)),2)),'r','LineWidth',2);
    
    legend({'BG -> BG','BG -> HF','HF -> BG','HF -> HF'},'FontSize',20)
    xlim([-0.5 1])
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', 18)
    ylim([-6 6]);
    vline(0,'k')
    hline(0,'k');
    box off
    title([phase ' Deconvolution'],'FontSize',20)
    xlabel('Time [sec]', 'FontSize',20)
    ylabel('Amplitude [microV]', 'FontSize',20);
    
%     export_fig(['/home/student/a/agert/Desktop/Into the Wild/LabvsWLFO/autoplot/2x2_ERP_Deconv_Unfold_' filtType '.eps'])
    
    % plot 2x2 design of NO deconvolution
    figure('units','normalized','outerposition',[0 0 1 1])
    hold on, plot(d2nd.times,squeeze(mean(mean(n2n_nodc(:,chan,:)),2)),'b--','LineWidth',2);
    hold on, plot(d2nd.times,squeeze(mean(mean(n2h_nodc(:,chan,:)),2)),'r--','LineWidth',2);
    hold on, plot(d2nd.times,squeeze(mean(mean(h2n_nodc(:,chan,:)),2)),'b','LineWidth',2);
    hold on, plot(d2nd.times,squeeze(mean(mean(h2h_nodc(:,chan,:)),2)),'r','LineWidth',2);
    
    legend({'BG -> BG','BG -> HF','HF -> BG','HF -> HF'},'FontSize',20)
    xlim([-0.5 1])
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', 18)
    ylim([-6 6])
    vline(0,'k')
    hline(0,'k')
    box off
    title([phase ' No Deconvolution'],'FontSize',20)
    xlabel('Time [sec]', 'FontSize',20)
    ylabel('Amplitude [microV]', 'FontSize',20);

    
    % within face fixation
    if strcmp(phase,'WLFO')
        figure('units','normalized','outerposition',[0 0 1 1])
        hold on,plot(d2nd.times,squeeze(mean(mean(n2h(:,chan,:)),2)),'b','LineWidth',2);
        hold on, plot(d2nd.times,squeeze(mean(mean(h2h(:,chan,:)),2)),'r--','LineWidth',2);
        hold on, plot(d2nd.times,squeeze(mean(mean(SB (:,chan,:)),2)),'r','LineWidth',2);
        
        
        legend({'BG -> HF','HF -> HF','Within Box'},'FontSize',20)
        xlim([-0.5 1])
        xt = get(gca, 'XTick');
        set(gca, 'FontSize', 18)
        ylim([-6 6])
        vline(0,'k')
        % vline(0.17,'--k')
        hline(0,'k')
        box off
        title([phase ' Samebox Effect'],'FontSize',20)
        xlabel('Time [sec]', 'FontSize',20)
        ylabel('Amplitude [microV]', 'FontSize',20);
        %     export_fig(['/home/student/a/agert/Desktop/Into the Wild/LabvsWLFO/autoplot/SB_ERP_' filtType '.png'])
        %     export_fig(['/home/student/a/agert/Desktop/Into the Wild/LabvsWLFO/autoplot/SB_ERP.eps'])
    end
    %% 1x2 ERP plots
    n2x=(n2h+n2n)./2;
    h2x=(h2h+h2n)./2;
    x2n=(h2n+n2n)./2;
    x2h=(h2h+n2h)./2;
    
    % plot 1x2 for type of previous fixation
    figure('units','normalized','outerposition',[0 0 1 1])
    hold on,  plot(d2nd.times,squeeze(mean(mean(n2x(:,chan,:)),2)),'k--','LineWidth',2);
    hold on, plot(d2nd.times,squeeze(mean(mean(h2x(:,chan,:)),2)),'k','LineWidth',2);
    
    
    legend({'BG -> X','HF -> X'},'FontSize',20)
    xlim([-0.5 1])
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', 18)
    ylim([-6 6])
    vline(0,'k')
    % vline(0.17,'--k')
    hline(0,'k')
    box off
    title([phase ' Main Effect Prev'],'FontSize',20)
    xlabel('Time [sec]', 'FontSize',20)
    ylabel('Amplitude [microV]', 'FontSize',20);
    %export_fig(['/home/student/a/agert/Desktop/Into the Wild/LabvsWLFO/autoplot/1x2_ERP_MainPrev_Unfold_' filtType '.eps'])
    
    % plot 1x2 for type of current fixation
    x2nCI=zeros(length(squeeze(mean(x2n(:,chan,:),2))),2); 
    x2hCI=zeros(length(squeeze(mean(x2h(:,chan,:),2))),2); 
    for timep=1:length(squeeze(mean(x2n(:,chan,:),2)))
        x2nCI(timep,:)= bootci(10000,@mean,mean(x2n(:,chan,timep),2));
        x2hCI(timep,:)= bootci(10000,@mean,mean(x2h(:,chan,timep),2));
    end
    
    figure('units','normalized','outerposition',[0 0 1 1])
%     hold on, plot(d2nd.times,squeeze(mean(mean(x2h(:,chan,:)),2)),'r','LineWidth',2);
%     
    x1 = [squeeze(d2nd.times), fliplr(squeeze(d2nd.times))];
    inBetween = [x2nCI(:,1)', fliplr(x2nCI(:,2)')];
    fill(x1, inBetween, [0,0.5,0.5]);
    hold on;
    plot(d2nd.times,squeeze(mean(mean(x2n(:,chan,:)),2)),'b','LineWidth',2);
%     hold on, plot(d2nd.times,x2nCI,'b--','LineWidth',2);
%     hold on, plot(d2nd.times,x2hCI,'r--','LineWidth',2);

    inBetween2 = [x2hCI(:,1)', fliplr(x2hCI(:,2)')];
    fill(x1, inBetween2, [1,0.5,0]);
    hold on;
    plot(d2nd.times,squeeze(mean(mean(x2h(:,chan,:)),2)),'r','LineWidth',2);

    legend({'X -> BG','X -> HF'},'FontSize',20)
    xlim([-0.5 1])
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', 18)
    ylim([-7 7])
    vline(0,'k')
    vline(0.12,'k--') %N170 start
    vline(0.2,'k--') %N170 stop
    vline(-0.3,'k--')
    % vline(0.17,'--k')
    hline(0,'k')
    
    box off
    title([phase ' Main Effect Curr'],'FontSize',20)
    xlabel('Time [sec]', 'FontSize',20)
    ylabel('Amplitude [microV]', 'FontSize',20);
    export_fig(['/home/student/a/agert/Desktop/Into the Wild/LabvsWLFO/autoplot/1x2_ERP_MainCurr_Unfold_' filtType '_' phase '_withCI.eps'])
   


%% Plot effect topos (timerange specified for N170)

effectData=squeeze(mean(x2h(:,:,:)))-squeeze(mean(x2n(:,:,:)));
tp1=find(d2nd.times>0.120, 1, 'first');
tp2=find(d2nd.times<0.200, 1, 'last');
topoData=squeeze(mean(effectData(:,tp1:tp2),2));
figure, topoplot(topoData',d2nd.chanlocs,'conv','on');
title(['N170 topo 120:200 ' phase])

%% viewing behavior
d2nd2 = subjectwise_getparam(d2nd,{{'sac_amplitude',[0.5 1 2.5 5 7.5 10]}},1);

uf_plotParam(d2nd2,'channel','Oz','plotSeparate','all','plotParam',{'fix_avgpos_y','fix_avgpos_x'})
uf_plotParam(d2nd2,'channel','Oz', 'plotSeparate','all','plotParam','sac_amplitude','deconv',1)
set(findobj(gca,'-property','LineWidth'),'LineWidth',2)

%% Component amplitudes
%P100: 80-130
%N170: 120-200


p100Start=find(d2nd.times>=0.080,1,'first');
p100End=find(d2nd.times<=0.130,1,'last');
n170Start=find(d2nd.times>=0.120,1,'first');
n170End=find(d2nd.times<=0.200,1,'last');

meanN2N.p100=squeeze(mean(ufERP.data(:,chan,p100Start:p100End,1),2));
meanN2H.p100=squeeze(mean(ufERP.data(:,chan,p100Start:p100End,2),2));
meanH2N.p100=squeeze(mean(ufERP.data(:,chan,p100Start:p100End,3),2));
meanH2H.p100=squeeze(mean(ufERP.data(:,chan,p100Start:p100End,4),2));

meanN2N.n170=squeeze(mean(ufERP.data(:,chan,n170Start:n170End,1),2));
meanN2H.n170=squeeze(mean(ufERP.data(:,chan,n170Start:n170End,2),2));
meanH2N.n170=squeeze(mean(ufERP.data(:,chan,n170Start:n170End,3),2));
meanH2H.n170=squeeze(mean(ufERP.data(:,chan,n170Start:n170End,4),2));


meanX2N.p100=squeeze(mean(mean(ufERP.data(:,chan,p100Start:p100End,[1,3]),2),4));
meanX2H.p100=squeeze(mean(mean(ufERP.data(:,chan,p100Start:p100End,[2,4]),2),4));

meanX2N.n170=squeeze(mean(mean(ufERP.data(:,chan,n170Start:n170End,[1,3]),2),4));
meanX2H.n170=squeeze(mean(mean(ufERP.data(:,chan,n170Start:n170End,[2,4]),2),4));


meanN2X.p100=squeeze(mean(mean(ufERP.data(:,chan,p100Start:p100End,[1,2]),2),4));
meanH2X.p100=squeeze(mean(mean(ufERP.data(:,chan,p100Start:p100End,[3,4]),2),4));

meanN2X.n170=squeeze(mean(mean(ufERP.data(:,chan,n170Start:n170End,[1,2]),2),4));
meanH2X.n170=squeeze(mean(mean(ufERP.data(:,chan,n170Start:n170End,[3,4]),2),4));



%%
p100=meanH2H.p100;
n170=meanH2H.n170;
[p100Amp,p100Time]=max(p100,[],2);
[n170Amp,n170Time]=min(n170,[],2);
p2pEff=p100Amp-n170Amp;
p100PeakTime=ufERP.times(p100Start+p100Time-1);
n170PeakTime=ufERP.times(n170Start+n170Time-1);
p100AmpCI = bootci(10000,@mean,p100Amp);
n170AmpCI = bootci(10000,@mean,n170Amp);
p2pAmpCI  = bootci(10000,@mean,p2pEff);
p100TimeCI = bootci(10000,@mean,p100PeakTime);
n170TimeCI = bootci(10000,@mean,n170PeakTime);

t=0;

if t==100
    fprintf(['P100 for ' phase ' H2H on average ' num2str(mean(p100Amp)) 'uV 95CI [' num2str(p100AmpCI(1)) ';' num2str(p100AmpCI(2))  '] at ' num2str(mean(p100PeakTime)) 'ms [' num2str(p100TimeCI(1)) ';' num2str(p100TimeCI(2))  '].\n'])
elseif t==170
    fprintf(['N170 for ' phase ' H2H on average ' num2str(mean(n170Amp)) 'uV 95CI [' num2str(n170AmpCI(1)) ';' num2str(n170AmpCI(2))  '] at ' num2str(mean(n170PeakTime)) 'ms [' num2str(n170TimeCI(1)) ';' num2str(n170TimeCI(2))  '].\n'])
elseif t==0
    fprintf(['P2P for ' phase ' H2H on average ' num2str(mean(p2pEff)) 'uV 95CI [' num2str(p2pAmpCI(1)) ';' num2str(p2pAmpCI(2))  '].\n'])
    
end


%% conditionwise differences

cond1=meanX2N;
cond2=meanX2H;

[cond1p100Amp,cond1p100Time]=max(cond1.p100,[],2);
[cond1n170Amp,cond1n170Time]=min(cond1.n170,[],2);
[cond2p100Amp,cond2p100Time]=max(cond2.p100,[],2);
[cond2n170Amp,cond2n170Time]=min(cond2.n170,[],2);
condDiffp100=cond2p100Amp-cond1p100Amp;
condDiffn170=cond2n170Amp-cond1n170Amp;
p100AmpCI = bootci(10000,@mean,condDiffp100);
n170AmpCI = bootci(10000,@mean,condDiffn170);

fprintf(['Diff P100 for ' phase ' on average ' num2str(mean(condDiffp100)) 'uV 95CI [' num2str(p100AmpCI(1)) ';' num2str(p100AmpCI(2))  '].\n'])

fprintf(['Diff N170 for ' phase ' on average ' num2str(mean(condDiffn170)) 'uV 95CI [' num2str(n170AmpCI(1)) ';' num2str(n170AmpCI(2))  '].\n'])


% [N2NAmp,N2NTime]=min(meanN2N,[],2);
% [N2HAmp,N2HTime]=min(meanN2H,[],2);
% [H2NAmp,H2NTime]=min(meanH2N,[],2);
% [H2HAmp,H2HTime]=min(meanH2H,[],2);
%
%
% [X2HAmp,X2HTime]=max(meanX2H,[],2);
% [X2NAmp,X2NTime]=max(meanX2N,[],2);
%
% X2NPeakTime=ufERP.times(n170Start+X2NTime);
% X2HPeakTime=ufERP.times(n170Start+X2HTime);

% X2NAmpCI = bootci(10000,@mean,X2NAmp);
% X2HAmpCI = bootci(10000,@mean,X2HAmp);
% X2NTimeCI = bootci(10000,@mean,X2NPeakTime);
% X2HTimeCI = bootci(10000,@mean,X2HPeakTime);

% fprintf(['N170 for ' phase ' N2N on average ' num2str(mean(N2NRange)) 'uV.\n'])
% fprintf(['N170 for ' phase ' N2H on average ' num2str(mean(N2HRange)) 'uV.\n'])
% fprintf(['N170 for ' phase ' H2N on average ' num2str(mean(H2NRange)) 'uV.\n'])
% fprintf(['N170 for ' phase ' H2H on average ' num2str(mean(H2HRange)) 'uV.\n'])

% fprintf(['N170 for ' phase ' X2N on average ' num2str(mean(X2NAmp)) 'uV 95CI [' num2str(X2NCI(1)) ' ' num2str(X2NAmpCI(2))  '] at ' num2str(mean(X2NPeakTime)) 'ms.\n'])
% fprintf(['N170 for ' phase ' X2H on average ' num2str(mean(X2HAmp)) 'uV 95CI [' num2str(X2HCI(1)) ' ' num2str(X2HAmpCI(2))  '] at ' num2str(mean(X2HPeakTime)) 'ms.\n'])