diff --git a/api/nlopt.c b/api/nlopt.c index ed2e6711d9155725181dfc32c01d05a975dd295f..bf7ee9f5ffb6b53961505fcdcba9fbc0fdbcf45d 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -106,7 +106,8 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) /*************************************************************************/ -nlopt_result nlopt_minimize( +/* same as nlopt_minimize, but xtol_abs is required to be non-NULL */ +static nlopt_result nlopt_minimize_( nlopt_algorithm algorithm, int n, nlopt_func f, void *f_data, const double *lb, const double *ub, /* bounds */ @@ -118,6 +119,8 @@ nlopt_result nlopt_minimize( { int i; nlopt_data d; + nlopt_stopping stop; + d.f = f; d.f_data = f_data; d.lb = lb; @@ -132,6 +135,18 @@ nlopt_result nlopt_minimize( if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i]) return NLOPT_INVALID_ARGS; + stop.n = n; + stop.fmin_max = (isnan(fmin_max) || (my_isinf(fmin_max) && fmin_max < 0)) + ? -MY_INF : fmin_max; + stop.ftol_rel = ftol_rel; + stop.ftol_abs = ftol_abs; + stop.xtol_rel = xtol_rel; + stop.xtol_abs = xtol_abs; + stop.nevals = 0; + stop.maxeval = maxeval; + stop.maxtime = maxtime; + stop.start = nlopt_seconds(); + switch (algorithm) { case NLOPT_GLOBAL_DIRECT: case NLOPT_GLOBAL_DIRECT_L: @@ -184,15 +199,18 @@ nlopt_result nlopt_minimize( else scale[i] = 0.01 * x[i] + 0.0001; } - iret = subplex(f_subplex, fmin, x, n, &d, xtol_rel, maxeval, - fmin_max, !my_isinf(fmin_max), scale); + iret = subplex(f_subplex, fmin, x, n, &d, &stop, scale); free(scale); switch (iret) { case -2: return NLOPT_INVALID_ARGS; + case -10: return NLOPT_MAXTIME_REACHED; case -1: return NLOPT_MAXEVAL_REACHED; case 0: return NLOPT_XTOL_REACHED; case 1: return NLOPT_SUCCESS; case 2: return NLOPT_FMIN_MAX_REACHED; + case 20: return NLOPT_FTOL_REACHED; + case -200: return NLOPT_OUT_OF_MEMORY; + default: return NLOPT_FAILURE; /* unknown return code */ } break; } @@ -234,3 +252,31 @@ nlopt_result nlopt_minimize( return NLOPT_SUCCESS; } + +nlopt_result nlopt_minimize( + nlopt_algorithm algorithm, + int n, nlopt_func f, void *f_data, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *fmin, /* out: minimum */ + double fmin_max, double ftol_rel, double ftol_abs, + double xtol_rel, const double *xtol_abs, + int maxeval, double maxtime) +{ + nlopt_result ret; + if (xtol_abs) + ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub, + x, fmin, fmin_max, ftol_rel, ftol_abs, + xtol_rel, xtol_abs, maxeval, maxtime); + else { + int i; + double *xtol = (double *) malloc(sizeof(double) * n); + if (!xtol) return NLOPT_OUT_OF_MEMORY; + for (i = 0; i < n; ++i) xtol[i] = -1; + ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub, + x, fmin, fmin_max, ftol_rel, ftol_abs, + xtol_rel, xtol, maxeval, maxtime); + free(xtol); + } + return ret; +} diff --git a/api/nlopt_minimize.3 b/api/nlopt_minimize.3 index 6c165f59005532013965fcee047ea35613e62d56..b8b59c0eb573eed79d2fa27c9f9f628eb9481514 100644 --- a/api/nlopt_minimize.3 +++ b/api/nlopt_minimize.3 @@ -222,7 +222,7 @@ some algorithms do not support all of the stopping criteria. .B fmin_max Stop when a function value less than or equal to .I fmin_max -is found. Set to Inf (see constraints section above) to disable. +is found. Set to -Inf or NaN (see constraints section above) to disable. .TP .B ftol_rel Relative tolerance on function value: stop when an optimization step diff --git a/subplex/Makefile.am b/subplex/Makefile.am index 8ad39b62c022f8dff7711e686a4597f1d6660119..bae7d9c8e4da5c04a6e13e4df2eec2f91ec5a69f 100644 --- a/subplex/Makefile.am +++ b/subplex/Makefile.am @@ -1,3 +1,5 @@ +AM_CPPFLAGS = -I$(top_srcdir)/util + noinst_LTLIBRARIES = libsubplex.la libsubplex_la_SOURCES = subplex.c subplex.h diff --git a/subplex/subplex.c b/subplex/subplex.c index 28c302ef39b5e29b2bf2d94a3c48229559742b3f..17867c6b7225d54064be4ef4615a4f42b2463b07 100644 --- a/subplex/subplex.c +++ b/subplex/subplex.c @@ -58,6 +58,7 @@ COMMENTS #include #include #include +#include #include "subplex.h" @@ -1265,7 +1266,7 @@ static int evalf_(D_fp f,void*fdata, integer *ns, integer *ips, doublereal *xs, */ static int simplx_(D_fp f, void *fdata, integer *n, doublereal *step, integer * - ns, integer *ips, integer *maxnfe, logical *cmode, doublereal *x, + ns, integer *ips, nlopt_stopping *stop, logical *cmode, doublereal *x, doublereal *fx, integer *nfe, doublereal *s, doublereal *fs, integer * iflag) { @@ -1384,6 +1385,7 @@ static int simplx_(D_fp f, void *fdata, integer *n, doublereal *step, integer * if (usubc_1.irepl > 0) { isubc_1.new__ = FALSE_; evalf_((D_fp)f,fdata, ns, &ips[1], &s[s_dim1 + 1], n, &x[1], &fs[1], nfe); + stop->nevals++; } else { fs[1] = *fx; } @@ -1392,6 +1394,7 @@ static int simplx_(D_fp f, void *fdata, integer *n, doublereal *step, integer * for (j = 2; j <= i__1; ++j) { evalf_((D_fp)f, fdata,ns, &ips[1], &s[j * s_dim1 + 1], n, &x[1], &fs[j], nfe); + stop->nevals++; /* L10: */ } il = 1; @@ -1413,6 +1416,7 @@ L20: goto L40; } evalf_((D_fp)f,fdata, ns, &ips[1], &s[itemp * s_dim1 + 1], n, &x[1], &fr, nfe); + stop->nevals++; if (fr < fs[il]) { /* expand */ @@ -1424,6 +1428,7 @@ L20: goto L40; } evalf_((D_fp)f,fdata, ns, &ips[1], &s[ih * s_dim1 + 1], n, &x[1], &fe, nfe); + stop->nevals++; if (fe < fr) { fs[ih] = fe; } else { @@ -1455,6 +1460,7 @@ L20: } evalf_((D_fp)f,fdata, ns, &ips[1], &s[itemp * s_dim1 + 1], n, &x[1], &fc, nfe); + stop->nevals++; /* Computing MIN */ d__1 = fr, d__2 = fs[ih]; if (fc < min(d__1,d__2)) { @@ -1476,6 +1482,7 @@ L20: } evalf_((D_fp)f,fdata, ns, &ips[1], &s[j * s_dim1 + 1], n, &x[1], &fs[j], nfe); + stop->nevals++; } /* L30: */ } @@ -1493,17 +1500,17 @@ L40: *fx = isubc_1.sfbest; } L50: - if (usubc_1.nfstop > 0 && *fx <= isubc_1.sfstop && usubc_1.nfxe >= - usubc_1.nfstop) { - *iflag = 2; - } else if (*nfe >= *maxnfe) { - *iflag = -1; - } else if (dist_(ns, &s[ih * s_dim1 + 1], &s[il * s_dim1 + 1]) <= tol || - small) { - *iflag = 0; - } else { - goto L20; - } + if (*fx < stop->fmin_max) + *iflag = 2; + else if (nlopt_stop_evals(stop)) + *iflag = -1; + else if (nlopt_stop_time(stop)) + *iflag = -10; + else if (dist_(ns, &s[ih * s_dim1 + 1], &s[il * s_dim1 + 1]) <= tol + || small) + *iflag = 0; + else + goto L20; /* end main loop, return best point */ @@ -1790,9 +1797,11 @@ static int setstp_(integer *nsubs, integer *n, doublereal *deltax, -lf2c -lm (in that order) */ -static int subplx_(D_fp f, void *fdata, integer *n, doublereal *tol, integer * - maxnfe, integer *mode, const doublereal *scale, doublereal *x, doublereal * - fx, integer *nfe, doublereal *work, integer *iwork, integer *iflag) +static int subplx_(D_fp f, void *fdata, integer *n, + nlopt_stopping *stop, integer *mode, + const doublereal *scale, doublereal *x, doublereal *fx, + integer *nfe, doublereal *work, integer *iwork, + integer *iflag) { /* Initialized data */ @@ -1814,7 +1823,7 @@ static int subplx_(D_fp f, void *fdata, integer *n, doublereal *tol, integer * static integer istptr; static doublereal scl, dum; static integer ins; - static doublereal sfx; + static doublereal sfx, sfx_old, *x_old; @@ -1913,6 +1922,8 @@ static int subplx_(D_fp f, void *fdata, integer *n, doublereal *tol, integer * --scale; --work; --iwork; + x_old = work; + work += *n; /* Function Body */ /* ----------------------------------------------------------- */ @@ -1924,12 +1935,6 @@ static int subplx_(D_fp f, void *fdata, integer *n, doublereal *tol, integer * if (*n < 1) { goto L120; } - if (*tol < 0.) { - goto L120; - } - if (*maxnfe < 1) { - goto L120; - } if (scale[1] >= 0.) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { @@ -2031,6 +2036,7 @@ static int subplx_(D_fp f, void *fdata, integer *n, doublereal *tol, integer * isubc_1.new__ = TRUE_; usubc_1.initx = TRUE_; evalf_((D_fp)f, fdata, &c__0, &iwork[1], &dum, n, &x[1], &sfx, nfe); + stop->nevals++; usubc_1.initx = FALSE_; } else { @@ -2071,12 +2077,12 @@ L40: ipptr = 1; /* simplex loop */ - + sfx_old = sfx; + memcpy(&x_old[1], &x[1], sizeof(doublereal) * *n); L60: ns = iwork[ins]; L70: - simplx_((D_fp)f, fdata, n, &work[istptr], &ns, &iwork[ipptr], maxnfe, &cmode, &x[ - 1], &sfx, nfe, &work[isptr], &work[ifsptr], iflag); + simplx_((D_fp)f, fdata, n, &work[istptr], &ns, &iwork[ipptr], stop, &cmode, &x[1], &sfx, nfe, &work[isptr], &work[ifsptr], iflag); cmode = FALSE_; if (*iflag != 0) { goto L110; @@ -2097,6 +2103,15 @@ L70: /* check termination */ + if (nlopt_stop_ftol(stop, sfx, sfx_old)) { + *iflag = 20; + goto L110; + } + if (nlopt_stop_x(stop, &x[1], &x_old[1])) { + *iflag = 0; + goto L110; + } + L90: istep = istptr; i__1 = *n; @@ -2106,7 +2121,7 @@ L90: d__1)) * usubc_1.psi; /* Computing MAX */ d__6 = (d__3 = x[i__], abs(d__3)); - if (max(d__4,d__5) / max(d__6,1.) > *tol) { + if (max(d__4,d__5) / max(d__6,1.) > stop->xtol_rel) { setstp_(&nsubs, n, &work[1], &work[istptr]); goto L40; } @@ -2147,47 +2162,49 @@ L120: fmin: (output) value of f at minimum x[n]: (input) starting guess position, (output) computed minimum fdata: data pointer passed to f + + old args: tol: relative error tolerance for x maxnfe: maximum number of function evaluations fmin_max, use_fmin_max: if use_fmin_max, stop when f <= fmin_max + + new args: nlopt_stopping *stop (stopping criteria) + scale[n]: (input) scale & initial stepsizes for components of x (if *scale < 0, |*scale| is used for all components) Return value: = -2 : invalid input - = -1 : maxnfe exceeded + = -10 : maxtime exceeded + = -1 : maxevals exceeded = 0 : tol satisfied = 1 : limit of machine precision = 2 : fstop reached (fstop usage is determined by values of options minf, nfstop, and irepl. default: f(x) not tested against fstop) + = 20 : ftol reached + = -200 : out of memory */ int subplex(subplex_func f, double *fmin, double *x, int n, void *fdata, - double tol, int maxnfe, - double fmin_max, int use_fmin_max, - const double *scale) + nlopt_stopping *stop, + const double *scale) { int mode = 0, *iwork, nsmax, nsmin, errflag, nfe; double *work; nsmax = min(5,n); nsmin = min(2,n); - work = (double*) malloc(sizeof(double) * (2*n + nsmax*(nsmax+4) + 1)); + work = (double*) malloc(sizeof(double) * (3*n + nsmax*(nsmax+4) + 1)); + if (!work) + return -200; iwork = (int*) malloc(sizeof(int) * (n + n/nsmin + 1)); - if (!work || !iwork) { - fprintf(stderr, "subplex: error, out of memory!\n"); - exit(EXIT_FAILURE); - } - - if (use_fmin_max) { /* stop when fmin_max is reached */ - subopt_(&n); - usubc_1.nfstop = 1; - usubc_1.fstop = fmin_max; - mode = 2; + if (!iwork) { + free(work); + return -200; } subplx_(f,fdata, &n, - &tol, &maxnfe, &mode, + stop, &mode, scale, x, fmin, &nfe, work, iwork, &errflag); diff --git a/subplex/subplex.h b/subplex/subplex.h index 0fccb7914c133aa51e91b739df1cba4d06af7c43..f9f5b2ff8fa6e0157973ed8e3f081e6955903d4e 100644 --- a/subplex/subplex.h +++ b/subplex/subplex.h @@ -1,6 +1,8 @@ #ifndef SUBPLEX_H #define SUBPLEX_H +#include "nlopt-util.h" + #ifdef __cplusplus extern "C" { @@ -9,8 +11,7 @@ extern "C" typedef double (*subplex_func)(int n, const double *x, void *func_data); extern int subplex(subplex_func f, double *fmin, double *x, int n, void *fdata, - double tol, int maxnfe, - double fmin_max, int use_fmin_max, + nlopt_stopping *stop, const double *scale); #ifdef __cplusplus diff --git a/test/testopt.cpp b/test/testopt.cpp index eb8e4db2cc8b2ead63e8ef3515a6ac16355adc0f..a723367d55cf89810499c2edc4ca86e7e0a792f7 100644 --- a/test/testopt.cpp +++ b/test/testopt.cpp @@ -16,7 +16,7 @@ #include "testfuncs.h" static nlopt_algorithm algorithm = NLOPT_GLOBAL_DIRECT; -static double ftol_rel = 0, ftol_abs = 0, xtol_rel = 0; +static double ftol_rel = 0, ftol_abs = 0, xtol_rel = 0, xtol_abs = 0, fmin_max = -HUGE_VAL; static int maxeval = 1000; static double maxtime = 0.0; @@ -40,7 +40,7 @@ static int test_function(int ifunc) { testfunc func; int i; - double *x, fmin, f0; + double *x, fmin, f0, *xtabs; nlopt_result ret; double start = nlopt_seconds(); @@ -56,9 +56,12 @@ static int test_function(int ifunc) nlopt_algorithm_name(algorithm)); return 0; } - x = (double *) malloc(sizeof(double) * func.n * 2); + x = (double *) malloc(sizeof(double) * func.n * 3); if (!x) { fprintf(stderr, "testopt: Out of memory!\n"); return 0; } + xtabs = x + func.n * 2; + for (i = 0; i < func.n; ++i) xtabs[i] = xtol_abs; + printf("-----------------------------------------------------------\n"); printf("Optimizing %s (%d dims) using %s algorithm\n", @@ -87,7 +90,7 @@ static int test_function(int ifunc) func.n, func.f, func.f_data, func.lb, func.ub, x, &fmin, - HUGE_VAL, ftol_rel, ftol_abs, xtol_rel, NULL, + 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); @@ -120,8 +123,14 @@ static void usage(FILE *f) " -r : use random seed for starting guesses\n" " -a : use optimization algorithm \n" " -o : use objective function \n" - " -e : use at most evals (default: %d)\n", - maxeval); + " -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" + " -X : absolute tolerance on x (default: disabled)\n" + " -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" + , maxeval); } int main(int argc, char **argv) @@ -134,7 +143,7 @@ int main(int argc, char **argv) if (argc <= 1) usage(stdout); - while ((c = getopt(argc, argv, "hLvra:o:e:")) != -1) + while ((c = getopt(argc, argv, "hLvr:a:o:e:t:x:X:f:F:m:")) != -1) switch (c) { case 'h': usage(stdout); @@ -147,7 +156,7 @@ int main(int argc, char **argv) testfuncs_verbose = 1; break; case 'r': - srand((unsigned) atoi(optarg)); + nlopt_srand((unsigned long) atoi(optarg)); break; case 'a': c = atoi(optarg); @@ -165,6 +174,24 @@ int main(int argc, char **argv) case 'e': maxeval = atoi(optarg); break; + case 't': + maxtime = atof(optarg); + break; + case 'x': + xtol_rel = atof(optarg); + break; + case 'X': + xtol_abs = atof(optarg); + break; + case 'f': + ftol_rel = atof(optarg); + break; + case 'F': + ftol_abs = atof(optarg); + break; + case 'm': + fmin_max = atof(optarg); + break; default: fprintf(stderr, "harminv: invalid argument -%c\n", c); usage(stderr); diff --git a/util/Makefile.am b/util/Makefile.am index 0bbe4b29573ece26a4804ec468a5d89fa509eac3..9d7214342b16e2b528b9e018cb614539afd6d689 100644 --- a/util/Makefile.am +++ b/util/Makefile.am @@ -1,2 +1,2 @@ noinst_LTLIBRARIES = libutil.la -libutil_la_SOURCES = mt19937ar.c timer.c nlopt-util.h +libutil_la_SOURCES = mt19937ar.c timer.c stop.c nlopt-util.h diff --git a/util/nlopt-util.h b/util/nlopt-util.h index 895db62f0a9346dc492e9dc76ad2a4f32fab1ee1..ea9d7cf309f43ec5489f9242a0948f12e9416e0f 100644 --- a/util/nlopt-util.h +++ b/util/nlopt-util.h @@ -14,6 +14,27 @@ extern unsigned long nlopt_time_seed(void); extern void nlopt_init_genrand(unsigned long s); extern double nlopt_urand(double a, double b); +/* stopping criteria */ +typedef struct { + int n; + double fmin_max; + double ftol_rel; + double ftol_abs; + double xtol_rel; + const double *xtol_abs; + int nevals, maxeval; + double maxtime, start; +} nlopt_stopping; +extern int nlopt_stop_f(const nlopt_stopping *stop, double f, double oldf); +extern int nlopt_stop_ftol(const nlopt_stopping *stop, double f, double oldf); +extern int nlopt_stop_x(const nlopt_stopping *stop, + const double *x, const double *oldx); +extern int nlopt_stop_xs(const nlopt_stopping *stop, + const double *xs, const double *oldxs, + const double *scale_min, const double *scale_max); +extern int nlopt_stop_evals(const nlopt_stopping *stop); +extern int nlopt_stop_time(const nlopt_stopping *stop); + #ifdef __cplusplus } /* extern "C" */ #endif /* __cplusplus */ diff --git a/util/stop.c b/util/stop.c new file mode 100644 index 0000000000000000000000000000000000000000..a48a83b5a0a99f9d13b5106b5184fc1b4d470d75 --- /dev/null +++ b/util/stop.c @@ -0,0 +1,59 @@ +#include +#include "nlopt-util.h" + +/* utility routines to implement the various stopping criteria */ + +static int relstop(double old, double new, double reltol, double abstol) +{ + return(fabs(new - old) < abstol + || fabs(new - old) < reltol * (fabs(new) + fabs(old)) * 0.5); +} + +int nlopt_stop_ftol(const nlopt_stopping *s, const double f, double oldf) +{ + return (relstop(oldf, f, s->ftol_rel, s->ftol_abs)); +} + +int nlopt_stop_f(const nlopt_stopping *s, const double f, double oldf) +{ + return (f <= s->fmin_max || nlopt_stop_ftol(s, f, oldf)); +} + +int nlopt_stop_x(const nlopt_stopping *s, const double *x, const double *oldx) +{ + int i; + for (i = 0; i < s->n; ++i) + if (relstop(oldx[i], x[i], s->xtol_rel, s->xtol_abs[i])) + return 1; + return 0; +} + +static double sc(double x, double smin, double smax) +{ + return smin + x * (smax - smin); +} + +/* some of the algorithms rescale x to a unit hypercube, so we need to + scale back before we can compare to the tolerances */ +int nlopt_stop_xs(const nlopt_stopping *s, + const double *xs, const double *oldxs, + const double *scale_min, const double *scale_max) +{ + int i; + for (i = 0; i < s->n; ++i) + if (relstop(sc(oldxs[i], scale_min[i], scale_max[i]), + sc(xs[i], scale_min[i], scale_max[i]), + s->xtol_rel, s->xtol_abs[i])) + return 1; + return 0; +} + +int nlopt_stop_evals(const nlopt_stopping *s) +{ + return (s->maxeval > 0 && s->nevals >= s->maxeval); +} + +int nlopt_stop_time(const nlopt_stopping *s) +{ + return (s->maxtime > 0 && nlopt_seconds() - s->start >= s->maxtime); +}