mma/mma.c

   1 #include <stdlib.h>
   2 #include <math.h>
   3 #include <string.h>
   4 #include <stdio.h>
   5
   6 #include "mma.h"
   7
   8 int mma_verbose = 0; /* > 0 for verbose output */
   9
  10 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  11 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  12
  13 /* magic minimum value for rho in MMA ... the 2002 paper says it should
  14    be a "fixed, strictly positive `small' number, e.g. 1e-5"
  15    ... grrr, I hate these magic numbers, which seem like they
  16    should depend on the objective function in some way */
  17 #define MMA_RHOMIN 1e-5
  18
  19 nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
  20                           const double *lb, const double *ub, /* bounds */
  21                           double *x, /* in: initial guess, out: minimizer */
  22                           double *minf,
  23                           nlopt_stopping *stop)
  24 {
  25      nlopt_result ret = NLOPT_SUCCESS;
  26      double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
  27      int j, k = 0;
  28
  29      sigma = (double *) malloc(sizeof(double) * 6*n);
  30      if (!sigma) return NLOPT_OUT_OF_MEMORY;
  31      dfdx = sigma + n;
  32      dfdx_cur = dfdx + n;
  33      xcur = dfdx_cur + n;
  34      xprev = xcur + n;
  35      xprevprev = xprev + n;
  36      for (j = 0; j < n; ++j)
  37           sigma[j] = 0.5 * (ub[j] - lb[j]);
  38      rho = 1.0;
  39
  40      fcur = *minf = f(n, x, dfdx, f_data);
  41      memcpy(xcur, x, sizeof(double) * n);
  42      stop->nevals++;
  43      while (1) { /* outer iterations */
  44           double fprev = fcur;
  45           if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
  46           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
  47           else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
  48           if (ret != NLOPT_SUCCESS) goto done;
  49           if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
  50           memcpy(xprev, xcur, sizeof(double) * n);
  51
  52           while (1) { /* inner iterations */
  53                double gcur = *minf, w = 0;
  54                for (j = 0; j < n; ++j) {
  55                     /* because we have no constraint functions, minimizing
  56                        the MMA approximate function can be done analytically */
  57                     double dx = -sigma[j]*sigma[j]*dfdx[j]
  58                          / (2*sigma[j]*fabs(dfdx[j]) + rho);
  59                     xcur[j] = x[j] + dx;
  60                     if (xcur[j] > x[j] + 0.9*sigma[j])
  61                          xcur[j] = x[j] + 0.9*sigma[j];
  62                     else if (xcur[j] < x[j] - 0.9*sigma[j])
  63                          xcur[j] = x[j] - 0.9*sigma[j];
  64                     if (xcur[j] > ub[j]) xcur[j] = ub[j];
  65                     else if (xcur[j] < lb[j]) xcur[j] = lb[j];
  66                     dx = xcur[j] - x[j];
  67                     gcur += (sigma[j]*sigma[j]*dfdx[j]*dx
  68                              + sigma[j]*fabs(dfdx[j])*dx*dx
  69                              + 0.5*rho*dx*dx) / (sigma[j]*sigma[j]-dx*dx);
  70                     w += 0.5*rho*dx*dx / (sigma[j]*sigma[j]-dx*dx);
  71                }
  72
  73                fcur = f(n, xcur, dfdx_cur, f_data);
  74                stop->nevals++;
  75                if (fcur < *minf) {
  76                     *minf = fcur;
  77                     memcpy(x, xcur, sizeof(double)*n);
  78                     memcpy(dfdx, dfdx_cur, sizeof(double)*n);
  79                }
  80                if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
  81                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
  82                else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
  83                if (ret != NLOPT_SUCCESS) goto done;
  84
  85                if (gcur >= fcur) break;
  86                rho = MIN(10*rho, 1.1 * (rho + (fcur - gcur) / w));
  87                if (mma_verbose)
  88                     printf("MMA inner iteration: rho -> %g\n", rho);
  89           }
  90
  91           if (nlopt_stop_ftol(stop, fcur, fprev))
  92                ret = NLOPT_FTOL_REACHED;
  93           if (nlopt_stop_x(stop, xcur, xprev))
  94                ret = NLOPT_XTOL_REACHED;
  95           if (ret != NLOPT_SUCCESS) goto done;
  96
  97           /* update rho and sigma for iteration k+1 */
  98           rho = MAX(0.1 * rho, MMA_RHOMIN);
  99           if (mma_verbose)
 100                printf("MMA outer iteration: rho -> %g\n", rho);
 101           if (k > 1)
 102                for (j = 0; j < n; ++j) {
 103                     double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
 104                     double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
 105                     sigma[j] *= gam;
 106                     sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
 107                     sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
 108                }
 109      }
 110
 111  done:
 112      free(sigma);
 113      return ret;
 114 }