Hi all. As you can see in my profile, I am writing code about Particle Swarm Optimization in SAS. A few days ago, Rick recommended that I vectorize the algorithm to make it more optimal and perform better, as well as being much more aesthetically pleasing and easier to understand.
Here is the original code 100% functional and without vectorizing:
/*creating a template so as to display our results*/
proc template;
/*the main structure of this template comes from the SAS documentation (https://documentation.sas.com/doc/en/pgmsascdc/v_021/imlug/imlug_tables_sect015.htm)*/
define table PSO_Results;
column Bestfun bx1 bx2 average gradx1 gradx2; /* names and order of columns */
define header topHeader; /* header at top of table */
text "Best results of the simulation";
end;
define header SpanHeader; /* define spanning header */
text "Best variables"; /* title of spanning header */
start = x1; /* span starts at second column */
end = x2; /* span ends at third column */
end;
end;
run;
/****************************************BASIC PSO****************************************/
proc iml;
/*****initialization*****/
/*function we want to minimize*/
start f(x);
y= x[1]**2 + x[2]**2 -100 ;
return(y);
finish f;
xmin=-100;
xmax=100;
/***parameters and variables***/
n=20; /*number of particules*/
p=2; /*number of dimensions*/
c1=2.05; /*acceleration factor*/
c2=2.05; /*acceleration factor*/
maxt=1000; /*maximum iteration*/
maxrun=100; /*trials of the algorithm*/
X=j(n,p+1,0); /*initialization of the particle's matrix*/
lb=j(1,p,xmin);
ub=j(1,p,xmax);
/*other useful matrixes*/
ffmin=j(maxt,maxrun,0);
fft=j(1,maxrun,0);
fff=j(1, maxrun,0);
rgbest=j(maxrun,p,0);
nb_t=j(2, maxrun,0);
do run=1 to maxrun;
nb_t[1,run]=run;
/***initialization of the pso***/
do i=1 to n;
do j=1 to p;
/*sample in a uniform distribution to initialize the particles position*/
call randgen(u,'uniform');
X[i,j]=lb[j]+u*(ub[j]-lb[j]);
end;
/*to fill the third column of X (which corresponds to the output of the function), we evaluate the fitting for each particule*/
X[i,p+1]=f(X[i,]);
end;
pbest=X; /*initialization of pbest so we can iterate then*/
V=0.1*X; /*initialization of the velocity (10%) */
/*we find the current minimum of f, given the current values of X (if we want to minimize the function)*/
titi=min(pbest[,p+1]);
/*then we iterate to assign the coordinates of this minimum to gbest*/
do i=1 to n;
if pbest[i,p+1]=titi then gbest=pbest[i,];
else e=0;
end;
* print X, pbest, gbest;
/**********************************/
t=1; /*resetting the iteration counter*/
tolerance=1; /*resetting the tolerance(stopping criteria)*/
do until(tolerance < (10**-(12))|t=maxt);
/*for all the particles in all dimensions : */
do i=1 to n;
do j=1 to p;
/*velocity updates*/
call randgen(u1,'uniform');
call randgen(u2,'uniform');
V[i,j]=V[i,j]+c1*u1*(pbest[i,j]-X[i,j])+c2*u2*(gbest[1,j]-X[i,j]);
/*position updates*/
X[i,j]= X[i,j]+V[i,j];
/*handling boundary violations*/
if X[i,j]<lb[j] then X[i,j]=lb[j];
else e=0;
if X[i,j]>ub[j] then X[i,j]=ub[j];
else e=0;
end;
end;
/*computing fitness for all the particules + updating pbest*/
do i=1 to n;
X[i,p+1]=f(X[i,]);
if X[i,p+1]<pbest[i,p+1] then pbest[i,]=X[i,];
else e=0;
end;
ffmin[t,run]=gbest[,p+1];
fft[run]=t;
/* updating gbest*/
titi=min(pbest[,p+1]);
do i=1 to n;
if pbest[i,p+1]=titi then gbest=pbest[i,];
else e=0;
end;
/*calculating tolerance*/
if t>100 then tolerance=abs(ffmin[(t-100),run]-gbest[,p+1]);
else e=0;
t=t+1;
end;
*print gbest;
nb_t[2,run]=t; /*storing the number of iterations it took the algorithm to converge*/
fff[run]=gbest[,p+1];
rgbest[run,]=gbest[,1:p];
end;
/*displaying our best results by putting them in a template*/
fff2=fff`;
bestfun=mean(fff2); /*taking the best run among all of the runs of the simulation*/
bx1=mean(rgbest[1]);
bx2=mean(rgbest[2]);
*print bestrun, bestfun, best_variables;
titi=nb_t`;
average=mean(titi[,2]);
best_variables= bx1||bx2;
opt={1 1};
call nlpfdd(fx,gradx, hfx, 'f',best_variables, opt);
toto=bestfun || bx1 || bx2 || average || gradx; /*we concatenate the main results*/
best_results= TableCreate({"Bestfun" "bx1" "bx2" "average" "gradx1" "gradx2"}, toto);
call TablePrint(best_results) template="PSO_Results";
quit;
Here all my progress made but I have several problems that I mention in the code and that I would like to tell you about here :
/*creating a template so as to display our results*/
proc template;
/*the main structure of this template comes from the SAS documentation (https://documentation.sas.com/doc/en/pgmsascdc/v_021/imlug/imlug_tables_sect015.htm)*/
define table PSO_Results;
column Bestfun bx1 bx2 average gradx1 gradx2; /* names and order of columns */
define header topHeader; /* header at top of table */
text "Best results of the simulation";
end;
define header SpanHeader; /* define spanning header */
text "Best variables"; /* title of spanning header */
start = x1; /* span starts at second column */
end = x2; /* span ends at third column */
end;
end;
run;
proc iml;
/* vectorization of f: R^2 --> R */
start f(x);
ros= 100*(x[,2]-x[,1]##2)##2+(x[,1]-1)##2;
return(ros);
finish f;
xmin=-5; /* min has to be >= 0 or else x^(1/4) and x^(3/4) are undefined */
xmax=10;
/***parameters and variables***/
n=20; /*number of particles*/
p=2;
maxt=1000; /*maximum iteration*/
maxrun=100; /*trials of the algorithm*/
X=j(n,p+1,0); /*initialization of the particle's matrix*/
lb=j(1,p,xmin);
ub=j(1,p,xmax);
c1=2.05; /*acceleration factor*/
c2=2.05; /*acceleration factor*/
phi=c1+c2;
chi= 2/abs(2-phi-sqrt(phi**2 -4*phi));
/*other useful matrixes*/
ffmin=j(maxt,maxrun,0);
fft=j(1,maxrun,0);
fff=j(1, maxrun,0);
rgbest=j(maxrun,p,0);
nb_t=j(2, maxrun,0);
call randseed(1234);
Y = j(n,P);/* allocate for random values */
do run=1 to maxrun; /*this line is not important in this example*/
/***initialization of the pso***/
call randgen(Y, 'uniform');
/*sample in a uniform distribution to initialize the particles position*/
X[,1:p] = lb + (ub-lb)#Y;
/*to fill the third column of X (which corresponds to the output of the function), we evaluate the fitting for each particule*/
X[,p+1]=f(X);
pbest=X; /*initialization of pbest so we can iterate then*/
V=0.1*X[,1:p];
/*then we iterate to assign the coordinates of this minimum to gbest*/
titi=min(pbest[,p+1]);
minIdx = loc(pbest[,p+1]=titi);
if ncol(minIdx)>0 then gbest = pbest[minIdx, ];
else e=0;
t=1; /*resetting the iteration counter*/
tolerance=1; /*resetting the tolerance(stopping criteria)*/
do until(tolerance < (10**-(12))|t=maxt);
/*for all the particles in all dimensions : */
call randseed(1234);
u1 = j(n,p);
u2 = j(n,p);
/*velocity updates*/
call randgen(u1,'uniform');
call randgen(u2,'uniform');
V=chi*(V+c1*u1#(pbest[,1:p]-X[,1:p]) + c2*u2#(gbest[,1:p]-X[,1:p]));
/*position updates*/
X[,1:p] = X[,1:p] + V;
/*handling boundary violations*/
X[,1:p]=X[,1:p]<>lb;
X[,1:p]=X[,1:p]><ub;
/*computing fitness for all the particules + updating pbest*/
X[,p+1]=f(X[,1:p]);
inferior_Idx=loc(X[,p+1]<pbest[,p+1]);
if ncol(inferior_Idx)>0 then pbest[inferior_Idx,]=X[inferior_Idx,];
else e=0;
ffmin[t,run]=gbest[,p+1];
fft[run]=t;
/* updating gbest*/
titi=min(pbest[,p+1]);
minIdx_2 = loc(pbest[,p+1]=titi);
if ncol(minIdx_2)>0 then
gbest = pbest[minIdx, ];
/*calculating tolerance*/
if t>100 then tolerance=abs(ffmin[(t-100),run]-gbest[,p+1]);
else e=0;
t=t+1;
end;
*print gbest;
nb_t[2,run]=t; /*storing the number of iterations it took the algorithm to converge*/
fff[run]=gbest[,p+1];
rgbest[run,]=gbest[,1:p];
end;
/*displaying our best results by putting them in a template*/
fff2=fff`;
bestfun=mean(fff2); /*taking the best run among all of the runs of the simulation*/
bx1=mean(rgbest[,1]);
bx2=mean(rgbest[,2]);
titi=nb_t`;
print titi;
average=mean(titi[,2]);
best_variables= bx1||bx2;
opt={1 1};
call nlpfdd(fx,gradx, hfx, 'f',best_variables, opt);
toto=bestfun || bx1 || bx2 || average || gradx; /*we concatenate the main results*/
best_results= TableCreate({"Bestfun" "bx1" "bx2" "average" "gradx1" "gradx2"}, toto);
call TablePrint(best_results) template="PSO_Results";
quit;
To begin with, the definition of bounds is complicated. In this case they work but only for 2 columns or p=2. I would like to generalize it but I don't know how.
Another problem lies in the fact of speed. Gbest is a 1x3 format vector. while X in this case is nx2. subtracting from each other is tricky. Actually what I want is for it to be subtracted line by line. I've thought about concatenating gbest n times so I can subtract it. But it seems to me that there has to be another way to do it.
At the moment everything is ready, although I still have many pending things. Like for example, how do you remove the do until loop in vector format?
Thank you so much everyone. I am learning more than ever with all your advice. And being a bit heavy.
All the best.
