api/nlopt.c

   1 #include <stdlib.h>
   2 #include <math.h>
   3 #include <string.h>
   4
   5 #include "nlopt.h"
   6 #include "nlopt-util.h"
   7 #include "config.h"
   8
   9 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  10
  11 /*************************************************************************/
  12
  13 #ifdef INFINITY
  14 #  define MY_INF INFINITY
  15 #else
  16 #  define MY_INF HUGE_VAL
  17 #endif
  18
  19 static int my_isinf(double x) {
  20      return fabs(x) >= HUGE_VAL * 0.99
  21 #ifdef HAVE_ISINF
  22           || isinf(x)
  23 #endif
  24           ;
  25 }
  26
  27 #ifndef HAVE_ISNAN
  28 static int my_isnan(double x) { return x != x; }
  29 #  define isnan my_isnan
  30 #endif
  31
  32 /*************************************************************************/
  33
  34 void nlopt_version(int *major, int *minor, int *bugfix)
  35 {
  36      *major = MAJOR_VERSION;
  37      *minor = MINOR_VERSION;
  38      *bugfix = BUGFIX_VERSION;
  39 }
  40
  41 /*************************************************************************/
  42
  43 static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = {
  44      "DIRECT (global)",
  45      "DIRECT-L (global)",
  46      "Subplex (local)",
  47      "StoGO (global)",
  48      "StoGO with randomized search (global)",
  49      "Low-storage BFGS (LBFGS) (local)"
  50 };
  51
  52 const char *nlopt_algorithm_name(nlopt_algorithm a)
  53 {
  54      if (a < 0 || a >= NLOPT_NUM_ALGORITHMS) return "UNKNOWN";
  55      return nlopt_algorithm_names[a];
  56 }
  57
  58 /*************************************************************************/
  59
  60 static int nlopt_srand_called = 0;
  61 void nlopt_srand(unsigned long seed) {
  62      nlopt_srand_called = 1;
  63      nlopt_init_genrand(seed);
  64 }
  65
  66 void nlopt_srand_time() {
  67      nlopt_srand(nlopt_time_seed());
  68 }
  69
  70 /*************************************************************************/
  71
  72 typedef struct {
  73      nlopt_func f;
  74      void *f_data;
  75      const double *lb, *ub, *x0;
  76      double *xtmp;
  77 } nlopt_data;
  78
  79 #define RECENTER 1 /* 0 to disable recentering */
  80
  81 /* for global-search algorithms that ignore the starting guess,
  82    but always check the center of the search box, we perform a
  83    coordinate transformation to put the initial guess x0 at the
  84    center of the box, and store the transformed x in xtmp. */
  85 static void recenter_x(int n, const double *x,
  86                        const double *lb, const double *ub,
  87                        const double *x0, double *xtmp)
  88 {
  89      int i;
  90      for (i = 0; i < n; ++i) {
  91 #if RECENTER
  92           /* Lagrange interpolating polynomial */
  93           double xm = 0.5 * (lb[i] + ub[i]);
  94           double dlu = 1. / (lb[i] - ub[i]);
  95           double dlm = 1. / (lb[i] - xm);
  96           double dum = 1. / (ub[i] - xm);
  97           double dxu = x[i] - ub[i];
  98           double dxl = x[i] - lb[i];
  99           double dxm = x[i] - xm;
 100           xtmp[i] = (lb[i] * (dxu * dlu) * (dxm * dlm)
 101                      - ub[i] * (dxl * dlu) * (dxm * dum)
 102                      + x0[i] * (dxl * dlm) * (dxu * dum));
 103 #else
 104           xtmp[i] = x[i];
 105 #endif
 106      }
 107 }
 108
 109 /* transform grad from df/dxtmp to df/dx */
 110 static void recenter_grad(int n, const double *x,
 111                           const double *lb, const double *ub,
 112                           const double *x0,
 113                           double *grad)
 114 {
 115 #if RECENTER
 116      if (grad) {
 117           int i;
 118           for (i = 0; i < n; ++i) {
 119                double xm = 0.5 * (lb[i] + ub[i]);
 120                double dlu = 1. / (lb[i] - ub[i]);
 121                double dlm = 1. / (lb[i] - xm);
 122                double dum = 1. / (ub[i] - xm);
 123                double dxu = x[i] - ub[i];
 124                double dxl = x[i] - lb[i];
 125                double dxm = x[i] - xm;
 126                grad[i] *= (lb[i] * dlu*dlm * (dxm + dxu)
 127                            - ub[i] * dum*dlu * (dxm + dxl)
 128                            + x0[i] * dlm*dum * (dxu + dxl));
 129           }
 130      }
 131 #endif
 132 }
 133
 134 static double f_recenter(int n, const double *x, double *grad, void *data_)
 135 {
 136      nlopt_data *data = (nlopt_data *) data_;
 137      double f;
 138      recenter_x(n, x, data->lb, data->ub, data->x0, data->xtmp);
 139      f = data->f(n, data->xtmp, grad, data->f_data);
 140      recenter_grad(n, x, data->lb, data->ub, data->x0, grad);
 141      return f;
 142 }
 143
 144 #include "subplex.h"
 145
 146 static double f_subplex(int n, const double *x, void *data_)
 147 {
 148      int i;
 149      nlopt_data *data = (nlopt_data *) data_;
 150      double f;
 151
 152      /* subplex does not support bound constraints, but it supports
 153         discontinuous objectives so we can just return Inf for invalid x */
 154      for (i = 0; i < n; ++i)
 155           if (x[i] < data->lb[i] || x[i] > data->ub[i])
 156                return MY_INF;
 157
 158      f = data->f(n, x, NULL, data->f_data);
 159      return (isnan(f) ? MY_INF : f);
 160 }
 161
 162 #include "direct.h"
 163
 164 static double f_direct(int n, const double *x, int *undefined, void *data_)
 165 {
 166      nlopt_data *data = (nlopt_data *) data_;
 167      double f;
 168      recenter_x(n, x, data->lb, data->ub, data->x0, data->xtmp);
 169      f = data->f(n, data->xtmp, NULL, data->f_data);
 170      *undefined = isnan(f) || my_isinf(f);
 171      return f;
 172 }
 173
 174 #include "stogo.h"
 175
 176 #include "l-bfgs-b.h"
 177
 178 /*************************************************************************/
 179
 180 /* same as nlopt_minimize, but xtol_abs is required to be non-NULL */
 181 static nlopt_result nlopt_minimize_(
 182      nlopt_algorithm algorithm,
 183      int n, nlopt_func f, void *f_data,
 184      const double *lb, const double *ub, /* bounds */
 185      double *x, /* in: initial guess, out: minimizer */
 186      double *fmin, /* out: minimum */
 187      double fmin_max, double ftol_rel, double ftol_abs,
 188      double xtol_rel, const double *xtol_abs,
 189      int maxeval, double maxtime)
 190 {
 191      int i;
 192      nlopt_data d;
 193      nlopt_stopping stop;
 194
 195      d.f = f;
 196      d.f_data = f_data;
 197      d.lb = lb;
 198      d.ub = ub;
 199      d.x0 = d.xtmp = NULL;
 200
 201      /* make sure rand generator is inited */
 202      if (!nlopt_srand_called)
 203           nlopt_srand_time(); /* default is non-deterministic */
 204
 205      /* check bound constraints */
 206      for (i = 0; i < n; ++i)
 207           if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
 208                return NLOPT_INVALID_ARGS;
 209
 210      stop.n = n;
 211      stop.fmin_max = (isnan(fmin_max) || (my_isinf(fmin_max) && fmin_max < 0))
 212           ? -MY_INF : fmin_max;
 213      stop.ftol_rel = ftol_rel;
 214      stop.ftol_abs = ftol_abs;
 215      stop.xtol_rel = xtol_rel;
 216      stop.xtol_abs = xtol_abs;
 217      stop.nevals = 0;
 218      stop.maxeval = maxeval;
 219      stop.maxtime = maxtime;
 220      stop.start = nlopt_seconds();
 221
 222      switch (algorithm) {
 223          case NLOPT_GLOBAL_DIRECT:
 224          case NLOPT_GLOBAL_DIRECT_L: {
 225               int iret;
 226               d.xtmp = (double *) malloc(sizeof(double) * n*2);
 227               if (!d.xtmp) return NLOPT_OUT_OF_MEMORY;
 228               memcpy(d.xtmp + n, x, sizeof(double) * n); d.x0 = d.xtmp + n;
 229               iret = direct_optimize(f_direct, &d, n, lb, ub, x, fmin,
 230                                      maxeval, 500, ftol_rel, ftol_abs,
 231                                      xtol_rel, xtol_rel,
 232                                      DIRECT_UNKNOWN_FGLOBAL, -1.0,
 233                                      NULL,
 234                                      algorithm == NLOPT_GLOBAL_DIRECT
 235                                      ? DIRECT_ORIGINAL
 236                                      : DIRECT_GABLONSKY);
 237               recenter_x(n, x, lb, ub, d.x0, x);
 238               free(d.xtmp);
 239               switch (iret) {
 240                   case DIRECT_INVALID_BOUNDS:
 241                   case DIRECT_MAXFEVAL_TOOBIG:
 242                   case DIRECT_INVALID_ARGS:
 243                        return NLOPT_INVALID_ARGS;
 244                   case DIRECT_INIT_FAILED:
 245                   case DIRECT_SAMPLEPOINTS_FAILED:
 246                   case DIRECT_SAMPLE_FAILED:
 247                        return NLOPT_FAILURE;
 248                   case DIRECT_MAXFEVAL_EXCEEDED:
 249                   case DIRECT_MAXITER_EXCEEDED:
 250                        return NLOPT_MAXEVAL_REACHED;
 251                   case DIRECT_GLOBAL_FOUND:
 252                        return NLOPT_SUCCESS;
 253                   case DIRECT_VOLTOL:
 254                   case DIRECT_SIGMATOL:
 255                        return NLOPT_XTOL_REACHED;
 256                   case DIRECT_OUT_OF_MEMORY:
 257                        return NLOPT_OUT_OF_MEMORY;
 258               }
 259               break;
 260          }
 261
 262          case NLOPT_GLOBAL_STOGO:
 263          case NLOPT_GLOBAL_STOGO_RANDOMIZED: {
 264               int iret;
 265               d.xtmp = (double *) malloc(sizeof(double) * n*2);
 266               if (!d.xtmp) return NLOPT_OUT_OF_MEMORY;
 267               memcpy(d.xtmp + n, x, sizeof(double) * n); d.x0 = d.xtmp + n;
 268               iret = stogo_minimize(n, f_recenter, &d, x, fmin, lb, ub, &stop,
 269                                     algorithm == NLOPT_GLOBAL_STOGO
 270                                     ? 0 : 2*n);
 271               recenter_x(n, x, lb, ub, d.x0, x);
 272               free(d.xtmp);
 273               if (!iret) return NLOPT_FAILURE;
 274               break;
 275          }
 276
 277          case NLOPT_LOCAL_SUBPLEX: {
 278               int iret;
 279               double *scale = (double *) malloc(sizeof(double) * n);
 280               if (!scale) return NLOPT_OUT_OF_MEMORY;
 281               for (i = 0; i < n; ++i) {
 282                    if (!my_isinf(ub[i]) && !my_isinf(lb[i]))
 283                         scale[i] = (ub[i] - lb[i]) * 0.01;
 284                    else if (!my_isinf(lb[i]) && x[i] > lb[i])
 285                         scale[i] = (x[i] - lb[i]) * 0.01;
 286                    else if (!my_isinf(ub[i]) && x[i] < ub[i])
 287                         scale[i] = (ub[i] - x[i]) * 0.01;
 288                    else
 289                         scale[i] = 0.01 * x[i] + 0.0001;
 290               }
 291               iret = subplex(f_subplex, fmin, x, n, &d, &stop, scale);
 292               free(scale);
 293               switch (iret) {
 294                   case -2: return NLOPT_INVALID_ARGS;
 295                   case -10: return NLOPT_MAXTIME_REACHED;
 296                   case -1: return NLOPT_MAXEVAL_REACHED;
 297                   case 0: return NLOPT_XTOL_REACHED;
 298                   case 1: return NLOPT_SUCCESS;
 299                   case 2: return NLOPT_FMIN_MAX_REACHED;
 300                   case 20: return NLOPT_FTOL_REACHED;
 301                   case -200: return NLOPT_OUT_OF_MEMORY;
 302                   default: return NLOPT_FAILURE; /* unknown return code */
 303               }
 304               break;
 305          }
 306
 307          case NLOPT_LOCAL_LBFGS: {
 308               int iret, *nbd = (int *) malloc(sizeof(int) * n);
 309               if (!nbd) return NLOPT_OUT_OF_MEMORY;
 310               for (i = 0; i < n; ++i) {
 311                    int linf = my_isinf(lb[i]) && lb[i] < 0;
 312                    int uinf = my_isinf(ub[i]) && ub[i] > 0;
 313                    nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
 314               }
 315               iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub,
 316                                      MIN(n, 5), 0.0, ftol_rel,
 317                                      xtol_abs ? *xtol_abs : xtol_rel,
 318                                      maxeval);
 319               free(nbd);
 320               if (iret <= 0) {
 321                    switch (iret) {
 322                        case -1: return NLOPT_INVALID_ARGS;
 323                        case -2: default: return NLOPT_FAILURE;
 324                    }
 325               }
 326               else {
 327                    *fmin = f(n, x, NULL, f_data);
 328                    switch (iret) {
 329                        case 5: return NLOPT_MAXEVAL_REACHED;
 330                        case 2: return NLOPT_XTOL_REACHED;
 331                        case 1: return NLOPT_FTOL_REACHED;
 332                        default: return NLOPT_SUCCESS;
 333                    }
 334               }
 335               break;
 336          }
 337
 338          default:
 339               return NLOPT_INVALID_ARGS;
 340      }
 341
 342      return NLOPT_SUCCESS;
 343 }
 344
 345 nlopt_result nlopt_minimize(
 346      nlopt_algorithm algorithm,
 347      int n, nlopt_func f, void *f_data,
 348      const double *lb, const double *ub, /* bounds */
 349      double *x, /* in: initial guess, out: minimizer */
 350      double *fmin, /* out: minimum */
 351      double fmin_max, double ftol_rel, double ftol_abs,
 352      double xtol_rel, const double *xtol_abs,
 353      int maxeval, double maxtime)
 354 {
 355      nlopt_result ret;
 356      if (xtol_abs)
 357           ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub,
 358                                 x, fmin, fmin_max, ftol_rel, ftol_abs,
 359                                 xtol_rel, xtol_abs, maxeval, maxtime);
 360      else {
 361           int i;
 362           double *xtol = (double *) malloc(sizeof(double) * n);
 363           if (!xtol) return NLOPT_OUT_OF_MEMORY;
 364           for (i = 0; i < n; ++i) xtol[i] = -1;
 365           ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub,
 366                                 x, fmin, fmin_max, ftol_rel, ftol_abs,
 367                                 xtol_rel, xtol, maxeval, maxtime);
 368           free(xtol);
 369      }
 370      return ret;
 371 }