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

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

diff --git a/api/nlopt.c b/api/nlopt.c
index 0ed5bcb8c96a1e6ad460218d7b44db751d5acd64..669dd9fe3d329ec059574688b2a349ab6d7629e9 100644 (file)
--- a/api/nlopt.c
+++ b/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;

diff --git a/cdirect/cdirect.c b/cdirect/cdirect.c
index 9d95f8f90655f74ca8dfabc7d36e4a544b3551fc..52b55fafbb9f8f612eaa83c903a0e0ee423b56ce 100644 (file)
--- a/cdirect/cdirect.c
+++ b/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;

diff --git a/test/testopt.cpp b/test/testopt.cpp
index e7c312f86f1b0c8dc0e07d988ffce162534178fd..1b7a7542d53da1a28619b344ce9ed78ca4652f9b 100644 (file)
--- a/test/testopt.cpp
+++ b/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: