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