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 ... in particular,
  17    note that rho is dimensionful (= dimensions of objective function) */
  18 #define MMA_RHOMIN 1e-5
  19
  20 /***********************************************************************/
  21 /* function for MMA's dual optimization of the approximant function */
  22
  23 typedef struct {
  24      int n; /* must be set on input to dimension of x */
  25      const double *x, *lb, *ub, *sigma, *dfdx; /* arrays of length n */
  26      const double *dfcdx; /* m-by-n array of fc gradients */
  27      double fval, rho; /* must be set on input */
  28      const double *fcval, *rhoc; /* arrays of length m */
  29      double *xcur; /* array of length n, output each time */
  30      double gval, wval, *gcval; /* output each time (array length m) */
  31 } dual_data;
  32
  33 static double dual_func(int m, const double *y, double *grad, void *d_)
  34 {
  35      dual_data *d = (dual_data *) d_;
  36      int n = d->n;
  37      const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma,
  38           *dfdx = d->dfdx;
  39      const double *dfcdx = d->dfcdx;
  40      double rho = d->rho, fval = d->fval;
  41      const double *rhoc = d->rhoc, *fcval = d->fcval;
  42      double *xcur = d->xcur;
  43      double *gcval = d->gcval;
  44      int i, j;
  45      double val;
  46
  47      val = d->gval = fval;
  48      d->wval = 0;
  49      for (i = 0; i < m; ++i) val += y[i] * (gcval[i] = fcval[i]);
  50      if (grad) { for (i = 0; i < m; ++i) grad[i] = fcval[i]; }
  51
  52      for (j = 0; j < n; ++j) {
  53           int x_pinned;
  54           double a, b, dx, denom, denominv, ca, cb, sigma2;
  55
  56           /* first, compute xcur[j] for y */
  57           a = dfdx[j];
  58           b = fabs(dfdx[j]) * sigma[j] + 0.5 * rho;
  59           for (i = 0; i < m; ++i) {
  60                a += dfcdx[i*n + j] * y[i];
  61                b += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i];
  62           }
  63           sigma2 = sigma[j] * sigma[j];
  64           dx = -a * sigma2 / b;
  65           xcur[j] = x[j] + dx;
  66           if (xcur[j] > ub[j]) xcur[j] = ub[j];
  67           if (xcur[j] < lb[j]) xcur[j] = lb[j];
  68           if (xcur[j] > x[j] + 0.9*sigma[j]) xcur[j] = x[j] + 0.9*sigma[j];
  69           if (xcur[j] < x[j] - 0.9*sigma[j]) xcur[j] = x[j] - 0.9*sigma[j];
  70           x_pinned = xcur[j] != x[j] + dx;
  71           dx = xcur[j] - x[j];
  72
  73           /* function value: */
  74           denom = sigma2 - dx*dx;
  75           denominv = 1.0 / denom;
  76           ca = sigma2 * dx;
  77           cb = dx*dx;
  78           val += (a*ca + b*cb) * denominv;
  79
  80           /* update gval, wval, gcval */
  81           d->gval += (dfdx[j] * ca + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * cb)
  82                * denominv;
  83           d->wval += 0.5 * cb * denominv;
  84           for (i = 0; i < m; ++i)
  85                gcval[i] += (dfcdx[i*n+j] * ca
  86                             + (fabs(dfcdx[i*n+j])*sigma[j] + 0.5*rhoc[j]) * cb)
  87                     * denominv;
  88
  89           /* gradient */
  90           if (grad)
  91                for (i = 0; i < m; ++i) {
  92                     double dady, dbdy;
  93                     dady = dfcdx[i*n + j];
  94                     dbdy = fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i];
  95                     grad[i] += (dady*ca + dbdy*cb) * denominv;
  96                     if (!x_pinned) { /* dx/dy nonzero */
  97                          int k;
  98                          double dxdy;
  99                          dxdy = (a*dbdy - dady*b) * (sigma2 / (b*b));
 100                          grad[i] +=
 101                               ((sigma2 * dfdx[j] + 2 * (fabs(dfdx[j])*sigma[j]
 102                                                         + 0.5*rho) * dx)
 103                                * denom
 104                                + 2*dx * (dfdx[j]*ca + (fabs(dfdx[j])*sigma[j]
 105                                                        + 0.5*rho) * cb))
 106                               * (denominv*denominv * dxdy);
 107                          for (k = 0; k < m; ++k)
 108                               grad[i] +=
 109                                    ((sigma2 * dfcdx[k*n + j]
 110                                      + 2 * (fabs(dfcdx[k*n + j])*sigma[j]
 111                                             + 0.5*rhoc[k]) * dx)
 112                                     * denom
 113                                     + 2*dx * (dfcdx[k*n + j]*ca
 114                                               + (fabs(dfcdx[k*n + j])*sigma[j]
 115                                                  + 0.5*rhoc[k]) * cb))
 116                                    * (denominv*denominv * dxdy) * y[k];
 117                     }
 118                }
 119      }
 120
 121      /* negate because we are maximizing */
 122      if (grad) { for (i = 0; i < m; ++i) grad[i] = -grad[i]; }
 123      return -val;
 124 }
 125
 126 /***********************************************************************/
 127
 128 nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
 129                           int m, nlopt_func fc,
 130                           void *fc_data_, ptrdiff_t fc_data_size,
 131                           const double *lb, const double *ub, /* bounds */
 132                           double *x, /* in: initial guess, out: minimizer */
 133                           double *minf,
 134                           nlopt_stopping *stop,
 135                           nlopt_algorithm dual_alg,
 136                           double dual_tolrel, int dual_maxeval)
 137 {
 138      nlopt_result ret = NLOPT_SUCCESS;
 139      double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
 140      double *dfcdx, *dfcdx_cur, *fcval, *rhoc, *gcval, *y, *dual_lb, *dual_ub;
 141      int i, j, k = 0;
 142      char *fc_data = (char *) fc_data_;
 143      dual_data dd;
 144      int feasible;
 145
 146      sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*6));
 147      if (!sigma) return NLOPT_OUT_OF_MEMORY;
 148      dfdx = sigma + n;
 149      dfdx_cur = dfdx + n;
 150      xcur = dfdx_cur + n;
 151      xprev = xcur + n;
 152      xprevprev = xprev + n;
 153      fcval = xprevprev + n;
 154      rhoc = fcval + m;
 155      gcval = rhoc + m;
 156      dual_lb = gcval + m;
 157      dual_ub = dual_lb + m;
 158      y = dual_ub + m;
 159      dfcdx = y + m;
 160      dfcdx_cur = dfcdx + m*n;
 161
 162      dd.n = n;
 163      dd.x = x;
 164      dd.lb = lb;
 165      dd.ub = ub;
 166      dd.sigma = sigma;
 167      dd.dfdx = dfdx;
 168      dd.dfcdx = dfcdx;
 169      dd.fcval = fcval;
 170      dd.rhoc = rhoc;
 171      dd.xcur = xcur;
 172      dd.gcval = gcval;
 173
 174      for (j = 0; j < n; ++j) {
 175           if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
 176                sigma[j] = 1.0; /* arbitrary default */
 177           else
 178                sigma[j] = 0.5 * (ub[j] - lb[j]);
 179      }
 180      rho = 1.0;
 181      for (i = 0; i < m; ++i) {
 182           rhoc[i] = 1.0;
 183           dual_lb[i] = y[i] = 0.0;
 184           dual_ub[i] = HUGE_VAL;
 185      }
 186
 187      dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
 188      stop->nevals++;
 189      feasible = 1;
 190      for (i = 0; i < m; ++i) {
 191           fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_data_size * i);
 192           feasible = feasible && (fcval[i] <= 0);
 193      }
 194      memcpy(xcur, x, sizeof(double) * n);
 195
 196      if (!feasible) { ret = NLOPT_FAILURE; goto done; } /* TODO: handle this */
 197
 198      while (1) { /* outer iterations */
 199           double fprev = fcur;
 200           if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 201           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 202           else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 203           if (ret != NLOPT_SUCCESS) goto done;
 204           if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
 205           memcpy(xprev, xcur, sizeof(double) * n);
 206
 207           while (1) { /* inner iterations */
 208                double min_dual;
 209                int feasible_cur, inner_done;
 210
 211                /* solve dual problem */
 212                dd.rho = rho;
 213                nlopt_minimize(dual_alg, m, dual_func, &dd,
 214                               dual_lb, dual_ub, y, &min_dual,
 215                               -HUGE_VAL, dual_tolrel, 0, 0, NULL, dual_maxeval,
 216                               stop->maxtime - (nlopt_seconds() - stop->start));
 217                dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 218
 219                fcur = f(n, xcur, dfdx_cur, f_data);
 220                stop->nevals++;
 221                inner_done = dd.gval >= fcur;
 222                feasible_cur = 1;
 223                for (i = 0; i < m; ++i) {
 224                     fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n,
 225                                       fc_data + fc_data_size * i);
 226                     feasible_cur = feasible_cur && (fcval_cur[i] <= 0);
 227                     inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]);
 228                }
 229
 230                if (fcur < *minf && (feasible_cur || !feasible)) {
 231                     feasible = feasible_cur;
 232                     dd.fval = *minf = fcur;
 233                     memcpy(fcval, fcval_cur, sizeof(double)*m);
 234                     memcpy(x, xcur, sizeof(double)*n);
 235                     memcpy(dfdx, dfdx_cur, sizeof(double)*n);
 236                     memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
 237                }
 238                if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 239                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 240                else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 241                if (ret != NLOPT_SUCCESS) goto done;
 242
 243                if (inner_done) break;
 244
 245                if (fcur > dd.gval)
 246                     rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
 247                for (i = 0; i < m; ++i)
 248                     if (fcval_cur[i] > dd.gcval[i])
 249                          rhoc[i] =
 250                               MIN(10*rhoc[i],
 251                                   1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
 252                                          / dd.wval));
 253
 254                if (mma_verbose)
 255                     printf("MMA inner iteration: rho -> %g\n", rho);
 256           }
 257
 258           if (nlopt_stop_ftol(stop, fcur, fprev))
 259                ret = NLOPT_FTOL_REACHED;
 260           if (nlopt_stop_x(stop, xcur, xprev))
 261                ret = NLOPT_XTOL_REACHED;
 262           if (ret != NLOPT_SUCCESS) goto done;
 263
 264           /* update rho and sigma for iteration k+1 */
 265           rho = MAX(0.1 * rho, MMA_RHOMIN);
 266           if (mma_verbose)
 267                printf("MMA outer iteration: rho -> %g\n", rho);
 268           for (i = 0; i < m; ++i)
 269                rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN);
 270           if (k > 1)
 271                for (j = 0; j < n; ++j) {
 272                     double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
 273                     double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
 274                     sigma[j] *= gam;
 275                     sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
 276                     sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
 277                }
 278      }
 279
 280  done:
 281      free(sigma);
 282      return ret;
 283 }