%-----------------------------------------------------------------------
%
%   false_neighbers_Kennel.m
%                                  written by N. Chino on May 1, 2017.
%                                renamed by N. Chino on April 1, 2018.
% 
%   This is a renamed version of the program, false_neighbors_general3.m,
%  written on May 1, 2017.  Subprogram is falnei2.m.
%
%  This program determines an acceptable minimum embedding dimension by 
%  looking at the behavior of near neighbers under changes in embedding
%  dimension from d to d+1.  This method was proposed by Kennel et al.
%  (1992).  It is assumed, in this program, that the appropriate time lag
%  (delay increment) is known.  One method to estimate the delay increment
%  is to apply the minimum mutual information criterion (MMI) proposed by
%  Fraser and Swinney (1986).
%----------------------------------------------------------------------
%
%load data/chino/Henon-map.dat
%load data/chino/Henon-map10000.dat
%x=Henon_map10000';
load data/chino/Henon-map3000.dat
xorig=Henon_map3000';
x=xorig';
%Nx=size(x,2);
%
%load data/chino/Langford-torus.dat
%x=Langford_torus;
%
%load data/chino/logistic-map.dat
%x=logistic_map;
%
%load data/chino/lorenz-butterfly.dat
%x=lorenz_butterfly;
%
Nx=size(x,1);
disp('1. Size of the time series, Nx=');
Nx
% Specify the time lag as one of the embedding parameters.
%lag=1;
%lag=2; 
%lag=8;
lag=6;
%lag=10;
disp('2. Time lag, lag=');
lag
Rtol=10; Atol=2.0;
disp('3. Value of the first criterion, Rtol=');
Rtol
disp('4. Value of the second criterion, Atol=');
Atol
% Rtol=25; Atol=1.0;
%
maxdim=10;
for dim=1:maxdim
  Ny=Nx-(dim-1)*lag;
  for i=1:Ny
    j=1:dim;
    y(i,j)=x(i+(j-1)*lag);
  end;
  yrow=size(y,1);
%
  dmp1=dim+1;
  Nyp1=Nx-(dmp1-1)*lag;
  for i=1:Nyp1
    j=1:dmp1;
    yp1(i,j)=x(i+(j-1)*lag);
  end;
  yp1row=size(yp1,1);
%
  for i=1:Nyp1
    [rd1,rdp1,innp1]=falnei2(x,y,yrow,yp1,yp1row,i);
    Rdnn1(i)=rd1;
    Rdnnp1(i)=rdp1;
    Inn1p1(i)=innp1;
  end;
%
  nofn1(dim)=0;
  for i=1:Nyp1
   crit1(i)=abs(x(i+dim*lag)-x(Inn1p1(i)+dim*lag))/Rdnn1(i);
   if crit1(i)>Rtol nofn1(dim)=nofn1(dim)+1;
   end;
  end;
  %
  xy=x(1:Ny);
  Ra=std(xy); 
%
  nofn2(dim)=0;
  for i=1:Nyp1
   crit2(i)=Rdnnp1(i)/Ra;
   if crit2(i)>Atol nofn2(dim)=nofn2(dim)+1;
   end;
  end;
  %
  clear Ny i j yrow dmp1 Nyp1 yp1 yp1row
  clear rd1 rdp1 innp1 Rdnn1 Rdnnp1 Inn1p1
  clear crit1 crit2 nofn1(dim) nofn2(dim)
end;
%
figure(1);
title('Number of false nearest neighbors by criterion 1, R_{tol}=10');
%title('Number of false nearest neighbors by criterion 1, R_tol=25');
plot(nofn1,'k');
ylabel('false neighbors count');
xlabel('Fig.1 Change in false neighbors with R_{tol}=10');
%xlabel('Fig.1 Change in false neighbors with R_tol=25');
pause;
%
figure(2);
title('Number of false nearest neighbors by criterion 2, A_{tol}=2');
%title('Number of false nearest neighbors by criterion 2, A_tol=1');
plot(nofn2,'k');
ylabel('false neighbors count');
xlabel('Fig.2 Change in false neighbors with A_{tol}=2');
%xlabel('Fig.2 Change in false neighbors with A_tol=1');
% axis([0 10 0 50]);
    
    

