From: stevenj Date: Sat, 25 Aug 2007 16:01:46 +0000 (-0400) Subject: recenter coords so that starting guess is utilized by global routines X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=3e978154f3ffa5faf3e575bd7931b2ea4af6836f;p=nlopt.git recenter coords so that starting guess is utilized by global routines darcs-hash:20070825160146-c8de0-27bd7e71a2e8ddcc6ab8e73c21f05e9e98a5dc8e.gz --- diff --git a/api/nlopt.c b/api/nlopt.c index 529c9eb..7b0e097 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -1,5 +1,6 @@ #include #include +#include #include "nlopt.h" #include "nlopt-util.h" @@ -71,9 +72,75 @@ void nlopt_srand_time() { typedef struct { nlopt_func f; void *f_data; - const double *lb, *ub; + const double *lb, *ub, *x0; + double *xtmp; } 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_) @@ -97,12 +164,15 @@ static double f_subplex(int n, const double *x, void *data_) static double f_direct(int n, const double *x, int *undefined, void *data_) { nlopt_data *data = (nlopt_data *) data_; - double f = data->f(n, x, NULL, data->f_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); *undefined = isnan(f) || my_isinf(f); return f; } #include "stogo.h" + #include "l-bfgs-b.h" /*************************************************************************/ @@ -126,6 +196,7 @@ 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) @@ -150,15 +221,22 @@ static nlopt_result nlopt_minimize_( switch (algorithm) { case NLOPT_GLOBAL_DIRECT: - case NLOPT_GLOBAL_DIRECT_L: - switch (direct_optimize(f_direct, &d, n, lb, ub, x, fmin, - maxeval, 500, ftol_rel, ftol_abs, - xtol_rel, xtol_rel, - DIRECT_UNKNOWN_FGLOBAL, -1.0, - NULL, - algorithm == NLOPT_GLOBAL_DIRECT - ? DIRECT_ORIGINAL - : DIRECT_GABLONSKY)) { + case NLOPT_GLOBAL_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, 500, ftol_rel, ftol_abs, + xtol_rel, xtol_rel, + DIRECT_UNKNOWN_FGLOBAL, -1.0, + NULL, + algorithm == NLOPT_GLOBAL_DIRECT + ? DIRECT_ORIGINAL + : DIRECT_GABLONSKY); + recenter_x(n, x, lb, ub, d.x0, x); + free(d.xtmp); + switch (iret) { case DIRECT_INVALID_BOUNDS: case DIRECT_MAXFEVAL_TOOBIG: case DIRECT_INVALID_ARGS: @@ -179,15 +257,23 @@ static nlopt_result nlopt_minimize_( return NLOPT_OUT_OF_MEMORY; } break; + } case NLOPT_GLOBAL_STOGO: - case NLOPT_GLOBAL_STOGO_RANDOMIZED: - if (!stogo_minimize(n, f, f_data, x, fmin, lb, ub, - maxeval, maxtime, - algorithm == NLOPT_GLOBAL_STOGO - ? 0 : 2*n)) - return NLOPT_FAILURE; + 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, + maxeval, maxtime, + algorithm == NLOPT_GLOBAL_STOGO + ? 0 : 2*n); + recenter_x(n, x, lb, ub, d.x0, x); + free(d.xtmp); + if (!iret) return NLOPT_FAILURE; break; + } case NLOPT_LOCAL_SUBPLEX: { int iret; diff --git a/api/nlopt_minimize.3 b/api/nlopt_minimize.3 index 312b658..a217805 100644 --- a/api/nlopt_minimize.3 +++ b/api/nlopt_minimize.3 @@ -42,10 +42,13 @@ design variables using the specified .IR algorithm . The minimum function value found is returned in .IR fmin , -with the corresponding design variable values stored in the array +with the corresponding design variable values returned in the array .I x of length .IR n . +The input values in +.I x +should be a starting guess for the optimum. The inputs .I lb and @@ -67,9 +70,7 @@ require the gradient (derivatives) of the function to be supplied via .IR f , and other algorithms do not require derivatives. Some of the algorithms attempt to find a global minimum within the given bounds, -and others find only a local minimum (with the initial value of -.I x -as a starting guess). +and others find only a local minimum. .PP The .B nlopt_minimize diff --git a/stogo/global.cc b/stogo/global.cc index 0b6c3bf..2f61703 100644 --- a/stogo/global.cc +++ b/stogo/global.cc @@ -64,10 +64,8 @@ void Global::FillRegular(RTBox SampleBox, RTBox box) { RVector m(dim), x(dim); if (det_pnts>0) { - // Add midpoint box.Midpoint(m) ; - tmpTrial.xvals=m ; tmpTrial.objval=DBL_MAX ; - SampleBox.AddTrial(tmpTrial) ; + tmpTrial.objval=DBL_MAX ; // Add the rest i=1 ; flag=1 ; dir=0 ; x=m ; @@ -83,6 +81,9 @@ void Global::FillRegular(RTBox SampleBox, RTBox box) { } i++ ; } + // Add midpoint + tmpTrial.xvals=m ; + SampleBox.AddTrial(tmpTrial) ; } } @@ -107,8 +108,8 @@ double Global::NewtonTest(RTBox box, int axis, RCRVector x_av, int *noutside) { double maxgrad=0 ; // Create sampling points - FillRegular(SampleBox, box); FillRandom(SampleBox, box); + FillRegular(SampleBox, box); // Perform the actual sampling while ( !SampleBox.EmptyBox() ) { diff --git a/test/testopt.cpp b/test/testopt.cpp index a723367..e7c312f 100644 --- a/test/testopt.cpp +++ b/test/testopt.cpp @@ -17,8 +17,9 @@ 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; +static int maxeval = 1000, iterations = 1; static double maxtime = 0.0; +static double xinit_tol = -1; static void listalgs(FILE *f) { @@ -39,7 +40,7 @@ static void listfuncs(FILE *f) static int test_function(int ifunc) { testfunc func; - int i; + int i, iter; double *x, fmin, f0, *xtabs; nlopt_result ret; double start = nlopt_seconds(); @@ -66,9 +67,26 @@ static int test_function(int ifunc) printf("-----------------------------------------------------------\n"); printf("Optimizing %s (%d dims) using %s algorithm\n", func.name, func.n, nlopt_algorithm_name(algorithm)); + printf("lower bounds at lb = ["); + for (i = 0; i < func.n; ++i) printf(" %g", func.lb[i]); + printf("]\n"); + printf("upper bounds at ub = ["); + for (i = 0; i < func.n; ++i) printf(" %g", func.ub[i]); + printf("]\n"); + + printf("Starting guess x = ["); - for (i = 0; i < func.n; ++i) - printf(" %g", x[i] = nlopt_urand(func.lb[i], func.ub[i])); + for (i = 0; i < func.n; ++i) { + 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 + x[i] = nlopt_urand(-xinit_tol, xinit_tol) + + (1 + nlopt_urand(-xinit_tol, xinit_tol)) * func.xmin[i]; + printf(" %g", x[i]); + } printf("]\n"); f0 = func.f(func.n, x, x + func.n, func.f_data); printf("Starting function value = %g\n", f0); @@ -85,29 +103,31 @@ static int test_function(int ifunc) } } - testfuncs_counter = 0; - ret = nlopt_minimize(algorithm, - func.n, func.f, func.f_data, - func.lb, func.ub, - x, &fmin, - fmin_max, ftol_rel, ftol_abs, xtol_rel, xtabs, - maxeval, maxtime); - printf("finished after %g seconds.\n", nlopt_seconds() - start); - printf("return code %d from nlopt_minimize\n", ret); - if (ret < 0) { - fprintf(stderr, "testopt: error in nlopt_minimize\n"); - return 0; + for (iter = 0; iter < iterations; ++iter) { + testfuncs_counter = 0; + ret = nlopt_minimize(algorithm, + func.n, func.f, func.f_data, + func.lb, func.ub, + x, &fmin, + fmin_max, ftol_rel, ftol_abs, xtol_rel, xtabs, + maxeval, maxtime); + printf("finished after %g seconds.\n", nlopt_seconds() - start); + printf("return code %d from nlopt_minimize\n", ret); + if (ret < 0) { + fprintf(stderr, "testopt: error in nlopt_minimize\n"); + return 0; + } + printf("Found minimum f = %g after %d evaluations.\n", + fmin, testfuncs_counter); + printf("Minimum at x = ["); + for (i = 0; i < func.n; ++i) printf(" %g", x[i]); + printf("]\n"); + printf("|f - fmin| = %g, |f - fmin| / |fmin| = %e\n", + fabs(fmin - func.fmin), fabs(fmin - func.fmin) / fabs(func.fmin)); } - printf("Found minimum f = %g after %d evaluations.\n", - fmin, testfuncs_counter); - printf("Minimum at x = ["); - for (i = 0; i < func.n; ++i) printf(" %g", x[i]); - printf("]\n"); printf("vs. global minimum f = %g at x = [", func.fmin); for (i = 0; i < func.n; ++i) printf(" %g", func.xmin[i]); printf("]\n"); - printf("|f - fmin| = %g, |f - fmin| / |fmin| = %e\n", - fabs(fmin - func.fmin), fabs(fmin - func.fmin) / fabs(func.fmin)); free(x); return 1; @@ -120,9 +140,9 @@ static void usage(FILE *f) " -h : print this help\n" " -L : list available algorithms and objective functions\n" " -v : verbose mode\n" - " -r : use random seed for starting guesses\n" " -a : use optimization algorithm \n" " -o : use objective function \n" + " -0 : starting guess within + (1+) * optimum\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" @@ -130,6 +150,8 @@ static void usage(FILE *f) " -f : relative tolerance on f (default: disabled)\n" " -F : absolute tolerance on f (default: disabled)\n" " -m : minimize f until is reached (default: disabled)\n" + " -i : iterate optimization times (default: 1)\n" + " -r : use random seed for starting guesses\n" , maxeval); } @@ -143,7 +165,7 @@ int main(int argc, char **argv) if (argc <= 1) usage(stdout); - while ((c = getopt(argc, argv, "hLvr:a:o:e:t:x:X:f:F:m:")) != -1) + while ((c = getopt(argc, argv, "hLv0:r:a:o:i:e:t:x:X:f:F:m:")) != -1) switch (c) { case 'h': usage(stdout); @@ -174,6 +196,9 @@ int main(int argc, char **argv) case 'e': maxeval = atoi(optarg); break; + case 'i': + iterations = atoi(optarg); + break; case 't': maxtime = atof(optarg); break; @@ -192,6 +217,9 @@ int main(int argc, char **argv) case 'm': fmin_max = atof(optarg); break; + case '0': + xinit_tol = atof(optarg); + break; default: fprintf(stderr, "harminv: invalid argument -%c\n", c); usage(stderr);