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