45 #ifndef O2SCL_MIN_QUAD_GOLDEN_H 46 #define O2SCL_MIN_QUAD_GOLDEN_H 53 #include <gsl/gsl_min.h> 54 #include <o2scl/min.h> 56 #ifndef DOXYGEN_NO_O2NS 108 #ifndef DOXYGEN_INTERNAL 121 double prev_stored_step;
131 double xlower,
double *flower,
double xupper,
134 *flower=func(xlower);
135 *fupper=func(xupper);
136 *fminimum=func(xminimum);
161 golden_mean=0.3819660112501052;
162 golden_ratio=1.6180339887498950;
168 int set(func_t &func,
double xmin,
double lower,
double upper) {
182 double lower,
double fl,
double upper,
186 std::string tmp=((std::string)
"Invalid interval (lower > upper) ")+
187 "in min_quad_golden::set_with_values().";
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().";
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().";
217 prev_stored_step=0.0;
233 double quad_step_size=prev_stored_step;
236 double x_eval, f_eval;
238 double x_midpoint=0.5*(x_l+x_u);
239 double tol=rel_err_val*fabs(x_m)+abs_err_val;
241 double dbl_eps=std::numeric_limits<double>::epsilon();
243 if (fabs(stored_step)-tol>-2.0*dbl_eps) {
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;
253 if (fabs(c2)>dbl_eps) {
258 quad_step_size=c1/c2;
263 quad_step_size=stored_step;
266 prev_stored_step=stored_step;
267 stored_step=step_size;
270 x_trial=x_m+quad_step_size;
272 if (fabs(quad_step_size)<fabs(0.5*prev_stored_step) &&
273 x_trial>x_l && x_trial<x_u) {
276 step_size=quad_step_size;
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);
283 }
else if ((x_small!=x_prev_small && x_small<x_m &&
285 (x_small!=x_prev_small && x_small>x_m &&
289 double outside_interval, inside_interval;
292 outside_interval=x_l-x_m;
293 inside_interval=x_u-x_m;
295 outside_interval=x_u-x_m;
296 inside_interval=x_l-x_m;
299 if (fabs(inside_interval) <= tol) {
301 double tmp=outside_interval;
302 outside_interval=inside_interval;
307 double step=inside_interval;
310 if (fabs(outside_interval)<fabs(inside_interval)) {
311 scale_factor=0.5*sqrt (-outside_interval/inside_interval);
313 scale_factor=(5.0/11.0)*
314 (0.1-inside_interval/outside_interval);
318 step_size=scale_factor*step;
326 if (x_m<x_midpoint) {
333 step_size=golden_mean*step;
338 if (fabs(step_size)>tol) {
339 x_eval=x_m+step_size;
341 x_eval=x_m+(step_size >= 0 ? +1.0 : -1.0)*fabs(tol);
346 f_eval=(*uf)(x_eval);
359 x_prev_small=x_small;
360 f_prev_small=f_small;
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;
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) {
411 int rx=
set(func,
x2,
x1,
x3);
414 "min_quad_golden::min_bkt().",rx);
425 "quad_golden::min_bkt().",status,this->
err_nonconv);
429 status=gsl_min_test_interval(x_lower,x_upper,this->
tol_abs,
435 std::string s=
"Function gsl_min_test_interval() failed in ";
436 s+=
"min_quad_golden::min_bkt().";
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)-
447 this->
print_iter(x2,f_minimum,iter,fabs(x_lower-x_upper)-
452 this->
print_iter(x2,f_minimum,iter,fabs(x_lower-x_upper),
453 this->
tol_abs,
"min_quad_golden");
459 if (iter>=this->ntrial) {
461 "iterations in min_quad_golden::min_bkt().",
472 virtual const char *
type() {
return "min_quad_golden"; }
476 #ifndef DOXYGEN_NO_O2NS One-dimensional bracketing minimization [abstract base].
double tol_rel
The tolerance for the minimum function value.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Minimization of a function using the safeguarded step-length algorithm of Gill and Murray [GSL]...
int verbose
Output control.
virtual const char * type()
Return string denoting type ("min_quad_golden")
invalid argument supplied by user
int ntrial
Maximum number of iterations.
double tol_abs
The tolerance for the location of the minimum.
exceeded max number of iterations
bool err_nonconv
If true, call the error handler if the routine does not "converge".
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
double f_upper
Value at upper bound.
int last_ntrial
The number of iterations used in the most recent minimization.
iteration has not converged
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.
double f_lower
Value at lower bound.
double x_upper
Upper bound.
static const double x3[11]
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
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.
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.
static const double x2[5]
static const double x1[5]
int iterate()
Perform an iteration.