%
%--------------------------------------------------------------
% ECP2012: A real Genetic Algorithm starting from the historical one
%---------------------------------------------------------------
Npop=200;
L=6;
Ngen=200;
pc=0.6;
pm=0.2;
sigma=0.1;
%
% initialisation
A=-5+10*rand(Npop,L+1);
%
for i=1:Ngen
% evaluation
for i=1:Npop
  A(i,L+1)=JJ(A(i,1:L));
end
% selection
% 1st step: sort elements
[s,p]=sort(-A(:,L+1));
A=A(p,:);
% 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  
% crossover  
Acrois=[];
  for k=1:Npop/2
    u1=floor(Npop*rand()+1);
    u2=floor(Npop*rand()+1);
    alpha=rand();
    if (rand()<pc)
      uenf1=alpha*Asel(u1,:)+(1-alpha)*Asel(u2,:);  
      uenf2=alpha*Asel(u2,:)+(1-alpha)*Asel(u1,:);  
    else
      uenf1=Asel(u1,:);
      uenf2=Asel(u2,:);
    end   
    Acrois=[Acrois;uenf1;uenf2];
   end 
% mutation   
Amut=Acrois;
  for k=1:Npop
      if (rand()<pm) 
         Amut(k,:)=Amut(k,:)+sigma*randn(1,L+1);
    end
  end
A=Amut;
end  
%
for i=1:Npop
  A(i,L+1)=JJ(A(i,1:L));
end
[s,p]=sort(-A(:,L+1));
A=A(p,:);
disp('best point:');
disp(A(Npop,1:L))
disp('corresponding value:');
disp(A(Npop,L+1))


