supported more stopping criteria in subplex (still a little too pessimistic)

darcs-hash:20070824211557-c8de0-f0f5566b654f5fd185b9459652ef0ea638d58707.gz

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;
+}

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

subplex/Makefile.am

+AM_CPPFLAGS = -I$(top_srcdir)/util
+
 noinst_LTLIBRARIES = libsubplex.la
 libsubplex_la_SOURCES = subplex.c subplex.h
 

subplex/subplex.c

@@ -58,6 +58,7 @@ COMMENTS
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include <string.h>
 
 #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);

subplex/subplex.h

 #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

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 <s> : use random seed <s> for starting guesses\n"
 	  " -a <n> : use optimization algorithm <n>\n"
 	  " -o <n> : use objective function <n>\n"
-	  " -e <n> : use at most <n> evals (default: %d)\n",
-	  maxeval);
+	  " -e <n> : use at most <n> evals (default: %d, 0 to disable)\n"
+	  " -t <t> : use at most <t> seconds (default: disabled)\n"
+	  " -x <t> : relative tolerance <t> on x (default: disabled)\n"
+	  " -X <t> : absolute tolerance <t> on x (default: disabled)\n"
+	  " -f <t> : relative tolerance <t> on f (default: disabled)\n"
+	  " -F <t> : absolute tolerance <t> on f (default: disabled)\n"
+	  " -m <m> : minimize f until <m> 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);

util/Makefile.am

 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

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 */

util/stop.c

+#include <math.h>
+#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);
+}