//
//------------------------------------------------------
// 
disp('ES with CSA adaptation for Rastrigin');
//
//-----------------------------------------------------
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,indbest]=selection(Xpar,Xoff);
[mu,n]=size(Xpar);
[lambda,n]=size(Xoff);
n=n-1;
//Xtot=[Xpar;Xoff];          // selection plus
Xtot=Xoff;                   // selection comma
[Fsort,ind]=sort(-Xtot(:,n+1));
indbest=ind(1:mu);
Xparnew=Xtot(indbest,:);
endfunction
//
//-----------------------------------------------------
//-----------------------------------------------------
//---------main program -------------------------------
//
xset('window',0);xbasc();
disp('real valued ES for minimizing black box function');
//
// ----------------------------------------------------
//
n=2;                          // number of variables of the optimization problem
mu=5;                           //parent number       
lambda=50;                       // offspring number 
siginit=1;                       // initial mutation strength
Ngen=20;
// ------- parameters for CSA adaptation-------------
c=1/sqrt(n);
chiN=sqrt(n)*(1-1/(4*n)+1/(21*n^2));
D=sqrt(n);
//
// ----------------------------------------------------
//
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;
V=zeros(1,n);
Xpar=evaluation(Xpar);
Y=1/mu*sum(Xpar(:,1:n),1);
Nevalf=mu;
//
// ---- MAIN LOOP : BEGINNING---------------------------
//------------------------------------------------------
//
for gen=1:(Ngen-1)
//
W=[];
  for i=1:lambda
    w=sigma*rand(1,n,'normal');  // mutation       
    W=[W;w];
    xoff=Y+w;   
    Xoff(i,1:n)=xoff;
end
//
Xoff=evaluation(Xoff);          // evaluation 
Nevalf=Nevalf+lambda;
//
[Xpar,ind]=selection(Xpar,Xoff); // selection
//
fmin=[fmin,Xpar(1,$)];
bestpar=Xpar(1,1:n);
Ntot=[Ntot,Nevalf];
//
// ---  adaptativity of sigma-----------------
//
Z=1/mu*sum(W(ind,:),1);
Y=Y+Z;
V=(1-c)*V+sqrt(c*(2-c))*sqrt(mu)/sigma*Z;
sigma=sigma*exp((norm(V)-chiN)/(D*chiN));
//
//------------------------------------------------------
//-------MAIN LOOP : ENDING-----------------------------
//
end
//-----------------------------------------------------
//-----------------------------------------------------
//
//
//
plot2d(Ntot,fmin);
disp('minimum obtained:');disp(bestpar(1:n));
disp('corresponding value by f:');disp(fmin($));
disp('evaluation number:');disp(Nevalf);
//
