api/nlopt.c

   1 /* Copyright (c) 2007-2008 Massachusetts Institute of Technology
   2  *
   3  * Permission is hereby granted, free of charge, to any person obtaining
   4  * a copy of this software and associated documentation files (the
   5  * "Software"), to deal in the Software without restriction, including
   6  * without limitation the rights to use, copy, modify, merge, publish,
   7  * distribute, sublicense, and/or sell copies of the Software, and to
   8  * permit persons to whom the Software is furnished to do so, subject to
   9  * the following conditions:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  16  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  17  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  18  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  19  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  20  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  21  */
  22
  23 #include <stdlib.h>
  24 #include <math.h>
  25 #include <string.h>
  26 #include <float.h>
  27
  28 #include "nlopt.h"
  29 #include "nlopt-util.h"
  30 #include "config.h"
  31
  32 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  33
  34 /*************************************************************************/
  35
  36 #ifdef INFINITY
  37 #  define MY_INF INFINITY
  38 #else
  39 #  define MY_INF HUGE_VAL
  40 #endif
  41
  42 int nlopt_isinf(double x) {
  43      return fabs(x) >= HUGE_VAL * 0.99
  44 #ifdef HAVE_ISINF
  45           || isinf(x)
  46 #endif
  47           ;
  48 }
  49
  50 #ifndef HAVE_ISNAN
  51 static int my_isnan(double x) { return x != x; }
  52 #  define isnan my_isnan
  53 #endif
  54
  55 /*************************************************************************/
  56
  57 void nlopt_version(int *major, int *minor, int *bugfix)
  58 {
  59      *major = MAJOR_VERSION;
  60      *minor = MINOR_VERSION;
  61      *bugfix = BUGFIX_VERSION;
  62 }
  63
  64 /*************************************************************************/
  65
  66 static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = {
  67      "DIRECT (global, no-derivative)",
  68      "DIRECT-L (global, no-derivative)",
  69      "Randomized DIRECT-L (global, no-derivative)",
  70      "Unscaled DIRECT (global, no-derivative)",
  71      "Unscaled DIRECT-L (global, no-derivative)",
  72      "Unscaled Randomized DIRECT-L (global, no-derivative)",
  73      "Original DIRECT version (global, no-derivative)",
  74      "Original DIRECT-L version (global, no-derivative)",
  75      "Subplex (local, no-derivative)",
  76 #ifdef WITH_CXX
  77      "StoGO (global, derivative-based)",
  78      "StoGO with randomized search (global, derivative-based)",
  79 #else
  80      "StoGO (NOT COMPILED)",
  81      "StoGO randomized (NOT COMPILED)",
  82 #endif
  83 #ifdef WITH_NOCEDAL_LBFGS
  84      "original NON-FREE L-BFGS code by Nocedal et al. (local, deriv.-based)",
  85 #else
  86      "original NON-FREE L-BFGS code by Nocedal et al. (NOT COMPILED)",
  87 #endif
  88      "Limited-memory BFGS (L-BFGS) (local, derivative-based)",
  89      "Principal-axis, praxis (local, no-derivative)",
  90      "Limited-memory variable-metric, rank 1 (local, derivative-based)",
  91      "Limited-memory variable-metric, rank 2 (local, derivative-based)",
  92      "Truncated Newton (local, derivative-based)",
  93      "Truncated Newton with restarting (local, derivative-based)",
  94      "Preconditioned truncated Newton (local, derivative-based)",
  95      "Preconditioned truncated Newton with restarting (local, derivative-based)",
  96      "Controlled random search (CRS2) with local mutation (global, no-derivative)",
  97      "Multi-level single-linkage (MLSL), random (global, no-derivative)",
  98      "Multi-level single-linkage (MLSL), random (global, derivative)",
  99      "Multi-level single-linkage (MLSL), quasi-random (global, no-derivative)",
 100      "Multi-level single-linkage (MLSL), quasi-random (global, derivative)",
 101      "Method of Moving Asymptotes (MMA) (local, derivative)"
 102 };
 103
 104 const char *nlopt_algorithm_name(nlopt_algorithm a)
 105 {
 106      if (a < 0 || a >= NLOPT_NUM_ALGORITHMS) return "UNKNOWN";
 107      return nlopt_algorithm_names[a];
 108 }
 109
 110 /*************************************************************************/
 111
 112 static int nlopt_srand_called = 0;
 113 void nlopt_srand(unsigned long seed) {
 114      nlopt_srand_called = 1;
 115      nlopt_init_genrand(seed);
 116 }
 117
 118 void nlopt_srand_time() {
 119      nlopt_srand(nlopt_time_seed());
 120 }
 121
 122 /*************************************************************************/
 123
 124 typedef struct {
 125      nlopt_func f;
 126      void *f_data;
 127      const double *lb, *ub;
 128 } nlopt_data;
 129
 130 #include "subplex.h"
 131 #include "praxis.h"
 132
 133 static double f_subplex(int n, const double *x, void *data_)
 134 {
 135      int i;
 136      nlopt_data *data = (nlopt_data *) data_;
 137      double f;
 138
 139      /* subplex does not support bound constraints, but it supports
 140         discontinuous objectives so we can just return Inf for invalid x */
 141      for (i = 0; i < n; ++i)
 142           if (x[i] < data->lb[i] || x[i] > data->ub[i])
 143                return MY_INF;
 144
 145      f = data->f(n, x, NULL, data->f_data);
 146      return (isnan(f) ? MY_INF : f);
 147 }
 148
 149 #include "direct.h"
 150
 151 static double f_direct(int n, const double *x, int *undefined, void *data_)
 152 {
 153      nlopt_data *data = (nlopt_data *) data_;
 154      double f;
 155      f = data->f(n, x, NULL, data->f_data);
 156      *undefined = isnan(f) || nlopt_isinf(f);
 157      return f;
 158 }
 159
 160 #ifdef WITH_CXX
 161 #  include "stogo.h"
 162 #endif
 163
 164 #include "cdirect.h"
 165
 166 #ifdef WITH_NOCEDAL
 167 #  include "l-bfgs-b.h"
 168 #endif
 169
 170 #include "luksan.h"
 171
 172 #include "crs.h"
 173
 174 #include "mlsl.h"
 175 #include "mma.h"
 176
 177 /*************************************************************************/
 178
 179 /* for "hybrid" algorithms that combine global search with some
 180    local search algorithm, most of the time we anticipate that people
 181    will stick with the default local search algorithm, so we
 182    don't add this as a parameter to nlopt_minimize.  Instead, the user
 183    can call nlopt_{set/get}_hybrid_local_algorithm to get/set the defaults. */
 184
 185 /* default local-search algorithms */
 186 static nlopt_algorithm local_search_alg_deriv = NLOPT_LD_LBFGS;
 187 static nlopt_algorithm local_search_alg_nonderiv = NLOPT_LN_SUBPLEX;
 188
 189 static int local_search_maxeval = -1; /* no maximum by default */
 190
 191 void nlopt_get_local_search_algorithm(nlopt_algorithm *deriv,
 192                                       nlopt_algorithm *nonderiv,
 193                                       int *maxeval)
 194 {
 195      *deriv = local_search_alg_deriv;
 196      *nonderiv = local_search_alg_nonderiv;
 197      *maxeval = local_search_maxeval;
 198 }
 199
 200 void nlopt_set_local_search_algorithm(nlopt_algorithm deriv,
 201                                       nlopt_algorithm nonderiv,
 202                                       int maxeval)
 203 {
 204      local_search_alg_deriv = deriv;
 205      local_search_alg_nonderiv = nonderiv;
 206      local_search_maxeval = maxeval;
 207 }
 208
 209 /*************************************************************************/
 210
 211 /* same as nlopt_minimize, but xtol_abs is required to be non-NULL */
 212 static nlopt_result nlopt_minimize_(
 213      nlopt_algorithm algorithm,
 214      int n, nlopt_func f, void *f_data,
 215      int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size,
 216      const double *lb, const double *ub, /* bounds */
 217      double *x, /* in: initial guess, out: minimizer */
 218      double *minf, /* out: minimum */
 219      double minf_max, double ftol_rel, double ftol_abs,
 220      double xtol_rel, const double *xtol_abs,
 221      int maxeval, double maxtime)
 222 {
 223      int i;
 224      nlopt_data d;
 225      nlopt_stopping stop;
 226
 227      /* some basic argument checks */
 228      if (!minf || !f || n < 0 || m < 0
 229           || (m > 0 && !fc)) return NLOPT_INVALID_ARGS;
 230      if (n == 0) { /* trivial case: no degrees of freedom */
 231           *minf = f(n, x, NULL, f_data);
 232           return NLOPT_SUCCESS;
 233      }
 234      else if (n < 0 || !lb || !ub || !x)
 235           return NLOPT_INVALID_ARGS;
 236
 237      /* nonlinear constraints are only supported with MMA */
 238      if (m != 0 && algorithm != NLOPT_LD_MMA) return NLOPT_INVALID_ARGS;
 239
 240      d.f = f;
 241      d.f_data = f_data;
 242      d.lb = lb;
 243      d.ub = ub;
 244
 245      /* make sure rand generator is inited */
 246      if (!nlopt_srand_called)
 247           nlopt_srand_time(); /* default is non-deterministic */
 248
 249      /* check bound constraints */
 250      for (i = 0; i < n; ++i)
 251           if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
 252                return NLOPT_INVALID_ARGS;
 253
 254      stop.n = n;
 255      stop.minf_max = (isnan(minf_max)
 256                       || (nlopt_isinf(minf_max) && minf_max < 0))
 257           ? -MY_INF : minf_max;
 258      stop.ftol_rel = ftol_rel;
 259      stop.ftol_abs = ftol_abs;
 260      stop.xtol_rel = xtol_rel;
 261      stop.xtol_abs = xtol_abs;
 262      stop.nevals = 0;
 263      stop.maxeval = maxeval;
 264      stop.maxtime = maxtime;
 265      stop.start = nlopt_seconds();
 266
 267      switch (algorithm) {
 268          case NLOPT_GN_DIRECT:
 269          case NLOPT_GN_DIRECT_L:
 270          case NLOPT_GN_DIRECT_L_RAND:
 271               return cdirect(n, f, f_data, lb, ub, x, minf, &stop, 0.0,
 272                              (algorithm != NLOPT_GN_DIRECT)
 273                              + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
 274                              + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
 275
 276          case NLOPT_GN_DIRECT_NOSCAL:
 277          case NLOPT_GN_DIRECT_L_NOSCAL:
 278          case NLOPT_GN_DIRECT_L_RAND_NOSCAL:
 279               return cdirect_unscaled(n, f, f_data, lb, ub, x, minf,
 280                                       &stop, 0.0,
 281                                       (algorithm != NLOPT_GN_DIRECT)
 282                                       + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
 283                                       + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
 284
 285          case NLOPT_GN_ORIG_DIRECT:
 286          case NLOPT_GN_ORIG_DIRECT_L:
 287               switch (direct_optimize(f_direct, &d, n, lb, ub, x, minf,
 288                                       maxeval, -1, 0.0, 0.0,
 289                                       pow(xtol_rel, (double) n), -1.0,
 290                                       stop.minf_max, 0.0,
 291                                       NULL,
 292                                       algorithm == NLOPT_GN_ORIG_DIRECT
 293                                       ? DIRECT_ORIGINAL
 294                                       : DIRECT_GABLONSKY)) {
 295                   case DIRECT_INVALID_BOUNDS:
 296                   case DIRECT_MAXFEVAL_TOOBIG:
 297                   case DIRECT_INVALID_ARGS:
 298                        return NLOPT_INVALID_ARGS;
 299                   case DIRECT_INIT_FAILED:
 300                   case DIRECT_SAMPLEPOINTS_FAILED:
 301                   case DIRECT_SAMPLE_FAILED:
 302                        return NLOPT_FAILURE;
 303                   case DIRECT_MAXFEVAL_EXCEEDED:
 304                   case DIRECT_MAXITER_EXCEEDED:
 305                        return NLOPT_MAXEVAL_REACHED;
 306                   case DIRECT_GLOBAL_FOUND:
 307                        return NLOPT_MINF_MAX_REACHED;
 308                   case DIRECT_VOLTOL:
 309                   case DIRECT_SIGMATOL:
 310                        return NLOPT_XTOL_REACHED;
 311                   case DIRECT_OUT_OF_MEMORY:
 312                        return NLOPT_OUT_OF_MEMORY;
 313               break;
 314          }
 315
 316          case NLOPT_GD_STOGO:
 317          case NLOPT_GD_STOGO_RAND:
 318 #ifdef WITH_CXX
 319               if (!stogo_minimize(n, f, f_data, x, minf, lb, ub, &stop,
 320                                   algorithm == NLOPT_GD_STOGO
 321                                   ? 0 : 2*n))
 322                    return NLOPT_FAILURE;
 323               break;
 324 #else
 325               return NLOPT_FAILURE;
 326 #endif
 327
 328          case NLOPT_LN_SUBPLEX: {
 329               int iret;
 330               double *scale = (double *) malloc(sizeof(double) * n);
 331               if (!scale) return NLOPT_OUT_OF_MEMORY;
 332               for (i = 0; i < n; ++i) {
 333                    if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i]))
 334                         scale[i] = (ub[i] - lb[i]) * 0.01;
 335                    else if (!nlopt_isinf(lb[i]) && x[i] > lb[i])
 336                         scale[i] = (x[i] - lb[i]) * 0.01;
 337                    else if (!nlopt_isinf(ub[i]) && x[i] < ub[i])
 338                         scale[i] = (ub[i] - x[i]) * 0.01;
 339                    else
 340                         scale[i] = 0.01 * x[i] + 0.0001;
 341               }
 342               iret = nlopt_subplex(f_subplex, minf, x, n, &d, &stop, scale);
 343               free(scale);
 344               switch (iret) {
 345                   case -2: return NLOPT_INVALID_ARGS;
 346                   case -10: return NLOPT_MAXTIME_REACHED;
 347                   case -1: return NLOPT_MAXEVAL_REACHED;
 348                   case 0: return NLOPT_XTOL_REACHED;
 349                   case 1: return NLOPT_SUCCESS;
 350                   case 2: return NLOPT_MINF_MAX_REACHED;
 351                   case 20: return NLOPT_FTOL_REACHED;
 352                   case -200: return NLOPT_OUT_OF_MEMORY;
 353                   default: return NLOPT_FAILURE; /* unknown return code */
 354               }
 355               break;
 356          }
 357
 358          case NLOPT_LN_PRAXIS: {
 359               double h0 = HUGE_VAL;
 360               for (i = 0; i < n; ++i) {
 361                    if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i]))
 362                         h0 = MIN(h0, (ub[i] - lb[i]) * 0.01);
 363                    else if (!nlopt_isinf(lb[i]) && x[i] > lb[i])
 364                         h0 = MIN(h0, (x[i] - lb[i]) * 0.01);
 365                    else if (!nlopt_isinf(ub[i]) && x[i] < ub[i])
 366                         h0 = MIN(h0, (ub[i] - x[i]) * 0.01);
 367                    else
 368                         h0 = MIN(h0, 0.01 * x[i] + 0.0001);
 369               }
 370               return praxis_(0.0, DBL_EPSILON, h0, n, x, f_subplex, &d,
 371                              &stop, minf);
 372          }
 373
 374 #ifdef WITH_NOCEDAL
 375          case NLOPT_LD_LBFGS_NOCEDAL: {
 376               int iret, *nbd = (int *) malloc(sizeof(int) * n);
 377               if (!nbd) return NLOPT_OUT_OF_MEMORY;
 378               for (i = 0; i < n; ++i) {
 379                    int linf = nlopt_isinf(lb[i]) && lb[i] < 0;
 380                    int uinf = nlopt_isinf(ub[i]) && ub[i] > 0;
 381                    nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
 382               }
 383               iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub,
 384                                      MIN(n, 5), 0.0, ftol_rel,
 385                                      xtol_abs ? *xtol_abs : xtol_rel,
 386                                      maxeval);
 387               free(nbd);
 388               if (iret <= 0) {
 389                    switch (iret) {
 390                        case -1: return NLOPT_INVALID_ARGS;
 391                        case -2: default: return NLOPT_FAILURE;
 392                    }
 393               }
 394               else {
 395                    *minf = f(n, x, NULL, f_data);
 396                    switch (iret) {
 397                        case 5: return NLOPT_MAXEVAL_REACHED;
 398                        case 2: return NLOPT_XTOL_REACHED;
 399                        case 1: return NLOPT_FTOL_REACHED;
 400                        default: return NLOPT_SUCCESS;
 401                    }
 402               }
 403               break;
 404          }
 405 #endif
 406
 407          case NLOPT_LD_LBFGS:
 408               return luksan_plis(n, f, f_data, lb, ub, x, minf, &stop);
 409
 410          case NLOPT_LD_VAR1:
 411          case NLOPT_LD_VAR2:
 412               return luksan_plip(n, f, f_data, lb, ub, x, minf, &stop,
 413                    algorithm == NLOPT_LD_VAR1 ? 1 : 2);
 414
 415          case NLOPT_LD_TNEWTON:
 416          case NLOPT_LD_TNEWTON_RESTART:
 417          case NLOPT_LD_TNEWTON_PRECOND:
 418          case NLOPT_LD_TNEWTON_PRECOND_RESTART:
 419               return luksan_pnet(n, f, f_data, lb, ub, x, minf, &stop,
 420                                  1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
 421                                  1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
 422
 423          case NLOPT_GN_CRS2_LM:
 424               return crs_minimize(n, f, f_data, lb, ub, x, minf, &stop, 0);
 425
 426          case NLOPT_GN_MLSL:
 427          case NLOPT_GD_MLSL:
 428          case NLOPT_GN_MLSL_LDS:
 429          case NLOPT_GD_MLSL_LDS:
 430               return mlsl_minimize(n, f, f_data, lb, ub, x, minf, &stop,
 431                                    (algorithm == NLOPT_GN_MLSL ||
 432                                     algorithm == NLOPT_GN_MLSL_LDS)
 433                                    ? local_search_alg_nonderiv
 434                                    : local_search_alg_deriv,
 435                                    local_search_maxeval,
 436                                    algorithm >= NLOPT_GN_MLSL_LDS);
 437
 438          case NLOPT_LD_MMA:
 439               return mma_minimize(n, f, f_data, m, fc, fc_data, fc_datum_size,
 440                                   lb, ub, x, minf, &stop,
 441                                   local_search_alg_deriv, 1e-12, 100000);
 442
 443          default:
 444               return NLOPT_INVALID_ARGS;
 445      }
 446
 447      return NLOPT_SUCCESS;
 448 }
 449
 450 nlopt_result nlopt_minimize_constrained(
 451      nlopt_algorithm algorithm,
 452      int n, nlopt_func f, void *f_data,
 453      int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size,
 454      const double *lb, const double *ub, /* bounds */
 455      double *x, /* in: initial guess, out: minimizer */
 456      double *minf, /* out: minimum */
 457      double minf_max, double ftol_rel, double ftol_abs,
 458      double xtol_rel, const double *xtol_abs,
 459      int maxeval, double maxtime)
 460 {
 461      nlopt_result ret;
 462      if (xtol_abs)
 463           ret = nlopt_minimize_(algorithm, n, f, f_data,
 464                                 m, fc, fc_data, fc_datum_size, lb, ub,
 465                                 x, minf, minf_max, ftol_rel, ftol_abs,
 466                                 xtol_rel, xtol_abs, maxeval, maxtime);
 467      else {
 468           int i;
 469           double *xtol = (double *) malloc(sizeof(double) * n);
 470           if (!xtol) return NLOPT_OUT_OF_MEMORY;
 471           for (i = 0; i < n; ++i) xtol[i] = -1;
 472           ret = nlopt_minimize_(algorithm, n, f, f_data,
 473                                 m, fc, fc_data, fc_datum_size, lb, ub,
 474                                 x, minf, minf_max, ftol_rel, ftol_abs,
 475                                 xtol_rel, xtol, maxeval, maxtime);
 476           free(xtol);
 477      }
 478      return ret;
 479 }
 480
 481 nlopt_result nlopt_minimize(
 482      nlopt_algorithm algorithm,
 483      int n, nlopt_func f, void *f_data,
 484      const double *lb, const double *ub, /* bounds */
 485      double *x, /* in: initial guess, out: minimizer */
 486      double *minf, /* out: minimum */
 487      double minf_max, double ftol_rel, double ftol_abs,
 488      double xtol_rel, const double *xtol_abs,
 489      int maxeval, double maxtime)
 490 {
 491      return nlopt_minimize_constrained(
 492           algorithm, n, f, f_data, 0, NULL, NULL, 0,
 493           lb, ub, x, minf, minf_max, ftol_rel, ftol_abs,
 494           xtol_rel, xtol_abs, maxeval, maxtime);
 495 }