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

function [pop1,f]=DE(m,n,gen,a,b,F,CR,functie)
    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,:));     
    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,:));
            if (f2<f(i)) then 
                pop1(i,:)=pop2(i,:);
                f(i)=f2;
            end
        end
        [valmin,imin]=min(f);
        disp(valmin);  // print the value of the best element 
     end 
     // plot
     if (n==1) then
       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*");
     end
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) with probability 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));  // mutation
    for j=1:n                                    // binomial crossover
        u=grand(1,1,"unf",0,1);
        if (u<1-CR) then
            trialVector(1,j)=pop(i,j);
        end
    end
    disp
endfunction
// Sphere
function y=DE_f1(x)
    n=size(x,2);
    y=0;
    for j=1:n
        y=y+x(j)^2;
    end
endfunction
// Griewank
function y=DE_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

// Differential Evolution
// Parameters:
//   m:  population size
//   n: problem size
//   gen:  numar de generatii
//   [a,b]: bounds
//   F: scale factor
//   CR: crossover probability
//   function: function to be minimized
// Exemplu de apel: DE(10,1,100,-50,50,0.5,0.55,DE_f1)
DE(10,10,50,-50,50,0.5,0.5,DE_f2)
