min.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 #ifndef O2SCL_MINIMIZE_H
24 #define O2SCL_MINIMIZE_H
25 
26 // For fabs()
27 #include <cmath>
28 
29 #include <o2scl/err_hnd.h>
30 #include <o2scl/funct.h>
31 
32 /** \file min.h
33  \brief One-dimensional minimization base class and associated functions
34 */
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief One-dimensional minimization [abstract base]
41  */
42  template<class func_t=funct, class dfunc_t=func_t> class min_base {
43 
44  public:
45 
46  min_base() {
47  verbose=0;
48  ntrial=100;
49  tol_rel=1.0e-4;
50  tol_abs=1.0e-4;
51  last_ntrial=0;
52  bracket_iter=20;
53  err_nonconv=true;
54  }
55 
56  virtual ~min_base() {}
57 
58  /// Output control
59  int verbose;
60 
61  /// Maximum number of iterations
62  int ntrial;
63 
64  /// The tolerance for the minimum function value
65  double tol_rel;
66 
67  /// The tolerance for the location of the minimum
68  double tol_abs;
69 
70  /// The number of iterations used in the most recent minimization
72 
73  /** \brief The number of iterations for automatically
74  bracketing a minimum (default 20)
75  */
77 
78  /// If true, call the error handler if the routine does not "converge"
80 
81  /** \brief Print out iteration information.
82 
83  Depending on the value of the variable \ref verbose, this
84  prints out the iteration information. If verbose=0, then no
85  information is printed, while if verbose>1, then after each
86  iteration, the present values of x and y are output to
87  std::cout along with the iteration number. If verbose>=2 then
88  each iteration waits for a character.
89  */
90  virtual int print_iter(double x, double y, int iter, double value=0.0,
91  double limit=0.0, std::string comment="") {
92 
93  if (verbose<=0) return success;
94 
95  char ch;
96 
97  std::cout << comment << " Iteration: " << iter << std::endl;
98  if (x<0) std::cout << x << " ";
99  else std::cout << " " << x << " ";
100  if (y<0) std::cout << y << " ";
101  else std::cout << " " << y << " ";
102  if (value<0) std::cout << value << " ";
103  else std::cout << " " << value << " ";
104  if (limit<0) std::cout << limit << std::endl;
105  else std::cout << " " << limit << std::endl;
106  if (verbose>1) {
107  std::cout << "Press a key and type enter to continue. ";
108  std::cin >> ch;
109  }
110 
111  return success;
112  }
113 
114  /** \brief Calculate the minimum \c min of \c func w.r.t 'x'.
115 
116  If this is not overloaded, it attempts to bracket the
117  minimum using bracket() and then calls min_bkt() with
118  the newly bracketed minimum.
119  */
120  virtual int min(double &x, double &fmin, func_t &func)=0;
121 
122  /** \brief Calculate the minimum \c min of \c func with x2
123  bracketed between x1 and x3
124 
125  If this is not overloaded, it ignores the bracket and calls min().
126  */
127  virtual int min_bkt(double &x2, double x1, double x3, double &fmin,
128  func_t &func)=0;
129 
130  /** \brief Calculate the minimum \c min of \c func with
131  derivative \c dfunc w.r.t 'x'.
132 
133  If this is not overloaded, it attempts to bracket the
134  minimum using bracket() and then calls min_bkt_de() with
135  the newly bracketed minimum.
136  */
137  virtual int min_de(double &x, double &fmin, func_t &func,
138  dfunc_t &df)=0;
139 
140  /** \brief Given interval <tt>(ax,bx)</tt>, attempt to bracket a
141  minimum for function \c func.
142 
143  Upon success, <tt>fa=func(ax)</tt>, <tt>fb=func(bx)</tt>, and
144  <tt>fc=func(cx)</tt> with <tt>fb<fa</tt>, <tt>fb<fc</tt> and
145  <tt>ax<bx<cx</tt>. The initial values of \c cx, \c fa,
146  \c fb, and \c fc are all ignored.
147 
148  The number of iterations is controlled by \ref bracket_iter.
149 
150  \note This algorithm can fail if there's a minimum which has a
151  much smaller size than \f$ bx-ax \f$, or if the function has the
152  same value at \c ax, \c bx, and the midpoint <tt>(ax+bx)/2</tt>.
153 
154  \future Improve this algorithm with the golden ratio
155  method in gsl/min/bracketing.c?
156  */
157  virtual int bracket(double &ax, double &bx, double &cx, double &fa,
158  double &fb, double &fc, func_t &func) {
159 
160  double x=ax, x2=bx, x3=(ax+bx)/2.0;
161  double fx, fx3, fx2;
162  int i=0;
163 
164  bool done=false;
165  while(done==false && i<bracket_iter) {
166  fx=func(x);
167  fx2=func(x2);
168  fx3=func(x3);
169 
170  if (verbose>0) {
171  std::cout << "Function min::bracket(), Iteration: "
172  << i << std::endl;
173  std::cout << " " << x << " " << x3 << " " << x2 << std::endl;
174  std::cout << " " << fx << " " << fx3 << " " << fx2 << std::endl;
175  if (verbose>1) {
176  char ch;
177  std::cout << "Press a key and type enter to continue. ";
178  std::cin >> ch;
179  }
180  }
181 
182  if (fx3>=fx2 && fx3>=fx) {
183  // If the middle value is higher than the endpoints,
184  // try again with one side or the other
185  if (fx2>fx) {
186  x2=x3;
187  x3=(x+x2)/2.0;
188  } else {
189  x=x3;
190  x3=(x+x2)/2.0;
191  }
192  } else if (fx<=fx3 && fx3<=fx2) {
193  // If we're increasing, extend to the left
194  x3=x;
195  x=x2-2.0*(x2-x);
196  } else if (fx3<fx2 && fx3<fx) {
197  // If we've succeeded, copy the results over
198  done=true;
199  ax=x;
200  bx=x3;
201  cx=x2;
202  fa=fx;
203  fb=fx3;
204  fc=fx2;
205  } else {
206  // Otherwise we're decreasing, extend to the right
207  x3=x2;
208  x2=x+2.0*(x2-x);
209  }
210  i++;
211  }
212 
213  if (done==false) {
214  O2SCL_ERR("Too many iterations in min::bracket().",
215  exc_emaxiter);
216  }
217 
218  return 0;
219  }
220 
221  /// Return string denoting type ("min")
222  virtual const char *type() { return "min"; }
223 
224  };
225 
226  /** \brief One-dimensional bracketing minimization [abstract base]
227  */
228  template<class func_t, class dfunc_t=func_t>
229  class min_bkt_base : public min_base<func_t,dfunc_t> {
230 
231  public:
232 
233  min_bkt_base() {
234  bracket_iter=20;
235  }
236 
237  virtual ~min_bkt_base() {}
238 
239  /** \brief The number of iterations for automatically
240  bracketing a minimum (default 20)
241  */
243 
244  /** \brief Calculate the minimum \c min of \c func w.r.t 'x'.
245 
246  If this is not overloaded, it attempts to bracket the
247  minimum using bracket() and then calls min_bkt() with
248  the newly bracketed minimum.
249  */
250  virtual int min(double &x, double &fmin, func_t &func) {
251  double xl, xr, f, fl, fr;
252  xl=x*0.9;
253  xr=x*1.1;
254  if (this->bracket(xl,x,xr,fl,f,fr,func)!=0) {
255  O2SCL_CONV_RET("Failed to bracket in min().",exc_efailed,
256  this->err_nonconv);
257  }
258  return min_bkt(x,xl,xr,fmin,func);
259  }
260 
261  /** \brief Calculate the minimum \c min of \c func with x2
262  bracketed between x1 and x3
263 
264  If this is not overloaded, it ignores the bracket and calls min().
265  */
266  virtual int min_bkt(double &x2, double x1, double x3, double &fmin,
267  func_t &func)=0;
268 
269  /** \brief Calculate the minimum \c min of \c func with
270  derivative \c dfunc w.r.t 'x'.
271 
272  If this is not overloaded, it attempts to bracket the
273  minimum using bracket() and then calls min_bkt_de() with
274  the newly bracketed minimum.
275  */
276  virtual int min_de(double &x, double &fmin, func_t &func,
277  dfunc_t &df) {
278  double xl, xr, f, fl, fr;
279  xl=x*0.9;
280  xr=x*1.1;
281  if (this->bracket(xl,x,xr,fl,f,fr,func)!=0) {
282  O2SCL_CONV_RET("Failed to bracket in min_de().",exc_efailed,
283  this->err_nonconv);
284  }
285  return min_bkt(x,xl,xr,fmin,func);
286  }
287 
288  /// Return string denoting type ("min_bkt")
289  virtual const char *type() { return "min_bkt"; }
290 
291  };
292 
293  /** \brief One-dimensional minimization using derivatives [abstract base]
294 
295  At the moment there are no minimizers of this type implemented in
296  \o2.
297 
298  \future Create a version of \ref o2scl::mmin_conf which
299  implements a minimizer with this interface.
300  */
301  template<class func_t, class dfunc_t=func_t>
302  class min_de_base : public min_base<func_t,dfunc_t> {
303 
304  public:
305 
306  min_de_base() {
307  }
308 
309  virtual ~min_de_base() {}
310 
311  /** \brief Calculate the minimum \c min of \c func w.r.t 'x'.
312 
313  If this is not overloaded, it attempts to bracket the
314  minimum using bracket() and then calls min_bkt() with
315  the newly bracketed minimum.
316  */
317  virtual int min(double &x, double &fmin, func_t &func)=0;
318 
319  /** \brief Calculate the minimum \c min of \c func with x2
320  bracketed between x1 and x3
321 
322  If this is not overloaded, it ignores the bracket and calls min().
323  */
324  virtual int min_bkt(double &x2, double x1, double x3, double &fmin,
325  func_t &func)=0;
326 
327  /** \brief Calculate the minimum \c min of \c func with
328  derivative \c dfunc w.r.t 'x'.
329 
330  If this is not overloaded, it attempts to bracket the
331  minimum using bracket() and then calls min_bkt_de() with
332  the newly bracketed minimum.
333  */
334  virtual int min_de(double &x, double &fmin, func_t &func,
335  dfunc_t &df)=0;
336 
337  /// Return string denoting type ("min_de")
338  virtual const char *type() { return "min_de"; }
339 
340  };
341 
342  /// \name Helper functions for constrained min in src/min/min.h
343  //@{
344  /** \brief Constrain \c x to be within \c width
345  of the value given by \c center
346 
347  Defining \f$ c= \f$ \c center, \f$ w= \f$ \c width, \f$ h= \f$
348  \c height, this returns the value \f$ h (1+|x-c-w|/w) \f$ if \f$
349  x>c+w \f$ and \f$ h (1+|x-c+w|/w) \f$ if \f$ x<c-w \f$ . The
350  value near \f$ x=c-w \f$ or \f$ x=c+w \f$ is \f$ h \f$ (the
351  value of the function exactly at these points is zero) and the
352  value at \f$ x=c-2w \f$ or \f$ x=c+2w \f$ is \f$ 2 h \f$ .
353 
354  It is important to note that, for large distances of \c x
355  from \c center, this only scales linearly. If you are trying to
356  constrain a function which decreases more than linearly by
357  making \c x far from \c center, then a minimizer may
358  ignore this constraint.
359  */
360  inline double constraint(double x, double center, double width,
361  double height) {
362  double ret=0.0;
363  if (x>center+width) {
364  ret=height*(1.0+fabs(x-center-width)/width);
365  } else if (x<center-width) {
366  ret=height*(1.0+fabs(x-center+width)/width);
367  }
368  return ret;
369  }
370 
371  /** \brief Constrain \c x to be within \c width of the value given
372  by \c center
373 
374  Defining \f$ c= \f$ \c center, \f$ w= \f$ \c width, \f$ h= \f$
375  \c height, \f$ t= \f$ \c tightness, and \f$ \ell= \f$ \c
376  exp_arg_limit, this returns the value
377  \f[
378  h \left(\frac{x-c}{w}\right)^2 \left[
379  1+ e^{t\left(x-c+w\right)\left(c+w-x\right)/w^2}
380  \right]^{-1}
381  \f]
382 
383  This function is continuous and differentiable. Note that if
384  \f$ x=c \f$ , then the function returns zero.
385 
386  The exponential is handled gracefully by assuming that anything
387  smaller than \f$ \exp(-\ell) \f$ is zero. This creates a small
388  discontinuity which can be removed with the sufficiently large
389  value of \f$ \ell \f$.
390 
391  It is important to note that, for large distances of \c x from
392  \c center, this scales quadratically. If you are trying to
393  constrain a function which decreases faster than quadratically
394  by making \c x far from \c center, then a minimizer may ignore
395  this constraint.
396 
397  In the limit \f$ t \rightarrow \infty \f$, this function
398  converges towards the squared value of \ref constraint(), except
399  exactly at the points \f$ x=c-w \f$ and \f$ x=c+w \f$.
400  */
401  inline double cont_constraint(double x, double center, double width,
402  double height, double tightness=40.0,
403  double exp_arg_limit=50.0) {
404  double ret, wid2=width*width;
405  double arg=tightness/wid2*(x-center+width)*(center+width-x);
406  if (arg<-exp_arg_limit) {
407  ret=(x-center)*(x-center)/wid2;
408  } else {
409  ret=(x-center)*(x-center)/wid2/(1.0+exp(arg));
410  }
411  return ret*height;
412  }
413 
414  /** \brief Constrain \c x to be greater than the value given by \c
415  center
416 
417  Defining \f$ c= \f$ \c center, \f$ w= \f$ \c width, \f$ h= \f$
418  \c height, this returns \f$ h(1+|x-c|/w) \f$ if \f$ x<c \f$ and
419  zero otherwise. The value at \f$ x=c \f$ is \f$ h \f$ , while
420  the value at \f$ x=c-w \f$ is \f$ 2 h \f$ .
421 
422  It is important to note that, for large distances of \c x from
423  \c center, this only scales linearly. If you are trying to
424  constrain a function which decreases more than linearly by
425  making \c x far from \c center, then a minimizer may ignore this
426  constraint.
427  */
428  inline double lower_bound(double x, double center, double width,
429  double height) {
430  double ret=0.0;
431  if (x<center) ret=height*(1.0+fabs(x-center)/width);
432  return ret;
433  }
434 
435  /** \brief Constrain \c x to be greater than the value given by \c
436  center
437 
438  Defining \f$ c= \f$ \c center, \f$ w= \f$ \c width, \f$ h= \f$
439  \c height, \f$ t= \f$ \c tightness, and \f$ \ell= \f$ \c
440  exp_arg_limit, this returns \f$ h(c-x+w)/(w+w\exp(t(x-c)/w)) \f$
441  and has the advantage of being a continuous and differentiable
442  function. The value of the function exactly at \f$ x=c \f$ is
443  \f$ h/2 \f$, but for \f$ x \f$ just below \f$ c \f$ the function
444  is \f$ h \f$ and just above \f$ c \f$ the function is quite
445  small.
446 
447  The exponential is handled gracefully by assuming that anything
448  smaller than \f$ \exp(-\ell) \f$ is zero. This creates a small
449  discontinuity which can be removed with the sufficiently large
450  value of \f$ \ell \f$.
451 
452  It is important to note that, for large distances of \c x
453  from \c center, this only scales linearly. If you are trying to
454  constrain a function which decreases more than linearly by
455  making \c x far from \c center, then a minimizer may
456  ingore this constraint.
457 
458  In the limit \f$ t \rightarrow \infty \f$, this function
459  converges towards \ref lower_bound(), except exactly at the
460  point \f$ x=c \f$.
461  */
462  inline double cont_lower_bound(double x, double center, double width,
463  double height, double tightness=40.0,
464  double exp_arg_limit=50.0) {
465  double ret, arg=tightness*(x-center)/width;
466  if (arg>exp_arg_limit) {
467  ret=0.0;
468  } else if (arg<-exp_arg_limit) {
469  ret=height*(center-x+width)/width;
470  } else {
471  ret=height*(center-x+width)/width/(1.0+exp(arg));
472  }
473  return ret;
474  }
475  //@}
476 
477 #ifndef DOXYGEN_NO_O2NS
478 }
479 #endif
480 
481 #endif
o2scl::min_bkt_base::min_de
virtual int min_de(double &x, double &fmin, func_t &func, dfunc_t &df)
Calculate the minimum min of func with derivative dfunc w.r.t 'x'.
Definition: min.h:276
o2scl::min_de_base::type
virtual const char * type()
Return string denoting type ("min_de")
Definition: min.h:338
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl::min_base::tol_abs
double tol_abs
The tolerance for the location of the minimum.
Definition: min.h:68
o2scl::min_bkt_base::bracket_iter
int bracket_iter
The number of iterations for automatically bracketing a minimum (default 20)
Definition: min.h:242
o2scl::min_base::type
virtual const char * type()
Return string denoting type ("min")
Definition: min.h:222
o2scl::min_base::ntrial
int ntrial
Maximum number of iterations.
Definition: min.h:62
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::min_base::min_de
virtual int min_de(double &x, double &fmin, func_t &func, dfunc_t &df)=0
Calculate the minimum min of func with derivative dfunc w.r.t 'x'.
o2scl::exc_emaxiter
@ exc_emaxiter
exceeded max number of iterations
Definition: err_hnd.h:73
o2scl::min_base::bracket_iter
int bracket_iter
The number of iterations for automatically bracketing a minimum (default 20)
Definition: min.h:76
o2scl::min_bkt_base::min_bkt
virtual int min_bkt(double &x2, double x1, double x3, double &fmin, func_t &func)=0
Calculate the minimum min of func with x2 bracketed between x1 and x3.
o2scl_inte_qng_coeffs::x2
static const double x2[5]
Definition: inte_qng_gsl.h:66
o2scl::min_bkt_base::type
virtual const char * type()
Return string denoting type ("min_bkt")
Definition: min.h:289
o2scl::min_base::print_iter
virtual int print_iter(double x, double y, int iter, double value=0.0, double limit=0.0, std::string comment="")
Print out iteration information.
Definition: min.h:90
o2scl::min_de_base
One-dimensional minimization using derivatives [abstract base].
Definition: min.h:302
o2scl::min_base::last_ntrial
int last_ntrial
The number of iterations used in the most recent minimization.
Definition: min.h:71
o2scl::success
@ success
Success.
Definition: err_hnd.h:47
o2scl_inte_qng_coeffs::x1
static const double x1[5]
Definition: inte_qng_gsl.h:48
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::min_de_base::min_de
virtual int min_de(double &x, double &fmin, func_t &func, dfunc_t &df)=0
Calculate the minimum min of func with derivative dfunc w.r.t 'x'.
o2scl::min_bkt_base
One-dimensional bracketing minimization [abstract base].
Definition: min.h:229
o2scl::min_bkt_base::min
virtual int min(double &x, double &fmin, func_t &func)
Calculate the minimum min of func w.r.t 'x'.
Definition: min.h:250
o2scl::min_de_base::min_bkt
virtual int min_bkt(double &x2, double x1, double x3, double &fmin, func_t &func)=0
Calculate the minimum min of func with x2 bracketed between x1 and x3.
o2scl::min_de_base::min
virtual int min(double &x, double &fmin, func_t &func)=0
Calculate the minimum min of func w.r.t 'x'.
o2scl::min_base::tol_rel
double tol_rel
The tolerance for the minimum function value.
Definition: min.h:65
o2scl::min_base
One-dimensional minimization [abstract base].
Definition: min.h:42
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::cont_constraint
double cont_constraint(double x, double center, double width, double height, double tightness=40.0, double exp_arg_limit=50.0)
Constrain x to be within width of the value given by center.
Definition: min.h:401
o2scl::min_base::min
virtual int min(double &x, double &fmin, func_t &func)=0
Calculate the minimum min of func w.r.t 'x'.
o2scl::constraint
double constraint(double x, double center, double width, double height)
Constrain x to be within width of the value given by center.
Definition: min.h:360
o2scl::min_base::min_bkt
virtual int min_bkt(double &x2, double x1, double x3, double &fmin, func_t &func)=0
Calculate the minimum min of func with x2 bracketed between x1 and x3.
o2scl::min_base::bracket
virtual int bracket(double &ax, double &bx, double &cx, double &fa, double &fb, double &fc, func_t &func)
Given interval (ax,bx), attempt to bracket a minimum for function func.
Definition: min.h:157
o2scl::lower_bound
double lower_bound(double x, double center, double width, double height)
Constrain x to be greater than the value given by center.
Definition: min.h:428
o2scl::min_base::err_nonconv
bool err_nonconv
If true, call the error handler if the routine does not "converge".
Definition: min.h:79
o2scl::min_base::verbose
int verbose
Output control.
Definition: min.h:59
o2scl_inte_qng_coeffs::x3
static const double x3[11]
Definition: inte_qng_gsl.h:94
o2scl::cont_lower_bound
double cont_lower_bound(double x, double center, double width, double height, double tightness=40.0, double exp_arg_limit=50.0)
Constrain x to be greater than the value given by center.
Definition: min.h:462

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