function [ parameters ] = getswjparameters( samples, isInTrial, blinks, samplerate )
%-------------------------------------------------------------------
%
%  FUNCTION getswjparameters.m
%
%   From a set of eye movement data it calculates the necessary parameters
%   to calculate the SWJ index in the swjs.m function
%
%   Distinctive Features of Saccadic Intrusions and Microsaccades in
%   Progressive Supranuclear Palsy. Otero-Millan, Leigh, Serra, Troncoso,
%   Macknik, Martinez-Conde. Journal of Neuroscience (under review).
%
%   Input:
%   
%   All the parameters are cell arrays. Each element of the cell
%   corresponds to one individual subject
%   
%   samples{i}(:,1)     timestamps of the recording
%   samples{i}(:,2)     horizontal position of the left eye   
%   samples{i}(:,3)     vertical position of the left eye  
%   samples{i}(:,4)     horizontal position of the right eye     
%   samples{i}(:,5)     vertical position of the right eye  
%   isInTrial{i}        binary vector indicating for each sample if it
%                       belons to a trial (1) or not (0)
%   blinks{i}           binary vector indicating for each sample if it
%                       belons to a blink (1) or not (0)
%   samplerate           samplerate of the recordings
%
% Output:
%  paramters            parameters necessary to calculate the swj index
%     parameters.DM     Difference in direction mixed gaussian means [270 270]
%     parameters.DS     Difference in direction mixed gaussian stds [30 7]
%     parameters.DR     Difference in direction mixed gaussian ratios [0.4 0.6]
%     
%     parameters.RM     Magnitude similarity index mixed gaussian means [0 0]
%     parameters.RS     Magnitude similarity index mixed gaussian stds [0.39 0.16]
%     parameters.RR     Magnitude similarity index mixed gaussian ratios [0.4 0.6]
%     parameters.ISIP   Intersaccadic interval ex-gaussian mean, std and tau [125 60 180]
%     parameters.swjindexth     Threshold for the swj index [0.0014]


swj_mag1 = [];
swj_mag2 = [];
swj_dir1 = [];
swj_dir2 = [];
swj_isi  = [];

allsamples = samples;
for isubj=1:length(samples)
    [lsac, rsac] = saccades( allsamples{isubj}, samplerate, [], [], [], [] , isInTrial{isubj}, blinks{isubj} );
    
    samples = allsamples{isubj};
    if ( ~isempty (lsac) )
        mag = zeros(size(lsac,1),1);
        for i = 1:size(lsac,1)
            idx = lsac(i,1):lsac(i,2);
            dX = max(samples(idx,2)) - min(samples(idx,2));
            dY = max(samples(idx,3)) - min(samples(idx,3));
            mag(i)  = sqrt( dX.^2 + dY.^2);
        end
        dx = samples(lsac(:,2),2) - samples(lsac(:,1),2);
        dy = samples(lsac(:,2),3) - samples(lsac(:,1),3);
        directions	= atan2( dx, dy );
        dir = mod( (directions * (180 / pi)) + 360, 360);
        pairsidx = find(cos( deg2rad(dir(1:end-1))- deg2rad(dir(2:end))) > 0) ;
        mag1 = mag(pairsidx);
        mag2 = mag(pairsidx+1);
        dir1 = dir(pairsidx);
        dir2 = dir(pairsidx+1);
        isi = (lsac(pairsidx+1,1) - lsac(pairsidx,2))*1000/samplerate;
        swj_mag1 = cat(1,swj_mag1,mag1);
        swj_mag2 = cat(1,swj_mag2,mag2);
        swj_dir1 = cat(1,swj_dir1,dir1);
        swj_dir2 = cat(1,swj_dir2,dir2);
        swj_isi  = cat(1,swj_isi,isi);
    end
    if ( ~isempty ( rsac) )
        mag = zeros(size(rsac,1),1);
        for i = 1:size(rsac,1)
            idx = rsac(i,1):rsac(i,2);
            dX = max(samples(idx,4)) - min(samples(idx,4));
            dY = max(samples(idx,5)) - min(samples(idx,5));
            mag(i)  = sqrt( dX.^2 + dY.^2);
        end
        dx = samples(rsac(:,2),4) - samples(rsac(:,1),4);
        dy = samples(rsac(:,2),5) - samples(rsac(:,1),5);
        directions	= atan2( dx, dy );
        dir = mod( (directions * (180 / pi)) + 360, 360);
        pairsidx = find(cos( deg2rad(dir(1:end-1))- deg2rad(dir(2:end))) > 0 );
        mag1 = mag(pairsidx);
        mag2 = mag(pairsidx+1);
        dir1 = dir(pairsidx);
        dir2 = dir(pairsidx+1);
        isi = (rsac(pairsidx+1,1) - rsac(pairsidx,2))*1000/samplerate;
        swj_mag1 = cat(1,swj_mag1,mag1);
        swj_mag2 = cat(1,swj_mag2,mag2);
        swj_dir1 = cat(1,swj_dir1,dir1);
        swj_dir2 = cat(1,swj_dir2,dir2);
        swj_isi  = cat(1,swj_isi,isi);
    end
    
end

[ parameters ] = calculateswjparameters( swj_mag1, swj_mag2, swj_dir1, swj_dir2, swj_isi );

end
    
function [ parameters ] =  calculateswjparameters( swj_mag1, swj_mag2, swj_dir1, swj_dir2, swj_isi )
relmag = (swj_mag2 - swj_mag1) ./ (swj_mag2 + swj_mag1);
dirdif = mod(acos(cos((swj_dir1 - swj_dir2)*pi/180))*180/pi .* sign((swj_dir1 - swj_dir2))+90,360);
isi = swj_isi;

[relmagfitParams] = fitdata( relmag, 'Mix-2-Gaussian'  );
[dirdiffitParams] = fitdata( dirdif, 'Mix-2-Gaussian'  );
[isifitParams] = fitdata( isi, 'Ex-gaussian' );


isi = swj_isi;

parameters.DM = dirdiffitParams([1 2]);
parameters.DS = dirdiffitParams([3 4]);
parameters.DR = dirdiffitParams([5 6]);

parameters.RM = relmagfitParams([1 2]);
parameters.RS = relmagfitParams([3 4]);
parameters.RR = relmagfitParams([5 6]);

parameters.ISIP = isifitParams;

f1 = 1-normcdf( abs(dirdif-parameters.DM(1)), 0, parameters.DS(1))*parameters.DR(1)-normcdf( abs(dirdif-parameters.DM(2)), 0, parameters.DS(2))*parameters.DR(2);
f2 = 1-normcdf( abs(relmag-parameters.RM(1)), 0, parameters.RS(1))*parameters.RR(1)-normcdf( abs(relmag-parameters.RM(2)), 0, parameters.RS(2))*parameters.RR(2);
f3 = ( double(isi>=200).*(1-swjexgausscdf( isi,parameters.ISIP))+double(isi<200).*(swjexgausscdf( isi, parameters.ISIP)));

pswj = f1.*f2.*f3;


[i] = kmeans(log(pswj), 2);
parameters.swjindexth = max(min(pswj(i==1)),min(pswj(i==2)));

end

function [fitParams] = fitdata( data, distribution )



%-- fitting options ---------------------------------------------------
switch(distribution)
    case 'Ex-gaussian'
        data = data(data<1e3 & data>0);
        [phat] = fminsearch('swjMLE', [ 150,36,100],optimset('MaxFunEvals',500),'swjexgausspdf',data);
        fitParams = phat;
    case 'Mix-2-Gaussian'
        [mu, s, t] = fit_mix_gaussian(data,2);
        fitParams = [mu s t];
end
end