function [pop,pbest,gbest]=initializepso(num, bounds, evalFN,evalOps,options)
% function [pop]=initializepso(populationSize, variableBounds,evalFN,
%                           evalOps,options)
%    initializega creates a matrix of random numbers with 
%    a number of rows equal to the populationSize and a number
%    columns equal to the number of rows in bounds plus 1 for
%    the f(x) value which is found by applying the evalFN.
%    This is used by the ga to create the population if it
%    is not supplied.
%
% pop            - the initial, evaluated, random population 
% populatoinSize - the size of the population, i.e. the number to create
% variableBounds - a matrix which contains the bounds of each variable, i.e.
%                  [var1_high var1_low; var2_high var2_low; ....]
% evalFN         - the evaluation fn, usually the name of the .m file for 
%                  evaluation
% evalOps        - any options to be passed to the eval function defaults []
% options        - options to the initialize function, ie. 
%                  [type prec] where eps is the epsilon value 
%                  and the second option is 1 for float and 0 for binary, 
%                  prec is the precision of the variables defaults [1e-6 1]

%rand('seed',512341234)
D=size(bounds,1);
pop=zeros(num,2*D+1);
diff=bounds(:,2)-bounds(:,1);
for m=1:D
    pop(:,m)=bounds(m,1)+rand(num,1)*diff(m);
end
for m=D+1:2*D
    pop(:,m)=(rand(num,1)-0.5)*2*diff(m-D);
end
for m=1:num
    eval(['[sol,val]=' evalFN '(pop(m,1:2),1);']);
%    [sol,val]=fitness(pop(m,1:2),1);
    pop(m,2*D+1)=val;
end
pbest(:,1:D)=pop(:,1:D);
pbest(:,D+1)=pop(:,2*D+1);
[y,n]=max(pop(:,2*D+1));
gbest(1,1:D)=pop(n,1:D);
gbest(1,D+1)=pop(n,2*D+1);