////////////////////////////////
//Kriging model
/////////////////////////////////////
//
function y=f(x)   //    the function to model
n=length(x);
y=n+sum(x.^2-cos(2*%pi*x));                 
endfunction
/////////////////////////////////////////////////////////////
function x= annealing (f, x0,Nmax)   // optimization with simulated annealing
  Xmin=0.1*ones(x0);Xmax=5*ones(x0);
y0 = f(x0);
x = x0; y = y0; alpha = 0.9; T = 1; s = 0.1;icount=1;
Neval=[];
Fbest=[];
for k = 1 : 20
   T = alpha * T; 
   accept = 0;
   for l = 1 : Nmax/20
      x1 = max(Xmin,min(Xmax,x0 + s*(0.5*ones(x0) - rand (x0))));
      y1 = f(x1);
      icount=icount+1;
      dy = y1 - y0;
      if rand() < exp (- dy / T) then
         x0 = x1; y0 = y1; accept = accept + 1;
         if y0 < y then x = x0; y = y0; end
      end
   end
   if accept < 25/100*Nmax/20 then s = s / 2; end
   if accept > 75/100*Nmax/20 then s = 2 * s; end
   x0 = x; y0 = y;
   Neval=[Neval;icount];
   Fbest=[Fbest;y];
end
endfunction
//////////////////////////////////////////////
function y=c(X,Y,T)        // covariance 
  y=T(1)*exp(-sum((X-Y).^2./T(3:$)))+T(2)
endfunction 
//////////////////////////////////////////////
function A=cov(Xexact,T) // covariance matrix
A=zeros(N,N);
for i=1:N
  for j=1:N
    A(i,j)=c(Xexact(i,1:n),Xexact(j,1:n),T);
  end
end 
endfunction    
///////////////////////////////////////////////
function y=vrai(T)     // maximum likelihood function
  A=cov(Xexact,T);
  if cond(A)>1E10 then
     y=1E10
  else
      y=log(det(A))+Xexact(:,n+1)'*inv(A)*Xexact(:,n+1)
  end    
endfunction  
//////////////////////////////////////////////
function y=fkrige(X)     // kriging function
  b=zeros(N,1);
  for i=1:N
    b(i)=c(X,Xexact(i,1:n),T);
  end  
  omega=linsolve(A,-b);
  y=omega'*Xexact(:,n+1);
endfunction   
/////////////////////
n=evstr(x_dialog('valeur de n (traces pour n=1 ou 2)','2'));
N=evstr(x_dialog('valeur de N','20'));
Xlo=evstr(x_dialog('borne minimum domaine','-1'));
Xup=evstr(x_dialog('borne maximum domaine','1'));
////////////////////
//
Xexact=(Xup-Xlo)*rand(N,n+1)+Xlo; // position of points
for i=1:N
  Xexact(i,n+1)=f(Xexact(i,1:n));
end 
T0=[1,1,ones(1,n)]; 
T=annealing(vrai,T0,2000);       // optimal parameters
A=cov(Xexact,T)
disp('optimal parameters:');disp(T);
disp('condition number:');disp(cond(A));
//////////////////////////////////////////////////
if (n==2) then  // plot in the case n=2
  Nplot=100;
  x=linspace(Xlo,Xup,Nplot);
  y=linspace(Xlo,Xup,Nplot);
  zex=zeros(Nplot,Nplot);
  zkr=zeros(Nplot,Nplot);  
  for i=1:Nplot
    for j=1:Nplot
      zex(i,j)=f([x(i),y(j)]);        // exact
      zkr(i,j)=fkrige([x(i),y(j)]);  // model
    end
  end
  xset('window',0)
  clf()
  xsetech([0,0,1.0,0.5])
  contour2d(x,y,zex,[0:0.01:0.1,0.2:1,1:10])
  xtitle('exact function');
  xsetech([0,0.5,1.0,0.5])
  contour2d(x,y,zkr,[0:0.01:0.1,0.2:1,1:10])
  plot2d(Xexact(:,1),Xexact(:,2),-1)
  xtitle('kriging function');
  xset('window',1)
  clf()
  surf(x,y,zex)
  xtitle('exact function');
  xset('window',2)
  clf()
  surf(x,y,zkr)
  xtitle('kriging function');
 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)=f(x(i))         // exact
      zkr(i)=fkrige(x(i))  // model
  end
  xset('window',0)
  clf()
  plot2d(x,zex)
  plot2d(x,zkr,2)
  plot2d(Xexact(:,1),Xexact(:,2),-1)
  xtitle('exact function (black) vs kriging (blue)');
 end      
  
    
 

 
