min.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 #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=funct11, 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  /** \brief Constrain \c x to be within \c width
343  of the value given by \c center
344 
345  Defining \f$ c= \f$ \c center, \f$ w= \f$ \c width, \f$ h= \f$
346  \c height, this returns the value \f$ h (1+|x-c-w|/w) \f$ if \f$
347  x>c+w \f$ and \f$ h (1+|x-c+w|/w) \f$ if \f$ x<c-w \f$ . The
348  value near \f$ x=c-w \f$ or \f$ x=c+w \f$ is \f$ h \f$ (the
349  value of the function exactly at these points is zero) and the
350  value at \f$ x=c-2w \f$ or \f$ x=c+2w \f$ is \f$ 2 h \f$ .
351 
352  It is important to note that, for large distances of \c x
353  from \c center, this only scales linearly. If you are trying to
354  constrain a function which decreases more than linearly by
355  making \c x far from \c center, then a minimizer may
356  ignore this constraint.
357  */
358  inline double constraint(double x, double center, double width,
359  double height) {
360  double ret=0.0;
361  if (x>center+width) {
362  ret=height*(1.0+fabs(x-center-width)/width);
363  } else if (x<center-width) {
364  ret=height*(1.0+fabs(x-center+width)/width);
365  }
366  return ret;
367  }
368 
369  /** \brief Constrain \c x to be within \c width of the value given
370  by \c center
371 
372  Defining \f$ c= \f$ \c center, \f$ w= \f$ \c width, \f$ h= \f$
373  \c height, \f$ t= \f$ \c tightness, and \f$ \ell= \f$ \c
374  exp_arg_limit, this returns the value
375  \f[
376  h \left(\frac{x-c}{w}\right)^2 \left[
377  1+ e^{t\left(x-c+w\right)\left(c+w-x\right)/w^2}
378  \right]^{-1}
379  \f]
380 
381  This function is continuous and differentiable. Note that if
382  \f$ x=c \f$ , then the function returns zero.
383 
384  The exponential is handled gracefully by assuming that anything
385  smaller than \f$ \exp(-\ell) \f$ is zero. This creates a small
386  discontinuity which can be removed with the sufficiently large
387  value of \f$ \ell \f$.
388 
389  It is important to note that, for large distances of \c x from
390  \c center, this scales quadratically. If you are trying to
391  constrain a function which decreases faster than quadratically
392  by making \c x far from \c center, then a minimizer may ignore
393  this constraint.
394 
395  In the limit \f$ t \rightarrow \infty \f$, this function
396  converges towards the squared value of \ref constraint(), except
397  exactly at the points \f$ x=c-w \f$ and \f$ x=c+w \f$.
398  */
399  inline double cont_constraint(double x, double center, double width,
400  double height, double tightness=40.0,
401  double exp_arg_limit=50.0) {
402  double ret, wid2=width*width;
403  double arg=tightness/wid2*(x-center+width)*(center+width-x);
404  if (arg<-exp_arg_limit) {
405  ret=(x-center)*(x-center)/wid2;
406  } else {
407  ret=(x-center)*(x-center)/wid2/(1.0+exp(arg));
408  }
409  return ret*height;
410  }
411 
412  /** \brief Constrain \c x to be greater than the value given by \c
413  center
414 
415  Defining \f$ c= \f$ \c center, \f$ w= \f$ \c width, \f$ h= \f$
416  \c height, this returns \f$ h(1+|x-c|/w) \f$ if \f$ x<c \f$ and
417  zero otherwise. The value at \f$ x=c \f$ is \f$ h \f$ , while
418  the value at \f$ x=c-w \f$ is \f$ 2 h \f$ .
419 
420  It is important to note that, for large distances of \c x from
421  \c center, this only scales linearly. If you are trying to
422  constrain a function which decreases more than linearly by
423  making \c x far from \c center, then a minimizer may ignore this
424  constraint.
425  */
426  inline double lower_bound(double x, double center, double width,
427  double height) {
428  double ret=0.0;
429  if (x<center) ret=height*(1.0+fabs(x-center)/width);
430  return ret;
431  }
432 
433  /** \brief Constrain \c x to be greater than the value given by \c
434  center
435 
436  Defining \f$ c= \f$ \c center, \f$ w= \f$ \c width, \f$ h= \f$
437  \c height, \f$ t= \f$ \c tightness, and \f$ \ell= \f$ \c
438  exp_arg_limit, this returns \f$ h(c-x+w)/(w+w\exp(t(x-c)/w)) \f$
439  and has the advantage of being a continuous and differentiable
440  function. The value of the function exactly at \f$ x=c \f$ is
441  \f$ h/2 \f$, but for \f$ x \f$ just below \f$ c \f$ the function
442  is \f$ h \f$ and just above \f$ c \f$ the function is quite
443  small.
444 
445  The exponential is handled gracefully by assuming that anything
446  smaller than \f$ \exp(-\ell) \f$ is zero. This creates a small
447  discontinuity which can be removed with the sufficiently large
448  value of \f$ \ell \f$.
449 
450  It is important to note that, for large distances of \c x
451  from \c center, this only scales linearly. If you are trying to
452  constrain a function which decreases more than linearly by
453  making \c x far from \c center, then a minimizer may
454  ingore this constraint.
455 
456  In the limit \f$ t \rightarrow \infty \f$, this function
457  converges towards \ref lower_bound(), except exactly at the
458  point \f$ x=c \f$.
459  */
460  inline double cont_lower_bound(double x, double center, double width,
461  double height, double tightness=40.0,
462  double exp_arg_limit=50.0) {
463  double ret, arg=tightness*(x-center)/width;
464  if (arg>exp_arg_limit) {
465  ret=0.0;
466  } else if (arg<-exp_arg_limit) {
467  ret=height*(center-x+width)/width;
468  } else {
469  ret=height*(center-x+width)/width/(1.0+exp(arg));
470  }
471  return ret;
472  }
473 
474 #ifndef DOXYGEN_NO_O2NS
475 }
476 #endif
477 
478 #endif
int bracket_iter
The number of iterations for automatically bracketing a minimum (default 20)
Definition: min.h:76
One-dimensional bracketing minimization [abstract base].
Definition: min.h:229
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:460
virtual const char * type()
Return string denoting type ("min_de")
Definition: min.h:338
double tol_rel
The tolerance for the minimum function value.
Definition: min.h:65
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
One-dimensional minimization using derivatives [abstract base].
Definition: min.h:302
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 &#39;x&#39;.
int verbose
Output control.
Definition: min.h:59
int ntrial
Maximum number of iterations.
Definition: min.h:62
double tol_abs
The tolerance for the location of the minimum.
Definition: min.h:68
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:399
exceeded max number of iterations
Definition: err_hnd.h:73
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:358
bool err_nonconv
If true, call the error handler if the routine does not "converge".
Definition: min.h:79
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
virtual int min(double &x, double &fmin, func_t &func)
Calculate the minimum min of func w.r.t &#39;x&#39;.
Definition: min.h:250
int last_ntrial
The number of iterations used in the most recent minimization.
Definition: min.h:71
generic failure
Definition: err_hnd.h:61
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.
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
int bracket_iter
The number of iterations for automatically bracketing a minimum (default 20)
Definition: min.h:242
static const double x3[11]
Definition: inte_qng_gsl.h:94
virtual int min(double &x, double &fmin, func_t &func)=0
Calculate the minimum min of func w.r.t &#39;x&#39;.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
One-dimensional minimization [abstract base].
Definition: min.h:42
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 &#39;x&#39;.
Definition: min.h:276
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:426
static const double x2[5]
Definition: inte_qng_gsl.h:66
static const double x1[5]
Definition: inte_qng_gsl.h:48
virtual const char * type()
Return string denoting type ("min")
Definition: min.h:222
virtual const char * type()
Return string denoting type ("min_bkt")
Definition: min.h:289
Success.
Definition: err_hnd.h:47

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