From 44d918d93b1bf4c0856531b551941f4a94b3397d Mon Sep 17 00:00:00 2001 From: stevenj Date: Wed, 29 Aug 2007 21:07:20 -0400 Subject: [PATCH] got rid of buggy f_recenter (Lagrange interpolation isn't guaranteed to be monotonic) darcs-hash:20070830010720-c8de0-19c1e18baffbe16a121a8c968f2c0d2a31f27058.gz --- api/nlopt.c | 116 +++++++--------------------------------------- cdirect/cdirect.c | 11 ++--- test/testopt.cpp | 18 +++++-- 3 files changed, 33 insertions(+), 112 deletions(-) diff --git a/api/nlopt.c b/api/nlopt.c index 0ed5bcb..669dd9f 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -75,75 +75,9 @@ void nlopt_srand_time() { typedef struct { nlopt_func f; void *f_data; - const double *lb, *ub, *x0; - double *xtmp; + const double *lb, *ub; } nlopt_data; -#define RECENTER 1 /* 0 to disable recentering */ - -/* for global-search algorithms that ignore the starting guess, - but always check the center of the search box, we perform a - coordinate transformation to put the initial guess x0 at the - center of the box, and store the transformed x in xtmp. */ -static void recenter_x(int n, const double *x, - const double *lb, const double *ub, - const double *x0, double *xtmp) -{ - int i; - for (i = 0; i < n; ++i) { -#if RECENTER - /* Lagrange interpolating polynomial */ - double xm = 0.5 * (lb[i] + ub[i]); - double dlu = 1. / (lb[i] - ub[i]); - double dlm = 1. / (lb[i] - xm); - double dum = 1. / (ub[i] - xm); - double dxu = x[i] - ub[i]; - double dxl = x[i] - lb[i]; - double dxm = x[i] - xm; - xtmp[i] = (lb[i] * (dxu * dlu) * (dxm * dlm) - - ub[i] * (dxl * dlu) * (dxm * dum) - + x0[i] * (dxl * dlm) * (dxu * dum)); -#else - xtmp[i] = x[i]; -#endif - } -} - -/* transform grad from df/dxtmp to df/dx */ -static void recenter_grad(int n, const double *x, - const double *lb, const double *ub, - const double *x0, - double *grad) -{ -#if RECENTER - if (grad) { - int i; - for (i = 0; i < n; ++i) { - double xm = 0.5 * (lb[i] + ub[i]); - double dlu = 1. / (lb[i] - ub[i]); - double dlm = 1. / (lb[i] - xm); - double dum = 1. / (ub[i] - xm); - double dxu = x[i] - ub[i]; - double dxl = x[i] - lb[i]; - double dxm = x[i] - xm; - grad[i] *= (lb[i] * dlu*dlm * (dxm + dxu) - - ub[i] * dum*dlu * (dxm + dxl) - + x0[i] * dlm*dum * (dxu + dxl)); - } - } -#endif -} - -static double f_recenter(int n, const double *x, double *grad, void *data_) -{ - nlopt_data *data = (nlopt_data *) data_; - double f; - recenter_x(n, x, data->lb, data->ub, data->x0, data->xtmp); - f = data->f(n, data->xtmp, grad, data->f_data); - recenter_grad(n, x, data->lb, data->ub, data->x0, grad); - return f; -} - #include "subplex.h" static double f_subplex(int n, const double *x, void *data_) @@ -168,8 +102,7 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) { nlopt_data *data = (nlopt_data *) data_; double f; - recenter_x(n, x, data->lb, data->ub, data->x0, data->xtmp); - f = data->f(n, data->xtmp, NULL, data->f_data); + f = data->f(n, x, NULL, data->f_data); *undefined = isnan(f) || my_isinf(f); return f; } @@ -205,7 +138,6 @@ static nlopt_result nlopt_minimize_( d.f_data = f_data; d.lb = lb; d.ub = ub; - d.x0 = d.xtmp = NULL; /* make sure rand generator is inited */ if (!nlopt_srand_called) @@ -235,25 +167,17 @@ static nlopt_result nlopt_minimize_( return cdirect(n, f, f_data, lb, ub, x, fmin, &stop, 0.0, (algorithm == NLOPT_GLOBAL_DIRECT ? 0 : 1) + 10 * (algorithm == NLOPT_GLOBAL_DIRECT_L_RANDOMIZED ? 2 : 0)); - + case NLOPT_GLOBAL_ORIG_DIRECT: case NLOPT_GLOBAL_ORIG_DIRECT_L: - { - int iret; - d.xtmp = (double *) malloc(sizeof(double) * n*2); - if (!d.xtmp) return NLOPT_OUT_OF_MEMORY; - memcpy(d.xtmp + n, x, sizeof(double) * n); d.x0 = d.xtmp + n; - iret = direct_optimize(f_direct, &d, n, lb, ub, x, fmin, - maxeval, -1, 0.0, 0.0, - pow(xtol_rel, (double) n), -1.0, - stop.fmin_max, 0.0, - NULL, - algorithm == NLOPT_GLOBAL_ORIG_DIRECT - ? DIRECT_ORIGINAL - : DIRECT_GABLONSKY); - recenter_x(n, x, lb, ub, d.x0, x); - free(d.xtmp); - switch (iret) { + switch (direct_optimize(f_direct, &d, n, lb, ub, x, fmin, + maxeval, -1, 0.0, 0.0, + pow(xtol_rel, (double) n), -1.0, + stop.fmin_max, 0.0, + NULL, + algorithm == NLOPT_GLOBAL_ORIG_DIRECT + ? DIRECT_ORIGINAL + : DIRECT_GABLONSKY)) { case DIRECT_INVALID_BOUNDS: case DIRECT_MAXFEVAL_TOOBIG: case DIRECT_INVALID_ARGS: @@ -272,24 +196,16 @@ static nlopt_result nlopt_minimize_( return NLOPT_XTOL_REACHED; case DIRECT_OUT_OF_MEMORY: return NLOPT_OUT_OF_MEMORY; - } break; } case NLOPT_GLOBAL_STOGO: - case NLOPT_GLOBAL_STOGO_RANDOMIZED: { - int iret; - d.xtmp = (double *) malloc(sizeof(double) * n*2); - if (!d.xtmp) return NLOPT_OUT_OF_MEMORY; - memcpy(d.xtmp + n, x, sizeof(double) * n); d.x0 = d.xtmp + n; - iret = stogo_minimize(n, f_recenter, &d, x, fmin, lb, ub, &stop, - algorithm == NLOPT_GLOBAL_STOGO - ? 0 : 2*n); - recenter_x(n, x, lb, ub, d.x0, x); - free(d.xtmp); - if (!iret) return NLOPT_FAILURE; + case NLOPT_GLOBAL_STOGO_RANDOMIZED: + if (!stogo_minimize(n, f, f_data, x, fmin, lb, ub, &stop, + algorithm == NLOPT_GLOBAL_STOGO + ? 0 : 2*n)) + return NLOPT_FAILURE; break; - } case NLOPT_LOCAL_SUBPLEX: { int iret; diff --git a/cdirect/cdirect.c b/cdirect/cdirect.c index 9d95f8f..52b55fa 100644 --- a/cdirect/cdirect.c +++ b/cdirect/cdirect.c @@ -429,7 +429,7 @@ nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data, double magic_eps, int which_alg) { params p; - int i, x_center = 1; + int i; double *rnew; nlopt_result ret = NLOPT_OUT_OF_MEMORY; @@ -443,7 +443,7 @@ nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data, p.f = f; p.f_data = f_data; p.xmin = x; - p.fmin = f(n, x, NULL, f_data); stop->nevals++; + p.fmin = HUGE_VAL; p.work = 0; p.iwork = 0; p.hull = 0; @@ -461,15 +461,10 @@ nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data, if (!(rnew = (double *) malloc(sizeof(double) * p.L))) goto done; for (i = 0; i < n; ++i) { rnew[2+i] = 0.5 * (lb[i] + ub[i]); - x_center = x_center - && (fabs(rnew[2+i]-x[i]) < 1e-13*(1+fabs(x[i]))); rnew[2+n+i] = ub[i] - lb[i]; } rnew[0] = rect_diameter(n, rnew+2+n, &p); - if (x_center) - rnew[1] = p.fmin; /* avoid computing f(center) twice */ - else - rnew[1] = function_eval(rnew+2, &p); + rnew[1] = function_eval(rnew+2, &p); if (!rb_tree_insert(&p.rtree, rnew)) { free(rnew); goto done; diff --git a/test/testopt.cpp b/test/testopt.cpp index e7c312f..1b7a754 100644 --- a/test/testopt.cpp +++ b/test/testopt.cpp @@ -17,7 +17,7 @@ static nlopt_algorithm algorithm = NLOPT_GLOBAL_DIRECT; static double ftol_rel = 0, ftol_abs = 0, xtol_rel = 0, xtol_abs = 0, fmin_max = -HUGE_VAL; -static int maxeval = 1000, iterations = 1; +static int maxeval = 1000, iterations = 1, center_start = 0; static double maxtime = 0.0; static double xinit_tol = -1; @@ -77,14 +77,19 @@ static int test_function(int ifunc) printf("Starting guess x = ["); for (i = 0; i < func.n; ++i) { - if (xinit_tol < 0) { /* random starting point near center of box */ + if (center_start) + x[i] = (func.ub[i] + func.lb[i]) * 0.5; + else if (xinit_tol < 0) { /* random starting point near center of box */ double dx = (func.ub[i] - func.lb[i]) * 0.25; double xm = 0.5 * (func.ub[i] + func.lb[i]); x[i] = nlopt_urand(xm - dx, xm + dx); } - else + else { x[i] = nlopt_urand(-xinit_tol, xinit_tol) + (1 + nlopt_urand(-xinit_tol, xinit_tol)) * func.xmin[i]; + if (x[i] > func.ub[i]) x[i] = func.ub[i]; + else if (x[i] < func.lb[i]) x[i] = func.lb[i]; + } printf(" %g", x[i]); } printf("]\n"); @@ -143,6 +148,7 @@ static void usage(FILE *f) " -a : use optimization algorithm \n" " -o : use objective function \n" " -0 : starting guess within + (1+) * optimum\n" + " -c : starting guess at center of cell\n" " -e : use at most evals (default: %d, 0 to disable)\n" " -t : use at most seconds (default: disabled)\n" " -x : relative tolerance on x (default: disabled)\n" @@ -165,7 +171,7 @@ int main(int argc, char **argv) if (argc <= 1) usage(stdout); - while ((c = getopt(argc, argv, "hLv0:r:a:o:i:e:t:x:X:f:F:m:")) != -1) + while ((c = getopt(argc, argv, "hLvc0:r:a:o:i:e:t:x:X:f:F:m:")) != -1) switch (c) { case 'h': usage(stdout); @@ -217,7 +223,11 @@ int main(int argc, char **argv) case 'm': fmin_max = atof(optarg); break; + case 'c': + center_start = 1; + break; case '0': + center_start = 0; xinit_tol = atof(optarg); break; default: -- 2.30.2