initial (untested) stab at ~full MMA

darcs-hash:20080815003121-c8de0-638a9763e951bdab04146be55d5fc80125b8e5df.gz

diff --git a/api/nlopt.c b/api/nlopt.c
index 71570bdfe118cb6ac6d4f55e1cbcf290a394dfb7..fb58dd10648e624fa01254ceff436ab745971256 100644 (file)
--- a/api/nlopt.c
+++ b/api/nlopt.c
@@ -17,7 +17,7 @@
 #  define MY_INF HUGE_VAL
 #endif
 
-static int my_isinf(double x) {
+int nlopt_isinf(double x) {
      return fabs(x) >= HUGE_VAL * 0.99
 #ifdef HAVE_ISINF
          || isinf(x)
@@ -131,7 +131,7 @@ static double f_direct(int n, const double *x, int *undefined, void *data_)
      nlopt_data *data = (nlopt_data *) data_;
      double f;
      f = data->f(n, x, NULL, data->f_data);
-     *undefined = isnan(f) || my_isinf(f);
+     *undefined = isnan(f) || nlopt_isinf(f);
      return f;
 }
 
@@ -202,7 +202,12 @@ static nlopt_result nlopt_minimize_(
      nlopt_stopping stop;
 
      /* some basic argument checks */
-     if (n <= 0 || !f || !lb || !ub || !x || !minf)
+     if (!minf || !f) return NLOPT_INVALID_ARGS;
+     if (n == 0) { /* trivial case: no degrees of freedom */
+         *minf = f(n, x, NULL, f_data);
+         return NLOPT_SUCCESS;
+     }
+     else if (n < 0 || !lb || !ub || !x)
          return NLOPT_INVALID_ARGS;
 
      d.f = f;
@@ -220,7 +225,8 @@ static nlopt_result nlopt_minimize_(
               return NLOPT_INVALID_ARGS;
 
      stop.n = n;
-     stop.minf_max = (isnan(minf_max) || (my_isinf(minf_max) && minf_max < 0))
+     stop.minf_max = (isnan(minf_max) 
+                     || (nlopt_isinf(minf_max) && minf_max < 0))
          ? -MY_INF : minf_max;
      stop.ftol_rel = ftol_rel;
      stop.ftol_abs = ftol_abs;
@@ -297,11 +303,11 @@ static nlopt_result nlopt_minimize_(
              double *scale = (double *) malloc(sizeof(double) * n);
              if (!scale) return NLOPT_OUT_OF_MEMORY;
              for (i = 0; i < n; ++i) {
-                  if (!my_isinf(ub[i]) && !my_isinf(lb[i]))
+                  if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i]))
                        scale[i] = (ub[i] - lb[i]) * 0.01;
-                  else if (!my_isinf(lb[i]) && x[i] > lb[i])
+                  else if (!nlopt_isinf(lb[i]) && x[i] > lb[i])
                        scale[i] = (x[i] - lb[i]) * 0.01;
-                  else if (!my_isinf(ub[i]) && x[i] < ub[i])
+                  else if (!nlopt_isinf(ub[i]) && x[i] < ub[i])
                        scale[i] = (ub[i] - x[i]) * 0.01;
                   else
                        scale[i] = 0.01 * x[i] + 0.0001;
@@ -325,11 +331,11 @@ static nlopt_result nlopt_minimize_(
         case NLOPT_LN_PRAXIS: {
              double h0 = HUGE_VAL;
              for (i = 0; i < n; ++i) {
-                  if (!my_isinf(ub[i]) && !my_isinf(lb[i]))
+                  if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i]))
                        h0 = MIN(h0, (ub[i] - lb[i]) * 0.01);
-                  else if (!my_isinf(lb[i]) && x[i] > lb[i])
+                  else if (!nlopt_isinf(lb[i]) && x[i] > lb[i])
                        h0 = MIN(h0, (x[i] - lb[i]) * 0.01);
-                  else if (!my_isinf(ub[i]) && x[i] < ub[i])
+                  else if (!nlopt_isinf(ub[i]) && x[i] < ub[i])
                        h0 = MIN(h0, (ub[i] - x[i]) * 0.01);
                   else
                        h0 = MIN(h0, 0.01 * x[i] + 0.0001);
@@ -343,8 +349,8 @@ static nlopt_result nlopt_minimize_(
              int iret, *nbd = (int *) malloc(sizeof(int) * n);
              if (!nbd) return NLOPT_OUT_OF_MEMORY;
              for (i = 0; i < n; ++i) {
-                  int linf = my_isinf(lb[i]) && lb[i] < 0;
-                  int uinf = my_isinf(ub[i]) && ub[i] > 0;
+                  int linf = nlopt_isinf(lb[i]) && lb[i] < 0;
+                  int uinf = nlopt_isinf(ub[i]) && ub[i] > 0;
                   nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
              }
              iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub,

diff --git a/mma/mma.c b/mma/mma.c
index 891d83dd1df7a77e79a56e55892dfb35dd8a7465..280ee3f363141fe30f540de60bb678c881b18500 100644 (file)
--- a/mma/mma.c
+++ b/mma/mma.c
@@ -13,33 +13,188 @@ int mma_verbose = 0; /* > 0 for verbose output */
 /* magic minimum value for rho in MMA ... the 2002 paper says it should
    be a "fixed, strictly positive `small' number, e.g. 1e-5"
    ... grrr, I hate these magic numbers, which seem like they
-   should depend on the objective function in some way */
+   should depend on the objective function in some way ... in particular,
+   note that rho is dimensionful (= dimensions of objective function) */
 #define MMA_RHOMIN 1e-5
 
+/***********************************************************************/
+/* function for MMA's dual optimization of the approximant function */
+
+typedef struct {
+     int n; /* must be set on input to dimension of x */
+     const double *x, *lb, *ub, *sigma, *dfdx; /* arrays of length n */
+     const double *dfcdx; /* m-by-n array of fc gradients */
+     double fval, rho; /* must be set on input */
+     const double *fcval, *rhoc; /* arrays of length m */
+     double *xcur; /* array of length n, output each time */
+     double gval, wval, *gcval; /* output each time (array length m) */
+} dual_data;
+
+static double dual_func(int m, const double *y, double *grad, void *d_)
+{
+     dual_data *d = (dual_data *) d_;
+     int n = d->n;
+     const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma, 
+         *dfdx = d->dfdx;
+     const double *dfcdx = d->dfcdx;
+     double rho = d->rho, fval = d->fval;
+     const double *rhoc = d->rhoc, *fcval = d->fcval;
+     double *xcur = d->xcur;
+     double *gcval = d->gcval;
+     int i, j;
+     double val;
+
+     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;
+
+         /* first, compute xcur[j] for y */
+         a = dfdx[j];
+         b = 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];
+         }
+         sigma2 = sigma[j] * sigma[j];
+         dx = -a * sigma2 / b;
+         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;
+         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;
+
+         /* update gval, wval, gcval */
+         d->gval += (dfdx[j] * ca + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * cb)
+              * denominv;
+         d->wval += 0.5 * cb * 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)
+                   * 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]; }
+     return -val;
+}
+
+/***********************************************************************/
+
 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,
                          const double *lb, const double *ub, /* bounds */
                          double *x, /* in: initial guess, out: minimizer */
                          double *minf,
-                         nlopt_stopping *stop)
+                         nlopt_stopping *stop,
+                         nlopt_algorithm dual_alg, 
+                         double dual_tolrel, int dual_maxeval)
 {
      nlopt_result ret = NLOPT_SUCCESS;
      double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
-     int j, k = 0;
+     double *dfcdx, *dfcdx_cur, *fcval, *rhoc, *gcval, *y, *dual_lb, *dual_ub;
+     int i, j, k = 0;
+     char *fc_data = (char *) fc_data_;
+     dual_data dd;
+     int feasible;
      
-     sigma = (double *) malloc(sizeof(double) * 6*n);
+     sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*6));
      if (!sigma) return NLOPT_OUT_OF_MEMORY;
      dfdx = sigma + n;
      dfdx_cur = dfdx + n;
      xcur = dfdx_cur + n;
      xprev = xcur + n;
      xprevprev = xprev + n;
-     for (j = 0; j < n; ++j)
-         sigma[j] = 0.5 * (ub[j] - lb[j]);
+     fcval = xprevprev + n;
+     rhoc = fcval + m;
+     gcval = rhoc + m;
+     dual_lb = gcval + m;
+     dual_ub = dual_lb + m;
+     y = dual_ub + m;
+     dfcdx = y + m;
+     dfcdx_cur = dfcdx + m*n;
+
+     dd.n = n;
+     dd.x = x;
+     dd.lb = lb;
+     dd.ub = ub;
+     dd.sigma = sigma;
+     dd.dfdx = dfdx;
+     dd.dfcdx = dfcdx;
+     dd.fcval = fcval;
+     dd.rhoc = rhoc;
+     dd.xcur = xcur;
+     dd.gcval = gcval;
+
+     for (j = 0; j < n; ++j) {
+         if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
+              sigma[j] = 1.0; /* arbitrary default */
+         else
+              sigma[j] = 0.5 * (ub[j] - lb[j]);
+     }
      rho = 1.0;
+     for (i = 0; i < m; ++i) {
+         rhoc[i] = 1.0;
+         dual_lb[i] = y[i] = 0.0;
+         dual_ub[i] = HUGE_VAL;
+     }
 
-     fcur = *minf = f(n, x, dfdx, f_data);
-     memcpy(xcur, x, sizeof(double) * n);
+     dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
      stop->nevals++;
+     feasible = 1;
+     for (i = 0; i < m; ++i) {
+         fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_data_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 */
          double fprev = fcur;
          if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
@@ -50,40 +205,52 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
          memcpy(xprev, xcur, sizeof(double) * n);
 
          while (1) { /* inner iterations */
-              double gcur = *minf, w = 0;
-              for (j = 0; j < n; ++j) {
-                   /* because we have no constraint functions, minimizing
-                      the MMA approximate function can be done analytically */
-                   double dx = -sigma[j]*sigma[j]*dfdx[j]
-                        / (2*sigma[j]*fabs(dfdx[j]) + rho);
-                   xcur[j] = x[j] + dx;
-                   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];
-                   if (xcur[j] > ub[j]) xcur[j] = ub[j];
-                   else if (xcur[j] < lb[j]) xcur[j] = lb[j];
-                   dx = xcur[j] - x[j];
-                   gcur += (sigma[j]*sigma[j]*dfdx[j]*dx
-                            + sigma[j]*fabs(dfdx[j])*dx*dx
-                            + 0.5*rho*dx*dx) / (sigma[j]*sigma[j]-dx*dx);
-                   w += 0.5*rho*dx*dx / (sigma[j]*sigma[j]-dx*dx);
-              }
+              double min_dual;
+              int feasible_cur, inner_done;
+
+              /* 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));
+              dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 
               fcur = f(n, xcur, dfdx_cur, f_data);
               stop->nevals++;
-              if (fcur < *minf) {
-                   *minf = fcur;
+              inner_done = dd.gval >= fcur;
+              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);
+                   feasible_cur = feasible_cur && (fcval_cur[i] <= 0);
+                   inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]);
+              }
+
+              if (fcur < *minf && (feasible_cur || !feasible)) {
+                   feasible = feasible_cur;
+                   dd.fval = *minf = fcur;
+                   memcpy(fcval, fcval_cur, sizeof(double)*m);
                    memcpy(x, xcur, sizeof(double)*n);
                    memcpy(dfdx, dfdx_cur, sizeof(double)*n);
+                   memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
               }
               if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
               else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
               else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
               if (ret != NLOPT_SUCCESS) goto done;
 
-              if (gcur >= fcur) break;
-              rho = MIN(10*rho, 1.1 * (rho + (fcur - gcur) / w));
+              if (inner_done) break;
+
+              if (fcur > dd.gval)
+                   rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
+              for (i = 0; i < m; ++i)
+                   if (fcval_cur[i] > dd.gcval[i])
+                        rhoc[i] = 
+                             MIN(10*rhoc[i], 
+                                 1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i]) 
+                                        / dd.wval));
+              
               if (mma_verbose)
                    printf("MMA inner iteration: rho -> %g\n", rho);
          }
@@ -98,6 +265,8 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
          rho = MAX(0.1 * rho, MMA_RHOMIN);
          if (mma_verbose)
               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 (j = 0; j < n; ++j) {
                    double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);

diff --git a/util/nlopt-util.h b/util/nlopt-util.h
index e0eb13485c2c70a7f438015f9817f7c6b5e1c332..fffb750b2ce419159a38b36aaeeea5a18eb5372b 100644 (file)
--- a/util/nlopt-util.h
+++ b/util/nlopt-util.h
@@ -6,6 +6,8 @@ extern "C"
 {
 #endif /* __cplusplus */
 
+int nlopt_isinf(double x);
+
 /* seconds timer */
 extern double nlopt_seconds(void);
 extern unsigned long nlopt_time_seed(void);