// PSO
// Parameters:
//   m:  number of elements (particles)
//   n: size of each element (particles)
//   gen:  number of generations
//   [a,b]: domain of variables
//   gamma: constriction factor
//   r1: random factor
//   r2: random factor
//   functie: function to be minimized
// Exemplu de apel:[xbest,fbest]=PSO(10,1,100,-50,50,0.1,0.5,0.5,@PSO_f1)

function [xBest,fbest]=PSO(m,n,gen,a,b,gama,r1,r2,functie)
x=initpop(m,n,a,b); //population initialization
best=x; // best position of each element
v=zeros(m,n); // initialization of velocities
f=evalpop(functie,x);  // evaluation of the initial population
fb=f;
indBest=findBest(f);   // find the best element of the population
fbest=f(1,indBest); 
xBest=x(indBest,:);
for k=1:gen           
   for i=1:m           
      // update of velocity
      v(i,:)=gama*(v(i,:)+r1*grand(1,n,'unf',0,1).*(best(i,:)-x(i,:))+ r2*grand(1,n,'unf',0,1).*(xBest(1,:)-x(i,:)));
      // change the position
      x(i,:)=x(i,:)+v(i,:);
      // correction of elements which are outside the definition domain
      x(i,:)=repair(x(i,:),a,b);
      // evaluation of the new position
      f(1,i)=functie(x(i,:));
      // update the personal best
      if f(1,i)<fb(1,i)
          best(i,:)=x(i,:);
          fb(1,i)=f(1,i);
      end    
      // update the global best
      if f(1,i)<fbest
         indBest=i;
         fbest=f(1,i);
         xBest=x(indBest,:);
      end           
   end
end 
 // plot the best position of each particle
 
 if n==1
      xx=a:0.1:b; y=zeros(1,size(xx,2));
      for j=1:size(xx,2)
         y(j)=functie(xx(j));
      end   
      plot(xx,y,'b',best,fb,'r*');
 end     
endfunction


function pop = 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

function y=repair(x,a,b)
y=zeros(1,size(x,2));
for j=1:size(x,2)
    if (x(1,j)<a) or (x(1,j)>b)
        y(1,j)=grand(1,1,'unf',a,b);
    else
        y(1,j)=x(1,j);
    end
end    
endfunction


function f = evalpop(functie,x)
for i=1:size(x,1)
    f(1,i)=functie(x(i,:));
end
endfunction

function indBest=findBest(f)
indBest=1;
for i=2:size(f,2)
    if f(1,i)<f(1,indBest)
        indBest=i;
    end
end
endfunction

function y=f1(x)
n=size(x,2);
y=0;
for j=1:n
   y=y+(x(j)-10)^2;
end
endfunction

function y=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

// [xBest,fbest]=PSO(m,n,gen,a,b,gamma,r1,r2,functie)
[xbest,fbest]=PSO(10,1,100,-50,50,1,0.5,0.5,f2)








