braindynamicsusyd / braintrak Goto Github PK

Should be able to do feather.interactive_fit
If feather.latest ~= 1, then it should re-run interactive fit on feather.subrange(feather.latest)

The gui should show the original and new chisq values

Should be mcmc.fit.single rather than mcmc.fit_single

% detecteyemovements() - detect saccades & fixations in eye tracking data.
% Saccade detection is based on the algorithm by
% Engbert & Mergenthaler (2006). Saccades are defined as
% (monocular or binocular) outliers in 2D velocity space.
% Velocity thresholds for saccade detection are determined
% adaptively as a multiple of the (median-based) SD of all
% data samples in the epoch. Fixations are defined as
% intervals in-between saccades. Eye movements can be added
% as new events to EEGLAB's event structure. For various
% other options, see below.
%
% Usage:
% >> EEG = detecteyemovements(EEG,left_eye_xy,right_eye_xy,vfac,mindur,...
% degperpixel,smooth,globalthresh,clusterdist,clustermode,
% plotfig,writesac,writefix)
%
% Required inputs:
% EEG - [string] EEG struct, also containing synchronized eye
% tracking data (see pop_importeyetracker)
% left_eye_xy - [vector of two channel indices],
% specifying channel indices of X- (first value) and
% Y-component (second value) of left eye. Leave empty []
% if the left eye was not recorded.
% right_eye_xy - [vector of two channel indices],
% specifying channel indices of X- (first value) and
% Y-component (second value) of right eye. Leave empty []
% if the right eye was not recorded.
% vfac - [double] velocity factor ("lambda") to determine
% the velocity threshold for saccade detection
% (cf. Engbert & Mergenthaler, 2006)
% mindur - [integer] minimum saccade duration (in samples)
% (cf. Engbert & Mergenthaler, 2006)
% degperpixel - [double] visual angle of one screen pixel
% if this value is left empty [], saccade characteristics
% are reported in the original data metric (pixel?)
% instead of in degrees of visual angle
% smooth - [0/1] if set to 1, the raw data is smoothed over a
% 5-sample window to suppress noise
% noise. Recommended for high native ET sampling rates.
% globalthresh - [0/1]. Use the same thresholds for all epochs?
% 0: Adaptive velocity thresholds are computed
% individually for each data epoch.
% 1: Adaptive velocity thresholds are first computed for
% each epoch, but then the mean thresholds are applied to
% each epochs (i.e. same detection parameters are used for
% all epochs). Setting is irrelevant if the input data is
% still continuous (= only one data epoch).
% clusterdist - [integer] value in sampling points that defines the
% minimum allowed fixation duration between two saccades.
% If the off- and onsets of two temp. adjacent sacc. are
% closer together than 'clusterdist' samples, these
% saccades are regarded as a "cluster" and treated
% according to the 'clustermode' setting (see below).
% clusterdist is irrelevant if clustermode == 1.
% clustermode - [1,2,3,4]. Integer between 1 and 4.
% 1: keep all saccades, do nothing
% 2: keep only first saccade of each cluster
% 3: keep only largest sacc. of each cluster
% 4: combine all movements into one (longer) saccade
% this new saccade is defined as the movement that
% occurs between the onset of the 1st saccade in the
% cluster and the offset of the last sacc. in cluster
% WARNING: CLUSTERMODE 4 is experimental and untested!
% plotfig - [0/1] Show a figure with eye movement properties?
% 0: do not plot a figure.
% 1: plot a figure displaying properties of detected
% saccades & fixations
% writesac - [0/1]: Add saccades to EEG.event?
% 0: detect saccades, but do not store them in EEG.event.
% 1: add detected saccades as new events to EEG.event.
% writefix - [0/1]: Add fixations to EEG.event?
% 0: detect fixations, but do not add them to EEG.event.
% 1: add detected fixations as new events to EEG.event.
% Note: It is recommended to first test the parameters of
% saccade/fixation detection without adding events.
% For this, set writesac and writefix to 0.
%
% Outputs:
% EEG - EEG structure. If writesac or writefix were set to 1,
% the EEG structure (EEG.event/EEG.urevent/EEG.epoch)
% will contain additional "saccade" and "fixation" events
% with their respective properties
%
% See also: vecvel, velthresh, microsacc_plugin, binsacc, saccpar,
% mergesacc, addevents
%
%
% MAJOR UPDATE 03/2018: bad ET intervals, as marked by "bad_ET" events
% in EEG.event (see function pop_rej_contin.m) are excluded from eye movement
% detection. Bad intervals are now ignored when estimating velocity thresholds
% and saccade/fixation events overlapping with these intervals are removed
%
% An example call of the function might look like this:
% >> EEG = detecteyemovements(EEG,[],[33 34],6,4,0.037,1,0,25,4,1,1,0)
%
% In this example, the eye position data for the right eye is stored in
% channels 33 (horiz.) and 34 (vertical). The left eye was not recorded.
% The velocity threshold is set to 6 times the (median-based)
% SD of all velocity samples in the epoch. The minimum duration of
% saccades to be detected is 4 samples. In the experiment, one screen
% pixel corresponded to 0.037 degrees of visual angle.
% The raw data is smoothed prior to saccade detection (smooth: 1).
% Adaptive velocity thresholds (X and Y-threshold for each eye) are
% determined individually for each data epoch (globalthresh: 0). For saccades
% separated by fixations of less than 25 samples, only the first saccade
% is kept (clusterdist: 25, clustermode: 2). A figure with the
% results is plotted. Detected saccades are stored as new events in
% EEG.event, but fixations are not stored.
%
% The eye movement detection is based on:
%
% Engbert, R., & Kliegl, R. (2003). Microsaccades uncover the orientation
% of covert attention. Vision Research, Vol. 43, 1035-1045
%
% ...as well as...
%
% Engbert, R., & Mergenthaler, K. (2006). Microsaccades are triggered by
% low retinal image slip, PNAS, Vol. 103 (18), 7192-7197
%
% Author: od
% Copyright (C) 2009-2018 Olaf Dimigen, HU Berlin
% [email protected]

% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, 51 Franklin Street, Boston, MA 02110-1301, USA

function EEG = detecteyemovements(EEG,left_eye_xy,right_eye_xy,vfac,mindur,degperpixel,smooth,globalthresh,clusterdist,clustermode,plotfig,writesac,writefix)

allsac = [];
allfix = [];

% data of which eye is available?
ldata = false;
if length(left_eye_xy) == 2, ldata = true; end
rdata = false;
if length(right_eye_xy) == 2, rdata = true; end

if length(left_eye_xy) == 1 || length(right_eye_xy) == 1
error('%s(): For each recorded eye, horizontal (X) and vertical (Y) gaze channel must be specified.',mfilename);
end

nsample = size(EEG.data,2);
nepochs = size(EEG.data,3);
badepochs = zeros(nepochs,1);
nbadsmp = 0;

% preallocate storage for sacc. detection thresholds
l_msdx = NaN(nepochs,1);
l_msdy = NaN(nepochs,1);
r_msdx = NaN(nepochs,1);
r_msdy = NaN(nepochs,1);

% warn message if back-to-back saccades are detected
clusterwarning = false;

% critical bugfix 2013-10-01, by OD
% due to bug in third-party function "smoothdata":
% function smoothdata() was removed from the toolbox
% function vecvel() was updated to incorporate different levels of smoothing:
% options:
% - smoothlevel 0: no smoothing, simple diff()
% - smoothlevel 1: 3-point window
% - smoothlevel 2: 5-point window
if smooth
smoothlevel = 2; % 5-point smoothing
else
smoothlevel = 0; % no smoothing
end

%% screen feedback

% detection feedback
if ~ldata && ~rdata
error('%s(): Please correctly specify channels containing eye tracking data.\nFor each recorded eye, both horiz. (X) and vert. (Y) channel must be specified.',mfilename)
else
% summary of detection parameters
fprintf('\n--------------------------------------------------------------------')
fprintf('\nDetecting saccades after Engbert & Mergenthaler (2006)\n')
fprintf('\nVelocity threshold factor (vfac): %.2f SD',vfac);
fprintf('\nMinimum saccade duration (mindur): %.2f samples (%.2f ms)',mindur,mindur1000/EEG.srate);
if ~isempty(degperpixel) | isnan(degperpixel) % bugfix 2016-11-12 by OD: added case if degperpixel = NaN (from GUI input)
fprintf('\nVisual angle per screen pixel: %f°',degperpixel);
metric = 'deg';
else
fprintf('\nWARNING: No input provided for degperpixel!\nSpatial saccade properties are given in original metric (pixel?)');
degperpixel = 1;
metric = 'pix';
end
if nepochs < 2, fprintf('\n-- Using continuous data.'), else fprintf('\n-- Using epoched data.'); end
if ldata && rdata, fprintf('\n-- Using binocular data:'); else fprintf('\n-- Using monocular data:'); end
if ldata
fprintf('\n\tLeft horiz.: "%s"',EEG.chanlocs(left_eye_xy(1)).labels);
fprintf('\n\tLeft verti.: "%s"',EEG.chanlocs(left_eye_xy(2)).labels);
end
if rdata
fprintf('\n\tRight horiz.: "%s"',EEG.chanlocs(right_eye_xy(1)).labels);
fprintf('\n\tRight verti.: "%s"',EEG.chanlocs(right_eye_xy(2)).labels);
end
if smooth, fprintf('\n-- Raw data is smoothed in 5-sample window.'); else fprintf('\n-- Raw data is not smoothed.'); end
if nepochs > 1
if globalthresh
fprintf('\n-- Velocity thresholds computed globally across all %i epochs',nepochs);
else
fprintf('\n-- Velocity thresholds computed individually for each epoch');
end
end
fprintf('\n-- Treatment of saccade clusters:');
switch clustermode
case 1
fprintf('\n\tAll saccades are kept');
case 2
fprintf('\n\tSaccades separated by fixations < %i ms are clustered.',clusterdist(1000/EEG.srate));
fprintf('\n\tFirst sacc. of each cluster is kept.');
case 3
fprintf('\n\tSaccades separated by fixations < %i ms are clustered.',clusterdist*(1000/EEG.srate));
fprintf('\n\tLargest sacc. of each cluster is kept.');
case 4
fprintf('\n\tSaccades separated by fixations < %i ms are clustered.',clusterdist*(1000/EEG.srate));
fprintf('\n\tClusters are combined into one saccade.');
otherwise
error('%s(): Unknown input for clustermode, should be: 1,2,3,4.',mfilename)
end
if plotfig, fprintf('\n-- A figure with eye movement properties is plotted.'); end
if writesac && writefix, fprintf('\n-- Saccades and fixations will be added to EEG.event.');
elseif writesac, fprintf('\n-- Saccades will be added to EEG.event.');
elseif writefix, fprintf('\n-- Fixations will be added to EEG.event.');
else, fprintf('\n-- Saccades and fixations are detected, but NOT stored anywhere.');
end
end

%% get "bad_ET" intervals from EEG.event structure
badvector = zeros(1,size(EEG.data,2)*size(EEG.data,3));
ix_badETevent = find(ismember({EEG.event.type},'bad_ET'));
if ~isempty(ix_badETevent)

fprintf('\n--------------------------------------------------------------------')
fprintf('\nFound \"bad_ET\" events in EEG.event.\nThese intervals will be ignored for saccade detection!')
fprintf('\n--------------------------------------------------------------------')
bad_lat     = [EEG.event(ix_badETevent).latency];
bad_dur     = [EEG.event(ix_badETevent).duration];
bad_ET      = [bad_lat; bad_dur]';
bad_ET(:,3) = bad_ET(:,1)+bad_ET(:,2)-1;

% create long vector (as long as EEG) indicating bad samples
for j = 1:size(bad_ET,1)
    badvector(bad_ET(j,1):bad_ET(j,3)) = 1;
end

end
% reshape "badvector" to 3D if data is already epoched
badvector = reshape(badvector,1,size(EEG.data,2),size(EEG.data,3));

%% pre-compute saccade velocity thresholds for all epochs
% this enables the option to use the same threshold for all epochs, not
% provided by the original Engbert & Mergenthaler implementation
% new 03/2018: exclude bad ET intervals from velocity estimation
fprintf('\nComputing adaptive velocity thresholds...')
% if any(badvector(:))
% fprintf('\n-- Found %i "bad_ET" events marking bad eye-tracking intervals in EEG.event.',length(ix_badETevent))
% fprintf('\n-- These intervals (%.2f%% of data) will be ignored when computing velocity thresholds.',(sum(badvector)/length(badvector(:))*100))
% end
for e=1:nepochs
ix_goodET = ~badvector(:,:,e); % get index of non-bad ET samples
if ldata
l = EEG.data([left_eye_xy(1) left_eye_xy(2)],ix_goodET,e)';
vl = vecvel(l,EEG.srate,smoothlevel);
[l_msdx(e), l_msdy(e)] = velthresh(vl);
end
if rdata
r = EEG.data([right_eye_xy(1) right_eye_xy(2)],ix_goodET,e)';
vr = vecvel(r,EEG.srate,smoothlevel);
[r_msdx(e), r_msdy(e)] = velthresh(vr);
end
end

%% detect saccades & fixations
for e=1:nepochs

%% saccades of left eye
sac = [];
if ldata
    l = EEG.data([left_eye_xy(1) left_eye_xy(2)],:,e)'; % don't exclude bad_ET samples here (otherwise discont. high-velocity jumps in time series)
    % bad/missing samples in eye track?
    badsmp = sum(sum(l<=0)); if badsmp > 0, badepochs(e) = 1; nbadsmp = nbadsmp + badsmp; end
    vl = vecvel(l,EEG.srate,smoothlevel); % get eye velocities
    % detect monocular saccades
    if globalthresh % use precomputed velocity thresholds (mean of all epochs)
        sacL = microsacc_plugin(l,vl,vfac,mindur,mean(l_msdx),mean(l_msdy));
    else % compute velocity thresholds from this epoch only
        sacL = microsacc_plugin(l,vl,vfac,mindur,l_msdx(e),l_msdy(e));        
    end
end
%% saccades of right eye
if rdata
    r = EEG.data([right_eye_xy(1) right_eye_xy(2)],:,e)';
    % bad/missing samples in eye track?
    badsmp = sum(sum(r<=0)); if badsmp > 0, badepochs(e) = 1; nbadsmp = nbadsmp + badsmp; end        
    vr = vecvel(r,EEG.srate,smoothlevel); % get eye velocities
    % detect monocular saccades
    if globalthresh % use precomputed velocity thresholds (mean of all epochs)
        sacR = microsacc_plugin(r,vr,vfac,mindur,mean(r_msdx),mean(r_msdy));
    else 
        % compute velocity thresholds from this epoch only
        sacR = microsacc_plugin(r,vr,vfac,mindur,r_msdx(e),r_msdy(e));             
    end
end

%% binocular saccades
if ldata && rdata
    [sac, tmp, tmp] = binsacc(sacL,sacR); clear tmp
    sac = saccpar(sac); % average saccade characteristics of both eyes
    sac = mergesacc(sac,(l+r)./2,clusterdist,clustermode); % merge nearby saccades (e.g. glissades)
elseif ldata
    sac = sacL; clear sacL;
    sac = saccpar([sac sac]);
    sac = mergesacc(sac,l,clusterdist,clustermode);
elseif rdata
    sac = sacR; clear sacR;
    sac = saccpar([sac sac]);
    sac = mergesacc(sac,r,clusterdist,clustermode);    
end

%% update various saccade metrics
if ~isempty(sac)
    
    % define saccade duration as difference between saccade offset and
    % saccade onset sample. In saccpar(), monocular saccade durations
    % of both eyes are averaged, leading to uneven values (e.g.: 10.5
    % samples) different from the difference between onset and offset
    % values (which are the monocular extremes).
    % Instead: use difference between offset and onset
    sac(:,3) = sac(:,2)-sac(:,1)+1;
    
    % report saccade velocity/distance/amplitude as visual angles
    sac(:,[5 6 8]) = sac(:,[5 6 8]) .* degperpixel;
    
    % report saccade angles in degree rather than radians
    sac(:,[7 9]) = sac(:,[7 9]) * 180/pi;
    
    % add index of corresp. data epoch
    sac(:,10) = e;
    
    % store screen location for start and end of saccade
    if ldata && rdata
        gazexy = (l+r)./2; % binoc. recordings: average across eyes
    elseif ldata
        gazexy = l;
    elseif rdata
        gazexy = r;
    end
    % get position immediatly before sacc. onset and after sacc. offset
    startsmp = sac(:,1)-1; endsmp = sac(:,2)+1;
    if startsmp(1) < 1, startsmp(1,1) = 1; end
    if endsmp(end) > size(gazexy,2), endsmp(end) = size(gazexy,2); end
    sac(:,11) = gazexy(startsmp,1);
    sac(:,12) = gazexy(startsmp,2);
    sac(:,13) = gazexy(endsmp  ,1);
    sac(:,14) = gazexy(endsmp  ,2);
end

% columns of [sac]:
% 1: saccade onset (sample)
% 2: saccade offset (sample)
% 3: duration (samples)
% 4: delay between eyes (samples)
% 5: vpeak (peak velocity)
% 6: saccade "distance" (eucly. dist. between start and endpoint)
% 7: saccade angle (based on saccade "distance")
% 8: saccade "amplitude" (eucly. dist. of min/max in full saccade trajectory)
% 9: saccade angle (based on saccade "amplitude")
%10: index of corresponding data epoch (1 in case of contin. data)
%11: horizontal (x) gaze position before start of saccade (pixel)
%12: vertial (y) gaze position before start of saccade (pixel)
%13: horizontal (x) gaze position after end of saccade (pixel)
%14: vertical (y) gaze position after end of saccade (pixel)


%% remove saccades that occured during "bad_ET" intervals [*]
% Delete all saccades whos onset or offset occurs during a 
% bad_ET interval (signal jumps to blinks are otherwise detected as 
% saccades)
%
% [*] Note: It is not clear/trivial how to treat blinks in the context 
% of EM detection, since we do not know what really happened during a
% blink or loss of the signal. For example, is a long fixation with a 
% blink in the middle really two fixations? Or should it be treated as 
% one long fixation? Or should all fixations that are ended or started 
% by blinks or are interrupted by them be completely removed?
% The solution here is preliminary: I first remove all saccades that 
% started or ended during a blink interval. The remaining saccades are 
% used to define fixatinos. Finally, all fixations are removed that 
% overlap with "bad_ET" intervals in EEG.event. So this approach 
% removes any eye movement event that overlaps with a bad_ET interval. 
% OD, 2018-03-08
    
%% delete saccades starting/ending during "bad_ET" intervals
badETsmp = find(badvector(:,:,e));
if ~isempty(sac)
    ix_fakesac = find(ismember(sac(:,1),badETsmp) | ismember(sac(:,2),badETsmp));
    sac(ix_fakesac,:) = [];
%         if ~isempty(ix_fakesac)
%             fprintf('\n\n-- Removed %i saccades that occured during "bad_ET" intervals',length(ix_fakesac))
%         end
end
   
%% get fixations
nsac = size(sac,1);
fix = [];
if ~isempty(sac)
    if nsac > 1
        for f = 1:nsac-1
            
            fix(f,1) = sac(f,2)+1;
            fix(f,2) = sac(f+1,1)-1;

            % catch special case: if Engbert algorithms are applied 
            % without any saccade clustering [>> mergesacc()] there can
            % be back-to-back saccades with intervening "fixations" of
            % zero sample duration. Catch this by setting the duration
            % of these fixations to one sample
            if fix(f,1) > fix(f,2)
                fix(f,1) = fix(f,2); % 1-sample fixation
                clusterwarning = true;
            end
            
        end
    end
    % if epoch does not begin with saccade, add first fixation
    if sac(1,1) > 1
        fix = [[1 sac(1,1)-1]; fix];
    end
    % if epoch does not end with saccade, add last fixation
    if sac(end,2) < nsample
        fix = [fix;[sac(end,2)+1 nsample]];
    end
    
    % add more fixation properties
    for f = 1:size(fix,1)
        
        % fixation duration (in samples!)
        fix(f,3) = fix(f,2)-fix(f,1)+1;
        
        % mean fix. position: left eye
        if ldata
            fix(f,4) = mean( l(fix(f,1):fix(f,2),1) );
            fix(f,5) = mean( l(fix(f,1):fix(f,2),2) );
        else
            fix(f,4) = NaN;
            fix(f,5) = NaN;
        end
        % mean fix. position: right eye
        if rdata
            fix(f,6) = mean( r(fix(f,1):fix(f,2),1) );
            fix(f,7) = mean( r(fix(f,1):fix(f,2),2) );
        else
            fix(f,6) = NaN;
            fix(f,7) = NaN;
        end
        % binocular fixation position
        fix(f,8) = nanmean(fix(f,[4 6]));
        fix(f,9) = nanmean(fix(f,[5 7]));
    end
    fix(:,10) = e; % add index of corresp. data epoch
    
    % recompute latencies of eye movement events 
    % (only necessary for epoched data)
    offset = (e-1)*nsample;
    sac(:,[1 2]) = sac(:,[1 2])+offset;
    % special case: single sacc. lasts entire epoch
    if ~isempty(fix) 
        fix(:,[1 2]) = fix(:,[1 2])+offset;
    end
end

% columns of [fix]:
% 1: fixation onset (sample)
% 2: fixation offset (sample)
% 3: duration (samples)
% 4: mean fix position (L eye X)
% 5: mean fix position (L eye Y)
% 6: mean fix position (R eye X)
% 7: mean fix position (R eye Y)
% 8: mean fix position (L/R average X)
% 9: mean fix position (L/R average Y)
%10: index of corresponding data epoch (1 in case of contin. data)

%% remove fixations overlapping with "bad_ET" intervals
% this includes fixations interrupted by a blink (!)
% development note: loop is slow, need to implement more efficiently
if ~isempty(fix) & any(badvector(:))
    badfix = false(size(fix,1),1); % preallocate logical
    
    % go tru fixations, check whether "bad"
    for p = 1:size(fix,1)
        fixsmp = fix(p,1):fix(p,2);
        if any(ismember(fixsmp,badETsmp))
            badfix(p) = true;
        end
    end
    % remove fixations overlapping with "bad_ET" intervals
    fix(badfix,:) = [];
    %         if any(badfix)
    %             fprintf('\n-- Removed %i fixations that overlapped with "bad_ET" intervals',sum(badfix))
    %         end
end
  
% slow, but simple:
allsac = [allsac;sac];
allfix = [allfix;fix];

end % epoch loop

%% remove artificial eye movements caused by boundaries (data breaks)
% Remove all EMs whose onset is detected in temporal proximity to boundary.
% Otherwise, data breaks will likely result in additional fake sacc./fix.
% Applies only if eye movements are detected in the continuous data
if nepochs == 1

ix_bnd = find(cellfun(@(x) strcmp(x,'boundary'),{EEG.event.type})); % bug fix: now robust against numeric types

% any data breaks?
% (due to manual rejections or function pop_rej_eyecontin)
if ~isempty(ix_bnd)
    
    % minimum distance from data break in milliseconds (hard-coded)
    BOUNDDIST_MS = 50;
    BOUNDDIST    = round(BOUNDDIST_MS / (1000/EEG.srate)); 
    % data break latencies
    bound_lats   = round([EEG.event(ix_bnd).latency]); 

    % mark all samples close to data break
    boundvector  = zeros(1,EEG.pnts);
    for b = 1:length(bound_lats)
        lowr = bound_lats(b)-BOUNDDIST;
        uppr = bound_lats(b)+BOUNDDIST;
        if lowr <= 0, lowr = 1; end
        if uppr > EEG.pnts, uppr = EEG.pnts; end
        boundvector(lowr:uppr) = 1;
    end
    nearboundsmp = find(boundvector);
    
    % option 1: event onset is close to boundary
    % fakesac = find(ismember(allsac(:,1),nearboundsmp));
    % fakefix = find(ismember(allfix(:,1),nearboundsmp));
    
    % option 2: event on- or offset is close to boundary
    fakesac = find(ismember(allsac(:,1),nearboundsmp) | ismember(allsac(:,2),nearboundsmp));
    fakefix = find(ismember(allfix(:,1),nearboundsmp) | ismember(allfix(:,2),nearboundsmp));

    allsac(fakesac,:) = [];
    allfix(fakefix,:) = [];
    fprintf('\n--------------------------------------------------------------------');
    fprintf('\nFound %i data breaks (boundary events) in the continuous data',length(ix_bnd));
    fprintf('\nRemoving eye movements that might be artifacts of data breaks:');
    fprintf('\nRemoved %i saccades  < %i ms away from a boundary',length(fakesac),BOUNDDIST_MS);
    fprintf('\nRemoved %i fixations < %i ms away from a boundary',length(fakefix),BOUNDDIST_MS);
end

end

if clusterwarning && clustermode == 1
fprintf('\n--------------------------------------------------------------------');
fprintf('\n\n*** WARNING! ***\nDetected pairs of saccades that ended/started on adjacent data samples.');
fprintf('\nYou should probably use saccade clustering! (clustermode: 2,3,4)');
end

%% user feedback if bad eye-tracking data (values <= 0) but no bad_ET events
if sum(badepochs)>0 & isempty(ix_badETevent)
warning('\nI found bad or missing data (i.e. values <= 0) in the ET channels but no "bad_ET" events!');
fprintf('\nDetails:');
fprintf('\n-- %i of %i epochs contained gaze position values <= 0',sum(badepochs),nepochs);
fprintf('\n-- Total number of bad samples <= 0: %i',nbadsmp);
fprintf('\n-- Did you detect or reject bad intervals with out-of-range values (''Reject data based on eyetrack'')?');
fprintf('\n-- Did you subtract a baseline from eye channels (explaining negative values)?');
fprintf('\n-- Note that unmarked blinks or missing data in the eye-track will...:')
fprintf('\n ...distort the velocity threshold for saccade detection!');
fprintf('\n ...be erroneously detected as additional saccades/fixations!');
end

%% user feedback: saccade & fixation detection
fprintf('\n--------------------------------------------------------------------');
fprintf('\nVelocity thresholds used:'); if nepochs > 1, fprintf(' (mean across epochs):'); end;
if ldata, fprintf('\n\tLeft eye. Horiz.: %.2f %s/s. Vert.: %.2f %s/s',mean(l_msdx(e)vfacdegperpixel),metric,mean(l_msdyvfacdegperpixel),metric); end
if rdata, fprintf('\n\tRight eye. Horiz.: %.2f %s/s. Vert.: %.2f %s/s',mean(r_msdx(e)vfacdegperpixel),metric,mean(r_msdyvfacdegperpixel),metric); end
fprintf('\n--------------------------------------------------------------------')
if ~isempty(allsac)
fprintf('\n%i saccades detected:',size(allsac,1));
fprintf('\n\tMedian amplitude: %.2f %s',median(allsac(:,6)),metric);
fprintf('\n\tMedian duration: %.2f ms',median(allsac(:,3))*1000/EEG.srate);
fprintf('\n\tMedian peak veloc. %.2f %s/s',median(allsac(:,5)),metric);
end
if ~isempty(allfix)
fprintf('\n%i fixations detected:',size(allfix,1));
fprintf('\n\tMedian duration: %.2f ms',median(allfix(:,3))*1000/EEG.srate);
if ldata,fprintf('\n\tMedian fix. pos. left eye: Horiz.: %.2f px. Vert.: %.2f px',median(allfix(:,4)),median(allfix(:,5))); end
if rdata,fprintf('\n\tMedian fix. pos. right eye: Horiz.: %.2f px. Vert.: %.2f px',median(allfix(:,6)),median(allfix(:,7))); end
end

%% plot figure with eye movements properties
if plotfig
if ldata && rdata
fprintf('\n--------------------------------------------------------------------')
fprintf('\nPlotting eye movement properties...\nFor plotting, fixation locations are averaged across eyes.')
ploteyemovements(allsac(:,6),allsac(:,5),allsac(:,7),allfix(:,3)*1000/EEG.srate,mean(allfix(:,[4 6]),2),mean(allfix(:,[5 7]),2),metric);
elseif ldata
fprintf('\n--------------------------------------------------------------------')
fprintf('\nPlotting eye movement properties of left eye...')
ploteyemovements(allsac(:,6),allsac(:,5),allsac(:,7),allfix(:,3)*1000/EEG.srate,allfix(:,4),allfix(:,5),metric);
else
fprintf('\n--------------------------------------------------------------------')
fprintf('\nPlotting eye movement properties of right eye...')
ploteyemovements(allsac(:,6),allsac(:,5),allsac(:,7),allfix(:,3)*1000/EEG.srate,allfix(:,6),allfix(:,7),metric);
end
end

%% write eye movements to EEG.event & EEG.urevent
if writefix || writesac

% check: are there already eye movements events in EEG.event?
em_types = {'saccade','fixation','L_saccade','R_saccade','L_fixation','R_fixation'};
if any(cellfun(@(x) any(strcmp(x,em_types)),{EEG.event.type})) % updated [v.0.337]
    fprintf('\n--------------------------------------------------------------------')
    fprintf('\n*** WARNING! ***:\nFound existing eye movement events in EEG.event!');
    fprintf('\nNew events will be added to the already existing events!');
end

% write saccades
if writesac
    fprintf('\n--------------------------------------------------------------------')
    fprintf('\nAdding %i saccades to EEG.event...\n',size(allsac,1))
    EEG = addevents(EEG,allsac(:,[1 3 5 6 7 10 11:14]),{'latency','duration','sac_vmax','sac_amplitude','sac_angle','epoch','sac_startpos_x','sac_startpos_y','sac_endpos_x','sac_endpos_y'},'saccade');
end

% write fixations
if writefix
    fprintf('\n--------------------------------------------------------------------')
    fprintf('\nAdding %i fixations to EEG.event...\n',size(allfix,1))
    if ldata && rdata % binocular
        EEG = addevents(EEG,allfix(:,[1 3 8 9 10]),{'latency','duration','fix_avgpos_x','fix_avgpos_y','epoch'},'fixation');
        %EEG = addevents(EEG,allfix(:,[1 3 4 5 6 7 10]),{'latency','duration','fixposition_ly','fixposition_lx','fixposition_ry','fixposition_rx','epoch'},'fixation');
    elseif ldata % left eye recorded
        EEG = addevents(EEG,allfix(:,[1 3 4 5 10]),{'latency','duration','fix_avgpos_x','fix_avgpos_y','epoch'},'fixation');        
    elseif rdata % right eye recorded
        EEG = addevents(EEG,allfix(:,[1 3 6 7 10]),{'latency','duration','fix_avgpos_x','fix_avgpos_y','epoch'},'fixation');
    end
end
fprintf('--------------------------------------------------------------------')

end
fprintf('\nDone.\n\n')
end

Develop ways to summarize and compare trajectories

Should be able to pass a feather object to mcmc.fit
This will continue fitting using the same options as in the existing feather. In fact, it could look like

mcmc.fit(feather,target_f,target_P)

which would automatically reuse the model, npts_per_fit, and the previous ICs and prior

If fit_cluster decides to skip the fit, the old value for chisq is copied over. In the existing data sets, the data has been post-processed such that the AIC etc. is NaN and chisq is computed again based on the fixed parameters (see converter.m in a5ef029 in romesh-code). This needs to be integrated into fit_cluster

If a parallel chain is run, the diagnostic output should also be generated for the final chain, not just the workers. Otherwise it cannot be directly compared to the output from the single threaded routine

There should be some way to mark data as bad manually. This would mean storing the fit/don't fit criteria within the input file (rather than computing within fit_cluster.m)

Analysis routines expect model.weights but it is not a default property

Will need significant changes to support multiple spectra etc.

Maybe mcmc.data.load_subject or something descriptive like that

First time running the code and when trying

run fit_kaggle.m

it throws up error

Undefined function or variable 'matlabpool_cluster'.

some cursory googling suggests to me this isn't a built-in matlab function. Perhaps it's a custom utils function not contained in this repo??