Annotation of /trunk/gadget/pso.cc

Parent Directory | Revision Log

---

Revision 22 - (view) (download)

1 :	agomez	22
2 :			#include "gadget.h"
3 :			#include "optinfo.h"
4 :			#include "mathfunc.h"
5 :			#include "doublevector.h"
6 :			#include "intvector.h"
7 :			#include "errorhandler.h"
8 :			#include "ecosystem.h"
9 :			#include "global.h"
10 :			#include "pso.h"
11 :
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 :			//==============================================================
23 :			//// calulate swarm size based on dimensionality
24 :			int OptInfoPso::pso_calc_swarm_size(int dim) {
25 :			int size = 10. + 2. * sqrt(dim);
26 :			return (size > PSO_MAX_SIZE ? PSO_MAX_SIZE : size);
27 :			}
28 :
29 :			void OptInfoPso::OptimiseLikelihood() {
30 :
31 :			handle.logMessage(LOGINFO, "\nStarting PSO optimisation algorithm\n");
32 :			double oldf, newf, steplength, tmp;
33 :			int i, offset,rchange,rcheck,rnumber;
34 :			int nvars = EcoSystem->numOptVariables();
35 :			DoubleVector x(nvars);
36 :			DoubleVector trialx(nvars);
37 :			DoubleVector bestx(nvars);
38 :			struct Best {double bestf; int index;};
39 :			struct Best best;
40 :			double bestf;
41 :			DoubleVector lowerb(nvars);
42 :			DoubleVector upperb(nvars);
43 :			DoubleVector init(nvars);
44 :			//DoubleVector initialstep(nvars, rho);
45 :			IntVector param(nvars, 0);
46 :			IntVector lbound(nvars, 0);
47 :			IntVector rbounds(nvars, 0);
48 :			IntVector trapped(nvars, 0);
49 :
50 :			/**
51 :			* Internal variables for PSO
52 :			*/
53 :			DoubleMatrix pos(size,nvars,0.0); // position matrix
54 :			DoubleMatrix vel(size,nvars,0.0); // velocity matrix
55 :			DoubleMatrix pos_b(size,nvars,0.0); // best position matrix
56 :			DoubleVector fit(size); // particle fitness vector
57 :			DoubleVector fit_b(size); // best fitness vector
58 :			DoubleMatrix pos_nb(size,nvars,0.0); // what is the best informed position for each particle
59 :			IntMatrix comm(size,size,0); // communications:who informs who rows : those who inform cols : those who are informed
60 :			int improved; // whether solution->error was improved during the last iteration
61 :			int d, step;
62 :			double a, b;
63 :			double rho1, rho2; //random numbers (coefficients)
64 :			double w; //current omega
65 :			Inform_fun inform_fun;
66 :			Calc_inertia_fun calc_inertia_fun; // inertia weight update function
67 :			int psosteps=0;
68 :
69 :			if (this->seed!=0)
70 :			srand(this->seed);
71 :
72 :			handle.logMessage(LOGINFO,"Starting PSO with particles ",size,"\n");
73 :			switch (nhood_strategy)
74 :			{
75 :			case PSO_NHOOD_GLOBAL :
76 :			// comm matrix not used
77 :			handle.logMessage(LOGINFO,"PSO will use global communication\n");
78 :			inform_fun = &OptInfoPso::inform_global;
79 :			break;
80 :			case PSO_NHOOD_RING :
81 :			handle.logMessage(LOGINFO,"PSO will use ring communication\n");
82 :			init_comm_ring(comm);
83 :			inform_fun = &OptInfoPso::inform_ring;
84 :			break;
85 :			case PSO_NHOOD_RANDOM :
86 :			handle.logMessage(LOGINFO,"PSO will use random communication\n");
87 :			init_comm_random(comm);
88 :			inform_fun = &OptInfoPso::inform_random;
89 :			break;
90 :			}
91 :			// SELECT APPROPRIATE INERTIA WEIGHT UPDATE FUNCTION
92 :			switch (w_strategy)
93 :			{
94 :			case PSO_W_CONST :
95 :			handle.logMessage(LOGINFO,"PSO will use constant inertia\n");
96 :			calc_inertia_fun = &OptInfoPso::calc_inertia_const;
97 :			break;
98 :			case PSO_W_LIN_DEC :
99 :			handle.logMessage(LOGINFO,"PSO will use linear decreasing inertia\n");
100 :			calc_inertia_fun = &OptInfoPso::calc_inertia_lin_dec;
101 :			break;
102 :			}
103 :			//
104 :			//initialize omega using standard value
105 :			//
106 :			w = PSO_INERTIA;
107 :			handle.logMessage(LOGINFO,"PSO initial inertia ",w,"\n");
108 :			handle.logMessage(LOGINFO,"PSO c1",c1,"\n");
109 :			handle.logMessage(LOGINFO,"PSO c2",c2,"\n");
110 :			handle.logMessage(LOGINFO,"PSO goal",goal,"\n");
111 :
112 :			EcoSystem->scaleVariables();
113 :			#ifdef _OPENMP
114 :			int numThr = omp_get_max_threads ( );
115 :			#pragma omp parallel for
116 :			for (i=0;i<numThr;i++)
117 :			EcoSystems[i]->scaleVariables();
118 :			#endif
119 :			EcoSystem->getOptScaledValues(x);
120 :			EcoSystem->getOptLowerBounds(lowerb);
121 :			EcoSystem->getOptUpperBounds(upperb);
122 :			EcoSystem->getOptInitialValues(init);
123 :
124 :			for (i = 0; i < nvars; i++) {
125 :			// Scaling the bounds, because the parameters are scaled
126 :			lowerb[i] = lowerb[i] / init[i];
127 :			upperb[i] = upperb[i] / init[i];
128 :			if (lowerb[i] > upperb[i]) {
129 :			tmp = lowerb[i];
130 :			lowerb[i] = upperb[i];
131 :			upperb[i] = tmp;
132 :			}
133 :			bestx[i] = x[i];
134 :			trialx[i] = x[i];
135 :			param[i] = i;
136 :			}
137 :
138 :			bestf = EcoSystem->SimulateAndUpdate(trialx);
139 :			if (bestf != bestf) { //check for NaN
140 :			handle.logMessage(LOGINFO, "Error starting PSO optimisation with f(x) = infinity");
141 :			converge = -1;
142 :			iters = 1;
143 :			return;
144 :			}
145 :
146 :			offset = EcoSystem->getFuncEval(); //number of function evaluations done before loop
147 :
148 :			iters = 0;
149 :
150 :			//
151 :			// Initialize the particles
152 :			//
153 :			i=0;
154 :			handle.logMessage(LOGINFO,"Initialising PSO ",size,"particles\n");
155 :			//
156 :			// Initialize the fist particle with the best known position
157 :			//
158 :			for (d=0; d<nvars; d++) {
159 :			// generate two numbers within the specified range and around trialx
160 :			a = trialx[d] ;
161 :			// initialize position
162 :			pos[i][d] = a;
163 :			pos_b[i][d] = a;
164 :			pos_nb[i][d] = a;
165 :			a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX);
166 :			b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX);
167 :			vel[i][d] = (a-b) / 2.;
168 :			}
169 :			fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
170 :			fit_b[i] = fit[i]; // this is also the personal best
171 :			if (fit[i] < bestf) {
172 :			bestf = fit[i];
173 :			best.bestf=bestf;
174 :			bestx=pos[i];
175 :			best.index=i;
176 :			}
177 :			//
178 :			// Now other particles with a random position in the solution space
179 :			//
180 :			best.bestf=bestf;
181 :			best.index=0;
182 :			handle.logMessage(LOGINFO,"before Best.bestf=",best.bestf);
183 :			#pragma omp declare reduction(bestomp: struct Best : omp_out = omp_in.bestf < omp_out.bestf ? omp_in : omp_out)
184 :			#pragma omp parallel for default(shared) private(i,d,a,b) reduction(bestomp:best) firstprivate(bestf)
185 :			for (i=1; i<size; i++) {
186 :			// for each dimension
187 :			ostringstream spos,svel;
188 :			spos << "Initial position of particle " << i<<" [" ;
189 :			svel << "Initial velocity of particle " << i<<" [" ;
190 :
191 :			for (d=0; d<nvars; d++) {
192 :			// generate two numbers within the specified range and around trialx
193 :			a = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX);
194 :			b = lowerb[d] + (upperb[d] - lowerb[d]) * (rand()*1.0/RAND_MAX);
195 :			// initialize position
196 :			pos[i][d] = a;
197 :			// best position is trialx
198 :			pos_b[i][d] = a;
199 :			pos_nb[i][d] = a;
200 :			// initialize velocity
201 :			vel[i][d] = (a-b) / 2.;
202 :			spos << pos[i][d] << " ";
203 :			svel << vel[i][d] << " ";
204 :			}
205 :			#ifdef _OPENMP
206 :			fit[i] = EcoSystems[omp_get_thread_num ()]->SimulateAndUpdate(pos[i]);
207 :			best.bestf = fit[i];
208 :			#else
209 :			fit[i] = EcoSystem->SimulateAndUpdate(pos[i]);
210 :			#endif
211 :			fit_b[i] = fit[i]; // this is also the personal best
212 :			// update gbest??
213 :			//#pragma omp critical (best1)
214 :			if (fit[i] < bestf) {
215 :			// update best fitness
216 :			handle.logMessage(LOGINFO,"iInside Best.bestf=",best.bestf);
217 :			// copy particle pos to gbest vector
218 :			//memmove((void *)&bestx, (void *)&pos[i],
219 :			// sizeof(double) * nvars);
220 :			best.index=i;
221 :			bestf=best.bestf;
222 :			//bestx=pos[i];
223 :			}
224 :			spos<<" ]";
225 :			svel<<" ]";
226 :			handle.logMessage(LOGDEBUG, spos.str().c_str());
227 :			handle.logMessage(LOGDEBUG, svel.str().c_str());
228 :
229 :
230 :			}
231 :			//
232 :			// Compare the reduction with the best value
233 :			//
234 :
235 :			if (best.bestf < bestf){
236 :			bestx=pos[best.index];
237 :			bestf=best.bestf;
238 :			}
239 :
240 :
241 :			//
242 :			// RUN ALGORITHM
243 :			//
244 :			psosteps=0;
245 :			while (1) {
246 :
247 :			handle.logMessage(LOGINFO, "PSO optimisation after", psosteps*size,"\n");
248 :			if (isZero(bestf)) {
249 :			iters = EcoSystem->getFuncEval() - offset;
250 :			handle.logMessage(LOGINFO, "Error in PSO optimisation after", iters, "function evaluations, f(x) = 0");
251 :			converge = -1;
252 :			return;
253 :			}
254 :			/* if too many function evaluations occur, terminate the algorithm */
255 :
256 :			iters = EcoSystem->getFuncEval() - offset;
257 :			if (psosteps > psoiter) {
258 :			handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n");
259 :			handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
260 :			handle.logMessage(LOGINFO, "The optimisation stopped because the maximum number of PSO steps was reached");
261 :			handle.logMessage(LOGINFO, "was reached and NOT because an optimum was found for this run");
262 :
263 :			bestf = EcoSystem->SimulateAndUpdate(bestx);
264 :			for (i = 0; i < nvars; i++)
265 :			bestx[i] = bestx[i] * init[i];
266 :			handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", bestf);
267 :			EcoSystem->storeVariables(bestf, bestx);
268 :			return;
269 :			}
270 :
271 :			if (w_strategy)
272 :			w = (this->*calc_inertia_fun)(psosteps);
273 :			// check optimization goal
274 :			if (bestf <= goal) {
275 :			handle.logMessage(LOGINFO, "\nStopping PSO optimisation algorithm\n");
276 :			handle.logMessage(LOGINFO, "Goal achieved!!! @ step ", psosteps);
277 :			bestf = EcoSystem->SimulateAndUpdate(bestx);
278 :			handle.logMessage(LOGINFO, "\nPSO finished with a likelihood score of", bestf);
279 :			break;
280 :			}
281 :
282 :			// update pos_nb matrix (find best of neighborhood for all particles)
283 :			(this->*inform_fun)(comm, pos_nb, pos_b, fit_b, bestx, improved);
284 :
285 :			// the value of improved was just used; reset it
286 :			improved = 0;
287 :
288 :			// update all particles
289 :			best.bestf=bestf;
290 :			#pragma omp parallel for default(shared) private(i,d,rho1,rho2) reduction(bestomp:best) firstprivate(bestf)
291 :			for (i=0; i<size; i++) {
292 :			// for each dimension
293 :			handle.logMessage(LOGDEBUG, "In parallel",omp_get_thread_num(),"\n");
294 :			ostringstream spos,svel;
295 :			spos << "New position of particle " << i<<" [" ;
296 :			svel << "New velocity of particle " << i<<" [" ;
297 :			for (d=0; d<nvars; d++) {
298 :			// calculate stochastic coefficients
299 :			rho1 = c1 * (rand()*1.0)/RAND_MAX;
300 :			rho2 = c2 * (rand()*1.0)/RAND_MAX;
301 :			// update velocity
302 :			vel[i][d] = w * vel[i][d] + \
303 :			rho1 * (pos_b[i][d] - pos[i][d]) + \
304 :			rho2 * (pos_nb[i][d] - pos[i][d]);
305 :			// update position
306 :			pos[i][d] += vel[i][d];
307 :			spos << pos[i][d] <<" ";
308 :			svel << vel[i][d] <<" ";
309 :			// clamp position within bounds?
310 :			if (clamp_pos) {
311 :			if (pos[i][d] < lowerb[d]) {
312 :			pos[i][d] = lowerb[d];
313 :			vel[i][d] = 0;
314 :			} else if (pos[i][d] > upperb[d]) {
315 :			pos[i][d] = upperb[d];
316 :			vel[i][d] = 0;
317 :			}
318 :			} else {
319 :			// enforce periodic boundary conditions
320 :			if (pos[i][d] < lowerb[d]) {
321 :
322 :			pos[i][d] = upperb[d] - fmod(lowerb[d] - pos[i][d],
323 :			upperb[d] - lowerb[d]);
324 :			vel[i][d] = 0;
325 :
326 :			} else if (pos[i][d] > upperb[d]) {
327 :
328 :			pos[i][d] = lowerb[d] + fmod(pos[i][d] - upperb[d],
329 :			upperb[d] - lowerb[d]);
330 :			vel[i][d] = 0;
331 :
332 :
333 :			}
334 :			}
335 :
336 :			} //end for d
337 :			spos<< "]" << endl;
338 :			svel<< "]" << endl;
339 :			handle.logMessage(LOGDEBUG,spos.str().c_str());
340 :			handle.logMessage(LOGDEBUG,svel.str().c_str());
341 :
342 :			// update particle fitness
343 :			fit[i] = EcoSystems[omp_get_thread_num()]->SimulateAndUpdate(pos[i]);
344 :			// update personal best position?
345 :			best.bestf = fit[i];
346 :			if (fit[i] < fit_b[i]) {
347 :			fit_b[i] = fit[i];
348 :			// copy contents of pos[i] to pos_b[i]
349 :			//memmove((void *)&pos_b[i], (void *)&pos[i],
350 :			// sizeof(double) * nvars);
351 :			pos_b[i]=pos[i];
352 :			}
353 :			// update gbest??
354 :			handle.logMessage(LOGDEBUG, "Particle ",i);
355 :			handle.logMessage(LOGDEBUG, "The likelihood score is", bestf, "at the point");
356 :			handle.logMessage(LOGDEBUG, "Its likelihood score is", fit[i], "at the point");
357 :			//#pragma omp critical (best2)
358 :			if (fit[i] < bestf) {
359 :			improved = 1;
360 :			// update best fitness
361 :			best.index=i;
362 :			bestf=fit[i];
363 :			// copy particle pos to gbest vector
364 :			//memmove((void *)&bestx, (void *)&pos[i],
365 :			// sizeof(double) * nvars);
366 :			// for (d = 0; d < nvars; d++)
367 :			// bestx[d] = pos[i][d] * init[d];
368 :			// iters = EcoSystem->getFuncEval() - offset;
369 :			// handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations");
370 :			// handle.logMessage(LOGINFO, "The likelihood score is", best.bestf, "at the point");
371 :			// EcoSystem->storeVariables(best.bestf, bestx);
372 :			// EcoSystem->writeBestValues();
373 :			// bestx=pos[i];
374 :
375 :			}
376 :			} // end i
377 :			if (best.bestf < bestf){
378 :			bestf=best.bestf;
379 :			bestx=pos[best.index];
380 :			for (d = 0; d < nvars; d++)
381 :			bestx[d] = bestx[d] * init[d];
382 :			iters = EcoSystem->getFuncEval() - offset; //To change. It is incorrect when is parallel.
383 :			handle.logMessage(LOGINFO, "\nNew optimum found after", iters, "function evaluations");
384 :			handle.logMessage(LOGINFO, "The likelihood score is", bestf, "at the point");
385 :			EcoSystem->storeVariables(bestf, bestx);
386 :			EcoSystem->writeBestValues();
387 :			bestx=pos[best.index];
388 :			}
389 :
390 :			psosteps++;
391 :
392 :			} // while (newf < bestf)
393 :
394 :			iters = EcoSystem->getFuncEval() - offset;
395 :			handle.logMessage(LOGINFO, "Existing PSO after ", iters, "function evaluations ...");
396 :			for (d = 0; d < nvars; d++)
397 :			bestx[d] = bestx[d] * init[d];
398 :			EcoSystem->storeVariables(bestf, bestx);
399 :			EcoSystem->writeBestValues();
400 :
401 :
402 :			}
403 :
404 :			void OptInfoPso::OptimiseLikelihoodREP() {
405 :			OptimiseLikelihood();
406 :
407 :			}
408 :
409 :			void OptInfoPso::OptimiseLikelihoodOMP() {
410 :			OptimiseLikelihood();
411 :			}
412 :
413 :			/**
414 :			* \brief calculate constant inertia weight equal to w_max
415 :			*/
416 :			double OptInfoPso::calc_inertia_const(int step) {
417 :			return w_max;
418 :			}
419 :			/**
420 :			* \brief calculate linearly decreasing inertia weight
421 :			*/
422 :			double OptInfoPso::calc_inertia_lin_dec(int step) {
423 :
424 :			int dec_stage = 3 * psoiter / 4;
425 :			if (step <= dec_stage)
426 :			return w_min + (w_max - w_min) * \
427 :			(dec_stage - step) / dec_stage;
428 :			else
429 :			return w_min;
430 :			}
431 :
432 :
433 :
434 :			/**
435 :			* \brief NEIGHBORHOOD (COMM) MATRIX STRATEGIES. global neighborhood
436 :			*/
437 :			void OptInfoPso::inform_global(IntMatrix& comm, DoubleMatrix& pos_nb,
438 :			DoubleMatrix& pos_b, DoubleVector& fit_b,
439 :			DoubleVector& gbest,int improved)
440 :			{
441 :
442 :			int i;
443 :			// all particles have the same attractor (gbest)
444 :			// copy the contents of gbest to pos_nb
445 :			for (i=0; i<size; i++)
446 :			//memmove((void *)&pos_nb[i][0], &gbest[0],
447 :			// sizeof(double) * gbest.Size());
448 :			pos_nb[i]=gbest;
449 :			}
450 :
451 :
452 :			/**
453 :			* \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
454 :			*/
455 :			void OptInfoPso::inform(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, int improved)
456 :			{
457 :			int i, j;
458 :			int b_n; // best neighbor in terms of fitness
459 :
460 :			// for each particle
461 :			for (j=0; j<size; j++) {
462 :			b_n = j; // self is best
463 :			// who is the best informer??
464 :			for (i=0; i<size; i++)
465 :			// the i^th particle informs the j^th particle
466 :			if (comm[i][j] && fit_b[i] < fit_b[b_n])
467 :			// found a better informer for j^th particle
468 :			b_n = i;
469 :			// copy pos_b of b_n^th particle to pos_nb[j]
470 :			//memmove((void *)&pos_nb[j][0],
471 :			// (void *)&pos_b[b_n][0],
472 :			// sizeof(double) * pos_nb.Ncol(j));
473 :			pos_nb[j]=pos_b[b_n];
474 :			}
475 :			}
476 :
477 :			/**
478 :			* \brief topology initialization :: this is a static (i.e. fixed) topology
479 :			*/
480 :			void OptInfoPso::init_comm_ring(IntMatrix& comm) {
481 :			int i;
482 :			// reset array
483 :			comm.setToZero();
484 :
485 :			// choose informers
486 :			for (i=0; i<size; i++) {
487 :			// set diagonal to 1
488 :			comm[i][i] = 1;
489 :			if (i==0) {
490 :			// look right
491 :			comm[i][1] = 1;
492 :			// look left
493 :			comm[i][size-1] = 1;
494 :			} else if (i==size-1) {
495 :			// look right
496 :			comm[i][0] = 1;
497 :			// look left
498 :			comm[i][i-1] = 1;
499 :			} else {
500 :			// look right
501 :			comm[i][i+1] = 1;
502 :			// look left
503 :			comm[i][i-1] = 1;
504 :			}
505 :
506 :			}
507 :			// Print Matrix
508 :			handle.logMessage(LOGDEBUG, "\nPSO ring communication\n");
509 :			for (i=0;i<size;i++){
510 :			string message="[";
511 :			for (int j=0;j<size;j++)
512 :			{
513 :			std::ostringstream ss;
514 :			ss<< comm[i][j];
515 :			message=message+" "+ss.str();
516 :			}
517 :			message=message+"]\n";
518 :			handle.logMessage(LOGDEBUG, message.c_str());
519 :			}
520 :
521 :			}
522 :
523 :
524 :
525 :
526 :			void OptInfoPso::inform_ring(IntMatrix& comm, DoubleMatrix& pos_nb,
527 :			DoubleMatrix& pos_b, DoubleVector& fit_b,
528 :			DoubleVector& gbest, int improved)
529 :			{
530 :
531 :			// update pos_nb matrix
532 :			inform(comm, pos_nb, pos_b, fit_b, improved);
533 :
534 :			}
535 :
536 :			/**
537 :			* \brief random neighborhood topology
538 :			*/
539 :			void OptInfoPso::init_comm_random(IntMatrix& comm) {
540 :
541 :			int i, j, k;
542 :			// reset array
543 :			comm.setToZero();
544 :
545 :			// choose informers
546 :			for (i=0; i<size; i++) {
547 :			// each particle informs itself
548 :			comm[i][i] = 1;
549 :			// choose kappa (on average) informers for each particle
550 :			for (k=0; k<nhood_size; k++) {
551 :			// generate a random index
552 :			//j = gsl_rng_uniform_int(settings->rng, settings->size);
553 :			j=rand()*size/RAND_MAX;
554 :			// particle i informs particle j
555 :			comm[i] [j] = 1;
556 :			}
557 :			}
558 :			}
559 :
560 :
561 :
562 :			void OptInfoPso::inform_random(IntMatrix& comm, DoubleMatrix& pos_nb,
563 :			DoubleMatrix& pos_b, DoubleVector& fit_b,
564 :			DoubleVector& gbest, int improved)
565 :			{
566 :
567 :
568 :			// regenerate connectivity??
569 :			if (!improved)
570 :			init_comm_random(comm);
571 :			inform(comm, pos_nb, pos_b, fit_b, improved);
572 :
573 :			}
574 :

---

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