src/api/optimize.c

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