Npop=60; L=2; Ngen=100; pc=0.6; pm=0.2;xmin=-5;xmax=5;
sigma=(xmax-xmin)/100;
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=[];oldbestval=0;
for gen=1:Ngen
// evaluation  + elitism  //////////////////////////////////////
for i=1:Npop
  A(i,L+1)=J(A(i,1:L));
end
//  elitism
[bestval,pos]=min(A(:,L+1)); if (bestval>oldbestval)&(gen>1) then
    u=int(Npop*rand()+1); A(u,:)=bestelem;
end
// selection
// 1st step: sort elements
[s,p]=sort(A(:,L+1)); A=A(p,:); val=[val,A(Npop,L+1)];
bestelem=A(Npop,:); oldbestval=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]=sort(A(:,L+1)); A=A(p,:); disp(A(Npop,1:L)); plot2d(val)
