// Differential Evolution
// Parametri:
//   m:  numarul elemente in populatie
//   n: dimensiunea fiecarui element
//   gen:  numar de generatii
//   [a,b]: domeniul de definitie al functiei de optimizat (minimizat)
//   F: factor de scalare
//   CR: probabilitate incrucisare
//   selectie: numele functiei de selectie
//   functie: numele functiei de minimizat
// 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); //initializarea populatiei
    pop2=zeros(m,n);     //alocare spatiu pentru populatia de copii obtinuta dupa mutatie si incrucisare
    f=zeros(1,m);             //alocare spatiu pentru valorile elementelor populatiei 
    fnew=zeros(1,m);
    for i=1:m
         f(i)=functie(pop1(i,:));     
    end
    [valmin,imin]=min(f);
    for k=1:gen
        // generare populatie de urmasi prin recombinare si mutatie  
        for i=1:m
            pop2(i,:)=DE_Rand1Bin(i,pop1,F,CR);
            //pop2(i,:)=DE_Best1Bin(i,imin,pop1,F,CR);
            f2(i)=functie(pop2(i,:));
        end
        for i=1:m
            j=closest(pop2(i,:),pop1)
            if f2(i)<f(j) then
                pop1(j,:)=pop2(i,:)
                f(j)=f2(i)
            end
        end
        [valmin,imin]=min(f);
        disp(valmin);  // afisarea valorii celui mai bun element din populatie
     end 
     // reprezentarea grafica a functiei si a elementelor populatiei
     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

function d=distance(x,y)
    n=size(x)
    d=(x-y)*(x-y)'
endfunction

function imin=closest(x,pop)
    m=size(pop,1)
    imin=1
    dmin=distance(x,pop(imin,:))
    for i=2:m
        d=distance(x,pop(i,:))
        if  d<dmin then
            dmin=d
            imin=i
        end
    end
endfunction

// Differential Evolution:  initializarea populatiei (m vectori avand cate n componente)
// fiecare componenta se genereaza uniform aleator in intervalul [a,b]
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: construire candidat nou
// Pas 1:  Mutatie:   trialVector=pop(r1)+F*(pop(r2)-pop(r3))  
// Pas 2: Incrucisare: 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));  // 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

// Differential Evolution: construire candidat nou
// Pas 1:  Mutatie:   trialVector=pop(ibest)+F*(pop(r2)-pop(r3))  
// Pas 2: Incrucisare: 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



// Strategie evolutiva:  functie test (functia patratica)
function y=DE_f1(x)
    n=size(x,2);
    y=0;
    for j=1:n
        y=y+x(j)^2;
    end
endfunction

// Strategie evolutiva: functie test multi-modala (multe minime locale)
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
// Parametri:
//   m:  numarul elemente in populatie
//   n: dimensiunea fiecarui element
//   gen:  numar de generatii
//   [a,b]: domeniul de definitie al functiei de optimizat (minimizat)
//   F: factor de scalare
//   CR: probabilitate incrucisare
//   selectie: numele functiei de selectie
//   functie: numele functiei de minimizat
// Exemplu de apel: DE(10,1,100,-50,50,0.5,0.55,DE_f1)
DE(30,1,500,-50,50,0.5,0.5,DE_f2)
