anneal_gsl.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, 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_funct11 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  See also a multi-threaded version of this class in \ref
126  anneal_mt.
127 
128  \comment
129  \future There's x0, old_x, new_x, best_x, and x? There's probably
130  some duplication here which could be avoided.
131  9/19/14: Some better documentation now, and it looks like
132  all four have some utility.
133  \endcomment
134 
135  \future Implement a more general simulated annealing routine
136  which would allow the solution of discrete problems like the
137  Traveling Salesman problem.
138  \comment
139  This could probably be done just by creating a parent abstract
140  base class which doesn't have any reference to step() and
141  step_vec.
142  \endcomment
143 
144  */
145  template<class func_t=multi_funct11,
147  class rng_t=rng_gsl> class anneal_gsl :
148  public anneal_base<func_t,vec_t,rng_t> {
149 
150  public:
151 
153 
154  anneal_gsl() {
155  boltz=1.0;
156  step_vec.resize(1);
157  step_vec[0]=1.0;
158  T_start=1.0;
159  T_dec=1.5;
160  step_dec=1.5;
161  min_step_ratio=100.0;
162  }
163 
164  virtual ~anneal_gsl() {
165  }
166 
167  /** \brief Calculate the minimum \c fmin of \c func w.r.t the
168  array \c x0 of size \c nvar.
169  */
170  virtual int mmin(size_t nvar, vec_t &x0, double &fmin,
171  func_t &func) {
172 
173  if (nvar==0) {
174  O2SCL_ERR2("Tried to minimize over zero variables ",
175  " in anneal_gsl::mmin().",exc_einval);
176  }
177 
178  fmin=0.0;
179 
180  step_norm=1.0;
181 
182  allocate(nvar);
183 
184  double E, new_E, best_E, T, old_E;
185  int i, iter=0;
186  size_t j;
187 
188  for(j=0;j<nvar;j++) {
189  x[j]=x0[j];
190  best_x[j]=x0[j];
191  }
192 
193  E=func(nvar,x);
194  best_E=E;
195 
196  // Setup initial temperature and step sizes
197  start(nvar,T);
198 
199  bool done=false;
200 
201  while (!done) {
202 
203  // Copy old value of x for next() function
204  for(j=0;j<nvar;j++) old_x[j]=x[j];
205  old_E=E;
206 
207  size_t nmoves=0;
208 
209  for (i=0;i<this->ntrial;++i) {
210  for (j=0;j<nvar;j++) new_x[j]=x[j];
211 
212  step(new_x,nvar);
213  new_E=func(nvar,new_x);
214 
215  // Store best value obtained so far
216  if(new_E<=best_E){
217  for(j=0;j<nvar;j++) best_x[j]=new_x[j];
218  best_E=new_E;
219  }
220 
221  // Take the crucial step: see if the new point is accepted
222  // or not, as determined by the Boltzmann probability
223  if (new_E<E) {
224  for(j=0;j<nvar;j++) x[j]=new_x[j];
225  E=new_E;
226  nmoves++;
227  } else {
228  double r=this->dist();
229  if (r < exp(-(new_E-E)/(boltz*T))) {
230  for(j=0;j<nvar;j++) x[j]=new_x[j];
231  E=new_E;
232  nmoves++;
233  }
234  }
235 
236  }
237 
238  if (this->verbose>0) {
239  this->print_iter(nvar,best_x,best_E,iter,T,"anneal_gsl");
240  iter++;
241  }
242 
243  // See if we're finished and proceed to the next step
244  next(nvar,old_x,old_E,x,E,T,nmoves,best_E,done);
245 
246  }
247 
248  for(j=0;j<nvar;j++) x0[j]=best_x[j];
249  fmin=best_E;
250 
251  return 0;
252  }
253 
254  /// Return string denoting type ("anneal_gsl")
255  virtual const char *type() { return "anneal_gsl"; }
256 
257  /** \brief Boltzmann factor (default 1.0).
258  */
259  double boltz;
260 
261  /// Set the step sizes
262  template<class vec2_t> int set_step(size_t nv, vec2_t &stepv) {
263  if (nv>0) {
264  step_vec.resize(nv);
265  for(size_t i=0;i<nv;i++) step_vec[i]=stepv[i];
266  }
267  return 0;
268  }
269 
270  /// Initial temperature (default 1.0)
271  double T_start;
272 
273  /// Factor to decrease temperature by (default 1.5)
274  double T_dec;
275 
276  /// Factor to decrease step size by (default 1.5)
277  double step_dec;
278 
279  /// Ratio between minimum step size and \ref tol_abs (default 100.0)
281 
282  /** \brief Copy constructor
283  */
287 
288  boltz=ag.boltz;
289  T_start=ag.T_start;
290  T_dec=ag.T_dec;
291  step_dec=ag.step_dec;
292  min_step_ratio=ag.min_step_ratio;
293  step_vec=ag.step_vec;
294  x=ag.x;
295  new_x=ag.new_x;
296  best_x=ag.best_x;
297  old_x=ag.old_x;
298 
299  }
300 
301  /** \brief Copy constructor from operator=
302  */
305  if (this != &ag) {
307  boltz=ag.boltz;
308  T_start=ag.T_start;
309  T_dec=ag.T_dec;
310  step_dec=ag.step_dec;
311  min_step_ratio=ag.min_step_ratio;
312  step_vec=ag.step_vec;
313  x=ag.x;
314  new_x=ag.new_x;
315  best_x=ag.best_x;
316  old_x=ag.old_x;
317  }
318  return *this;
319  }
320 
321 #ifndef DOXYGEN_INTERNAL
322 
323  protected:
324 
325  /// \name Storage for points in parameter space
326  //@{
327  /// Current point
328  ubvector x;
329  /// Proposed next point
330  ubvector new_x;
331  /// Optimum point over all iterations
332  ubvector best_x;
333  /// Last point from previous iteration
334  ubvector old_x;
335  //@}
336 
337  /// Normalization for step
338  double step_norm;
339 
340  /// Vector of step sizes
341  ubvector step_vec;
342 
343  /// Determine how to change the minimization for the next iteration
344  virtual int next(size_t nvar, vec_t &x_old, double min_old,
345  vec_t &x_new, double min_new, double &T,
346  size_t n_moves, double best_E, bool &finished) {
347 
348  if (T/T_dec<this->tol_abs) {
349  finished=true;
350  return success;
351  }
352  if (n_moves==0) {
353  // If we haven't made any progress, shrink the step by
354  // decreasing step_norm
355  step_norm/=step_dec;
356  // Also reset x to best value so far
357  for(size_t i=0;i<nvar;i++) {
358  x_new[i]=best_x[i];
359  }
360  min_new=best_E;
361  }
362  T/=T_dec;
363  return success;
364  }
365 
366  /// Setup initial temperature and stepsize
367  virtual int start(size_t nvar, double &T) {
368  T=T_start;
369  return success;
370  }
371 
372  /** \brief Allocate memory for a minimizer over \c n dimensions
373  with stepsize \c step and Boltzmann factor \c boltz_factor
374  */
375  virtual int allocate(size_t n, double boltz_factor=1.0) {
376  x.resize(n);
377  old_x.resize(n);
378  new_x.resize(n);
379  best_x.resize(n);
380  boltz=boltz_factor;
381  return 0;
382  }
383 
384  /** \brief Make a step to a new attempted minimum
385  */
386  virtual int step(vec_t &sx, int nvar) {
387  size_t nstep=step_vec.size();
388  for(int i=0;i<nvar;i++) {
389  double u=this->dist();
390 
391  // Construct the step in the ith direction
392  double step_i=step_norm*step_vec[i%nstep];
393  // Fix if the step is too small
394  if (step_i<this->tol_abs*min_step_ratio) {
395  step_i=this->tol_abs*min_step_ratio;
396  }
397 
398  sx[i]=(2.0*u-1.0)*step_i+sx[i];
399  }
400  return 0;
401  }
402 
403 #endif
404 
405  };
406 
407 #ifndef DOXYGEN_NO_O2NS
408 }
409 #endif
410 
411 #endif
int set_step(size_t nv, vec2_t &stepv)
Set the step sizes.
Definition: anneal_gsl.h:262
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
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, double best_E, bool &finished)
Determine how to change the minimization for the next iteration.
Definition: anneal_gsl.h:344
ubvector best_x
Optimum point over all iterations.
Definition: anneal_gsl.h:332
ubvector step_vec
Vector of step sizes.
Definition: anneal_gsl.h:341
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::prob_dens_uniform dist
The random distribution object.
Definition: anneal.h:123
ubvector new_x
Proposed next point.
Definition: anneal_gsl.h:330
double T_dec
Factor to decrease temperature by (default 1.5)
Definition: anneal_gsl.h:274
virtual int allocate(size_t n, double boltz_factor=1.0)
Allocate memory for a minimizer over n dimensions with stepsize step and Boltzmann factor boltz_facto...
Definition: anneal_gsl.h:375
double boltz
Boltzmann factor (default 1.0).
Definition: anneal_gsl.h:259
Simulated annealing base.
Definition: anneal.h:67
double T_start
Initial temperature (default 1.0)
Definition: anneal_gsl.h:271
ubvector x
Current point.
Definition: anneal_gsl.h:328
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
virtual int start(size_t nvar, double &T)
Setup initial temperature and stepsize.
Definition: anneal_gsl.h:367
double min_step_ratio
Ratio between minimum step size and tol_abs (default 100.0)
Definition: anneal_gsl.h:280
ubvector old_x
Last point from previous iteration.
Definition: anneal_gsl.h:334
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:170
double tol_abs
The independent variable tolerance.
Definition: mmin.h:203
virtual int step(vec_t &sx, int nvar)
Make a step to a new attempted minimum.
Definition: anneal_gsl.h:386
Multidimensional minimization by simulated annealing (GSL)
Definition: anneal_gsl.h:147
Random number generator (GSL)
Definition: rng_gsl.h:55
virtual int print_iter(size_t nv, vec_t &x, double y, int iter, double tptr, std::string comment)
Print out iteration information.
Definition: anneal.h:98
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct11
Multi-dimensional function typedef.
Definition: multi_funct.h:45
virtual const char * type()
Return string denoting type ("anneal_gsl")
Definition: anneal_gsl.h:255
double step_dec
Factor to decrease step size by (default 1.5)
Definition: anneal_gsl.h:277
Success.
Definition: err_hnd.h:47
double step_norm
Normalization for step.
Definition: anneal_gsl.h:338
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
int ntrial
Maximum number of iterations.
Definition: mmin.h:197

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