diff options
Diffstat (limited to 'gsl-1.9/multimin/linear_minimize.c')
| -rw-r--r-- | gsl-1.9/multimin/linear_minimize.c | 247 |
1 files changed, 247 insertions, 0 deletions
diff --git a/gsl-1.9/multimin/linear_minimize.c b/gsl-1.9/multimin/linear_minimize.c new file mode 100644 index 0000000..60677d6 --- /dev/null +++ b/gsl-1.9/multimin/linear_minimize.c @@ -0,0 +1,247 @@ +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_poly.h> + +/* Find a minimum in x=[0,1] of the interpolating quadratic through + * (0,f0) (1,f1) with derivative fp0 at x=0. The interpolating + * polynomial is q(x) = f0 + fp0 * z + (f1-f0-fp0) * z^2 + */ + +static double +interp_quad (double f0, double fp0, double f1, double zl, double zh) +{ + double fl = f0 + zl*(fp0 + zl*(f1 - f0 -fp0)); + double fh = f0 + zh*(fp0 + zh*(f1 - f0 -fp0)); + double c = 2 * (f1 - f0 - fp0); /* curvature */ + + double zmin = zl, fmin = fl; + + if (fh < fmin) { zmin = zh; fmin = fh; } + + if (c > 0) /* positive curvature required for a minimum */ + { + double z = -fp0 / c; /* location of minimum */ + if (z > zl && z < zh) { + double f = f0 + z*(fp0 + z*(f1 - f0 -fp0)); + if (f < fmin) { zmin = z; fmin = f; }; + } + } + + return zmin; +} + +/* Find a minimum in x=[0,1] of the interpolating cubic through + * (0,f0) (1,f1) with derivatives fp0 at x=0 and fp1 at x=1. + * + * The interpolating polynomial is: + * + * c(x) = f0 + fp0 * z + eta * z^2 + xi * z^3 + * + * where eta=3*(f1-f0)-2*fp0-fp1, xi=fp0+fp1-2*(f1-f0). + */ + +static double +cubic (double c0, double c1, double c2, double c3, double z) +{ + return c0 + z * (c1 + z * (c2 + z * c3)); +} + +static void +check_extremum (double c0, double c1, double c2, double c3, double z, + double *zmin, double *fmin) +{ + /* could make an early return by testing curvature >0 for minimum */ + + double y = cubic (c0, c1, c2, c3, z); + + if (y < *fmin) + { + *zmin = z; /* accepted new point*/ + *fmin = y; + } +} + +static double +interp_cubic (double f0, double fp0, double f1, double fp1, double zl, double zh) +{ + double eta = 3 * (f1 - f0) - 2 * fp0 - fp1; + double xi = fp0 + fp1 - 2 * (f1 - f0); + double c0 = f0, c1 = fp0, c2 = eta, c3 = xi; + double zmin, fmin; + double z0, z1; + + zmin = zl; fmin = cubic(c0, c1, c2, c3, zl); + check_extremum (c0, c1, c2, c3, zh, &zmin, &fmin); + + { + int n = gsl_poly_solve_quadratic (3 * c3, 2 * c2, c1, &z0, &z1); + + if (n == 2) /* found 2 roots */ + { + if (z0 > zl && z0 < zh) + check_extremum (c0, c1, c2, c3, z0, &zmin, &fmin); + if (z1 > zl && z1 < zh) + check_extremum (c0, c1, c2, c3, z1, &zmin, &fmin); + } + else if (n == 1) /* found 1 root */ + { + if (z0 > zl && z0 < zh) + check_extremum (c0, c1, c2, c3, z0, &zmin, &fmin); + } + } + + return zmin; +} + + +static double +interpolate (double a, double fa, double fpa, + double b, double fb, double fpb, double xmin, double xmax, + int order) +{ + /* Map [a,b] to [0,1] */ + double z, alpha, zmin, zmax; + + zmin = (xmin - a) / (b - a); + zmax = (xmax - a) / (b - a); + + if (zmin > zmax) + { + double tmp = zmin; + zmin = zmax; + zmax = tmp; + }; + + if (order > 2 && GSL_IS_REAL(fpb)) { + z = interp_cubic (fa, fpa * (b - a), fb, fpb * (b - a), zmin, zmax); + } else { + z = interp_quad (fa, fpa * (b - a), fb, zmin, zmax); + } + + alpha = a + z * (b - a); + + return alpha; +} + +/* recommended values from Fletcher are + rho = 0.01, sigma = 0.1, tau1 = 9, tau2 = 0.05, tau3 = 0.5 */ + +static int +minimize (gsl_function_fdf * fn, double rho, double sigma, + double tau1, double tau2, double tau3, + int order, double alpha1, double *alpha_new) +{ + double f0, fp0, falpha, falpha_prev, fpalpha, fpalpha_prev, delta, + alpha_next; + double alpha = alpha1, alpha_prev = 0.0; + double a, b, fa, fb, fpa, fpb; + const size_t bracket_iters = 100, section_iters = 100; + size_t i = 0; + + GSL_FN_FDF_EVAL_F_DF (fn, 0.0, &f0, &fp0); + falpha_prev = f0; + fpalpha_prev = fp0; + + /* Avoid uninitialized variables morning */ + a = 0.0; b = alpha; + fa = f0; fb = 0.0; + fpa = fp0; fpb = 0.0; + + /* Begin bracketing */ + + while (i++ < bracket_iters) + { + falpha = GSL_FN_FDF_EVAL_F (fn, alpha); + + /* Fletcher's rho test */ + + if (falpha > f0 + alpha * rho * fp0 || falpha >= falpha_prev) + { + a = alpha_prev; fa = falpha_prev; fpa = fpalpha_prev; + b = alpha; fb = falpha; fpb = GSL_NAN; + break; /* goto sectioning */ + } + + fpalpha = GSL_FN_FDF_EVAL_DF (fn, alpha); + + /* Fletcher's sigma test */ + + if (fabs (fpalpha) <= -sigma * fp0) + { + *alpha_new = alpha; + return GSL_SUCCESS; + } + + if (fpalpha >= 0) + { + a = alpha; fa = falpha; fpa = fpalpha; + b = alpha_prev; fb = falpha_prev; fpb = fpalpha_prev; + break; /* goto sectioning */ + } + + delta = alpha - alpha_prev; + + { + double lower = alpha + delta; + double upper = alpha + tau1 * delta; + + alpha_next = interpolate (alpha_prev, falpha_prev, fpalpha_prev, + alpha, falpha, fpalpha, lower, upper, order); + + } + + alpha_prev = alpha; + falpha_prev = falpha; + fpalpha_prev = fpalpha; + alpha = alpha_next; + } + + /* Sectioning of bracket [a,b] */ + + while (i++ < section_iters) + { + delta = b - a; + + { + double lower = a + tau2 * delta; + double upper = b - tau3 * delta; + + alpha = interpolate (a, fa, fpa, b, fb, fpb, lower, upper, order); + } + + falpha = GSL_FN_FDF_EVAL_F (fn, alpha); + + if ((a-alpha)*fpa <= GSL_DBL_EPSILON) { + /* roundoff prevents progress */ + return GSL_ENOPROG; + }; + + if (falpha > f0 + rho * alpha * fp0 || falpha >= fa) + { + /* a_next = a; */ + b = alpha; fb = falpha; fpb = GSL_NAN; + } + else + { + fpalpha = GSL_FN_FDF_EVAL_DF (fn, alpha); + + if (fabs(fpalpha) <= -sigma * fp0) + { + *alpha_new = alpha; + return GSL_SUCCESS; /* terminate */ + } + + if ( ((b-a) >= 0 && fpalpha >= 0) || ((b-a) <=0 && fpalpha <= 0)) + { + b = a; fb = fa; fpb = fpa; + a = alpha; fa = falpha; fpa = fpalpha; + } + else + { + a = alpha; fa = falpha; fpa = fpalpha; + } + } + } + + return GSL_SUCCESS; +} |