WildLab/code/unfold/unfold_lab.m

function [ufresult]=unfold_lab(sub_num)
%%
cfg = struct();
cfg.subject = num2str(sub_num);
cfg.filtering='causal';
cfg.mainfolder = ['/net/store/nbp/projects/IntoTheWild/Daten/EEG/Lab/VP' cfg.subject '/preprocessed/' cfg.filtering];
tmp = load(fullfile(cfg.mainfolder,['4_ITW_Lab_subj' cfg.subject '_cleaningTimes.mat']));
cfg.cleantimes = tmp.rej(:,1:2);

cfg.file = ['3_ITW_Lab_subj' cfg.subject '_channelrejTriggersXensor.set'];
tmp = pop_loadset('filename',['2_ITW_Lab_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_Lab_subj' cfg.subject '_ICAcleancont.mat']));
cfg.badcomponents = tmp.comps_to_rej;

cfg.timewin=[-0.5 1];


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


%% Load ICA and clean
addpath('/net/store/nbp/projects/EEG/blind_spot/amica')
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

idxs = ~ismember({cfg.urchanlocs.labels},alldel);

%interpolate missing channels
EEG= pop_interp(EEG,cfg.urchanlocs(idxs),'spherical');


%% adjust for trigger delay
% as calculated from photosensor
if sub_num<36
    trigger_delay=0.0102; %10.2ms
else
    trigger_delay=0.0348; %34.8ms
end

for evt=1:length(EEG.event)
    if str2num(EEG.event(evt).type)<=5 %only for stim onset [1 2 3 4 5]
        assert(str2num(EEG.event(evt).type)>0)
        EEG.event(evt).latency = EEG.event(evt).latency + (trigger_delay*EEG.srate);
    end
end
EEG = eeg_checkset(EEG);

%% adjust cleaning times to be between 1 and length of data
[a,b]=find(cfg.cleantimes(:,1)<1);
if ~isempty(a)
    cfg.cleantimes(a,1)=1;
end

[a,b]=find(cfg.cleantimes(:,1)>length(EEG.data));
if ~isempty(a)
    cfg.cleantimes(a,:)=[];
end


[a,b]=find(cfg.cleantimes(:,2)>length(EEG.data));
if ~isempty(a)
    cfg.cleantimes(a,2)=length(EEG.data);
end
numstim=0;
%% for each fixation we need to save ALL the necessary information in the EEG structure
% for e = 1:length(EEG.event)
for e = length(EEG.event):-1:1
    
    % 1-5 stimulus
    % 180 stim offset
    % 200 exp start
    % 255 exp stop
    
    % if the trigger is a fixation
    if ismember(EEG.event(e).type, ...
            {'L_fixation' 'R_fixation'})

        EEG.event(e).urtype = EEG.event(e).type;
        
        % if it is a meaningful fixation
        
        EEG.event(e).type = 'fixation';
        
        continue

    end

    % if the trigger was the beginning of a trial
    if str2num(EEG.event(e).type) <= 5
        EEG.event(e).urtype = EEG.event(e).type;
        EEG.event(e).type = 'stimulus';
    end
    
    
    
    % recode all fix that are not Face or None to "OtherHeads"
    if strcmp(EEG.event(e).type, 'stimulus')
        isStim=1;
    else
        isStim=0;
    end
    
    
    
    if isStim==1
        
        numstim=numstim+1;
        
        if EEG.event(e).humanface==1
            EEG.event(e).curr='HF';
            
        elseif EEG.event(e).car==1
            EEG.event(e).curr='Car';
        elseif EEG.event(e).obj==1
            EEG.event(e).curr='Obj';
        elseif EEG.event(e).none==1
            EEG.event(e).curr='None';
        else
            %first in a trial has no prevprev, therefore we define it as
            %'outside' as outside will not be analysed (easiest solution)
            EEG.event(e).curr='outside';
        end
        
        if EEG.event(e).prevhumanface==1
            EEG.event(e).prev='HF';
        elseif EEG.event(e).prevcar==1
            EEG.event(e).prev='Car';
        elseif EEG.event(e).prevobj==1
            EEG.event(e).prev='Obj';
        elseif EEG.event(e).prevnone==1
            EEG.event(e).prev='BlockBeginning';
            
        end
        
        
        if EEG.event(e).prevprevhumanface==1
            EEG.event(e).prevprev='HF';
        elseif EEG.event(e).prevprevnone==1
            EEG.event(e).prevprev='None';
        else
            %first in a trial has no prevprev, therefore we define it as
            %'outside' as outside will not be analysed (easiest solution)
            EEG.event(e).prevprev='BlockBeginning';
        end
        
    end
  
end


%% save information about saccaded and fixations
if sum(strcmp({EEG.event.type},'L_saccade'))>0
    ixSac = find(strcmp({EEG.event.type},'L_saccade'));
else
    ixSac = find(strcmp({EEG.event.type},'R_saccade'));
end

ixFix = find(strcmp({EEG.event.type},'fixation'));

%find blinks
ixBlk=find(cellfun(@(x)strcmp(x,'R_blink'),{EEG.event.type}));
if isempty(ixBlk)
    ixBlk=find(cellfun(@(x)strcmp(x,'L_blink'),{EEG.event.type}));
end


latSac = [EEG.event(ixSac).latency];
latFix = [EEG.event(ixFix).latency];
latBlk = [EEG.event(ixBlk).latency];

endSac = latSac+[EEG.event(ixSac).duration];
endBlk = latBlk+[EEG.event(ixBlk).duration];

% put the velocity for blink sac/data loss sections to NaN
for ixS=1:length(ixSac)
    
    blkStarAftertSacStart=find(latBlk>=latSac(ixS),1,'first');
    blkStopBeforeSacStop=find(endBlk<=endSac(ixS),1,'last');
    if blkStarAftertSacStart== blkStopBeforeSacStop
        EEG.event(ixSac(ixS)).sac_vmax=NaN;
    end
end



for e = 1:length(latFix)
    % for each fixation, find the closest saccade before the event and copy
    % the saccade information
    
    % when  Sac --> Fix, then lat(Fix) - lat(Sac) should be positive
    ix = find((latFix(e) - latSac) > 0,1,'last');
    
    if ~isempty(ix)
        for fn = {'sac_endpos_x','sac_endpos_y','sac_amplitude','sac_startpos_x','sac_startpos_y','sac_vmax'}
            EEG.event(ixFix(e)).(fn{1})= EEG.event(ixSac(ix)).(fn{1});
        end
    end
    
end

% remove all fixation with saccade velocity of 0 (= blink saccades)
for ixF=1:length(ixFix)
    if EEG.event(ixFix(ixF)).sac_vmax==0
        EEG.event(ixFix(ixF)).type='blinkFix';
    end
end

ixFix = find(strcmp({EEG.event.type},'fixation'));
latFix =  [EEG.event(ixFix).latency];

% remove all fixations close to a blink/data loss section
if ~isempty(ixBlk)
    ixBlk=fliplr(ixBlk);
    for iB=ixBlk
        blkStop=EEG.event(iB).latency+EEG.event(iB).duration;
        delFix=find(abs(latFix-blkStop)<50);
        
        if ~isempty(delFix)
            for dF=1:length(delFix)
                EEG.event(ixFix(delFix(dF))).type='blinkFix';
            end
        end
    end
end

for e = 1:length(latSac)
    % for each fixation, find the closest saccade before the event and copy
    % the saccade information
    
    % when  Sac --> Fix, then lat(Fix) - lat(Sac) should be positive
    ix = find((latSac(e) - latFix) > 0,1,'last');
    if ~isempty(ix) && ~isempty( EEG.event(ixFix(ix)).prev)
        for fn = {'prev','curr'}
            EEG.event(ixSac(e)).(fn{1})= EEG.event(ixFix(ix)).(fn{1});
            
        end
        EEG.event(ixSac(e)).type='saccade';
    end
    
end


%% Remove large saccades
fprintf('\n\n %u large saccades (>40) will be removed. \n\n\n', length(find([EEG.event.sac_amplitude]>40)))
EEG.event(cellfun(@(x)(x>40),{EEG.event.sac_amplitude})) = [];



%% Config
EEG.unfold = [];
cfgDesign = [];

cfgDesign.splinespacing = 'quantile';

cfgDesign.eventtypes = {
    {'fixation'}   % fixation onset
    {'stimulus'}}; % stimulus onset

if sub_num==51 %51 has only 2 relevant sac before fix - maybe due to stopping recording between trials
    cfgDesign.formula = {
        'y~1 + spl(fix_avgpos_x,5) + spl(fix_avgpos_y,5)'
        'y~1 + cat(curr) + cat(prev)+ humanface:prevhumanface'};
else
    cfgDesign.formula = {
        'y~1 + spl(fix_avgpos_x,5) + spl(fix_avgpos_y,5) +  spl(sac_amplitude,5)'
        'y~1 + cat(curr) + cat(prev)+ humanface:prevhumanface'};
end


cfgDesign.codingschema = 'effects';
cfgDesign.mashup=1;

% if codingschema=effects -> reference cat is the last one
% if codingschema=reference -> reference cat is the first one
if strcmp(cfgDesign.codingschema, 'effects')
    cfgDesign.categorical = {
        'curr',{'Car','HF','Obj'};
        'prev',{'BlockBeginning','Car','HF','Obj'};
        };
elseif strcmp(cfgDesign.codingschema, 'reference')
    cfgDesign.categorical = {
        'curr',{'Obj','HF','Car','BlockBeginning'};
        'prev',{'Obj','HF','Car','BlockBeginning'};
        };
end

cfgFit = [];
cfgFit.precondition = 1;
cfgFit.lsmriterations = 1500;
% if subset of channels is wanted:
% elecfind = @(x)find(strcmp(x,{EEG.chanlocs.labels}));
% cfgFit.channel = [elecfind('PO8')];
cfgFit.channel  = 1:length(EEG.chanlocs);
if cfg.regularize
    cfgFit.method='glmnet';
    cfgFit.glmnetalpha = 0;
end

EEG = uf_designmat(EEG,cfgDesign);

% if we want a mashup design, exchange everything <1 with 0, in all
% uninteresting parameters
if cfgDesign.mashup
    parampos=find(ismember(EEG.unfold.colnames,{'prev_Car','curr_Car','prev_BlockBeginning'}));
    
    for p=1:length(parampos)
        pos=find(EEG.unfold.X(:,parampos(p))<0);
        EEG.unfold.X(pos,parampos(p))=0;
    end
    
    %the interaction is basically only the multiplicatiom of curr and prev
    interactpos=find(ismember(EEG.unfold.colnames,{'humanface:prevhumanface'}));
    mainpos =find(ismember(EEG.unfold.colnames,{'curr_HF','prev_HF'}));
    EEG.unfold.X(:,interactpos)=EEG.unfold.X(:,mainpos(1)).*EEG.unfold.X(:,mainpos(2));
end


%% Add cyclical regressor
% if cfgDesign.cyclicalAngle
%     EEG = scenes_add_cyclical(EEG);
% end


%% Timeshift the design matrix
for t=1:size(cfg.timewin,1)
    cfgTimeshift = [];
    cfgTimeshift.timelimits = cfg.timewin(t,:);
    
    EEG = uf_timeexpandDesignmat(EEG,cfgTimeshift);
    
    EEG = uf_continuousArtifactExclude(EEG,'winrej',cfg.cleantimes);
    
    
    %% Fit it iteratively
    EEG= uf_glmfit(EEG,cfgFit);
 
    
    %% Make a massive uni-variate fit without de-convolution
    EEGepoch = uf_epoch(EEG,'winrej',cfg.cleantimes,'timelimits',cfgTimeshift.timelimits);
    EEGepoch= uf_glmfit_nodc(EEGepoch);

    
    %% Save the data
    
    ufresult = uf_condense(EEGepoch);
    
    save(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/Lab/unfold_Lab/' cfg.filtering '/ufresult_subj' num2str(sub_num) '_allElec_mashup_' num2str(cfg.timewin(t,1)*1000) '-' num2str(cfg.timewin(t,2)*1000) '_regularized' num2str(cfg.regularize) '.mat'],'ufresult')
    
    
    disp('successfully saved')
end