PaGMO  1.1.5
sga.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 <boost/random/variate_generator.hpp>
28 #include <boost/random/normal_distribution.hpp>
29 #include <boost/math/special_functions/round.hpp>
30 #include <string>
31 #include <vector>
32 #include <algorithm>
33 
34 #include "../exceptions.h"
35 #include "../population.h"
36 #include "../types.h"
37 #include "base.h"
38 #include "sga.h"
39 #include "../problem/base_stochastic.h"
40 
41 namespace pagmo { namespace algorithm {
42 
44 
59 sga::sga(int gen, const double &cr, const double &m, int elitism, mutation::type mut, double width, selection::type sel, crossover::type cro)
60  :base(),m_gen(gen),m_cr(cr),m_m(m),m_elitism(elitism),m_mut(mut,width),m_sel(sel),m_cro(cro)
61 {
62  if (gen < 0) {
63  pagmo_throw(value_error,"number of generations must be nonnegative");
64  }
65  if (cr > 1 || cr < 0) {
66  pagmo_throw(value_error,"crossover probability must be in the [0,1] range");
67  }
68  if (m < 0 || m > 1) {
69  pagmo_throw(value_error,"mutation probability must be in the [0,1] range");
70  }
71  if (elitism < 1) {
72  pagmo_throw(value_error,"elitisim must be greater than zero");
73  }
74  if (width <0 || width >1) {
75  pagmo_throw(value_error,"mutation width must be in the [0,1] range");
76  }
77 
78 }
79 
82 {
83  return base_ptr(new sga(*this));
84 }
85 
87 
94 void sga::evolve(population &pop) const
95 {
96  // Let's store some useful variables.
97  const problem::base &prob = pop.problem();
98  const problem::base::size_type D = prob.get_dimension(), Di = prob.get_i_dimension(), prob_c_dimension = prob.get_c_dimension(), prob_f_dimension = prob.get_f_dimension();
99  const decision_vector &lb = prob.get_lb(), &ub = prob.get_ub();
100  const population::size_type NP = pop.size();
101  const problem::base::size_type Dc = D - Di;
102 
103 
104  //We perform some checks to determine wether the problem/population are suitable for SGA
105  if ( prob_c_dimension != 0 ) {
106  pagmo_throw(value_error,"The problem is not box constrained and SGA is not suitable to solve it");
107  }
108 
109  if ( prob_f_dimension != 1 ) {
110  pagmo_throw(value_error,"The problem is not single objective and SGA is not suitable to solve it");
111  }
112 
113  if (NP < 5) {
114  pagmo_throw(value_error,"for SGA at least 5 individuals in the population are needed");
115  }
116 
117  // Get out if there is nothing to do.
118  if (m_gen == 0) {
119  return;
120  }
121  // Some vectors used during evolution are allocated here.
122  decision_vector dummy(D,0); //used for initialisation purposes
123  std::vector<decision_vector > X(NP,dummy), Xnew(NP,dummy);
124 
125  std::vector<fitness_vector > fit(NP); //fitness
126 
127  fitness_vector bestfit;
128  decision_vector bestX(D,0);
129 
130  std::vector<double> selectionfitness(NP), cumsum(NP), cumsumTemp(NP);
131  std::vector <int> selection(NP);
132 
133  int tempID;
134  std::vector<int> fitnessID(NP);
135 
136  // Initialise the chromosomes and their fitness to that of the initial deme
137  for (pagmo::population::size_type i = 0; i<NP; i++ ) {
138  X[i] = pop.get_individual(i).cur_x;
139  fit[i] = pop.get_individual(i).cur_f;
140  }
141 
142  // Find the best member and store in bestX and bestfit
143  double bestidx = pop.get_best_idx();
144  bestX = pop.get_individual(bestidx).cur_x;
145  bestfit = pop.get_individual(bestidx).cur_f;
146 
147 
148  // Main SGA loop
149  for (int j = 0; j<m_gen; j++) {
150 
151  switch (m_sel) {
152  case selection::BEST20: { //selects the best 20% and puts multiple copies in Xnew
153  //Sort the individuals according to their fitness
154  for (pagmo::population::size_type i=0; i<NP; i++) fitnessID[i]=i;
155  for (pagmo::population::size_type i=0; i < (NP-1); ++i) {
156  for (pagmo::population::size_type j=i+1; j<NP; ++j) {
157  if ( prob.compare_fitness(fit[j],fit[i]) ) {
158  //swap fitness values
159  fit[i].swap(fit[j]);
160  //swap id's
161  tempID = fitnessID[i];
162  fitnessID[i] = fitnessID[j];
163  fitnessID[j] = tempID;
164  }
165  }
166  }
167  int best20 = NP/5;
168  for (pagmo::population::size_type i=0; i<NP; ++i) {
169  selection[i] = fitnessID[i % best20];
170  }
171  break;
172  }
173 
174  case selection::ROULETTE: {
175  //We scale all fitness values from 0 (worst) to absolute value of the best fitness
176  fitness_vector worstfit=fit[0];
177  for (pagmo::population::size_type i = 1; i < NP;i++) {
178  if (prob.compare_fitness(worstfit,fit[i])) worstfit=fit[i];
179  }
180 
181  for (pagmo::population::size_type i = 0; i < NP; i++) {
182  selectionfitness[i] = fabs(worstfit[0] - fit[i][0]);
183  }
184 
185  // We build and normalise the cumulative sum
186  cumsumTemp[0] = selectionfitness[0];
187  for (pagmo::population::size_type i = 1; i< NP; i++) {
188  cumsumTemp[i] = cumsumTemp[i - 1] + selectionfitness[i];
189  }
190  for (pagmo::population::size_type i = 0; i < NP; i++) {
191  cumsum[i] = cumsumTemp[i]/cumsumTemp[NP-1];
192  }
193 
194  //we throw a dice and pick up the corresponding index
195  double r2;
196  for (pagmo::population::size_type i = 0; i < NP; i++) {
197  r2 = m_drng();
198  for (pagmo::population::size_type j = 0; j < NP; j++) {
199  if (cumsum[j] > r2) {
200  selection[i]=j;
201  break;
202  }
203  }
204  }
205  break;
206  }
207  }
208 
209  //Xnew stores the new selected generation of chromosomes
210  for (pagmo::population::size_type i = 0; i < NP; i++) {
211  Xnew[i]=X[selection[i]];
212  }
213 
214  //2 - Crossover
215  {
216  int r1,L;
217  decision_vector member1,member2;
218 
219  for (pagmo::population::size_type i=0; i< NP; i++) {
220  //for each chromosome selected i.e. in Xnew
221  member1 = Xnew[i];
222  //we select a mating patner different from the self (i.e. no masturbation)
223  do {
224  r1 = boost::uniform_int<int>(0,NP - 1)(m_urng);
225  } while ( r1 == boost::numeric_cast<int>(i) );
226  member2 = Xnew[r1];
227  //and we operate crossover
228  switch (m_cro) {
229  //0 - binomial crossover
230  case crossover::BINOMIAL: {
231  size_t n = boost::uniform_int<int>(0,D-1)(m_urng);
232  for (size_t L = 0; L < D; ++L) { /* perform D binomial trials */
233  if ((m_drng() < m_cr) || L + 1 == D) { /* change at least one parameter */
234  member1[n] = member2[n];
235  }
236  n = (n+1)%D;
237  }
238  break; }
239  //1 - exponential crossover
240  case crossover::EXPONENTIAL: {
241  size_t n = boost::uniform_int<int>(0,D-1)(m_urng);
242  L = 0;
243  do {
244  member1[n] = member2[n];
245  n = (n+1) % D;
246  L++;
247  } while ( (m_drng() < m_cr) && (L < boost::numeric_cast<int>(D)) );
248  break; }
249  }
250  Xnew[i] = member1;
251 
252  } }
253 
254  //3 - Mutation
255  switch (m_mut.m_type) {
256  case mutation::GAUSSIAN: {
257  boost::normal_distribution<double> dist;
258  boost::variate_generator<boost::lagged_fibonacci607 &, boost::normal_distribution<double> > delta(m_drng,dist);
259  for (pagmo::problem::base::size_type k = 0; k < Dc;k++) { //for each continuous variable
260  double std = (ub[k]-lb[k]) * m_mut.m_width;
261  for (pagmo::population::size_type i = 0; i < NP;i++) { //for each individual
262  if (m_drng() < m_m) {
263  double mean = Xnew[i][k];
264  double tmp = (delta() * std + mean);
265  if ( (tmp < ub[k]) && (tmp > lb[k]) ) Xnew[i][k] = tmp;
266  }
267  }
268  }
269  for (pagmo::problem::base::size_type k = Dc; k < D;k++) { //for each integer variable
270  double std = (ub[k]-lb[k]) * m_mut.m_width;
271  for (pagmo::population::size_type i = 0; i < NP;i++) { //for each individual
272  if (m_drng() < m_m) {
273  double mean = Xnew[i][k];
274  double tmp = boost::math::iround(delta() * std + mean);
275  if ( (tmp < ub[k]) && (tmp > lb[k]) ) Xnew[i][k] = tmp;
276  }
277  }
278  }
279  break;
280  }
281  case mutation::RANDOM: {
282  for (pagmo::population::size_type i = 0; i < NP;i++) {
283  for (pagmo::problem::base::size_type j = 0; j < Dc;j++) { //for each continuous variable
284  if (m_drng() < m_m) {
285  Xnew[i][j] = boost::uniform_real<double>(lb[j],ub[j])(m_drng);
286  }
287  }
288  for (pagmo::problem::base::size_type j = Dc; j < D;j++) {//for each integer variable
289  if (m_drng() < m_m) {
290  Xnew[i][j] = boost::uniform_int<int>(lb[j],ub[j])(m_urng);
291  }
292  }
293  }
294  break;
295  }
296  }
297 
298  // If the problem is a stochastic optimization chage the seed and re-evaluate taking care to update also best and local bests
299  try
300  {
301  //4 - Evaluate the new population (stochastic problem)
302  dynamic_cast<const pagmo::problem::base_stochastic &>(prob).set_seed(m_urng());
303  pop.clear(); // Removes memory based on different seeds (champion and best_x, best_f, best_c)
304 
305  // We re-evaluate the best individual (for elitism)
306  prob.objfun(bestfit,bestX);
307  // Re-evaluate wrt new seed the particle position and memory
308  for (pagmo::population::size_type i = 0; i < NP;i++) {
309  // We evaluate here the new individual fitness
310  prob.objfun(fit[i],Xnew[i]);
311  // We update the velocity (in case coupling with PSO via archipelago)
312  //dummy = Xnew[i];
313  //std::transform(dummy.begin(), dummy.end(), pop.get_individual(i).cur_x.begin(), dummy.begin(),std::minus<double>());
315  pop.push_back(Xnew[i]);
316  //pop.set_v(i,dummy);
317  if (prob.compare_fitness(fit[i], bestfit)) {
318  bestfit = fit[i];
319  bestX = Xnew[i];
320  }
321  }
322  }
323  catch (const std::bad_cast& e)
324  {
325  //4 - Evaluate the new population (deterministic problem)
326  for (pagmo::population::size_type i = 0; i < NP;i++) {
327  prob.objfun(fit[i],Xnew[i]);
328  dummy = Xnew[i];
329  std::transform(dummy.begin(), dummy.end(), pop.get_individual(i).cur_x.begin(), dummy.begin(),std::minus<double>());
330  //updates x and v (cache avoids to recompute the objective function and constraints)
331  pop.set_x(i,Xnew[i]);
332  pop.set_v(i,dummy);
333  if (prob.compare_fitness(fit[i], bestfit)) {
334  bestfit = fit[i];
335  bestX = Xnew[i];
336  }
337  }
338  }
339 
340  //5 - Reinsert best individual every m_elitism generations
341  if (j % m_elitism == 0) {
342  int worst=0;
343  for (pagmo::population::size_type i = 1; i < NP;i++) {
344  if ( prob.compare_fitness(fit[worst],fit[i]) ) worst=i;
345  }
346  Xnew[worst] = bestX;
347  fit[worst] = bestfit;
348  dummy = Xnew[worst];
349  std::transform(dummy.begin(), dummy.end(), pop.get_individual(worst).cur_x.begin(), dummy.begin(),std::minus<double>());
350  //updates x and v (cache avoids to recompute the objective function)
351  pop.set_x(worst,Xnew[worst]);
352  pop.set_v(worst,dummy);
353  }
354  X = Xnew;
355  } // end of main SGA loop
356 }
357 
359 std::string sga::get_name() const
360 {
361  return "A Simple Genetic Algorithm";
362 }
363 
365 
368 std::string sga::human_readable_extra() const
369 {
370  std::ostringstream s;
371  s << "gen:" << m_gen << ' ';
372  s << "CR:" << m_cr << ' ';
373  s << "M:" << m_m << ' ';
374  s << "elitism:" << m_elitism << ' ';
375  s << "mutation:";
376  switch (m_mut.m_type) {
377  case mutation::RANDOM: {
378  s << "RANDOM ";
379  break;
380  }
381  case mutation::GAUSSIAN: {
382  s << "GAUSSIAN (" << m_mut.m_width << ") ";
383  break;
384  }
385  }
386  s << "selection:";
387  switch (m_sel) {
388  case selection::BEST20: {
389  s << "BEST20 ";
390  break;
391  }
392  case selection::ROULETTE: {
393  s << "ROULETTE ";
394  break;
395  }
396  }
397  s << "crossover:";
398  switch (m_cro) {
399  case crossover::EXPONENTIAL: {
400  s << "EXPONENTIAL ";
401  break;
402  }
403  case crossover::BINOMIAL: {
404  s << "BINOMIAL ";
405  break;
406  }
407  }
408 
409  return s.str();
410 }
411 
412 }} //namespaces
413 
414 BOOST_CLASS_EXPORT_IMPLEMENT(pagmo::algorithm::sga)
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
type
Selection type, best 20% or roulette.
Definition: sga.h:59
fitness_vector cur_f
Current fitness vector.
Definition: population.h:89
const individual_type & get_individual(const size_type &) const
Get constant reference to individual at position n.
Definition: population.cpp:277
Selection info.
Definition: sga.h:57
Base algorithm class.
type
Crossover type, binomial or exponential.
Definition: sga.h:90
STL namespace.
Base problem class.
Definition: problem/base.h:148
Population class.
Definition: population.h:70
size_type get_dimension() const
Return global dimension.
void evolve(population &) const
Evolve implementation.
Definition: sga.cpp:94
bool compare_fitness(const fitness_vector &, const fitness_vector &) const
Compare fitness vectors.
double m_width
Mutation width.
Definition: sga.h:76
fitness_vector objfun(const decision_vector &) const
Return fitness of pagmo::decision_vector.
size_type get_i_dimension() const
Return integer dimension.
std::vector< double > fitness_vector
Fitness vector type.
Definition: types.h:42
c_size_type get_c_dimension() const
Return global constraints dimension.
const decision_vector & get_ub() const
Upper bounds getter.
Base Stochastic Optimization Problem.
type
Mutation type, gaussian or random.
Definition: sga.h:64
base_ptr clone() const
Clone method.
Definition: sga.cpp:81
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.
sga(int gen=1, const double &cr=.95, const double &m=.02, int elitism=1, mutation::type mut=mutation::GAUSSIAN, double width=0.1, selection::type sel=selection::ROULETTE, crossover::type cro=crossover::EXPONENTIAL)
Constructor.
Definition: sga.cpp:59
f_size_type get_f_dimension() const
Return fitness dimension.
const decision_vector & get_lb() const
Lower bounds getter.
rng_double m_drng
Random number generator for double-precision floating point values.
std::string get_name() const
Algorithm name.
Definition: sga.cpp:359
The Simple Genetic Algorithm (SGA)
Definition: sga.h:53
std::string human_readable_extra() const
Extra human readable algorithm info.
Definition: sga.cpp:368
type m_type
Mutation type.
Definition: sga.h:74
decision_vector::size_type size_type
Problem's size type: the same as pagmo::decision_vector's size type.
Definition: problem/base.h:160