min_brent_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 /* min/brent.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
26  *
27  * This program is free software; you can redistribute it and/or modify
28  * it under the terms of the GNU General Public License as published by
29  * the Free Software Foundation; either version 3 of the License, or (at
30  * your option) any later version.
31  *
32  * This program is distributed in the hope that it will be useful, but
33  * WITHOUT ANY WARRANTY; without even the implied warranty of
34  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35  * General Public License for more details.
36  *
37  * You should have received a copy of the GNU General Public License
38  * along with this program; if not, write to the Free Software
39  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
40  * 02110-1301, USA.
41  */
42 #ifndef O2SCL_MIN_BRENT_GSL_H
43 #define O2SCL_MIN_BRENT_GSL_H
44 
45 /** \file min_brent_gsl.h
46  \brief File defining \ref o2scl::min_brent_gsl
47 */
48 #include <limits>
49 
50 // For gsl_min_test_interval()
51 #include <gsl/gsl_min.h>
52 
53 #include <o2scl/min.h>
54 
55 
56 #ifndef DOXYGEN_NO_O2NS
57 namespace o2scl {
58 #endif
59 
60  /** \brief One-dimensional minimization using Brent's method (GSL)
61 
62  The minimization in the function min_bkt() is complete when the
63  bracketed interval is smaller than
64  \f$ \mathrm{tol} = \mathrm{tol\_abs} + \mathrm{tol\_rel} \cdot
65  \mathrm{min} \f$, where \f$ \mathrm{min} =
66  \mathrm{min}(|\mathrm{lower}|,|\mathrm{upper}|) \f$.
67 
68  Note that this algorithm requires that the initial guess
69  already brackets the minimum, i.e. \f$ x_1 < x_2 < x_3 \f$,
70  \f$ f(x_1) > f(x_2) \f$ and \f$ f(x_3) > f(x_2) \f$. This is
71  different from \ref min_cern, where the initial value
72  of the first parameter to min_cern::min_bkt() is
73  ignored.
74 
75  The set functions throw an error if the initial bracket is not
76  correctly specified. The function \ref iterate() never calls the
77  error handler. The function min_bkt() calls the error handler if
78  the tolerances are negative or if the number of iterations is
79  insufficient to give the specified tolerance.
80 
81  Setting \ref min_base::err_nonconv to false will force
82  min_bkt() not to call the error handler when the number of
83  iterations is insufficient.
84 
85  Note that if there are more than 1 local minima in the specified
86  interval, there is no guarantee that this method will find the
87  global minimum.
88 
89  See also \ref root_brent_gsl for a similar algorithm
90  applied as a solver rather than a minimizer.
91 
92  \note There was a bug in this minimizer which was fixed for
93  GSL-1.11 which has also been fixed here.
94  */
95  template<class func_t=funct11> class min_brent_gsl :
96  public min_bkt_base<func_t> {
97 
98 #ifndef DOXYGEN_INTERNAL
99 
100  protected:
101 
102  /// The function
103  func_t *uf;
104 
105  /// \name Temporary storage
106  //@{
107  double d, e, v, w, f_v, f_w;
108  //@}
109 
110  /// Compute the function values at the various points
111  int compute_f_values(func_t &func, double xminimum, double *fminimum,
112  double xlower, double *flower, double xupper,
113  double *fupper) {
114 
115  *flower=func(xlower);
116  *fupper=func(xupper);
117  *fminimum=func(xminimum);
118 
119  return success;
120  }
121 
122 #endif
123 
124  public:
125 
126  /// Location of minimum
127  double x_minimum;
128  /// Lower bound
129  double x_lower;
130  /// Upper bound
131  double x_upper;
132  /// Minimum value
133  double f_minimum;
134  /// Value at lower bound
135  double f_lower;
136  /// Value at upper bound
137  double f_upper;
138 
139  min_brent_gsl() {
140  }
141 
142  virtual ~min_brent_gsl() {}
143 
144  /// Set the function and the initial bracketing interval
145  int set(func_t &func, double xmin, double lower, double upper) {
146 
147  double fmin, fl, fu;
148  int status=compute_f_values(func,lower,&fl,xmin,&fmin,upper,&fu);
149 
150  status=set_with_values(func,xmin,fmin,lower,fl,upper,fu);
151 
152  return status;
153  }
154 
155  /** \brief Set the function, the initial bracketing interval,
156  and the corresponding function values.
157  */
158  int set_with_values(func_t &func, double xmin, double fmin,
159  double lower, double fl, double upper,
160  double fu) {
161 
162  if (lower > upper) {
163  std::string tmp=((std::string)"Invalid interval (lower > upper) ")+
164  "in min_brent_gsl::set_with_values().";
165  O2SCL_ERR(tmp.c_str(),exc_einval);
166  }
167  if (xmin>=upper || xmin<=lower) {
168  std::string tmp=((std::string)"'xmin' was not inside interval ")+
169  "in min_brent_gsl::set_with_values().";
170  O2SCL_ERR(tmp.c_str(),exc_einval);
171  }
172  if (fmin>=fl || fmin>=fu) {
173  std::string tmp=((std::string)"Endpoints don't enclose minimum ")+
174  "in min_brent_gsl::set_with_values().";
175  O2SCL_ERR(tmp.c_str(),exc_einval);
176  }
177 
178  x_lower=lower;
179  x_minimum=xmin;
180  x_upper=upper;
181  f_lower=fl;
182  f_minimum=fmin;
183  f_upper=fu;
184 
185  uf=&func;
186 
187  /* golden=(3-sqrt(5))/2 */
188  const double golden=0.3819660;
189 
190  v=x_lower+golden*(x_upper-x_lower);
191  w=v;
192  d=0;
193  e=0;
194  f_v=func(v);
195  f_w=f_v;
196 
197  return success;
198  }
199 
200  /** \brief Perform an iteration
201 
202  \future It looks like x_left and x_right can be removed. Also,
203  it would be great to replace the one-letter variable names
204  with something more meaningful.
205  */
206  int iterate() {
207 
208  const double x_left=x_lower;
209  const double x_right=x_upper;
210 
211  const double z=x_minimum;
212  double u, f_u;
213  const double f_z=f_minimum;
214 
215  /* golden=(3-sqrt(5))/2 */
216  const double golden=0.3819660;
217 
218  const double w_lower=(z-x_left);
219  const double w_upper=(x_right-z);
220 
221  double sqrt_dbl_eps=sqrt(std::numeric_limits<double>::epsilon());
222 
223  const double tolerance=sqrt_dbl_eps*fabs(z);
224 
225  double p=0, q=0, r=0;
226 
227  const double midpoint=0.5*(x_left+x_right);
228 
229  if (fabs (e) > tolerance) {
230  /* fit parabola */
231 
232  r=(z-w)*(f_z-f_v);
233  q=(z-v)*(f_z-f_w);
234  p=(z-v)*q-(z-w)*r;
235  q=2*(q-r);
236 
237  if (q > 0) {
238  p=-p;
239  } else {
240  q=-q;
241  }
242 
243  r=e;
244  e=d;
245  }
246 
247  if (fabs (p) < fabs (0.5*q*r) && p < q*w_lower &&
248  p < q*w_upper) {
249 
250  double t2=2*tolerance;
251 
252  d=p / q;
253  u=z+d;
254 
255  if ((u-x_left) < t2 || (x_right-u) < t2) {
256  d=(z < midpoint) ? tolerance : -tolerance ;
257  }
258 
259  } else {
260 
261  e=(z < midpoint) ? x_right-z : -(z-x_left) ;
262  d=golden*e;
263 
264  }
265 
266 
267  if (fabs (d) >= tolerance) {
268  u=z+d;
269  } else {
270  u=z+((d > 0) ? tolerance : -tolerance);
271  }
272 
273  f_u=(*uf)(u);
274 
275  if (f_u<=f_z) {
276 
277  if (u<z) {
278 
279  x_upper=z;
280  f_upper=f_z;
281 
282  } else {
283 
284  x_lower=z;
285  f_lower=f_z;
286 
287  }
288 
289  v=w;
290  f_v=f_w;
291  w=z;
292  f_w=f_z;
293  x_minimum=u;
294  f_minimum=f_u;
295 
296  return success;
297 
298  } else {
299 
300  if (u<z) {
301 
302  x_lower=u;
303  f_lower=f_u;
304 
305  } else {
306 
307  x_upper=u;
308  f_upper=f_u;
309 
310  }
311 
312  if (f_u<=f_w || w==z) {
313 
314  v=w;
315  f_v=f_w;
316  w=u;
317  f_w=f_u;
318  return success;
319 
320  } else if (f_u<=f_v || v==z || v==w) {
321 
322  v=u;
323  f_v=f_u;
324  return success;
325  }
326 
327  }
328 
329  return success;
330  }
331 
332  /** \brief Calculate the minimum \c fmin of \c func
333  with \c x2 bracketed between \c x1 and \c x3.
334 
335  Note that this algorithm requires that the initial guess
336  already brackets the minimum, i.e. \f$ x_1 < x_2 < x_3 \f$,
337  \f$ f(x_1) > f(x_2) \f$ and \f$ f(x_3) > f(x_2) \f$. This is
338  different from \ref min_cern, where the initial value
339  of the first parameter to min_cern::min_bkt() is
340  ignored.
341  */
342  virtual int min_bkt(double &x2, double x1, double x3, double &fmin,
343  func_t &func) {
344 
345  int status, iter=0;
346 
347  int rx=set(func,x2,x1,x3);
348  if (rx!=0) {
349  O2SCL_ERR2("Function set() failed in ",
350  "min_brent_gsl::min_bkt().",rx);
351  }
352 
353  do {
354  iter++;
355  status=iterate();
356  x2=x_minimum;
357  if (status) {
358  // This should never actually happen in the current
359  // version, but we retain it for now
360  O2SCL_CONV2_RET("Function iterate() failed in gsl_min_",
361  "brent::min_bkt().",status,this->err_nonconv);
362  break;
363  }
364 
365  status=gsl_min_test_interval(x_lower,x_upper,this->tol_abs,
366  this->tol_rel);
367  if (status>0) {
368  // The function gsl_min_test_interval() fails if the
369  // tolerances are negative or if the lower bound is larger
370  // than the upper bound
371  std::string s="Function gsl_min_test_interval() failed in ";
372  s+="min_brent_gsl::min_bkt().";
373  O2SCL_ERR(s.c_str(),status);
374  }
375 
376  if (this->verbose>0) {
377  if (x_lower*x_upper<0.0) {
378  if (x_lower<x_upper) {
379  this->print_iter(x2,f_minimum,iter,fabs(x_lower-x_upper)-
380  this->tol_rel*x_lower,this->tol_abs,
381  "min_brent_gsl");
382  } else {
383  this->print_iter(x2,f_minimum,iter,fabs(x_lower-x_upper)-
384  this->tol_rel*x_upper,this->tol_abs,
385  "min_brent_gsl");
386  }
387  } else {
388  this->print_iter(x2,f_minimum,iter,fabs(x_lower-x_upper),
389  this->tol_abs,"min_brent_gsl");
390  }
391  }
392 
393  } while (status == gsl_continue && iter<this->ntrial);
394 
395  if (iter>=this->ntrial) {
396  O2SCL_CONV2_RET("Exceeded maximum number of ",
397  "iterations in min_brent_gsl::min_bkt().",
399  }
400 
401  this->last_ntrial=iter;
402  fmin=f_minimum;
403 
404  return status;
405  }
406 
407  /// Return string denoting type ("min_brent_gsl")
408  virtual const char *type() { return "min_brent_gsl"; }
409 
410  };
411 
412 #ifndef DOXYGEN_NO_O2NS
413 }
414 #endif
415 
416 #endif
double f_minimum
Minimum value.
One-dimensional bracketing minimization [abstract base].
Definition: min.h:229
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
double x_upper
Upper bound.
int verbose
Output control.
Definition: min.h:59
double f_lower
Value at lower bound.
invalid argument supplied by user
Definition: err_hnd.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
int iterate()
Perform an iteration.
exceeded max number of iterations
Definition: err_hnd.h:73
int set_with_values(func_t &func, double xmin, double fmin, double lower, double fl, double upper, double fu)
Set the function, the initial bracketing interval, and the corresponding function values...
bool err_nonconv
If true, call the error handler if the routine does not "converge".
Definition: min.h:79
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:298
int last_ntrial
The number of iterations used in the most recent minimization.
Definition: min.h:71
iteration has not converged
Definition: err_hnd.h:51
One-dimensional minimization using Brent&#39;s method (GSL)
Definition: min_brent_gsl.h:95
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 const char * type()
Return string denoting type ("min_brent_gsl")
static const double x3[11]
Definition: inte_qng_gsl.h:94
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
double x_lower
Lower bound.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
double x_minimum
Location of minimum.
double f_upper
Value at upper bound.
func_t * uf
The function.
int compute_f_values(func_t &func, double xminimum, double *fminimum, double xlower, double *flower, double xupper, double *fupper)
Compute the function values at the various points.
static const double x2[5]
Definition: inte_qng_gsl.h:66
static const double x1[5]
Definition: inte_qng_gsl.h:48
Success.
Definition: err_hnd.h:47
virtual int min_bkt(double &x2, double x1, double x3, double &fmin, func_t &func)
Calculate the minimum fmin of func with x2 bracketed between x1 and x3.

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