api/nlopt.c

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