// parameters
function y=J(x)
  y=sum(x.*x-cos(2*%pi*x))+length(x)
//  y=x^2;
endfunction
n=2;xmin=-5;xmax=5;Vmax=1;
Npop=6;w=0.7;c1=0.7;c2=0.7;
Ngen=20;
//
//initialisation
PSO=zeros(Npop,3*n+1);       // size Npop*(3n)
PSO(:,1:n)=xmin+(xmax-xmin)*rand(Npop,n);
PSO(:,(n+1:2*n))=-Vmax+2*Vmax*rand(Npop,n);
PSO(:,(2*n+1:3*n))=PSO(:,1:n);
PSO(:,3*n+1)=1E10;
pg=zeros(1,n+1);pg(n+1)=1E10;pgval=[];
for i=1:Ngen
  val=[];
  for j=1:Npop
    valj=J(PSO(j,1:n))
    val=[val;valj] // the evaluation of the swarm
    // computation of the new pi
    if (valj<PSO(j,3*n+1)) then
      PSO(j,2*n+1:3*n)=PSO(j,1:n);
      PSO(j,3*n+1)=valj;
    end 
  end
  // computation of pg
  [a,b]=gsort(val);ind=b($);
  best=[PSO(ind,1:n),val(ind)];
  if(val(ind)<pg(n+1)) then 
    pg=best;
  end
 // disp(pg);
  pgval=[pgval,pg(n+1)];
  // new velocitiy
PSO(:,(n+1:2*n))=w*PSO(:,(n+1:2*n))+c1*rand(Npop,n).*(PSO(:,2*n+1:3*n)-PSO(:,1:n))+c2*rand(Npop,n).*(ones(Npop,1)*pg(1:n)-PSO(:,1:n));
  // new value of x
PSO(:,1:n)=PSO(:,1:n)+PSO(:,(n+1:2*n)); 
end
clf()
disp('best elemet found for Rastrigin function (position+value):');disp(pg)
disp('number of evaluations:');disp(Ngen*Npop)
plot(1:Ngen,pgval)

