function distance = statedistance(data,nSkip,nMaxDt)
% STATEDISTANCE Compute distance between sampled states.
%   STATEDISTANCE(DATA,NSKIP,NMAXDT) returns the euclidean distance between all 
%   states in the data set that are at least NSKIP samples apart and that may be 
%   tracked forward in time for NMAXDT samples. DATA should be formated such 
%   that the rows correspond to samples of the state and the columns are the 
%   states at that sample increment. 
%
%   NSKIP is the number of data samples to skip both forward and backward in 
%   time to avoid sampling nearest neighbors on the same cycle. 
%
%   NMAXDT is the number of time samples that nearest neighbors will be tracked
%   forward in time.
%
%   Note that this method does not work for data sets that are large, i.e. on
%   the order of magnitude of more than ~8000 time samples (m > ~8000) because 
%   the output requires too much memory. As an alternative, it might be
%   desirable to use the function FINDNEARESTNEIGHBORS
%
%   Input
%       DATA     | m x n | double | m is # of samples and n is # of states
%       NSKIP    | 1 x 1 | double | integer, # of samples to skip to avoid nn
%       NMAXDT   | 1 x 1 | double | integer, # samples will be tracked forward 
%
%   Output
%       DISTANCE | (m-nMaxDt-nSkip-1)*(m-nMaxDt-nSkip)/2 x 1 | double | [1]
%
%       [1] DISTANCE is a column vector that can be reconstructed to yield 
%           the statedistance matrix given in doc/statedistance.pdf
%
%   Example 
%
%   Reference
%
%   See also COMPUTELYAPUNOVSPECTRUM, DOC/STATEDISTANCE.PDF.  

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

  nSamples = size(data,1);    % Total # of time samples
  nValid   = nSamples-nMaxDt; % # of samples that may be tracked forward NMAXDT
  
% Find distance to every other state in the data set
  distance = zeros((nValid-nSkip-1)*(nValid-nSkip)/2,1); % Pre-allocate

  for i = 1:(nValid-nSkip-1)
      j = (i+nSkip+1):nValid; % Indicies we'd like to compare against
      vectorDifferences = ones(length(j),1)*data(i,:)-data(j,:);
      k = (i-1)*(nValid-nSkip-i/2)+(j-nSkip-i);
      distance(k) = sum(vectorDifferences.^2,2);      
  end;  
  
% Convert to Euclid distance
  distance = sqrt(distance); 
  
% Number of samples to skip for indexing k 
%  i                                                  = (i-1)*(nValid-nSkip)-(i-1)*i/2
%  i                                                  = (i-1)*(nValid-nSkip-i/2)
%  1 0                                                = 0*(nValid-nSkip)-0
%  2 nValid-nSkip-1                                   = 1*(nValid-nSkip)-1
%  3 nValid-nSkip-1 + nValid-nSkip-2                  = 2*(nValid-nSkip)-3
%  4 nValid-nSkip-1 + nValid-nSkip-2 + nValid-nSkip-3 = 3*(nValid-nSkip)-6

%{
% Old slower method
  for i = 1:(nValid-nSkip-1)
      for j = (i+nSkip+1):(nValid)
          vectorDifference = data(i,:)-data(j,:);
          k = (i-1)*(nValid-nSkip)-sum(1:(i-1))+(j-nSkip-i);
          distance(k) = vectorDifference*vectorDifference';
          %stateDistanceCheck(i,j) = distance(k);
      end;
  end;
  distance = sqrt(distance);
%}

end % END statedistance