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