//
//------------------------------------------------------
// hybrid method (ES + steepest descent) for Rastrigin
disp('hybrid method (ES + SD) for Rastrigin');
// 
//
//-----------------------------------------------------
clear();
//
function y=rastrigin(x)  // the function to optimize
n=prod(size(x));
y=n+sum(x.^2-cos(2*%pi*x));
endfunction
//
//-----------------------------------------------------
function y=rastrigingrad(x)                     // the gradient of the function to optimize
y=2*x+2*%pi*sin(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);
[mu,n]=size(Xpar);
[lambda,n]=size(Xoff);
n=n-1;
Xtot=[Xpar;Xoff];
[Fsort,ind]=sort(-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
//
//-----------------------------------------------------
function z=dotprod(x,y);                        // computes the dot product of x and y
z=sum(x.*y);
endfunction
//-----------------------------------------------------
function d=descentdirection(f,x,fx,gx);         // descent direction: gradient
d=-gx;
endfunction
//----------------------------------------------------
function [xnew,fnew,itback]=backtracking(f,x,fx,gx,d,tau,bet,alphainit);// line search by backtracking until Armijo condition
alpha=alphainit;xnew=x+alpha*d;
fnew=f(xnew);
itback=1;
while(fnew>fx+bet*alpha*dotprod(gx,d))
  alpha=tau*alpha;
  xnew=x+alpha*d;
  fnew=f(xnew);
  itback=itback+1;
end
endfunction
//-------------------------------------------------
//                                           // steepest descent method with Armijo condition
function [xbest,fbest,itfct,itgrad]=steepest(x,fx,epsilongrad,itermax,deltagain,tau,bet,alphainit)
itgrad=0;itfct=0;
gx=rastrigingrad(x);itgrad=itgrad+1;
//
delta=1;
while (norm(gx)>epsilongrad)&(delta>0)&(itgrad<itermax)
  d=descentdirection(rastrigin,x,fx,gx);  
  [x,fxnew,itback]=backtracking(rastrigin,x,fx,gx,d,tau,bet,alphainit);
  delta=(fx-fxnew)/(itback+prod(size(x)))-deltagain;
  fx=fxnew;
  gx=rastrigingrad(x);
  itgrad=itgrad+1;
  itfct=itfct+itback;
end
xbest=x;fbest=fx;
endfunction
//
//-------------------------------------------------// clustering
//
function [Xclust,ind,Matr,Numb]=clustering(Xpar,Nclust);
[Npop,n2]=size(Xpar);n=n2-1;
bari=zeros(Nclust,n);
Matr=zeros(Nclust,Nclust);
Numb=zeros(1,Nclust);
Xclust=zeros(Nclust,n+1);
Xparclust=[Xpar,1+int(Nclust*rand(Npop,1))];

inon=1;
while inon==1;
 inon=0;
 for i=1:Nclust
  index=find(Xparclust(:,$)==i);
  if (length(index)>0) then
  bari(i,:)=1/length(index)*sum(Xparclust(index,1:n),'r');
  end
 end
 for i=1:Npop
   dist0= norm(Xparclust(i,1:n)-bari(Xparclust(i,$),:));
   for j=1:Nclust
      dist(j)=norm(Xparclust(i,1:n)-bari(j,:));
      if dist(j)<dist0 then 
          inon=1;
          Xparclust(i,$)=j;dist0=dist(j);
       end
    end   
 end     
end
 for i=1:Nclust
  index=find(Xparclust(:,$)==i);
  [minval,indic]=mini(Xparclust(index,n+1)); 
  ind(i)=index(indic(1));
  if (length(index)>0) then
    Xclust(i,:)=Xpar(index(indic(1)),:);
  end
  Numb(i)=length(index);
 end 
for i=1:Nclust
 for j=1:Nclust
    Matr(i,j)=norm(bari(i,:)-bari(j,:));
 end
end
endfunction
//
//
//-------------------------------------------dynamic results displays

//
function result(Ntot,fmin,co)
nl=length(Ntot);
if (nl>1) then
xset('window',0);
plot2d([Ntot($-1),Ntot($)],[fmin($-1),fmin($)],co);
xtitle('convergence history','Neval','fmin');
end
endfunction
//
//-----------------------------------------------------
//-----------------------------------------------------
//---------main program -------------------------------
//
xset('window',0);xbasc();
// ----------------------------------------------------
//
n=6;                          // number of variables of the optimization problem
mu=5;                           //parent number       
lambda=35;                       // offspring number 
siginit=0.5;                       // initial mutation strength
alpha=1.03;                       // 1/5 rule coefficient
Ngen=60;
//
Nclustmax=0;                   // maximal cluster number
itgradmax=10;                  // maximal steepest descent iteration
valmin=3;                      // minimal number of consecutive GA iterations
CVRlim=1;                      // coefficient of variation of the population threshold
//
//
epsilongrad=1E-10;               // stopping criterion for gradient value
tau=0.3;                        // backtraking coefficient
bet=1E-4;                       // Armijo coefficient
alphainit=0.1;                  // first step length of steepest descent 
//
// ----------------------------------------------------
//
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;Ntot=[];
mineval=0;newval=zeros(mu,1);itfct=0;fmin=[];
CVold=1;valgen=1;Nevalg=0;
//
//-----------------------------------------------------
//-----------------------------------------------------
//
sigma=siginit;
Xpar=evaluation(Xpar);
Nevalf=mu;
Nclust=Nclustmax;
gen=0;
//
//     MAIN LOOP : BEGINNING
//
//
for gen=1:Ngen
  for i=1:lambda
    u=int(mu*rand())+1;
    xpar=Xpar(u,:);
    xoff=mutate(xpar,sigma);              // ES mutation
    Xoff(i,:)=xoff;
  end
//
oldval=Xoff(:,$);
Xoff=evaluation(Xoff);                    // evaluation  of the offsprings
Nevalf=Nevalf+lambda;
//
if (gen==1) then Nold=Nevalf;oldmean=mean(Xoff(:,$));end
//
Xpar=selection(Xpar,Xoff);                // selection lambda+mu
bestpar=Xpar(1,1:n);
fmin=[fmin,Xpar(1,$)];
Ntot=[Ntot,Nevalf];
result(Ntot,fmin,1);                      // dynamic result display
//
newval=Xoff(:,$);
ind=find(oldval>newval);
taux=length(ind)/lambda;
if (taux<0.2) then sigma=sigma/alpha; else sigma=alpha*sigma;end
//
// ---------------------------------------------------------------------------
//
CVnew=mean(Xoff(:,$))/st_deviation(Xoff(:,$));
CVR=CVnew/CVold;
CVold=CVnew;
//  
if (CVR>CVRlim)&(Nclust>0)&(gen-valgen>valmin) then           // shift to steepest descent method
   newmean=mean(Xoff(:,$));Nnew=Nevalf;
   deltagain=(oldmean-newmean)/(Nnew-Nold);
   oldmean=newmean;Nold=Nnew;
//  
   [Xclust,ind,Matr,Numb]=clustering(Xoff,Nclust);
   Xoffnew=[];
   for i=1:Nclust
      if (Numb(i)>0) then
        [xnew,fnew,itfct,itgrad]=steepest(Xclust(i,1:n),Xclust(i,$),epsilongrad,itgradmax,deltagain,tau,bet,alphainit);
        Nevalf=Nevalf+itfct+itgrad*n;
        Nevalg=Nevalg+itgrad;
        Xoffnew=[Xoffnew;xnew,fnew];
      end
   end  
   Xpar=selection(Xpar,Xoffnew);                // selection lambda+mu
   bestpar=Xpar(1,1:n);
   fmin=[fmin,Xpar(1,$)];
   Ntot=[Ntot,Nevalf];
   result(Ntot,fmin,2);                      // dynamic result display
   valgen=gen;
end
//------------------------------------------------------------------------------
//
//     MAIN LOOP : ENDING
//
end
//-----------------------------------------------------
//
//
//
disp('minimum obtained:');disp(bestpar(1:n));
disp('corresponding value by f:');disp(fmin($));
disp('evaluation number:');disp(Nevalf);
//
