%-----------------------------------------------------------------------
%
%   minimum_mutual_information_criterion_dim2.m
%                                  written by N. Chino on March 18, 2018
%                            revised a bit by N. Chino on March 29, 2018
%
%  This program determines an acceptable minimum mutual information sug-
%  gested by Fraser and Swinney (1986), by computing the specified mutual
%  information values corresponding to the discrete delay increments, xi's
%  on a certain dimensional embedding space.
%    In this program the descrete delay increment, xi, varies from 1 to 
%  100, and the embedding dimension is two.
%
%------------------------------------------------------------------------
%
%load data/chino/Henon-map5000.dat
%x=Henon_map5000;
load data/chino/Henon-map3000.dat
x=Henon_map3000;
%
%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;
%
disp('1. size of the original time series, x');
Nx=size(x,1)
%
disp('2. maximum number of embedding dimension');
dim=2
%
% Specify the maximum number of embedding dimension, dim, as the other
% embedding parameters.
%
disp('3. maximum number of discrete delay increment');
lag=50
%
for glag=1:lag
    %
% Compute the NT by dim coordinate matrix y from the original single time
% series, x, where NT=Nx-(dim-1)*lag.  Therefore, NT varies according to
% the two embedding parameter values, lag and dim.
%
Ny=Nx-(dim-1)*glag;
for i=1:Ny
  for j=1:dim
    y(i,j)=x(i+(j-1)*glag);
  end
end
disp('4. total number of reconstructed two-dimensional data');
yrow=size(y,1)
% 
% ymin and ymax are row vectors of the minima and maxima of the original
%   time series.
%
disp('5. minimum value of each of the two dimensions');
ymin=min(y)
disp('6. maximum value of each of the two dimensions');
ymax=max(y)
disp('7. range of the reconstruced data on each of the two dimensions');
rg=ymax-ymin
%
% Specify the number of classes of the coordinates in nine coordinate
% dimensions.
%Nc=[10,10,10,10,10,10,10]
disp('8. original number of classes in each of the two dimensions');
Nc=[4,4]
Ncscalar=4;
%
% Round rg./Nc to nearest 2 decimal digits.
% Here, the size of the decimal digits must be altered according to the 
% absolute values of the coordinates.
disp('9. width of each of the classes on each of the two dimensions');
len=round(rg./Nc,2)
%
disp('10. lower bound of the first class in each of the two dimensions');
limo=ymin-len./2
%
for j=1:dim
    limc(j,1)=limo(j);
    Ncmp2=Nc(j)+2;
    for k=2:Ncmp2
        limc(j,k)=limc(j,1)+len(j)*(k-1);
    end
end
%
disp('11. bounds of the final 5 classes in each of the two dimensions');
limc
%
Ncsp1=Ncscalar+1;
Ncsp2=Ncscalar+2;
% Define the two-dimensional array of zeros, fc, as iniial values
% of frequencies which fall in its specified cells.
fc=zeros(Ncsp1,Ncsp1);
%
% print out the total number of cells of the two-dimensional array
% for the frequency table.
disp('12. total number of cells of the frequency table');
nocells=Ncsp1^dim
%
% Identify the classes into which each point, i, falls on each of the
% two dimensions, and preserve them, as jc(i,j), i=1,Ny, j=1,dim.,
% for later use.
for i=1:Ny
    for j=1:dim
        Ncmp1=Nc(j)+1;
        for k=1:Ncsp1
            kp1=k+1;
            lb=limc(j,k);
            ub=limc(j,kp1);
            yy=y(i,j);
            if (yy >= lb) && (yy < ub)
                jc(i,j)=k;
            end
        end
    end
    j1=jc(i,1); j2=jc(i,2);
    fc(j1,j2)=fc(j1,j2)+1;
end
%
% test outputs of the frequency tables
%
%disp('frequency table for the first and second array');
%fc(:,:)
%
fcsum=0;
for j1=1:Ncsp1
    for j2=1:Ncsp1
        fcsum=fcsum+fc(j1,j2);
    end
end
%
disp('13. sum of the frequencies of the two-dimensional array, fc');
fcsum
%
% Compute the sum of fcsum along with each of the two dimensions
%
dim=1;
for k=1:Ncsp1
    fcd(dim,k)=0;
    for j2=1:Ncsp1
        fcd(dim,k)=fcd(dim,k)+fc(k,j2);
    end
end
%
dim=2;
for k=1:Ncsp1
    fcd(dim,k)=0;
    for j1=1:Ncsp1
        fcd(dim,k)=fcd(dim,k)+fc(j1,k);
    end
end
%
disp('14. sums of the freq. of each of the two classes on two dimensions'); 
fcd
%
pcd=fcd./Ny;
disp('15. estimated prob. into which each point falls each class on each dimension');
pcd
%
pjoint=fc./Ny;
%
for j1=1:Ncsp1
    for j2=1:Ncsp1
        if pjoint(j1,j2)>=1 j1,j2
        end
    end
end
%
M(glag)=0;
for i=1:Ny
    for j=1:dim
        js(j)=jc(i,j);
    end
    js1=js(1);js2=js(2);
    pjs=pjoint(js1,js2);
    arglog=pjs/(pcd(1,js1)*pcd(2,js2));
    Ms=pjs*log(arglog);
    M(glag)=M(glag)+Ms;
end
%
end
%
disp('16. Minimum mutual information values');
M
%
figure(1);
plot(M,'k');
ylabel('minimum mutual information');
xlabel('delay increment on the two-dimensional space');
pause;
%
Mmax=M(1);
for glag=1:lag
    Mp(glag)=100*M(glag)/Mmax;
end
%
disp('17. The minimum mutual information values');
M
%
figure(2);
plot(Mp,'k');
ylabel('per cent of the minimum mutual information');
xlabel('delay increment per cent on the two-dimensional space');
