function [divergenceMean,divergenceSd] = computedivergence(data,nNearest,nMaxDt,averagingMethod,phase)
% COMPUTEDIVERGENCE Calculates the divegence of nearest neighbors.
%   COMPUTEDIVERGENCE(DATA,NNEAREST,NMAXDT) returns the average normalized 
%   expansion distance of NNEAREST neighbors from time t to time t+dt*NMAXDT.
%   The expansion distances are normalized by the initial distance between the
%   nearest neighbors.
%
%   COMPUTEDIVERGENCE(DATA,NNEAREST,NMAXDT,AVERAGINGMETHOD) where
%   AVERAGINGMETHOD is either 'standard' (default) or 'rosenstein'. The
%   'standard' is as described above. For 'rosenstein', the natural log of the 
%   expansion distance is taken before averaging. After averaging, the
%   exponential is used to transform the units back to divergence. Finally, the 
%   result is normalized by the initial divergence to equate to the 'standard' 
%   method. This was proposed by Rosenstein to be more robust to differences in 
%   initial separation, but depends on the assumption that expansion is 
%   exponential with divergence time.
%
%   NMAXDT is an integer of data samples forward expansion should be tracked.
%
%   [DIVERGENCEMEAN,DIVERGENCESD] = COMPUTEDIVERGENCE(...) returns the standard
%   deviation with respect to divergence time.
% 
%   Input
%       DATA            | m x n | double | m is # samples and n is # states
%       NNEAREST        | 1 x 1 | double | integer # of nearest neighbors
%       NMAXDT          | 1 x 1 | double | integer # of samples to follow
%       AVERAGINGMETHOD |       | string | optional input argument (see above) 
%
%   Output
%       DIVERGENCEMEAN | (nMaxDt+1) x 1 | double | 
%       DIVERGENCESD   | (nMaxDt+1) x 1 | double | 
%   
%   Example 
%      lorenzex; % Run the example script
%
%   Reference
%     [1] Rosenstein, M. T., Collins, J. J., & De Luca, C. J., 1993, "A 
%         practical method for calculating largest Lyapunov exponents from small 
%         data sets," Physica D: Nonlinear Phenomena, Vol. 65, pp. 117-134.
%
%   See also LORENZEX, COMPUTELYAPUNOVSPECTRUM.  

% Created by: Kevin P. Granata Musculoskeletal Biomechanics Laboratory
%             http://www.biomechanics.esm.vt.edu
% Last Revision: 19 October 2007

% Check that both nNearest and nMaxDt are integers
  if (fix(nNearest) ~= nNearest)
      error('Input argument nNearest should be an integer.');
  elseif (fix(nMaxDt) ~= nMaxDt)
      error('Input argument nMaxDt should be an integer.');
  end;

% Handle optional argument averagingMethod
  if ~exist('averagingMethod','var')
      averagingMethod = 'standard';
  end;
  
% Mean-zero data (Required for auto-correlation)
  data = data-repmat(mean(data,1),size(data,1),1);

% Determine the number of samples to skip to avoid sampling nearest neighbors on
% same cycle. To ensure nSkip is sufficiently large, skip +/- first minimum of 
% autoCorrelation, i.e. +/- nSkip.
  autoCorrelation  = xcorr(data(:,1));
  autoCorrelation  = autoCorrelation(ceil(length(autoCorrelation)/2):end);
  nSkip = find(diff(autoCorrelation)>0,1); % Find first min of autoCorrelation

% Find nearest neighbors
  nearestIndex = findnearestneighbors(data,nNearest,nSkip,nMaxDt);

% Track nearest neighbors
  stateDistance = tracknearestneighbors(data, nearestIndex, nMaxDt);

% Determine which averaging method to apply
  switch lower(averagingMethod)
      case 'standard'
          % Normalize stateDistance by initial stateDistance, excluding any time 
          % samples where initial stateDistance is zero for any neighbors.
          keepStates      = any(stateDistance(:,1,:) ~= 0,3);
          stateDistance   = stateDistance(keepStates,:,:);
          initialDistance = repmat(stateDistance(:,1,:),[1 nMaxDt+1 1]);
          stateDistance   = stateDistance./initialDistance;
          
          % Mean & SD of stateDistance./initialDistance over nearest neighbors & 
          % divergence time
          divergenceMean = mean(mean(stateDistance,3),1)';
          divergenceSd   = std(mean(stateDistance,3),[],1)';
          
      case 'rosenstein'
          % THIS APPROACH NEEDS TO BE CHECKED 6-19-2007, Jesse Norris
          
          % Rosenstein approach of taking the natural log before averaging to 
          % avoid problems with various inital distances. Excluding any time 
          % samples where stateDistance is zero for any of the nearest neighbors 
          % at any separation time.
          keepStates    = any(any(stateDistance(:,:,:) ~= 0,3),2);
          stateDistance = stateDistance(keepStates,:,:);
          
          % exp(mean(ln(stateDistance))) over nearest neighbors/divergence time
          divergenceMean = exp(mean(mean(log(stateDistance),3),1));        
          divergenceSd   = exp(std(mean(log(stateDistance),3),[],1));                
          
          % divide divergence mean & std by it's intial value to produce similar
          % normalization to the standard approach
          divergenceSd   = divergenceSd'/divergenceMean(1); 
          divergenceMean = divergenceMean'/divergenceMean(1);
          
      case 'phase'
          divergenceMean = zeros(nMaxDt+1,length(phase));
          divergenceSd   = zeros(nMaxDt+1,length(phase));
          
          for i = 1:length(phase)
              iPhase = phase{i}(phase{i} < (size(data,1)-nMaxDt));
              keepStates = any(any(stateDistance(iPhase,:,:) ~= 0,3),2);
              stateDistanceTemp = stateDistance(keepStates,:,:);              
              
              % exp(mean(ln(stateDistance))) over nearest neighbors/divergence time
              divergenceMean(:,i) = exp(mean(mean(log(stateDistanceTemp),3),1))';        
              divergenceSd(:,i)   = exp(std(mean(log(stateDistanceTemp),3),[],1))';                
              
              % divide divergence mean & std by it's intial value to produce similar
              % normalization to the standard approach
              divergenceSd(:,i)   = divergenceSd(:,i)/divergenceMean(1,i); 
              divergenceMean(:,i) = divergenceMean(:,i)/divergenceMean(1,i);
              
              clear stateDistanceTemp;
          end          
          
      otherwise
          error('Unknown input argument for ''aveargingMethod''');
  end;
  
end % END computedivergence