|
|
|
1 |
#include "gadget.h" //All the required standard header files are in here |
#include "gadget.h" //All the required standard header files are in here |
2 |
#include "optinfo.h" |
#include "optinfo.h" |
3 |
#include "mathfunc.h" |
#include "mathfunc.h" |
7 |
#include "ecosystem.h" |
#include "ecosystem.h" |
8 |
#include "global.h" |
#include "global.h" |
9 |
#include "pso.h" |
#include "pso.h" |
10 |
|
#include <float.h> |
11 |
|
|
12 |
#ifdef _OPENMP |
#ifdef _OPENMP |
13 |
#include "omp.h" |
#include "omp.h" |
18 |
extern Ecosystem** EcoSystems; |
extern Ecosystem** EcoSystems; |
19 |
#endif |
#endif |
20 |
|
|
|
|
|
21 |
//============================================================== |
//============================================================== |
22 |
//// calulate swarm size based on dimensionality |
//// calulate swarm size based on dimensionality |
23 |
int OptInfoPso::pso_calc_swarm_size(int dim) { |
int OptInfoPso::pso_calc_swarm_size(int dim) { |
28 |
void OptInfoPso::OptimiseLikelihood() { |
void OptInfoPso::OptimiseLikelihood() { |
29 |
|
|
30 |
handle.logMessage(LOGINFO, "\nStarting PSO optimisation algorithm\n"); |
handle.logMessage(LOGINFO, "\nStarting PSO optimisation algorithm\n"); |
31 |
double oldf, newf, steplength, tmp; |
double bestf, tmp; |
32 |
int i, offset,rchange,rcheck,rnumber; |
int i, offset,rchange,rcheck,rnumber; |
33 |
int nvars = EcoSystem->numOptVariables(); |
int nvars = EcoSystem->numOptVariables(); |
34 |
DoubleVector x(nvars); |
DoubleVector x(nvars); |
35 |
DoubleVector trialx(nvars); |
// DoubleVector trialx(nvars); |
36 |
DoubleVector bestx(nvars); |
DoubleVector bestx(nvars); |
|
struct Best {double bestf; int index;}; |
|
|
struct Best best; |
|
|
double bestf; |
|
37 |
DoubleVector lowerb(nvars); |
DoubleVector lowerb(nvars); |
38 |
DoubleVector upperb(nvars); |
DoubleVector upperb(nvars); |
39 |
DoubleVector init(nvars); |
DoubleVector init(nvars); |
46 |
/** |
/** |
47 |
* Internal variables for PSO |
* Internal variables for PSO |
48 |
*/ |
*/ |
49 |
|
// printf("nvars:%d\n",nvars); |
50 |
|
if (size == 0) size = pso_calc_swarm_size(nvars); |
51 |
DoubleMatrix pos(size,nvars,0.0); // position matrix |
DoubleMatrix pos(size,nvars,0.0); // position matrix |
52 |
DoubleMatrix vel(size,nvars,0.0); // velocity matrix |
DoubleMatrix vel(size,nvars,0.0); // velocity matrix |
53 |
DoubleMatrix pos_b(size,nvars,0.0); // best position matrix |
DoubleMatrix pos_b(size,nvars,0.0); // best position matrix |
64 |
Calc_inertia_fun calc_inertia_fun; // inertia weight update function |
Calc_inertia_fun calc_inertia_fun; // inertia weight update function |
65 |
int psosteps=0; |
int psosteps=0; |
66 |
|
|
67 |
if (this->seed!=0) |
if (seed == 0) seed = 1234; |
|
srand(this->seed); |
|
68 |
|
|
69 |
handle.logMessage(LOGINFO,"Starting PSO with particles ",size,"\n"); |
handle.logMessage(LOGINFO,"Starting PSO with particles ",size,"\n"); |
70 |
switch (nhood_strategy) |
|
71 |
{ |
// SELECT APPROPRIATE NHOOD UPDATE FUNCTION |
72 |
|
switch (nhood_strategy) { |
73 |
case PSO_NHOOD_GLOBAL : |
case PSO_NHOOD_GLOBAL : |
74 |
// comm matrix not used |
// comm matrix not used |
75 |
handle.logMessage(LOGINFO,"PSO will use global communication\n"); |
handle.logMessage(LOGINFO,"PSO will use global communication\n"); |
87 |
break; |
break; |
88 |
} |
} |
89 |
// SELECT APPROPRIATE INERTIA WEIGHT UPDATE FUNCTION |
// SELECT APPROPRIATE INERTIA WEIGHT UPDATE FUNCTION |
90 |
switch (w_strategy) |
switch (w_strategy) { |
|
{ |
|
91 |
case PSO_W_CONST : |
case PSO_W_CONST : |
92 |
handle.logMessage(LOGINFO,"PSO will use constant inertia\n"); |
handle.logMessage(LOGINFO,"PSO will use constant inertia\n"); |
93 |
calc_inertia_fun = &OptInfoPso::calc_inertia_const; |
calc_inertia_fun = &OptInfoPso::calc_inertia_const; |
106 |
handle.logMessage(LOGINFO,"PSO c2",c2,"\n"); |
handle.logMessage(LOGINFO,"PSO c2",c2,"\n"); |
107 |
handle.logMessage(LOGINFO,"PSO goal",goal,"\n"); |
handle.logMessage(LOGINFO,"PSO goal",goal,"\n"); |
108 |
|
|
109 |
|
// if (scale) |
110 |
EcoSystem->scaleVariables(); |
EcoSystem->scaleVariables(); |
111 |
#ifdef _OPENMP |
//#ifdef _OPENMP |
112 |
int numThr = omp_get_max_threads ( ); |
// int numThr = omp_get_max_threads ( ); |
113 |
#pragma omp parallel for |
//#pragma omp parallel for |
114 |
for (i=0;i<numThr;i++) |
// for (i=0;i<numThr;i++) |
115 |
EcoSystems[i]->scaleVariables(); |
// EcoSystems[i]->scaleVariables(); |
116 |
#endif |
//#endif |
117 |
EcoSystem->getOptScaledValues(x); |
EcoSystem->getOptScaledValues(x); |
118 |
EcoSystem->getOptLowerBounds(lowerb); |
EcoSystem->getOptLowerBounds(lowerb); |
119 |
EcoSystem->getOptUpperBounds(upperb); |
EcoSystem->getOptUpperBounds(upperb); |
120 |
EcoSystem->getOptInitialValues(init); |
EcoSystem->getOptInitialValues(init); |
121 |
|
|
122 |
|
for (i = 0; i < nvars; ++i) { |
123 |
|
bestx[i] = x[i]; |
124 |
|
param[i] = i; |
125 |
|
} |
126 |
|
|
127 |
|
//se escala? |
128 |
|
// if (scale) |
129 |
for (i = 0; i < nvars; i++) { |
for (i = 0; i < nvars; i++) { |
130 |
// Scaling the bounds, because the parameters are scaled |
// Scaling the bounds, because the parameters are scaled |
131 |
lowerb[i] = lowerb[i] / init[i]; |
lowerb[i] = lowerb[i] / init[i]; |
135 |
lowerb[i] = upperb[i]; |
lowerb[i] = upperb[i]; |
136 |
upperb[i] = tmp; |
upperb[i] = tmp; |
137 |
} |
} |
|
bestx[i] = x[i]; |
|
|
trialx[i] = x[i]; |
|
|
param[i] = i; |
|
138 |
} |
} |
139 |
|
|
140 |
bestf = EcoSystem->SimulateAndUpdate(trialx); |
bestf = EcoSystem->SimulateAndUpdate(x); |
141 |
if (bestf != bestf) { //check for NaN |
if (bestf != bestf) { //check for NaN |
142 |
handle.logMessage(LOGINFO, "Error starting PSO optimisation with f(x) = infinity"); |
handle.logMessage(LOGINFO, "Error starting PSO optimisation with f(x) = infinity"); |
143 |
converge = -1; |
converge = -1; |
146 |
} |
} |
147 |
|
|
148 |
offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop |
offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop |
149 |
|
// newf = bestf; |
150 |
|
// oldf = bestf; |
151 |
|
|
152 |
iters = 0; |
iters = 0; |
153 |
|
|
156 |
// |
// |
157 |
i=0; |
i=0; |
158 |
handle.logMessage(LOGINFO,"Initialising PSO ",size,"particles\n"); |
handle.logMessage(LOGINFO,"Initialising PSO ",size,"particles\n"); |
159 |
|
|
160 |
|
// SWARM INITIALIZATION |
161 |
// |
// |
162 |
// Initialize the fist particle with the best known position |
// Initialize the fist particle with the best known position |
163 |
// |
// |
164 |
for (d=0; d<nvars; d++) { |
for (d=0; d<nvars; d++) { |
165 |
// generate two numbers within the specified range and around trialx |
// generate two numbers within the specified range and around trialx |
166 |
a = trialx[d] ; |
a = x[d]; |
167 |
// initialize position |
// initialize position |
168 |
pos[i][d] = a; |
pos[i][d] = a; |
169 |
pos_b[i][d] = a; |
pos_b[i][d] = a; |
170 |
pos_nb[i][d] = a; |
pos_nb[i][d] = a; |
171 |
a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX); |
// a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand() * 1.0 / RAND_MAX); |
172 |
b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX); |
b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX); |
173 |
vel[i][d] = (a-b) / 2.; |
vel[i][d] = (a-b) / 2.; |
174 |
} |
} |
175 |
fit[i] = EcoSystem->SimulateAndUpdate(pos[i]); |
fit[i] = EcoSystem->SimulateAndUpdate(pos[i]); |
176 |
fit_b[i] = fit[i]; // this is also the personal best |
fit_b[i] = fit[i]; // this is also the personal best |
177 |
if (fit[i] < bestf) { |
if (fit[i] < bestf) { |
178 |
bestf = fit[i]; |
bestf = fit[i]; |
|
best.bestf=bestf; |
|
179 |
bestx=pos[i]; |
bestx=pos[i]; |
|
best.index=i; |
|
180 |
} |
} |
181 |
// |
|
182 |
// Now other particles with a random position in the solution space |
// Now other particles with a random position in the solution space |
183 |
// |
// |
184 |
best.bestf=bestf; |
//#pragma omp parallel for default(shared) private(i,d,a,b) |
|
best.index=0; |
|
|
handle.logMessage(LOGINFO,"before Best.bestf=",best.bestf); |
|
|
#pragma omp declare reduction(bestomp: struct Best : omp_out = omp_in.bestf < omp_out.bestf ? omp_in : omp_out) |
|
|
#pragma omp parallel for default(shared) private(i,d,a,b) reduction(bestomp:best) firstprivate(bestf) |
|
185 |
for (i=1; i<size; i++) { |
for (i=1; i<size; i++) { |
|
// for each dimension |
|
186 |
ostringstream spos,svel; |
ostringstream spos,svel; |
187 |
spos << "Initial position of particle " << i<<" [" ; |
spos << "Initial position of particle " << i<<" [" ; |
188 |
svel << "Initial velocity of particle " << i<<" [" ; |
svel << "Initial velocity of particle " << i<<" [" ; |
189 |
|
// for each dimension |
190 |
for (d=0; d<nvars; d++) { |
for (d=0; d<nvars; d++) { |
191 |
// generate two numbers within the specified range and around trialx |
// generate two numbers within the specified range and around trialx |
192 |
a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX); |
a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX); |
193 |
b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX); |
b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX); |
194 |
// initialize position |
// initialize position |
195 |
pos[i][d] = a; |
pos[i][d] = a; |
196 |
// best position is trialx |
// best position is trialx |
201 |
spos << pos[i][d] << " "; |
spos << pos[i][d] << " "; |
202 |
svel << vel[i][d] << " "; |
svel << vel[i][d] << " "; |
203 |
} |
} |
204 |
#ifdef _OPENMP |
//#ifdef _OPENMP |
205 |
fit[i] = EcoSystems[omp_get_thread_num ()]->SimulateAndUpdate(pos[i]); |
// fit[i] = EcoSystems[omp_get_thread_num ()]->SimulateAndUpdate(pos[i]); |
206 |
best.bestf = fit[i]; |
//#else |
|
#else |
|
207 |
fit[i] = EcoSystem->SimulateAndUpdate(pos[i]); |
fit[i] = EcoSystem->SimulateAndUpdate(pos[i]); |
208 |
#endif |
//#endif |
209 |
fit_b[i] = fit[i]; // this is also the personal best |
fit_b[i] = fit[i]; // this is also the personal best |
210 |
// update gbest?? |
// update gbest?? |
211 |
//#pragma omp critical (best1) |
//#pragma omp critical (best1) |
212 |
if (fit[i] < bestf) { |
if (fit[i] < bestf) { |
213 |
// update best fitness |
// update best fitness |
214 |
handle.logMessage(LOGINFO,"iInside Best.bestf=",best.bestf); |
bestf = fit[i]; |
215 |
// copy particle pos to gbest vector |
// copy particle pos to gbest vector |
216 |
//memmove((void *)&bestx, (void *)&pos[i], |
//memmove((void *)&bestx, (void *)&pos[i], |
217 |
// sizeof(double) * nvars); |
// sizeof(double) * nvars); |
218 |
best.index=i; |
bestx = pos[i]; |
|
bestf=best.bestf; |
|
|
//bestx=pos[i]; |
|
219 |
} |
} |
220 |
spos<<" ]"; |
spos<<" ]"; |
221 |
svel<<" ]"; |
svel<<" ]"; |
222 |
handle.logMessage(LOGDEBUG, spos.str().c_str()); |
handle.logMessage(LOGDEBUG, spos.str().c_str()); |
223 |
handle.logMessage(LOGDEBUG, svel.str().c_str()); |
handle.logMessage(LOGDEBUG, svel.str().c_str()); |
224 |
|
|
|
|
|
|
} |
|
|
// |
|
|
// Compare the reduction with the best value |
|
|
// |
|
|
|
|
|
if (best.bestf < bestf){ |
|
|
bestx=pos[best.index]; |
|
|
bestf=best.bestf; |
|
225 |
} |
} |
226 |
|
|
|
|
|
227 |
// |
// |
228 |
// RUN ALGORITHM |
// RUN ALGORITHM |
229 |
// |
// |
246 |
handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of PSO steps was reached"); |
handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of PSO steps was reached"); |
247 |
handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run"); |
handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run"); |
248 |
|
|
249 |
bestf = EcoSystem->SimulateAndUpdate(bestx); |
score = EcoSystem->SimulateAndUpdate(bestx); |
250 |
for (i = 0; i < nvars; i++) |
for (i = 0; i < nvars; i++) |
251 |
bestx[i] = bestx[i] * init[i]; |
bestx[i] = bestx[i] * init[i]; |
252 |
handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", bestf); |
handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", score); |
253 |
EcoSystem->storeVariables(bestf, bestx); |
EcoSystem->storeVariables(score, bestx); |
254 |
return; |
return; |
255 |
} |
} |
256 |
|
// update inertia weight |
257 |
if (w_strategy) |
if (w_strategy) w = (this->*calc_inertia_fun)(psosteps); |
|
w = (this->*calc_inertia_fun)(psosteps); |
|
258 |
// check optimization goal |
// check optimization goal |
259 |
if (bestf <= goal) { |
if (bestf <= goal) { |
260 |
handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n"); |
handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n"); |
261 |
handle.logMessage(LOGINFO, "Goal achieved!!! @ step ", psosteps); |
handle.logMessage(LOGINFO, "Goal achieved!!! @ step ", psosteps); |
262 |
bestf = EcoSystem->SimulateAndUpdate(bestx); |
score = EcoSystem->SimulateAndUpdate(bestx); |
263 |
handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", bestf); |
handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", score); |
264 |
break; |
break; |
265 |
} |
} |
266 |
|
|
271 |
improved = 0; |
improved = 0; |
272 |
|
|
273 |
// update all particles |
// update all particles |
274 |
best.bestf=bestf; |
//#pragma omp parallel for private(d,rho1,rho2) default(shared) |
|
#pragma omp parallel for default(shared) private(i,d,rho1,rho2) reduction(bestomp:best) firstprivate(bestf) |
|
275 |
for (i=0; i<size; i++) { |
for (i=0; i<size; i++) { |
276 |
// for each dimension |
// handle.logMessage(LOGDEBUG, "In parallel",omp_get_thread_num(),"\n"); |
|
handle.logMessage(LOGDEBUG, "In parallel",omp_get_thread_num(),"\n"); |
|
277 |
ostringstream spos,svel; |
ostringstream spos,svel; |
278 |
spos << "New position of particle " << i<<" [" ; |
spos << "New position of particle " << i<<" [" ; |
279 |
svel << "New velocity of particle " << i<<" [" ; |
svel << "New velocity of particle " << i<<" [" ; |
280 |
|
// for each dimension |
281 |
for (d=0; d<nvars; d++) { |
for (d=0; d<nvars; d++) { |
282 |
// calculate stochastic coefficients |
// calculate stochastic coefficients |
283 |
rho1 = c1 * (rand()*1.0)/RAND_MAX; |
rho1 = c1 * (rand_r(&seed) * 1.0) / RAND_MAX; |
284 |
rho2 = c2 * (rand()*1.0)/RAND_MAX; |
rho2 = c2 * (rand_r(&seed) * 1.0) / RAND_MAX; |
285 |
// update velocity |
// update velocity |
286 |
vel[i][d] = w * vel[i][d] + \ |
vel[i][d] = w * vel[i][d] + rho1 * (pos_b[i][d] - pos[i][d]) + rho2 * (pos_nb[i][d] - pos[i][d]); |
|
rho1 * (pos_b[i][d] - pos[i][d]) + \ |
|
|
rho2 * (pos_nb[i][d] - pos[i][d]); |
|
287 |
// update position |
// update position |
288 |
pos[i][d] += vel[i][d]; |
pos[i][d] += vel[i][d]; |
289 |
spos << pos[i][d] <<" "; |
spos << pos[i][d] <<" "; |
293 |
if (pos[i][d] < lowerb[d]) { |
if (pos[i][d] < lowerb[d]) { |
294 |
pos[i][d] = lowerb[d]; |
pos[i][d] = lowerb[d]; |
295 |
vel[i][d] = 0; |
vel[i][d] = 0; |
296 |
} else if (pos[i][d] > upperb[d]) { |
} |
297 |
|
else if (pos[i][d] > upperb[d]) { |
298 |
pos[i][d] = upperb[d]; |
pos[i][d] = upperb[d]; |
299 |
vel[i][d] = 0; |
vel[i][d] = 0; |
300 |
} |
} |
301 |
} else { |
} |
302 |
|
else { |
303 |
// enforce periodic boundary conditions |
// enforce periodic boundary conditions |
304 |
if (pos[i][d] < lowerb[d]) { |
if (pos[i][d] < lowerb[d]) { |
305 |
|
pos[i][d] = upperb[d] - fmod(lowerb[d] - pos[i][d], upperb[d] - lowerb[d]); |
|
pos[i][d] = upperb[d] - fmod(lowerb[d] - pos[i][d], |
|
|
upperb[d] - lowerb[d]); |
|
306 |
vel[i][d] = 0; |
vel[i][d] = 0; |
307 |
|
} |
308 |
} else if (pos[i][d] > upperb[d]) { |
else if (pos[i][d] > upperb[d]) { |
309 |
|
pos[i][d] = lowerb[d] + fmod(pos[i][d] - upperb[d], upperb[d] - lowerb[d]); |
|
pos[i][d] = lowerb[d] + fmod(pos[i][d] - upperb[d], |
|
|
upperb[d] - lowerb[d]); |
|
310 |
vel[i][d] = 0; |
vel[i][d] = 0; |
311 |
|
|
|
|
|
312 |
} |
} |
313 |
} |
} |
314 |
|
} //end for nvars |
315 |
|
|
|
} //end for d |
|
316 |
spos<< "]" << endl; |
spos<< "]" << endl; |
317 |
svel<< "]" << endl; |
svel<< "]" << endl; |
318 |
handle.logMessage(LOGDEBUG,spos.str().c_str()); |
handle.logMessage(LOGDEBUG,spos.str().c_str()); |
319 |
handle.logMessage(LOGDEBUG,svel.str().c_str()); |
handle.logMessage(LOGDEBUG,svel.str().c_str()); |
320 |
|
|
321 |
// update particle fitness |
// update particle fitness |
322 |
fit[i] = EcoSystems[omp_get_thread_num()]->SimulateAndUpdate(pos[i]); |
// fit[i] = EcoSystems[omp_get_thread_num()]->SimulateAndUpdate(pos[i]); |
323 |
|
fit[i] = EcoSystem->SimulateAndUpdate(pos[i]); |
324 |
// update personal best position? |
// update personal best position? |
|
best.bestf = fit[i]; |
|
325 |
if (fit[i] < fit_b[i]) { |
if (fit[i] < fit_b[i]) { |
326 |
fit_b[i] = fit[i]; |
fit_b[i] = fit[i]; |
|
// copy contents of pos[i] to pos_b[i] |
|
|
//memmove((void *)&pos_b[i], (void *)&pos[i], |
|
|
// sizeof(double) * nvars); |
|
327 |
pos_b[i]=pos[i]; |
pos_b[i]=pos[i]; |
328 |
} |
} |
329 |
// update gbest?? |
// update gbest?? |
334 |
if (fit[i] < bestf) { |
if (fit[i] < bestf) { |
335 |
improved = 1; |
improved = 1; |
336 |
// update best fitness |
// update best fitness |
|
best.index=i; |
|
337 |
bestf=fit[i]; |
bestf=fit[i]; |
|
// copy particle pos to gbest vector |
|
|
//memmove((void *)&bestx, (void *)&pos[i], |
|
|
// sizeof(double) * nvars); |
|
|
// for (d = 0; d < nvars; d++) |
|
|
// bestx[d] = pos[i][d] * init[d]; |
|
|
// iters = EcoSystem->getFuncEval() - offset; |
|
|
// handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations"); |
|
|
// handle.logMessage(LOGINFO, "The likelihood score is", best.bestf, "at the point"); |
|
|
// EcoSystem->storeVariables(best.bestf, bestx); |
|
|
// EcoSystem->writeBestValues(); |
|
|
// bestx=pos[i]; |
|
|
|
|
|
} |
|
|
} // end i |
|
|
if (best.bestf < bestf){ |
|
|
bestf=best.bestf; |
|
|
bestx=pos[best.index]; |
|
338 |
for (d = 0; d < nvars; d++) |
for (d = 0; d < nvars; d++) |
339 |
bestx[d] = bestx[d] * init[d]; |
bestx[d] = pos[i][d] * init[d]; |
340 |
iters = EcoSystem->getFuncEval() - offset; //To change. It is incorrect when is parallel. |
iters = EcoSystem->getFuncEval() - offset; |
341 |
handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations"); |
handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations"); |
342 |
handle.logMessage(LOGINFO, "The likelihood score is", bestf, "at the point"); |
handle.logMessage(LOGINFO, "The likelihood score is", bestf, "at the point"); |
343 |
EcoSystem->storeVariables(bestf, bestx); |
EcoSystem->storeVariables(bestf, bestx); |
344 |
EcoSystem->writeBestValues(); |
EcoSystem->writeBestValues(); |
345 |
bestx=pos[best.index]; |
bestx = pos[i]; |
346 |
|
|
347 |
} |
} |
348 |
|
} // end i |
349 |
|
|
350 |
psosteps++; |
psosteps++; |
351 |
|
|
358 |
EcoSystem->storeVariables(bestf, bestx); |
EcoSystem->storeVariables(bestf, bestx); |
359 |
EcoSystem->writeBestValues(); |
EcoSystem->writeBestValues(); |
360 |
|
|
|
|
|
361 |
} |
} |
362 |
|
|
363 |
|
#ifdef _OPENMP |
364 |
void OptInfoPso::OptimiseLikelihoodREP() { |
void OptInfoPso::OptimiseLikelihoodREP() { |
365 |
OptimiseLikelihood(); |
OptimiseLikelihood(); |
366 |
|
|
367 |
} |
} |
368 |
|
|
369 |
void OptInfoPso::OptimiseLikelihoodOMP() { |
void OptInfoPso::OptimiseLikelihoodOMP() { |
370 |
OptimiseLikelihood(); |
|
371 |
|
handle.logMessage(LOGINFO, "\nStarting PSO optimisation algorithm\n"); |
372 |
|
double tmp; |
373 |
|
int i, offset, rchange, rcheck, rnumber; |
374 |
|
int nvars = EcoSystem->numOptVariables(); |
375 |
|
DoubleVector x(nvars); |
376 |
|
// DoubleVector trialx(nvars); |
377 |
|
DoubleVector bestx(nvars); |
378 |
|
DoubleVector scalex(nvars); |
379 |
|
DoubleVector lowerb(nvars); |
380 |
|
DoubleVector upperb(nvars); |
381 |
|
DoubleVector init(nvars); |
382 |
|
//DoubleVector initialstep(nvars, rho); |
383 |
|
IntVector param(nvars, 0); |
384 |
|
IntVector lbound(nvars, 0); |
385 |
|
IntVector rbounds(nvars, 0); |
386 |
|
IntVector trapped(nvars, 0); |
387 |
|
struct Best { double bestf; int index;}; |
388 |
|
struct Best best; |
389 |
|
|
390 |
|
/** |
391 |
|
* Internal variables for PSO |
392 |
|
*/ |
393 |
|
// printf("nvars:%d\n",nvars); |
394 |
|
if (size == 0) size = pso_calc_swarm_size(nvars); |
395 |
|
DoubleMatrix pos(size, nvars, 0.0); // position matrix |
396 |
|
DoubleMatrix vel(size, nvars, 0.0); // velocity matrix |
397 |
|
DoubleMatrix pos_b(size, nvars, 0.0); // best position matrix |
398 |
|
DoubleVector fit(size); // particle fitness vector |
399 |
|
DoubleVector fit_b(size); // best fitness vector |
400 |
|
DoubleMatrix pos_nb(size, nvars, 0.0); // what is the best informed position for each particle |
401 |
|
IntMatrix comm(size, size, 0); // communications:who informs who rows : those who inform cols : those who are informed |
402 |
|
int improved; // whether solution->error was improved during the last iteration |
403 |
|
int d, step; |
404 |
|
double a, b; |
405 |
|
double rho1, rho2; //random numbers (coefficients) |
406 |
|
double w; //current omega |
407 |
|
Inform_fun inform_fun; |
408 |
|
Calc_inertia_fun calc_inertia_fun; // inertia weight update function |
409 |
|
int psosteps = 0; //!= iters? |
410 |
|
|
411 |
|
// if (seed == 0) seed = 1234; |
412 |
|
|
413 |
|
handle.logMessage(LOGINFO, "Starting PSO with particles ", size, "\n"); |
414 |
|
|
415 |
|
// SELECT APPROPRIATE NHOOD UPDATE FUNCTION |
416 |
|
switch (nhood_strategy) { |
417 |
|
case PSO_NHOOD_GLOBAL: |
418 |
|
// comm matrix not used |
419 |
|
handle.logMessage(LOGINFO, "PSO will use global communication\n"); |
420 |
|
inform_fun = &OptInfoPso::inform_global; |
421 |
|
break; |
422 |
|
case PSO_NHOOD_RING: |
423 |
|
handle.logMessage(LOGINFO, "PSO will use ring communication\n"); |
424 |
|
init_comm_ring(comm); |
425 |
|
inform_fun = &OptInfoPso::inform_ring; |
426 |
|
break; |
427 |
|
case PSO_NHOOD_RANDOM: |
428 |
|
handle.logMessage(LOGINFO, "PSO will use random communication\n"); |
429 |
|
init_comm_random(comm); |
430 |
|
inform_fun = &OptInfoPso::inform_random; |
431 |
|
break; |
432 |
|
} |
433 |
|
// SELECT APPROPRIATE INERTIA WEIGHT UPDATE FUNCTION |
434 |
|
switch (w_strategy) { |
435 |
|
case PSO_W_CONST: |
436 |
|
handle.logMessage(LOGINFO, "PSO will use constant inertia\n"); |
437 |
|
calc_inertia_fun = &OptInfoPso::calc_inertia_const; |
438 |
|
break; |
439 |
|
case PSO_W_LIN_DEC: |
440 |
|
handle.logMessage(LOGINFO, "PSO will use linear decreasing inertia\n"); |
441 |
|
calc_inertia_fun = &OptInfoPso::calc_inertia_lin_dec; |
442 |
|
break; |
443 |
|
} |
444 |
|
// |
445 |
|
//initialize omega using standard value |
446 |
|
// |
447 |
|
w = PSO_INERTIA; |
448 |
|
handle.logMessage(LOGINFO, "PSO initial inertia ", w, "\n"); |
449 |
|
handle.logMessage(LOGINFO, "PSO c1", c1, "\n"); |
450 |
|
handle.logMessage(LOGINFO, "PSO c2", c2, "\n"); |
451 |
|
handle.logMessage(LOGINFO, "PSO goal", goal, "\n"); |
452 |
|
|
453 |
|
EcoSystem->scaleVariables(); |
454 |
|
int numThr = omp_get_max_threads ( ); |
455 |
|
// if (scale) { |
456 |
|
#pragma omp parallel for |
457 |
|
for (i=0;i<numThr;i++) |
458 |
|
EcoSystems[i]->scaleVariables(); |
459 |
|
// } |
460 |
|
//#endif |
461 |
|
EcoSystem->getOptScaledValues(x); |
462 |
|
EcoSystem->getOptLowerBounds(lowerb); |
463 |
|
EcoSystem->getOptUpperBounds(upperb); |
464 |
|
EcoSystem->getOptInitialValues(init); |
465 |
|
|
466 |
|
|
467 |
|
for (i = 0; i < nvars; ++i) { |
468 |
|
bestx[i] = x[i]; |
469 |
|
param[i] = i; |
470 |
|
} |
471 |
|
//se escala? |
472 |
|
// if (scale) |
473 |
|
for (i = 0; i < nvars; i++) { |
474 |
|
scalex[i] = x[i]; |
475 |
|
// Scaling the bounds, because the parameters are scaled |
476 |
|
lowerb[i] = lowerb[i] / init[i]; |
477 |
|
upperb[i] = upperb[i] / init[i]; |
478 |
|
if (lowerb[i] > upperb[i]) { |
479 |
|
tmp = lowerb[i]; |
480 |
|
lowerb[i] = upperb[i]; |
481 |
|
upperb[i] = tmp; |
482 |
|
} |
483 |
|
} |
484 |
|
|
485 |
|
best.bestf = EcoSystem->SimulateAndUpdate(x); |
486 |
|
best.index=0; |
487 |
|
if (best.bestf != best.bestf) { //check for NaN |
488 |
|
handle.logMessage(LOGINFO, "Error starting PSO optimisation with f(x) = infinity"); |
489 |
|
converge = -1; |
490 |
|
iters = 1; |
491 |
|
return; |
492 |
|
} |
493 |
|
|
494 |
|
offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop |
495 |
|
// newf = bestf; |
496 |
|
// oldf = bestf; |
497 |
|
|
498 |
|
iters = 0; |
499 |
|
|
500 |
|
// |
501 |
|
// Initialize the particles |
502 |
|
// |
503 |
|
i = 0; |
504 |
|
handle.logMessage(LOGINFO, "Initialising PSO ", size, "particles\n"); |
505 |
|
|
506 |
|
// SWARM INITIALIZATION |
507 |
|
// |
508 |
|
// Initialize the fist particle with the best known position |
509 |
|
// |
510 |
|
// best.bestf=bestf; |
511 |
|
|
512 |
|
for (d = 0; d < nvars; d++) { |
513 |
|
// generate two numbers within the specified range and around trialx |
514 |
|
a = x[d]; |
515 |
|
// initialize position |
516 |
|
pos[i][d] = a; |
517 |
|
pos_b[i][d] = a; |
518 |
|
pos_nb[i][d] = a; |
519 |
|
// a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand() * 1.0 / RAND_MAX); |
520 |
|
b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX); |
521 |
|
vel[i][d] = (a - b) / 2.; |
522 |
|
} |
523 |
|
fit[i] = EcoSystem->SimulateAndUpdate(pos[i]); |
524 |
|
fit_b[i] = fit[i]; // this is also the personal best |
525 |
|
if (fit[i] < best.bestf) { |
526 |
|
best.bestf = fit[i]; |
527 |
|
// best.bestf=bestf; |
528 |
|
bestx = pos[i]; |
529 |
|
best.index=i; |
530 |
} |
} |
531 |
|
|
532 |
|
// Now other particles with a random position in the solution space |
533 |
|
// |
534 |
|
// handle.logMessage(LOGINFO, "ini 1 best.bestf=",best.bestf, "\n"); |
535 |
|
|
536 |
|
#pragma omp declare reduction(bestomp: Best : omp_out = omp_in.bestf < omp_out.bestf ? omp_in : omp_out) initializer(omp_priv=omp_orig) |
537 |
|
#pragma omp parallel for default(shared) private(i,d,a,b) reduction(bestomp:best) //firstprivate(bestf) |
538 |
|
for (i = 1; i < size; i++) { |
539 |
|
ostringstream spos, svel; |
540 |
|
spos << "Initial position of particle " << i << " ["; |
541 |
|
svel << "Initial velocity of particle " << i << " ["; |
542 |
|
// for each dimension |
543 |
|
for (d = 0; d < nvars; d++) { |
544 |
|
// generate two numbers within the specified range and around trialx |
545 |
|
a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX); |
546 |
|
b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX); |
547 |
|
// initialize position |
548 |
|
pos[i][d] = a; |
549 |
|
// best position is trialx |
550 |
|
pos_b[i][d] = a; |
551 |
|
pos_nb[i][d] = a; |
552 |
|
// initialize velocity |
553 |
|
vel[i][d] = (a - b) / 2.; |
554 |
|
spos << pos[i][d] << " "; |
555 |
|
svel << vel[i][d] << " "; |
556 |
|
} |
557 |
|
//#ifdef _OPENMP |
558 |
|
fit[i] = EcoSystems[omp_get_thread_num ()]->SimulateAndUpdate(pos[i]); |
559 |
|
// best.bestf = fit[i]; |
560 |
|
//#else |
561 |
|
// fit[i] = EcoSystem->SimulateAndUpdate(pos[i]); |
562 |
|
//#endif |
563 |
|
fit_b[i] = fit[i]; // this is also the personal best |
564 |
|
// update gbest?? |
565 |
|
//#pragma omp critical (best1) |
566 |
|
if (fit[i] < best.bestf) { |
567 |
|
// update best fitness |
568 |
|
// bestf = fit[i]; |
569 |
|
// copy particle pos to gbest vector |
570 |
|
//memmove((void *)&bestx, (void *)&pos[i], |
571 |
|
// sizeof(double) * nvars); |
572 |
|
// bestx = pos[i]; |
573 |
|
best.index=i; |
574 |
|
// bestf=best.bestf; |
575 |
|
best.bestf = fit[i]; |
576 |
|
} |
577 |
|
spos << " ]"; |
578 |
|
svel << " ]"; |
579 |
|
handle.logMessage(LOGDEBUG, spos.str().c_str()); |
580 |
|
handle.logMessage(LOGDEBUG, svel.str().c_str()); |
581 |
|
|
582 |
|
} |
583 |
|
|
584 |
|
// |
585 |
|
// Compare the reduction with the best value |
586 |
|
// |
587 |
|
|
588 |
|
// if (best.bestf < bestf){ |
589 |
|
// bestx=pos[best.index]; |
590 |
|
// bestf=best.bestf; |
591 |
|
// } |
592 |
|
|
593 |
|
// |
594 |
|
// RUN ALGORITHM |
595 |
|
// |
596 |
|
// handle.logMessage(LOGINFO, "ini 2 best.bestf=",best.bestf, "\n"); |
597 |
|
psosteps = 0; |
598 |
|
while (1) { |
599 |
|
|
600 |
|
handle.logMessage(LOGINFO, "PSO optimisation after", psosteps * size, "\n"); |
601 |
|
if (isZero(best.bestf)) { |
602 |
|
iters = (EcoSystems[0]->getFuncEval() * numThr) + EcoSystem->getFuncEval() - offset; |
603 |
|
handle.logMessage(LOGINFO, "Error in PSO optimisation after", iters, "function evaluations, f(x) = 0"); |
604 |
|
converge = -1; |
605 |
|
return; |
606 |
|
} |
607 |
|
/* if too many function evaluations occur, terminate the algorithm */ |
608 |
|
|
609 |
|
iters = (EcoSystems[0]->getFuncEval() * numThr) + EcoSystem->getFuncEval() - offset; |
610 |
|
if (psosteps > psoiter) { |
611 |
|
handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n"); |
612 |
|
handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations"); |
613 |
|
handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of PSO steps was reached"); |
614 |
|
handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run"); |
615 |
|
|
616 |
|
score = EcoSystem->SimulateAndUpdate(bestx); |
617 |
|
for (i = 0; i < nvars; i++) |
618 |
|
bestx[i] = bestx[i] * init[i]; |
619 |
|
handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", score); |
620 |
|
EcoSystem->storeVariables(score, bestx); |
621 |
|
return; |
622 |
|
} |
623 |
|
// update inertia weight |
624 |
|
if (w_strategy) w = (this->*calc_inertia_fun)(psosteps); |
625 |
|
// check optimization goal |
626 |
|
if (best.bestf <= goal) { |
627 |
|
handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n"); |
628 |
|
handle.logMessage(LOGINFO, "Goal achieved!!! @ step ", psosteps); |
629 |
|
score = EcoSystem->SimulateAndUpdate(bestx); |
630 |
|
handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", score); |
631 |
|
break; |
632 |
|
} |
633 |
|
|
634 |
|
// update pos_nb matrix (find best of neighborhood for all particles) |
635 |
|
(this->*inform_fun)(comm, pos_nb, pos_b, fit_b, bestx, improved); |
636 |
|
|
637 |
|
// the value of improved was just used; reset it |
638 |
|
improved = 0; |
639 |
|
|
640 |
|
// update all particles |
641 |
|
// best.bestf=bestf; |
642 |
|
//#pragma omp declare reduction(bestomp2: Best : omp_out = omp_in.bestf < omp_out.bestf ? omp_in : omp_out) initializer(omp_priv=omp_orig) |
643 |
|
#pragma omp parallel for default(shared) private(i,d,rho1,rho2) reduction(bestomp:best) //firstprivate(best) |
644 |
|
for (i = 0; i < size; i++) { |
645 |
|
// handle.logMessage(LOGDEBUG, "In parallel",omp_get_thread_num(),"\n"); |
646 |
|
ostringstream spos, svel; |
647 |
|
spos << "New position of particle " << i << " ["; |
648 |
|
svel << "New velocity of particle " << i << " ["; |
649 |
|
// for each dimension |
650 |
|
for (d = 0; d < nvars; d++) { |
651 |
|
// calculate stochastic coefficients |
652 |
|
rho1 = c1 * (rand_r(&seed) * 1.0) / RAND_MAX; |
653 |
|
rho2 = c2 * (rand_r(&seed) * 1.0) / RAND_MAX; |
654 |
|
// update velocity |
655 |
|
vel[i][d] = w * vel[i][d] + rho1 * (pos_b[i][d] - pos[i][d]) + rho2 * (pos_nb[i][d] - pos[i][d]); |
656 |
|
// update position |
657 |
|
pos[i][d] += vel[i][d]; |
658 |
|
spos << pos[i][d] << " "; |
659 |
|
svel << vel[i][d] << " "; |
660 |
|
// clamp position within bounds? |
661 |
|
if (clamp_pos) { |
662 |
|
if (pos[i][d] < lowerb[d]) { |
663 |
|
pos[i][d] = lowerb[d]; |
664 |
|
vel[i][d] = 0; |
665 |
|
} |
666 |
|
else if (pos[i][d] > upperb[d]) { |
667 |
|
pos[i][d] = upperb[d]; |
668 |
|
vel[i][d] = 0; |
669 |
|
} |
670 |
|
} |
671 |
|
else { |
672 |
|
// enforce periodic boundary conditions |
673 |
|
if (pos[i][d] < lowerb[d]) { |
674 |
|
pos[i][d] = upperb[d] - fmod(lowerb[d] - pos[i][d], upperb[d] - lowerb[d]); |
675 |
|
vel[i][d] = 0; |
676 |
|
} |
677 |
|
else if (pos[i][d] > upperb[d]) { |
678 |
|
pos[i][d] = lowerb[d] + fmod(pos[i][d] - upperb[d], upperb[d] - lowerb[d]); |
679 |
|
vel[i][d] = 0; |
680 |
|
|
681 |
|
} |
682 |
|
} |
683 |
|
} //end for nvars |
684 |
|
|
685 |
|
spos << "]" << endl; |
686 |
|
svel << "]" << endl; |
687 |
|
handle.logMessage(LOGDEBUG, spos.str().c_str()); |
688 |
|
handle.logMessage(LOGDEBUG, svel.str().c_str()); |
689 |
|
|
690 |
|
// update particle fitness |
691 |
|
fit[i] = EcoSystems[omp_get_thread_num()]->SimulateAndUpdate(pos[i]); |
692 |
|
// fit[i] = EcoSystem->SimulateAndUpdate(pos[i]); |
693 |
|
// update personal best position? |
694 |
|
// best.bestf = fit[i]; |
695 |
|
if (fit[i] < fit_b[i]) { |
696 |
|
fit_b[i] = fit[i]; |
697 |
|
pos_b[i] = pos[i]; |
698 |
|
} |
699 |
|
// update gbest?? |
700 |
|
handle.logMessage(LOGDEBUG, "Particle ", i); |
701 |
|
handle.logMessage(LOGDEBUG, "The likelihood score is", best.bestf, "at the point"); |
702 |
|
handle.logMessage(LOGDEBUG, "Its likelihood score is", fit[i], "at the point"); |
703 |
|
//#pragma omp critical (best2) |
704 |
|
if (fit[i] < best.bestf) { |
705 |
|
improved = 1; |
706 |
|
// update best fitness |
707 |
|
best.index=i; |
708 |
|
best.bestf = fit[i]; |
709 |
|
for (d = 0; d < nvars; d++) |
710 |
|
bestx[d] = pos[i][d] * init[d]; |
711 |
|
iters = (EcoSystems[0]->getFuncEval() * numThr) + EcoSystem->getFuncEval() - offset; |
712 |
|
handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations"); |
713 |
|
handle.logMessage(LOGINFO, "The likelihood score is", best.bestf, "at the point"); |
714 |
|
EcoSystem->storeVariables(best.bestf, bestx); |
715 |
|
EcoSystem->writeBestValues(); |
716 |
|
bestx = pos[i]; |
717 |
|
|
718 |
|
} |
719 |
|
} // end i |
720 |
|
// if (best.bestf < bestf){ |
721 |
|
// bestf=best.bestf; |
722 |
|
// bestx=pos[best.index]; |
723 |
|
// for (d = 0; d < nvars; d++) |
724 |
|
// bestx[d] = bestx[d] * init[d]; |
725 |
|
// iters = EcoSystem->getFuncEval() - offset; //To change. It is incorrect when is parallel. |
726 |
|
// handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations"); |
727 |
|
// handle.logMessage(LOGINFO, "The likelihood score is", bestf, "at the point"); |
728 |
|
// EcoSystem->storeVariables(bestf, bestx); |
729 |
|
// EcoSystem->writeBestValues(); |
730 |
|
// bestx=pos[best.index]; |
731 |
|
// } |
732 |
|
// handle.logMessage(LOGINFO, "best.bestf=",best.bestf, "\n"); |
733 |
|
psosteps += numThr; |
734 |
|
|
735 |
|
} // while (newf < bestf) |
736 |
|
|
737 |
|
iters = (EcoSystems[0]->getFuncEval() * numThr) + EcoSystem->getFuncEval() - offset; |
738 |
|
handle.logMessage(LOGINFO, "Existing PSO after ", iters, "function evaluations ..."); |
739 |
|
for (d = 0; d < nvars; d++) |
740 |
|
bestx[d] = bestx[d] * init[d]; |
741 |
|
EcoSystem->storeVariables(best.bestf, bestx); |
742 |
|
EcoSystem->writeBestValues(); |
743 |
|
|
744 |
|
} |
745 |
|
#endif |
746 |
|
|
747 |
/** |
/** |
748 |
* \brief calculate constant inertia weight equal to w_max |
* \brief calculate constant inertia weight equal to w_max |
749 |
*/ |
*/ |
757 |
|
|
758 |
int dec_stage = 3 * psoiter / 4; |
int dec_stage = 3 * psoiter / 4; |
759 |
if (step <= dec_stage) |
if (step <= dec_stage) |
760 |
return w_min + (w_max - w_min) * \ |
return w_min + (w_max - w_min) * (dec_stage - step) / dec_stage; |
|
(dec_stage - step) / dec_stage; |
|
761 |
else |
else |
762 |
return w_min; |
return w_min; |
763 |
} |
} |
764 |
|
|
|
|
|
|
|
|
765 |
/** |
/** |
766 |
* \brief NEIGHBORHOOD (COMM) MATRIX STRATEGIES. global neighborhood |
* \brief NEIGHBORHOOD (COMM) MATRIX STRATEGIES. global neighborhood |
767 |
*/ |
*/ |
768 |
void OptInfoPso::inform_global(IntMatrix& comm, DoubleMatrix& pos_nb, |
void OptInfoPso::inform_global(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, DoubleVector& gbest, int improved) { |
|
DoubleMatrix& pos_b, DoubleVector& fit_b, |
|
|
DoubleVector& gbest,int improved) |
|
|
{ |
|
769 |
|
|
770 |
int i; |
int i; |
771 |
// all particles have the same attractor (gbest) |
// all particles have the same attractor (gbest) |
776 |
pos_nb[i]=gbest; |
pos_nb[i]=gbest; |
777 |
} |
} |
778 |
|
|
|
|
|
779 |
/** |
/** |
780 |
* \brief general inform function :: according to the connectivity matrix COMM, it copies the best position (from pos_b) of the informers of each particle to the pos_nb matrix |
* \brief general inform function :: according to the connectivity matrix COMM, it copies the best position (from pos_b) of the informers of each particle to the pos_nb matrix |
781 |
*/ |
*/ |
782 |
void OptInfoPso::inform(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, int improved) |
void OptInfoPso::inform(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, int improved) { |
|
{ |
|
783 |
int i, j; |
int i, j; |
784 |
int b_n; // best neighbor in terms of fitness |
int b_n; // best neighbor in terms of fitness |
785 |
|
|
817 |
comm[i][1] = 1; |
comm[i][1] = 1; |
818 |
// look left |
// look left |
819 |
comm[i][size-1] = 1; |
comm[i][size-1] = 1; |
820 |
} else if (i==size-1) { |
} |
821 |
|
else if (i == size - 1) { |
822 |
// look right |
// look right |
823 |
comm[i][0] = 1; |
comm[i][0] = 1; |
824 |
// look left |
// look left |
825 |
comm[i][i-1] = 1; |
comm[i][i-1] = 1; |
826 |
} else { |
} |
827 |
|
else { |
828 |
// look right |
// look right |
829 |
comm[i][i+1] = 1; |
comm[i][i+1] = 1; |
830 |
// look left |
// look left |
836 |
handle.logMessage(LOGDEBUG, "\nPSO ring communication\n"); |
handle.logMessage(LOGDEBUG, "\nPSO ring communication\n"); |
837 |
for (i=0;i<size;i++){ |
for (i=0;i<size;i++){ |
838 |
string message="["; |
string message="["; |
839 |
for (int j=0;j<size;j++) |
for (int j = 0; j < size; j++) { |
|
{ |
|
840 |
std::ostringstream ss; |
std::ostringstream ss; |
841 |
ss<< comm[i][j]; |
ss<< comm[i][j]; |
842 |
message=message+" "+ss.str(); |
message=message+" "+ss.str(); |
847 |
|
|
848 |
} |
} |
849 |
|
|
850 |
|
void OptInfoPso::inform_ring(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, DoubleVector& gbest, int improved) { |
|
|
|
|
|
|
|
void OptInfoPso::inform_ring(IntMatrix& comm, DoubleMatrix& pos_nb, |
|
|
DoubleMatrix& pos_b, DoubleVector& fit_b, |
|
|
DoubleVector& gbest, int improved) |
|
|
{ |
|
851 |
|
|
852 |
// update pos_nb matrix |
// update pos_nb matrix |
853 |
inform(comm, pos_nb, pos_b, fit_b, improved); |
inform(comm, pos_nb, pos_b, fit_b, improved); |
871 |
for (k=0; k<nhood_size; k++) { |
for (k=0; k<nhood_size; k++) { |
872 |
// generate a random index |
// generate a random index |
873 |
//j = gsl_rng_uniform_int(settings->rng, settings->size); |
//j = gsl_rng_uniform_int(settings->rng, settings->size); |
874 |
j=rand()*size/RAND_MAX; |
j = rand_r(&this->seed) * size / RAND_MAX; |
875 |
// particle i informs particle j |
// particle i informs particle j |
876 |
comm[i] [j] = 1; |
comm[i] [j] = 1; |
877 |
} |
} |
878 |
} |
} |
879 |
} |
} |
880 |
|
|
881 |
|
void OptInfoPso::inform_random(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, DoubleVector& gbest, int improved) { |
|
|
|
|
void OptInfoPso::inform_random(IntMatrix& comm, DoubleMatrix& pos_nb, |
|
|
DoubleMatrix& pos_b, DoubleVector& fit_b, |
|
|
DoubleVector& gbest, int improved) |
|
|
{ |
|
|
|
|
882 |
|
|
883 |
// regenerate connectivity?? |
// regenerate connectivity?? |
884 |
if (!improved) |
if (!improved) init_comm_random(comm); |
|
init_comm_random(comm); |
|
885 |
inform(comm, pos_nb, pos_b, fit_b, improved); |
inform(comm, pos_nb, pos_b, fit_b, improved); |
886 |
|
|
887 |
} |
} |