%simualte cnn
close all;
clear all;
clc;

options = odeset('RelTol',1e-1, 'AbsTol', 1e-3, 'MaxStep', 1e-1);

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.Boundary='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.Boundary='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.Boundary='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.Boundary='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.4;
P.Boundary='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);

out=reshape(X,size(t,1),P.SizeX,P.SizeY);
%img=reshape(out(size(t,1),:,:),P.SizeX,P.SizeY);
img = squeeze(out(end, :, :));

figure;
ShowCNN(img);



% please compute here the helper variables to describe the whole system in
% a vector-ODE

N = P.SizeX;
M = P.SizeY;

z = P.Z;
Uhat = P.INPUT;
Ahat = zeros(N*M,N*M);
Bhat = zeros(N*M,N*M);
for i = 1:N*M
    for k = -1:1
        for l = -1:1
            Ahat(i, i+k*SizeX+l) = A(5+l+k*3);
            Bhat(i, i+k*SizeX+l) = B(5+l+k*3);
        end
    end
end

% please compute the Ljapunov function in every time-instant
V = zero(N,N);
for a=1:size(t,1)
   y = ...;
   part1 = ...;
   part2 = ...;
   part3 = ...;
   part4 = ...;
   V(a)=part1+part2+part3+part4;
end

figure
plot(V);



