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 nlopt_result nlopt_minimize(
 110      nlopt_algorithm algorithm,
 111      int n, nlopt_func f, void *f_data,
 112      const double *lb, const double *ub, /* bounds */
 113      double *x, /* in: initial guess, out: minimizer */
 114      double *fmin, /* out: minimum */
 115      double fmin_max, double ftol_rel, double ftol_abs,
 116      double xtol_rel, const double *xtol_abs,
 117      int maxeval, double maxtime)
 118 {
 119      int i;
 120      nlopt_data d;
 121      d.f = f;
 122      d.f_data = f_data;
 123      d.lb = lb;
 124      d.ub = ub;
 125
 126      /* make sure rand generator is inited */
 127      if (!nlopt_srand_called)
 128           nlopt_srand_time(); /* default is non-deterministic */
 129
 130      /* check bound constraints */
 131      for (i = 0; i < n; ++i)
 132           if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
 133                return NLOPT_INVALID_ARGS;
 134
 135      switch (algorithm) {
 136          case NLOPT_GLOBAL_DIRECT:
 137          case NLOPT_GLOBAL_DIRECT_L:
 138               switch (direct_optimize(f_direct, &d, n, lb, ub, x, fmin,
 139                                       maxeval, 500, ftol_rel, ftol_abs,
 140                                       xtol_rel, xtol_rel,
 141                                       DIRECT_UNKNOWN_FGLOBAL, -1.0,
 142                                       NULL,
 143                                       algorithm == NLOPT_GLOBAL_DIRECT
 144                                       ? DIRECT_ORIGINAL
 145                                       : DIRECT_GABLONSKY)) {
 146                   case DIRECT_INVALID_BOUNDS:
 147                   case DIRECT_MAXFEVAL_TOOBIG:
 148                   case DIRECT_INVALID_ARGS:
 149                        return NLOPT_INVALID_ARGS;
 150                   case DIRECT_INIT_FAILED:
 151                   case DIRECT_SAMPLEPOINTS_FAILED:
 152                   case DIRECT_SAMPLE_FAILED:
 153                        return NLOPT_FAILURE;
 154                   case DIRECT_MAXFEVAL_EXCEEDED:
 155                   case DIRECT_MAXITER_EXCEEDED:
 156                        return NLOPT_MAXEVAL_REACHED;
 157                   case DIRECT_GLOBAL_FOUND:
 158                        return NLOPT_SUCCESS;
 159                   case DIRECT_VOLTOL:
 160                   case DIRECT_SIGMATOL:
 161                        return NLOPT_XTOL_REACHED;
 162                   case DIRECT_OUT_OF_MEMORY:
 163                        return NLOPT_OUT_OF_MEMORY;
 164               }
 165               break;
 166
 167          case NLOPT_GLOBAL_STOGO:
 168               if (!stogo_minimize(n, f, f_data, x, fmin, lb, ub,
 169                                   maxeval, maxtime))
 170                    return NLOPT_FAILURE;
 171               break;
 172
 173          case NLOPT_LOCAL_SUBPLEX: {
 174               int iret;
 175               double *scale = (double *) malloc(sizeof(double) * n);
 176               if (!scale) return NLOPT_OUT_OF_MEMORY;
 177               for (i = 0; i < n; ++i) {
 178                    if (!my_isinf(ub[i]) && !my_isinf(lb[i]))
 179                         scale[i] = (ub[i] - lb[i]) * 0.01;
 180                    else if (!my_isinf(lb[i]) && x[i] > lb[i])
 181                         scale[i] = (x[i] - lb[i]) * 0.01;
 182                    else if (!my_isinf(ub[i]) && x[i] < ub[i])
 183                         scale[i] = (ub[i] - x[i]) * 0.01;
 184                    else
 185                         scale[i] = 0.01 * x[i] + 0.0001;
 186               }
 187               iret = subplex(f_subplex, fmin, x, n, &d, xtol_rel, maxeval,
 188                              fmin_max, !my_isinf(fmin_max), scale);
 189               free(scale);
 190               switch (iret) {
 191                   case -2: return NLOPT_INVALID_ARGS;
 192                   case -1: return NLOPT_MAXEVAL_REACHED;
 193                   case 0: return NLOPT_XTOL_REACHED;
 194                   case 1: return NLOPT_SUCCESS;
 195                   case 2: return NLOPT_FMIN_MAX_REACHED;
 196               }
 197               break;
 198          }
 199
 200          case NLOPT_LOCAL_LBFGS: {
 201               int iret, *nbd = (int *) malloc(sizeof(int) * n);
 202               if (!nbd) return NLOPT_OUT_OF_MEMORY;
 203               for (i = 0; i < n; ++i) {
 204                    int linf = my_isinf(lb[i]) && lb[i] < 0;
 205                    int uinf = my_isinf(ub[i]) && ub[i] > 0;
 206                    nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
 207               }
 208               iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub,
 209                                      MIN(n, 5), 0.0, ftol_rel,
 210                                      xtol_abs ? *xtol_abs : xtol_rel,
 211                                      maxeval);
 212               free(nbd);
 213               if (iret <= 0) {
 214                    switch (iret) {
 215                        case -1: return NLOPT_INVALID_ARGS;
 216                        case -2: default: return NLOPT_FAILURE;
 217                    }
 218               }
 219               else {
 220                    *fmin = f(n, x, NULL, f_data);
 221                    switch (iret) {
 222                        case 5: return NLOPT_MAXEVAL_REACHED;
 223                        case 2: return NLOPT_XTOL_REACHED;
 224                        case 1: return NLOPT_FTOL_REACHED;
 225                        default: return NLOPT_SUCCESS;
 226                    }
 227               }
 228               break;
 229          }
 230
 231          default:
 232               return NLOPT_INVALID_ARGS;
 233      }
 234
 235      return NLOPT_SUCCESS;
 236 }