Skip to content
Snippets Groups Projects
Commit e6af344c authored by Stefen Guzik's avatar Stefen Guzik
Browse files

Reworked Brent 1D minimizer
  -- Starts with known parameter brackets for Curves
  -- Allows for min at interval extent
  -- Tolerance at maximum for gsl_brent
parent 4157fae1
No related branches found
No related tags found
No related merge requests found
......@@ -2874,11 +2874,9 @@ static Vertex *VERTEX;
static double min1d(double (*funct) (double), double *xmin)
{
// we should think about the tolerance more carefully...
double ax = 1.e-15, bx = 1.e-12, cx = 1.e-11, fa, fb, fx, tol = 1.e-4;
mnbrak(&ax, &bx, &cx, &fa, &fx, &fb, funct);
//Msg::Info("--MIN1D : ax %12.5E bx %12.5E cx %12.5E",ax,bx,cx);
return (brent(ax, bx, cx, funct, tol, xmin));
// 0. for tolerance allows for maximum as code in gsl_brent
return (brent(CURVE->ubeg, 0.5*(CURVE->ubeg + CURVE->uend), CURVE->uend,
funct, 0., xmin));
}
static void projectPS(int N, double x[], double res[])
......@@ -2900,10 +2898,6 @@ static void projectPS(int N, double x[], double res[])
static double projectPC(double u)
{
if(u < CURVE->ubeg)
u = CURVE->ubeg;
if(u < CURVE->ubeg)
u = CURVE->ubeg;
Vertex c = InterpolateCurve(CURVE, u, 0);
return sqrt(SQU(c.Pos.X - VERTEX->Pos.X) +
SQU(c.Pos.Y - VERTEX->Pos.Y) +
......
......@@ -20,6 +20,29 @@ double fn1(double x, void *params)
return val;
}
int gsl_min_fminimizer_set_with_values_MOD(gsl_min_fminimizer *s,
gsl_function *f,
double x_minimum, double f_minimum,
double x_lower, double f_lower,
double x_upper, double f_upper)
{
s->function = f;
s->x_minimum = x_minimum;
s->x_lower = x_lower;
s->x_upper = x_upper;
if (x_lower > x_upper)
{
GSL_ERROR ("invalid interval (lower > upper)", GSL_EINVAL);
}
s->f_lower = f_lower;
s->f_upper = f_upper;
s->f_minimum = f_minimum;
return (s->type->set) (s->state, s->function,
x_minimum, f_minimum,
x_lower, f_lower,
x_upper, f_upper);
}
#define MAXITER 100
// Returns the minimum betwen ax and cx to a given tolerance tol using
......@@ -28,9 +51,14 @@ double fn1(double x, void *params)
double brent(double ax, double bx, double cx,
double (*f) (double), double tol, double *xmin)
{
// The solver can stall at the following internal tolerance in brent - so
// this is about the lowest possible tolerance.
if(tol == 0.) tol = 10.*GSL_SQRT_DBL_EPSILON;
const double tolsq = tol*tol;
int status;
int iter = 0;
double a, b, c; // a < b < c
double fa, fb, fc;
const gsl_min_fminimizer_type *T;
gsl_min_fminimizer *s;
gsl_function F;
......@@ -42,38 +70,56 @@ double brent(double ax, double bx, double cx,
c = cx;
}
else {
a = ax;
c = cx;
}
// if a-b < tol, return func(a)
if(fabs(c - a) < tol) {
*xmin = ax;
return (f(*xmin));
a = cx;
c = ax;
}
nrfunc = f;
F.function = &fn1;
F.params = 0;
fa = f(a);
fb = f(b);
fc = f(c);
if(fb > fa || fb > fc) {
if(a < c) {
b = a;
fb = fa;
}
else {
b = c;
fb = fc;
}
}
// if a-c < tol, return func(b)
if(gsl_min_test_interval(a, c, tolsq, tol) == GSL_SUCCESS) {
*xmin = b;
return (fb);
}
T = gsl_min_fminimizer_brent;
s = gsl_min_fminimizer_alloc(T);
gsl_min_fminimizer_set(s, &F, b, a, c);
// gsl_min_fminimizer_set(s, &F, b, a, c); // Standard
// This mod avoids some error checks so we can have a minimum at an extent.
gsl_min_fminimizer_set_with_values_MOD(s, &F, b, fb, a, fa, c, fc);
do {
iter++;
status = gsl_min_fminimizer_iterate(s);
if(status)
break; // solver problem
b = gsl_min_fminimizer_minimum(s);
// this is deprecated: we should use gsl_min_fminimizer_x_minimum(s) instead
a = gsl_min_fminimizer_x_lower(s);
c = gsl_min_fminimizer_x_upper(s);
// we should think about the tolerance more carefully...
status = gsl_min_test_interval(a, c, tol, tol);
// fb = gsl_min_fminimizer_f_minimum(s);
// fa = gsl_min_fminimizer_f_lower(s);
// fc = gsl_min_fminimizer_f_upper(s);
status = gsl_min_test_interval(a, c, tolsq, tol);
// ... in case curve becomes too flat - not required?
// std::fabs(fa - fc) < (std::fabs(fb) + tol)*tol
}
while(status == GSL_CONTINUE && iter < MAXITER);
b = gsl_min_fminimizer_x_minimum(s);
if(status != GSL_SUCCESS)
Msg::Error("MIN1D not converged: f(%g) = %g", b, fn1(b, NULL));
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment