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) return NLOPT_INVALID_ARGS;
 207      if (n == 0) { /* trivial case: no degrees of freedom */
 208           *minf = f(n, x, NULL, f_data);
 209           return NLOPT_SUCCESS;
 210      }
 211      else if (n < 0 || !lb || !ub || !x)
 212           return NLOPT_INVALID_ARGS;
 213
 214      /* nonlinear constraints are only supported with MMA */
 215      if (m != 0 && algorithm != NLOPT_LD_MMA) return NLOPT_INVALID_ARGS;
 216
 217      d.f = f;
 218      d.f_data = f_data;
 219      d.lb = lb;
 220      d.ub = ub;
 221
 222      /* make sure rand generator is inited */
 223      if (!nlopt_srand_called)
 224           nlopt_srand_time(); /* default is non-deterministic */
 225
 226      /* check bound constraints */
 227      for (i = 0; i < n; ++i)
 228           if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
 229                return NLOPT_INVALID_ARGS;
 230
 231      stop.n = n;
 232      stop.minf_max = (isnan(minf_max)
 233                       || (nlopt_isinf(minf_max) && minf_max < 0))
 234           ? -MY_INF : minf_max;
 235      stop.ftol_rel = ftol_rel;
 236      stop.ftol_abs = ftol_abs;
 237      stop.xtol_rel = xtol_rel;
 238      stop.xtol_abs = xtol_abs;
 239      stop.nevals = 0;
 240      stop.maxeval = maxeval;
 241      stop.maxtime = maxtime;
 242      stop.start = nlopt_seconds();
 243
 244      switch (algorithm) {
 245          case NLOPT_GN_DIRECT:
 246          case NLOPT_GN_DIRECT_L:
 247          case NLOPT_GN_DIRECT_L_RAND:
 248               return cdirect(n, f, f_data, lb, ub, x, minf, &stop, 0.0,
 249                              (algorithm != NLOPT_GN_DIRECT)
 250                              + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
 251                              + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
 252
 253          case NLOPT_GN_DIRECT_NOSCAL:
 254          case NLOPT_GN_DIRECT_L_NOSCAL:
 255          case NLOPT_GN_DIRECT_L_RAND_NOSCAL:
 256               return cdirect_unscaled(n, f, f_data, lb, ub, x, minf,
 257                                       &stop, 0.0,
 258                                       (algorithm != NLOPT_GN_DIRECT)
 259                                       + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
 260                                       + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
 261
 262          case NLOPT_GN_ORIG_DIRECT:
 263          case NLOPT_GN_ORIG_DIRECT_L:
 264               switch (direct_optimize(f_direct, &d, n, lb, ub, x, minf,
 265                                       maxeval, -1, 0.0, 0.0,
 266                                       pow(xtol_rel, (double) n), -1.0,
 267                                       stop.minf_max, 0.0,
 268                                       NULL,
 269                                       algorithm == NLOPT_GN_ORIG_DIRECT
 270                                       ? DIRECT_ORIGINAL
 271                                       : DIRECT_GABLONSKY)) {
 272                   case DIRECT_INVALID_BOUNDS:
 273                   case DIRECT_MAXFEVAL_TOOBIG:
 274                   case DIRECT_INVALID_ARGS:
 275                        return NLOPT_INVALID_ARGS;
 276                   case DIRECT_INIT_FAILED:
 277                   case DIRECT_SAMPLEPOINTS_FAILED:
 278                   case DIRECT_SAMPLE_FAILED:
 279                        return NLOPT_FAILURE;
 280                   case DIRECT_MAXFEVAL_EXCEEDED:
 281                   case DIRECT_MAXITER_EXCEEDED:
 282                        return NLOPT_MAXEVAL_REACHED;
 283                   case DIRECT_GLOBAL_FOUND:
 284                        return NLOPT_MINF_MAX_REACHED;
 285                   case DIRECT_VOLTOL:
 286                   case DIRECT_SIGMATOL:
 287                        return NLOPT_XTOL_REACHED;
 288                   case DIRECT_OUT_OF_MEMORY:
 289                        return NLOPT_OUT_OF_MEMORY;
 290               break;
 291          }
 292
 293          case NLOPT_GD_STOGO:
 294          case NLOPT_GD_STOGO_RAND:
 295 #ifdef WITH_CXX
 296               if (!stogo_minimize(n, f, f_data, x, minf, lb, ub, &stop,
 297                                   algorithm == NLOPT_GD_STOGO
 298                                   ? 0 : 2*n))
 299                    return NLOPT_FAILURE;
 300               break;
 301 #else
 302               return NLOPT_FAILURE;
 303 #endif
 304
 305          case NLOPT_LN_SUBPLEX: {
 306               int iret;
 307               double *scale = (double *) malloc(sizeof(double) * n);
 308               if (!scale) return NLOPT_OUT_OF_MEMORY;
 309               for (i = 0; i < n; ++i) {
 310                    if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i]))
 311                         scale[i] = (ub[i] - lb[i]) * 0.01;
 312                    else if (!nlopt_isinf(lb[i]) && x[i] > lb[i])
 313                         scale[i] = (x[i] - lb[i]) * 0.01;
 314                    else if (!nlopt_isinf(ub[i]) && x[i] < ub[i])
 315                         scale[i] = (ub[i] - x[i]) * 0.01;
 316                    else
 317                         scale[i] = 0.01 * x[i] + 0.0001;
 318               }
 319               iret = nlopt_subplex(f_subplex, minf, x, n, &d, &stop, scale);
 320               free(scale);
 321               switch (iret) {
 322                   case -2: return NLOPT_INVALID_ARGS;
 323                   case -10: return NLOPT_MAXTIME_REACHED;
 324                   case -1: return NLOPT_MAXEVAL_REACHED;
 325                   case 0: return NLOPT_XTOL_REACHED;
 326                   case 1: return NLOPT_SUCCESS;
 327                   case 2: return NLOPT_MINF_MAX_REACHED;
 328                   case 20: return NLOPT_FTOL_REACHED;
 329                   case -200: return NLOPT_OUT_OF_MEMORY;
 330                   default: return NLOPT_FAILURE; /* unknown return code */
 331               }
 332               break;
 333          }
 334
 335          case NLOPT_LN_PRAXIS: {
 336               double h0 = HUGE_VAL;
 337               for (i = 0; i < n; ++i) {
 338                    if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i]))
 339                         h0 = MIN(h0, (ub[i] - lb[i]) * 0.01);
 340                    else if (!nlopt_isinf(lb[i]) && x[i] > lb[i])
 341                         h0 = MIN(h0, (x[i] - lb[i]) * 0.01);
 342                    else if (!nlopt_isinf(ub[i]) && x[i] < ub[i])
 343                         h0 = MIN(h0, (ub[i] - x[i]) * 0.01);
 344                    else
 345                         h0 = MIN(h0, 0.01 * x[i] + 0.0001);
 346               }
 347               return praxis_(0.0, DBL_EPSILON, h0, n, x, f_subplex, &d,
 348                              &stop, minf);
 349          }
 350
 351 #ifdef WITH_NOCEDAL
 352          case NLOPT_LD_LBFGS_NOCEDAL: {
 353               int iret, *nbd = (int *) malloc(sizeof(int) * n);
 354               if (!nbd) return NLOPT_OUT_OF_MEMORY;
 355               for (i = 0; i < n; ++i) {
 356                    int linf = nlopt_isinf(lb[i]) && lb[i] < 0;
 357                    int uinf = nlopt_isinf(ub[i]) && ub[i] > 0;
 358                    nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
 359               }
 360               iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub,
 361                                      MIN(n, 5), 0.0, ftol_rel,
 362                                      xtol_abs ? *xtol_abs : xtol_rel,
 363                                      maxeval);
 364               free(nbd);
 365               if (iret <= 0) {
 366                    switch (iret) {
 367                        case -1: return NLOPT_INVALID_ARGS;
 368                        case -2: default: return NLOPT_FAILURE;
 369                    }
 370               }
 371               else {
 372                    *minf = f(n, x, NULL, f_data);
 373                    switch (iret) {
 374                        case 5: return NLOPT_MAXEVAL_REACHED;
 375                        case 2: return NLOPT_XTOL_REACHED;
 376                        case 1: return NLOPT_FTOL_REACHED;
 377                        default: return NLOPT_SUCCESS;
 378                    }
 379               }
 380               break;
 381          }
 382 #endif
 383
 384          case NLOPT_LD_LBFGS:
 385               return luksan_plis(n, f, f_data, lb, ub, x, minf, &stop);
 386
 387          case NLOPT_LD_VAR1:
 388          case NLOPT_LD_VAR2:
 389               return luksan_plip(n, f, f_data, lb, ub, x, minf, &stop,
 390                    algorithm == NLOPT_LD_VAR1 ? 1 : 2);
 391
 392          case NLOPT_LD_TNEWTON:
 393          case NLOPT_LD_TNEWTON_RESTART:
 394          case NLOPT_LD_TNEWTON_PRECOND:
 395          case NLOPT_LD_TNEWTON_PRECOND_RESTART:
 396               return luksan_pnet(n, f, f_data, lb, ub, x, minf, &stop,
 397                                  1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
 398                                  1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
 399
 400          case NLOPT_GN_CRS2_LM:
 401               return crs_minimize(n, f, f_data, lb, ub, x, minf, &stop, 0);
 402
 403          case NLOPT_GN_MLSL:
 404          case NLOPT_GD_MLSL:
 405          case NLOPT_GN_MLSL_LDS:
 406          case NLOPT_GD_MLSL_LDS:
 407               return mlsl_minimize(n, f, f_data, lb, ub, x, minf, &stop,
 408                                    (algorithm == NLOPT_GN_MLSL ||
 409                                     algorithm == NLOPT_GN_MLSL_LDS)
 410                                    ? local_search_alg_nonderiv
 411                                    : local_search_alg_deriv,
 412                                    local_search_maxeval,
 413                                    algorithm >= NLOPT_GN_MLSL_LDS);
 414
 415          case NLOPT_LD_MMA:
 416               return mma_minimize(n, f, f_data, m, fc, fc_data, fc_datum_size,
 417                                   lb, ub, x, minf, &stop,
 418                                   local_search_alg_deriv, 1e-12, 100000);
 419
 420          default:
 421               return NLOPT_INVALID_ARGS;
 422      }
 423
 424      return NLOPT_SUCCESS;
 425 }
 426
 427 nlopt_result nlopt_minimize_c(
 428      nlopt_algorithm algorithm,
 429      int n, nlopt_func f, void *f_data,
 430      int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size,
 431      const double *lb, const double *ub, /* bounds */
 432      double *x, /* in: initial guess, out: minimizer */
 433      double *minf, /* out: minimum */
 434      double minf_max, double ftol_rel, double ftol_abs,
 435      double xtol_rel, const double *xtol_abs,
 436      int maxeval, double maxtime)
 437 {
 438      nlopt_result ret;
 439      if (xtol_abs)
 440           ret = nlopt_minimize_(algorithm, n, f, f_data,
 441                                 m, fc, fc_data, fc_datum_size, lb, ub,
 442                                 x, minf, minf_max, ftol_rel, ftol_abs,
 443                                 xtol_rel, xtol_abs, maxeval, maxtime);
 444      else {
 445           int i;
 446           double *xtol = (double *) malloc(sizeof(double) * n);
 447           if (!xtol) return NLOPT_OUT_OF_MEMORY;
 448           for (i = 0; i < n; ++i) xtol[i] = -1;
 449           ret = nlopt_minimize_(algorithm, n, f, f_data,
 450                                 m, fc, fc_data, fc_datum_size, lb, ub,
 451                                 x, minf, minf_max, ftol_rel, ftol_abs,
 452                                 xtol_rel, xtol, maxeval, maxtime);
 453           free(xtol);
 454      }
 455      return ret;
 456 }
 457
 458 nlopt_result nlopt_minimize(
 459      nlopt_algorithm algorithm,
 460      int n, nlopt_func f, void *f_data,
 461      const double *lb, const double *ub, /* bounds */
 462      double *x, /* in: initial guess, out: minimizer */
 463      double *minf, /* out: minimum */
 464      double minf_max, double ftol_rel, double ftol_abs,
 465      double xtol_rel, const double *xtol_abs,
 466      int maxeval, double maxtime)
 467 {
 468      return nlopt_minimize_c(algorithm, n, f, f_data, 0, NULL, NULL, 0,
 469                              lb, ub, x, minf, minf_max, ftol_rel, ftol_abs,
 470                              xtol_rel, xtol_abs, maxeval, maxtime);
 471 }