% loads unfold results, adds ufresult.data and calculates the
% SE over trials for each subject

function se = unfold_calculateSE(subj,phase)
cd '/net/store/nbp/users/agert/itw'
init_itw

addpath('/net/store/nbp/users/agert/unfold')
init_unfold

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

filtType='causal';
% time ranges for P100 and N170
t_p100=[0.07 0.150];
t_n170=[0.120 0.200];


%% LOADING
codingscheme = 'mashup';


%% load ufresults
ufPath = fullfile(projectFolder,['Daten/EEG/' phase '/unfold_' phase '/' filtType '/']);
ufPath = fullfile(ufPath,['ufresult_subj' num2str(subj) ,'_allElec_' codingscheme '_-500-1000_regularized0.mat']);

fprintf('\n \n loading subject %i \n',subj)

d = load(ufPath);
ufresult = d.ufresult;
clear d; % make some space

if strcmp(phase,'WLFO')
    ufam = uf_predictContinuous(ufresult,'predictAt',{{'sac_amplitude',[0.5 1 2.5 5 7.5 10]},{'fix_avgpos_y',linspace(500,1500,5)},{'fix_avgpos_x',linspace(900,3000,5)},{'samebox_1',1}});
    ufam = uf_addmarginal(ufam);
elseif strcmp(phase,'Lab')
    ufam=ufresult;
end


%% load EEG data
cfg = struct();
cfg.subject = num2str(subj);
cfg.filtering='causal';
cfg.mainfolder = ['/net/store/nbp/projects/IntoTheWild/Daten/EEG/' phase '/VP' cfg.subject '/preprocessed/' cfg.filtering];
cfg.file = ['3_ITW_' phase '_subj' cfg.subject '_channelrejTriggersXensor.set'];
tmp = pop_loadset('filename',['2_ITW_' phase '_subj' cfg.subject '_bandpass_resample_deblank.set'],'filepath',cfg.mainfolder,'loadmode','info');
tmp.chanlocs(find(strcmp({tmp.chanlocs.labels},'CzAmp2'))).labels = 'Iz';
tmp=pop_chanedit(tmp, 'lookup','/net/store/nbp/projects/IntoTheWild/EEG_analysis/eeglab13_5_4b/plugins/dipfit2.3/standard_BESA/standard-10-5-cap385.elp');
cfg.urchanlocs= tmp.chanlocs;
cfg.avgref=1;
cfg.saccade=1;
cfg.regularize=0;
tmp = load(fullfile(cfg.mainfolder,['6_ITW_' phase '_subj' cfg.subject '_ICAcleancont.mat']));
cfg.badcomponents = tmp.comps_to_rej;
cfg.timewin=[-0.5 1];

%load EEG
EEG = pop_loadset('filename',cfg.file,'filepath',cfg.mainfolder);
EEG =  pop_select(EEG,'nochannel',{'AUX1','AUX2'});

% Load ICA and clean
mod = loadmodout12(fullfile(cfg.mainfolder,'amica'));
EEG.icaweights = mod.W;
EEG.icasphere = mod.S;
EEG.icawinv = [];EEG.icaact = [];EEG.icachansind = [];
EEG = eeg_checkset(EEG);
EEG = pop_subcomp(EEG,cfg.badcomponents);

% interpolate missing channels
% ET channels (Gaze etc.) not interpolated
alldel = {'CzAmp2' 'BIP1' 'BIP2' 'BIP3' 'BIP4' 'BIP5' 'BIP6' 'BIP7' 'BIP8' 'AUX1' 'AUX2' 'AUX3' 'AUX4' 'AUX5' 'AUX6' 'AUX7' 'AUX8','TIME', 'L_GAZE_X', 'L_GAZE_Y', 'L_AREA','R_GAZE_X', 'R_GAZE_Y', 'R_AREA', 'INPUT', 'L-GAZE-X', 'L-GAZE-Y', 'L-AREA','R-GAZE-X', 'R-GAZE-Y', 'R-AREA', 'INPUT'};
%remove VEOG from EEG, as there is no corresponding location
EEG = pop_select(EEG, 'nochannel', {'VEOG'});

if cfg.avgref
    %calculate avg ref
    EEG = pop_reref( EEG, []); %Participants’ averages were then re-referenced to a common average reference. (Rossion & Caharel, 2011)
end

%interpolate missing channels
idxs = ~ismember({cfg.urchanlocs.labels},alldel);
EEG= pop_interp(EEG,cfg.urchanlocs(idxs),'spherical');


%% add the EEG data to the unfold data
ufresult.data=EEG.data;

clear EEG mod; % make some space


%% find P100 and N170 peaks

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

P100Start=find(ufresult.times>t_p100(1),1,'first');
P100Stop=find(ufresult.times<t_p100(2),1,'last');
N170Start=find(ufresult.times>t_n170(1),1,'first');
N170Stop=find(ufresult.times<t_n170(2),1,'last');

paramNames=ufresult.unfold.colnames;
if strcmp(phase,'WLFO')
    [~,intPos]   =find(ismember(paramNames,'(Intercept)'));
elseif strcmp(phase,'Lab')
    [~,intPos]   =find(ismember(paramNames,'2_(Intercept)'));
end
[~,currHFPos]=find(ismember(paramNames,'curr_HF'));

ufFace=mean(ufam.beta(chan,:,intPos)+ufresult.beta(chan,:,currHFPos));
ufObj =mean(ufam.beta(chan,:,intPos)-ufresult.beta(chan,:,currHFPos));
clear ufam;% make some space

[~,P100Face]=max(ufFace(P100Start:P100Stop));
P100Face=P100Start+P100Face-1;
[~,N170Face]=min(ufFace(N170Start:N170Stop));
N170Face=N170Start-1+N170Face; %-1 as correction for shift by one position

[~,P100Obj]=max(ufObj(P100Start:P100Stop));
P100Obj=P100Start+P100Obj-1;
[~,N170Obj]=min(ufObj(N170Start:N170Stop));
N170Obj=N170Start-1+N170Obj; %-1 as correction for shift by one position


%% calculate SE
contrast = zeros(7,size(ufresult.unfold.beta_dc,2),size(ufresult.unfold.beta_dc,3));

%y = b0==N1 + bHF==N1 - (b0==P1 + bHF==P1)

contrast(1, N170Face,intPos)    = 1;
contrast(1, N170Face,currHFPos) = 1;

contrast(2, P100Face,intPos)   = (1);
contrast(2, P100Face,currHFPos)= (1);

%y = b0==N1 - bHF==N1 - (y = b0==P1 - bHF==P1)
contrast(3, N170Obj,intPos)    = 1;
contrast(3, N170Obj,currHFPos) = -1;

contrast(4, P100Obj,intPos)   = (1);
contrast(4, P100Obj,currHFPos)= (-1);

contrast=uf_contrast_addmarginal(ufresult,contrast);
%y =  bHF==N1+ bHF==N1  - bHF==P1  - bHF==P1   - b0==P1 + b0==P1 + b0==N1 - b0==N1
contrast(5, :,:) =   contrast(1, :,:)-contrast(2, :,:); %N1Face - P1Face
contrast(6, :,:) =   contrast(3, :,:)-contrast(4, :,:); %N1Obje - P1Obje
contrast(7, :,:) =   contrast(5, :,:)-contrast(6, :,:); %P2PFace - P2Pobj

tic
[se] = uf_se(ufresult,'channels', chan,'contrast',contrast);

toc

save(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/' phase '/unfold_' phase '/' filtType 'unfoldSE_subj' num2str(subj) '.mat'],'se')
end

%end