anneal_gsl.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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 /* siman/siman.c
24  *
25  * Copyright (C) 2007 Brian Gough
26  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Mark Galassi
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 3 of the License, or (at
31  * your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful, but
34  * WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
36  * General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
41  * 02110-1301, USA.
42  */
43 #ifndef O2SCL_ANNEAL_GSL_H
44 #define O2SCL_ANNEAL_GSL_H
45 
46 /** \file anneal_gsl.h
47  \brief File defining \ref o2scl::anneal_gsl
48 */
49 #include <random>
50 
51 #include <boost/numeric/ublas/vector.hpp>
52 #include <boost/numeric/ublas/matrix.hpp>
53 
54 #include <o2scl/anneal.h>
55 #include <o2scl/multi_funct.h>
56 
57 #ifndef DOXYGEN_NO_O2NS
58 namespace o2scl {
59 #endif
60 
61  /** \brief Multidimensional minimization by simulated annealing (GSL)
62 
63  This class is a modification of simulated annealing as
64  implemented in GSL in the function \c gsl_siman_solve(). It acts
65  as a generic multidimensional minimizer for any function given a
66  generic temperature schedule specified by the user.
67 
68  There are a large variety of strategies for choosing the
69  temperature evolution. To offer the user the largest possible
70  flexibility, the temperature evolution is controlled by a
71  the virtual functions start() and next() which can be freely
72  changed by creating a child class which overwrites these
73  functions.
74 
75  The simulated annealing algorithm proposes a displacement of one
76  coordinate of the previous point by
77  \f[
78  x_{i,\mathrm{new}} = \mathrm{step\_size}_i
79  (2 u_i - 1) + x_{i,\mathrm{old}}
80  \f]
81  where the \f$u_i\f$ are random numbers between 0 and 1. The
82  displacement is accepted or rejected based on the Metropolis
83  method. The random number generator is set in the parent,
84  anneal.
85 
86  The default behavior is as follows: Initially, the step sizes
87  are chosen to be 1.0 (or whatever was recently specified in \ref
88  set_step() ) and the temperature to be \ref T_start (default
89  1.0). Each iteration decreases the temperature by a factor of
90  \ref T_dec (default 1.5) for each step, and the minimizer is
91  finished when the next decrease would bring the temperature
92  below \ref o2scl::mmin_base::tol_abs. If none of the
93  mmin_base::ntrial steps in a particular iteration changes the
94  value of the minimum, and the step sizes are greater than \ref
95  min_step_ratio (default 100) times \ref
96  o2scl::mmin_base::tol_abs, then the step sizes are decreased by
97  a factor of \ref step_dec (default 1.5) for the next iteration.
98 
99  If \ref o2scl::mmin_base::verbose is greater than zero, then
100  \ref mmin() will print out information and/or request a keypress
101  after the function iterations for each temperature.
102 
103  An example demonstrating the usage of this class is given in
104  <tt>examples/ex_anneal.cpp</tt> and in the \ref ex_anneal_sect .
105 
106  \comment
107  The form of the user-specified function is as in \ref
108  multi_funct has a "function value" which is the value of the
109  function (given in the third argument as a number of type \c
110  double), and a "return value" (the integer return value). The
111  initial function evaluation which is performed at the
112  user-specified initial guess must give 0 as the return value. If
113  subsequent function evaluations have a non-zero return value,
114  then the resulting point is ignored and a new point is selected.
115 
116  This class thus can sometimes handle constrained minimization
117  problems. If the user ensures that the function's return value
118  is non-zero when the function is evaluated outside the allowed
119  region, the minimizer will not accept any steps which take the
120  minimizer outside the allowed region. Note that this should be
121  done with care, however, as this approach may cause convergence
122  problems with sufficiently difficult functions or constraints.
123  \endcomment
124 
125  \comment
126  \future There's x0, old_x, new_x, and x? There's probably
127  some duplication here which could be avoided.
128  9/19/14: Some better documentation now, and it looks like
129  all four have some utility.
130  \endcomment
131 
132  \future Implement a more general simulated annealing routine
133  which would allow the solution of discrete problems like the
134  Traveling Salesman problem.
135  \comment
136  This could probably be done just by creating a parent abstract
137  base class which doesn't have any reference to step() and
138  step_vec.
139  \endcomment
140 
141  */
142  template<class func_t=multi_funct,
144  class rng_t=rng_gsl> class anneal_gsl :
145  public anneal_base<func_t,vec_t,rng_t> {
146 
147  public:
148 
150 
151  anneal_gsl() {
152  boltz=1.0;
153  step_vec.resize(1);
154  step_vec[0]=1.0;
155  T_start=1.0;
156  T_dec=1.5;
157  step_dec=1.5;
158  min_step_ratio=100.0;
159  }
160 
161  virtual ~anneal_gsl() {
162  }
163 
164  /** \brief Calculate the minimum \c fmin of \c func w.r.t the
165  array \c x0 of size \c nvar.
166  */
167  virtual int mmin(size_t nvar, vec_t &x0, double &fmin,
168  func_t &func) {
169 
170  if (nvar==0) {
171  O2SCL_ERR2("Tried to minimize over zero variables ",
172  " in anneal_gsl::mmin().",exc_einval);
173  }
174 
175  fmin=0.0;
176 
177  step_norm=1.0;
178 
179  double E, new_E, T, old_E;
180  int i, iter=0;
181  size_t j;
182 
183  // Current point
184  vec_t x(nvar);
185  // Proposed next point
186  vec_t new_x(nvar);
187  // Last point from previous iteration
188  vec_t old_x(nvar);
189 
190  for(j=0;j<nvar;j++) {
191  x[j]=x0[j];
192  }
193 
194  E=func(nvar,x);
195 
196  // We use x0 and fmin to store the best point found
197  fmin=E;
198 
199  // Setup initial temperature and step sizes
200  start(nvar,T);
201 
202  bool done=false;
203 
204  while (!done) {
205 
206  // Copy old value of x for next() function
207  for(j=0;j<nvar;j++) old_x[j]=x[j];
208  old_E=E;
209 
210  size_t nmoves=0;
211 
212  for (i=0;i<this->ntrial;++i) {
213  for (j=0;j<nvar;j++) new_x[j]=x[j];
214 
215  step(new_x,nvar);
216  new_E=func(nvar,new_x);
217 
218  // Store best value obtained so far
219  if(new_E<=fmin) {
220  for(j=0;j<nvar;j++) x0[j]=new_x[j];
221  fmin=new_E;
222  }
223 
224  // Take the crucial step: see if the new point is accepted
225  // or not, as determined by the Boltzmann probability
226  if (new_E<E) {
227  for(j=0;j<nvar;j++) x[j]=new_x[j];
228  E=new_E;
229  nmoves++;
230  } else {
231  double r=this->dist();
232  if (r < exp(-(new_E-E)/(boltz*T))) {
233  for(j=0;j<nvar;j++) x[j]=new_x[j];
234  E=new_E;
235  nmoves++;
236  }
237  }
238 
239  }
240 
241  if (this->verbose>0) {
242  this->print_iter(nvar,x0,fmin,iter,T,"anneal_gsl");
243  iter++;
244  if (this->verbose>1) {
245  char ch;
246  std::cin >> ch;
247  }
248  }
249 
250  // See if we're finished and proceed to the next step
251  next(nvar,old_x,old_E,x,E,T,nmoves,x0,fmin,done);
252 
253  }
254 
255  return 0;
256  }
257 
258  /// Return string denoting type ("anneal_gsl")
259  virtual const char *type() { return "anneal_gsl"; }
260 
261  /** \brief Boltzmann factor (default 1.0).
262  */
263  double boltz;
264 
265  /// Set the step sizes
266  template<class vec2_t> int set_step(size_t nv, vec2_t &stepv) {
267  if (nv>0) {
268  step_vec.resize(nv);
269  for(size_t i=0;i<nv;i++) step_vec[i]=stepv[i];
270  }
271  return 0;
272  }
273 
274  /// Initial temperature (default 1.0)
275  double T_start;
276 
277  /// Factor to decrease temperature by (default 1.5)
278  double T_dec;
279 
280  /// Factor to decrease step size by (default 1.5)
281  double step_dec;
282 
283  /// Ratio between minimum step size and \ref tol_abs (default 100.0)
285 
286  /** \brief Copy constructor
287  */
291 
292  boltz=ag.boltz;
293  T_start=ag.T_start;
294  T_dec=ag.T_dec;
295  step_dec=ag.step_dec;
297  step_vec=ag.step_vec;
298  }
299 
300  /** \brief Copy constructor from operator=
301  */
304  if (this != &ag) {
306  boltz=ag.boltz;
307  T_start=ag.T_start;
308  T_dec=ag.T_dec;
309  step_dec=ag.step_dec;
311  step_vec=ag.step_vec;
312  }
313  return *this;
314  }
315 
316 #ifndef DOXYGEN_INTERNAL
317 
318  protected:
319 
320  /// \name Storage for points in parameter space
321  //@{
322  //@}
323 
324  /// Normalization for step
325  double step_norm;
326 
327  /// Vector of step sizes
329 
330  /// Determine how to change the minimization for the next iteration
331  virtual int next(size_t nvar, vec_t &x_old, double min_old,
332  vec_t &x_new, double min_new, double &T,
333  size_t n_moves, vec_t &best_x, double best_E,
334  bool &finished) {
335 
336  if (T/T_dec<this->tol_abs) {
337  finished=true;
338  return success;
339  }
340  if (n_moves==0) {
341  // If we haven't made any progress, shrink the step by
342  // decreasing step_norm
344  // Also reset x to best value so far
345  for(size_t i=0;i<nvar;i++) {
346  x_new[i]=best_x[i];
347  }
348  min_new=best_E;
349  }
350  T/=T_dec;
351  return success;
352  }
353 
354  /// Setup initial temperature and stepsize
355  virtual int start(size_t nvar, double &T) {
356  T=T_start;
357  return success;
358  }
359 
360  /** \brief Make a step to a new attempted minimum
361  */
362  virtual int step(vec_t &sx, int nvar) {
363  size_t nstep=step_vec.size();
364  for(int i=0;i<nvar;i++) {
365  double u=this->dist();
366 
367  // Construct the step in the ith direction
368  double step_i=step_norm*step_vec[i%nstep];
369  // Fix if the step is too small
370  if (step_i<this->tol_abs*min_step_ratio) {
371  step_i=this->tol_abs*min_step_ratio;
372  }
373 
374  sx[i]=(2.0*u-1.0)*step_i+sx[i];
375  }
376  return 0;
377  }
378 
379 #endif
380 
381  };
382 
383 #ifndef DOXYGEN_NO_O2NS
384 }
385 #endif
386 
387 #endif
o2scl::anneal_gsl::set_step
int set_step(size_t nv, vec2_t &stepv)
Set the step sizes.
Definition: anneal_gsl.h:266
o2scl::rng_gsl
Random number generator (GSL)
Definition: rng_gsl.h:55
o2scl::anneal_base::operator=
anneal_base< func_t, vec_t, rng_t > & operator=(const anneal_base< func_t, vec_t, rng_t > &ab)
Copy constructor from operator=.
Definition: anneal.h:142
boost::numeric::ublas::vector< double >
o2scl::anneal_gsl::type
virtual const char * type()
Return string denoting type ("anneal_gsl")
Definition: anneal_gsl.h:259
o2scl::anneal_base< multi_funct, boost::numeric::ublas::vector< double >, rng_gsl >::dist
o2scl::prob_dens_uniform dist
The random distribution object.
Definition: anneal.h:123
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
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::anneal_gsl::min_step_ratio
double min_step_ratio
Ratio between minimum step size and tol_abs (default 100.0)
Definition: anneal_gsl.h:284
o2scl::anneal_gsl::T_dec
double T_dec
Factor to decrease temperature by (default 1.5)
Definition: anneal_gsl.h:278
o2scl::mmin_base< multi_funct, multi_funct, boost::numeric::ublas::vector< double > >::ntrial
int ntrial
Maximum number of iterations.
Definition: mmin.h:197
o2scl::anneal_gsl
Multidimensional minimization by simulated annealing (GSL)
Definition: anneal_gsl.h:144
o2scl::anneal_base
Simulated annealing base.
Definition: anneal.h:67
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::success
@ success
Success.
Definition: err_hnd.h:47
o2scl::anneal_gsl::step_vec
ubvector step_vec
Vector of step sizes.
Definition: anneal_gsl.h:328
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::mmin_base< multi_funct, multi_funct, boost::numeric::ublas::vector< double > >::verbose
int verbose
Output control.
Definition: mmin.h:194
o2scl::anneal_gsl::boltz
double boltz
Boltzmann factor (default 1.0).
Definition: anneal_gsl.h:263
o2scl::anneal_gsl::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 x0 of size nvar.
Definition: anneal_gsl.h:167
o2scl::anneal_gsl::next
virtual int next(size_t nvar, vec_t &x_old, double min_old, vec_t &x_new, double min_new, double &T, size_t n_moves, vec_t &best_x, double best_E, bool &finished)
Determine how to change the minimization for the next iteration.
Definition: anneal_gsl.h:331
o2scl::anneal_gsl::step_dec
double step_dec
Factor to decrease step size by (default 1.5)
Definition: anneal_gsl.h:281
o2scl::anneal_base< multi_funct, boost::numeric::ublas::vector< double >, rng_gsl >::print_iter
virtual int print_iter(size_t nv, boost::numeric::ublas::vector< double > &x, double y, int iter, double tptr, std::string comment)
Print out iteration information.
Definition: anneal.h:98
o2scl::anneal_gsl::step
virtual int step(vec_t &sx, int nvar)
Make a step to a new attempted minimum.
Definition: anneal_gsl.h:362
o2scl::anneal_gsl::start
virtual int start(size_t nvar, double &T)
Setup initial temperature and stepsize.
Definition: anneal_gsl.h:355
o2scl::anneal_gsl::step_norm
double step_norm
Normalization for step.
Definition: anneal_gsl.h:325
o2scl::mmin_base< multi_funct, multi_funct, boost::numeric::ublas::vector< double > >::tol_abs
double tol_abs
The independent variable tolerance.
Definition: mmin.h:203
o2scl::anneal_gsl::T_start
double T_start
Initial temperature (default 1.0)
Definition: anneal_gsl.h:275

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