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