Annotation of /trunk/gadget/pso.cc

Parent Directory | Revision Log

---

Revision 28 - (view) (download)

1 :	agomez	22	#include "gadget.h"
2 :			#include "optinfo.h"
3 :			#include "mathfunc.h"
4 :			#include "doublevector.h"
5 :			#include "intvector.h"
6 :			#include "errorhandler.h"
7 :			#include "ecosystem.h"
8 :			#include "global.h"
9 :			#include "pso.h"
10 :	ulcessvp	28	#include <float.h>
11 :	agomez	22
12 :			#ifdef _OPENMP
13 :			#include "omp.h"
14 :			#endif
15 :
16 :			extern Ecosystem* EcoSystem;
17 :			#ifdef _OPENMP
18 :			extern Ecosystem** EcoSystems;
19 :			#endif
20 :
21 :			//==============================================================
22 :			//// calulate swarm size based on dimensionality
23 :			int OptInfoPso::pso_calc_swarm_size(int dim) {
24 :			int size = 10. + 2. * sqrt(dim);
25 :	ulcessvp	28	return (size > PSO_MAX_SIZE ? PSO_MAX_SIZE : size);
26 :			}
27 :	agomez	22
28 :			void OptInfoPso::OptimiseLikelihood() {
29 :
30 :	ulcessvp	28	handle.logMessage(LOGINFO, "\nStarting PSO optimisation algorithm\n");
31 :			double bestf, tmp;
32 :			int i, offset, rchange, rcheck, rnumber;
33 :			int nvars = EcoSystem->numOptVariables();
34 :			DoubleVector x(nvars);
35 :			// DoubleVector trialx(nvars);
36 :			DoubleVector bestx(nvars);
37 :			DoubleVector lowerb(nvars);
38 :			DoubleVector upperb(nvars);
39 :			DoubleVector init(nvars);
40 :			//DoubleVector initialstep(nvars, rho);
41 :			IntVector param(nvars, 0);
42 :			IntVector lbound(nvars, 0);
43 :			IntVector rbounds(nvars, 0);
44 :			IntVector trapped(nvars, 0);
45 :	agomez	22
46 :	ulcessvp	28	/**
47 :			* 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
52 :			DoubleMatrix vel(size, nvars, 0.0); // velocity matrix
53 :			DoubleMatrix pos_b(size, nvars, 0.0); // best position matrix
54 :			DoubleVector fit(size); // particle fitness vector
55 :			DoubleVector fit_b(size); // best fitness vector
56 :			DoubleMatrix pos_nb(size, nvars, 0.0); // what is the best informed position for each particle
57 :			IntMatrix comm(size, size, 0); // communications:who informs who rows : those who inform cols : those who are informed
58 :			int improved; // whether solution->error was improved during the last iteration
59 :			int d, step;
60 :			double a, b;
61 :			double rho1, rho2; //random numbers (coefficients)
62 :			double w; //current omega
63 :			Inform_fun inform_fun;
64 :			Calc_inertia_fun calc_inertia_fun; // inertia weight update function
65 :			int psosteps = 0;
66 :	agomez	22
67 :	ulcessvp	28	if (seed == 0) seed = 1234;
68 :	agomez	22
69 :	ulcessvp	28	handle.logMessage(LOGINFO, "Starting PSO with particles ", size, "\n");
70 :	agomez	22
71 :	ulcessvp	28	// SELECT APPROPRIATE NHOOD UPDATE FUNCTION
72 :			switch (nhood_strategy) {
73 :			case PSO_NHOOD_GLOBAL:
74 :			// comm matrix not used
75 :			handle.logMessage(LOGINFO, "PSO will use global communication\n");
76 :			inform_fun = &OptInfoPso::inform_global;
77 :			break;
78 :			case PSO_NHOOD_RING:
79 :			handle.logMessage(LOGINFO, "PSO will use ring communication\n");
80 :			init_comm_ring(comm);
81 :			inform_fun = &OptInfoPso::inform_ring;
82 :			break;
83 :			case PSO_NHOOD_RANDOM:
84 :			handle.logMessage(LOGINFO, "PSO will use random communication\n");
85 :			init_comm_random(comm);
86 :			inform_fun = &OptInfoPso::inform_random;
87 :			break;
88 :			}
89 :			// SELECT APPROPRIATE INERTIA WEIGHT UPDATE FUNCTION
90 :			switch (w_strategy) {
91 :			case PSO_W_CONST:
92 :			handle.logMessage(LOGINFO, "PSO will use constant inertia\n");
93 :			calc_inertia_fun = &OptInfoPso::calc_inertia_const;
94 :			break;
95 :			case PSO_W_LIN_DEC:
96 :			handle.logMessage(LOGINFO, "PSO will use linear decreasing inertia\n");
97 :			calc_inertia_fun = &OptInfoPso::calc_inertia_lin_dec;
98 :			break;
99 :			}
100 :			//
101 :			//initialize omega using standard value
102 :			//
103 :			w = PSO_INERTIA;
104 :			handle.logMessage(LOGINFO, "PSO initial inertia ", w, "\n");
105 :			handle.logMessage(LOGINFO, "PSO c1", c1, "\n");
106 :			handle.logMessage(LOGINFO, "PSO c2", c2, "\n");
107 :			handle.logMessage(LOGINFO, "PSO goal", goal, "\n");
108 :	agomez	22
109 :	ulcessvp	28	// if (scale)
110 :			EcoSystem->scaleVariables();
111 :			//#ifdef _OPENMP
112 :			// int numThr = omp_get_max_threads ( );
113 :			//#pragma omp parallel for
114 :			// for (i=0;i<numThr;i++)
115 :			// EcoSystems[i]->scaleVariables();
116 :			//#endif
117 :			EcoSystem->getOptScaledValues(x);
118 :			EcoSystem->getOptLowerBounds(lowerb);
119 :			EcoSystem->getOptUpperBounds(upperb);
120 :			EcoSystem->getOptInitialValues(init);
121 :
122 :			for (i = 0; i < nvars; ++i) {
123 :			bestx[i] = x[i];
124 :			param[i] = i;
125 :	agomez	22	}
126 :
127 :	ulcessvp	28	//se escala?
128 :			// if (scale)
129 :			for (i = 0; i < nvars; i++) {
130 :			// Scaling the bounds, because the parameters are scaled
131 :			lowerb[i] = lowerb[i] / init[i];
132 :			upperb[i] = upperb[i] / init[i];
133 :			if (lowerb[i] > upperb[i]) {
134 :			tmp = lowerb[i];
135 :			lowerb[i] = upperb[i];
136 :			upperb[i] = tmp;
137 :			}
138 :			}
139 :	agomez	22
140 :	ulcessvp	28	bestf = EcoSystem->SimulateAndUpdate(x);
141 :			if (bestf != bestf) { //check for NaN
142 :			handle.logMessage(LOGINFO, "Error starting PSO optimisation with f(x) = infinity");
143 :			converge = -1;
144 :			iters = 1;
145 :			return;
146 :			}
147 :	agomez	22
148 :	ulcessvp	28	offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop
149 :			// newf = bestf;
150 :			// oldf = bestf;
151 :
152 :			iters = 0;
153 :
154 :	agomez	22	//
155 :	ulcessvp	28	// Initialize the particles
156 :			//
157 :			i = 0;
158 :			handle.logMessage(LOGINFO, "Initialising PSO ", size, "particles\n");
159 :
160 :			// SWARM INITIALIZATION
161 :			//
162 :			// Initialize the fist particle with the best known position
163 :			//
164 :			for (d = 0; d < nvars; d++) {
165 :			// generate two numbers within the specified range and around trialx
166 :			a = x[d];
167 :			// initialize position
168 :			pos[i][d] = a;
169 :			pos_b[i][d] = a;
170 :			pos_nb[i][d] = a;
171 :			// a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand() * 1.0 / RAND_MAX);
172 :			b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX);
173 :			vel[i][d] = (a - b) / 2.;
174 :			}
175 :			fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
176 :			fit_b[i] = fit[i]; // this is also the personal best
177 :			if (fit[i] < bestf) {
178 :			bestf = fit[i];
179 :			bestx = pos[i];
180 :			}
181 :
182 :	agomez	22	// Now other particles with a random position in the solution space
183 :			//
184 :	ulcessvp	28	//#pragma omp parallel for default(shared) private(i,d,a,b)
185 :			for (i = 1; i < size; i++) {
186 :			ostringstream spos, svel;
187 :			spos << "Initial position of particle " << i << " [";
188 :			svel << "Initial velocity of particle " << i << " [";
189 :			// for each dimension
190 :			for (d = 0; d < nvars; d++) {
191 :			// generate two numbers within the specified range and around trialx
192 :			a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX);
193 :			b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand_r(&seed) * 1.0 / RAND_MAX);
194 :			// initialize position
195 :			pos[i][d] = a;
196 :			// best position is trialx
197 :			pos_b[i][d] = a;
198 :			pos_nb[i][d] = a;
199 :			// initialize velocity
200 :			vel[i][d] = (a - b) / 2.;
201 :	agomez	22	spos << pos[i][d] << " ";
202 :			svel << vel[i][d] << " ";
203 :	ulcessvp	28	}
204 :			//#ifdef _OPENMP
205 :			// fit[i] = EcoSystems[omp_get_thread_num ()]->SimulateAndUpdate(pos[i]);
206 :			//#else
207 :			fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
208 :			//#endif
209 :			fit_b[i] = fit[i]; // this is also the personal best
210 :			// update gbest??
211 :	agomez	22	//#pragma omp critical (best1)
212 :	ulcessvp	28	if (fit[i] < bestf) {
213 :			// update best fitness
214 :			bestf = fit[i];
215 :			// copy particle pos to gbest vector
216 :			//memmove((void *)&bestx, (void *)&pos[i],
217 :			// sizeof(double) * nvars);
218 :			bestx = pos[i];
219 :			}
220 :			spos << " ]";
221 :			svel << " ]";
222 :	agomez	22	handle.logMessage(LOGDEBUG, spos.str().c_str());
223 :			handle.logMessage(LOGDEBUG, svel.str().c_str());
224 :
225 :	ulcessvp	28	}
226 :
227 :	agomez	22	//
228 :	ulcessvp	28	// RUN ALGORITHM
229 :	agomez	22	//
230 :	ulcessvp	28	psosteps = 0;
231 :			while (1) {
232 :	agomez	22
233 :	ulcessvp	28	handle.logMessage(LOGINFO, "PSO optimisation after", psosteps * size, "\n");
234 :			if (isZero(bestf)) {
235 :			iters = EcoSystem->getFuncEval() - offset;
236 :			handle.logMessage(LOGINFO, "Error in PSO optimisation after", iters, "function evaluations, f(x) = 0");
237 :			converge = -1;
238 :			return;
239 :			}
240 :			/* if too many function evaluations occur, terminate the algorithm */
241 :	agomez	22
242 :	ulcessvp	28	iters = EcoSystem->getFuncEval() - offset;
243 :			if (psosteps > psoiter) {
244 :			handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n");
245 :			handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
246 :			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");
248 :	agomez	22
249 :	ulcessvp	28	score = EcoSystem->SimulateAndUpdate(bestx);
250 :			for (i = 0; i < nvars; i++)
251 :			bestx[i] = bestx[i] * init[i];
252 :			handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", score);
253 :			EcoSystem->storeVariables(score, bestx);
254 :			return;
255 :			}
256 :			// update inertia weight
257 :			if (w_strategy) w = (this->*calc_inertia_fun)(psosteps);
258 :			// check optimization goal
259 :			if (bestf <= goal) {
260 :			handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n");
261 :			handle.logMessage(LOGINFO, "Goal achieved!!! @ step ", psosteps);
262 :			score = EcoSystem->SimulateAndUpdate(bestx);
263 :			handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", score);
264 :			break;
265 :			}
266 :	agomez	22
267 :	ulcessvp	28	// update pos_nb matrix (find best of neighborhood for all particles)
268 :			(this->*inform_fun)(comm, pos_nb, pos_b, fit_b, bestx, improved);
269 :	agomez	22
270 :	ulcessvp	28	// the value of improved was just used; reset it
271 :			improved = 0;
272 :	agomez	22
273 :	ulcessvp	28	// update all particles
274 :			//#pragma omp parallel for private(d,rho1,rho2) default(shared)
275 :			for (i = 0; i < size; i++) {
276 :			// handle.logMessage(LOGDEBUG, "In parallel",omp_get_thread_num(),"\n");
277 :			ostringstream spos, svel;
278 :			spos << "New position of particle " << i << " [";
279 :			svel << "New velocity of particle " << i << " [";
280 :			// for each dimension
281 :			for (d = 0; d < nvars; d++) {
282 :			// calculate stochastic coefficients
283 :			rho1 = c1 * (rand_r(&seed) * 1.0) / RAND_MAX;
284 :			rho2 = c2 * (rand_r(&seed) * 1.0) / RAND_MAX;
285 :			// update velocity
286 :			vel[i][d] = w * vel[i][d] + rho1 * (pos_b[i][d] - pos[i][d]) + rho2 * (pos_nb[i][d] - pos[i][d]);
287 :			// update position
288 :			pos[i][d] += vel[i][d];
289 :			spos << pos[i][d] << " ";
290 :			svel << vel[i][d] << " ";
291 :			// clamp position within bounds?
292 :			if (clamp_pos) {
293 :			if (pos[i][d] < lowerb[d]) {
294 :			pos[i][d] = lowerb[d];
295 :			vel[i][d] = 0;
296 :			}
297 :			else if (pos[i][d] > upperb[d]) {
298 :			pos[i][d] = upperb[d];
299 :			vel[i][d] = 0;
300 :			}
301 :			}
302 :			else {
303 :			// enforce periodic boundary conditions
304 :			if (pos[i][d] < lowerb[d]) {
305 :			pos[i][d] = upperb[d] - fmod(lowerb[d] - pos[i][d], upperb[d] - lowerb[d]);
306 :			vel[i][d] = 0;
307 :			}
308 :			else if (pos[i][d] > upperb[d]) {
309 :			pos[i][d] = lowerb[d] + fmod(pos[i][d] - upperb[d], upperb[d] - lowerb[d]);
310 :			vel[i][d] = 0;
311 :	agomez	22
312 :	ulcessvp	28	}
313 :			}
314 :			} //end for nvars
315 :	agomez	22
316 :	ulcessvp	28	spos << "]" << endl;
317 :			svel << "]" << endl;
318 :			handle.logMessage(LOGDEBUG, spos.str().c_str());
319 :			handle.logMessage(LOGDEBUG, svel.str().c_str());
320 :	agomez	22
321 :	ulcessvp	28	// update particle fitness
322 :			// fit[i] = EcoSystems[omp_get_thread_num()]->SimulateAndUpdate(pos[i]);
323 :			fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
324 :			// update personal best position?
325 :			if (fit[i] < fit_b[i]) {
326 :			fit_b[i] = fit[i];
327 :			pos_b[i] = pos[i];
328 :			}
329 :			// update gbest??
330 :			handle.logMessage(LOGDEBUG, "Particle ", i);
331 :			handle.logMessage(LOGDEBUG, "The likelihood score is", bestf, "at the point");
332 :			handle.logMessage(LOGDEBUG, "Its likelihood score is", fit[i], "at the point");
333 :			//#pragma omp critical (best2)
334 :			if (fit[i] < bestf) {
335 :			improved = 1;
336 :			// update best fitness
337 :			bestf = fit[i];
338 :			for (d = 0; d < nvars; d++)
339 :			bestx[d] = pos[i][d] * init[d];
340 :			iters = EcoSystem->getFuncEval() - offset;
341 :			handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations");
342 :			handle.logMessage(LOGINFO, "The likelihood score is", bestf, "at the point");
343 :			EcoSystem->storeVariables(bestf, bestx);
344 :			EcoSystem->writeBestValues();
345 :			bestx = pos[i];
346 :	agomez	22
347 :	ulcessvp	28	}
348 :			} // end i
349 :	agomez	22
350 :	ulcessvp	28	psosteps++;
351 :	agomez	22
352 :			} // while (newf < bestf)
353 :
354 :			iters = EcoSystem->getFuncEval() - offset;
355 :			handle.logMessage(LOGINFO, "Existing PSO after ", iters, "function evaluations ...");
356 :	ulcessvp	28	for (d = 0; d < nvars; d++)
357 :			bestx[d] = bestx[d] * init[d];
358 :	agomez	22	EcoSystem->storeVariables(bestf, bestx);
359 :			EcoSystem->writeBestValues();
360 :
361 :			}
362 :
363 :	ulcessvp	28	#ifdef _OPENMP
364 :	agomez	22	void OptInfoPso::OptimiseLikelihoodREP() {
365 :	ulcessvp	28	OptimiseLikelihood();
366 :	agomez	22
367 :			}
368 :
369 :			void OptInfoPso::OptimiseLikelihoodOMP() {
370 :
371 :	ulcessvp	28	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 :	agomez	22	/**
748 :			* \brief calculate constant inertia weight equal to w_max
749 :			*/
750 :			double OptInfoPso::calc_inertia_const(int step) {
751 :	ulcessvp	28	return w_max;
752 :	agomez	22	}
753 :			/**
754 :			* \brief calculate linearly decreasing inertia weight
755 :			*/
756 :			double OptInfoPso::calc_inertia_lin_dec(int step) {
757 :
758 :			int dec_stage = 3 * psoiter / 4;
759 :			if (step <= dec_stage)
760 :	ulcessvp	28	return w_min + (w_max - w_min) * (dec_stage - step) / dec_stage;
761 :	agomez	22	else
762 :	ulcessvp	28	return w_min;
763 :	agomez	22	}
764 :
765 :			/**
766 :			* \brief NEIGHBORHOOD (COMM) MATRIX STRATEGIES. global neighborhood
767 :			*/
768 :	ulcessvp	28	void OptInfoPso::inform_global(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, DoubleVector& gbest, int improved) {
769 :	agomez	22
770 :			int i;
771 :	ulcessvp	28	// all particles have the same attractor (gbest)
772 :			// copy the contents of gbest to pos_nb
773 :			for (i = 0; i < size; i++)
774 :			//memmove((void *)&pos_nb[i][0], &gbest[0],
775 :	agomez	22	// sizeof(double) * gbest.Size());
776 :	ulcessvp	28	pos_nb[i] = gbest;
777 :	agomez	22	}
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
781 :			*/
782 :	ulcessvp	28	void OptInfoPso::inform(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, int improved) {
783 :	agomez	22	int i, j;
784 :			int b_n; // best neighbor in terms of fitness
785 :
786 :	ulcessvp	28	// for each particle
787 :			for (j = 0; j < size; j++) {
788 :			b_n = j; // self is best
789 :			// who is the best informer??
790 :			for (i = 0; i < size; i++)
791 :			// the i^th particle informs the j^th particle
792 :			if (comm[i][j] && fit_b[i] < fit_b[b_n])
793 :			// found a better informer for j^th particle
794 :			b_n = i;
795 :			// copy pos_b of b_n^th particle to pos_nb[j]
796 :			//memmove((void *)&pos_nb[j][0],
797 :			// (void *)&pos_b[b_n][0],
798 :			// sizeof(double) * pos_nb.Ncol(j));
799 :			pos_nb[j] = pos_b[b_n];
800 :			}
801 :	agomez	22	}
802 :
803 :			/**
804 :			* \brief topology initialization :: this is a static (i.e. fixed) topology
805 :			*/
806 :			void OptInfoPso::init_comm_ring(IntMatrix& comm) {
807 :			int i;
808 :	ulcessvp	28	// reset array
809 :			comm.setToZero();
810 :	agomez	22
811 :	ulcessvp	28	// choose informers
812 :			for (i = 0; i < size; i++) {
813 :			// set diagonal to 1
814 :			comm[i][i] = 1;
815 :			if (i == 0) {
816 :			// look right
817 :			comm[i][1] = 1;
818 :			// look left
819 :			comm[i][size - 1] = 1;
820 :			}
821 :			else if (i == size - 1) {
822 :			// look right
823 :			comm[i][0] = 1;
824 :			// look left
825 :			comm[i][i - 1] = 1;
826 :			}
827 :			else {
828 :			// look right
829 :			comm[i][i + 1] = 1;
830 :			// look left
831 :			comm[i][i - 1] = 1;
832 :			}
833 :	agomez	22
834 :			}
835 :	ulcessvp	28	// Print Matrix
836 :	agomez	22	handle.logMessage(LOGDEBUG, "\nPSO ring communication\n");
837 :	ulcessvp	28	for (i = 0; i < size; i++) {
838 :			string message = "[";
839 :			for (int j = 0; j < size; j++) {
840 :			std::ostringstream ss;
841 :			ss << comm[i][j];
842 :			message = message + " " + ss.str();
843 :			}
844 :			message = message + "]\n";
845 :			handle.logMessage(LOGDEBUG, message.c_str());
846 :	agomez	22	}
847 :
848 :			}
849 :
850 :	ulcessvp	28	void OptInfoPso::inform_ring(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, DoubleVector& gbest, int improved) {
851 :	agomez	22
852 :	ulcessvp	28	// update pos_nb matrix
853 :	agomez	22	inform(comm, pos_nb, pos_b, fit_b, improved);
854 :
855 :			}
856 :
857 :			/**
858 :			* \brief random neighborhood topology
859 :			*/
860 :			void OptInfoPso::init_comm_random(IntMatrix& comm) {
861 :
862 :			int i, j, k;
863 :	ulcessvp	28	// reset array
864 :	agomez	22	comm.setToZero();
865 :
866 :	ulcessvp	28	// choose informers
867 :			for (i = 0; i < size; i++) {
868 :			// each particle informs itself
869 :			comm[i][i] = 1;
870 :			// choose kappa (on average) informers for each particle
871 :			for (k = 0; k < nhood_size; k++) {
872 :			// generate a random index
873 :			//j = gsl_rng_uniform_int(settings->rng, settings->size);
874 :			j = rand_r(&this->seed) * size / RAND_MAX;
875 :			// particle i informs particle j
876 :			comm[i][j] = 1;
877 :			}
878 :	agomez	22	}
879 :			}
880 :
881 :	ulcessvp	28	void OptInfoPso::inform_random(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, DoubleVector& gbest, int improved) {
882 :	agomez	22
883 :	ulcessvp	28	// regenerate connectivity??
884 :			if (!improved) init_comm_random(comm);
885 :	agomez	22	inform(comm, pos_nb, pos_b, fit_b, improved);
886 :
887 :			}
888 :	ulcessvp	28

---

root@forge.cesga.es	ViewVC Help
Powered by ViewVC 1.0.0