PaGMO  1.1.5
pso_generational_racing.cpp
1 /*****************************************************************************
2  * Copyright (C) 2004-2015 The PaGMO development team, *
3  * Advanced Concepts Team (ACT), European Space Agency (ESA) *
4  * *
5  * https://github.com/esa/pagmo *
6  * *
7  * act@esa.int *
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  * This program is distributed in the hope that it will be useful, *
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
17  * GNU General Public License for more details. *
18  * *
19  * You should have received a copy of the GNU General Public License *
20  * along with this program; if not, write to the *
21  * Free Software Foundation, Inc., *
22  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
23  *****************************************************************************/
24 
25 #include <boost/random/uniform_int.hpp>
26 #include <boost/random/uniform_real.hpp>
27 #include <vector>
28 #include <cmath>
29 #include <iostream>
30 
31 #include "pso_generational_racing.h"
32 #include "../problem/base_stochastic.h"
33 
34 namespace pagmo { namespace algorithm {
35 
36 using namespace util::racing;
37 
39 
58 pso_generational_racing::pso_generational_racing(int gen, double omega, double eta1, double eta2, double vcoeff, int variant, int neighb_type, int neighb_param, unsigned int nr_eval_per_x, unsigned int max_fevals): base(), m_gen(gen), m_omega(omega), m_eta1(eta1), m_eta2(eta2), m_vcoeff(vcoeff), m_variant(variant), m_neighb_type(neighb_type), m_neighb_param(neighb_param), m_nr_eval_per_x(nr_eval_per_x), m_fevals(0), m_max_fevals(max_fevals) {
59  if (gen < 0) {
60  pagmo_throw(value_error,"number of generations must be nonnegative");
61  }
62 
63  if (m_omega < 0 || m_omega > 1) {
64  if( variant < 5 )
65  // variants using Inertia weight
66  pagmo_throw(value_error,"the particles' inertia must be in the [0,1] range");
67  else
68  // variants using Constriction coefficient
69  pagmo_throw(value_error,"the particles' constriction coefficient must be in the [0,1] range");
70  }
71 
72  if (eta1 < 0 || eta2 < 0 || eta1 > 4 || eta2 > 4) {
73  pagmo_throw(value_error,"the eta parameters must be in the [0,4] range");
74  }
75 
76  if (vcoeff <= 0 || vcoeff > 1) {
77  pagmo_throw(value_error,"fraction of variables' range in which velocity may vary should be in ]0,1]");
78  }
79 
80  if (variant < 1 || variant > 6) {
81  pagmo_throw(value_error,"algorithm variant must be one of 1 ... 6");
82  }
83 
84  if (neighb_type < 1 || neighb_type > 4) {
85  pagmo_throw(value_error,"swarm topology variant must be one of 1 ... 4");
86  }
87  // And neighb param???
88 }
89 
90 
93 {
94  return base_ptr(new pso_generational_racing(*this));
95 }
96 
97 // Construct the racing structure from lists of decision vectors tailored to
98 // this pso_gen implementation. If only a single list of decision vectors is
99 // required, simply leave the second list empty.
100 void pso_generational_racing::racing__construct_race_environment( util::racing::race_pop & race_structure, const problem::base& prob, const std::vector<decision_vector> &x_list1, const std::vector<decision_vector> &x_list2 ) const
101 {
102  util::racing::racing_population lbpop( prob );
103  for( unsigned int p = 0; p < x_list1.size(); p++ ){
104  lbpop.push_back_noeval( x_list1[p] );
105  }
106  for( unsigned int p = 0; p < x_list2.size(); p++ ){
107  lbpop.push_back_noeval( x_list2[p] );
108  }
109  race_structure.register_population( lbpop );
110 }
111 
112 // Run several kinds of racing encountered in pso_gen:
113 //
114 // 1) Between two individuals: Their indices should be set as idx1 and idx2
115 // respectively.
116 // 2) Among the whole population: In this case set idx1 and idx2 both to -1.
117 //
118 // A tuple containing the index of the winner, and the number of fitness
119 // evaluation consumed will be returned.
120 std::pair<population::size_type, unsigned int> pso_generational_racing::racing__race_for_winner( util::racing::race_pop &race_structure, int idx1, int idx2, unsigned int max_fevals) const
121 {
122  // Check if the provided indices are valid
123  if(!(idx1<0 && idx2<0) && !(idx1>=0 && idx2>=0)){
124  pagmo_throw(value_error, "The pair must either be both -1 (race all) or both explicitly specified.");
125  }
126 
127  // If this is true, the max_fevals need to be scaled down accordingly, as
128  // it has a different meaning when being passed in.
129  bool use_data_count_termination = false;
130 
131  std::vector<population::size_type> active_set;
132  // Run race on all the individuals
133  if(idx1 < 0){
134  active_set.resize(race_structure.size());
135  for(population::size_type i = 0; i < race_structure.size(); i++){
136  active_set[i] = i;
137  }
138  if(use_data_count_termination){
139  max_fevals /= race_structure.size();
140  }
141  }
142  // Race between the two specified individuals
143  else{
144  active_set.resize(2);
145  active_set[0] = idx1;
146  active_set[1] = idx2;
147  if(use_data_count_termination){
148  max_fevals /= 2;
149  }
150  }
151  std::pair<std::vector<population::size_type>, unsigned int> res;
152  if(use_data_count_termination){
153  res = race_structure.run(1, 0, max_fevals, 0.001, active_set, race_pop::MAX_DATA_COUNT, true, false);
154  }
155  else{
156  res = race_structure.run(1, 0, max_fevals, 0.001, active_set, race_pop::MAX_BUDGET, true, false);
157  }
158  return std::make_pair(res.first[0], res.second);
159 }
160 
161 // Compute the averaged fitness value of a stochastic problem
162 void pso_generational_racing::particle__average_fitness(const problem::base &prob, fitness_vector &f, const decision_vector &x) const
163 {
164  f = fitness_vector(prob.get_f_dimension(), 0);
165  for(unsigned int i = 0; i < m_nr_eval_per_x; i++){
166  dynamic_cast<const pagmo::problem::base_stochastic &>(prob).set_seed(m_urng());
167  fitness_vector f_tmp = prob.objfun(x);
168  for(unsigned int j = 0; j < f.size(); j++){
169  f[j] += f_tmp[j] / m_nr_eval_per_x;
170  }
171  }
172 }
173 
174 // Compute the averaged fitness value of a decision vector in a stochastic problem,
175 void pso_generational_racing::particle__average_fitness_set_best(const problem::base &prob, std::vector<fitness_vector> &F, population::size_type& best_idx, fitness_vector& best_fit, const std::vector<decision_vector> &X) const
176 {
177  F = std::vector<fitness_vector>(X.size(), fitness_vector(prob.get_f_dimension()));
178  for(unsigned int i = 0; i < X.size(); i++){
179  particle__average_fitness(prob, F[i], X[i]);
180  if(i==0){
181  best_idx = i;
182  best_fit = F[i];
183  }
184  else{
185  if(prob.compare_fitness(F[i], best_fit)){
186  best_idx = i;
187  best_fit = F[i];
188  }
189  }
190  }
191 }
192 
194 
216 {
217  // Let's store some useful variables.
218  const problem::base &prob = pop.problem();
219  const problem::base::size_type D = prob.get_dimension(), prob_i_dimension = prob.get_i_dimension(), prob_c_dimension = prob.get_c_dimension(), prob_f_dimension = prob.get_f_dimension();
220  const problem::base::size_type Dc = D - prob_i_dimension;
221  const decision_vector &lb = prob.get_lb(), &ub = prob.get_ub();
222  const population::size_type swarm_size = pop.size();
223 
224 
225  //We perform some checks to determine wether the problem/population are suitable for PSO
226  if( Dc == 0 ){
227  pagmo_throw(value_error,"There is no continuous part in the problem decision vector for PSO to optimise");
228  }
229 
230  if( prob_c_dimension != 0 ){
231  pagmo_throw(value_error,"The problem is not box constrained and PSO is not suitable to solve it");
232  }
233 
234  if( prob_f_dimension != 1 ){
235  pagmo_throw(value_error,"The problem is not single objective and PSO is not suitable to solve it");
236  }
237 
238  // Get out if there is nothing to do.
239  if (swarm_size == 0 || m_gen == 0) {
240  return;
241  }
242 
243  // If problem is not stochastic, no point to use a racing based PSO
244  try{
245  dynamic_cast<const pagmo::problem::base_stochastic &>(prob).set_seed(m_urng());
246  }
247  catch (const std::bad_cast& e){
248  // TODO: Should throw, but current algorithm's serialization test does not supply stochastic problem
249  // pagmo_throw(value_error, "The problem is not stochastic and racing based PSO is not suitable to solve it");
250  std::cout << "The problem is not stochastic and racing based PSO is not suitable to solve it" << std::endl;
251  return;
252  }
253 
254  // Some vectors used during evolution are allocated here.
255  std::vector<double> dummy(Dc,0); // used for initialisation purposes
256 
257  std::vector<decision_vector> X(swarm_size,dummy); // particles' current positions
258  std::vector<fitness_vector> fit(swarm_size); // particles' current fitness values
259 
260  std::vector<decision_vector > V(swarm_size,dummy); // particles' velocities
261 
262  std::vector<decision_vector> lbX(swarm_size,dummy); // particles' previous best positions
263  std::vector<fitness_vector> lbfit(swarm_size); // particles' fitness values at their previous best positions
264 
265 
266  std::vector< std::vector<int> > neighb(swarm_size); // swarm topology (iterators over indexes of each particle's neighbors in the swarm)
267 
268  decision_vector best_neighb(Dc); // search space position of particles' best neighbor
269  fitness_vector best_fit; // fitness at the best found search space position (tracked only when using topologies 1 or 4)
270  population::size_type best_idx = 0; // index into the best solution in the swarm (best_fit == lbfit[best_idx]) (as with best_fit, this is only tracked when using topologies 1 or 4)
271  bool best_fit_improved; // flag indicating whether the best solution's fitness improved (tracked only when using topologies 1 or 4)
272 
273  decision_vector minv(Dc), maxv(Dc); // Maximum and minumum velocity allowed
274 
275  double vwidth; // Temporary variable
276  double new_x; // Temporary variable
277 
278  population::size_type p; // for iterating over particles
279  population::size_type n; // for iterating over particles's neighbours
280  problem::base::size_type d; // for iterating over problem dimensions
281 
282 
283  // Initialise the minimum and maximum velocity
284  for( d = 0; d < Dc; d++ ){
285  vwidth = ( ub[d] - lb[d] ) * m_vcoeff;
286  minv[d] = -1.0 * vwidth;
287  maxv[d] = vwidth;
288  }
289 
290  // Copy the particle positions, their velocities and their fitness
291  for( p = 0; p < swarm_size; p++ ){
292  X[p] = pop.get_individual(p).cur_x;
293  V[p] = pop.get_individual(p).cur_v;
294  fit[p] = pop.get_individual(p).cur_f;
295  }
296 
297  // Initialise particles' previous best positions
298  for( p = 0; p < swarm_size; p++ ){
299  lbX[p] = pop.get_individual(p).best_x;
300  lbfit[p] = pop.get_individual(p).best_f;
301  }
302 
303  // Initialize the Swarm's topology
304  switch( m_neighb_type ){
305  case 1: initialize_topology__gbest( pop, best_neighb, best_fit, neighb ); break;
306  case 3: initialize_topology__von( neighb ); break;
307  case 4: initialize_topology__adaptive_random( neighb );
308  best_idx = pop.get_best_idx();
309  best_fit = lbfit[ best_idx ]; // need to track improvements in best found fitness, to know when to rewire
310  break;
311  case 2:
312  default: initialize_topology__lbest( neighb );
313  }
314 
315  m_fevals = 0;
316 
317  // Initialize the seed to be used to construct different race_pop instances
318  // to be used to do racing in different contexts: Within the new X, or
319  // within X + lbX, so that past evaluation data can be reused.
320  unsigned int racing_seed = m_urng();
321  util::racing::race_pop race_lbX(racing_seed);
322  util::racing::race_pop race_lbX_and_X(racing_seed);
323 
324  racing__construct_race_environment(race_lbX, pop.problem(), lbX, std::vector<decision_vector>());
325 
326  if( m_neighb_type == 1 ){
327 
328  // Using the gbest (fully connected) topology, race all the individuals
329  std::pair<population::size_type, unsigned int> res
330  = racing__race_for_winner(race_lbX, -1, -1, pop.size() * m_nr_eval_per_x);
331  best_idx = res.first;
332  m_fevals += res.second;
333 
334  // Set lbfit to be the averaged fitness values from racing
335  lbfit = race_lbX.get_mean_fitness();
336  best_neighb = lbX[ best_idx ];
337  best_fit = lbfit[ best_idx ];
338  }
339 
340  // auxiliary varibables specific to the Fully Informed Particle Swarm variant
341  double acceleration_coefficient = m_eta1 + m_eta2;
342  double sum_forces;
343 
344  double r1 = 0.0;
345  double r2 = 0.0;
346 
347  bool forced_terminate = false;
348 
349  // Invoke racing in a neighbourhood topology to determine the best
350  // neighbour. Otherwise, use the averaged fitness values.
351  bool racing_opt__race_neighbourhood = true;
352 
353  /* --- Main PSO loop ---
354  */
355  // For each generation
356  int g = 0;
357  while( g < m_gen && m_fevals < m_max_fevals ){
358  g++;
359 
360  // Initialize a new list of internal seeds for use in racing
361  unsigned int cur_racing_seed = m_urng();
362  race_lbX.set_seed(cur_racing_seed);
363  race_lbX_and_X.set_seed(cur_racing_seed);
364 
365  // Update Velocity
366  for( p = 0; p < swarm_size; p++ ){
367 
368  // identify the current particle's best neighbour. not needed if
369  // m_neighb_type == 1 (gbest): best_neighb directly tracked in this
370  // function . not needed if m_variant == 6 (FIPS): all neighbours
371  // are considered, no need to identify the best one
372  if( m_neighb_type != 1 && m_variant != 6){
373  if( !racing_opt__race_neighbourhood ){
374  best_neighb = particle__get_best_neighbor( p, neighb, lbX, lbfit, prob );
375  }
376  else{
377  best_neighb = particle__racing_get_best_neighbor( p, neighb, lbX, race_lbX );
378  // If the swarm has not been completely processed but
379  // racing has exhausted the allowable evaluation budget, we
380  // have to terminate. The final feval count is capped at
381  // its maximum. This is still fair as the information
382  // gained from the exceeded fevals will not used at all
383  // within the evolution.
384  if(m_fevals > m_max_fevals){
385  m_fevals = m_max_fevals;
386  forced_terminate = true;
387  break;
388  }
389  }
390  }
391  /*-------PSO canonical (with inertia weight) ---------------------------------------------*/
392  /*-------Original algorithm used in PaGMO paper-------------------------------------------*/
393  if( m_variant == 1 ){
394  for( d = 0; d < Dc; d++ ){
395  r1 = m_drng();
396  r2 = m_drng();
397  V[p][d] = m_omega * V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) + m_eta2 * r2 * (best_neighb[d] - X[p][d]);
398  }
399  }
400 
401  /*-------PSO canonical (with inertia weight) ---------------------------------------------*/
402  /*-------and with equal random weights of social and cognitive components-----------------*/
403  /*-------Check with Rastrigin-------------------------------------------------------------*/
404  else if( m_variant == 2 ){
405  for( d = 0; d < Dc; d++ ){
406  r1 = m_drng();
407  V[p][d] = m_omega * V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) + m_eta2 * r1 * (best_neighb[d] - X[p][d]);
408  }
409  }
410 
411  /*-------PSO variant (commonly mistaken in literature for the canonical)----------------*/
412  /*-------Same random number for all components------------------------------------------*/
413  else if( m_variant == 3 ){
414  r1 = m_drng();
415  r2 = m_drng();
416  for( d = 0; d < Dc; d++ ){
417  V[p][d] = m_omega * V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) + m_eta2 * r2 * (best_neighb[d] - X[p][d]);
418  }
419  }
420 
421  /*-------PSO variant (commonly mistaken in literature for the canonical)----------------*/
422  /*-------Same random number for all components------------------------------------------*/
423  /*-------and with equal random weights of social and cognitive components---------------*/
424  else if( m_variant == 4 ){
425  r1 = m_drng();
426  for( d = 0; d < Dc; d++ ){
427  V[p][d] = m_omega * V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) + m_eta2 * r1 * (best_neighb[d] - X[p][d]);
428  }
429  }
430 
431  /*-------PSO variant with constriction coefficients------------------------------------*/
432  /* ''Clerc's analysis of the iterative system led him to propose a strategy for the
433  * placement of "constriction coefficients" on the terms of the formulas; these
434  * coefficients controlled the convergence of the particle and allowed an elegant and
435  * well-explained method for preventing explosion, ensuring convergence, and
436  * eliminating the arbitrary Vmax parameter. The analysis also takes the guesswork
437  * out of setting the values of phi_1 and phi_2.''
438  * ''this is the canonical particle swarm algorithm of today.''
439  * [Poli et al., 2007] http://dx.doi.org/10.1007/s11721-007-0002-0
440  * [Clerc and Kennedy, 2002] http://dx.doi.org/10.1109/4235.985692
441  *
442  * This being the canonical PSO of today, this variant is set as the default in PaGMO.
443  *-------------------------------------------------------------------------------------*/
444  else if( m_variant == 5 ){
445  for( d = 0; d < Dc; d++ ){
446  r1 = m_drng();
447  r2 = m_drng();
448  V[p][d] = m_omega * ( V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) + m_eta2 * r2 * (best_neighb[d] - X[p][d]) );
449  }
450  }
451 
452  /*-------Fully Informed Particle Swarm-------------------------------------------------*/
453  /* ''Whereas in the traditional algorithm each particle is affected by its own
454  * previous performance and the single best success found in its neighborhood, in
455  * Mendes' fully informed particle swarm (FIPS), the particle is affected by all its
456  * neighbors, sometimes with no influence from its own previous success.''
457  * ''With good parameters, FIPS appears to find better solutions in fewer iterations
458  * than the canonical algorithm, but it is much more dependent on the population topology.''
459  * [Poli et al., 2007] http://dx.doi.org/10.1007/s11721-007-0002-0
460  * [Mendes et al., 2004] http://dx.doi.org/10.1109/TEVC.2004.826074
461  *-------------------------------------------------------------------------------------*/
462  else if( m_variant == 6 ){
463  for( d = 0; d < Dc; d++ ){
464  sum_forces = 0.0;
465  for( n = 0; n < neighb[p].size(); n++ )
466  sum_forces += m_drng() * acceleration_coefficient * ( lbX[ neighb[p][n] ][d] - X[p][d] );
467 
468  V[p][d] = m_omega * ( V[p][d] + sum_forces / neighb[p].size() );
469  }
470  }
471  }
472 
473  if(forced_terminate){
474  break;
475  }
476 
477  // Update Position
478  for( p = 0; p < swarm_size; p++ ){
479  // We now check that the velocity does not exceed the maximum allowed per component
480  // and we perform the position update and the feasibility correction
481  for( d = 0; d < Dc; d++ ){
482 
483  if( V[p][d] > maxv[d] )
484  V[p][d] = maxv[d];
485 
486  else if( V[p][d] < minv[d] )
487  V[p][d] = minv[d];
488 
489  // update position
490  new_x = X[p][d] + V[p][d];
491 
492  // feasibility correction
493  // (velocity updated to that which would have taken the previous position
494  // to the newly corrected feasible position)
495  if( new_x < lb[d] ){
496  new_x = lb[d];
497  V[p][d] = 0.0;
498 // new_x = boost::uniform_real<double>(lb[d],ub[d])(m_drng);
499 // V[p][d] = new_x - X[p][d];
500 // V[p][d] = 0;
501  }
502  else if( new_x > ub[d] ){
503  new_x = ub[d];
504  V[p][d] = 0.0;
505 // new_x = boost::uniform_real<double>(lb[d],ub[d])(m_drng);
506 // V[p][d] = new_x - X[p][d];
507 // V[p][d] = 0;
508  }
509  X[p][d] = new_x;
510  }
511  }
512  std::vector<bool> local_bests_improved = std::vector<bool>(lbX.size(), false);
513 
514  // Problem is stochastic, change the seed
515  dynamic_cast<const pagmo::problem::base_stochastic &>(prob).set_seed(m_urng());
516 
517  // Stochastic problem being solved with the assistance of racing
518 
519  racing__construct_race_environment(race_lbX_and_X, pop.problem(), lbX, X);
520  race_lbX_and_X.inherit_memory(race_lbX);
521 
522  for( p = 0; p < swarm_size && !forced_terminate; p++ ){
523  std::pair<population::size_type, unsigned int> res =
524  racing__race_for_winner(race_lbX_and_X, p, p + swarm_size, 2 * m_nr_eval_per_x);
525  m_fevals += res.second;
526  if (m_fevals > m_max_fevals){
527  forced_terminate = true;
528  break;
529  }
530  if(res.first == p+swarm_size){
531  local_bests_improved[p] = true;
532  }
533  }
534 
535  // If we run out of budget in the previous loop, the following
536  // operations are not meaningful as they may be out-of-date
537  if(forced_terminate){
538  break;
539  }
540 
541  pop.clear();
542 
543  // Update lbfit and fit to hold the averaged fitness values
544  std::vector<fitness_vector> averaged_fitness =
545  race_lbX_and_X.get_mean_fitness();
546  for(unsigned int i = 0; i < swarm_size; i++){
547  lbfit[i] = averaged_fitness[i];
548  fit[i] = averaged_fitness[i+swarm_size];
549  }
550 
551  // We update the particles memory if a better point has been reached
552  best_fit_improved = false;
553  // The decisions here (whether to update local best or not) are
554  // based on the results of racing performed previously.
555  for( p = 0; p < swarm_size; p++ ){
556  // Use racing to update local bests
557  if(local_bests_improved[p]){
558  // current position better than previous best update the
559  // particle's previous best position
560  lbfit[p] = fit[p];
561  lbX[p] = X[p];
562 
563  // case of improvement not detected below when
564  // m_neighb_type == 4
565  if( p == best_idx )
566  best_fit_improved = true;
567  }
568  pop.push_back(lbX[p]);
569  pop.set_v(p,V[p]);
570 
571  // Skipping set_x of X[p] intentionally, so that this information
572  // is not lost when m_gen is 1
573  // pop.set_x(p,X[p]);
574  }
575 
576  // Construct the new race environment with the updated lbX. Pass on the
577  // memory to the next generation -- some new memory may come from the
578  // previous race between lbX and X. NOTE: This is only possible if the
579  // ground seed of the race doesn't change in the next iteration,
580  // otherwise the transferred cache will be cleared anyway.
581  racing__construct_race_environment(race_lbX, pop.problem(), lbX, std::vector<decision_vector>());
582  race_lbX.inherit_memory(race_lbX_and_X);
583 
584  // update the best position observed so far by any particle in the swarm
585  // (only performed if swarm topology is gbest or random varying)
586  if( m_neighb_type == 1 ){
587  // Race to get the best in lbX
588  std::pair<population::size_type, unsigned int> res =
589  racing__race_for_winner(race_lbX, -1, -1, swarm_size * m_nr_eval_per_x);
590 
591  // Set fitness to the averaged fitness from race
592  lbfit = race_lbX.get_mean_fitness();
593 
594  m_fevals += res.second;
595 
596  p = res.first;
597  if( p != best_idx ){
598  best_idx = p;
599  best_fit_improved = true;
600  }
601  best_neighb = lbX[ best_idx ];
602  best_fit = lbfit[ best_idx ];
603  }
604  if( m_neighb_type == 4 ){
605  // Improvement -> A decision vector other than the one
606  // previously identified as being the top solution now has the
607  // best fitness
608  best_neighb = lbX[ best_idx ];
609  best_fit = lbfit[ best_idx ];
610  for( p = 0; p < swarm_size; p++ ){
611  if( prob.compare_fitness( lbfit[p], best_fit ) ){
612  best_idx = p;
613  best_neighb = lbX[p];
614  best_fit = lbfit[p];
615  best_fit_improved = true;
616  }
617  }
618  }
619 
620  // reset swarm topology if no improvement was observed in the best
621  // found fitness value
622  if( m_neighb_type == 4 && !best_fit_improved )
623  {
624  initialize_topology__adaptive_random( neighb );
625  }
626 
627  } // end of main PSO loop
628  std::cout << "PSO terminated: gen = " << g << ", incurred fevals = " << m_fevals << std::endl;
629 }
630 
631 
642 decision_vector pso_generational_racing::particle__get_best_neighbor( population::size_type pidx, std::vector< std::vector<int> > &neighb, const std::vector<decision_vector> &lbX, const std::vector<fitness_vector> &lbfit, const problem::base &prob) const
643 {
644  population::size_type nidx, bnidx; // neighbour index; best neighbour index
645 
646  switch( m_neighb_type ){
647  case 1: // { gbest }
648  // ERROR: execution should not reach this point, as the global best position is not tracked using the neighb vector
649  pagmo_throw(value_error,"particle__get_best_neighbor() invoked while using a gbest swarm topology");
650  break;
651  case 2: // { lbest }
652  case 3: // { von }
653  case 4: // { adaptive random }
654  default:
655  // iterate over indexes of the particle's neighbours, and identify the best
656  bnidx = neighb[pidx][0];
657  for( nidx = 1; nidx < neighb[pidx].size(); nidx++ )
658  if( prob.compare_fitness( lbfit[ neighb[pidx][nidx] ], lbfit[ bnidx ] ) )
659  bnidx = neighb[pidx][nidx];
660  return lbX[bnidx];
661  }
662 }
663 
665 decision_vector pso_generational_racing::particle__racing_get_best_neighbor( population::size_type pidx, std::vector< std::vector<int> > &neighb, const std::vector<decision_vector> &lbX, util::racing::race_pop& race) const
666 {
667  std::vector<population::size_type> active_indices;
668  for(population::size_type nidx = 0; nidx < neighb[pidx].size(); nidx++){
669  bool repeated = false;
670  for(unsigned int i = 0; i < active_indices.size(); i++){
671  // During the intialization of random adaptive neighbourhood no
672  // check is performed to avoid repeated indices
673  if((unsigned int)neighb[pidx][nidx] == active_indices[i]){
674  repeated = true;
675  break;
676  }
677  }
678  if(!repeated)
679  active_indices.push_back(neighb[pidx][nidx]);
680  }
681 
682  std::pair<std::vector<population::size_type>, unsigned int> race_res = race.run(2, 0, neighb[pidx].size() * m_nr_eval_per_x, 0.01, active_indices, race_pop::MAX_BUDGET, true, false);
683  //std::pair<std::vector<population::size_type>, unsigned int> race_res = race.run(1, 0, m_nr_eval_per_x, 0.01, active_indices, race_pop::MAX_DATA_COUNT, true, false);
684  std::vector<population::size_type> winners = race_res.first;
685  m_fevals += race_res.second;
686  return lbX[winners[0]];
687 }
688 
708 void pso_generational_racing::initialize_topology__gbest( const population &pop, decision_vector &gbX, fitness_vector &gbfit, std::vector< std::vector<int> > &neighb ) const
709 {
710  // The best position already visited by the swarm will be tracked in pso::evolve() as particles are evaluated.
711  // Here we define the initial values of the variables that will do that tracking.
712  gbX = pop.champion().x;
713  gbfit = pop.champion().f;
714 
715  /* The usage of a gbest swarm topology along with a FIPS (fully informed particle swarm) velocity update formula
716  * is discouraged. However, because a user might still configure such a setup, we must ensure FIPS has access to
717  * the list of indices of particles' neighbours:
718  */
719  if( m_variant == 6 ){
720  unsigned int i;
721  for( i = 0; i < neighb.size(); i++ )
722  neighb[0].push_back( i );
723  for( i = 1; i < neighb.size(); i++ )
724  neighb[i] = neighb[0];
725  }
726 }
727 
728 
741 void pso_generational_racing::initialize_topology__lbest( std::vector< std::vector<int> > &neighb ) const
742 {
743  int swarm_size = neighb.size();
744  int pidx; // for iterating over particles
745  int nidx, j; // for iterating over particles' neighbours
746 
747  int radius = m_neighb_param / 2;
748 
749  for( pidx = 0; pidx < swarm_size; pidx++ ){
750  for( j = -radius; j <= radius; j++ ){
751  if( j == 0 ) j++;
752  nidx = (pidx + j) % swarm_size;
753  if( nidx < 0 ) nidx = swarm_size + nidx;
754  neighb[pidx].push_back( nidx );
755  }
756  }
757 }
758 
759 
769 const int vonNeumann_neighb_diff[4][2] = { {-1,0}, {1,0}, {0,-1}, {0,1} };
770 
788 void pso_generational_racing::initialize_topology__von( std::vector< std::vector<int> > &neighb ) const
789 {
790  int swarm_size = neighb.size();
791  int cols, rows; // lattice structure
792  int pidx, nidx; // particle and neighbour indices, in the swarm and neighbourhood vectors
793  int p_x, p_y; // particle's coordinates in the lattice
794  int n_x, n_y; // particle neighbor's coordinates in the lattice
795 
796  rows = std::sqrt( swarm_size );
797  while( swarm_size % rows != 0 )
798  rows -= 1;
799  cols = swarm_size / rows;
800 
801  for( pidx = 0; pidx < swarm_size; pidx++ ){
802  p_x = pidx % cols;
803  p_y = pidx / cols;
804 
805  for( nidx = 0; nidx < 4; nidx++ ){
806  n_x = ( p_x + vonNeumann_neighb_diff[nidx][0] ) % cols; if( n_x < 0 ) n_x = cols + n_x; // sign of remainder(%) in a division when at least one of the operands is negative is compiler implementation specific. The 'if' here ensures the same behaviour across compilers
807  n_y = ( p_y + vonNeumann_neighb_diff[nidx][1] ) % rows; if( n_y < 0 ) n_y = rows + n_y;
808 
809  neighb[pidx].push_back( n_y * cols + n_x );
810  }
811  }
812 }
813 
814 
837 void pso_generational_racing::initialize_topology__adaptive_random( std::vector< std::vector<int> > &neighb ) const
838 {
839  int swarm_size = neighb.size();
840  int pidx; // for iterating over particles
841  int nidx, j; // for iterating over particles being connected to
842 
843  // clear previously defined topology
844  for( pidx = 0; pidx < swarm_size; pidx++ )
845  neighb[pidx].clear();
846 
847  // define new topology
848  for( pidx = 0; pidx < swarm_size; pidx++ ){
849 
850  // the particle always connects back to itself, thus guaranteeing a minimum indegree of 1
851  neighb[pidx].push_back( pidx );
852 
853  for( j = 1; j < m_neighb_param; j++ ){
854  nidx = m_drng() * swarm_size;
855  neighb[nidx].push_back( pidx );
856  // No check performed to see whether pidx is already in neighb[nidx],
857  // leading to a performance penalty in particle__get_best_neighbor() when it occurs.
858  }
859  }
860 }
861 
862 
865 {
866  return "PSO - Generational with racing";
867 }
868 
869 
871 
875 {
876  std::ostringstream s;
877  s << "gen:" << m_gen << ' ';
878  s << "omega:" << m_omega << ' ';
879  s << "eta1:" << m_eta1 << ' ';
880  s << "eta2:" << m_eta2 << ' ';
881  s << "variant:" << m_variant << ' ';
882  s << "topology:" << m_neighb_type << ' ';
883  if( m_neighb_type == 2 || m_neighb_type == 4 )
884  s << "topology param.:" << m_neighb_param << ' ';
885  return s.str();
886 }
887 
888 }} //namespaces
889 
890 BOOST_CLASS_EXPORT_IMPLEMENT(pagmo::algorithm::pso_generational_racing)
boost::shared_ptr< base > base_ptr
Alias for shared pointer to base algorithm.
Root PaGMO namespace.
std::vector< double > decision_vector
Decision vector type.
Definition: types.h:40
Particle Swarm optimization generational with racing.
Fixed number of iterations.
Definition: race_pop.h:70
fitness_vector cur_f
Current fitness vector.
Definition: population.h:89
pso_generational_racing(int gen=1, double omega=0.7298, double eta1=2.05, double eta2=2.05, double vcoeff=0.5, int variant=5, int neighb_type=2, int neighb_param=4, unsigned int nr_eval_per_x=5, unsigned int max_fevals=std::numeric_limits< unsigned int >::max())
Constructor.
const individual_type & get_individual(const size_type &) const
Get constant reference to individual at position n.
Definition: population.cpp:277
void register_population(const population &)
Update the population on which the race will run.
Definition: race_pop.cpp:41
Base algorithm class.
const int vonNeumann_neighb_diff[4][2]
Von Neumann neighborhood (increments on particles' lattice coordinates that produce the coordinates o...
Definition: pso.cpp:474
Base problem class.
Definition: problem/base.h:148
Population class.
Definition: population.h:70
Fixed number of function evaluations.
Definition: race_pop.h:69
size_type get_dimension() const
Return global dimension.
bool compare_fitness(const fitness_vector &, const fitness_vector &) const
Compare fitness vectors.
std::string human_readable_extra() const
Extra human readable algorithm info.
void evolve(population &) const
Evolve implementation.
size_type get_i_dimension() const
Return integer dimension.
fitness_vector best_f
Best fitness vector so far.
Definition: population.h:95
std::vector< double > fitness_vector
Fitness vector type.
Definition: types.h:42
c_size_type get_c_dimension() const
Return global constraints dimension.
decision_vector cur_v
Current velocity vector.
Definition: population.h:85
const decision_vector & get_ub() const
Upper bounds getter.
Base Stochastic Optimization Problem.
container_type::size_type size_type
Population size type.
Definition: population.h:192
decision_vector cur_x
Current decision vector.
Definition: population.h:83
rng_uint32 m_urng
Random number generator for unsigned integer values.
f_size_type get_f_dimension() const
Return fitness dimension.
std::string get_name() const
Algorithm name.
const decision_vector & get_lb() const
Lower bounds getter.
rng_double m_drng
Random number generator for double-precision floating point values.
decision_vector best_x
Best decision vector so far.
Definition: population.h:91
decision_vector::size_type size_type
Problem's size type: the same as pagmo::decision_vector's size type.
Definition: problem/base.h:160
Racing mechanism for a population.
Definition: race_pop.h:60
void set_seed(unsigned int)
Set a new racing seed.
Definition: race_pop.cpp:695
std::vector< fitness_vector > get_mean_fitness(const std::vector< population::size_type > &active_set=std::vector< population::size_type >()) const
Returns mean fitness of the individuals based on past evaluation data.
Definition: race_pop.cpp:532
void inherit_memory(const race_pop &)
Inherits the memory of another race_pop object.
Definition: race_pop.cpp:653