api/optimize.c

   1 /* Copyright (c) 2007-2010 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 #include "slsqp.h"
  64
  65 /*********************************************************************/
  66
  67 static double f_bound(int n, const double *x, void *data_)
  68 {
  69      int i;
  70      nlopt_opt data = (nlopt_opt) data_;
  71      double f;
  72
  73      /* some methods do not support bound constraints, but support
  74         discontinuous objectives so we can just return Inf for invalid x */
  75      for (i = 0; i < n; ++i)
  76           if (x[i] < data->lb[i] || x[i] > data->ub[i])
  77                return HUGE_VAL;
  78
  79      f = data->f((unsigned) n, x, NULL, data->f_data);
  80      return (isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f);
  81 }
  82
  83 static double f_noderiv(int n, const double *x, void *data_)
  84 {
  85      nlopt_opt data = (nlopt_opt) data_;
  86      return data->f((unsigned) n, x, NULL, data->f_data);
  87 }
  88
  89 static double f_direct(int n, const double *x, int *undefined, void *data_)
  90 {
  91      nlopt_opt data = (nlopt_opt) data_;
  92      double f;
  93      unsigned i, j;
  94      f = data->f((unsigned) n, x, NULL, data->f_data);
  95      *undefined = isnan(f) || nlopt_isinf(f);
  96      if (nlopt_get_force_stop(data)) return f;
  97      for (i = 0; i < data->m && !*undefined; ++i) {
  98           nlopt_eval_constraint(data->work, NULL, data->fc+i, (unsigned) n, x);
  99           if (nlopt_get_force_stop(data)) return f;
 100           for (j = 0; j < data->fc[i].m; ++j)
 101                if (data->work[j] > 0)
 102                     *undefined = 1;
 103      }
 104      return f;
 105 }
 106
 107 /*********************************************************************/
 108
 109 /* get min(dx) for algorithms requiring a scalar initial step size */
 110 static nlopt_result initial_step(nlopt_opt opt, const double *x, double *step)
 111 {
 112      unsigned freedx = 0, i;
 113
 114      if (!opt->dx) {
 115           freedx = 1;
 116           if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
 117                return NLOPT_OUT_OF_MEMORY;
 118      }
 119
 120      *step = HUGE_VAL;
 121      for (i = 0; i < opt->n; ++i)
 122           if (*step > fabs(opt->dx[i]))
 123                *step = fabs(opt->dx[i]);
 124
 125      if (freedx) { free(opt->dx); opt->dx = NULL; }
 126      return NLOPT_SUCCESS;
 127 }
 128
 129 /*********************************************************************/
 130
 131 /* return true if [lb,ub] is finite in every dimension (n dimensions) */
 132 static int finite_domain(unsigned n, const double *lb, const double *ub)
 133 {
 134      unsigned i;
 135      for (i = 0; i < n; ++i)
 136           if (nlopt_isinf(ub[i] - lb[i])) return 0;
 137      return 1;
 138 }
 139
 140 /*********************************************************************/
 141 /* wrapper functions, only for derivative-free methods, that
 142    eliminate dimensions with lb == ub.   (The gradient-based methods
 143    should handle this case directly, since they operate on much
 144    larger vectors where I am loathe to make copies unnecessarily.) */
 145
 146 typedef struct {
 147      nlopt_func f;
 148      nlopt_mfunc mf;
 149      void *f_data;
 150      unsigned n; /* true dimension */
 151      double *x; /* scratch vector of length n */
 152      const double *lb, *ub; /* bounds, of length n */
 153 } elimdim_data;
 154
 155 static void *elimdim_makedata(nlopt_func f, nlopt_mfunc mf, void *f_data,
 156                               unsigned n, double *x, const double *lb,
 157                               const double *ub)
 158 {
 159      elimdim_data *d = (elimdim_data *) malloc(sizeof(elimdim_data));
 160      if (!d) return NULL;
 161      d->f = f; d->mf = mf; d->f_data = f_data; d->n = n; d->x = x;
 162      d->lb = lb; d->ub = ub;
 163      return d;
 164 }
 165
 166 static double elimdim_func(unsigned n0, const double *x0, double *grad, void *d_)
 167 {
 168      elimdim_data *d = (elimdim_data *) d_;
 169      double *x = d->x;
 170      const double *lb = d->lb, *ub = d->ub;
 171      unsigned n = d->n, i, j;
 172
 173      (void) n0; /* unused */
 174      (void) grad; /* assert: grad == NULL */
 175      for (i = j = 0; i < n; ++i) {
 176           if (lb[i] == ub[i])
 177                x[i] = lb[i];
 178           else /* assert: j < n0 */
 179                x[i] = x0[j++];
 180      }
 181      return d->f(n, x, NULL, d->f_data);
 182 }
 183
 184
 185 static void elimdim_mfunc(unsigned m, double *result,
 186                           unsigned n0, const double *x0, double *grad, void *d_)
 187 {
 188      elimdim_data *d = (elimdim_data *) d_;
 189      double *x = d->x;
 190      const double *lb = d->lb, *ub = d->ub;
 191      unsigned n = d->n, i, j;
 192
 193      (void) n0; /* unused */
 194      (void) grad; /* assert: grad == NULL */
 195      for (i = j = 0; i < n; ++i) {
 196           if (lb[i] == ub[i])
 197                x[i] = lb[i];
 198           else /* assert: j < n0 */
 199                x[i] = x0[j++];
 200      }
 201      d->mf(m, result, n, x, NULL, d->f_data);
 202 }
 203
 204 /* compute the eliminated dimension: number of dims with lb[i] != ub[i] */
 205 static unsigned elimdim_dimension(unsigned n, const double *lb, const double *ub)
 206 {
 207      unsigned n0 = 0, i;
 208      for (i = 0; i < n; ++i) n0 += lb[i] != ub[i] ? 1U : 0;
 209      return n0;
 210 }
 211
 212 /* modify v to "shrunk" version, with dimensions for lb[i] == ub[i] elim'ed */
 213 static void elimdim_shrink(unsigned n, double *v,
 214                            const double *lb, const double *ub)
 215 {
 216      unsigned i, j;
 217      if (v)
 218           for (i = j = 0; i < n; ++i)
 219                if (lb[i] != ub[i])
 220                     v[j++] = v[i];
 221 }
 222
 223 /* inverse of elimdim_shrink */
 224 static void elimdim_expand(unsigned n, double *v,
 225                            const double *lb, const double *ub)
 226 {
 227      unsigned i, j;
 228      if (v && n > 0) {
 229           j = elimdim_dimension(n, lb, ub) - 1;
 230           for (i = n - 1; i > 0; --i) {
 231                if (lb[i] != ub[i])
 232                     v[i] = v[j--];
 233                else
 234                     v[i] = lb[i];
 235           }
 236           if (lb[0] == ub[0])
 237                v[0] = lb[0];
 238      }
 239 }
 240
 241 /* given opt, create a new opt with equal-constraint dimensions eliminated */
 242 static nlopt_opt elimdim_create(nlopt_opt opt)
 243 {
 244      nlopt_opt opt0 = nlopt_copy(opt);
 245      double *x;
 246      unsigned i;
 247
 248      if (!opt0) return NULL;
 249      x = (double *) malloc(sizeof(double) * opt->n);
 250      if (opt->n && !x) { nlopt_destroy(opt0); return NULL; }
 251
 252      opt0->n = elimdim_dimension(opt->n, opt->lb, opt->ub);
 253      elimdim_shrink(opt->n, opt0->lb, opt->lb, opt->ub);
 254      elimdim_shrink(opt->n, opt0->ub, opt->lb, opt->ub);
 255      elimdim_shrink(opt->n, opt0->xtol_abs, opt->lb, opt->ub);
 256      elimdim_shrink(opt->n, opt0->dx, opt->lb, opt->ub);
 257
 258      opt0->munge_on_destroy = opt0->munge_on_copy = NULL;
 259
 260      opt0->f = elimdim_func;
 261      opt0->f_data = elimdim_makedata(opt->f, NULL, opt->f_data,
 262                                      opt->n, x, opt->lb, opt->ub);
 263      if (!opt0->f_data) goto bad;
 264
 265      for (i = 0; i < opt->m; ++i) {
 266           opt0->fc[i].f = elimdim_func;
 267           opt0->fc[i].mf = elimdim_mfunc;
 268           opt0->fc[i].f_data = elimdim_makedata(opt->fc[i].f, opt->fc[i].mf,
 269                                                 opt->fc[i].f_data,
 270                                                 opt->n, x, opt->lb, opt->ub);
 271           if (!opt0->fc[i].f_data) goto bad;
 272      }
 273
 274      for (i = 0; i < opt->p; ++i) {
 275           opt0->h[i].f = elimdim_func;
 276           opt0->h[i].mf = elimdim_mfunc;
 277           opt0->h[i].f_data = elimdim_makedata(opt->h[i].f, opt->h[i].mf,
 278                                                opt->h[i].f_data,
 279                                                opt->n, x, opt->lb, opt->ub);
 280           if (!opt0->h[i].f_data) goto bad;
 281      }
 282
 283      return opt0;
 284 bad:
 285      free(x);
 286      nlopt_destroy(opt0);
 287      return NULL;
 288 }
 289
 290 /* like nlopt_destroy, but also frees elimdim_data */
 291 static void elimdim_destroy(nlopt_opt opt)
 292 {
 293      unsigned i;
 294      if (!opt) return;
 295
 296      free(((elimdim_data*) opt->f_data)->x);
 297      free(opt->f_data); opt->f_data = NULL;
 298
 299      for (i = 0; i < opt->m; ++i) {
 300           free(opt->fc[i].f_data);
 301           opt->fc[i].f_data = NULL;
 302      }
 303      for (i = 0; i < opt->p; ++i) {
 304           free(opt->h[i].f_data);
 305           opt->h[i].f_data = NULL;
 306      }
 307
 308      nlopt_destroy(opt);
 309 }
 310
 311 /* return whether to use elimdim wrapping. */
 312 static int elimdim_wrapcheck(nlopt_opt opt)
 313 {
 314      if (!opt) return 0;
 315      if (elimdim_dimension(opt->n, opt->lb, opt->ub) == opt->n) return 0;
 316      switch (opt->algorithm) {
 317          case NLOPT_GN_DIRECT:
 318          case NLOPT_GN_DIRECT_L:
 319          case NLOPT_GN_DIRECT_L_RAND:
 320          case NLOPT_GN_DIRECT_NOSCAL:
 321          case NLOPT_GN_DIRECT_L_NOSCAL:
 322          case NLOPT_GN_DIRECT_L_RAND_NOSCAL:
 323          case NLOPT_GN_ORIG_DIRECT:
 324          case NLOPT_GN_ORIG_DIRECT_L:
 325          case NLOPT_LN_PRAXIS:
 326          case NLOPT_LN_COBYLA:
 327          case NLOPT_LN_NEWUOA:
 328          case NLOPT_LN_NEWUOA_BOUND:
 329          case NLOPT_LN_BOBYQA:
 330          case NLOPT_LN_NELDERMEAD:
 331          case NLOPT_LN_SBPLX:
 332          case NLOPT_GN_ISRES:
 333               return 1;
 334
 335          default: return 0;
 336      }
 337 }
 338
 339 /*********************************************************************/
 340
 341 #define POP(defaultpop) (opt->stochastic_population > 0 ?               \
 342                          opt->stochastic_population :                   \
 343                          (nlopt_stochastic_population > 0 ?             \
 344                           nlopt_stochastic_population : (defaultpop)))
 345
 346 /* unlike nlopt_optimize() below, only handles minimization case */
 347 static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf)
 348 {
 349      const double *lb, *ub;
 350      nlopt_algorithm algorithm;
 351      nlopt_func f; void *f_data;
 352      unsigned n, i;
 353      int ni;
 354      nlopt_stopping stop;
 355
 356      if (!opt || !x || !minf || !opt->f
 357          || opt->maximize) return NLOPT_INVALID_ARGS;
 358
 359      /* reset stopping flag */
 360      nlopt_set_force_stop(opt, 0);
 361      opt->force_stop_child = NULL;
 362
 363      /* copy a few params to local vars for convenience */
 364      n = opt->n;
 365      ni = (int) n; /* most of the subroutines take "int" arg */
 366      lb = opt->lb; ub = opt->ub;
 367      algorithm = opt->algorithm;
 368      f = opt->f; f_data = opt->f_data;
 369
 370      if (n == 0) { /* trivial case: no degrees of freedom */
 371           *minf = opt->f(n, x, NULL, opt->f_data);
 372           return NLOPT_SUCCESS;
 373      }
 374
 375      *minf = HUGE_VAL;
 376
 377      /* make sure rand generator is inited */
 378      nlopt_srand_time_default(); /* default is non-deterministic */
 379
 380      /* check bound constraints */
 381      for (i = 0; i < n; ++i)
 382           if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
 383                return NLOPT_INVALID_ARGS;
 384
 385      stop.n = n;
 386      stop.minf_max = opt->stopval;
 387      stop.ftol_rel = opt->ftol_rel;
 388      stop.ftol_abs = opt->ftol_abs;
 389      stop.xtol_rel = opt->xtol_rel;
 390      stop.xtol_abs = opt->xtol_abs;
 391      stop.nevals = 0;
 392      stop.maxeval = opt->maxeval;
 393      stop.maxtime = opt->maxtime;
 394      stop.start = nlopt_seconds();
 395      stop.force_stop = &(opt->force_stop);
 396
 397      switch (algorithm) {
 398          case NLOPT_GN_DIRECT:
 399          case NLOPT_GN_DIRECT_L:
 400          case NLOPT_GN_DIRECT_L_RAND:
 401               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
 402               return cdirect(ni, f, f_data,
 403                              lb, ub, x, minf, &stop, 0.0,
 404                              (algorithm != NLOPT_GN_DIRECT)
 405                              + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND
 406                                     ? 2 : (algorithm != NLOPT_GN_DIRECT))
 407                              + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND
 408                                     ? 1 : (algorithm != NLOPT_GN_DIRECT)));
 409
 410          case NLOPT_GN_DIRECT_NOSCAL:
 411          case NLOPT_GN_DIRECT_L_NOSCAL:
 412          case NLOPT_GN_DIRECT_L_RAND_NOSCAL:
 413               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
 414               return cdirect_unscaled(ni, f, f_data, lb, ub, x, minf,
 415                                       &stop, 0.0,
 416                                       (algorithm != NLOPT_GN_DIRECT)
 417                                       + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
 418                                       + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
 419
 420          case NLOPT_GN_ORIG_DIRECT:
 421          case NLOPT_GN_ORIG_DIRECT_L: {
 422               direct_return_code dret;
 423               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
 424               opt->work = (double*) malloc(sizeof(double) *
 425                                            nlopt_max_constraint_dim(opt->m,
 426                                                                     opt->fc));
 427               if (!opt->work) return NLOPT_OUT_OF_MEMORY;
 428               dret = direct_optimize(f_direct, opt, ni, lb, ub, x, minf,
 429                                      stop.maxeval, -1,
 430                                      stop.start, stop.maxtime,
 431                                      0.0, 0.0,
 432                                      pow(stop.xtol_rel, (double) n), -1.0,
 433                                      stop.force_stop,
 434                                      stop.minf_max, 0.0,
 435                                      NULL,
 436                                      algorithm == NLOPT_GN_ORIG_DIRECT
 437                                      ? DIRECT_ORIGINAL
 438                                      : DIRECT_GABLONSKY);
 439               free(opt->work); opt->work = NULL;
 440               switch (dret) {
 441                   case DIRECT_INVALID_BOUNDS:
 442                   case DIRECT_MAXFEVAL_TOOBIG:
 443                   case DIRECT_INVALID_ARGS:
 444                        return NLOPT_INVALID_ARGS;
 445                   case DIRECT_INIT_FAILED:
 446                   case DIRECT_SAMPLEPOINTS_FAILED:
 447                   case DIRECT_SAMPLE_FAILED:
 448                        return NLOPT_FAILURE;
 449                   case DIRECT_MAXFEVAL_EXCEEDED:
 450                   case DIRECT_MAXITER_EXCEEDED:
 451                        return NLOPT_MAXEVAL_REACHED;
 452                   case DIRECT_MAXTIME_EXCEEDED:
 453                        return NLOPT_MAXTIME_REACHED;
 454                   case DIRECT_GLOBAL_FOUND:
 455                        return NLOPT_MINF_MAX_REACHED;
 456                   case DIRECT_VOLTOL:
 457                   case DIRECT_SIGMATOL:
 458                        return NLOPT_XTOL_REACHED;
 459                   case DIRECT_OUT_OF_MEMORY:
 460                        return NLOPT_OUT_OF_MEMORY;
 461                   case DIRECT_FORCED_STOP:
 462                        return NLOPT_FORCED_STOP;
 463               }
 464               break;
 465          }
 466
 467          case NLOPT_GD_STOGO:
 468          case NLOPT_GD_STOGO_RAND:
 469 #ifdef WITH_CXX
 470               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
 471               if (!stogo_minimize(ni, f, f_data, x, minf, lb, ub, &stop,
 472                                   algorithm == NLOPT_GD_STOGO
 473                                   ? 0 : (int) POP(2*n)))
 474                    return NLOPT_FAILURE;
 475               break;
 476 #else
 477               return NLOPT_INVALID_ARGS;
 478 #endif
 479
 480 #if 0
 481               /* lacking a free/open-source license, we no longer use
 482                  Rowan's code, and instead use by "sbplx" re-implementation */
 483          case NLOPT_LN_SUBPLEX: {
 484               int iret, freedx = 0;
 485               if (!opt->dx) {
 486                    freedx = 1;
 487                    if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
 488                         return NLOPT_OUT_OF_MEMORY;
 489               }
 490               iret = nlopt_subplex(f_bound, minf, x, n, opt, &stop, opt->dx);
 491               if (freedx) { free(opt->dx); opt->dx = NULL; }
 492               switch (iret) {
 493                   case -2: return NLOPT_INVALID_ARGS;
 494                   case -20: return NLOPT_FORCED_STOP;
 495                   case -10: return NLOPT_MAXTIME_REACHED;
 496                   case -1: return NLOPT_MAXEVAL_REACHED;
 497                   case 0: return NLOPT_XTOL_REACHED;
 498                   case 1: return NLOPT_SUCCESS;
 499                   case 2: return NLOPT_MINF_MAX_REACHED;
 500                   case 20: return NLOPT_FTOL_REACHED;
 501                   case -200: return NLOPT_OUT_OF_MEMORY;
 502                   default: return NLOPT_FAILURE; /* unknown return code */
 503               }
 504               break;
 505          }
 506 #endif
 507
 508          case NLOPT_LN_PRAXIS: {
 509               double step;
 510               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
 511                    return NLOPT_OUT_OF_MEMORY;
 512               return praxis_(0.0, DBL_EPSILON,
 513                              step, ni, x, f_bound, opt, &stop, minf);
 514          }
 515
 516 #ifdef WITH_NOCEDAL
 517          case NLOPT_LD_LBFGS_NOCEDAL: {
 518               int iret, *nbd = (int *) malloc(sizeof(int) * n);
 519               if (!nbd) return NLOPT_OUT_OF_MEMORY;
 520               for (i = 0; i < n; ++i) {
 521                    int linf = nlopt_isinf(lb[i]) && lb[i] < 0;
 522                    int uinf = nlopt_isinf(ub[i]) && ub[i] > 0;
 523                    nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
 524               }
 525               iret = lbfgsb_minimize(ni, f, f_data, x, nbd, lb, ub,
 526                                      ni < 5 ? ni : 5, 0.0, stop.ftol_rel,
 527                                      stop.xtol_abs[0] > 0 ? stop.xtol_abs[0]
 528                                      : stop.xtol_rel,
 529                                      stop.maxeval);
 530               free(nbd);
 531               if (iret <= 0) {
 532                    switch (iret) {
 533                        case -1: return NLOPT_INVALID_ARGS;
 534                        case -2: default: return NLOPT_FAILURE;
 535                    }
 536               }
 537               else {
 538                    *minf = f(n, x, NULL, f_data);
 539                    switch (iret) {
 540                        case 5: return NLOPT_MAXEVAL_REACHED;
 541                        case 2: return NLOPT_XTOL_REACHED;
 542                        case 1: return NLOPT_FTOL_REACHED;
 543                        default: return NLOPT_SUCCESS;
 544                    }
 545               }
 546               break;
 547          }
 548 #endif
 549
 550          case NLOPT_LD_LBFGS:
 551               return luksan_plis(ni, f, f_data, lb, ub, x, minf,
 552                                  &stop, opt->vector_storage);
 553
 554          case NLOPT_LD_VAR1:
 555          case NLOPT_LD_VAR2:
 556               return luksan_plip(ni, f, f_data, lb, ub, x, minf,
 557                                  &stop, opt->vector_storage,
 558                                  algorithm == NLOPT_LD_VAR1 ? 1 : 2);
 559
 560          case NLOPT_LD_TNEWTON:
 561          case NLOPT_LD_TNEWTON_RESTART:
 562          case NLOPT_LD_TNEWTON_PRECOND:
 563          case NLOPT_LD_TNEWTON_PRECOND_RESTART:
 564               return luksan_pnet(ni, f, f_data, lb, ub, x, minf,
 565                                  &stop, opt->vector_storage,
 566                                  1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
 567                                  1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
 568
 569          case NLOPT_GN_CRS2_LM:
 570               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
 571               return crs_minimize(ni, f, f_data, lb, ub, x, minf, &stop,
 572                                   (int) POP(0), 0);
 573
 574          case NLOPT_G_MLSL:
 575          case NLOPT_G_MLSL_LDS:
 576          case NLOPT_GN_MLSL:
 577          case NLOPT_GD_MLSL:
 578          case NLOPT_GN_MLSL_LDS:
 579          case NLOPT_GD_MLSL_LDS: {
 580               nlopt_opt local_opt = opt->local_opt;
 581               nlopt_result ret;
 582               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
 583               if (!local_opt && (algorithm == NLOPT_G_MLSL
 584                                  || algorithm == NLOPT_G_MLSL_LDS))
 585                    return NLOPT_INVALID_ARGS;
 586               if (!local_opt) { /* default */
 587                    nlopt_algorithm local_alg = (algorithm == NLOPT_GN_MLSL ||
 588                                                 algorithm == NLOPT_GN_MLSL_LDS)
 589                         ? nlopt_local_search_alg_nonderiv
 590                         : nlopt_local_search_alg_deriv;
 591                    /* don't call MLSL recursively! */
 592                    if (local_alg >= NLOPT_GN_MLSL
 593                        && local_alg <= NLOPT_GD_MLSL_LDS)
 594                         local_alg = (algorithm == NLOPT_GN_MLSL ||
 595                                      algorithm == NLOPT_GN_MLSL_LDS)
 596                              ? NLOPT_LN_COBYLA : NLOPT_LD_MMA;
 597                    local_opt = nlopt_create(local_alg, n);
 598                    if (!local_opt) return NLOPT_FAILURE;
 599                    nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
 600                    nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
 601                    nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
 602                    nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
 603                    nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
 604               }
 605               if (opt->dx) nlopt_set_initial_step(local_opt, opt->dx);
 606               for (i = 0; i < n && stop.xtol_abs[i] > 0; ++i) ;
 607               if (local_opt->ftol_rel <= 0 && local_opt->ftol_abs <= 0 &&
 608                   local_opt->xtol_rel <= 0 && i < n) {
 609                    /* it is not sensible to call MLSL without *some*
 610                       nonzero tolerance for the local search */
 611                    nlopt_set_ftol_rel(local_opt, 1e-15);
 612                    nlopt_set_xtol_rel(local_opt, 1e-7);
 613               }
 614               opt->force_stop_child = local_opt;
 615               ret = mlsl_minimize(ni, f, f_data, lb, ub, x, minf, &stop,
 616                                   local_opt, (int) POP(0),
 617                                   algorithm >= NLOPT_GN_MLSL_LDS &&
 618                                   algorithm != NLOPT_G_MLSL);
 619               opt->force_stop_child = NULL;
 620               if (!opt->local_opt) nlopt_destroy(local_opt);
 621               return ret;
 622          }
 623
 624          case NLOPT_LD_MMA: {
 625               nlopt_opt dual_opt;
 626               nlopt_result ret;
 627 #define LO(param, def) (opt->local_opt ? opt->local_opt->param : (def))
 628               dual_opt = nlopt_create(LO(algorithm,
 629                                          nlopt_local_search_alg_deriv),
 630                                       nlopt_count_constraints(opt->m,
 631                                                               opt->fc));
 632               if (!dual_opt) return NLOPT_FAILURE;
 633               nlopt_set_ftol_rel(dual_opt, LO(ftol_rel, 1e-12));
 634               nlopt_set_ftol_abs(dual_opt, LO(ftol_abs, 0.0));
 635               nlopt_set_maxeval(dual_opt, LO(maxeval, 100000));
 636 #undef LO
 637
 638               ret = mma_minimize(n, f, f_data, opt->m, opt->fc,
 639                                  lb, ub, x, minf, &stop, dual_opt);
 640               nlopt_destroy(dual_opt);
 641               return ret;
 642          }
 643
 644          case NLOPT_LN_COBYLA: {
 645               nlopt_result ret;
 646               int freedx = 0;
 647               if (!opt->dx) {
 648                    freedx = 1;
 649                    if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
 650                         return NLOPT_OUT_OF_MEMORY;
 651               }
 652               return cobyla_minimize(n, f, f_data,
 653                                      opt->m, opt->fc,
 654                                      opt->p, opt->h,
 655                                      lb, ub, x, minf, &stop,
 656                                      opt->dx);
 657               if (freedx) { free(opt->dx); opt->dx = NULL; }
 658               return ret;
 659          }
 660
 661          case NLOPT_LN_NEWUOA: {
 662               double step;
 663               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
 664                    return NLOPT_OUT_OF_MEMORY;
 665               return newuoa(ni, 2*n+1, x, 0, 0, step,
 666                             &stop, minf, f_noderiv, opt);
 667          }
 668
 669          case NLOPT_LN_NEWUOA_BOUND: {
 670               double step;
 671               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
 672                    return NLOPT_OUT_OF_MEMORY;
 673               return newuoa(ni, 2*n+1, x, lb, ub, step,
 674                             &stop, minf, f_noderiv, opt);
 675          }
 676
 677          case NLOPT_LN_BOBYQA: {
 678               nlopt_result ret;
 679               int freedx = 0;
 680               if (!opt->dx) {
 681                    freedx = 1;
 682                    if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
 683                         return NLOPT_OUT_OF_MEMORY;
 684               }
 685               ret = bobyqa(ni, 2*n+1, x, lb, ub, opt->dx,
 686                            &stop, minf, opt->f, opt->f_data);
 687               if (freedx) { free(opt->dx); opt->dx = NULL; }
 688               return ret;
 689          }
 690
 691          case NLOPT_LN_NELDERMEAD:
 692          case NLOPT_LN_SBPLX:
 693          {
 694               nlopt_result ret;
 695               int freedx = 0;
 696               if (!opt->dx) {
 697                    freedx = 1;
 698                    if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
 699                         return NLOPT_OUT_OF_MEMORY;
 700               }
 701               if (algorithm == NLOPT_LN_NELDERMEAD)
 702                    ret= nldrmd_minimize(ni,f,f_data,lb,ub,x,minf,opt->dx,&stop);
 703               else
 704                    ret= sbplx_minimize(ni,f,f_data,lb,ub,x,minf,opt->dx,&stop);
 705               if (freedx) { free(opt->dx); opt->dx = NULL; }
 706               return ret;
 707          }
 708
 709          case NLOPT_AUGLAG:
 710          case NLOPT_AUGLAG_EQ:
 711          case NLOPT_LN_AUGLAG:
 712          case NLOPT_LN_AUGLAG_EQ:
 713          case NLOPT_LD_AUGLAG:
 714          case NLOPT_LD_AUGLAG_EQ: {
 715               nlopt_opt local_opt = opt->local_opt;
 716               nlopt_result ret;
 717               if ((algorithm == NLOPT_AUGLAG || algorithm == NLOPT_AUGLAG_EQ)
 718                   && !local_opt)
 719                    return NLOPT_INVALID_ARGS;
 720               if (!local_opt) { /* default */
 721                    local_opt = nlopt_create(
 722                         algorithm == NLOPT_LN_AUGLAG ||
 723                         algorithm == NLOPT_LN_AUGLAG_EQ
 724                         ? nlopt_local_search_alg_nonderiv
 725                         : nlopt_local_search_alg_deriv, n);
 726                    if (!local_opt) return NLOPT_FAILURE;
 727                    nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
 728                    nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
 729                    nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
 730                    nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
 731                    nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
 732               }
 733               if (opt->dx) nlopt_set_initial_step(local_opt, opt->dx);
 734               opt->force_stop_child = local_opt;
 735               ret = auglag_minimize(ni, f, f_data,
 736                                     opt->m, opt->fc,
 737                                     opt->p, opt->h,
 738                                     lb, ub, x, minf, &stop,
 739                                     local_opt,
 740                                     algorithm == NLOPT_AUGLAG_EQ
 741                                     || algorithm == NLOPT_LN_AUGLAG_EQ
 742                                     || algorithm == NLOPT_LD_AUGLAG_EQ);
 743               opt->force_stop_child = NULL;
 744               if (!opt->local_opt) nlopt_destroy(local_opt);
 745               return ret;
 746          }
 747
 748          case NLOPT_GN_ISRES:
 749               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
 750               return isres_minimize(ni, f, f_data,
 751                                     (int) (opt->m), opt->fc,
 752                                     (int) (opt->p), opt->h,
 753                                     lb, ub, x, minf, &stop,
 754                                     (int) POP(0));
 755
 756          case NLOPT_LD_SLSQP:
 757               return nlopt_slsqp(n, f, f_data,
 758                                  opt->m, opt->fc,
 759                                  opt->p, opt->h,
 760                                  lb, ub, x, minf, &stop);
 761
 762          default:
 763               return NLOPT_INVALID_ARGS;
 764      }
 765
 766      return NLOPT_SUCCESS; /* never reached */
 767 }
 768
 769 /*********************************************************************/
 770
 771 typedef struct {
 772      nlopt_func f;
 773      void *f_data;
 774 } f_max_data;
 775
 776 /* wrapper for maximizing: just flip the sign of f and grad */
 777 static double f_max(unsigned n, const double *x, double *grad, void *data)
 778 {
 779      f_max_data *d = (f_max_data *) data;
 780      double val = d->f(n, x, grad, d->f_data);
 781      if (grad) {
 782           unsigned i;
 783           for (i = 0; i < n; ++i)
 784                grad[i] = -grad[i];
 785      }
 786      return -val;
 787 }
 788
 789 nlopt_result
 790 NLOPT_STDCALL nlopt_optimize(nlopt_opt opt, double *x, double *opt_f)
 791 {
 792      nlopt_func f; void *f_data;
 793      f_max_data fmd;
 794      int maximize;
 795      nlopt_result ret;
 796
 797      if (!opt || !opt_f || !opt->f) return NLOPT_INVALID_ARGS;
 798      f = opt->f; f_data = opt->f_data;
 799
 800      /* for maximizing, just minimize the f_max wrapper, which
 801         flips the sign of everything */
 802      if ((maximize = opt->maximize)) {
 803           fmd.f = f; fmd.f_data = f_data;
 804           opt->f = f_max; opt->f_data = &fmd;
 805           opt->stopval = -opt->stopval;
 806           opt->maximize = 0;
 807      }
 808
 809      { /* possibly eliminate lb == ub dimensions for some algorithms */
 810           nlopt_opt elim_opt = opt;
 811           if (elimdim_wrapcheck(opt)) {
 812                elim_opt = elimdim_create(opt);
 813                if (!elim_opt) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
 814                elimdim_shrink(opt->n, x, opt->lb, opt->ub);
 815           }
 816
 817           ret = nlopt_optimize_(elim_opt, x, opt_f);
 818
 819           if (elim_opt != opt) {
 820                elimdim_destroy(elim_opt);
 821                elimdim_expand(opt->n, x, opt->lb, opt->ub);
 822           }
 823      }
 824
 825 done:
 826      if (maximize) { /* restore original signs */
 827           opt->maximize = maximize;
 828           opt->stopval = -opt->stopval;
 829           opt->f = f; opt->f_data = f_data;
 830           *opt_f = -*opt_f;
 831      }
 832
 833      return ret;
 834 }
 835
 836 /*********************************************************************/
 837
 838 nlopt_result nlopt_optimize_limited(nlopt_opt opt, double *x, double *minf,
 839                                     int maxeval, double maxtime)
 840 {
 841      int save_maxeval;
 842      double save_maxtime;
 843      nlopt_result ret;
 844
 845      if (!opt) return NLOPT_INVALID_ARGS;
 846
 847      save_maxeval = nlopt_get_maxeval(opt);
 848      save_maxtime = nlopt_get_maxtime(opt);
 849
 850      /* override opt limits if maxeval and/or maxtime are more stringent */
 851      if (save_maxeval <= 0 || (maxeval > 0 && maxeval < save_maxeval))
 852           nlopt_set_maxeval(opt, maxeval);
 853      if (save_maxtime <= 0 || (maxtime > 0 && maxtime < save_maxtime))
 854           nlopt_set_maxtime(opt, maxtime);
 855
 856      ret = nlopt_optimize(opt, x, minf);
 857
 858      nlopt_set_maxeval(opt, save_maxeval);
 859      nlopt_set_maxtime(opt, save_maxtime);
 860
 861      return ret;
 862 }
 863
 864 /*********************************************************************/