min_quad_golden.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 /*--------------------------------------------------------------------------*/
24 /* */
25 /* quad_golden.c */
26 /* */
27 /* Copyright (C) 2007 James Howse */
28 /* Copyright (C) 2009 Brian Gough */
29 /* */
30 /* This program is free software; you can redistribute it and/or modify */
31 /* it under the terms of the GNU General Public License as published by */
32 /* the Free Software Foundation; either version 3 of the License, or (at */
33 /* your option) any later version. */
34 /* */
35 /* This program is distributed in the hope that it will be useful, but */
36 /* WITHOUT ANY WARRANTY; without even the implied warranty of */
37 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */
38 /* General Public License for more details. */
39 /* */
40 /* You should have received a copy of the GNU General Public License */
41 /* along with this program; if not, write to the Free Software */
42 /* Foundation, Inc., 51 Franklin Street, Fifth Floor, */
43 /* Boston, MA 02110-1301, USA. */
44 /*--------------------------------------------------------------------------*/
45 #ifndef O2SCL_MIN_QUAD_GOLDEN_H
46 #define O2SCL_MIN_QUAD_GOLDEN_H
47 
48 /** \file min_quad_golden.h
49  \brief File defining \ref o2scl::min_quad_golden
50 */
51 
52 #include <limits>
53 #include <gsl/gsl_min.h>
54 #include <o2scl/min.h>
55 
56 #ifndef DOXYGEN_NO_O2NS
57 namespace o2scl {
58 #endif
59 
60  /** \brief Minimization of a function using the safeguarded step-length
61  algorithm of Gill and Murray [GSL]
62 
63  This class is unfinished.
64 
65  Documentation from GSL:
66  \verbatim
67 
68  This algorithm performs univariate minimization (i.e., line
69  search). It requires only objective function values g(x) to
70  compute the minimum. The algorithm maintains an interval of
71  uncertainty [a,b] and a point x in the interval [a,b] such that
72  a < x < b, and g(a) > g(x) and g(x) < g(b). The algorithm also
73  maintains the three points with the smallest objective values x,
74  v and w such that g(x) < g(v) < g(w). The algorithm terminates
75  when max( x-a, b-x ) < 2(r |x| + t) where r and t are small
76  positive reals. At a given iteration, the algorithm first fits a
77  quadratic through the three points (x, g(x)), (v, g(v)) and (w,
78  g(w)) and computes the location of the minimum u of the
79  resulting quadratic. If u is in the interval [a,b] then g(u) is
80  computed. If u is not in the interval [a,b], and either v < x
81  and w < x, or v > x and w > x (i.e., the quadratic is
82  extrapolating), then a point u' is computed using a safeguarding
83  procedure and g(u') is computed. If u is not in the interval
84  [a,b], and the quadratic is not extrapolating, then a point u''
85  is computed using approximate golden section and g(u'') is
86  computed. After evaluating g() at the appropriate new point, a,
87  b, x, v, and w are updated and the next iteration is performed.
88  The algorithm is based on work presented in the following
89  references.
90 
91  Algorithms for Minimization without derivatives
92  Richard Brent
93  Prentice-Hall Inc., Englewood Cliffs, NJ, 1973
94 
95  Safeguarded Steplength Algorithms for Optimization using Descent Methods
96  Philip E. Gill and Walter Murray
97  Division of Numerical Analysis and Computing
98  National Physical Laboratory, Teddington, United Kingdom
99  NPL Report NAC 37, August 1974
100  \endverbatim
101 
102  \future Take common elements of this and min_brent and
103  move to a generic GSL minimizer type?
104  */
105  template<class func_t=funct> class min_quad_golden :
106  public min_bkt_base<func_t> {
107 
108 #ifndef DOXYGEN_INTERNAL
109 
110  protected:
111 
112  /// The function
113  func_t *uf;
114 
115  double x_prev_small;
116  double f_prev_small;
117  double x_small;
118  double f_small;
119  double step_size;
120  double stored_step;
121  double prev_stored_step;
122  size_t num_iter;
123 
124  double rel_err_val;
125  double abs_err_val;
126  double golden_mean;
127  double golden_ratio;
128 
129  /// Compute the function values at the various points
130  int compute_f_values(func_t &func, double xminimum, double *fminimum,
131  double xlower, double *flower, double xupper,
132  double *fupper) {
133 
134  *flower=func(xlower);
135  *fupper=func(xupper);
136  *fminimum=func(xminimum);
137 
138  return success;
139  }
140 
141 #endif
142 
143  public:
144 
145  /// Location of minimum
146  double x_minimum;
147  /// Lower bound
148  double x_lower;
149  /// Upper bound
150  double x_upper;
151  /// Minimum value
152  double f_minimum;
153  /// Value at lower bound
154  double f_lower;
155  /// Value at upper bound
156  double f_upper;
157 
158  min_quad_golden() {
159  rel_err_val=1.0e-06;
160  abs_err_val=1.0e-10;
161  golden_mean=0.3819660112501052;
162  golden_ratio=1.6180339887498950;
163  }
164 
165  virtual ~min_quad_golden() {}
166 
167  /// Set the function and the initial bracketing interval
168  int set(func_t &func, double xmin, double lower, double upper) {
169 
170  double fmin, fl, fu;
171  int status=compute_f_values(func,lower,&fl,xmin,&fmin,upper,&fu);
172 
173  status=set_with_values(func,xmin,fmin,lower,fl,upper,fu);
174 
175  return status;
176  }
177 
178  /** \brief Set the function, the initial bracketing interval,
179  and the corresponding function values.
180  */
181  int set_with_values(func_t &func, double xmin, double fmin,
182  double lower, double fl, double upper,
183  double fu) {
184 
185  if (lower > upper) {
186  std::string tmp=((std::string)"Invalid interval (lower > upper) ")+
187  "in min_quad_golden::set_with_values().";
188  O2SCL_ERR(tmp.c_str(),exc_einval);
189  }
190  if (xmin>=upper || xmin<=lower) {
191  std::string tmp=((std::string)"'xmin' was not inside interval ")+
192  "in min_quad_golden::set_with_values().";
193  O2SCL_ERR(tmp.c_str(),exc_einval);
194  }
195  if (fmin>= fl || fmin>=fu) {
196  std::string tmp=((std::string)"Endpoints don't enclose minimum ")+
197  "in min_quad_golden::set_with_values().";
198  O2SCL_ERR(tmp.c_str(),exc_einval);
199  }
200 
201  x_lower=lower;
202  x_minimum=xmin;
203  x_upper=upper;
204  f_lower=fl;
205  f_minimum=fmin;
206  f_upper=fu;
207 
208  uf=&func;
209 
210  x_prev_small=xmin;
211  x_small=xmin;
212  f_prev_small=fmin;
213  f_small=fmin;
214 
215  step_size=0.0;
216  stored_step=0.0;
217  prev_stored_step=0.0;
218  num_iter=0;
219 
220  return success;
221  }
222 
223  /** \brief Perform an iteration
224  */
225  int iterate() {
226 
227  const double x_m=x_minimum;
228  const double f_m=f_minimum;
229 
230  const double x_l=x_lower;
231  const double x_u=x_upper;
232 
233  double quad_step_size=prev_stored_step;
234 
235  double x_trial;
236  double x_eval, f_eval;
237 
238  double x_midpoint=0.5*(x_l+x_u);
239  double tol=rel_err_val*fabs(x_m)+abs_err_val;
240 
241  double dbl_eps=std::numeric_limits<double>::epsilon();
242 
243  if (fabs(stored_step)-tol>-2.0*dbl_eps) {
244 
245  // [GSL] Fit quadratic
246  double c3=(x_m-x_small)*(f_m-f_prev_small);
247  double c2=(x_m-x_prev_small)*(f_m-f_small);
248  double c1=(x_m-x_prev_small)*c2-(x_m-x_small)*c3;
249 
250  c2=2.0*(c2-c3);
251 
252  // [GSL] if (c2!=0)
253  if (fabs(c2)>dbl_eps) {
254  if (c2>0.0) {
255  c1=-c1;
256  }
257  c2=fabs(c2);
258  quad_step_size=c1/c2;
259  } else {
260  // [GSL] Handle case where c2 is nearly 0. Insure that the
261  // line search will NOT take a quadratic interpolation step
262  // in this iteration
263  quad_step_size=stored_step;
264  }
265 
266  prev_stored_step=stored_step;
267  stored_step=step_size;
268  }
269 
270  x_trial=x_m+quad_step_size;
271 
272  if (fabs(quad_step_size)<fabs(0.5*prev_stored_step) &&
273  x_trial>x_l && x_trial<x_u) {
274 
275  // [GSL] Take quadratic interpolation step
276  step_size=quad_step_size;
277 
278  // [GSL] Do not evaluate function too close to x_l or x_u
279  if ((x_trial-x_l)<2.0*tol || (x_u-x_trial)<2.0*tol) {
280  step_size=(x_midpoint >= x_m ? +1.0 : -1.0)*fabs(tol);
281  }
282 
283  } else if ((x_small!=x_prev_small && x_small<x_m &&
284  x_prev_small<x_m) ||
285  (x_small!=x_prev_small && x_small>x_m &&
286  x_prev_small>x_m)) {
287 
288  // [GSL] Take safeguarded function comparison step
289  double outside_interval, inside_interval;
290 
291  if (x_small<x_m) {
292  outside_interval=x_l-x_m;
293  inside_interval=x_u-x_m;
294  } else {
295  outside_interval=x_u-x_m;
296  inside_interval=x_l-x_m;
297  }
298 
299  if (fabs(inside_interval) <= tol) {
300  // [GSL] Swap inside and outside intervals
301  double tmp=outside_interval;
302  outside_interval=inside_interval;
303  inside_interval=tmp;
304  }
305 
306  {
307  double step=inside_interval;
308  double scale_factor;
309 
310  if (fabs(outside_interval)<fabs(inside_interval)) {
311  scale_factor=0.5*sqrt (-outside_interval/inside_interval);
312  } else {
313  scale_factor=(5.0/11.0)*
314  (0.1-inside_interval/outside_interval);
315  }
316 
317  stored_step=step;
318  step_size=scale_factor*step;
319  }
320 
321  } else {
322 
323  // [GSL] Take golden section step
324  double step;
325 
326  if (x_m<x_midpoint) {
327  step=x_u-x_m;
328  } else {
329  step=x_l-x_m;
330  }
331 
332  stored_step=step;
333  step_size=golden_mean*step;
334 
335  }
336 
337  // [GSL] Do not evaluate function too close to x_minimum
338  if (fabs(step_size)>tol) {
339  x_eval=x_m+step_size;
340  } else {
341  x_eval=x_m+(step_size >= 0 ? +1.0 : -1.0)*fabs(tol);
342  }
343 
344 
345  // [GSL] Evaluate function at the new point x_eval
346  f_eval=(*uf)(x_eval);
347 
348  // [GSL] Update {x,f}_lower, {x,f}_upper, {x,f}_prev_small,
349  // {x,f}_small, and {x,f}_minimum
350  if (f_eval <= f_m) {
351  if (x_eval<x_m) {
352  x_upper=x_m;
353  f_upper=f_m;
354  } else {
355  x_lower=x_m;
356  f_upper=f_m;
357  }
358 
359  x_prev_small=x_small;
360  f_prev_small=f_small;
361 
362  x_small=x_m;
363  f_small=f_m;
364 
365  x_minimum=x_eval;
366  f_minimum=f_eval;
367  } else {
368  if (x_eval<x_m) {
369  x_lower=x_eval;
370  f_lower=f_eval;
371  } else {
372  x_upper=x_eval;
373  f_upper=f_eval;
374  }
375 
376  if (f_eval <= f_small || fabs(x_small-x_m)<2.0*dbl_eps) {
377  x_prev_small=x_small;
378  f_prev_small=f_small;
379 
380  x_small=x_eval;
381  f_small=f_eval;
382  } else if (f_eval <= f_prev_small ||
383  fabs(x_prev_small-x_m)<2.0*dbl_eps ||
384  fabs(x_prev_small-x_small)<2.0*dbl_eps) {
385  x_prev_small=x_eval;
386  f_prev_small=f_eval;
387  }
388  }
389 
390  // Update for next iteration
391  num_iter++;
392 
393  return success;
394  }
395 
396  /** \brief Calculate the minimum \c fmin of \c func
397  with \c x2 bracketed between \c x1 and \c x3.
398 
399  Note that this algorithm requires that the initial guess
400  already brackets the minimum, i.e. \f$ x_1<x_2<x_3 \f$,
401  \f$ f(x_1)>f(x_2) \f$ and \f$ f(x_3)>f(x_2) \f$. This is
402  different from \ref min_cern, where the initial value
403  of the first parameter to min_cern::min_bkt() is
404  ignored.
405  */
406  virtual int min_bkt(double &x2, double x1, double x3, double &fmin,
407  func_t &func) {
408 
409  int status, iter=0;
410 
411  int rx=set(func,x2,x1,x3);
412  if (rx!=0) {
413  O2SCL_ERR2("Function set() failed in ",
414  "min_quad_golden::min_bkt().",rx);
415  }
416 
417  do {
418  iter++;
419  status=iterate();
420  x2=x_minimum;
421  if (status) {
422  // This should never actually happen in the current
423  // version, but we retain it for now
424  O2SCL_CONV2_RET("Function iterate() failed in gsl_min_",
425  "quad_golden::min_bkt().",status,this->err_nonconv);
426  break;
427  }
428 
429  status=gsl_min_test_interval(x_lower,x_upper,this->tol_abs,
430  this->tol_rel);
431  if (status>0) {
432  // The function gsl_min_test_interval() fails if the
433  // tolerances are negative or if the lower bound is larger
434  // than the upper bound
435  std::string s="Function gsl_min_test_interval() failed in ";
436  s+="min_quad_golden::min_bkt().";
437  O2SCL_ERR(s.c_str(),status);
438  }
439 
440  if (this->verbose>0) {
441  if (x_lower*x_upper<0.0) {
442  if (x_lower<x_upper) {
443  this->print_iter(x2,f_minimum,iter,fabs(x_lower-x_upper)-
444  this->tol_rel*x_lower,this->tol_abs,
445  "min_quad_golden");
446  } else {
447  this->print_iter(x2,f_minimum,iter,fabs(x_lower-x_upper)-
448  this->tol_rel*x_upper,this->tol_abs,
449  "min_quad_golden");
450  }
451  } else {
452  this->print_iter(x2,f_minimum,iter,fabs(x_lower-x_upper),
453  this->tol_abs,"min_quad_golden");
454  }
455  }
456 
457  } while (status == gsl_continue && iter<this->ntrial);
458 
459  if (iter>=this->ntrial) {
460  O2SCL_CONV2_RET("Exceeded maximum number of ",
461  "iterations in min_quad_golden::min_bkt().",
463  }
464 
465  this->last_ntrial=iter;
466  fmin=f_minimum;
467 
468  return status;
469  }
470 
471  /// Return string denoting type ("min_quad_golden")
472  virtual const char *type() { return "min_quad_golden"; }
473 
474  };
475 
476 #ifndef DOXYGEN_NO_O2NS
477 }
478 #endif
479 
480 #endif
o2scl::min_quad_golden::compute_f_values
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.
Definition: min_quad_golden.h:130
O2SCL_CONV2_RET
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:303
o2scl::min_quad_golden::f_minimum
double f_minimum
Minimum value.
Definition: min_quad_golden.h:152
o2scl::min_quad_golden::set
int set(func_t &func, double xmin, double lower, double upper)
Set the function and the initial bracketing interval.
Definition: min_quad_golden.h:168
o2scl::min_quad_golden::min_bkt
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.
Definition: min_quad_golden.h:406
o2scl::min_base< funct, funct >::tol_abs
double tol_abs
The tolerance for the location of the minimum.
Definition: min.h:68
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::min_quad_golden::type
virtual const char * type()
Return string denoting type ("min_quad_golden")
Definition: min_quad_golden.h:472
o2scl::min_base< funct, funct >::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_quad_golden::x_upper
double x_upper
Upper bound.
Definition: min_quad_golden.h:150
o2scl::exc_emaxiter
@ exc_emaxiter
exceeded max number of iterations
Definition: err_hnd.h:73
o2scl::min_quad_golden::set_with_values
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.
Definition: min_quad_golden.h:181
o2scl_inte_qng_coeffs::x2
static const double x2[5]
Definition: inte_qng_gsl.h:66
o2scl::min_base< funct, funct >::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_quad_golden::iterate
int iterate()
Perform an iteration.
Definition: min_quad_golden.h:225
o2scl::min_base< funct, funct >::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::min_quad_golden::f_upper
double f_upper
Value at upper bound.
Definition: min_quad_golden.h:156
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::min_bkt_base
One-dimensional bracketing minimization [abstract base].
Definition: min.h:229
o2scl::gsl_continue
@ gsl_continue
iteration has not converged
Definition: err_hnd.h:51
o2scl::min_base< funct, funct >::tol_rel
double tol_rel
The tolerance for the minimum function value.
Definition: min.h:65
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::min_quad_golden::f_lower
double f_lower
Value at lower bound.
Definition: min_quad_golden.h:154
o2scl::min_quad_golden::uf
func_t * uf
The function.
Definition: min_quad_golden.h:113
o2scl::min_quad_golden::x_minimum
double x_minimum
Location of minimum.
Definition: min_quad_golden.h:146
o2scl::min_quad_golden::x_lower
double x_lower
Lower bound.
Definition: min_quad_golden.h:148
o2scl::min_quad_golden
Minimization of a function using the safeguarded step-length algorithm of Gill and Murray [GSL].
Definition: min_quad_golden.h:105
o2scl::min_base< funct, funct >::err_nonconv
bool err_nonconv
If true, call the error handler if the routine does not "converge".
Definition: min.h:79
o2scl::min_base< funct, funct >::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

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