api/optimize.c

   1 /* Copyright (c) 2007-2009 Massachusetts Institute of Technology
   2  *
   3  * Permission is hereby granted, free of charge, to any person obtaining
   4  * a copy of this software and associated documentation files (the
   5  * "Software"), to deal in the Software without restriction, including
   6  * without limitation the rights to use, copy, modify, merge, publish,
   7  * distribute, sublicense, and/or sell copies of the Software, and to
   8  * permit persons to whom the Software is furnished to do so, subject to
   9  * the following conditions:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  16  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  17  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  18  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  19  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  20  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  21  */
  22
  23 #include <stdlib.h>
  24 #include <math.h>
  25 #include <float.h>
  26
  27 #include "nlopt-internal.h"
  28
  29 /*********************************************************************/
  30
  31 #ifndef HAVE_ISNAN
  32 static int my_isnan(double x) { return x != x; }
  33 #  define isnan my_isnan
  34 #endif
  35
  36 /*********************************************************************/
  37
  38 #include "praxis.h"
  39 #include "direct.h"
  40
  41 #ifdef WITH_CXX
  42 #  include "stogo.h"
  43 #endif
  44
  45 #include "cdirect.h"
  46
  47 #ifdef WITH_NOCEDAL
  48 #  include "l-bfgs-b.h"
  49 #endif
  50
  51 #include "luksan.h"
  52
  53 #include "crs.h"
  54
  55 #include "mlsl.h"
  56 #include "mma.h"
  57 #include "cobyla.h"
  58 #include "newuoa.h"
  59 #include "neldermead.h"
  60 #include "auglag.h"
  61 #include "bobyqa.h"
  62 #include "isres.h"
  63
  64 /*********************************************************************/
  65
  66 typedef struct {
  67      nlopt_func f;
  68      void *f_data;
  69      const double *lb, *ub;
  70 } nlopt_data;
  71
  72 #include "praxis.h"
  73
  74 static double f_bound(int n, const double *x, void *data_)
  75 {
  76      int i;
  77      nlopt_data *data = (nlopt_data *) data_;
  78      double f;
  79
  80      /* some methods do not support bound constraints, but support
  81         discontinuous objectives so we can just return Inf for invalid x */
  82      for (i = 0; i < n; ++i)
  83           if (x[i] < data->lb[i] || x[i] > data->ub[i])
  84                return HUGE_VAL;
  85
  86      f = data->f(n, x, NULL, data->f_data);
  87      return (isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f);
  88 }
  89
  90 static double f_noderiv(int n, const double *x, void *data_)
  91 {
  92      nlopt_data *data = (nlopt_data *) data_;
  93      return data->f(n, x, NULL, data->f_data);
  94 }
  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;
 100      f = data->f(n, x, NULL, data->f_data);
 101      *undefined = isnan(f) || nlopt_isinf(f);
 102      return f;
 103 }
 104
 105 /*********************************************************************/
 106
 107 /* get min(dx) for algorithms requiring a scalar initial step size */
 108 static nlopt_result initial_step(nlopt_opt opt, const double *x, double *step)
 109 {
 110      int freedx = 0, i;
 111
 112      if (!opt->dx) {
 113           freedx = 1;
 114           if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
 115                return NLOPT_OUT_OF_MEMORY;
 116      }
 117
 118      *step = HUGE_VAL;
 119      for (i = 0; i < opt->n; ++i)
 120           if (*step > fabs(opt->dx[i]))
 121                *step = fabs(opt->dx[i]);
 122
 123      if (freedx) { free(opt->dx); opt->dx = NULL; }
 124      return NLOPT_SUCCESS;
 125 }
 126
 127 /*********************************************************************/
 128
 129 #define POP(defaultpop) (opt->stochastic_population > 0 ?               \
 130                          opt->stochastic_population :                   \
 131                          (nlopt_stochastic_population > 0 ?             \
 132                           nlopt_stochastic_population : (defaultpop)))
 133
 134 nlopt_result nlopt_optimize(nlopt_opt opt, double *x, double *minf)
 135 {
 136      const double *lb, *ub;
 137      nlopt_algorithm algorithm;
 138      nlopt_func f; void *f_data;
 139      int n, i;
 140      nlopt_data d;
 141      nlopt_stopping stop;
 142
 143      if (!opt || !x || !minf || !opt->f) return NLOPT_INVALID_ARGS;
 144
 145      /* copy a few params to local vars for convenience */
 146      n = opt->n;
 147      lb = opt->lb; ub = opt->ub;
 148      algorithm = opt->algorithm;
 149      f = opt->f; f_data = opt->f_data;
 150
 151      if (n == 0) { /* trivial case: no degrees of freedom */
 152           *minf = f(n, x, NULL, f_data);
 153           return NLOPT_SUCCESS;
 154      }
 155
 156      *minf = HUGE_VAL;
 157
 158      d.f = f;
 159      d.f_data = f_data;
 160      d.lb = lb;
 161      d.ub = ub;
 162
 163      /* make sure rand generator is inited */
 164      if (!nlopt_srand_called)
 165           nlopt_srand_time(); /* default is non-deterministic */
 166
 167      /* check bound constraints */
 168      for (i = 0; i < n; ++i)
 169           if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
 170                return NLOPT_INVALID_ARGS;
 171
 172      stop.n = n;
 173      stop.minf_max = opt->minf_max;
 174      stop.ftol_rel = opt->ftol_rel;
 175      stop.ftol_abs = opt->ftol_abs;
 176      stop.xtol_rel = opt->xtol_rel;
 177      stop.xtol_abs = opt->xtol_abs;
 178      stop.nevals = 0;
 179      stop.maxeval = opt->maxeval;
 180      stop.maxtime = opt->maxtime;
 181      stop.start = nlopt_seconds();
 182
 183      switch (algorithm) {
 184          case NLOPT_GN_DIRECT:
 185          case NLOPT_GN_DIRECT_L:
 186          case NLOPT_GN_DIRECT_L_RAND:
 187               return cdirect(n, f, f_data,
 188                              lb, ub, x, minf, &stop, 0.0,
 189                              (algorithm != NLOPT_GN_DIRECT)
 190                              + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND
 191                                     ? 2 : (algorithm != NLOPT_GN_DIRECT))
 192                              + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND
 193                                     ? 1 : (algorithm != NLOPT_GN_DIRECT)));
 194
 195          case NLOPT_GN_DIRECT_NOSCAL:
 196          case NLOPT_GN_DIRECT_L_NOSCAL:
 197          case NLOPT_GN_DIRECT_L_RAND_NOSCAL:
 198               return cdirect_unscaled(n, f, f_data, lb, ub, x, minf,
 199                                       &stop, 0.0,
 200                                       (algorithm != NLOPT_GN_DIRECT)
 201                                       + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
 202                                       + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
 203
 204          case NLOPT_GN_ORIG_DIRECT:
 205          case NLOPT_GN_ORIG_DIRECT_L:
 206               switch (direct_optimize(f_direct, &d, n, lb, ub, x, minf,
 207                                       stop.maxeval, -1, 0.0, 0.0,
 208                                       pow(stop.xtol_rel, (double) n), -1.0,
 209                                       stop.minf_max, 0.0,
 210                                       NULL,
 211                                       algorithm == NLOPT_GN_ORIG_DIRECT
 212                                       ? DIRECT_ORIGINAL
 213                                       : DIRECT_GABLONSKY)) {
 214                   case DIRECT_INVALID_BOUNDS:
 215                   case DIRECT_MAXFEVAL_TOOBIG:
 216                   case DIRECT_INVALID_ARGS:
 217                        return NLOPT_INVALID_ARGS;
 218                   case DIRECT_INIT_FAILED:
 219                   case DIRECT_SAMPLEPOINTS_FAILED:
 220                   case DIRECT_SAMPLE_FAILED:
 221                        return NLOPT_FAILURE;
 222                   case DIRECT_MAXFEVAL_EXCEEDED:
 223                   case DIRECT_MAXITER_EXCEEDED:
 224                        return NLOPT_MAXEVAL_REACHED;
 225                   case DIRECT_GLOBAL_FOUND:
 226                        return NLOPT_MINF_MAX_REACHED;
 227                   case DIRECT_VOLTOL:
 228                   case DIRECT_SIGMATOL:
 229                        return NLOPT_XTOL_REACHED;
 230                   case DIRECT_OUT_OF_MEMORY:
 231                        return NLOPT_OUT_OF_MEMORY;
 232               break;
 233          }
 234
 235          case NLOPT_GD_STOGO:
 236          case NLOPT_GD_STOGO_RAND:
 237 #ifdef WITH_CXX
 238               if (!stogo_minimize(n, f, f_data, x, minf, lb, ub, &stop,
 239                                   algorithm == NLOPT_GD_STOGO
 240                                   ? 0 : POP(2*n)))
 241                    return NLOPT_FAILURE;
 242               break;
 243 #else
 244               return NLOPT_FAILURE;
 245 #endif
 246
 247 #if 0
 248               /* lacking a free/open-source license, we no longer use
 249                  Rowan's code, and instead use by "sbplx" re-implementation */
 250          case NLOPT_LN_SUBPLEX: {
 251               int iret, freedx = 0;
 252               if (!opt->dx) {
 253                    freedx = 1;
 254                    if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
 255                         return NLOPT_OUT_OF_MEMORY;
 256               }
 257               iret = nlopt_subplex(f_bound, minf, x, n, &d, &stop, opt->dx);
 258               if (freedx) { free(opt->dx); opt->dx = NULL; }
 259               switch (iret) {
 260                   case -2: return NLOPT_INVALID_ARGS;
 261                   case -10: return NLOPT_MAXTIME_REACHED;
 262                   case -1: return NLOPT_MAXEVAL_REACHED;
 263                   case 0: return NLOPT_XTOL_REACHED;
 264                   case 1: return NLOPT_SUCCESS;
 265                   case 2: return NLOPT_MINF_MAX_REACHED;
 266                   case 20: return NLOPT_FTOL_REACHED;
 267                   case -200: return NLOPT_OUT_OF_MEMORY;
 268                   default: return NLOPT_FAILURE; /* unknown return code */
 269               }
 270               break;
 271          }
 272 #endif
 273
 274          case NLOPT_LN_PRAXIS: {
 275               double step;
 276               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
 277                    return NLOPT_OUT_OF_MEMORY;
 278               return praxis_(0.0, DBL_EPSILON,
 279                              step, n, x, f_bound, &d, &stop, minf);
 280          }
 281
 282 #ifdef WITH_NOCEDAL
 283          case NLOPT_LD_LBFGS_NOCEDAL: {
 284               int iret, *nbd = (int *) malloc(sizeof(int) * n);
 285               if (!nbd) return NLOPT_OUT_OF_MEMORY;
 286               for (i = 0; i < n; ++i) {
 287                    int linf = nlopt_isinf(lb[i]) && lb[i] < 0;
 288                    int uinf = nlopt_isinf(ub[i]) && ub[i] > 0;
 289                    nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
 290               }
 291               iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub,
 292                                      MIN(n, 5), 0.0, stop.ftol_rel,
 293                                      stop.xtol_abs[0] > 0 ? stop.xtol_abs[0]
 294                                      : stop.xtol_rel,
 295                                      stop.maxeval);
 296               free(nbd);
 297               if (iret <= 0) {
 298                    switch (iret) {
 299                        case -1: return NLOPT_INVALID_ARGS;
 300                        case -2: default: return NLOPT_FAILURE;
 301                    }
 302               }
 303               else {
 304                    *minf = f(n, x, NULL, f_data);
 305                    switch (iret) {
 306                        case 5: return NLOPT_MAXEVAL_REACHED;
 307                        case 2: return NLOPT_XTOL_REACHED;
 308                        case 1: return NLOPT_FTOL_REACHED;
 309                        default: return NLOPT_SUCCESS;
 310                    }
 311               }
 312               break;
 313          }
 314 #endif
 315
 316          case NLOPT_LD_LBFGS:
 317               return luksan_plis(n, f, f_data, lb, ub, x, minf, &stop);
 318
 319          case NLOPT_LD_VAR1:
 320          case NLOPT_LD_VAR2:
 321               return luksan_plip(n, f, f_data, lb, ub, x, minf, &stop,
 322                    algorithm == NLOPT_LD_VAR1 ? 1 : 2);
 323
 324          case NLOPT_LD_TNEWTON:
 325          case NLOPT_LD_TNEWTON_RESTART:
 326          case NLOPT_LD_TNEWTON_PRECOND:
 327          case NLOPT_LD_TNEWTON_PRECOND_RESTART:
 328               return luksan_pnet(n, f, f_data, lb, ub, x, minf, &stop,
 329                                  1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
 330                                  1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
 331
 332          case NLOPT_GN_CRS2_LM:
 333               return crs_minimize(n, f, f_data, lb, ub, x, minf, &stop,
 334                                   POP(0), 0);
 335
 336          case NLOPT_GN_MLSL:
 337          case NLOPT_GD_MLSL:
 338          case NLOPT_GN_MLSL_LDS:
 339          case NLOPT_GD_MLSL_LDS: {
 340               nlopt_opt local_opt = opt->local_opt;
 341               nlopt_result ret;
 342               if (!local_opt) { /* default */
 343                    local_opt = nlopt_create((algorithm == NLOPT_GN_MLSL ||
 344                                              algorithm == NLOPT_GN_MLSL_LDS)
 345                                             ? nlopt_local_search_alg_nonderiv
 346                                             : nlopt_local_search_alg_deriv, n);
 347                    if (!local_opt) return NLOPT_FAILURE;
 348                    nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
 349                    nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
 350                    nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
 351                    nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
 352                    nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
 353               }
 354               for (i = 0; i < n && stop.xtol_abs[i] > 0; ++i) ;
 355               if (local_opt->ftol_rel <= 0 && local_opt->ftol_abs <= 0 &&
 356                   local_opt->xtol_rel <= 0 && i < n) {
 357                    /* it is not sensible to call MLSL without *some*
 358                       nonzero tolerance for the local search */
 359                    nlopt_set_ftol_rel(local_opt, 1e-15);
 360                    nlopt_set_xtol_rel(local_opt, 1e-7);
 361               }
 362               ret = mlsl_minimize(n, f, f_data, lb, ub, x, minf, &stop,
 363                                   local_opt, POP(0),
 364                                   algorithm >= NLOPT_GN_MLSL_LDS);
 365               if (!opt->local_opt) nlopt_destroy(local_opt);
 366               return ret;
 367          }
 368
 369          case NLOPT_LD_MMA: {
 370               nlopt_opt dual_opt;
 371               nlopt_result ret;
 372 #define LO(param, def) (opt->local_opt ? opt->local_opt->param : (def))
 373               dual_opt = nlopt_create(LO(algorithm,
 374                                          nlopt_local_search_alg_deriv),
 375                                       opt->m);
 376               if (!dual_opt) return NLOPT_FAILURE;
 377               nlopt_set_ftol_rel(dual_opt, LO(ftol_rel, 1e-12));
 378               nlopt_set_ftol_abs(dual_opt, LO(ftol_abs, 0.0));
 379               nlopt_set_maxeval(dual_opt, LO(maxeval, 100000));
 380 #undef LO
 381
 382               ret = mma_minimize(n, f, f_data, opt->m, opt->fc,
 383                                  lb, ub, x, minf, &stop, dual_opt);
 384               nlopt_destroy(dual_opt);
 385               return ret;
 386          }
 387
 388          case NLOPT_LN_COBYLA: {
 389               double step;
 390               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
 391                    return NLOPT_OUT_OF_MEMORY;
 392               return cobyla_minimize(n, f, f_data,
 393                                      opt->m, opt->fc,
 394                                      lb, ub, x, minf, &stop,
 395                                      step);
 396          }
 397
 398          case NLOPT_LN_NEWUOA: {
 399               double step;
 400               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
 401                    return NLOPT_OUT_OF_MEMORY;
 402               return newuoa(n, 2*n+1, x, 0, 0, step,
 403                             &stop, minf, f_noderiv, &d);
 404          }
 405
 406          case NLOPT_LN_NEWUOA_BOUND: {
 407               double step;
 408               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
 409                    return NLOPT_OUT_OF_MEMORY;
 410               return newuoa(n, 2*n+1, x, lb, ub, step,
 411                             &stop, minf, f_noderiv, &d);
 412          }
 413
 414          case NLOPT_LN_BOBYQA: {
 415               double step;
 416               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
 417                    return NLOPT_OUT_OF_MEMORY;
 418               return bobyqa(n, 2*n+1, x, lb, ub, step,
 419                             &stop, minf, f_noderiv, &d);
 420          }
 421
 422          case NLOPT_LN_NELDERMEAD:
 423          case NLOPT_LN_SBPLX:
 424          {
 425               nlopt_result ret;
 426               int freedx = 0;
 427               if (!opt->dx) {
 428                    freedx = 1;
 429                    if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
 430                         return NLOPT_OUT_OF_MEMORY;
 431               }
 432               if (algorithm == NLOPT_LN_NELDERMEAD)
 433                    ret= nldrmd_minimize(n,f,f_data,lb,ub,x,minf,opt->dx,&stop);
 434               else
 435                    ret= sbplx_minimize(n,f,f_data,lb,ub,x,minf,opt->dx,&stop);
 436               if (freedx) { free(opt->dx); opt->dx = NULL; }
 437               return ret;
 438          }
 439
 440          case NLOPT_LN_AUGLAG:
 441          case NLOPT_LN_AUGLAG_EQ:
 442          case NLOPT_LD_AUGLAG:
 443          case NLOPT_LD_AUGLAG_EQ: {
 444               nlopt_opt local_opt = opt->local_opt;
 445               nlopt_result ret;
 446               if (!local_opt) { /* default */
 447                    local_opt = nlopt_create(
 448                         algorithm == NLOPT_LN_AUGLAG ||
 449                         algorithm == NLOPT_LN_AUGLAG_EQ
 450                         ? nlopt_local_search_alg_nonderiv
 451                         : nlopt_local_search_alg_deriv, n);
 452                    if (!local_opt) return NLOPT_FAILURE;
 453                    nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
 454                    nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
 455                    nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
 456                    nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
 457               }
 458               ret = auglag_minimize(n, f, f_data,
 459                                     opt->m, opt->fc,
 460                                     opt->p, opt->h,
 461                                     lb, ub, x, minf, &stop,
 462                                     local_opt,
 463                                     algorithm == NLOPT_LN_AUGLAG_EQ
 464                                     || algorithm == NLOPT_LD_AUGLAG_EQ);
 465               if (!opt->local_opt) nlopt_destroy(local_opt);
 466               return ret;
 467          }
 468
 469          case NLOPT_GN_ISRES:
 470               return isres_minimize(n, f, f_data,
 471                                     opt->m, opt->fc,
 472                                     opt->p, opt->h,
 473                                     lb, ub, x, minf, &stop,
 474                                     POP(0));
 475
 476          default:
 477               return NLOPT_INVALID_ARGS;
 478      }
 479
 480      return NLOPT_SUCCESS; /* never reached */
 481 }
 482
 483 /*********************************************************************/
 484
 485 nlopt_result nlopt_optimize_limited(nlopt_opt opt, double *x, double *minf,
 486                                     int maxeval, double maxtime)
 487 {
 488      int save_maxeval;
 489      double save_maxtime;
 490      nlopt_result ret;
 491
 492      if (!opt) return NLOPT_INVALID_ARGS;
 493
 494      save_maxeval = nlopt_get_maxeval(opt);
 495      save_maxtime = nlopt_get_maxtime(opt);
 496
 497      /* override opt limits if maxeval and/or maxtime are more stringent */
 498      if (save_maxeval <= 0 || (maxeval > 0 && maxeval < save_maxeval))
 499           nlopt_set_maxeval(opt, maxeval);
 500      if (save_maxtime <= 0 || (maxtime > 0 && maxtime < save_maxtime))
 501           nlopt_set_maxtime(opt, maxtime);
 502
 503      ret = nlopt_optimize(opt, x, minf);
 504
 505      nlopt_set_maxeval(opt, save_maxeval);
 506      nlopt_set_maxtime(opt, save_maxtime);
 507
 508      return ret;
 509 }
 510
 511 /*********************************************************************/