%simualte cnn
clear all;
clc;
close all;
%just for sure

options = odeset('RelTol',1.0E-0,'AbsTol', 1.0E-0,'MaxStep', 1.0E-0);

Max_time = 10;


%{
%EDGE
A=[0 0 0; 0 1 0; 0 0 0];
B=[-1 -1 -1; -1 8 -1; -1 -1 -1];
P.Z=-1;
P.Bounadry='Periodic';
P.Input=ImgToCNN('Recall1.bmp');
P.SizeX=size(P.Input,1);
P.SizeY=size(P.Input,2);
InitState=zeros(P.SizeX,P.SizeY);
%}

%{
%Dilation
A=[0 0 0; 0 0 0; 0 0 0];
B=[0 1 0; 1 1 1; 0 1 0];
P.Z=4; 
P.Bounadry='Periodic';
P.Input=ImgToCNN('HOLLOW.bmp');
P.SizeX=size(P.Input,1);
P.SizeY=size(P.Input,2);
InitState=zeros(P.SizeX,P.SizeY);
%}


%{
%EROSION
A=[0 0 0; 0 0 0; 0 0 0];
B=[0 1 0; 1 1 1; 0 1 0];
P.Z=-4;  
P.Bounadry='Periodic';
P.Input=ImgToCNN('Recall1.bmp');
P.SizeX=size(P.Input,1);
P.SizeY=size(P.Input,2);
InitState=zeros(P.SizeX,P.SizeY);
%}

%{
%RECALL
A=[0.5 0.5 0.5; 0.5 4.0 0.5; 0.5 0.5 0.5];
B=[0 0 0; 0 0.4 0; 0 0 0];
P.Z=2.5; 
P.Bounadry='Constant';
P.BounValue=-1;
P.Input=ImgToCNN('Recall1.bmp');
P.SizeX=size(P.Input,1);
P.SizeY=size(P.Input,2);
InitState=ImgToCNN('Recall2.bmp');
%}


%THRES
A=[0 0 0; 0 2 0; 0 0 0];
B=[0 0 0; 0 0 0; 0 0 0];
P.Z=-0.8;
P.Bounadry='Periodic';
InitState=ImgToCNN('avergra2.bmp');
P.SizeX=size(InitState,1);
P.SizeY=size(InitState,2);
P.Input=zeros(P.SizeX,P.SizeY);


P.A=reshape(A,1,9);
P.B=reshape(B,1,9);

[t,X] = ode45(@cell_equation_boundary,[0 Max_time],InitState,options,P);

%show the image
out=reshape(X,size(t,1),P.SizeX,P.SizeY);
%cut off the border
img=reshape(out(size(t,1),:,:),P.SizeX,P.SizeY);
ShowCNN(img);


%Ahat, Bhat
z=ones(P.SizeX*P.SizeY,1)*P.Z;
Uhat=reshape(P.Input, P.SizeX*P.SizeY,1);
Ahat=zeros(P.SizeX*P.SizeY,P.SizeX*P.SizeY);
Bhat=zeros(P.SizeX*P.SizeY,P.SizeX*P.SizeY);

for a=1:(P.SizeX*P.SizeY)
   for k=-1:1
       for l=-1:1
           if ((a+l+k*P.SizeX)>0) && (a+l+k*P.SizeX<(P.SizeX*P.SizeY+1))
           Ahat(a,a+l+k*P.SizeX)= P.A(5+l+k*3);
           Bhat(a,a+l+k*P.SizeX)= P.B(5+l+k*3);
           end
       end
   end    
end



%ljapunov
V=zeros(size(t,1),1);
for a=1:size(t,1)
    y=nonlinearity1(X(a,:));
    part1=-0.5*y*Ahat*y';
    part2=inverseintegrate_nonlinearity(y);
    part3=-y*Bhat*Uhat;
    part4=-y*z;
   
    V(a)=part1+part2+part3+part4;
end
figure
plot(V);