From: stevenj Date: Wed, 29 Aug 2007 23:41:11 +0000 (-0400) Subject: support calling both original and new DIRECT code, make direct limits more dynamic... X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=8097fc2c6df53379aacad22f1af223d5d62e64db;p=nlopt.git support calling both original and new DIRECT code, make direct limits more dynamic, eps != ftol in direct, fixed bugs in convex hull code for identical points darcs-hash:20070829234111-c8de0-2ce2867d6444a75b3434ec6b971c0aafefd3d6b2.gz --- diff --git a/api/nlopt.c b/api/nlopt.c index 1e89fad..00ac174 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -43,6 +43,7 @@ void nlopt_version(int *major, int *minor, int *bugfix) static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = { "DIRECT (global)", "DIRECT-L (global)", + "Randomized DIRECT-L (global)", "Original DIRECT version (global)", "Original DIRECT-L version (global)", "Subplex (local)", @@ -230,8 +231,10 @@ static nlopt_result nlopt_minimize_( switch (algorithm) { case NLOPT_GLOBAL_DIRECT: case NLOPT_GLOBAL_DIRECT_L: - return cdirect(n, f, f_data, lb, ub, x, fmin, &stop, 1e-4, - algorithm == NLOPT_GLOBAL_DIRECT ? 0 : 1); + case NLOPT_GLOBAL_DIRECT_L_RANDOMIZED: + 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: @@ -241,9 +244,9 @@ static nlopt_result nlopt_minimize_( 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, 1e-4, 1e-4, - xtol_rel, xtol_rel, - DIRECT_UNKNOWN_FGLOBAL, -1.0, + maxeval, 999, 0.0, 0.0, + pow(xtol_rel, (double) n), -1.0, + stop.fmin_max, 0.0, NULL, algorithm == NLOPT_GLOBAL_ORIG_DIRECT ? DIRECT_ORIGINAL @@ -263,7 +266,7 @@ static nlopt_result nlopt_minimize_( case DIRECT_MAXITER_EXCEEDED: return NLOPT_MAXEVAL_REACHED; case DIRECT_GLOBAL_FOUND: - return NLOPT_SUCCESS; + return NLOPT_FMIN_MAX_REACHED; case DIRECT_VOLTOL: case DIRECT_SIGMATOL: return NLOPT_XTOL_REACHED; diff --git a/api/nlopt.h b/api/nlopt.h index 516931a..8f9e30d 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -15,6 +15,7 @@ typedef enum { NLOPT_GLOBAL_DIRECT = 0, NLOPT_GLOBAL_DIRECT_L, + NLOPT_GLOBAL_DIRECT_L_RANDOMIZED, NLOPT_GLOBAL_ORIG_DIRECT, NLOPT_GLOBAL_ORIG_DIRECT_L, diff --git a/cdirect/cdirect.c b/cdirect/cdirect.c index bd79997..aa21b96 100644 --- a/cdirect/cdirect.c +++ b/cdirect/cdirect.c @@ -245,7 +245,10 @@ static int convex_hull(rb_tree *t, double **hull) yminmin = n->k[1]; xmax = nmax->k[0]; - hull[nhull++] = n->k; + do { /* include any duplicate points at (xmin,yminmin) */ + hull[nhull++] = n->k; + n = rb_tree_succ(n); + } while (n && n->k[0] == xmin && n->k[1] == yminmin); if (xmin == xmax) return nhull; /* set nmax = min mode with x == xmax */ @@ -327,7 +330,12 @@ static int convex_hull(rb_tree *t, double **hull) } hull[nhull++] = k; } - hull[nhull++] = nmax->k; + + do { /* include any duplicate points at (xmax,ymaxmin) */ + hull[nhull++] = nmax->k; + nmax = rb_tree_succ(nmax); + } while (nmax && nmax->k[0] == xmax && n->k[1] == ymaxmin); + return nhull; } @@ -359,16 +367,19 @@ static nlopt_result divide_good_rects(params *p) divisions: for (i = 0; i < nhull; ++i) { double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K; - if (i > 0) - K1 = (hull[i][1] - hull[i-1][1]) / (hull[i][0] - hull[i-1][0]); - if (i < nhull-1) - K1 = (hull[i][1] - hull[i+1][1]) / (hull[i][0] - hull[i+1][0]); + int im, ip; + for (im = i-1; im >= 0 && hull[im][0] == hull[i][0]; --im); + for (ip = i+1; ip < nhull && hull[ip][0] == hull[i][0]; ++ip); + if (im >= 0) + K1 = (hull[i][1] - hull[im][1]) / (hull[i][0] - hull[im][0]); + if (ip < nhull) + K1 = (hull[i][1] - hull[ip][1]) / (hull[i][0] - hull[ip][0]); K = MAX(K1, K2); if (hull[i][1] - K * hull[i][0] <= p->fmin - magic_eps * fabs(p->fmin)) { /* "potentially optimal" rectangle, so subdivide */ - divided_some = 1; nlopt_result ret = divide_rect(hull[i], p); + divided_some = 1; if (ret != NLOPT_SUCCESS) return ret; xtol_reached = xtol_reached && small(hull[i] + 2+n, p); } @@ -423,8 +434,8 @@ nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data, nlopt_result ret = NLOPT_OUT_OF_MEMORY; p.magic_eps = magic_eps; - p.which_diam = which_alg % 2; - p.which_div = 0; + p.which_diam = which_alg % 10; + p.which_div = (which_alg / 10) % 10; p.lb = lb; p.ub = ub; p.stop = stop; p.n = n; diff --git a/direct/DIRect.c b/direct/DIRect.c index 60d43aa..994ee05 100644 --- a/direct/DIRect.c +++ b/direct/DIRect.c @@ -54,7 +54,7 @@ doublereal d__1; /* constants (FIXME: change to variable?) */ - integer MAXFUNC = *maxf <= 0 ? 100050 : *maxf + 50; + integer MAXFUNC = *maxf <= 0 ? 101000 : (*maxf + 1000 + *maxf / 2); const integer MAXDEEP = 1000; const integer MAXDIV = 5000; @@ -99,6 +99,7 @@ MY_ALLOC(c__, doublereal, MAXFUNC * (*n)); MY_ALLOC(length, integer, MAXFUNC * (*n)); MY_ALLOC(f, doublereal, MAXFUNC * 2); + if (*maxf <= 0) *maxf = MAXFUNC - 1000; MY_ALLOC(s, integer, MAXDIV * 2); MY_ALLOC(w, doublereal, (*n)); diff --git a/direct/direct.h b/direct/direct.h index 4ceda34..8b6f916 100644 --- a/direct/direct.h +++ b/direct/direct.h @@ -44,7 +44,7 @@ extern direct_return_code direct_optimize( double *x, double *fmin, int max_feval, int max_iter, - double reltol, double abstol, + double magic_eps, double magic_eps_abs, double volume_reltol, double sigma_reltol, double fglobal, diff --git a/direct/direct_wrap.c b/direct/direct_wrap.c index 406a44d..3a26838 100644 --- a/direct/direct_wrap.c +++ b/direct/direct_wrap.c @@ -20,8 +20,12 @@ x: an array of length dimension, set to optimum variables upon return fmin: on return, set to minimum f value + magic_eps, magic_eps_abs: Jones' "magic" epsilon parameter, and + also an absolute version of the same + (not multipled by fmin). Jones suggests + setting this to 1e-4, but 0 also works... + max_feval, max_iter: maximum number of function evaluations & DIRECT iters - reltol, abstol: relative and absolute tolerances (0 if none) volume_reltol: relative tolerance on hypercube volume (0 if none) sigma_reltol: relative tolerance on hypercube "measure" (??) (0 if none) @@ -44,7 +48,7 @@ direct_return_code direct_optimize( double *x, double *fmin, int max_feval, int max_iter, - double reltol, double abstol, + double magic_eps, double magic_eps_abs, double volume_reltol, double sigma_reltol, double fglobal, @@ -80,7 +84,7 @@ direct_return_code direct_optimize( u[i] = upper_bounds[i]; } - direct_direct_(f, x, &dimension, &reltol, abstol, + direct_direct_(f, x, &dimension, &magic_eps, magic_eps_abs, &max_feval, &max_iter, fmin, l, u,