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