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      "Low-storage BFGS (LBFGS) (local, derivative-based)",
  62      "Principal-axis, praxis (local, no-derivative)",
  63      "Limited-memory variable-metric, rank 1 (local, derivative-based)",
  64      "Limited-memory variable-metric, rank 2 (local, derivative-based)",
  65      "Truncated Newton (local, derivative-based)",
  66      "Truncated Newton with restarting (local, derivative-based)",
  67      "Preconditioned truncated Newton (local, derivative-based)",
  68      "Preconditioned truncated Newton with restarting (local, derivative-based)",
  69      "Controlled random search (CRS2) with local mutation (global, no-derivative)",
  70      "Multi-level single-linkage (MLSL), random (global, no-derivative)",
  71      "Multi-level single-linkage (MLSL), random (global, derivative)",
  72      "Multi-level single-linkage (MLSL), quasi-random (global, no-derivative)",
  73      "Multi-level single-linkage (MLSL), quasi-random (global, derivative)"
  74 };
  75
  76 const char *nlopt_algorithm_name(nlopt_algorithm a)
  77 {
  78      if (a < 0 || a >= NLOPT_NUM_ALGORITHMS) return "UNKNOWN";
  79      return nlopt_algorithm_names[a];
  80 }
  81
  82 /*************************************************************************/
  83
  84 static int nlopt_srand_called = 0;
  85 void nlopt_srand(unsigned long seed) {
  86      nlopt_srand_called = 1;
  87      nlopt_init_genrand(seed);
  88 }
  89
  90 void nlopt_srand_time() {
  91      nlopt_srand(nlopt_time_seed());
  92 }
  93
  94 /*************************************************************************/
  95
  96 typedef struct {
  97      nlopt_func f;
  98      void *f_data;
  99      const double *lb, *ub;
 100 } nlopt_data;
 101
 102 #include "subplex.h"
 103 #include "praxis.h"
 104
 105 static double f_subplex(int n, const double *x, void *data_)
 106 {
 107      int i;
 108      nlopt_data *data = (nlopt_data *) data_;
 109      double f;
 110
 111      /* subplex does not support bound constraints, but it supports
 112         discontinuous objectives so we can just return Inf for invalid x */
 113      for (i = 0; i < n; ++i)
 114           if (x[i] < data->lb[i] || x[i] > data->ub[i])
 115                return MY_INF;
 116
 117      f = data->f(n, x, NULL, data->f_data);
 118      return (isnan(f) ? MY_INF : f);
 119 }
 120
 121 #include "direct.h"
 122
 123 static double f_direct(int n, const double *x, int *undefined, void *data_)
 124 {
 125      nlopt_data *data = (nlopt_data *) data_;
 126      double f;
 127      f = data->f(n, x, NULL, data->f_data);
 128      *undefined = isnan(f) || my_isinf(f);
 129      return f;
 130 }
 131
 132 #ifdef WITH_CXX
 133 #  include "stogo.h"
 134 #endif
 135
 136 #include "cdirect.h"
 137
 138 #include "luksan.h"
 139
 140 #include "crs.h"
 141
 142 /*************************************************************************/
 143
 144 /* for "hybrid" algorithms that combine global search with some
 145    local search algorithm, most of the time we anticipate that people
 146    will stick with the default local search algorithm, so we
 147    don't add this as a parameter to nlopt_minimize.  Instead, the user
 148    can call nlopt_{set/get}_hybrid_local_algorithm to get/set the defaults. */
 149
 150 /* default local-search algorithms */
 151 static nlopt_algorithm local_search_alg_deriv = NLOPT_LD_LBFGS;
 152 static nlopt_algorithm local_search_alg_nonderiv = NLOPT_LN_SUBPLEX;
 153
 154 static int local_search_maxeval = -1; /* no maximum by default */
 155
 156 void nlopt_get_local_search_algorithm(nlopt_algorithm *deriv,
 157                                       nlopt_algorithm *nonderiv,
 158                                       int *maxeval)
 159 {
 160      *deriv = local_search_alg_deriv;
 161      *nonderiv = local_search_alg_nonderiv;
 162      *maxeval = local_search_maxeval;
 163 }
 164
 165 void nlopt_set_local_search_algorithm(nlopt_algorithm deriv,
 166                                       nlopt_algorithm nonderiv,
 167                                       int maxeval)
 168 {
 169      local_search_alg_deriv = deriv;
 170      local_search_alg_nonderiv = nonderiv;
 171      local_search_maxeval = maxeval;
 172 }
 173
 174 /*************************************************************************/
 175
 176 /* same as nlopt_minimize, but xtol_abs is required to be non-NULL */
 177 static nlopt_result nlopt_minimize_(
 178      nlopt_algorithm algorithm,
 179      int n, nlopt_func f, void *f_data,
 180      const double *lb, const double *ub, /* bounds */
 181      double *x, /* in: initial guess, out: minimizer */
 182      double *minf, /* out: minimum */
 183      double minf_max, double ftol_rel, double ftol_abs,
 184      double xtol_rel, const double *xtol_abs,
 185      int maxeval, double maxtime)
 186 {
 187      int i;
 188      nlopt_data d;
 189      nlopt_stopping stop;
 190
 191      /* some basic argument checks */
 192      if (n <= 0 || !f || !lb || !ub || !x || !minf)
 193           return NLOPT_INVALID_ARGS;
 194
 195      d.f = f;
 196      d.f_data = f_data;
 197      d.lb = lb;
 198      d.ub = ub;
 199
 200      /* make sure rand generator is inited */
 201      if (!nlopt_srand_called)
 202           nlopt_srand_time(); /* default is non-deterministic */
 203
 204      /* check bound constraints */
 205      for (i = 0; i < n; ++i)
 206           if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
 207                return NLOPT_INVALID_ARGS;
 208
 209      stop.n = n;
 210      stop.minf_max = (isnan(minf_max) || (my_isinf(minf_max) && minf_max < 0))
 211           ? -MY_INF : minf_max;
 212      stop.ftol_rel = ftol_rel;
 213      stop.ftol_abs = ftol_abs;
 214      stop.xtol_rel = xtol_rel;
 215      stop.xtol_abs = xtol_abs;
 216      stop.nevals = 0;
 217      stop.maxeval = maxeval;
 218      stop.maxtime = maxtime;
 219      stop.start = nlopt_seconds();
 220
 221      switch (algorithm) {
 222          case NLOPT_GN_DIRECT:
 223          case NLOPT_GN_DIRECT_L:
 224          case NLOPT_GN_DIRECT_L_RAND:
 225               return cdirect(n, f, f_data, lb, ub, x, minf, &stop, 0.0,
 226                              (algorithm != NLOPT_GN_DIRECT)
 227                              + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
 228                              + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
 229
 230          case NLOPT_GN_DIRECT_NOSCAL:
 231          case NLOPT_GN_DIRECT_L_NOSCAL:
 232          case NLOPT_GN_DIRECT_L_RAND_NOSCAL:
 233               return cdirect_unscaled(n, f, f_data, lb, ub, x, minf,
 234                                       &stop, 0.0,
 235                                       (algorithm != NLOPT_GN_DIRECT)
 236                                       + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
 237                                       + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
 238
 239          case NLOPT_GN_ORIG_DIRECT:
 240          case NLOPT_GN_ORIG_DIRECT_L:
 241               switch (direct_optimize(f_direct, &d, n, lb, ub, x, minf,
 242                                       maxeval, -1, 0.0, 0.0,
 243                                       pow(xtol_rel, (double) n), -1.0,
 244                                       stop.minf_max, 0.0,
 245                                       NULL,
 246                                       algorithm == NLOPT_GN_ORIG_DIRECT
 247                                       ? DIRECT_ORIGINAL
 248                                       : DIRECT_GABLONSKY)) {
 249                   case DIRECT_INVALID_BOUNDS:
 250                   case DIRECT_MAXFEVAL_TOOBIG:
 251                   case DIRECT_INVALID_ARGS:
 252                        return NLOPT_INVALID_ARGS;
 253                   case DIRECT_INIT_FAILED:
 254                   case DIRECT_SAMPLEPOINTS_FAILED:
 255                   case DIRECT_SAMPLE_FAILED:
 256                        return NLOPT_FAILURE;
 257                   case DIRECT_MAXFEVAL_EXCEEDED:
 258                   case DIRECT_MAXITER_EXCEEDED:
 259                        return NLOPT_MAXEVAL_REACHED;
 260                   case DIRECT_GLOBAL_FOUND:
 261                        return NLOPT_MINF_MAX_REACHED;
 262                   case DIRECT_VOLTOL:
 263                   case DIRECT_SIGMATOL:
 264                        return NLOPT_XTOL_REACHED;
 265                   case DIRECT_OUT_OF_MEMORY:
 266                        return NLOPT_OUT_OF_MEMORY;
 267               break;
 268          }
 269
 270          case NLOPT_GD_STOGO:
 271          case NLOPT_GD_STOGO_RAND:
 272 #ifdef WITH_CXX
 273               if (!stogo_minimize(n, f, f_data, x, minf, lb, ub, &stop,
 274                                   algorithm == NLOPT_GD_STOGO
 275                                   ? 0 : 2*n))
 276                    return NLOPT_FAILURE;
 277               break;
 278 #else
 279               return NLOPT_FAILURE;
 280 #endif
 281
 282          case NLOPT_LN_SUBPLEX: {
 283               int iret;
 284               double *scale = (double *) malloc(sizeof(double) * n);
 285               if (!scale) return NLOPT_OUT_OF_MEMORY;
 286               for (i = 0; i < n; ++i) {
 287                    if (!my_isinf(ub[i]) && !my_isinf(lb[i]))
 288                         scale[i] = (ub[i] - lb[i]) * 0.01;
 289                    else if (!my_isinf(lb[i]) && x[i] > lb[i])
 290                         scale[i] = (x[i] - lb[i]) * 0.01;
 291                    else if (!my_isinf(ub[i]) && x[i] < ub[i])
 292                         scale[i] = (ub[i] - x[i]) * 0.01;
 293                    else
 294                         scale[i] = 0.01 * x[i] + 0.0001;
 295               }
 296               iret = nlopt_subplex(f_subplex, minf, x, n, &d, &stop, scale);
 297               free(scale);
 298               switch (iret) {
 299                   case -2: return NLOPT_INVALID_ARGS;
 300                   case -10: return NLOPT_MAXTIME_REACHED;
 301                   case -1: return NLOPT_MAXEVAL_REACHED;
 302                   case 0: return NLOPT_XTOL_REACHED;
 303                   case 1: return NLOPT_SUCCESS;
 304                   case 2: return NLOPT_MINF_MAX_REACHED;
 305                   case 20: return NLOPT_FTOL_REACHED;
 306                   case -200: return NLOPT_OUT_OF_MEMORY;
 307                   default: return NLOPT_FAILURE; /* unknown return code */
 308               }
 309               break;
 310          }
 311
 312          case NLOPT_LN_PRAXIS: {
 313               double h0 = HUGE_VAL;
 314               for (i = 0; i < n; ++i) {
 315                    if (!my_isinf(ub[i]) && !my_isinf(lb[i]))
 316                         h0 = MIN(h0, (ub[i] - lb[i]) * 0.01);
 317                    else if (!my_isinf(lb[i]) && x[i] > lb[i])
 318                         h0 = MIN(h0, (x[i] - lb[i]) * 0.01);
 319                    else if (!my_isinf(ub[i]) && x[i] < ub[i])
 320                         h0 = MIN(h0, (ub[i] - x[i]) * 0.01);
 321                    else
 322                         h0 = MIN(h0, 0.01 * x[i] + 0.0001);
 323               }
 324               return praxis_(0.0, DBL_EPSILON, h0, n, x, f_subplex, &d,
 325                              &stop, minf);
 326          }
 327
 328          case NLOPT_LD_LBFGS:
 329               return luksan_plis(n, f, f_data, lb, ub, x, minf, &stop);
 330
 331          case NLOPT_LD_VAR1:
 332          case NLOPT_LD_VAR2:
 333               return luksan_plip(n, f, f_data, lb, ub, x, minf, &stop,
 334                    algorithm == NLOPT_LD_VAR1 ? 1 : 2);
 335
 336          case NLOPT_LD_TNEWTON:
 337          case NLOPT_LD_TNEWTON_RESTART:
 338          case NLOPT_LD_TNEWTON_PRECOND:
 339          case NLOPT_LD_TNEWTON_PRECOND_RESTART:
 340               return luksan_pnet(n, f, f_data, lb, ub, x, minf, &stop,
 341                                  1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
 342                                  1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
 343
 344          case NLOPT_GN_CRS2_LM:
 345               return crs_minimize(n, f, f_data, lb, ub, x, minf, &stop, 0);
 346
 347          case NLOPT_GN_MLSL:
 348          case NLOPT_GD_MLSL:
 349          case NLOPT_GN_MLSL_LDS:
 350          case NLOPT_GD_MLSL_LDS:
 351               return mlsl_minimize(n, f, f_data, lb, ub, x, minf, &stop,
 352                                    (algorithm == NLOPT_GN_MLSL ||
 353                                     algorithm == NLOPT_GN_MLSL_LDS)
 354                                    ? local_search_alg_nonderiv
 355                                    : local_search_alg_deriv,
 356                                    local_search_maxeval,
 357                                    algorithm >= NLOPT_GN_MLSL_LDS);
 358
 359          default:
 360               return NLOPT_INVALID_ARGS;
 361      }
 362
 363      return NLOPT_SUCCESS;
 364 }
 365
 366 nlopt_result nlopt_minimize(
 367      nlopt_algorithm algorithm,
 368      int n, nlopt_func f, void *f_data,
 369      const double *lb, const double *ub, /* bounds */
 370      double *x, /* in: initial guess, out: minimizer */
 371      double *minf, /* out: minimum */
 372      double minf_max, double ftol_rel, double ftol_abs,
 373      double xtol_rel, const double *xtol_abs,
 374      int maxeval, double maxtime)
 375 {
 376      nlopt_result ret;
 377      if (xtol_abs)
 378           ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub,
 379                                 x, minf, minf_max, ftol_rel, ftol_abs,
 380                                 xtol_rel, xtol_abs, maxeval, maxtime);
 381      else {
 382           int i;
 383           double *xtol = (double *) malloc(sizeof(double) * n);
 384           if (!xtol) return NLOPT_OUT_OF_MEMORY;
 385           for (i = 0; i < n; ++i) xtol[i] = -1;
 386           ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub,
 387                                 x, minf, minf_max, ftol_rel, ftol_abs,
 388                                 xtol_rel, xtol, maxeval, maxtime);
 389           free(xtol);
 390      }
 391      return ret;
 392 }