//
//------------------------------------------------------
// hybrid method (GA + steepest descent) for Rastrigin
disp('hybrid method (GA + 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));
y=sum(x.^2-10*cos(2*%pi*x)+10);
//y=n+exp(1)-n*exp(-0.2*sqrt(1/n*sum(x.^2)))-exp(1/n*sum(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);
y=2*x+10*2*%pi*sin(2*%pi*x);
//n=prod(size(x));
//y1=-n*exp(-0.2*sqrt(1/n*sum(x.^2)))*((-0.2)*1/n*2*x)/(2*(sqrt(1/n*sum(x.^2))));
//y2=-exp(1/n*sum(cos(2*%pi*x)))*1/n*(-2*%pi*sin(2*%pi*x));
//y=y1+y2;
endfunction
//
//-----------------------------------------------------
//
//-----------------------------------------------------
function evalpop=evaluation(X);  // evaluation of the whole population
[Npop,n]=size(X);n=n-1;
evalpop=zeros(Npop,1);
for i=1:Npop;
  evalpop(i)=rastrigin(X(i,1:n));
end
endfunction
//
//-----------------------------------------------------
function [bestpar,bestval]=thebest(Xtot);
[mu,n]=size(Xtot);n=n-1;
[Fsort,indice]=sort(-Xtot(:,n+1));
bestpar=Xtot(indice(1),1:n);
bestval=Xtot(indice(1),n+1);
endfunction
//
//-----------------------------------------------------
//
//
function [Xpar,bestpar]=evalplus1elitism(Xpar,Xoldpar,Npop,gen) // evaluation plus 1-elitism
val=evaluation(Xpar);    
Xpar(:,n+1)=val;
[minval,index]=min(Xoldpar(:,n+1));  
bestpar=Xoldpar(index(1),:);
[minval,index]=min(val);
bestnewpop=Xpar(index(1),:);
if (bestnewpop($)>bestpar($))&(gen>1) then
    nrand=int(Npop*rand())+1;
    Xpar(nrand,:)=bestpar;
else 
    bestpar=bestnewpop;
end
endfunction
//
//-----------------------------------------------------
//
function Xsel=rankselect(X)    // ranking selection
Xsel=X;
[Npop,n]=size(X);n=n-1;
[Fsort,ind]=sort(X(:,n+1));
X=X(ind,:);
p=1:Npop;roulette=cumsum(p)/sum(p);
for i=1:Npop;
  u=rand();
  index=find(roulette<u);Nselect=max(index)+1;
  Xsel(i,:)=X(Nselect,:);
end
endfunction
//
//-----------------------------------------------------
//
function xoff=blcross(xpar)    // blend crossover
[np,n]=size(xpar);n=n-1;
xoff=xpar;
for i=1:n
   u=rand();
   xoff(1,i)=u*xpar(1,i)+(1-u)*xpar(2,i);
   xoff(2,i)=u*xpar(2,i)+(1-u)*xpar(1,i);
end
endfunction
//
//-----------------------------------------------------
//
function xoff=mutate(xpar,ngen,Ngen,xmin,xmax) // non uniform mutation
[np,n]=size(xpar);n=n-1;
b=5;
xoff=xpar;
for i=1:n
  u1=rand();u2=rand();
  if (u1<1/2) then
    xoff(1,i)=xpar(1,i)+(xmax(i)-xpar(i))*u2*(1-(ngen-1)/Ngen)^b;
   else
    xoff(1,i)=xpar(1,i)-(xpar(i)-xmin(i))*u2*(1-(ngen-1)/Ngen)^b;
   end;
end;
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);
// ----------------------------------------------------
//
n=6;                          // number of variables of the optimization problem
Npop=40;                           //parent number       
//Ngen=80;
pc=0.9;
pm=0.6;
Nevalfmax=4000;
Nclustmax=4;                   // 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(Npop,1)*xmin;
Xmax=ones(Npop,1)*xmax;
u=rand(Npop,n);
pop=Xmin+(Xmax-Xmin).*u;  // random initialisation of the population
Xpar=[pop,rand(Npop,1)];  Xoldpar=Xpar;
Xoff=Xpar;Ntot=[];
mineval=0;newval=zeros(Npop,1);itfct=0;fmin=[];
CVold=1;valgen=1;Nevalg=0;Nevalf=0;gen=0;
//
Nclust=Nclustmax;
//----------------------------------------------------------
//     MAIN LOOP : BEGINNING
//
//
while(Nevalf<Nevalfmax)
gen=gen+1;
//for gen=1:(Ngen-1)
//
  Xpar=Xoff;
  [Xpar,bestpar]=evalplus1elitism(Xpar,Xoldpar,Npop,gen);  // evaluation plus 1-elitism
  Nevalf=Nevalf+Npop;
  fmin=[fmin,bestpar($)];
  Ntot=[Ntot,Nevalf];
  result(Ntot,fmin,1);                      // dynamic result display
  if (gen==1) then Nold=Nevalf;oldmean=mean(Xoff(:,$));end
 // 
 //-------------------------------------------------------------------------
  CVnew=mean(Xpar(:,$))/st_deviation(Xpar(:,$));
  CVR=CVnew/CVold;
  CVold=CVnew;
//
//  ---------------------// shift to steepest descent method   -------------------------------
//
if (CVR>CVRlim)&(Nclust>0)&(gen-valgen>valmin) then           
   newmean=mean(Xpar(:,$));Nnew=Nevalf;
   deltagain=(oldmean-newmean)/(Nnew-Nold);
   oldmean=newmean;Nold=Nnew;
//  
   [Xclust,ind,Matr,Numb]=clustering(Xpar,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;
        Xpar(ind(i),:)=[xnew,fnew];
        Xoffnew=[Xoffnew;xnew,fnew];
      end
   end  
   [bestpar,bestval]=thebest(Xoffnew);                // selection lambda+mu
   fmin=[fmin,bestval];
   Ntot=[Ntot,Nevalf];
   result(Ntot,fmin,2);                      // dynamic result display
   valgen=gen;
end
 //  
 //--------------------------------------------------------------------------
 // 
  Xoldpar=Xpar;
  Xpar=rankselect(Xpar); // rank selection
  for i=1:2:(Npop-1)
    u1=int(Npop*rand())+1;u2=int(Npop*rand())+1;u3=rand();
    xpar=[Xpar(u1,:);Xpar(u2,:)];Xoff(i:(i+1),:)=xpar;
    if (u3<pc) then
     xoff=blcross(xpar);   //crossover
     Xoff(i:(i+1),:)=xoff;
    for j=1:2
       u4=rand();
        if (u4<pm) then
          xoffl=mutate(xoff(j,:),Nevalf,Nevalfmax,xmin,xmax);           Xoff(i+j-1,:)=xoffl;
        end
    end
    end
  end
//
//------------------------------------------------------------------------------
//
//     MAIN LOOP : ENDING
//
end
//-----------------------------------------------------
//
//[Xpar,bestpar]=evalplus1elitism(Xoff,Xoldpar,Npop,gen);  // evaluation plus 1-elitism
fmin=[fmin,bestpar($)];
Ntot=[Ntot,Nevalf];
result(Ntot,fmin,1);                      // dynamic result display
disp('minimum obtained:');disp(bestpar(1:n));
disp('corresponding value by f:');disp(fmin($));
disp('evaluation number:');disp(Nevalf);
//

