WildLab/code/unfold/unfold_freeviewing.m

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

cfg.file = ['3_ITW_WLFO_subj' cfg.subject '_channelrejTriggersXensor.set'];
tmp = pop_loadset('filename',['2_ITW_WLFO_subj' cfg.subject '_bandpass_resample_deblank.set'],'filepath',cfg.mainfolder,'loadmode','info');
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');
tmp.chanlocs(find(strcmp({tmp.chanlocs.labels},'CzAmp2'))).labels = 'Iz';
cfg.urchanlocs= tmp.chanlocs;
cfg.avgref=1;
cfg.saccade=1;
cfg.regularize=0;

tmp = load(fullfile(cfg.mainfolder,['6_ITW_WLFO_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');
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

%% 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
    
    
    % if the trigger was the beginning of a trial
    if str2num(EEG.event(e).type) <= 171
        EEG.event(e).urtype = EEG.event(e).type;
        EEG.event(e).type = 'stimulus';
        continue
    end
    
    % 1-171 stimulus
    % 180 stim offset
    % 200 exp start
    % 255 exp stop
    
    % if the trigger is a fixation
    if ismember(EEG.event(e).type, ...
            {'HFtoHF'    'HFtoHH'    'HFtoNH'    'HFtoOL'    'HFtoNone'  'HFtoOS' ...
            'HHtoHF'    'HHtoHH'    'HHtoNH'    'HHtoOL'    'HHtoNone'   'HHtoOS' ...
            'NHtoHF'    'NHtoHH'    'NHtoNH'    'NHtoOL'    'NHtoNone'   'NHtoOS' ...
            'OLtoHF'    'OLtoHH'    'OLtoNH'    'OLtoOL'    'OLtoNone'   'OLtoOS' ...
            'NonetoHF'  'NonetoHH'  'NonetoNH'  'NonetoOL'  'NonetoNone' 'NonetoOS' ...
            'OStoHF'    'OStoHH'    'OStoNH'    'OStoOL'    'OStoNone'   'OStoOS' ...
            })
        
        
        EEG.event(e).urtype = EEG.event(e).type;
        
        % if it is a meaningful fixation
        toIDX=strfind(EEG.event(e).type,'to');
        
        if toIDX>0
            
            EEG.event(e).curr=EEG.event(e).type((toIDX+2):end);
            EEG.event(e).prev=EEG.event(e).type(1:toIDX-1);
            
        else error(['You forgot to code a fixation type:' EEG.event(e).type])
        end
        
        EEG.event(e).type = 'fixation';
        
        %         continue      end
        
        
    end
    
    % recode all fix that are not Face or None to "OtherHeads"
    if strcmp(EEG.event(e).type, 'fixation')
        isFix=1;
    else
        isFix=0;
    end
    
    if isFix==1
        if ismember(EEG.event(e).curr,{'OS'})==1
            EEG.event(e).curr='outside';
        elseif ismember(EEG.event(e).curr,{'HF','None'})==0
            EEG.event(e).curr='OtherHeads';
            
        end
        
        if ismember(EEG.event(e).prev,{'OS'})==1
            EEG.event(e).prev='outside';
        elseif ismember(EEG.event(e).prev,{'HF','None'})==0
            EEG.event(e).prev='OtherHeads';
        end
        
        if EEG.event(e).prevprevhumanface==1
            EEG.event(e).prevprev='HF';
        elseif EEG.event(e).prevprevnone==1
            EEG.event(e).prevprev='None';
        elseif EEG.event(e).prevprevoutside==1
            EEG.event(e).prevprev='outside';
        elseif EEG.event(e).prevprevhumanhead==1
            EEG.event(e).prevprev='OtherHeads';
        elseif EEG.event(e).prevprevnonhuman==1
            EEG.event(e).prevprev='OtherHeads';
        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='outside';
        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');
    
    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

% 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','samebox'}
            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

cfgDesign.formula = {
    'y~1 + cat(curr) + cat(prev)+ cat(samebox) + humanface:prevhumanface + spl(fix_avgpos_x,5) + spl(fix_avgpos_y,5) +  spl(sac_amplitude,5)',
    'y~1'};

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',{'outside','HF','OtherHeads','None'};
        'prev',{'outside','HF','OtherHeads','None'};
        'samebox',{1,0};
        };
elseif strcmp(cfgDesign.codingschema, 'reference')
    cfgDesign.categorical = {
        'curr',{'None','HF','OtherHeads','outside'};
        'prev',{'None','HF','OtherHeads','outside'};
        'samebox',{0,1};
        };
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,{'curr_OtherHeads','prev_OtherHeads','curr_outside','prev_outside'}));
    
    for p=1:length(parampos)
        pos=find(EEG.unfold.X(:,parampos(p))<0);
        EEG.unfold.X(pos,parampos(p))=0;
    end
    
    sameboxpos=find(strncmp(EEG.unfold.colnames,'same',4));
    sbpos0=find(EEG.unfold.X(:,sameboxpos)<0);
    EEG.unfold.X(sbpos0,sameboxpos)=0;
    
    sameboxnpos=find(strncmp(EEG.unfold.colnames,'sameboxn',8));
    psbpos0=find(EEG.unfold.X(:,sameboxnpos)<0);
    EEG.unfold.X(psbpos0,sameboxnpos)=0;
    
    %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


%% 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/WLFO/unfold_WLFO/' 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