function S = lyap_k(x, m, tau, r, delta_n)
%
% An implementation of the lyap_k function in the TISEAN
% package.
%
% arguments
%  x: observed time data
%  m: embedding dimensions
%  tau: time delay for embedding
%  r: radius of neighborhood
%  delta_n: time span
%

% length of the original observed data
%
lenx = length(x);

% we must first estimate the state space using embedding 
%
X = delay_reconstruct(x, tau, m);

% find neighbors
%
neighbors = [];
   
totalavg=0;

for i=1:length(X(:,1))
    for j=1:length(X(:,1))
    
        if (j < i-(m-1)*tau | j > i+(m-1)*tau)
            
            dim = [];
        
            for k=1:length(X(1,:))
                 dim = [dim, (X(i,k)-X(j,k))^2]; 
            end
        
            sumdim = sum(dim);
            distance = sumdim^(1/2);
            
            if (distance < r)
                neighbors = [neighbors, j];
            end
        end
    end

    % convert the indices to match with observed time series.
    % the current neighbor indices are for the state space
    %
    new_neighbors = [];

    for j=1:length(neighbors)        

        temp = neighbors(j) + ((m-1)*tau);     
    
        % make sure we didn't pick up neighbors too close to our
        % reference state
        %
        if (temp+delta_n < lenx)
            new_neighbors = [new_neighbors, temp];
        end
    end
    
    neighbors = new_neighbors;
    if (length(neighbors) > 0)
        avg = mean(abs(x(i+((m-1)*tau))-x([neighbors+delta_n])));
        totalavg = totalavg + log(avg); 
    end
      
end

S = totalavg/length(X);
