%-------------------------------------------------------------------------
% Computational modeling of fixation durations in scene viewing
% by Anke Cajar, Jochen Laubrock, & Ralf Engbert (2013)
%-------------------------------------------------------------------------
% Please cite:
% "Control of fixation duration during scene viewing by interaction of 
%  foveal and peripheral processing" 
%  Laubrock, J., Cajar, A., & Engbert, R. (in press). Journal of Vision.
%-------------------------------------------------------------------------
clear all;
close all;

%% model parameters (see journal article for more information)
N(1) = 9;
N(2) = 33;
N(3) = 22;
prob = 0.66;
t_sac = 253;
j_C = 2.61;
k_C = 3.15;
j_fLP = 2.58;
j_fHP = 3.44;
k_pLP = 2.53;
k_pHP = 3.06;
rho = 4.37;

%% no. of simulated fixation durations
RUNS = 1; 
% RUNS = 1: plot of single-trial simulations (see Fig. 6)
% RUNS > 1, e.g. RUNS = 1000: report of fixation duration distributions and mean values (see Figs. 7 and 8)

%% initialize variables 
if RUNS==1
    maxcondition = 1;
else
    maxcondition = 5;
end
L = 1000;
t = zeros(1,L+1);
term = L;
RT = zeros(1,RUNS);
m = zeros(maxcondition,2);

%% plot settings
figure('units', 'normalized', 'position',[0 0 1 1]);
col = {'r-','b-','g-','b--','g--'};
LINEWIDTH = 1.8;
FONTSIZE = 24;
MARKERSIZE = 15;

%% simulate fixation durations

for condition = 1:maxcondition
    n = zeros(L+1,3);
    for j = 1:RUNS 
                        
        dir = [1 1 1];          
        Nmax = [0 0 0];
        
        rate(1) = N(1)/t_sac;
        rate(2) = j_C*rate(1);
        rate(3) = k_C*rate(1);
        
        switch condition
            case 2  % foveal low-pass
                rate(2) = j_fLP*rate(1);
            case 3  % foveal high-pass
                rate(2) = j_fHP*rate(1);
            case 4  % peripheral low-pass
                rate(3) = k_pLP*rate(1);
            case 5  % peripheral high-pass
                rate(3) = k_pHP*rate(1);
        end
        k = 1;
        while k<L
            
            % 0. Foveal inhibition
            diff = n(k,2)-n(k,3); 
            if diff<0
                diff = 0;
            end
            inhib = rate(1)/(1+rho*diff/N(2));           
            
            % Gillespie Algorithm/minimal process method (Gillespie, 1978)
            % 1. Time step
            w(1) = inhib;
            for s = 2:3
                if n(k,s)==0 && Nmax(s)==1 
                    w(s) = 0;
                else
                    w(s) = rate(s);
                end
            end
            W = sum(w);
            dt = -1/W*log(1-rand);
            t(k+1) = t(k) + dt;
            
            % 2. Choose transition (linear selection algorithm)
            p = w/W;
            c = p(1);
            trans = 1;
            rp = rand;
            while c<rp
                trans = trans + 1;
                c = c + p(trans);
            end
            
            n(k+1,:) = n(k,:);
            
            randdir = dir(trans);
            if trans>1 && randdir~=0
                rn = rand(1);
                if rn>prob
                    randdir = randdir*-1; 
                end
            end

            n(k+1,trans) = n(k+1,trans)+randdir;
                        
            % 3. Check random-walk thresholds and orientations
            for s = 2:3
                if n(k+1,s)>=N(s)
                    n(k+1,s) = N(s);
                    Nmax(s) = 1;
                    dir(s) = -1;
                end

                if n(k+1,s)<=0
                    n(k+1,s) = 0;
                    if Nmax(s)==1
                        dir(s) = 0;
                    end
                end
            end
            
            % 4. Terminate simulation if timer reaches threshold
            if n(k,1)==N(1)
                term = k;
                k = L;
            end
            k = k+1;
        end
        RT(j) = t(term);
    end
    
    if RUNS==1
        fprintf(1,'RT = %.1f ms\n',RT);
    else
        fprintf(1,'condition %d: RT = %.1f ms (%.1f ms)\n',condition,mean(RT),median(RT));
        m(condition,1:2) = [mean(RT) std(RT)/sqrt(RUNS)];
    end
    
    %% the following code is for plotting the simulation results only
    if RUNS==1
        % plot activations
        n = n(1:term,:);
        t = t(1:term);
        t = [t; t];
        t = t(:);
        n1 = n(:,1)';
        n1 = [0 n1(1:end-1); n1];
        n1 = n1(:);
        n2 = n(:,2)';
        n2 = [0 n2(1:end-1); n2];
        n2 = n2(:);
        n3 = n(:,3)';
        n3 = [0 n3(1:end-1); n3];
        n3 = n3(:);
        n = [n1 n2 n3];
        plot(t,n,'LineWidth',LINEWIDTH)
        ylim([0 max(n(:))])
        set(gca,'FontSize',FONTSIZE,'LineWidth',LINEWIDTH)
        xlabel('Time t [ms]');
        ylabel('Activation');
        legend('Timer activation a_T(t)','Foveal activation a_F(t)','Peripheral activation a_P(t)', 'Location', 'NorthWest');
    else
        % plot fixation duration distributions
        subplot(1,2,1);
        [yd, xi] = ksdensity(RT);
        plot(xi,yd,char(col(condition)),'LineWidth',LINEWIDTH);
        hold on
        if condition==maxcondition
            axis square;
            set(gca,'FontSize',FONTSIZE,'LineWidth',LINEWIDTH);
            xlabel('Fixation duration [ms]');
            ylabel('Probability');
            xlim([0 1000])
            legend('Control','Foveal low-pass','Foveal high-pass','Peripheral low-pass', 'Peripheral high-pass');
        % plot mean fixation durations
            subplot(1,2,2);
            errorbar(1,m(1,1),m(1,2),'ro-', 'LineWidth',LINEWIDTH, 'MarkerSize', MARKERSIZE, 'MarkerFaceColor', 'r');
            hold on;
            errorbar([3 2]',[m(2,1) m(4,1)],[m(2,2) m(4,2)],'bd-', 'LineWidth',LINEWIDTH, 'MarkerSize', MARKERSIZE, 'MarkerFaceColor', 'b');
            errorbar([3 2]',[m(3,1) m(5,1)],[m(3,2) m(5,2)],'gs-', 'LineWidth',LINEWIDTH, 'MarkerSize', MARKERSIZE, 'MarkerFaceColor', 'g');
            axis square;
            set(gca,'FontSize',FONTSIZE,'LineWidth',LINEWIDTH, 'XTick',[1 2 3],'XTickLabel',{'Control','Periphery','Fovea'});
            xlabel('Mask region');
            ylabel('Fixation duration [ms]');
            legend('Control','Low-pass','High-pass', 'Location', 'NorthWest');
        end
    end
end