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