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