// Evolution strategy (mu+lambda)
//
// Parameters:
//   mu:  number of parents
//   lambda: number of children
//   n: problem size (number of componente of each population elements)
//   gen:  number of generation
//   [a,b]: range of values for the objective function variables
//          (in the case of a neural network is the initial range of the synaptic weights which is usually [-1,1])
//   sigma: perturbation parameter (standard deviation on the gaussian distribution used for random perturbation)
//   rho: the number of parents used in recombination
//   selection: selection function
//   function: objective function (in the case of a NN it is MSE - mean squared error corresponding to the training set)
//   nrHidden: number of hidden units
// Call example: SE_nn(20,20,100,-1,1,0.5,2,@SE_truncation,@SE_XOR,5);

function [pop,f,k,errMin]=SE_nn(mu,lambda,gen,a,b,sigma,rho,selection,objFunction,nrHidden)
n=4*nrHidden+1;
pop1=SE_initpop(mu,n,a,b); //population initialization
pop2=zeros(lambda,n); // initialization of offspring population (just for space allocation)
m=mu+lambda;
f=zeros(1,m);      //initialization of fitness values (just for space allocation)
for k=1:gen
   // creation of offspring by recombination and mutation  
   for i=1:lambda
      pop2(i,:)=SE_recombination(pop1,rho);
      pop2(i,:)=SE_mutation(pop2(i,:),sigma);
   end
   pop=[pop1;pop2];  // concatenation of the current population and the offspring population
   // evaluation of the population elements
   for i=1:m
      f(1,i)=objFunction(pop(i,:),nrHidden);     
   end   
   // selection of survivors
   pop1=selection(pop,f,mu); 
   disp(min(f));  
   if (min(f)==0)
       break;
   end;    
end 
 disp("Nr generations "); disp(k); 
 errMin=min(f);
 disp('Error:'); disp(errMin); 
endfunction

// ES:  population initialization (m vectors with n components)
// uniformly distributed random value 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

// ES: recombination
// rho randomly selected parents -> average
function child=SE_recombination(p, rho)
m=size(p,1); 
n=size(p,2);
child=zeros(1,n);
for k=1:rho
   ind=grand(1,1,"uin",1,m);  // selectare indice parinte
   child=child+p(ind,:); 
end
child=child/rho;
endfunction

// ES: mutation
// generate a new alement from the parent by applying a perturbation with a N(0,sigma) distribution
function child=SE_mutation(parent,sigma)
n=size(parent,2);
child=parent;
for i=1:n
   child(i)=child(i)+grand(1,1,"nor",0,sigma);
end   
endfunction

// tournament selection
function p2=SE_tournament(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
function p2=SE_truncation(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


function v=H(u)
if u<=0
    v=0;
else
    v=1;
end
endfunction

function err=SE_XOR(x,k)
b=[-1 -1 -1 -1; 0 0 1 1; 0 1 0 1];
d=[0 1 1 0];
w1=zeros(k,3); w2=zeros(1,k);
n=size(x,2); disp(n,"n=");
w1(:,1)=x(1:k)'; disp(w1(:,1),"w1=");
w1(:,2)=x(k+1:2*k)'; disp(w1(:,2),"w2=");
w1(:,3)=x(2*k+1:(3*k))'; disp(w1(:,3),"w3=");
w2=x(3*k+1:(4*k+1));
err=0;
for j=1:4
   x1=w1*b(:,j);
   for i=1:k
       y1(i)=H(x1(i));
   end    
//   disp(y1,"y1=");
   x2=[-1 y1']';
   y=H(w2*x2);
//   disp('y');
//   disp(y);
//   disp('d');
//   disp(d(j));
   err=err+(d(j)-y)^2;
end
disp(err,"err=");
endfunction
// SE_nn(mu,lambda,gen,a,b,sigma,rho,selection,function,nrHidden)

mean=0;
nrsucc=0;
for i=1:10
    [pop,f,k,errMin]=SE_nn(20,20,100,-1,1,0.5,2,SE_truncation,SE_XOR,1);
    if errMin==0
        nrsucc=nrsucc+1;
        mean=mean+k;
    end    
end
disp('nrsucc=',nrsucc);
disp('mean=',mean/nrsucc);











