function y=exactfunction(x)                //    the function to optimize
n=length(x);
y=n+sum(x.^2-cos(2*%pi*x));            // Rastrigin                           
endfunction
//
//----------------------------------------------------------------
//
function y=approxfunction(X,Xexact,Xexactsel,u)  // the approximated function by RBF interpolation
[N1,N2]=size(Xexact);n=N2-1;
[Nsize,N3]=size(Xexactsel);
y=0;
for i=1:Nsize 
   y=y+u(i)*phi(X-Xexactsel(i,1:n))          
end
endfunction
//
//----------------------------------------------------------------
//----------------------------------------------------------------
//
function y=phi(x)             // RBF function
// y=exp(-norm(x)^2);          // gaussian
 c=1;y=sqrt(c^2+norm(x)^2); // multiquadric
endfunction
//
//----------------------------------------------------------------
//
function [A,b]=construct(Xexact,Xexactsel) // construction of the matrix for interpolation
  //
eps=1E-6;
[N1,N2]=size(Xexact);n=N2-1;
[Nsize,N3]=size(Xexactsel);
A=zeros(N1,Nsize);b=zeros(N1,1);
for i=1:N1
  for j=1:Nsize
    A(i,j)=phi(Xexact(i,1:n)-Xexactsel(j,1:n));
  end
  b(i)=Xexact(i,n+1);
end
//A=A+eps*eye(A); // for condition reduction
endfunction
//
//----------------------------------------------------------------
function evalpop=evaluation(X);    // exact evaluation of the whole population
[Npop,N2]=size(X);n=N2-1;
evalpop=zeros(Npop,1);
for i=1:Npop;
  evalpop(i)=exactfunction(X(i,1:n));
end
endfunction
//
/////////////////////
n=evstr(x_dialog('valeur de n (traces pour n=1 ou 2)','2'));
N=evstr(x_dialog('nombre de centres','100'));
Nsize=evstr(x_dialog('nombre de degres de liberte','100'));
Xlo=evstr(x_dialog('borne minimum domaine','-1'));
Xup=evstr(x_dialog('borne maximum domaine','1'));
//------------------------------------------------------------
//
Xexact=(Xup-Xlo)*rand(N,n+1)+Xlo*ones(N,n+1);
evalpop=evaluation(Xexact);
Xexact(:,n+1)=evalpop;
//
uu=1:Nsize;
Xexactsel=Xexact(uu,:);
//
[A,b]=construct(Xexact,Xexactsel);
//u=inv(A'*A)*A'*b;  // case of least square minimization Nsize<N
u=inv(A)*b;
//
disp('condition number:');disp(cond(A));
//------------------------------------------------------------------------
///////////////////////////////////////////////////////////////
if (n==2) then  // plot in the case n=2
//
Ntest=40;
x=Xlo:((Xup-Xlo)/(Ntest-1)):Xup;
y=x;
z=zeros(Ntest,Ntest);
zapprox=zeros(Ntest,Ntest);
for i=1:Ntest
for j=1:Ntest
  z(i,j)=exactfunction([x(i),y(j)]);
  zapprox(i,j)=approxfunction([x(i),y(j)],Xexact,Xexactsel,u);
end
end
//
xset('window',0);
clf()
xsetech([0,0,1.0,0.5])
contour2d(x,y,z,[0:0.01:0.1,0.2:1,1:10],rect=[Xlo,Xlo,Xup,Xup]); 
xtitle('rastrigin exacte');
xsetech([0,0.5,1.0,0.5])
contour2d(x,y,zapprox,[0:0.01:0.1,0.2:1,1:10],rect=[Xlo,Xlo,Xup,Xup]);
plot2d(Xexact(:,1),Xexact(:,2),-1,rect=[Xlo,Xlo,Xup,Xup])
plot2d(Xexact(1:Nsize,1),Xexact(1:Nsize,2),-2,rect=[Xlo,Xlo,Xup,Xup])
xtitle('rastrigin reconstruit avec RBF');
//
xset('window',1);
clf()
surf(x,y,z); 
xtitle('rastrigin exacte');
//
xset('window',2);
clf()
surf(x,y,zapprox);
xtitle('rastrigin reconstruit avec RBF');
end
 ///////////////////////////////////////////////   
 if (n==1) then  // plot in the case n=1
  Nplot=100;
  x=linspace(Xlo,Xup,Nplot);
  zex=zeros(Nplot,1);
  zkr=zeros(Nplot,1);  
  for i=1:Nplot
      zex(i)=exactfunction(x(i))         // exact
      zkr(i)=approxfunction(x(i),Xexact,Xexactsel,u); // model
  end
  xset('window',0)
  clf()
  plot2d(x,zex)
  plot2d(x,zkr,2)
  plot2d(Xexact(:,1),Xexact(:,2),-1)
  xtitle('exact function (black) vs RBF (blue)');
 end     
