min_cern.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 /** \file min_cern.h
24  \brief File defining \ref o2scl::min_cern
25 */
26 #ifndef O2SCL_CERN_MINIMIZE_H
27 #define O2SCL_CERN_MINIMIZE_H
28 
29 #include <o2scl/funct.h>
30 #include <o2scl/min.h>
31 
32 #ifndef DOXYGEN_NO_O2NS
33 namespace o2scl {
34 #endif
35 
36  /** \brief One-dimensional minimization (CERNLIB)
37 
38  The golden section search is applied in the interval \f$ (a,b) \f$
39  using a fixed number \f$ n \f$ of function evaluations where
40  \f[
41  n=\left[ 2.08 \ln(|a-b|/\mathrm{tol\_abs})+\frac{1}{2}\right]+1
42  \f]
43 
44  The accuracy depends on the function. A choice of \f$
45  \mathrm{tol\_abs}>10^{-8} \f$ usually results in a relative error of
46  $x$ which is smaller than or of the order of \f$ \mathrm{tol\_abs}
47  \f$ .
48 
49  This routine strictly searches the interval \f$ (a,b) \f$ . If
50  the function is nowhere flat in this interval, then min_bkt()
51  will return either \f$ a \f$ or \f$ b \f$ and \ref min_type
52  is set to 1. Note that if there are more than 1 local minima
53  in the specified interval, there is no guarantee that this
54  method will find the global minimum.
55 
56  \note The number of function evaluations can be quite large if
57  mmin::tol_abs is sufficiently small. If mmin::tol_abs is
58  exactly zero, then the error handler will be called.
59 
60  \comment
61  If mmin::tol_abs is exactly zero, then \f$ 10^{-8} \f$
62  will be used instead.
63  [5/23/11 - Changed this consistent with the policy of
64  throwing on incorrect input.]
65  \endcomment
66 
67  Based on the CERNLIB routines RMINFC and DMINFC, which was
68  based on \ref Fletcher87, and \ref Krabs83 and is documented at
69  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d503/top.html
70  */
71 #ifdef DOXYGEN_NO_O2NS
72  template<class func_t=funct11> class min_cern :
73  public min_bkt_base
74 #else
75  template<class func_t=funct11> class min_cern :
76  public min_bkt_base<func_t>
77 #endif
78  {
79 
80 #ifndef DOXGYEN_INTERNAL
81 
82  protected:
83 
84  /// C analog of Fortran's "Nearest integer" function
85  inline int nint(double x) {
86  if (x<0.0) return ((int)(x-0.5));
87  else return ((int)(x+0.5));
88  }
89 
90 #endif
91 
92  public:
93 
94  min_cern() {
95  delta_set=false;
96  min_type=0;
97  }
98 
99  /** \brief Calculate the minimum \c min of \c func between
100  \c a and \c b
101 
102  The initial value of \c x is ignored.
103 
104  If there is no minimum in the given interval, then on exit
105  \c x will be equal to either \c a or \c b and \ref min_type
106  will be set to 1 instead of zero. The error handler is not
107  called, as this need not be interpreted as an error.
108  */
109  virtual int min_bkt(double &x, double a, double b, double &y,
110  func_t &func) {
111 
112  // The square root of 5
113  static const double w5=2.23606797749979;
114  static const double hv=(3.0-w5)/2.0, hw=(w5-1.0)/2.0;
115  double c, d, v=0.0, fv=0.0, w=0.0, fw=0.0, h;
116 
117  if (!delta_set) delta=10.0*this->tol_abs;
118 
119  int n=-1;
120  if (this->tol_abs==0.0) {
121  O2SCL_ERR("Tolerance zero in min_cern::min_bkt().",
122  exc_einval);
123  } else {
124  if (a!=b) n=nint(2.0*log(fabs((a-b)/this->tol_abs)));
125  }
126  this->last_ntrial=n;
127  c=a;
128  d=b;
129  if (a>b) {
130  c=b;
131  d=a;
132  }
133  bool llt=true, lge=true;
134 
135  // This is never an infinite loop so long as 'n' is finite. We
136  // should maybe check somehow that 'n' is not too large
137  while (true) {
138  h=d-c;
139  if (this->verbose>0) {
140  x=(c+d)/2.0;
141  y=func(x);
142  this->print_iter(x,y,this->last_ntrial-n,((double)n),this->tol_abs,
143  "min_cern");
144  }
145  if (n<0) {
146  x=(c+d)/2.0;
147  y=func(x);
148  if ((fabs(x-a)>delta && fabs(x-b)>delta)) min_type=0;
149  min_type=1;
150  return 0;
151  }
152  if (llt) {
153  v=c+hv*h;
154  fv=func(v);
155  }
156  if (lge) {
157  w=c+hw*h;
158  fw=func(w);
159  }
160  if (fv<fw) {
161  llt=true;
162  lge=false;
163  d=w;
164  w=v;
165  fw=fv;
166  } else {
167  llt=false;
168  lge=true;
169  c=v;
170  v=w;
171  fv=fw;
172  }
173  n--;
174  }
175 
176  O2SCL_ERR("Failed sanity check in min_cern::min_bkt().",
177  exc_esanity);
178  return exc_esanity;
179  }
180 
181  /** \brief Set the value of \f$ \delta \f$
182 
183  If this is not called before min_bkt() is used, then the
184  suggested value \f$ \delta=10 \mathrm{tol_abs} \f$ is used.
185  */
186  int set_delta(double d) {
187  if (d>0.0) {
188  delta=d;
189  delta_set=true;
190  }
191  return 0;
192  }
193 
194  /// Return string denoting type ("min_cern")
195  virtual const char *type() { return "min_cern"; }
196 
197  /// Type of minimum found
198  int min_type;
199 
200  protected:
201 
202 #ifndef DOXYGEN_INTERNAL
203 
204  /// The value of delta as specified by the user
205  double delta;
206 
207  /// True if the value of delta has been set
208  bool delta_set;
209 
210 #endif
211 
212  };
213 
214 #ifndef DOXYGEN_NO_O2NS
215 }
216 #endif
217 
218 #endif
One-dimensional bracketing minimization [abstract base].
Definition: min.h:229
int min_type
Type of minimum found.
Definition: min_cern.h:198
bool delta_set
True if the value of delta has been set.
Definition: min_cern.h:208
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
int verbose
Output control.
Definition: min.h:59
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
invalid argument supplied by user
Definition: err_hnd.h:59
double tol_abs
The tolerance for the location of the minimum.
Definition: min.h:68
int nint(double x)
C analog of Fortran&#39;s "Nearest integer" function.
Definition: min_cern.h:85
virtual const char * type()
Return string denoting type ("min_cern")
Definition: min_cern.h:195
int last_ntrial
The number of iterations used in the most recent minimization.
Definition: min.h:71
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
virtual int min_bkt(double &x, double a, double b, double &y, func_t &func)
Calculate the minimum min of func between a and b.
Definition: min_cern.h:109
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
One-dimensional minimization (CERNLIB)
Definition: min_cern.h:75
double delta
The value of delta as specified by the user.
Definition: min_cern.h:205
int set_delta(double d)
Set the value of .
Definition: min_cern.h:186

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