got rid of buggy f_recenter (Lagrange interpolation isn't guaranteed to be monotonic)

darcs-hash:20070830010720-c8de0-19c1e18baffbe16a121a8c968f2c0d2a31f27058.gz

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;

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;

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 <n> : use optimization algorithm <n>\n"
 	  " -o <n> : use objective function <n>\n"
 	  " -0 <x> : starting guess within <x> + (1+<x>) * optimum\n"
+	  "     -c : starting guess at center of cell\n"
 	  " -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"
@@ -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: