// In mathematical optimization, the Rastrigin function is a non-convex function 
// used as a performance test problem for optimization algorithms
//------------------------------------------------------
// real valued genetic algorithm 
// for the minimization of the rastrigin function 
// (with dynamic display for the case of 2 parameters)
//
//-----------------------------------------------------
// selection: proportionate with ranking method 
// crossover: blend 
// mutation: non uniform
// 1-elitism
//
function y=rastrigin(x,n)  // the function to optimize
y=n+sum(x.^2-cos(2*%pi*x));
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),n);
end
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]=gsort(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,:)=u*xpar(1,:)+(1-u)*xpar(2,:);
   xoff(2,:)=u*xpar(2,:)+(1-u)*xpar(1,:);
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
//
//---------main program ---------------------
//
clf()
clear all
disp('real valued GA for minimizing the rastrigin function');
disp('(with dynamic display for 2 parameters)');
//
xmin=-5;xmax=5;N=100;
x=xmin:((xmax-xmin)/(N-1)):xmax;
y=x;
z=zeros(N,N);
for i=1:N
for j=1:N
z(i,j)=rastrigin([x(i),y(j)],2);
end
end
//
n=evstr(x_dialog('parameter number of the rastrigin function','2')); 
Npop=evstr(x_dialog('population number','30'));     
Ngen=evstr(x_dialog('generation number','50')); 
pc=evstr(x_dialog('crossover probability','0.9')); 
pm=evstr(x_dialog('mutation probability','0.6')); 
//
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,zeros(Npop,1)];  
Xoff=Xpar;Xoldpar=Xpar;fmin=[];mineval=0;newval=zeros(Npop,1);Traj=[];
//
xset('window',2);
xset('window',3);
//
for gen=1:Ngen;
  Xpar=Xoff;
  [Xpar,bestpar]=evalplus1elitism(Xpar,Xoldpar,Npop,gen);  // evaluation plus 1-elitism
  Xoldpar=Xpar;
  fmin=[fmin,bestpar($)];
  //
  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;
    end
    for j=1:2
       u4=rand();
        if (u4<pm) then
          xoffl=mutate(xoff(j,:),gen,Ngen,xmin,xmax);  // mutation
          Xoff(i+j-1,:)=xoffl;
        end
    end
   end
//-----  dyanmic display if n=2 ----------------------------
   if (n==2) then
   xset('window',2);
   drawlater
   clf()
   Traj=[Traj;bestpar(1),bestpar(2)];
   contour2d(x,y,z,[0:0.01:0.1,0.2:1,1:10],rect=[-5.12,-5.12,5.12,5.12]);
   plot2d(Xpar(:,1),Xpar(:,2),-1,rect=[-5.12,-5.12,5.12,5.12]);  
   plot2d(Traj(:,1),Traj(:,2),5,rect=[-5.12,-5.12,5.12,5.12]);
   xtitle('trajectory display');
   drawnow
   xset('window',3);
   drawlater
   clf()
   contour2d(x,y,z,[0:0.01:0.1,0.2:1,1:10],rect=[-0.5,-0.5,0.5,0.5]);
   plot2d(Xpar(:,1),Xpar(:,2),-1,rect=[-0.5,-0.5,0.5,0.5]);   
   plot2d(Traj(:,1),Traj(:,2),5,rect=[-0.5,-0.5,0.5,0.5]);   
   xtitle('trajectory display');
   drawnow
   end
//-----  end of dyanmic display   ---------------------------
//
end
[Xpar,bestpar]=evalplus1elitism(Xoff,Xoldpar,Npop,gen);  // evaluation plus 1-elitism
fmin=[fmin,bestpar($)];
//
//----------- results displays --------------
//
xset('window',0);
clf();
plot2d(0:Npop:(Npop*Ngen),fmin);
xtitle('convergence history','Neval','fmin');
xset('window',1);
clf();
plot2d1('onl',0:Npop:(Npop*Ngen),fmin);
xtitle('convergence history (log scale)','Neval','fmin');
//
disp('minimum obtained:');disp(bestpar(1:n));
disp('corresponding value by f:');disp(bestpar($));
disp('evaluation number:');disp(Npop*Ngen);















