function demo_LM_RBF()
syms c % adaptable parameters in an MLP network
M=40;d=2;N=800
% initialization
c0=rand(1,M*(d+1)+M+1)*2-1;
% preparation of training data
% uniform sampling
% Substitute to the target function g
z=rand(N,2)*2*pi-pi;
y=g_sin(z(:,1),z(:,2));%+rand(N,1)*0.02-0.01;
[xx means_z]=kmeans(z,M);
c0(1:M*d)=reshape(means_z',M*d,1);
c0(M*d+1:M*d+M) = abs(c0(M*d+1:M*d+M)); 
plot3(z(:,1),z(:,2),y,'.');
h = g_hat2(z(:,1),z(:,2),c0,M);
% Levenberg-Marquardt method
options = optimoptions('fsolve','Algorithm', 'levenberg-marquardt')
% specification of a nonlinear system for MLP_learning
f = @(c)learning_RBF(c,z(:,1),z(:,2),y,M);
% apply fsolve
% Verification of c_zero
c_zero = fsolve(f,c0,options);
mean((learning_RBF(c_zero,z(:,1),z(:,2),y,M)).^2)
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)
figure
range=pi;
x1=-range:0.1:range;
x2=x1;
for i=1:length(x1)
 C(i,:)=g_hat2(x1(i)*ones(length(x2),1),x2',c_zero,M);
end
mesh(x1,x2,C);
end

% a nonlinear system for MLP_learning
function F = learning_RBF(c,x1,x2,y,M)
A = [x1 x2];
d=2;W=[];
%W=reshape(c(1:M*(d+1)),d+1,M);
for j=1:M
 W=[W c((j-1)*(d)+1:j*(d))'];
end
Y=W';
X=A;
K=size(Y,1);N=size(X,1);
A=sum(X.^2,2)*ones(1,K);
C=ones(N,1)*sum(Y.^2,2)';
B=X*Y';
D=sqrt(A-2*B+C);
sigma=c(M*d+1:M*d+M);
Dsigma=ones(N,1)*sigma;
H=-D./(Dsigma.^2);
V=exp(H);
hh = V*c(M*(d+1)+1:end-1)';
F=hh+ones(length(x1),1)*c(end);
F =F-y;
end
function h = g_sin(x1,x2)
C1=[1 1/2]'; %weight
C2=[-1/2 1]';
A = [x1 x2];
h = sin(A*C1); %activation function
h=h+cos(A*C2);
end

function h = g_hat2(x1,x2,c,M)
A = [x1 x2];
d=2;W=[];
%W=reshape(c(1:M*(d+1)),d+1,M);
for j=1:M
 W=[W c((j-1)*(d)+1:j*(d))'];
end
Y=W';
X=A;
K=size(Y,1);N=size(X,1);
A=sum(X.^2,2)*ones(1,K);
C=ones(N,1)*sum(Y.^2,2)';
B=X*Y';
D=sqrt(A-2*B+C);
sigma=c(M*d+1:M*d+M);
Dsigma=ones(N,1)*sigma;
H=-D./(Dsigma.^2);
V=exp(H);
hh = V*c(M*(d+1)+1:end-1)';
h=hh+ones(length(x1),1)*c(end);
end
