function [dE] = FNN(dataStream, tau, varargin);

%% [dE] = FNN(dataStream, tau);
%%
%%  Examples  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% [dE] = FNN(dataStream, tau, 'scan', 14, 'v');
%% [dE] = FNN(dataStream, tau, 'thresh', 15, 'order2', 4);
%% [dE] = FNN(dataStream, tau, 'plot', 'cutoff', 0.001);
%%  
%%  Input and output arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%  dataStream  (matrix )  N x 1 vector of time series data
%%  tau         (scalar )  estimated delay coordinate
%%
%%  dE          (scalar )  embedding dimension, dE
%%
%%  [argID,     (string )  See below. The values which are unambiguous can 
%%   value]     (varies )  be given without the preceeding argID.
%%
%% Valid argument IDs and corresponding values. Values which are
%% unambiguous (marked with '*') can be given without the preceeding argID.
%%
%%  'thresh'    (scalar )  threshold for Falseness calculation
%%  'order2'    (scalar )  threshold for high-dimensional correction
%%  'cutoff'    (scalar )  threshold for minimizing d_E
%%  'scan'      (scalar )  scan range of d_E for %FNNs
%%  'plot'      (*string)  draw %FNN(d_E) plot
%%  'v'         (*string)  verbose output; prints progress
%%
%%  Authors  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Erin Miller   (emiller@csss.santafe.edu)
%%  Nicky Grigg   (grigg@csss.santafe.edu)
%%  Aaron Clauset (clauset@csss.santafe.edu)
%%
%%  Santa Fe Institute, Complex Systems Summer School 2003
%%  date created : 06.30.03
%%  last modified: 07.01.03
%%  MATLAB file (not necessarily OCTAVE compatible)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p.plot     = 0;
p.thresh   = 10;
p.order2   = 2;
p.dECutoff = 0.001;
p.scan     = 0;
p.dEmax    = 14;
p.v        = 0;

% parse command-line parameters
i=1; 
while i<=length(varargin), 
  argok = 1; 
  if ischar(varargin{i}), 
    switch varargin{i},
     % FNN parameters
     case 'thresh',        p.thresh      = varargin{i+1}; i = i + 1;
     case 'order2',        p.order2      = varargin{i+1}; i = i + 1;
     case 'plot',          p.plot        = 1;
     case 'cutoff',        p.cutoff      = varargin{i+1}; i = i + 1;
     case 'scan',          p.dEmax       = varargin{i+1}; i = i + 1; p.scan = 1;
     case 'v',             p.v           = 1;
     otherwise argok=0; 
    end
  end
  if ~argok, 
    disp(['(FNNs) Ignoring invalid argument #' num2str(i+1)]); 
  end
  i = i+1; 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% check dimensionality of input matrix
[N, R] = size(dataStream);
if (R>1) dataStream = transpose(dataStream); end;

%% initial values for dE search
dE           = 1;
FalseNNs(dE) = 1.00;

%% search over range of dE for percent FNN
switch (p.scan)
 %% simple search for when FalseNNs(dE) < p.dECutoff   
 case 0
  while (FalseNNs(dE)>p.dECutoff)
   dE = dE + 1;
   n = floor(length(dataStream)/tau)*tau;
   M = zeros(n-(dE-1)*tau,dE);
   for k=1:dE
    M(:,k) = dataStream((k-1)*tau+1:n-(dE-k)*tau);
   end
 
   %% calculate Falseness for each pair of points
   warning off; t0 = clock;
   D        = zeros(size(M,1));
   for i=1:size(M,1)
        D(i,i+1:end) = transpose(abs(repmat(M(i,end),size(M,1)-i,1) - M(i+1:end,end)));
        D(i,i+1:end) = p.thresh<D(i,i+1:end)./ ...
            transpose(sqrt(sum((repmat(M(i,1:end-1),size(M,1)-i,1) - M(i+1:end,1:end-1)).^2,2))) & ...
            p.order2>D(i,i+1:end)./std(dataStream);
   end
 
   FalseNNs(dE) = sum(sum(D))./(prod(size(D))/2);
   if (p.v==1) fprintf('dE: %2.0f\t%%FNN: %10.8f\ttime: %6.4fs\n',dE, FalseNNs(dE), etime(clock,t0)); end;
  end;
 case 1
 %% scan full range of 1:p.dEmax, then find dE_opt after  
  for dE=2:p.dEmax
   dE = dE + 1;
   n = floor(length(dataStream)/tau)*tau;
   M = zeros(n-(dE-1)*tau,dE);
   for k=1:dE
    M(:,k) = dataStream((k-1)*tau+1:n-(dE-k)*tau);
   end
 
   %% calculate Falseness for each pair of points
   warning off; t0 = clock;
   D        = zeros(size(M,1));
   for i=1:size(M,1)
        D(i,i+1:end) = transpose(abs(repmat(M(i,end),size(M,1)-i,1) - M(i+1:end,end)));
        D(i,i+1:end) = p.thresh<D(i,i+1:end)./ ...
            transpose(sqrt(sum((repmat(M(i,1:end-1),size(M,1)-i,1) - M(i+1:end,1:end-1)).^2,2))) & ...
            p.order2>D(i,i+1:end)./std(dataStream);
   end
 
   FalseNNs(dE) = sum(sum(D))./(prod(size(D))/2);
   if (p.v==1) fprintf('dE: %2.0f\t%%FNN: %10.8f\ttime: %6.4fs\n',dE, FalseNNs(dE), etime(clock,t0)); end;
  end;
  dE = min(find(ceil(FalseNNs-p.dECutoff)==0))
end;

if (p.plot>0)
 plot(FalseNNs, 'bo:'); xlabel('Dimension, d_E'); ylabel('Percentage FNN'); axis([1 dE_max 0 1]);
 title(strcat('False Nearest Neighbors; d_E = ', int2str(dE)));
end;
