// Differential Evolution
// Parameters:
//   m:  population size
//   gen:  numar de generatii
//   [a,b]: bounds
//   F: scale factor
//   CR: crossover probability
//   functionn: function to be minimized
//  K= number of hidden units
// Exemplu de apel: DE(10,1,100,-50,50,0.5,0.55,SE_f1)

function [pop1,f,valmin,k]=DE(m,gen,a,b,F,CR,functie,K)
    n=4*K+1;
    pop1=DE_initpop(m,n,a,b); //initialization
    pop2=zeros(m,n);     //space allocation for the population of offsprings
    f=zeros(1,m);        //space allocation for the function values 
    for i=1:m
         f(i)=functie(pop1(i,:),K);     
    end
    [valmin,imin]=min(f);
    for k=1:gen
        // applying mutation and crossover 
        for i=1:m
            //pop2(i,:)=DE_Rand1Bin(i,pop1,F,CR);
            pop2(i,:)=DE_Best1Bin(i,imin,pop1,F,CR);
            f2=functie(pop2(i,:),K);
            if (f2<f(i)) then 
                pop1(i,:)=pop2(i,:);
                f(i)=f2;
            end
        end
        [valmin,imin]=min(f);
        if (valmin==0) then break
        end
    end    
    disp(valmin,"errmin=");  // print the value of the best element 
    disp(k,"number of generations=")
endfunction

// initialization of the population 
function pop = DE_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

// Differential Evolution: constructing a new candidate
// Step 1:  Mutation:   trialVector=pop(r1)+F*(pop(r2)-pop(r3))  
// Step 2: Crossover: trialVector(j)=pop(i,j) cu probabilitatea 1-CR
function trialVector=DE_Rand1Bin(i,pop, F, CR)
    m=size(pop,1); 
    n=size(pop,2);
    trialVector=zeros(1,n);
    r1=grand(1,1,"uin",1,m);
    r2=grand(1,1,"uin",1,m); 
    while(r2==r1)
        r2=grand(1,1,"uin",1,m);
    end
    r3=grand(1,1,"uin",1,m); 
    while((r3==r1)|(r3==r2))
        r3=grand(1,1,"uin",1,m);
    end    
    trialVector=pop(r1,:)+F*(pop(r2,:)-pop(r3,:)); 
    for j=1:n
        u=grand(1,1,"unf",0,1);
        if (u<1-CR) then
            trialVector(1,j)=pop(i,j);
        end
    end
endfunction

// DE/best/1/bin
// Step 1:  Mutation:   trialVector=pop(ibest)+F*(pop(r2)-pop(r3))  
// Step 2: Crossover: trialVector(j)=pop(i,j) cu probabilitatea 1-CR
function trialVector=DE_Best1Bin(i,ibest,pop, F, CR)
    m=size(pop,1); 
    n=size(pop,2);
    trialVector=zeros(1,n);
    r1=grand(1,1,"uin",1,m);
    r2=grand(1,1,"uin",1,m); 
    while(r2==r1)
        r2=grand(1,1,"uin",1,m);
    end
    trialVector=pop(ibest,:)+F*(pop(r1,:)-pop(r2,:));  // selectare indice parinte
    for j=1:n
        u=grand(1,1,"unf",0,1);
        if (u<1-CR) then
            trialVector(1,j)=pop(i,j);
        end
    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

// Differential Evolution
// Parameters:
//   m:  population size
//   gen:  numar de generatii
//   [a,b]: bounds
//   F: scale factor
//   CR: crossover probability
//   functionn: function to be minimized
//  K: number of hidden units
// call: DE(10,1,100,-50,50,0.5,0.55,SE_f1)
DE(10,50,-50,50,0.5,0.5,SE_XOR,3)

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