min_quad_golden.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 /*--------------------------------------------------------------------------*/
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=funct11> 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
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
Minimization of a function using the safeguarded step-length algorithm of Gill and Murray [GSL]...
int verbose
Output control.
Definition: min.h:59
virtual const char * type()
Return string denoting type ("min_quad_golden")
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
exceeded max number of iterations
Definition: err_hnd.h:73
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
double f_upper
Value at upper bound.
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
double x_lower
Lower bound.
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
double f_lower
Value at lower bound.
double x_upper
Upper bound.
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
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...
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
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.
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.
double f_minimum
Minimum value.
double x_minimum
Location of minimum.
func_t * uf
The function.
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
int iterate()
Perform an iteration.

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