// Genetic algorithm for the knapsack problem
// Input data:
//    m - population size (ex: 50)
//    n - number of objects (ex: 5)
//    gen - number of generations (ex: 100)
//    pm - mutation probability (ex: 0.01)
//    alfa - penalty coefficient used in constructing the fitness function  (ex: 0.8)
//    V - values of the objects (ex: [3,4,1,5,2])
//    W - weights of the objects (ex: [10,8,3,6,5])
//    C - knapsack capacity (ex: 20)
// Call:  knapsack(50,5,100,0.01,0.8, [3,4,1,5,2],[10,8,3,6,5],20); 
function pop=knapsack(m,n,gen,pm,alfa,V,W,C)
    pop=initpop(m,n); 
    for k=1:gen       
        copii=zeros(m,n);       
//        val=evalpop(pop,alfa,V, W, C);
//        prob=calculprob(val);          
//        indici1=ruleta(prob);          
//        indici2=ruleta(prob);
        for i=1:m/2
          index1=grand(1, 1, "uin", 1, m);
          index2=grand(1, 1, "uin", 1, m);
          if (better(pop(index1,:)',pop(index2,:)')==1) then parinte1=pop(index1,:);
          else parinte1=pop(index2,:);
          end    
          index1=grand(1, 1, "uin", 1, m);
          index2=grand(1, 1, "uin", 1, m);                      
          if (better(pop(index1,:)',pop(index2,:)')==1) then parinte2=pop(index1,:);
          else parinte2=pop(index2,:);
          end  
          [copil1,copil2]=incrucisare(parinte1,parinte2);
          copii(i,:)=copil1;
          copii(m-i+1,:)=copil2;
        end  
        for i=1:m                      
          if grand(1,1,"unf",0,1)<pm
             copii(i,:)=mutatie(copii(i,:));
          end
        end   
        pop=copii; disp(pop);         
    end   
    indBest=findBest(val);
    disp(pop(indBest,:),"Best element:");
    valBest=sum(V.*pop(indBest,:));
    weightBest=sum(W.*pop(indBest,:));
    disp(valBest,"Value of the selected objects:"); 
    disp(weightBest,"Weight of the selected objects:"); 
endfunction

function pop = initpop(m,n)
    for i=1:m
        for j=1:n
            pop(i,j)=grand(1,1,'uin',0,1);
        end
    end    
endfunction

function val=evalpop(pop,alfa,V, W, C)
[m,n]=size(pop);
val=[];
for i=1:m
   vrucsac=0; wrucsac=0;
   for j=1:n
      vrucsac=vrucsac+V(j)*pop(i,j);
      wrucsac=wrucsac+W(j)*pop(i,j);
   end   
   if wrucsac<=C
      val(i)=vrucsac;   
      // the constraint is satisfied
   else
      val(i)=(1-alfa)*vrucsac+alfa*(C-wrucsac); 
      // the constraints is not satisfied
   end   
end   
endfunction

function rez=better(x,y,V,W,C)
    cx=C-W*x
    cy=C-W*x
    vx=V*x
    vy=V*y
    if (cx>=0)&(cy>=0) then
        if (vx>=vy) then
           rez=1
        else rez=0    
        end    
    end
    if (cx>=0)&(cy<0) then
        rez=1
    end
    if (cx<0)&(cy>=0) then
        rez=0
    end
    if (cx<0)&(cy<0) then
        if (cx<cy) then rez=0
        else rez=1
        end    
    end
endfunction

function indBest=findBest(val)
    indBest=1;
    for i=2:length(val)
        if val(i)>val(indBest)
            indBest=i;
        end
    end
endfunction    

function prob=calculprob(val)
m=length(val);
prob=zeros(m);
s=0;
for i=1:m
   s=s+val(i);
end
for i=1:m
   prob(i)=val(i)/s;
end     
endfunction

function indici=ruleta(prob)
    indici=zeros(length(prob));
    for k=1:length(prob)
        i=1;
        s=prob(1);
        u=grand(1,1,"unf",0,1);
        while u>s
           i=i+1;
           s=s+prob(i);
        end
        indici(k)=i;
    end
endfunction

// crossover
function [copil1,copil2]=incrucisare(parinte1, parinte2)
    n=size(parinte1,2); 
    copil1=zeros(1,n);copil2=zeros(1,n);
    poz=grand(1,1,"uin",1,n-1); // generate the cut point
    for i=1:poz                 // fill in the first components of the offsprings   
      copil1(1,i)=parinte1(1,i);
      copil2(1,i)=parinte2(1,i);
   end
   for i=poz+1:n            // fill in the last components of the offsprings
      copil1(1,i)=parinte2(1,i);
      copil2(1,i)=parinte1(1,i);
   end
endfunction

function copil=mutatie(parinte)
    n=size(parinte,2);
    copil=parinte;
    poz=grand(1,1,"uin",1,n);
    copil(poz)=1-copil(poz);
endfunction
knapsack(50,5,100,0.01,0.8, [3,4,1,5,2],[10,8,3,6,5],20);


