function Xeval=evaluation(X);
  [Npop,n2]=size(X);n=n2-1;
  Xeval=X;
  for i=1:Npop;
    Xeval(i,$)=sum(X(i,1:n).^2-cos(2*%pi*X(i,1:n)))+n;
end
endfunction
///////////////////////////////
////////////////////////////////////////
function A=matricedistance(Xexact,h)
[N,n]=size(Xexact);
A=zeros(N,N);
for i=1:N
  for j=1:N
    A(i,j)=h(norm(Xexact(i,:)-Xexact(j,:)));
  end
end
endfunction  
//////////////////////////////////////////
function Xapprox=evaluationRBF(X,Xtot,h);
Xapprox=X;
[N,n2]=size(Xtot);
n=n2-1;
[Np,n3]=size(X);
A=matricedistance(Xtot(:,1:n),h);
w=linsolve(A,-Xtot(:,$));
for k=1:Np
  Jtilde=0;
  for i=1:N
    Jtilde=Jtilde+w(i)*h(norm(X(k,1:n)-Xtot(i,1:n)));
  end
  Xapprox(k,$)=Jtilde;
end  
endfunction
///////////////////////////////
function bestX=best(X);
[Npop,n2]=size(X);
[y,k]=sort(-X(:,$))
bestX=X(k(1),:);
endfunction
//////////////////////////////
function X=elitism(X,oldbestpar)
[Npop,n2]=size(X);
bestpar=best(X);
if(oldbestpar(:,$)<bestpar(:,$)) then
i=int(Npop*rand())+1;
X(i,:)=oldbestpar;
end
endfunction
//////////////////////////////// 
function Xsel=selection(X)  // selection par le rang
[Npop,n2]=size(X);
[y,k]=sort(X(:,$))  
X=X(k,:);  
p=1:Npop;
roulette=cumsum(p)/sum(p);
Xsel=[];
for i=1:Npop
  u=rand();
  w=find(u<roulette);
  n1=w(1);
  xsel=X(n1,:);
  Xsel=[Xsel;xsel];
end  
endfunction
////////////////////////////////// 
function xcrois=croisement(x)
// x: 2 lignes
xcrois=x;
u=rand();
xcrois(1,:)=u*x(1,:)+(1-u)*x(2,:);
xcrois(2,:)=(1-u)*x(1,:)+u*x(2,:);
endfunction
//////////////////////////////////
function xmut=mutation(x,gen,Ngen,xmin,xmax)
// x: 1 ligne
bet=5;
if (rand()<1/2) then
   xmut=x+rand()*(xmax-x)*(1-(gen-1)/Ngen)^bet;
else
   xmut=x-rand()*(x-xmin)*(1-(gen-1)/Ngen)^bet;
end
endfunction

//////////MAIN //////////////////////////
Npop=20;
Ngen=30;
n=2;
pc=0.9;
pm=0.6;
xmin=-1.5;xmax=1.5;
//
Ngeninit=2;
tau=2;
deff('y=h(x)','y=x');
//
Xtot=[];Nexacttot=[];
X=xmin+(xmax-xmin)*rand(Npop,n+1);
bestX=[];
for gen=1:Ngen
if (gen<Ngeninit) then 
  X=evaluation(X);   // evaluation
  Xtot=[Xtot;X];
else  
  X=evaluationRBF(X,Xtot,h);
  [y,k]=sort(-X(:,$)) ; 
  X=X(k,:);  
  Xexact=evaluation(X(1:tau,:));
  Xtot=[Xtot;Xexact];
end  
X=selection(X);    // selection
bestX=[bestX;best(Xtot)];
[n1,n2]=size(Xtot);
Nexacttot=[Nexacttot;n1];
//
oldbestpar=best(Xtot);
Xnew=[];
for j=1:Npop/2
   u1=int(Npop*rand())+1;
   u2=int(Npop*rand())+1;
   if (rand()<pc) then 
     Xnew=[Xnew;croisement(X([u1,u2],:))]; //croisement  
   else
     Xnew=[Xnew;X([u1,u2],:)];
   end  
end
for j=1:Npop
  if (rand()<pm) then
    Xnew(j,:)=mutation(Xnew(j,:),gen,Ngen,xmin,xmax);  // mutation
  end 
end  
X=Xnew;
end
//////////TRACE //////////////////////////
xbasc()
plot2d([1;Nexacttot],[Xtot(1,$);bestX(:,$)])
disp('meilleur individu:');disp(bestX($,1:n))
  
  
  
  
  
  
 
  
  
  
  
  
  
  
  




