function [phasedrift] = phasefunc(A, B, tau)
%
% Computes phase function at the specified lag (tau) and returns it.
%
% Usage:
%       phasefunc(A, B, tau)
%
% where A is row matrix holding spike/discrete time events,
%       B is another such row matrix (A itself for finding auto phase function),
%       tau is the lag at which the phase function is sought.
%
% Reference:
%       A Phase Function to Quantify Serial Dependence between Discrete Samples.
%       Ramana Dodla and Charles J. Wilson.
%       Biophys J. 2010 February 17; 98(4): L5-L7.
%       http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2820633/
%
% Example: The following code plots a curve in Fig. 2b of the above reference.
%
%    A = [0, 0.18595, 0.33955, 0.4997, 0.61495, 0.7413, 0.9024]
%    p = []; 
%    for i = 0:400
%      p=[p phasefunc(A, A, i * 0.001)]; 
%    end 
%    x=[0:1:400] * 0.001; 
%    plot(x, p, 'o')
%
% Author: Ramana Dodla, 2009
%
   [c] = phasedisp_c(A, B + tau); 
   thetadrift = mean(abs(c(:, 1)));
   norm = sqrt(2.0);
   phasedrift = thetadrift / norm;
end
%------- Compute c_i: -----------
function [c] = phasedisp_c(A, B)
   pairs = phasedisp_pairs(A, B);
   mg = mean(pairs(:, 1));
   md = mean(pairs(:, 2));
   f1 = mg * mg + md * md;
   f2 = sqrt(f1);
   c = (f1 - pairs(:, 1) * mg - pairs(:, 2) * md) / f2;
end
%------ Compute (gamma_i, delta_i): --
function [pairs] = phasedisp_pairs(A, B)
  i = 2; 
  j = 2; 
  pairs = [];
  T = max(A(1), B(1));
  while i <= length(A) && j <= length(B)
     if A(i) >= B(j)
        if A(i-1) >= B(j)
           T = A(i-1); 
	   j = j + 1;
        else
           g = (B(j) - T) / (A(i) - A(i-1));
           d = (B(j) - T) / (B(j) - B(j-1));
           pairs = [pairs; [g d]];
           T = B(j); 
	   j = j + 1;
        end;
     else
        if B(j-1) >= A(i)
           T = B(j-1);
           i = i+1;
        else
           g = (A(i) - T)/(A(i) - A(i-1));
           d = (A(i) - T)/(B(j) - B(j-1));
           pairs = [pairs; [g d]];
           T = A(i);
           i = i + 1;
        end;
     end
  end
end
