diff_evo.h
Go to the documentation of this file.
1  /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2010-2020, Edwin van Leeuwen and Andrew W. Steiner
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 require the gradient
47  of the function to be minimized.
48 
49  The method starts by initializing a the population of
50  points in the parameter space. The user can specify
51  a function which creates the initial population using
52  \ref set_init_function() . Alternatively, by default this
53  class chooses random points near the initial point
54  with a default step size of 0.01 in every direction.
55  The step size can be changed using \ref set_step() .
56 
57  After the initial population is created the algorithm will
58  repeat until a solution is found or the maximum number of
59  iterations is reached. Based on \ref Storn97.
60 
61  \note The constructor sets \ref o2scl::mmin_base::ntrial to 1000 .
62 
63  If the population converges prematurely, then \ref diff_evo::f
64  and \ref pop_size should be increased.
65 
66  \future AWS, 7/25/19, The code stops early if \ref nconv
67  generations go by without a better fit, but we could also
68  consider some code which terminates early if the minimum is
69  found to within a particular tolerance.
70 
71  \future This class probably has a simple parallelization like
72  \ref o2scl::anneal_para ?
73  */
74  template<class func_t=multi_funct,
76  class init_funct_t=mm_funct> class diff_evo :
77  public mmin_base<func_t,func_t,vec_t> {
78 
79  public:
80 
82 
83  /** \brief Population size (default 0)
84 
85  Should be at least 4. Typically between \f$ 5 d \f$ and \f$ 10
86  d \f$ where \f$ d \f$ is the dimensionality of the problem.
87 
88  If this is 0 (the default), then it is set by mmin to be
89  equal to \f$ 10 d \f$ .
90  */
91  size_t pop_size;
92 
93  /** \brief The number of generations without a better fit before we
94  assume that the algorithm has converged (default 25)
95  */
96  size_t nconv;
97 
98  /** \brief Differential weight (default 0.75)
99 
100  A parameter which controls the amplification of the
101  differential variation. Usually between 0 and 2.
102  */
103  double f;
104 
105  /** \brief Crossover probability (default 0.8)
106 
107  Usually between 0 and 1.
108  */
109  double cr;
110 
111  diff_evo() {
112  this->ntrial=1000;
113  f = 0.75;
114  cr = 0.8;
115  rand_init_funct = 0;
116  pop_size = 0;
117  nconv = 25;
118  step.resize(1);
119  step[0]=1.0e-2;
120  }
121 
122  virtual ~diff_evo() {
123  }
124 
125  /** \brief Set the function that is used to select the
126  initial population
127 
128  The init function is called in the beginning to fill
129  the population with random individuals, so it is best
130  to make this cover the part of the parameter space you
131  are interested in. The method will find solutions outside
132  this parameter space, but choosing a good init function will
133  help finding solutions faster.
134 
135  Note that this function stores a pointer to the user-specified
136  function so care must be taken to ensure the pointer is still
137  valid when the minimization is run.
138  */
139  virtual void set_init_function(init_funct_t &function) {
140  rand_init_funct = &function;
141  }
142 
143  /** \brief Calculate the minimum \c fmin of \c func w.r.t the
144  array \c x of size \c nvar.
145 
146  Initialize all agents x with random positions in the
147  search-space. Until a termination criterion is met (e.g.
148  number of iterations performed, or adequate fitness reached),
149  repeat the following: For each agent x in the population do:
150  Pick three agents a, b, and c from the population at random,
151  they must be distinct from each other as well as from agent x
152  Pick a random index {1, ..., n}, where the highest possible
153  value n is the dimensionality of the problem to be optimized.
154  Compute the agent's potentially new position y = [y1, ..., yn]
155  by iterating over each i {1, ..., n} as follows: Pick
156  ri~U(0,1) uniformly from the open range (0,1) If (i=R) or
157  (ri<CR) let yi = ai + F(bi - ci), otherwise let yi = xi If
158  (f(y) < f(x)) then replace the agent in the population with
159  the improved candidate solution, that is, set x = y in the
160  population.
161 
162  Pick the agent from the population that has the lowest fmin
163  and return it as the best found candidate solution.
164  */
165  virtual int mmin(size_t nvar, vec_t &x0, double &fmin, func_t &func) {
166 
167  // Keep track of number of generation without better solutions
168  size_t nconverged = 0;
169 
170  if (pop_size==0) {
171  // Automatically select pop_size based on on dimensionality.
172  pop_size = 10*nvar;
173  }
174 
175  initialize_population( nvar, x0 );
176 
177  fmins.resize(pop_size);
178 
179  // Set initial fmin
180  for (size_t x = 0; x < pop_size; ++x) {
181  vec_t agent_x;
182  agent_x.resize(nvar);
183  for (size_t i = 0; i < nvar; ++i) {
184  agent_x[i] = population[x*nvar+i];
185  }
186  double fmin_x = 0;
187  fmin_x=func(nvar,agent_x);
188  fmins[x]=fmin_x;
189  if (x==0) {
190  fmin = fmin_x;
191  for (size_t i = 0; i<nvar; ++i) {
192  x0[i] = agent_x[i];
193  }
194  } else if (fmin_x<fmin) {
195  fmin = fmin_x;
196  for (size_t i = 0; i<nvar; ++i) {
197  x0[i] = agent_x[i];
198  }
199  }
200  }
201 
202  int gen = 0;
203  while (gen < this->ntrial && nconverged <= nconv) {
204 
205  ++nconverged;
206  ++gen;
207 
208  // For each agent x in the population do:
209  for (size_t x = 0; x < pop_size; ++x) {
210 
211  std::vector<int> others;
212 
213  // Create a copy agent_x and agent_y of the current agent
214  // vector
215  vec_t agent_x, agent_y;
216  agent_x.resize(nvar);
217  agent_y.resize(nvar);
218  for (size_t i = 0; i < nvar; ++i) {
219  agent_x[i] = population[x*nvar+i];
220  agent_y[i] = population[x*nvar+i];
221  }
222 
223  // Pick three agents a, b, and c from the population at
224  // random, they must be distinct from each other as well as
225  // from agent x
226  others = pick_unique_agents( 3, x );
227 
228  // Pick a random index R in {1, ..., n}, where the highest
229  // possible value n is the dimensionality of the problem
230  // to be optimized.
231  size_t r = floor(gr.random()*nvar);
232 
233  for (size_t i = 0; i < nvar; ++i) {
234  // Pick ri~U(0,1) uniformly from the open range (0,1)
235  double ri = gr.random();
236  // If (i=R) or (ri<CR) let yi = ai + F(bi - ci), otherwise
237  // let yi = xi
238  if (i == r || ri < cr) {
239  agent_y[i] = population[others[0]*nvar+i] +
240  f*(population[others[1]*nvar+i]-
241  population[others[2]*nvar+i]);
242  }
243  }
244  // If (f(y) < f(x)) then replace the agent in the population
245  // with the improved candidate solution, that is, set x = y
246  // in the population
247  double fmin_y;
248 
249  fmin_y=func(nvar,agent_y);
250  if (fmin_y<fmins[x]) {
251  for (size_t i = 0; i < nvar; ++i) {
252  population[x*nvar+i] = agent_y[i];
253  fmins[x] = fmin_y;
254  }
255  if (fmin_y<fmin) {
256  fmin = fmin_y;
257  for (size_t i = 0; i<nvar; ++i) {
258  x0[i] = agent_y[i];
259  }
260  nconverged = 0;
261  }
262  }
263 
264  }
265  if (this->verbose > 0) {
266  this->print_iter( nvar, fmin, gen, x0 );
267  }
268  }
269 
270  this->last_ntrial=gen;
271 
272  if (gen>=this->ntrial) {
273  std::string str="Exceeded maximum number of iterations ("+
274  itos(this->ntrial)+") in diff_evo::mmin().";
275  O2SCL_CONV_RET(str.c_str(),exc_emaxiter,this->err_nonconv);
276  }
277 
278  return 0;
279  };
280 
281  /** \brief Print out iteration information.
282 
283  \comment
284  Depending on the value of the variable verbose, this prints
285  out the iteration information. If verbose=0, then no
286  information is printed, while if verbose>1, then after each
287  iteration, the present values of x and y are output to
288  std::cout along with the iteration number. If verbose>=2 then
289  each iteration waits for a character.
290  \endcomment
291  */
292  virtual void print_iter( size_t nvar, double fmin,
293  int iter, vec_t &best_fit ) {
294  std::cout << "Generation " << iter << std::endl;
295  std::cout << "Fmin: " << fmin << std::endl;
296  std::cout << "Parameters: ";
297  for (size_t i=0; i<nvar; ++i) {
298  std::cout << best_fit[i] << " ";
299  }
300  std::cout << std::endl;
301  std::cout << "Population: " << std::endl;
302  for (size_t i=0; i<pop_size; ++i ) {
303  std::cout << i << ": ";
304  for (size_t j = 0; j<nvar; ++j ) {
305  std::cout << population[i*nvar+j] << " ";
306  }
307  std::cout << "fmin: " << fmins[i] << std::endl;
308  }
309  if (this->verbose>1) {
310  char ch;
311  std::cin >> ch;
312  }
313  return;
314  }
315 
316  /// Set the step sizes for the default initialization
317  template<class vec2_t> int set_step(size_t nv, vec2_t &stepv) {
318  if (nv>0) {
319  step.resize(nv);
320  for(size_t i=0;i<nv;i++) step[i]=stepv[i];
321  }
322  return 0;
323  }
324 
325 #ifndef DOXYGEN_INTERNAL
326 
327  protected:
328 
329  /** \brief Vector containing the population.
330 
331  A long vector with all agents after each other
332  */
333  vec_t population;
334 
335  /// Vector that keeps track of the function values
337 
338  /** \brief Function that is used to produce random init variables
339 
340  This function is used to fill the population with random agents
341  */
342  init_funct_t *rand_init_funct;
343 
344  /// Random number generator
346 
347  /// Step size for initialization
348  std::vector<double> step;
349 
350  /** \brief Initialize a population of random agents
351  */
352  virtual int initialize_population( size_t nvar, vec_t &x0 ) {
353  population.resize(nvar*pop_size);
354  if (rand_init_funct==0) {
355  for(size_t i=0;i<pop_size;i++) {
356  for(size_t j=0;j<nvar;j++) {
357  double stepj=step[j%step.size()];
358  population[i*nvar+j]=x0[j]-stepj/2.0+stepj*gr.random();
359  }
360  }
361  } else {
362  for (size_t i = 0; i < pop_size; ++i) {
363  vec_t y(nvar);
364  (*rand_init_funct)( nvar, x0, y );
365  for (size_t j = 0; j < nvar; ++j) {
366  population[ i*nvar+j ] = y[j];
367 
368  }
369  }
370  }
371  return 0;
372  }
373 
374  /** \brief Pick number of unique agent id's
375 
376  Unique from x and each other
377 
378  Uses the Fisher-Yates algorithm.
379 
380  \future AWS, 7/25/19: GSL may have a better Fisher-Yates
381  implementation we should use here. Or is it in the
382  \ref o2scl::permutation class?
383  */
384  virtual std::vector<int> pick_unique_agents(int nr, size_t x) {
385  std::vector<int> ids;
386  std::vector<int> agents;
387  // Fill array with ids
388  for (size_t i=0; i<pop_size-1; ++i){
389  if (i<x) {
390  ids.push_back( i );
391  } else {
392  ids.push_back( i+1 );
393  }
394  }
395  // Shuffle according to Fisher-Yates
396  for (size_t i=ids.size()-1; i>ids.size()-nr-1; --i) {
397  int j = round(gr.random()*i);
398  std::swap( ids[i], ids[j] );
399  }
400  for (size_t i=ids.size()-1; i>ids.size()-nr-1; --i) {
401  agents.push_back( ids[i] );
402  }
403  return agents;
404  }
405 
406 #endif
407 
408 #ifndef DOXYGEN_INTERNAL
409 
410  private:
411 
416 
417 #endif
418 
419  };
420 
421 #ifndef DOXYGEN_NO_O2NS
422 }
423 #endif
424 
425 #endif
o2scl::rng_gsl
Random number generator (GSL)
Definition: rng_gsl.h:55
o2scl::diff_evo::pop_size
size_t pop_size
Population size (default 0)
Definition: diff_evo.h:91
boost::numeric::ublas::vector< double >
o2scl::diff_evo::initialize_population
virtual int initialize_population(size_t nvar, vec_t &x0)
Initialize a population of random agents.
Definition: diff_evo.h:352
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::mmin_base< multi_funct, multi_funct, boost::numeric::ublas::vector< double > >::ntrial
int ntrial
Maximum number of iterations.
Definition: mmin.h:197
o2scl::mmin_base< multi_funct, multi_funct, boost::numeric::ublas::vector< double > >::last_ntrial
int last_ntrial
The number of iterations for in the most recent minimization.
Definition: mmin.h:206
o2scl::exc_emaxiter
@ exc_emaxiter
exceeded max number of iterations
Definition: err_hnd.h:73
o2scl::diff_evo::set_step
int set_step(size_t nv, vec2_t &stepv)
Set the step sizes for the default initialization.
Definition: diff_evo.h:317
o2scl::diff_evo::population
vec_t population
Vector containing the population.
Definition: diff_evo.h:333
o2scl::diff_evo::cr
double cr
Crossover probability (default 0.8)
Definition: diff_evo.h:109
o2scl::rng_gsl::random
double random()
Return a random number in .
Definition: rng_gsl.h:82
o2scl::diff_evo::pick_unique_agents
virtual std::vector< int > pick_unique_agents(int nr, size_t x)
Pick number of unique agent id's.
Definition: diff_evo.h:384
o2scl::multi_funct
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct
Multi-dimensional function typedef in src/base/multi_funct.h.
Definition: multi_funct.h:45
o2scl::diff_evo::step
std::vector< double > step
Step size for initialization.
Definition: diff_evo.h:348
o2scl::diff_evo::nconv
size_t nconv
The number of generations without a better fit before we assume that the algorithm has converged (def...
Definition: diff_evo.h:96
o2scl::diff_evo::print_iter
virtual void print_iter(size_t nvar, double fmin, int iter, vec_t &best_fit)
Print out iteration information.
Definition: diff_evo.h:292
o2scl::mmin_base< multi_funct, multi_funct, boost::numeric::ublas::vector< double > >::verbose
int verbose
Output control.
Definition: mmin.h:194
o2scl::diff_evo::f
double f
Differential weight (default 0.75)
Definition: diff_evo.h:103
O2SCL_CONV_RET
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:297
o2scl::itos
std::string itos(int x)
Convert an integer to a string.
o2scl::mm_funct
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct
Array of multi-dimensional functions typedef in src/base/mm_funct.h.
Definition: mm_funct.h:43
o2scl::diff_evo::rand_init_funct
init_funct_t * rand_init_funct
Function that is used to produce random init variables.
Definition: diff_evo.h:342
o2scl::diff_evo::fmins
ubvector fmins
Vector that keeps track of the function values.
Definition: diff_evo.h:336
o2scl::mmin_base
Multidimensional minimization [abstract base].
Definition: mmin.h:164
o2scl::diff_evo::gr
rng_gsl gr
Random number generator.
Definition: diff_evo.h:345
o2scl::diff_evo
Multidimensional minimization by the differential evolution method.
Definition: diff_evo.h:76
o2scl::diff_evo::set_init_function
virtual void set_init_function(init_funct_t &function)
Set the function that is used to select the initial population.
Definition: diff_evo.h:139
o2scl::diff_evo::mmin
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:165

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