%%%============================================================
%%% Demo file for the Microsaccade Toolbox
%%% (R-language Version 0.9)
%%% Authors: Ralf Engbert, Petra Sinn, Konstantin Mergenthaler,
%%% and Hans A. Trukenbrod
%%% Date: February 20th, 2014
%%%
%%% Converted to Matlab by Anna-Katharina Hauperich, 20th June 2019
%%%============================================================
clearvars
close all
clc
%%------------------------------
%% 1. Detection of microsaccades
%%------------------------------
%% Set parameters
SAMPLING = 500;
MINDUR = 3;
VFAC = 5;

%% Read raw data (ASCII; subject 01, trial 005)
Part =1;
Trial = 5;
d = readtable(['./data/f0' num2str(Part) '.00' num2str(Trial) '.dat']);

%% Select epoch from trial, transform to matrix
idx = 3001:4500;
xl = table2array(d(idx,2:3));
xr = table2array(d(idx,4:5));

%% Apply running average
xls = smoothdataMSR(xl,1,2);
xrs = smoothdataMSR(xr,1,2);

%% Detection of microsaccades
msl = microsaccMSR(xl,VFAC,MINDUR,SAMPLING);
msr = microsaccMSR(xr,VFAC,MINDUR,SAMPLING);
sac = binsaccMSR(msl.table,msr.table);

%% (Binocular) microsaccades
N = sac.N(1);
bin = sac.bin;

%% Plot trajectory
set(0,'DefaultFigureColor',[1 1 1])
figure
subplot(1,2,1)
hold on
plot(xls(:,1),xls(:,2),'k')
MSXL = nan(size(xls(:,1)));
MSYL = nan(size(xls(:,2)));
for n = 1:N
    MSXL(bin(n,1):bin(n,2)) = xls(bin(n,1):bin(n,2),1);
    MSYL(bin(n,1):bin(n,2)) = xls(bin(n,1):bin(n,2),2);
    
end
plot(MSXL,MSYL,'r-')
title('Position')
xlabel('Horizontal position [deg]')
ylabel('Vertical position [deg]')
daspect([1 1 1])

%% Plot trajectory = 2D velocity space
subplot(1,2,2)
hold on
vls = vecvelMSR(xl,SAMPLING,2);
plot(vls(:,1),vls(:,2),'k')

MSVelXL = nan(size(xls(:,1)));
MSVelYL = nan(size(xls(:,2)));
for n = 1:N
    MSVelXL(bin(n,1):bin(n,2)) = vls(bin(n,1):bin(n,2),1);
    MSVelYL(bin(n,1):bin(n,2)) = vls(bin(n,1):bin(n,2),2);
end
plot(MSVelXL,MSVelYL,'r-')

phi = linspace(0,2*pi,300);
cx = msl.radius(1)*cos(phi);
cy = msl.radius(2)*sin(phi);
plot(cx,cy,'k:')
title('Velocity')
xlabel('Horizontal velocity [deg/sec]')
ylabel('Verticalvelocity [deg/sec]')
daspect([1 1 1])
print(['TrajectoryAndVel_P' num2str(Part) 'T' num2str(Trial) '_Matlab.png'], '-dpng', '-r400')

%%------------------------------------------
%% 2. Plot main sequence for microsaccades
%%------------------------------------------
%% Set parameters
SAMPLING = 500;
MINDUR = 3;
VFAC = 5.5;
ctab = [1 0 0;...
    0 0 1;...
    0 1 0; ...
    0 1 1;...
    1 0 1];
N = 0;
figure
hold on
count = 1;
for vp = 1:5
    for  trial = 1:5
        
        %% Read raw data (ASCII)
        filename = ['data/f0' num2str(vp) '.00' num2str(trial) '.dat'];
        d = readtable(filename);
        
        %% Select epoch from trial, transform to matrix
        xl = table2array(d(:,2:3));
        xr = table2array(d(:,4:5));
        
        %% Detection of microsaccades
        msl = microsaccMSR(xl,VFAC,MINDUR,SAMPLING);
        msr = microsaccMSR(xr,VFAC,MINDUR,SAMPLING);
        sac = binsaccMSR(msl.table,msr.table);
        tab = sacparMSR(sac);
        SizeTag = size(tab);
        N = N + SizeTag(1);
        plot(tab(:,8),tab(:,5),'o','Color', ctab(vp,:))
        ylabel('Peak velocity (deg/s)')
        xlabel('Amplitude (deg)')
        title('Main Sequence')
        xlim([0.02,2])
        ylim([5,200])
        set(gca, 'YScale', 'log')
        set(gca, 'XScale', 'log')
        count = count+1;
        if count>5
            count = 1;
        end
    end
end
print(['MainSequence_Matlab.png'], '-dpng', '-r400')

%%------------------------------------
%% 3. Generation of surrogate data
%%------------------------------------
%% Read raw data (ASCII; subject 01, trial 005)

d = readtable(['./data/f0' num2str(Part) '.00' num2str(Trial) '.dat']);
%% Select epoch from trial, transform to matrix
idx = 3001:4500;
xl = table2array(d(idx,2:3));
xr = table2array(d(idx,4:5));
%% Generation of surrogates
xls = smoothdataMSR(xl,1,500);
xlsur = surrogateMSR(xls,500);

%% Plot original and surrogate trajectories
figure
subplot(1,2,1)
hold on
plot(xls(:,1),xls(:,2),'k')
title('Original Trajectory')
xlabel('Horizontal position [deg]')
ylabel('Vertical position [deg]')
daspect([1 1 1])

subplot(1,2,2)
hold on
plot(xlsur(:,1),xlsur(:,2),'r')
title('Surrogate Trajectory')
xlabel('Horizontal position [deg]')
ylabel('Vertical position [deg]')
daspect([1 1 1])
print(['OriginalSurrogate_P' num2str(Part) 'T' num2str(Trial) '_Matlab.png'], '-dpng', '-r400')

%% Plot acf of original trajectory and surrogate tracetory

figure
subplot(2,1,1)
vls = vecvelMSR(xls,500,2);
autocorr(vls(:,1))
title('Original Data')
subplot(2,1,2)
hold on
for s = 1:20
    xlsur = surrogateMSR(xls,500);
    vlsur = vecvelMSR(xlsur,500,2);
    autocorr(vlsur(:,1))
end
title('Surrogate Data')

print(['OriginalSurrogateAutocorr_P' num2str(Part) 'T' num2str(Trial) '_Matlab.png'], '-dpng', '-r400')


%%-------------------------------------------------
%% 4. Surrogate analysis for different thresholds
%%-------------------------------------------------
%% range of detection thresholds
from=3;
to=8;
by=0.5;
vfac = from:by:to;

%% subjects and trials
vp = 5;
ntrials = 5;

%% read raw data
num = length(vfac);
ncol=3;
mstab = zeros(11,3);
for  v = 1:vp
    for  f = 1:ntrials
        disp(['... processing: vp =' num2str(v) ', trial =' num2str(f)]);
        sep='';
        filename = ['data/f0' num2str(v) '.00' num2str(f) '.dat'];
        d = readtable(filename);
        %% select epoch, transform to matrix
        xl = table2array(d(:,2:3));
        xr = table2array(d(:,4:5));
        Sizexr = size(xr);
        dur = Sizexr(1)/SAMPLING;
        for i = 1:num
            %% detect microsaccades
            VFAC=vfac(i);
            msl = microsaccMSR(xl,VFAC,3,500);
            msr = microsaccMSR(xr,VFAC,3,500);
            if ~isempty(msr.table) && ~isempty(msl.table)
                
                sac = binsaccMSR(msl.table,msr.table);
                N = sac.N(1)/dur;
            end
            %% computation of surrogate data
            xsl = surrogateMSR(xl,500);
            xsr = surrogateMSR(xr,500);
            msl = microsaccMSR(xsl,VFAC,3,500);
            msr = microsaccMSR(xsr,VFAC,3,500);
            if ~isempty(msr.table) && ~isempty(msl.table)
                sac = binsaccMSR(msl.table,msr.table);
                Nsur = sac.N(1)/dur;
                mstab(i,1) = vfac(i);
                mstab(i,2:3) = mstab(i,2:3) + [N,Nsur];
            end
        end
    end
end
mstab(:,2:3) = mstab(:,2:3)/(vp*ntrials);

%% Plot surrogate analysis
figure
hold on
plot(mstab(:,1),mstab(:,2),'-ko')
plot(mstab(:,1),mstab(:,3),'-ro')
plot(mstab(:,1),mstab(:,2)-mstab(:,3),'-bo')
title('Surrogate Analysis')
xlabel('Threshold multiplier \lambda')
ylabel('Rate [1/s]')
print(['SurrogateRate_P1toP' num2str(vp) '_Matlab.png'], '-dpng', '-r400')

%%-------------------------------------------------
%% 5. Random walk analysis
%%-------------------------------------------------

d = readtable(['./data/f0' num2str(Part) '.00' num2str(Trial) '.dat']);
xl = table2array(d(:,2:3));
xr = table2array(d(:,4:5));

%% Perform lagged distance analysis
rvl = lagdistMSR(xl);
vl = diff(xl);

%% Perform random shuffling
R1 = randperm(length(vl(:,1)));
R2 = randperm(length(vl(:,1)));

vls = [(vl(R1,1)),vl(R2,2)];
rxls = [cumsum(vls(:,1)),cumsum(vls(:,2))];
rvls = lagdistMSR(rxls);

%% Plot results
figure
hold on
plot(rvl(:,1),rvl(:,2),'k')
plot(rvls(:,1),rvls(:,2),'r')
title('Random Walk Analysis')
xlabel('lag')
ylabel('D(lag)')
set(gca, 'YScale', 'log')
set(gca, 'XScale', 'log')
print(['RandomWalk_P' num2str(Part) 'T' num2str(Trial) '_Matlab.png'], '-dpng', '-r400')

%%-------------------------------------------------
%% 6. Box-count analysis
%%-------------------------------------------------
dx = 0.01;
dt = 100;
msl =  microsaccMSR(xls,5,3,500);

%% Calculate box-count for drift epochs
xx = xls((msl.table(1,1)-dt):(msl.table(1,1)-1),:);
boxes = boxcountMSR(xx,dx);

%% Plot results
figure
hold on
plot(xx(:,1),xx(:,2),'k')
x_min = min(xx(:,1));
y_min = min(xx(:,2));
M = length(xx(:,1));

for   l = 1:M
    i = floor( (xx(l,1)-x_min)/dx ) + 1;
    j = floor( (xx(l,2)-y_min)/dx ) + 1;
    a = x_min+i*dx;
    b = y_min+j*dx;
    plot([a-dx,a],[b,b],'g')
    plot([a-dx,a],[b-dx,b-dx],'g')
    plot([a-dx,a-dx],[b-dx,b],'g')
    plot([a,a],[b-dx,b],'g')
end
print(['BoxCount_P' num2str(Part) 'T' num2str(Trial) '_Matlab.png'], '-dpng', '-r400')

%% EOF
license('inuse')
