diff_evo.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2010-2014, Edwin van Leeuwen
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_DIFF_EVO_H
24 #define O2SCL_DIFF_EVO_H
25 
26 /** \file diff_evo.h
27  \brief File defining \ref o2scl::diff_evo
28 */
29 
30 #include <vector>
31 #include <algorithm>
32 
33 #include <o2scl/rng_gsl.h>
34 #include <o2scl/mmin.h>
35 #include <o2scl/mm_funct.h>
36 
37 #ifndef DOXYGEN_NO_O2NS
38 namespace o2scl {
39 #endif
40 
41  /** \brief Multidimensional minimization by the differential
42  evolution method
43 
44  This class minimizes a function using differential evolution.
45  This method is a genetic algorithm and as such works well for
46  non continuous problems, since it does not rely on a gradient of
47  the function that is being mind.
48 
49  The method starts by initializing a random population of
50  candidate parameters. To do this the user needs to define a
51  function to create these random parameters, which can be
52  provided using \ref set_init_function().
53 
54  After the initial population is created the algorithm will
55  repeat a number of standard steps until a solution is found or
56  the maximum number of iterations is reached. Based on
57  \ref Storn97.
58 
59  \note The constructor sets \ref o2scl::mmin_base::ntrial to 1000 .
60 
61  If the population converges prematurely, then \ref diff_evo::f
62  and \ref pop_size should be increased.
63  */
64  template<class func_t=multi_funct11,
66  class init_funct_t=mm_funct11> class diff_evo :
67  public mmin_base<func_t,func_t,vec_t> {
68 
69  public:
70 
72 
73  /** \brief Population size (default 0)
74 
75  Should be at least 4. Typically between \f$ 5 d \f$ and \f$ 10
76  d \f$ where \f$ d \f$ is the dimensionality of the problem.
77 
78  If this is 0 (the default), then it is set by mmin to be
79  equal to \f$ 10 d \f$ .
80  */
81  size_t pop_size;
82 
83  /** \brief The number of generations without a better fit before we
84  assume that the algorithm has converged (default 25)
85  */
86  size_t nconv;
87 
88  /** \brief Differential weight (default 0.75)
89 
90  A parameter which controls the amplification of the
91  differential variation. Usually between 0 and 2.
92  */
93  double f;
94 
95  /** \brief Crossover probability (default 0.8)
96 
97  Usually between 0 and 1.
98  */
99  double cr;
100 
101  diff_evo() {
102  this->ntrial=1000;
103  f = 0.75;
104  cr = 0.8;
105  rand_init_funct = 0;
106  pop_size = 0;
107  nconv = 25;
108  }
109 
110  virtual ~diff_evo() {
111  }
112 
113  /** \brief Set the function that is used to produce random
114  init variables
115 
116  REQUIRED
117 
118  The init function is called in the beginning to fill
119  the population with random individuals, so it is best
120  to make this cover the part of the parameter space you
121  are interested in. The method will find solutions outside
122  this parameter space, but choosing a good init function will
123  help finding solutions faster.
124  */
125  virtual void set_init_function( init_funct_t &function ) {
126  rand_init_funct = &function;
127  }
128 
129  /** \brief Calculate the minimum \c fmin of \c func w.r.t the
130  array \c x of size \c nvar.
131 
132  Initialize all agents x with random positions in the
133  search-space. Until a termination criterion is met (e.g.
134  number of iterations performed, or adequate fitness reached),
135  repeat the following: For each agent x in the population do:
136  Pick three agents a, b, and c from the population at random,
137  they must be distinct from each other as well as from agent x
138  Pick a random index {1, ..., n}, where the highest possible
139  value n is the dimensionality of the problem to be optimized.
140  Compute the agent's potentially new position y = [y1, ..., yn]
141  by iterating over each i {1, ..., n} as follows: Pick
142  ri~U(0,1) uniformly from the open range (0,1) If (i=R) or
143  (ri<CR) let yi = ai + F(bi - ci), otherwise let yi = xi If
144  (f(y) < f(x)) then replace the agent in the population with
145  the improved candidate solution, that is, set x = y in the
146  population.
147 
148  Pick the agent from the population that has the lowest fmin
149  and return it as the best found candidate solution.
150  */
151  virtual int mmin(size_t nvar, vec_t &x0, double &fmin, func_t &func) {
152 
153  // Keep track of number of generation without better solutions
154  size_t nconverged = 0;
155 
156  if (pop_size==0) {
157  // Automatically select pop_size dependent on dimensionality.
158  pop_size = 10*nvar;
159  }
160 
161  initialize_population( nvar, x0 );
162 
163  fmins.resize(pop_size);
164 
165  // Set initial fmin
166  for (size_t x = 0; x < pop_size; ++x) {
167  vec_t agent_x;
168  agent_x.resize(nvar);
169  for (size_t i = 0; i < nvar; ++i) {
170  agent_x[i] = population[x*nvar+i];
171  }
172  double fmin_x = 0;
173  fmin_x=func(nvar,agent_x);
174  fmins[x]=fmin_x;
175  if (x==0) {
176  fmin = fmin_x;
177  for (size_t i = 0; i<nvar; ++i)
178  x0[i] = agent_x[i];
179  //x0 = agent_x;
180  } else if (fmin_x<fmin) {
181  fmin = fmin_x;
182  for (size_t i = 0; i<nvar; ++i)
183  x0[i] = agent_x[i];
184  //x0 = agent_x;
185  }
186  }
187 
188  int gen = 0;
189  while (gen < this->ntrial && nconverged <= nconv) {
190 
191  ++nconverged;
192  ++gen;
193 
194  // For each agent x in the population do:
195  for (size_t x = 0; x < pop_size; ++x) {
196 
197  std::vector<int> others;
198 
199  // Create a copy agent_x and agent_y of the current agent
200  // vector
201  vec_t agent_x, agent_y;
202  agent_x.resize(nvar);
203  agent_y.resize(nvar);
204  for (size_t i = 0; i < nvar; ++i) {
205  agent_x[i] = population[x*nvar+i];
206  agent_y[i] = population[x*nvar+i];
207  }
208 
209  // Pick three agents a, b, and c from the population at
210  // random, they must be distinct from each other as well as
211  // from agent x
212  others = pick_unique_agents( 3, x );
213 
214  // Pick a random index R in {1, ..., n}, where the highest
215  // possible value n is the dimensionality of the problem
216  // to be optimized.
217  size_t r = floor(gr.random()*nvar);
218 
219  for (size_t i = 0; i < nvar; ++i) {
220  // Pick ri~U(0,1) uniformly from the open range (0,1)
221  double ri = gr.random();
222  // If (i=R) or (ri<CR) let yi = ai + F(bi - ci), otherwise
223  // let yi = xi
224  if (i == r || ri < cr) {
225  agent_y[i] = population[others[0]*nvar+i] +
226  f*(population[others[1]*nvar+i]-
227  population[others[2]*nvar+i]);
228  }
229  }
230  // If (f(y) < f(x)) then replace the agent in the population
231  // with the improved candidate solution, that is, set x = y
232  // in the population
233  double fmin_y;
234 
235  fmin_y=func(nvar,agent_y);
236  if (fmin_y<fmins[x]) {
237  for (size_t i = 0; i < nvar; ++i) {
238  population[x*nvar+i] = agent_y[i];
239  fmins[x] = fmin_y;
240  }
241  if (fmin_y<fmin) {
242  fmin = fmin_y;
243  for (size_t i = 0; i<nvar; ++i) {
244  x0[i] = agent_y[i];
245  }
246  nconverged = 0;
247  }
248  }
249 
250  }
251  if (this->verbose > 0)
252  this->print_iter( nvar, fmin, gen, x0 );
253  }
254 
255  if(gen>=this->ntrial) {
256  std::string str="Exceeded maximum number of iterations ("+
257  itos(this->ntrial)+") in diff_evo::mmin().";
258  O2SCL_CONV_RET(str.c_str(),exc_emaxiter,this->err_nonconv);
259  }
260 
261  return 0;
262  };
263 
264  /** \brief Print out iteration information.
265 
266  \comment
267  Depending on the value of the variable verbose, this prints
268  out the iteration information. If verbose=0, then no
269  information is printed, while if verbose>1, then after each
270  iteration, the present values of x and y are output to
271  std::cout along with the iteration number. If verbose>=2 then
272  each iteration waits for a character.
273  \endcomment
274  */
275  virtual void print_iter( size_t nvar, double fmin,
276  int iter, vec_t &best_fit ) {
277  std::cout << "Generation " << iter << std::endl;
278  std::cout << "Fmin: " << fmin << std::endl;
279  std::cout << "Parameters: ";
280  for (size_t i=0; i<nvar; ++i) {
281  std::cout << best_fit[i] << " ";
282  }
283  std::cout << std::endl;
284  std::cout << "Population: " << std::endl;
285  for (size_t i=0; i<pop_size; ++i ) {
286  std::cout << i << ": ";
287  for (size_t j = 0; j<nvar; ++j ) {
288  std::cout << population[i*nvar+j] << " ";
289  }
290  std::cout << "fmin: " << fmins[i] << std::endl;
291  }
292  }
293 
294 #ifndef DOXYGEN_INTERNAL
295 
296  protected:
297 
298  /** \brief Vector containing the population.
299 
300  For now using one long vector with all agents after each other
301  */
302  vec_t population;
303 
304  /// Vector that keeps track of fmins values
305  ubvector fmins;
306 
307  /** \brief Function that is used to produce random init variables
308 
309  This function is used to fill the population with random agents
310  */
311  init_funct_t *rand_init_funct;
312 
313  /// Random number generator
315 
316  /** \brief Initialize a population of random agents
317  */
318  virtual int initialize_population( size_t nvar, vec_t &x0 ) {
319  if (rand_init_funct==0) {
320  O2SCL_ERR("No initialization function provided.",
321  exc_ebadfunc);
322 
323  }
324  population.resize(nvar*pop_size );
325  for (size_t i = 0; i < pop_size; ++i) {
326  vec_t y(nvar);
327  (*rand_init_funct)( nvar, x0, y );
328  for (size_t j = 0; j < nvar; ++j) {
329  population[ i*nvar+j ] = y[j];
330 
331  }
332  }
333  return 0;
334  }
335 
336  /** \brief Pick number of unique agent id's
337 
338  Unique from x and each other
339 
340  Uses the Fisher-Yates algorithm.
341 
342  */
343  virtual std::vector<int> pick_unique_agents( int nr, size_t x ) {
344  std::vector<int> ids;
345  std::vector<int> agents;
346  // Fill array with ids
347  for (size_t i=0; i<pop_size-1; ++i){
348  if (i<x) {
349  ids.push_back( i );
350  } else {
351  ids.push_back( i+1 );
352  }
353  }
354  // Shuffle according to Fisher-Yates
355  for (size_t i=ids.size()-1; i>ids.size()-nr-1; --i) {
356  int j = round(gr.random()*i);
357  std::swap( ids[i], ids[j] );
358  }
359  for (size_t i=ids.size()-1; i>ids.size()-nr-1; --i) {
360  agents.push_back( ids[i] );
361  }
362  return agents;
363  }
364 
365 #endif
366 
367 #ifndef DOXYGEN_INTERNAL
368 
369  private:
370 
374  (const diff_evo<func_t,vec_t,init_funct_t>&);
375 
376 #endif
377 
378  };
379 
380 #ifndef DOXYGEN_NO_O2NS
381 }
382 #endif
383 
384 #endif
size_t nconv
The number of generations without a better fit before we assume that the algorithm has converged (def...
Definition: diff_evo.h:86
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:292
Multidimensional minimization by the differential evolution method.
Definition: diff_evo.h:66
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct11
Array of multi-dimensional functions typedef.
Definition: mm_funct.h:43
exceeded max number of iterations
Definition: err_hnd.h:73
virtual int initialize_population(size_t nvar, vec_t &x0)
Initialize a population of random agents.
Definition: diff_evo.h:318
double random()
Return a random number in .
Definition: rng_gsl.h:82
size_t pop_size
Population size (default 0)
Definition: diff_evo.h:81
bool err_nonconv
If true, call the error handler if the routine does not "converge".
Definition: mmin.h:209
Multidimensional minimization [abstract base].
Definition: mmin.h:164
virtual void set_init_function(init_funct_t &function)
Set the function that is used to produce random init variables.
Definition: diff_evo.h:125
double cr
Crossover probability (default 0.8)
Definition: diff_evo.h:99
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
ubvector fmins
Vector that keeps track of fmins values.
Definition: diff_evo.h:305
init_funct_t * rand_init_funct
Function that is used to produce random init variables.
Definition: diff_evo.h:311
virtual int mmin(size_t nvar, vec_t &x0, double &fmin, func_t &func)
Calculate the minimum fmin of func w.r.t the array x of size nvar.
Definition: diff_evo.h:151
double f
Differential weight (default 0.75)
Definition: diff_evo.h:93
virtual void print_iter(size_t nvar, double fmin, int iter, vec_t &best_fit)
Print out iteration information.
Definition: diff_evo.h:275
virtual std::vector< int > pick_unique_agents(int nr, size_t x)
Pick number of unique agent id&#39;s.
Definition: diff_evo.h:343
Random number generator (GSL)
Definition: rng_gsl.h:55
problem with user-supplied function
Definition: err_hnd.h:69
vec_t population
Vector containing the population.
Definition: diff_evo.h:302
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct11
Multi-dimensional function typedef.
Definition: multi_funct.h:45
rng_gsl gr
Random number generator.
Definition: diff_evo.h:314
std::string itos(int x)
Convert an integer to a string.
int ntrial
Maximum number of iterations.
Definition: mmin.h:197

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).