recenter coords so that starting guess is utilized by global routines

[nlopt.git] / api / nlopt.c

diff --git a/api/nlopt.c b/api/nlopt.c
index 529c9eb2d057a1c6e17ac0e93314f58a0b75a493..7b0e09772447c84388ed016fae6f8179bed7faba 100644 (file)
--- a/api/nlopt.c
+++ b/api/nlopt.c
@@ -1,5 +1,6 @@
 #include <stdlib.h>
 #include <math.h>
+#include <string.h>
 
 #include "nlopt.h"
 #include "nlopt-util.h"
@@ -71,9 +72,75 @@ void nlopt_srand_time() {
 typedef struct {
      nlopt_func f;
      void *f_data;
-     const double *lb, *ub;
+     const double *lb, *ub, *x0;
+     double *xtmp;
 } 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_)
@@ -97,12 +164,15 @@ static double f_subplex(int n, const double *x, void *data_)
 static double f_direct(int n, const double *x, int *undefined, void *data_)
 {
      nlopt_data *data = (nlopt_data *) data_;
-     double f = data->f(n, x, NULL, data->f_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);
      *undefined = isnan(f) || my_isinf(f);
      return f;
 }
 
 #include "stogo.h"
+
 #include "l-bfgs-b.h"
 
 /*************************************************************************/
@@ -126,6 +196,7 @@ 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)
@@ -150,15 +221,22 @@ static nlopt_result nlopt_minimize_(
 
      switch (algorithm) {
         case NLOPT_GLOBAL_DIRECT:
-        case NLOPT_GLOBAL_DIRECT_L:
-             switch (direct_optimize(f_direct, &d, n, lb, ub, x, fmin,
-                                     maxeval, 500, ftol_rel, ftol_abs,
-                                     xtol_rel, xtol_rel,
-                                     DIRECT_UNKNOWN_FGLOBAL, -1.0,
-                                     NULL, 
-                                     algorithm == NLOPT_GLOBAL_DIRECT
-                                     ? DIRECT_ORIGINAL
-                                     : DIRECT_GABLONSKY)) {
+        case NLOPT_GLOBAL_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, 500, ftol_rel, ftol_abs,
+                                    xtol_rel, xtol_rel,
+                                    DIRECT_UNKNOWN_FGLOBAL, -1.0,
+                                    NULL, 
+                                    algorithm == NLOPT_GLOBAL_DIRECT
+                                    ? DIRECT_ORIGINAL
+                                    : DIRECT_GABLONSKY);
+             recenter_x(n, x, lb, ub, d.x0, x);
+             free(d.xtmp);
+             switch (iret) {
                  case DIRECT_INVALID_BOUNDS:
                  case DIRECT_MAXFEVAL_TOOBIG:
                  case DIRECT_INVALID_ARGS:
@@ -179,15 +257,23 @@ static nlopt_result nlopt_minimize_(
                       return NLOPT_OUT_OF_MEMORY;
              }
              break;
+        }
 
         case NLOPT_GLOBAL_STOGO:
-        case NLOPT_GLOBAL_STOGO_RANDOMIZED:
-             if (!stogo_minimize(n, f, f_data, x, fmin, lb, ub,
-                                 maxeval, maxtime,
-                                 algorithm == NLOPT_GLOBAL_STOGO
-                                 ? 0 : 2*n))
-                  return NLOPT_FAILURE;
+        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,
+                                   maxeval, maxtime,
+                                   algorithm == NLOPT_GLOBAL_STOGO
+                                   ? 0 : 2*n);
+             recenter_x(n, x, lb, ub, d.x0, x);
+             free(d.xtmp);
+             if (!iret) return NLOPT_FAILURE;
              break;
+        }
 
         case NLOPT_LOCAL_SUBPLEX: {
              int iret;