%% Solving a nonlinear system for LM learning, created by Jiann-Ming Wu
% Department of Applied Mathematics, National Dong Hwa University %
% Using Matlab (R)
% 2018 dec. 25
% learning is resolved by the Levenberg-Marquardt method %
function demo_LM_CNN_regressions()
syms c % adaptable parameters in an MLP network 
M=3;d=28;N=300
f_size=5; 
% initialization
c0_size=f_size*f_size*M; 
c0_size=c0_size+(d-f_size+1)^2 
c0=rand(1,c0_size)*2-1;
size(c0)
% preparation of training data
% uniform sampling
% Substitute to the target function g
z=rand(d,d,N);
A=z>0.5;
y=rand(N,1)*0.3-0.15;
size(y)
y_hat = g_hat_CNN(z,c0,M,f_size);
size(y_hat)

[trainImages,~,trainAngles] = digitTrain4DArrayData;
for n=1:N
    A(:,:,n)=squeeze(trainImages(:,:,1,n));
end
y=trainAngles(1:N);
% Levenberg-Marquardt method
options = optimoptions('fsolve','Algorithm', 'levenberg-marquardt') 
% specification of a nonlinear system for MLP_learning
f = @(c)learning_CNN(c,A,y,M,f_size);
% apply fsolve
% Verification of c_zero
c_zero = fsolve(f,c0,options);
BB=c_zero;
F = learning_CNN(BB,A,y,M,f_size);
mean((F).^2)
end
%mean((g_hat2(z(:,1),z(:,2),c_zero,M)-y).^2)
% testing phase
%z_test=rand(100,2)*2*pi-pi; %y_test=g_sin(z_test(:,1),z_test(:,2)); %mean((g_hat2(z_test(:,1),z_test(:,2),c_zero,M)-y_test).^2) end
function F = learning_CNN(BB,A,y,M,f_size) 
    B_size=M*f_size*f_size;
B=BB(1:B_size);
WW=BB(B_size+1:end);
%
for i=1:M
cc(i,:) = B((i-1)*f_size^2+1:i*f_size^2); 
for j=1:f_size
C(j,:,i)=cc(i,(j-1)*f_size+1:j*f_size); 
end
end
W_size=sqrt(length(WW)) 
for k=1:W_size
W(k,:)=WW((k-1)*W_size+1:k*W_size); 
end
N=size(A,3); 
for n=1:N
v(:,:,n)=zeros(size(A,2)-f_size+1,size(A,2)-f_size+1); 
for i=1:M
h(:,:,i)=convn(A(:,:,n),C(:,:,i),'valid');
v(:,:,n)=v(:,:,n)+tanh(h(:,:,i)); 
end
y_hat(n)=sum(sum(v(:,:,n).*W)); 
end
F=y_hat;
F =F-y;
end
function y_hat = g_hat_CNN(A,BB,M,f_size)
B_size=M*f_size*f_size; 
B=BB(1:B_size); 
WW=BB(B_size+1:end); %
for i=1:M
cc(i,:) = B((i-1)*f_size^2+1:i*f_size^2); 
for j=1:f_size
C(j,:,i)=cc(i,(j-1)*f_size+1:j*f_size); 
end
end
W_size=sqrt(length(WW)) 
for k=1:W_size
W(k,:)=WW((k-1)*W_size+1:k*W_size); 
end
N=size(A,3); 
for n=1:N
v(:,:,n)=zeros(size(A,2)-f_size+1,size(A,2)-f_size+1); 
for i=1:M
h(:,:,i)=convn(A(:,:,n),C(:,:,i),'valid');
v(:,:,n)=v(:,:,n)+tanh(h(:,:,i)); 
end
y_hat(n)=sum(sum(v(:,:,n).*W)); 
end
end
