//
// stratégie d'évolution
// variante: adaptation de sigma: loi du 1/5
// variante: stratégie plus ou comma
// choix de fonction: Rastrigin ou Rosenbrock
//
////////////////////////////////////////////////
n=2;   // n=2 si fonction rosen
lambda=50;
mu=5;
Ngen=200;
sigma=1;
/////////////////////////////////////////////////
function y=rast(x)
    y=sum(x.*x-cos(2*%pi*x))+length(x);
endfunction
///////////////////////////////////////////////////
/////////////////////////////////////////////////
function y=rosen(x)
    y=100*(x(2)-x(1)^2)^2+sum(x(1)-1).^2;
endfunction
///////////////////////////////////////////////////

Ap=10*rand(mu,n+1)-5*ones(mu,n+1)// initialisation
for i=1:mu
   Ap(i,n+1)=rosen(Ap(i,1:n))   // evaluation (fonction rosen ou rast) 
end
Ae=zeros(lambda,n+1)
BestJ=[]
BestX=[]
//
for i=1:Ngen
// 
  tau=0;   
  X=1/mu*sum(Ap(:,1:n),'r')  // croisement
  val=rosen(X);    
  for i=1:lambda
     Ae(i,1:n)=X+sigma*grand(1,n,'nor',0,1) // mutation
     Ae(i,n+1)=rosen(Ae(i,1:n))   // evaluation (fonction rosen ou rast) 
     if (Ae(i,n+1)<val) then
        tau=tau+1/lambda
     end    
  end
  Atot=[Ap;Ae]; // selection plus (enfants+ parents)
//  Atot=Ae        // selection virgule (enfants)
  [a,b]=gsort(-Atot(:,n+1))
  Atot=Atot(b,:)
  Ap=Atot(1:mu,:)  
  BestJ=[BestJ,Atot(1,n+1)]
  BestX=[BestX;Atot(1,1:n)]
  sigma=exp(1/3*(tau-1/5)/(1-1/5))*sigma
end 
figure(1)
clf() 
plot(1:Ngen,BestJ)
figure(2)
clf()
for j=1:n
 plot2d(1:Ngen,BestX(:,j),j)
end


  

