function [pop_natural,pop_costs] = pareto_GA( wfn, bounds, popsize, selection, mutrate, maxit )
% PARETO_GA - Implements a pareto genetic algorithm for multiple objective minimization (MOO) for *two* functions only
% 
% Inputs: 
%   wfn = a function handle accepting a double vector input and returning two columns of scalar output representing f_1 and f_2 function to minimize
%   bounds = [ number_of_variables x 2 ] input matrix where the ith row is the valid bounds for the ith variable
%   
% Outputs:
%   x_bin = [ number_of_samples x npar ] pareto encoding of each variable in x
% 
% Written by:
% -- 
% John L. Weatherwax                2006-08-28
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

% Parse the parameters of the genetic algorithm (if not provided defaults are used)
% 
if( nargin < 3 ) 
  popsize = 100; % the population size
end
if( nargin < 4 )
  selection = 0.5; % selection rate 
end
if( nargin < 5 ) 
  mutrate = 0.02; % the mutation rate
end
if( nargin < 6 )
  maxit = 20; 
end

% Check inputs:
npar = size(bounds,1); 
assert( size(bounds,2)==2, 'bounds input must have two columns' );

keep = floor( selection * popsize );

M = ceil( (popsize-keep)/2 ); % the number of matings 
odds = 1;
for ii=2:keep,
  odds = [odds ii*ones(1,ii)];
end
odds = keep+1-odds;
Nodds = length(odds);

% _________________________________________________________________________
% Create the initial population 
%
iga         = 0;                      % generation counter initialized 
pop         = rand(popsize,npar);     % random population of normalized continious variables
pop_natural = var_decode(pop,bounds); % Decode normalized variables into natural variables (the variables from the problem definition)
fout        = wfn( pop_natural );     % Find cost for each variable

%
% _________________________________________________________________________
while iga<maxit
  iga = iga+1; % increments generation counter
  cost = zeros(popsize,1);
  [g,ind] = sort(fout(:,1)); % sort on the first column or g = f_1
  pop = pop(ind,:); % sort chromosomes based on best values of f_1
  h = fout(ind,2); % extract values of f_2 ordered by thoes of f_1 
  ct=0; rank=1; q=0;
  while ct<popsize
    for ig=1:popsize
      if h(ig)<=min(h(1:ig))
        ct=ct+1;
        q(ct)=ig;
        cost(ct,1) = rank;
      end
      if rank==1 % save the most dominant results
        px=g(q); py=h(q); elite=length(q);
      end
      if ct==popsize; break; end
    end
    rank=rank+1;
  end
  pop=pop(q,:); % order the new population by the elite chromosomes
  
  %figure(1); clf; % plot new each time 
  figure(gcf); hold on; % plot on top of the results from problem 2 
  plot(px,py,'-g'); plot(g,h,'.k','MarkerSize',4); % plot (f_1,f_2) and the most dominant results 
  xlabel('objective function f_1'); ylabel('objective function f_2'); title('pareto front'); 
  %pause
  
  [cost,ind] = sort(cost);
  pop = pop(ind,:);
  
  % _________________________________________________________________________
  % tournament selection for determinining parents: 
  Ntourn=2;
  picks  = ceil(keep*rand(Ntourn,M));
  [c,pt] = min(cost(picks));
  for ib=1:M
    ma(ib)=picks(pt(ib),ib); % the mother 
  end
  picks  = ceil(keep*rand(Ntourn,M));
  [c,pt] = min(cost(picks));
  for ib=1:M
    pa(ib)=picks(pt(ib),ib); % the father 
  end
  
  % _________________________________________________________________________
  % Perform mating given the mother and father: 
  ix=1:2:keep; % index of mate #1
  xp=floor(rand(1,M)*npar); % crossover point
  r=rand(1,M); % mixing parameter
  xy = pop(ma+popsize*xp) - pop(pa+popsize*xp);
  % mix from the mother and the father 
  pop(keep+ix+popsize*xp) = pop(ma+popsize*xp) - r.*xy;
  % first offspring 
  pop(keep+ix+1+popsize*xp) = pop(pa+popsize*xp) + r.*xy;
  % second offspring
  for ic=1:M/2
    if xp(ic)<npar % perform crossover when last variable not selected
      pop(keep+ix(ic),:)   = [pop(ma(ic),1:xp(ic)) pop(pa(ic),xp(ic)+1:npar)];
      pop(keep+ix(ic)+1,:) = [pop(pa(ic),1:xp(ic)) pop(ma(ic),xp(ic)+1:npar)];
    end % if
  end % end ic
  
  %disp('the top 4 samples');
  %pop(1:4,:)
  
  % _________________________________________________________________________
  % Mutate the population 
  nmut = ceil((popsize-elite)*npar*mutrate); % total number of mutations 
  mrow = ceil(rand(1,nmut)*(popsize-elite));
  mcol = ceil(rand(1,nmut)*npar);
  mutindex = mrow+(mcol-1)*popsize;
  pop(mutindex) = rand(1,nmut);
  
  % _________________________________________________________________________
  % The new offspring and muted chromosomes are evaluated for cost
  row = sort(rem(mutindex,popsize));
  iq=1; rowmut(iq) = row(1);
  for ic=2:nmut
    if row(ic)>keep; break; end;
    if row(ic)>rowmut(iq)
      iq=iq+1; rowmut(ic)=row(ic);
    end
  end
  if rowmut(1)==0; rowmut=rowmut(2:length(rowmut)); end;
  fout(rowmut,:)         = wfn( var_decode(pop(rowmut,:),bounds) );
  fout(keep+1:popsize,:) = wfn( var_decode(pop(keep+1:popsize,:),bounds) ); 

  % _________________________________________________________________________
  % Stopping criteria
  if iga>maxit
    break;
  end
  
  %disp('the optimal solution at this timestep')
  %[iga cost(1) fout(1,:)]
  
end % iga

pop_natural = var_decode(pop,bounds);
pop_costs   = wfn( pop_natural );