api/nlopt.c

   1 #include <stdlib.h>
   2 #include <math.h>
   3
   4 #include "nlopt.h"
   5 #include "nlopt-util.h"
   6 #include "config.h"
   7
   8 #define MIN(a,b) ((a) < (b) ? (a) : (b))
   9
  10 /*************************************************************************/
  11
  12 #ifdef INFINITY
  13 #  define MY_INF INFINITY
  14 #else
  15 #  define MY_INF HUGE_VAL
  16 #endif
  17
  18 static int my_isinf(double x) {
  19      return fabs(x) >= HUGE_VAL * 0.99
  20 #ifdef HAVE_ISINF
  21           || isinf(x)
  22 #endif
  23           ;
  24 }
  25
  26 #ifndef HAVE_ISNAN
  27 static int my_isnan(double x) { return x != x; }
  28 #  define isnan my_isnan
  29 #endif
  30
  31 /*************************************************************************/
  32
  33 void nlopt_version(int *major, int *minor, int *bugfix)
  34 {
  35      *major = MAJOR_VERSION;
  36      *minor = MINOR_VERSION;
  37      *bugfix = BUGFIX_VERSION;
  38 }
  39
  40 /*************************************************************************/
  41
  42 static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = {
  43      "DIRECT (global)",
  44      "DIRECT-L (global)",
  45      "Subplex (local)",
  46      "StoGO (global)",
  47      "Low-storage BFGS (LBFGS) (local)"
  48 };
  49
  50 const char *nlopt_algorithm_name(nlopt_algorithm a)
  51 {
  52      if (a < 0 || a >= NLOPT_NUM_ALGORITHMS) return "UNKNOWN";
  53      return nlopt_algorithm_names[a];
  54 }
  55
  56 /*************************************************************************/
  57
  58 static int nlopt_srand_called = 0;
  59 void nlopt_srand(unsigned long seed) {
  60      nlopt_srand_called = 1;
  61      nlopt_init_genrand(seed);
  62 }
  63
  64 void nlopt_srand_time() {
  65      nlopt_srand(nlopt_time_seed());
  66 }
  67
  68 /*************************************************************************/
  69
  70 typedef struct {
  71      nlopt_func f;
  72      void *f_data;
  73      const double *lb, *ub;
  74 } nlopt_data;
  75
  76 #include "subplex.h"
  77
  78 static double f_subplex(int n, const double *x, void *data_)
  79 {
  80      int i;
  81      nlopt_data *data = (nlopt_data *) data_;
  82      double f;
  83
  84      /* subplex does not support bound constraints, but it supports
  85         discontinuous objectives so we can just return Inf for invalid x */
  86      for (i = 0; i < n; ++i)
  87           if (x[i] < data->lb[i] || x[i] > data->ub[i])
  88                return MY_INF;
  89
  90      f = data->f(n, x, NULL, data->f_data);
  91      return (isnan(f) ? MY_INF : f);
  92 }
  93
  94 #include "direct.h"
  95
  96 static double f_direct(int n, const double *x, int *undefined, void *data_)
  97 {
  98      nlopt_data *data = (nlopt_data *) data_;
  99      double f = data->f(n, x, NULL, data->f_data);
 100      *undefined = isnan(f) || my_isinf(f);
 101      return f;
 102 }
 103
 104 #include "stogo.h"
 105 #include "l-bfgs-b.h"
 106
 107 /*************************************************************************/
 108
 109 /* same as nlopt_minimize, but xtol_abs is required to be non-NULL */
 110 static nlopt_result nlopt_minimize_(
 111      nlopt_algorithm algorithm,
 112      int n, nlopt_func f, void *f_data,
 113      const double *lb, const double *ub, /* bounds */
 114      double *x, /* in: initial guess, out: minimizer */
 115      double *fmin, /* out: minimum */
 116      double fmin_max, double ftol_rel, double ftol_abs,
 117      double xtol_rel, const double *xtol_abs,
 118      int maxeval, double maxtime)
 119 {
 120      int i;
 121      nlopt_data d;
 122      nlopt_stopping stop;
 123
 124      d.f = f;
 125      d.f_data = f_data;
 126      d.lb = lb;
 127      d.ub = ub;
 128
 129      /* make sure rand generator is inited */
 130      if (!nlopt_srand_called)
 131           nlopt_srand_time(); /* default is non-deterministic */
 132
 133      /* check bound constraints */
 134      for (i = 0; i < n; ++i)
 135           if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
 136                return NLOPT_INVALID_ARGS;
 137
 138      stop.n = n;
 139      stop.fmin_max = (isnan(fmin_max) || (my_isinf(fmin_max) && fmin_max < 0))
 140           ? -MY_INF : fmin_max;
 141      stop.ftol_rel = ftol_rel;
 142      stop.ftol_abs = ftol_abs;
 143      stop.xtol_rel = xtol_rel;
 144      stop.xtol_abs = xtol_abs;
 145      stop.nevals = 0;
 146      stop.maxeval = maxeval;
 147      stop.maxtime = maxtime;
 148      stop.start = nlopt_seconds();
 149
 150      switch (algorithm) {
 151          case NLOPT_GLOBAL_DIRECT:
 152          case NLOPT_GLOBAL_DIRECT_L:
 153               switch (direct_optimize(f_direct, &d, n, lb, ub, x, fmin,
 154                                       maxeval, 500, ftol_rel, ftol_abs,
 155                                       xtol_rel, xtol_rel,
 156                                       DIRECT_UNKNOWN_FGLOBAL, -1.0,
 157                                       NULL,
 158                                       algorithm == NLOPT_GLOBAL_DIRECT
 159                                       ? DIRECT_ORIGINAL
 160                                       : DIRECT_GABLONSKY)) {
 161                   case DIRECT_INVALID_BOUNDS:
 162                   case DIRECT_MAXFEVAL_TOOBIG:
 163                   case DIRECT_INVALID_ARGS:
 164                        return NLOPT_INVALID_ARGS;
 165                   case DIRECT_INIT_FAILED:
 166                   case DIRECT_SAMPLEPOINTS_FAILED:
 167                   case DIRECT_SAMPLE_FAILED:
 168                        return NLOPT_FAILURE;
 169                   case DIRECT_MAXFEVAL_EXCEEDED:
 170                   case DIRECT_MAXITER_EXCEEDED:
 171                        return NLOPT_MAXEVAL_REACHED;
 172                   case DIRECT_GLOBAL_FOUND:
 173                        return NLOPT_SUCCESS;
 174                   case DIRECT_VOLTOL:
 175                   case DIRECT_SIGMATOL:
 176                        return NLOPT_XTOL_REACHED;
 177                   case DIRECT_OUT_OF_MEMORY:
 178                        return NLOPT_OUT_OF_MEMORY;
 179               }
 180               break;
 181
 182          case NLOPT_GLOBAL_STOGO:
 183               if (!stogo_minimize(n, f, f_data, x, fmin, lb, ub,
 184                                   maxeval, maxtime))
 185                    return NLOPT_FAILURE;
 186               break;
 187
 188          case NLOPT_LOCAL_SUBPLEX: {
 189               int iret;
 190               double *scale = (double *) malloc(sizeof(double) * n);
 191               if (!scale) return NLOPT_OUT_OF_MEMORY;
 192               for (i = 0; i < n; ++i) {
 193                    if (!my_isinf(ub[i]) && !my_isinf(lb[i]))
 194                         scale[i] = (ub[i] - lb[i]) * 0.01;
 195                    else if (!my_isinf(lb[i]) && x[i] > lb[i])
 196                         scale[i] = (x[i] - lb[i]) * 0.01;
 197                    else if (!my_isinf(ub[i]) && x[i] < ub[i])
 198                         scale[i] = (ub[i] - x[i]) * 0.01;
 199                    else
 200                         scale[i] = 0.01 * x[i] + 0.0001;
 201               }
 202               iret = subplex(f_subplex, fmin, x, n, &d, &stop, scale);
 203               free(scale);
 204               switch (iret) {
 205                   case -2: return NLOPT_INVALID_ARGS;
 206                   case -10: return NLOPT_MAXTIME_REACHED;
 207                   case -1: return NLOPT_MAXEVAL_REACHED;
 208                   case 0: return NLOPT_XTOL_REACHED;
 209                   case 1: return NLOPT_SUCCESS;
 210                   case 2: return NLOPT_FMIN_MAX_REACHED;
 211                   case 20: return NLOPT_FTOL_REACHED;
 212                   case -200: return NLOPT_OUT_OF_MEMORY;
 213                   default: return NLOPT_FAILURE; /* unknown return code */
 214               }
 215               break;
 216          }
 217
 218          case NLOPT_LOCAL_LBFGS: {
 219               int iret, *nbd = (int *) malloc(sizeof(int) * n);
 220               if (!nbd) return NLOPT_OUT_OF_MEMORY;
 221               for (i = 0; i < n; ++i) {
 222                    int linf = my_isinf(lb[i]) && lb[i] < 0;
 223                    int uinf = my_isinf(ub[i]) && ub[i] > 0;
 224                    nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
 225               }
 226               iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub,
 227                                      MIN(n, 5), 0.0, ftol_rel,
 228                                      xtol_abs ? *xtol_abs : xtol_rel,
 229                                      maxeval);
 230               free(nbd);
 231               if (iret <= 0) {
 232                    switch (iret) {
 233                        case -1: return NLOPT_INVALID_ARGS;
 234                        case -2: default: return NLOPT_FAILURE;
 235                    }
 236               }
 237               else {
 238                    *fmin = f(n, x, NULL, f_data);
 239                    switch (iret) {
 240                        case 5: return NLOPT_MAXEVAL_REACHED;
 241                        case 2: return NLOPT_XTOL_REACHED;
 242                        case 1: return NLOPT_FTOL_REACHED;
 243                        default: return NLOPT_SUCCESS;
 244                    }
 245               }
 246               break;
 247          }
 248
 249          default:
 250               return NLOPT_INVALID_ARGS;
 251      }
 252
 253      return NLOPT_SUCCESS;
 254 }
 255
 256 nlopt_result nlopt_minimize(
 257      nlopt_algorithm algorithm,
 258      int n, nlopt_func f, void *f_data,
 259      const double *lb, const double *ub, /* bounds */
 260      double *x, /* in: initial guess, out: minimizer */
 261      double *fmin, /* out: minimum */
 262      double fmin_max, double ftol_rel, double ftol_abs,
 263      double xtol_rel, const double *xtol_abs,
 264      int maxeval, double maxtime)
 265 {
 266      nlopt_result ret;
 267      if (xtol_abs)
 268           ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub,
 269                                 x, fmin, fmin_max, ftol_rel, ftol_abs,
 270                                 xtol_rel, xtol_abs, maxeval, maxtime);
 271      else {
 272           int i;
 273           double *xtol = (double *) malloc(sizeof(double) * n);
 274           if (!xtol) return NLOPT_OUT_OF_MEMORY;
 275           for (i = 0; i < n; ++i) xtol[i] = -1;
 276           ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub,
 277                                 x, fmin, fmin_max, ftol_rel, ftol_abs,
 278                                 xtol_rel, xtol, maxeval, maxtime);
 279           free(xtol);
 280      }
 281      return ret;
 282 }