function [x,endPop]=pso(bounds,w,c1,c2,termFns,termOps,evalFN,...
evalOps,startPop,pbest,gbest);
% PSO run a particle swarm algorithm
% [x,endPop,bestPop,trace]=pso(bounds,w,c1,c2,termFns,termOps,evalFn,evalOps,
%                               startPop,pbest,gbest);
%                                
% Output Arguments:
%   x            - the best solution found during the course of the run
%   endPop       - the final population 
%   bPop         - a trace of the best population
%   traceInfo    - a matrix of best and means of the ga for each generation
%
% Input Arguments:
%   bounds       - a matrix of upper and lower bounds on the variables
%   w            - the inertia weight
%   c1,c2        - the acceleration constants
%   termFns      - name of the .m termination function (['maxGenTerm'])
%   termOps      - options string to be passed to the termination function
%                  ([100]).
%   evalFn       - the name of the evaluation .m function
%   evalOps      - options to pass to the evaluation function ([NULL])
%   startPop     - a matrix of solutions that can be initialized
%                  from initialize.m
%   pbest        - the initial best solution of each particle
%   gbest        - the initial best solution of generation 1

%  This program is written by Li Jun on Dec,15,2003.All rights reserved

D    = size(bounds,1);
num  = size(startPop,1);
gen  = 1     ;
done = 0     ;
vnow = startPop(:,(D+1):2*D);
xnow = startPop(:,1:D);
xnew = xnow;
pb   = pbest;
gb   = gbest;
while(~done)
    %Caculate the fitness value of position xnow
    xnow=xnew;
    for m=1:num
    eval(['[sol,val]=' evalFN '(xnow(m,1:2),1);']);
    value(m,1)=val;
    end
    %replace the pbest
    for m=1:num
        if pb(m,D+1) < value(m,1)
            pb(m,1:D)=xnow(m,1:D);
            pb(m,D+1)=value(m,1);
        end
    end
    %replace the gbest
    [maxval,col]=max(value);
    if gb(D+1) < maxval
        gb=[xnow(col,:),maxval];
        bestg=gen;
    end
    %caculate the vnew and xnew
    vnow = w*vnow+c1*rand*(pb(:,1:D)-xnow)+c2*rand*([gb(1)*ones(num,1),...
           gb(2)*ones(num,1)]-xnow);
    xnew = xnow+vnow;
    if (100-gb(3)<1e-6) | (gen >= termOps)
        done=1;
    else
        gen=gen+1;
        w=w-0.01;    % w declines linearly
        fprintf(1,'%d ',gen);
    end
end
endPop=[xnow,value];
x=gb;
fprintf(1,'\n w= %f',w);
fprintf(1,'\n bestgen= %d',bestg);