//////////////////////////////////////////////
//
// Basic real GA
/////////////////////////////////////////////
Npop=50;
Ngen=30;
pc=0.6;
pm=0.2;
xmin=-5;
xmax=5;
sigma=(xmax-xmin)/100;
//
L=2;  // number of parameters
/////////////////////////////////////////////////
function y=J(x)   //
   y=sum(x.^2-cos(2*%pi*x))+length(x);  // rastrigin
endfunction
// initialisation ///////////////////////////////////
A= xmin+(xmax-xmin)*rand(Npop,L+1);
//
val=[];
for gen=1:Ngen
// evaluation   //////////////////////////////////////
for i=1:Npop
  A(i,L+1)=J(A(i,1:L));
end
[s,p]=gsort(A(:,L+1));
A=A(p,:);
val=[val,A(Npop,L+1)];
// 2st step: choose Npop elements with probability pi
//
Asel=[]; p=(1:Npop)/sum(1:Npop);ps=cumsum(p); for i=1:Npop
  u=rand();isel=1;
  while (u>ps(isel))
     isel=isel+1;
  end
Asel=[Asel;A(isel,:)]; end
// barycentric crossover
Acrois=[];
   for k=1:Npop/2
   u1=int(Npop*rand()+1);
   u2=int(Npop*rand()+1);
    if (rand()<pc) then
      alpha=rand();
      uenf1=alpha*Asel(u1,:)+(1-alpha)*Asel(u2,:);
      uenf2=(1-alpha)*Asel(u1,:)+alpha*Asel(u2,:);
    else
      uenf1=Asel(u1,1:L+1);uenf2=Asel(u2,1:L+1);
    end
    Acrois=[Acrois;uenf1;uenf2];
   end
// normal mutation
Amut=Acrois;
  for k=1:Npop
      if (rand()<pm) then
         epsilon=sigma*rand(1,L+1,'normal');Amut(k,:)=Amut(k,:)+epsilon;
      end
  end
A=Amut; end
 for i=1:Npop
  A(i,L+1)=J(A(i,1:L));
end;
[s,p]=gsort(A(:,L+1)); A=A(p,:); disp(A(Npop,1:L)); 
////////// history plot///////////////
clf()
plot2d(Npop:Npop:Ngen*Npop,val)
