//
//------------------------------------------------------
// 
disp('ES with CSA adaptation for can problem');
//
//-----------------------------------------------------
function y=can(x)       // the function to optimize: surface of the can
y=%pi*x(1)^2/4+%pi*x(1)*x(2);
endfunction
//
function y=volume(x);   // the volume of the can
  y=%pi*x(1)^2*x(2)/4;
endfunction  
//
function y=cantotdeath(x)    //  the death penalty
if (x(1)<1)|(volume(x)<300)  then 
  y=1E10;
else
  y=can(x)
end  
endfunction  
function y=cantotpenal(x)    //  dynamic penalization
if (x(1)<1)|(volume(x)<300)  then 
  y=1E10;
else
  y=can(x)+K*(max(0,(1-x(1))^2)+max(0,(300-volume(x))^2));
end  
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;
///-------------------------------------------------------------------------------------------- 
///    CHANGE HERE THE FUNCTION : cantotdeath or cantotpenal
///-------------------------------------------------------------------------------------------- 
  evalpop(i)=cantotdeath(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]=gsort(-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=10;                           //parent number       
lambda=100;                       // offspring number 
siginit=10;                       // initial mutation strength
Ngen=100;                        // number of generations
K=100;                            // initial penalization
tau=2;                        // geometric rate
// ------- parameters for CSA adaptation-------------
c=1/sqrt(n);
chiN=sqrt(n)*(1-1/(4*n)+1/(21*n^2));
D=sqrt(n);
//
// ----------------------------------------------------
//
xmin=1.*ones(1,n);
xmax=10.*ones(1,n);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;
//----------------------------------------------------
xset('window',1);
XV=linspace(1,10,100);YV=1200/%pi*ones(XV)./(XV.^2);
xset('pixmap',1);
//
// ---- MAIN LOOP : BEGINNING---------------------------
//------------------------------------------------------
//
for gen=1:(Ngen-1)
//
if (modulo(gen,10)==0) then 
  K=K*tau;  // Dynamic penalization
end  
//
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
//
clf();
plot2d(XV,YV,2,rect=[xmin,xmax]);
plot2d(Xpar(:,1),Xpar(:,2),-1,rect=[xmin,xmax]);
xset('wshow');
//xpause(100000)
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
//-----------------------------------------------------
//-----------------------------------------------------
//
//
dbest=(2400/%pi)^(1/3);hbest=1200/(%pi*dbest^2);
plot2d(dbest,hbest,-2);
//
xset('pixmap',0);
xset('window',0);
plot2d(Ntot,fmin,rect=[0,190,Nevalf,250]);
disp('best can obtained:');disp(bestpar(1:n));
disp('volume obtained:');disp(volume(bestpar(1:n)));
disp('corresponding surface:');disp(can(bestpar(1:n)));
disp('best possible can'); disp([dbest,hbest]);
disp('best possible surface:');disp(can([dbest,hbest]));
disp('number of evaluations:');disp(Nevalf);
//
