//
// multi objective optimization
// NSGA method without sharing
// comparaison of the exact and computed Pareto front on a simple example.
//
//////////////////////////////////////////////////////
function [J1,J2]=MO(x)   // 
   J1=(x(1)-1)^2+(x(2)-3)^2;
   J2=(x(1)-4)^2+(x(2)-2)^2;
endfunction 
/////////////////////////////////////////////////////
//
// calcul des fronts de Pareto et tracés
//////////////////////////////////////////////////////
//
function Xp=paretoplot(Pop,MOpop)
[Npop,b]=size(Pop);
Xp=[];JXp=[];
ind=(1:Npop);
count=1;
while(length(ind)>0)
Nfront=[];
  for i=ind;
  l=find(MOpop(ind,1)<=MOpop(i,1) & MOpop(ind,2)<=MOpop(i,2))  
  if (length(l)==1) then
    Xp=[Xp;[Pop(i,:),count]];
    JXp=[JXp;[MOpop(i,:),count]];
  else
    Nfront=[Nfront,i];
  end
end
ind=Nfront;
// exclude the current front
count=count+1;
end
clf()
for i=1:(count-1)
  ind=find(JXp(:,3)==i)
  [S1,S2]=gsort(JXp(ind,1));
  plot2d(JXp(ind(S2),1),JXp(ind(S2),2),i+1,rect=[0,0,30,30]);
  plot2d(JXp(ind(S2),1),JXp(ind(S2),2),-1,rect=[0,0,30,30]);
end
//
// front de Pareto exact
//
lambda=0:0.01:1;
J1p=10*(lambda-1)^2;
J2p=10*lambda^2;
plot2d(J1p, J2p)
endfunction
//
//////////////////////////////////////
// programme principal
//////////////////////////////////////
Npop=40;
L=2;                                  
Ngen=100;
pc=0.6;
pm=0.2;
xmin=-5;xmax=5;                      
sigma=(xmax-xmin)/100;
//
A= xmin+(xmax-xmin)*rand(Npop,L+1);
//
for gen=1:Ngen
//
// evaluation et calcul front de Pareto ///////////////
MOpop=[];
for i=1:Npop
[J1,J2]=MO(A(i,1:L));
MOpop=[MOpop;[J1,J2]];
end
drawlater
A=paretoplot(A(:,1:L),MOpop+1E-10*rand(Npop,2));
drawnow
sleep(1000)
//
///////////////////////////////////////////////////////
// selection
[s,p]=gsort(A(:,L+1));
A=A(p,:);
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  
//


