added nlopt_minimize_c function for constrained minimization (only via MMA for now...

darcs-hash:20080818204722-c8de0-3500ba9e001ac8d8f6a4593385c258ee6ed6ce4f.gz

diff --git a/api/nlopt.c b/api/nlopt.c
index b8019bf47c5b18f4bfb314bce1834f3ac9e43ed7..e10448fdc16072892bd912c6a777264bb6afc468 100644 (file)
--- a/api/nlopt.c
+++ b/api/nlopt.c
@@ -190,6 +190,7 @@ void nlopt_set_local_search_algorithm(nlopt_algorithm deriv,
 static nlopt_result nlopt_minimize_(
      nlopt_algorithm algorithm,
      int n, nlopt_func f, void *f_data,
+     int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size,
      const double *lb, const double *ub, /* bounds */
      double *x, /* in: initial guess, out: minimizer */
      double *minf, /* out: minimum */
@@ -210,6 +211,9 @@ static nlopt_result nlopt_minimize_(
      else if (n < 0 || !lb || !ub || !x)
          return NLOPT_INVALID_ARGS;
 
+     /* nonlinear constraints are only supported with MMA */
+     if (m != 0 && algorithm != NLOPT_LD_MMA) return NLOPT_INVALID_ARGS;
+
      d.f = f;
      d.f_data = f_data;
      d.lb = lb;
@@ -409,7 +413,7 @@ static nlopt_result nlopt_minimize_(
                                   algorithm >= NLOPT_GN_MLSL_LDS);
 
         case NLOPT_LD_MMA:
-             return mma_minimize(n, f, f_data, 0, NULL, NULL, 0,
+             return mma_minimize(n, f, f_data, m, fc, fc_data, fc_datum_size,
                                  lb, ub, x, minf, &stop,
                                  local_search_alg_deriv, 1e-8, 100000);
 
@@ -420,9 +424,10 @@ static nlopt_result nlopt_minimize_(
      return NLOPT_SUCCESS;
 }
 
-nlopt_result nlopt_minimize(
+nlopt_result nlopt_minimize_c(
      nlopt_algorithm algorithm,
      int n, nlopt_func f, void *f_data,
+     int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size,
      const double *lb, const double *ub, /* bounds */
      double *x, /* in: initial guess, out: minimizer */
      double *minf, /* out: minimum */
@@ -432,7 +437,8 @@ nlopt_result nlopt_minimize(
 {
      nlopt_result ret;
      if (xtol_abs)
-         ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub,
+         ret = nlopt_minimize_(algorithm, n, f, f_data,
+                               m, fc, fc_data, fc_datum_size, lb, ub,
                                x, minf, minf_max, ftol_rel, ftol_abs,
                                xtol_rel, xtol_abs, maxeval, maxtime);
      else {
@@ -440,10 +446,26 @@ nlopt_result nlopt_minimize(
          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,
+         ret = nlopt_minimize_(algorithm, n, f, f_data, 
+                               m, fc, fc_data, fc_datum_size, lb, ub,
                                x, minf, minf_max, ftol_rel, ftol_abs,
                                xtol_rel, xtol, maxeval, maxtime);
          free(xtol);
      }
      return ret;
 }
+
+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 *minf, /* out: minimum */
+     double minf_max, double ftol_rel, double ftol_abs,
+     double xtol_rel, const double *xtol_abs,
+     int maxeval, double maxtime)
+{
+     return nlopt_minimize_c(algorithm, n, f, f_data, 0, NULL, NULL, 0,
+                            lb, ub, x, minf, minf_max, ftol_rel, ftol_abs,
+                            xtol_rel, xtol_abs, maxeval, maxtime);
+}

diff --git a/api/nlopt.h b/api/nlopt.h
index edd04c41233f88f9d0bc92a45122765d0eac4a09..fa7629a6bbceeaed8bf25bb1268b7c2e3554242a 100644 (file)
--- a/api/nlopt.h
+++ b/api/nlopt.h
@@ -1,6 +1,8 @@
 #ifndef NLOPT_H
 #define NLOPT_H
 
+#include <stddef.h> /* for ptrdiff_t */
+
 #ifdef __cplusplus
 extern "C"
 {
@@ -89,6 +91,17 @@ extern nlopt_result nlopt_minimize(
      double xtol_rel, const double *xtol_abs,
      int maxeval, double maxtime);
 
+extern nlopt_result nlopt_minimize_c(
+     nlopt_algorithm algorithm,
+     int n, nlopt_func f, void *f_data,
+     int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size,
+     const double *lb, const double *ub, /* bounds */
+     double *x, /* in: initial guess, out: minimizer */
+     double *minf, /* out: minimum */
+     double minf_max, double ftol_rel, double ftol_abs,
+     double xtol_rel, const double *xtol_abs,
+     int maxeval, double maxtime);
+
 extern void nlopt_srand(unsigned long seed);
 extern void nlopt_srand_time(void);
 

diff --git a/mma/mma.c b/mma/mma.c
index 4c595c4b8fd2e58e707dba83c6374df32a4f0662..e91a979138f8f8269b36e04ae7e61626e3bbf2b3 100644 (file)
--- a/mma/mma.c
+++ b/mma/mma.c
@@ -30,6 +30,8 @@ typedef struct {
      double gval, wval, *gcval; /* output each time (array length m) */
 } dual_data;
 
+static double sqr(double x) { return x * x; }
+
 static double dual_func(int m, const double *y, double *grad, void *d_)
 {
      dual_data *d = (dual_data *) d_;
@@ -47,79 +49,56 @@ static double dual_func(int m, const double *y, double *grad, void *d_)
      val = d->gval = fval;
      d->wval = 0;
      for (i = 0; i < m; ++i) val += y[i] * (gcval[i] = fcval[i]);
-     if (grad) { for (i = 0; i < m; ++i) grad[i] = fcval[i]; }
 
      for (j = 0; j < n; ++j) {
-         int x_pinned;
-         double a, b, dx, denom, denominv, ca, cb, sigma2;
+         double u, v, dx, denominv, c, sigma2, dx2;
+
+         /* first, compute xcur[j] for y.  Because this objective is
+            separable, we can minimize over x analytically, and the minimum
+            dx is given by the solution of a quadratic equation:
+                    u dx^2 + 2 v sigma^2 dx + u sigma^2 = 0
+            where u and v are defined by the sums below.  Because of
+            the definitions, it is guaranteed that |u/v| <= sigma,
+            and it follows that the only dx solution with |dx| <= sigma
+            is given by:
+                    (v/u) sigma^2 (-1 + sqrt(1 - (u / v sigma)^2))
+             (which goes to zero as u -> 0). */
 
-         /* first, compute xcur[j] for y */
-         a = dfdx[j];
-         b = fabs(dfdx[j]) * sigma[j] + 0.5 * rho;
+         u = dfdx[j];
+         v = fabs(dfdx[j]) * sigma[j] + 0.5 * rho;
          for (i = 0; i < m; ++i) {
-              a += dfcdx[i*n + j] * y[i];
-              b += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i];
+              u += dfcdx[i*n + j] * y[i];
+              v += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i];
          }
-         sigma2 = sigma[j] * sigma[j];
-         dx = -a * sigma2 / b;
+         u *= (sigma2 = sqr(sigma[j]));
+         dx = u==0 ? 0 : (v/u)*sigma2 * (-1 + sqrt(1 - sqr(u/(v*sigma[j]))));
          xcur[j] = x[j] + dx;
          if (xcur[j] > ub[j]) xcur[j] = ub[j];
-         if (xcur[j] < lb[j]) xcur[j] = lb[j];
-         if (xcur[j] > x[j] + 0.9*sigma[j]) xcur[j] = x[j] + 0.9*sigma[j];
-         if (xcur[j] < x[j] - 0.9*sigma[j]) xcur[j] = x[j] - 0.9*sigma[j];
-         x_pinned = xcur[j] != x[j] + dx;
+         else if (xcur[j] < lb[j]) xcur[j] = lb[j];
+         if (xcur[j] > x[j]+0.9*sigma[j]) xcur[j] = x[j]+0.9*sigma[j];
+         else if (xcur[j] < x[j]-0.9*sigma[j]) xcur[j] = x[j]-0.9*sigma[j];
          dx = xcur[j] - x[j];
          
          /* function value: */
-         denom = sigma2 - dx*dx;
-         denominv = 1.0 / denom;
-         ca = sigma2 * dx;
-         cb = dx*dx;
-         val += (a*ca + b*cb) * denominv;
+         dx2 = dx * dx;
+         denominv = 1.0 / (sigma2 - dx2);
+         val += (u * dx + v * dx2) * denominv;
 
-         /* update gval, wval, gcval */
-         d->gval += (dfdx[j] * ca + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * cb)
+         /* update gval, wval, gcval (approximant functions) */
+         c = sigma2 * dx;
+         d->gval += (dfdx[j] * c + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * dx2)
               * denominv;
-         d->wval += 0.5 * cb * denominv;
+         d->wval += 0.5 * dx2 * denominv;
          for (i = 0; i < m; ++i)
-              gcval[i] += (dfcdx[i*n+j] * ca 
-                           + (fabs(dfcdx[i*n+j])*sigma[j] + 0.5*rhoc[j]) * cb)
+              gcval[i] += (dfcdx[i*n+j] * c + (fabs(dfcdx[i*n+j])*sigma[j] 
+                                               + 0.5*rhoc[j]) * dx2)
                    * denominv;
-
-         /* gradient */
-         if (grad)
-              for (i = 0; i < m; ++i) {
-                   double dady, dbdy;
-                   dady = dfcdx[i*n + j];
-                   dbdy = fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i];
-                   grad[i] += (dady*ca + dbdy*cb) * denominv;
-                   if (!x_pinned) { /* dx/dy nonzero */
-                        int k;
-                        double dxdy;
-                        dxdy = (a*dbdy - dady*b) * (sigma2 / (b*b));
-                        grad[i] += 
-                             ((sigma2 * dfdx[j] + 2 * (fabs(dfdx[j])*sigma[j]
-                                                       + 0.5*rho) * dx) 
-                              * denom
-                              + 2*dx * (dfdx[j]*ca + (fabs(dfdx[j])*sigma[j] 
-                                                      + 0.5*rho) * cb))
-                             * (denominv*denominv * dxdy);
-                        for (k = 0; k < m; ++k)
-                             grad[i] += 
-                                  ((sigma2 * dfcdx[k*n + j]
-                                    + 2 * (fabs(dfcdx[k*n + j])*sigma[j]
-                                           + 0.5*rhoc[k]) * dx) 
-                                   * denom
-                                   + 2*dx * (dfcdx[k*n + j]*ca 
-                                             + (fabs(dfcdx[k*n + j])*sigma[j] 
-                                                + 0.5*rhoc[k]) * cb))
-                                  * (denominv*denominv * dxdy) * y[k];
-                   }
-              }
      }
 
-     /* negate because we are maximizing */
-     if (grad) { for (i = 0; i < m; ++i) grad[i] = -grad[i]; }
+     /* gradient is easy to compute: since we are at a minimum x (dval/dx=0),
+       we only need the partial derivative with respect to y, and
+       we negate because we are maximizing: */
+     if (grad) for (i = 0; i < m; ++i) grad[i] = -gcval[i];
      return -val;
 }
 
@@ -127,7 +106,7 @@ static double dual_func(int m, const double *y, double *grad, void *d_)
 
 nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
                          int m, nlopt_func fc,
-                         void *fc_data_, ptrdiff_t fc_data_size,
+                         void *fc_data_, ptrdiff_t fc_datum_size,
                          const double *lb, const double *ub, /* bounds */
                          double *x, /* in: initial guess, out: minimizer */
                          double *minf,
@@ -188,13 +167,13 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
 
      dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
      stop->nevals++;
+     memcpy(xcur, x, sizeof(double) * n);
+
      feasible = 1;
      for (i = 0; i < m; ++i) {
-         fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_data_size * i);
+         fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_datum_size * i);
          feasible = feasible && (fcval[i] <= 0);
      }
-     memcpy(xcur, x, sizeof(double) * n);
-
      if (!feasible) { ret = NLOPT_FAILURE; goto done; } /* TODO: handle this */
 
      while (1) { /* outer iterations */
@@ -209,13 +188,21 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
          while (1) { /* inner iterations */
               double min_dual;
               int feasible_cur, inner_done;
+              nlopt_result reti;
 
               /* solve dual problem */
               dd.rho = rho;
-              nlopt_minimize(dual_alg, m, dual_func, &dd,
-                             dual_lb, dual_ub, y, &min_dual,
-                             -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval,
-                             stop->maxtime - (nlopt_seconds() - stop->start));
+              reti = nlopt_minimize(
+                   dual_alg, m, dual_func, &dd,
+                   dual_lb, dual_ub, y, &min_dual,
+                   -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval,
+                   stop->maxtime - (nlopt_seconds() - stop->start));
+              if (reti == NLOPT_FAILURE) reti = NLOPT_SUCCESS; /* hack */
+              if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
+                   ret = reti;
+                   goto done;
+              }
+
               dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 
               fcur = f(n, xcur, dfdx_cur, f_data);
@@ -224,7 +211,7 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
               feasible_cur = 1;
               for (i = 0; i < m; ++i) {
                    fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n, 
-                                     fc_data + fc_data_size * i);
+                                     fc_data + fc_datum_size * i);
                    feasible_cur = feasible_cur && (fcval_cur[i] <= 0);
                    inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]);
               }
@@ -255,6 +242,8 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
               
               if (mma_verbose)
                    printf("MMA inner iteration: rho -> %g\n", rho);
+              for (i = 0; i < mma_verbose; ++i)
+                   printf("                     rhoc[%d] -> %g\n", i,rhoc[i]);
          }
 
          if (nlopt_stop_ftol(stop, fcur, fprev))
@@ -269,14 +258,22 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
               printf("MMA outer iteration: rho -> %g\n", rho);
          for (i = 0; i < m; ++i)
               rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN);
-         if (k > 1)
+         for (i = 0; i < mma_verbose; ++i)
+              printf("                     rhoc[%d] -> %g\n", i, rhoc[i]);
+         if (k > 1) {
               for (j = 0; j < n; ++j) {
                    double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
                    double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
                    sigma[j] *= gam;
-                   sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
-                   sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
+                   if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
+                        sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
+                        sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
+                   }
               }
+              for (j = 0; j < mma_verbose; ++j)
+                   printf("                     sigma[%d] -> %g\n", 
+                          j, sigma[j]);
+         }
      }
 
  done:

diff --git a/mma/mma.h b/mma/mma.h
index 929a228f4dd8ff12421c201ca8e567ef872d6120..f6aac370a329b1a94b1573ed7edfc944a4fe4d0f 100644 (file)
--- a/mma/mma.h
+++ b/mma/mma.h
@@ -4,8 +4,6 @@
 #include "nlopt.h"
 #include "nlopt-util.h"
 
-#include <stddef.h> /* for ptrdiff_t */
-
 #ifdef __cplusplus
 extern "C"
 {
@@ -15,7 +13,7 @@ extern int mma_verbose;
 
 nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
                          int m, nlopt_func fc,
-                         void *fc_data_, ptrdiff_t fc_data_size,
+                         void *fc_data_, ptrdiff_t fc_datum_size,
                          const double *lb, const double *ub, /* bounds */
                          double *x, /* in: initial guess, out: minimizer */
                          double *minf,