// Evolution strategy (mu+lambda) - minimization problems
// Input data:
//   mu:  number of parents
//   lambda: number of children 
//   n: problem size (number of components of each population element)
//   gen:  number of generations
//   [a,b]: search range (domain of the objective function)
//   sigma: parameter controlling the amount of perturbation (standard deviation of the Gaussian distribution)
//   rho: number of parents used in the recombination
//   selectie: name of the selection function
//   functie: name of the objective function
// Call: SE(10,20,1,100,-50,50,0.1,5,@SE_trunchiere,@SE_f1)

function [pop,f]=SE(mu,lambda,n,gen,a,b,sigma,rho,selectie,functie)
    pop1=SE_initpop(mu,n,a,b); //population initialization
    pop2=zeros(lambda,n); //space for the population of children
    m=mu+lambda;
    f=zeros(1,m);      //space for the objective values of all elements (current population + children)
    for k=1:gen
        // construction of the children (recombination+mutation)  
        for i=1:lambda
            pop2(i,:)=SE_recombinare(pop1,rho);
            pop2(i,:)=SE_mutatie(pop2(i,:),sigma);
        end
        pop=[pop1;pop2];  // concatenation of all elements
        // evaluation
        for i=1:m
            f(1,i)=functie(pop(i,:));     
        end   
        // selection of survivors
        pop1=selectie(pop,f,mu); 
        disp(min(f));  // the value of the best element in the population
     end 
     disp(min(f)); 
     f=[];
     for i=1:mu
         f(i)=functie(pop1(i,:));     
     end
     // plot of the objective function + population
     x=a:0.1:b; y=zeros(1,size(x,2));
     for j=1:size(x,2)
         y(j)=functie(x(j));
     end   
     plot(x,y,"b",pop1',f,"r*");
endfunction

// Population initialization (uniformly random in [a,b])
function pop = SE_initpop(m,n,a,b)
    pop=zeros(m,n);
    for i=1:m
        for j=1:n
            pop(i,j)=grand(1,1,'unf',a,b);
        end   
    end   
endfunction

// Convex (intermediate) recombination
function copil=SE_recombinare(p, rho)
    m=size(p,1); 
    n=size(p,2);
    copil=zeros(1,n);
    for k=1:rho
        ind=grand(1,1,"uin",1,m);  // uniform selection of parents
        copil=copil+p(ind,:); 
    end
    copil=copil/rho;
endfunction

// Mutation using a N(0,sigma) perturbation 
function copil=SE_mutatie(parinte,sigma)
    n=size(parinte,2);
    copil=parinte;
    for i=1:n
        copil(i)=copil(i)+grand(1,1,"nor",0,sigma);
    end   
endfunction

// Tournament selection (survivors)
function p2=SE_turneu(p1,f1,mu)
m=size(p1,1);
p2=[];
sel=zeros(1,m);
for i=1:mu
   j1=grand(1,1,'uin',1,m);
   j2=grand(1,1,'uin',1,m);
   if f1(j1)<f1(j2)
      p2(i,:)=p1(j1,:);
   else   
      p2(i,:)=p1(j2,:);  
   end   
end   
endfunction

// Truncation selection (survivors)
function p2=SE_trunchiere(p1,f1,mu)
m=size(p1,1); 
n=size(p1,2);
p2=zeros(mu,n);
sel=zeros(1,m); 
for i=1:mu
   jmin=1;
   for j=2:m
      if sel(j)==0 
         if f1(j)<f1(jmin)
            jmin=j;
         end   
      end
   end   
   sel(jmin)=1;
   p2(i,:)=p1(jmin,:);
end  
endfunction

// Test function: sphere
function y=SE_f1(x)
    n=size(x,2);
    y=0;
    for j=1:n
        y=y+x(j)^2;
    end
endfunction

// Test function: Griewank
function y=SE_f2(x)
    n=size(x,2);
    s=0;p=1;
    for j=1:n
        s=s+x(j)^2;
        p=p*cos(x(j)/sqrt(j));   
    end   
    y=s/4000-p+1;
endfunction

// Test function: Ackley
function y = ackley(x)
n=size(x,2);
a = 20; b = 0.2; c = 2*%pi;
s1 = 0; s2 = 0;
for i=1:n;
   s1 = s1+x(i)^2;
   s2 = s2+cos(c*x(i));
end
y = -a*exp(-b*sqrt(1/n*s1))-exp(1/n*s2)+a+exp(1);

endfunction

SE(10,20,1,100,-50,50,0.1,5,SE_turneu,SE_f2)
