//
//------------------------------------------------------
//
// Evolution Strategy
// for the minimization of the rastrigin function 
//
//-----------------------------------------------------
//
function y=rastrigin(x)       // the function to optimize
n=prod(size(x));
y=n+sum(x.^2-cos(2*%pi*x));
endfunction
//
//-----------------------------------------------------
function Xpar=evaluation(X);  // evaluation of the whole population
Xpar=X;
[Npo,n]=size(X);n=n-1;
evalpop=zeros(Npo,1);
for i=1:Npo;
 evalpop(i)=rastrigin(X(i,1:n));
end
Xpar(:,n+1)=evalpop;
endfunction
//
//-----------------------------------------------------
function Xparnew=selection(Xpar,Xoff); // selection plus or comma
[mu,n]=size(Xpar);
[lambda,n]=size(Xoff);
n=n-1;
Xtot=[Xpar;Xoff];    // plus selection
//Xtot=[Xoff];      // comma selection
[Fsort,ind]=gsort(-Xtot(:,n+1));
Xtot=Xtot(ind,:);
Xparnew=Xtot(1:mu,:);
endfunction
//
//-----------------------------------------------------//
function xoff=mutate(xpar,sigma) // normal mutation
[np,n]=size(xpar);n=n-1;
xoff=xpar;
u=rand(1,n,'n');
xoff(1,1:n)=xpar(1,1:n)+sigma*u;
endfunction
//
//-----------------------------------------------------
//-----------------------------------------------------
//---------main program -------------------------------
//
disp('real valued ES for minimizing black box function');
//
// ----------------------------------------------------
//
n=evstr(x_dialog('parameter number of the rastrigin function','2')); 
mu=evstr(x_dialog('parent number','5')); 
lambda=evstr(x_dialog('offspring number','50'));    
Ngen=evstr(x_dialog('generation number','20')); 
// 
siginit=1;                       // initial mutation strength
alpha=1.25;                      // 1/5 rule coefficient
//
// ----------------------------------------------------
//
xmin=-5.*ones(1,n);
xmax=5.*ones(1,n);
Xmin=ones(mu,1)*xmin;
Xmax=ones(mu,1)*xmax;
u=rand(mu,n);
pop=Xmin+(Xmax-Xmin).*u;  // random initialisation of the population
Xpar=[pop,rand(mu,1)];  Xoldpar=Xpar;
Xoff=Xpar;min=[];Ntot=[];
mineval=0;newval=zeros(mu,1);itfct=0;fmin=[];
//
//-----------------------------------------------------
//-----------------------------------------------------
//
sigma=siginit;
Xpar=evaluation(Xpar);
Nevalf=mu;
//
//     MAIN LOOP : BEGINNING
//
//
for gen=1:(Ngen-1)
  //
  xcros=1/mu*sum(Xpar,'r');         // parent crossover
 // 
  for i=1:lambda
    u=int(mu*rand())+1;
    xpar=Xpar(u,:);
      xoff=mutate(xcros,sigma);        // parent crossover+mutation 
//    xoff=mutate(xpar,sigma);         // parent mutation, no crossover
    Xoff(i,:)=xoff;
  end
//
oldval=Xoff(:,$);
Xoff=evaluation(Xoff);                  // offspring evaluation 
Nevalf=Nevalf+lambda;
//
newval=Xoff(:,$);
ind=find(oldval>newval);
taux=length(ind)/lambda;
//
// ---  adaptativity of sigma-----------------
//
if (taux<0.2) then sigma=sigma/alpha; else sigma=alpha*sigma;end
//
//-----------------------------------------------------
//
Xpar=selection(Xpar,Xoff);              // selection plus or comma
fmin=[fmin,Xpar(1,$)];
bestpar=Xpar(1,1:n);
Ntot=[Ntot,Nevalf];
//
//     MAIN LOOP : ENDING
//
end
//-----------------------------------------------------
//-----------------------------------------------------
//
//
xset('window',0);
clf();
plot2d(Ntot,fmin)
xtitle('convergence history','Neval','fmin');
//
disp('minimum obtained:');disp(bestpar(1:n));
disp('corresponding value by f:');disp(fmin($));
disp('evaluation number:');disp(Nevalf);
//
