1 /* Copyright (c) 2007-2008 Massachusetts Institute of Technology
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:
11 * The above copyright notice and this permission notice shall be
12 * included in all copies or substantial portions of the Software.
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.
29 #include "nlopt-util.h"
31 #define MIN(a,b) ((a) < (b) ? (a) : (b))
33 /*************************************************************************/
35 int nlopt_isinf(double x) {
36 return fabs(x) >= HUGE_VAL * 0.99
44 static int my_isnan(double x) { return x != x; }
45 # define isnan my_isnan
48 /*************************************************************************/
50 void nlopt_version(int *major, int *minor, int *bugfix)
52 *major = MAJOR_VERSION;
53 *minor = MINOR_VERSION;
54 *bugfix = BUGFIX_VERSION;
57 /*************************************************************************/
59 static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = {
60 "DIRECT (global, no-derivative)",
61 "DIRECT-L (global, no-derivative)",
62 "Randomized DIRECT-L (global, no-derivative)",
63 "Unscaled DIRECT (global, no-derivative)",
64 "Unscaled DIRECT-L (global, no-derivative)",
65 "Unscaled Randomized DIRECT-L (global, no-derivative)",
66 "Original DIRECT version (global, no-derivative)",
67 "Original DIRECT-L version (global, no-derivative)",
69 "StoGO (global, derivative-based)",
70 "StoGO with randomized search (global, derivative-based)",
72 "StoGO (NOT COMPILED)",
73 "StoGO randomized (NOT COMPILED)",
75 #ifdef WITH_NOCEDAL_LBFGS
76 "original NON-FREE L-BFGS code by Nocedal et al. (local, deriv.-based)",
78 "original NON-FREE L-BFGS code by Nocedal et al. (NOT COMPILED)",
80 "Limited-memory BFGS (L-BFGS) (local, derivative-based)",
81 "Principal-axis, praxis (local, no-derivative)",
82 "Limited-memory variable-metric, rank 1 (local, derivative-based)",
83 "Limited-memory variable-metric, rank 2 (local, derivative-based)",
84 "Truncated Newton (local, derivative-based)",
85 "Truncated Newton with restarting (local, derivative-based)",
86 "Preconditioned truncated Newton (local, derivative-based)",
87 "Preconditioned truncated Newton with restarting (local, derivative-based)",
88 "Controlled random search (CRS2) with local mutation (global, no-derivative)",
89 "Multi-level single-linkage (MLSL), random (global, no-derivative)",
90 "Multi-level single-linkage (MLSL), random (global, derivative)",
91 "Multi-level single-linkage (MLSL), quasi-random (global, no-derivative)",
92 "Multi-level single-linkage (MLSL), quasi-random (global, derivative)",
93 "Method of Moving Asymptotes (MMA) (local, derivative)",
94 "COBYLA (Constrained Optimization BY Linear Approximations) (local, no-derivative)",
95 "NEWUOA unconstrained optimization via quadratic models (local, no-derivative)",
96 "Bound-constrained optimization via NEWUOA-based quadratic models (local, no-derivative)",
97 "Nelder-Mead simplex algorithm (local, no-derivative)",
98 "Sbplx variant of Nelder-Mead (re-implementation of Rowan's Subplex) (local, no-derivative)",
99 "Augmented Lagrangian method (local, no-derivative)",
100 "Augmented Lagrangian method (local, derivative)",
101 "Augmented Lagrangian method for equality constraints (local, no-derivative)",
102 "Augmented Lagrangian method for equality constraints (local, derivative)",
105 const char *nlopt_algorithm_name(nlopt_algorithm a)
107 if (a < 0 || a >= NLOPT_NUM_ALGORITHMS) return "UNKNOWN";
108 return nlopt_algorithm_names[a];
111 /*************************************************************************/
113 static int nlopt_srand_called = 0;
114 void nlopt_srand(unsigned long seed) {
115 nlopt_srand_called = 1;
116 nlopt_init_genrand(seed);
119 void nlopt_srand_time() {
120 nlopt_srand(nlopt_time_seed());
123 /*************************************************************************/
125 /* crude heuristics for initial step size of nonderivative algorithms */
126 static double initial_step(int n,
127 const double *lb, const double *ub, const double *x)
130 double step = HUGE_VAL;
131 for (i = 0; i < n; ++i) {
132 if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i])
133 && (ub[i] - lb[i]) * 0.25 < step && ub[i] > lb[i])
134 step = (ub[i] - lb[i]) * 0.25;
135 if (!nlopt_isinf(ub[i])
136 && ub[i] - x[i] < step && ub[i] > x[i])
138 if (!nlopt_isinf(lb[i])
139 && x[i] - lb[i] < step && x[i] > lb[i])
142 if (nlopt_isinf(step))
143 for (i = 0; i < n; ++i) {
144 if (!nlopt_isinf(ub[i])
145 && ub[i] - x[i] < step && ub[i] > x[i] + 1e-4)
147 if (!nlopt_isinf(lb[i])
148 && x[i] - lb[i] < step && x[i] > lb[i] + 1e-4)
151 if (nlopt_isinf(step)) {
153 for (i = 0; i < n; ++i)
154 if (fabs(x[i]) * 0.25 > step)
155 step = fabs(x[i]) * 0.25;
162 /*************************************************************************/
167 const double *lb, *ub;
172 static double f_bound(int n, const double *x, void *data_)
175 nlopt_data *data = (nlopt_data *) data_;
178 /* some methods do not support bound constraints, but support
179 discontinuous objectives so we can just return Inf for invalid x */
180 for (i = 0; i < n; ++i)
181 if (x[i] < data->lb[i] || x[i] > data->ub[i])
184 f = data->f(n, x, NULL, data->f_data);
185 return (isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f);
188 static double f_noderiv(int n, const double *x, void *data_)
190 nlopt_data *data = (nlopt_data *) data_;
191 return data->f(n, x, NULL, data->f_data);
196 static double f_direct(int n, const double *x, int *undefined, void *data_)
198 nlopt_data *data = (nlopt_data *) data_;
200 f = data->f(n, x, NULL, data->f_data);
201 *undefined = isnan(f) || nlopt_isinf(f);
212 # include "l-bfgs-b.h"
223 #include "neldermead.h"
226 #define AUGLAG_ALG(a) ((a) == NLOPT_LN_AUGLAG || \
227 (a) == NLOPT_LN_AUGLAG_EQ || \
228 (a) == NLOPT_LD_AUGLAG || \
229 (a) == NLOPT_LD_AUGLAG_EQ)
231 /*************************************************************************/
233 /* for "hybrid" algorithms that combine global search with some
234 local search algorithm, most of the time we anticipate that people
235 will stick with the default local search algorithm, so we
236 don't add this as a parameter to nlopt_minimize. Instead, the user
237 can call nlopt_{set/get}_hybrid_local_algorithm to get/set the defaults. */
239 /* default local-search algorithms */
240 static nlopt_algorithm local_search_alg_deriv = NLOPT_LD_MMA;
241 static nlopt_algorithm local_search_alg_nonderiv = NLOPT_LN_COBYLA;
243 static int local_search_maxeval = -1; /* no maximum by default */
245 void nlopt_get_local_search_algorithm(nlopt_algorithm *deriv,
246 nlopt_algorithm *nonderiv,
249 *deriv = local_search_alg_deriv;
250 *nonderiv = local_search_alg_nonderiv;
251 *maxeval = local_search_maxeval;
254 void nlopt_set_local_search_algorithm(nlopt_algorithm deriv,
255 nlopt_algorithm nonderiv,
258 local_search_alg_deriv = deriv;
259 local_search_alg_nonderiv = nonderiv;
260 local_search_maxeval = maxeval;
263 /*************************************************************************/
265 /* same as nlopt_minimize_econstrained,
266 but xtol_abs is required to be non-NULL */
267 static nlopt_result nlopt_minimize_(
268 nlopt_algorithm algorithm,
269 int n, nlopt_func f, void *f_data,
270 int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size,
271 int p, nlopt_func h, void *h_data, ptrdiff_t h_datum_size,
272 const double *lb, const double *ub, /* bounds */
273 double *x, /* in: initial guess, out: minimizer */
274 double *minf, /* out: minimum */
275 double minf_max, double ftol_rel, double ftol_abs,
276 double xtol_rel, const double *xtol_abs,
277 double htol_rel, double htol_abs,
278 int maxeval, double maxtime)
284 /* some basic argument checks */
285 if (!minf || !f || n < 0 || m < 0 || p < 0 || (p > 0 && !h)
286 || (m > 0 && !fc)) return NLOPT_INVALID_ARGS;
287 if (n == 0) { /* trivial case: no degrees of freedom */
288 *minf = f(n, x, NULL, f_data);
289 return NLOPT_SUCCESS;
291 else if (n < 0 || !lb || !ub || !x)
292 return NLOPT_INVALID_ARGS;
294 /* nonlinear constraints are only supported with some algorithms */
295 if (m != 0 && algorithm != NLOPT_LD_MMA && algorithm != NLOPT_LN_COBYLA
296 && !AUGLAG_ALG(algorithm))
297 return NLOPT_INVALID_ARGS;
299 /* nonlinear equality constraints (h(x) = 0) only via some algorithms */
300 if (p != 0 && !AUGLAG_ALG(algorithm))
301 return NLOPT_INVALID_ARGS;
308 /* make sure rand generator is inited */
309 if (!nlopt_srand_called)
310 nlopt_srand_time(); /* default is non-deterministic */
312 /* check bound constraints */
313 for (i = 0; i < n; ++i)
314 if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
315 return NLOPT_INVALID_ARGS;
318 stop.minf_max = (isnan(minf_max)
319 || (nlopt_isinf(minf_max) && minf_max < 0))
320 ? -HUGE_VAL : minf_max;
321 stop.ftol_rel = ftol_rel;
322 stop.ftol_abs = ftol_abs;
323 stop.xtol_rel = xtol_rel;
324 stop.xtol_abs = xtol_abs;
325 stop.htol_rel = htol_rel;
326 stop.htol_abs = htol_abs;
328 stop.maxeval = maxeval;
329 stop.maxtime = maxtime;
330 stop.start = nlopt_seconds();
333 case NLOPT_GN_DIRECT:
334 case NLOPT_GN_DIRECT_L:
335 case NLOPT_GN_DIRECT_L_RAND:
336 return cdirect(n, f, f_data, lb, ub, x, minf, &stop, 0.0,
337 (algorithm != NLOPT_GN_DIRECT)
338 + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
339 + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
341 case NLOPT_GN_DIRECT_NOSCAL:
342 case NLOPT_GN_DIRECT_L_NOSCAL:
343 case NLOPT_GN_DIRECT_L_RAND_NOSCAL:
344 return cdirect_unscaled(n, f, f_data, lb, ub, x, minf,
346 (algorithm != NLOPT_GN_DIRECT)
347 + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
348 + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
350 case NLOPT_GN_ORIG_DIRECT:
351 case NLOPT_GN_ORIG_DIRECT_L:
352 switch (direct_optimize(f_direct, &d, n, lb, ub, x, minf,
353 maxeval, -1, 0.0, 0.0,
354 pow(xtol_rel, (double) n), -1.0,
357 algorithm == NLOPT_GN_ORIG_DIRECT
359 : DIRECT_GABLONSKY)) {
360 case DIRECT_INVALID_BOUNDS:
361 case DIRECT_MAXFEVAL_TOOBIG:
362 case DIRECT_INVALID_ARGS:
363 return NLOPT_INVALID_ARGS;
364 case DIRECT_INIT_FAILED:
365 case DIRECT_SAMPLEPOINTS_FAILED:
366 case DIRECT_SAMPLE_FAILED:
367 return NLOPT_FAILURE;
368 case DIRECT_MAXFEVAL_EXCEEDED:
369 case DIRECT_MAXITER_EXCEEDED:
370 return NLOPT_MAXEVAL_REACHED;
371 case DIRECT_GLOBAL_FOUND:
372 return NLOPT_MINF_MAX_REACHED;
374 case DIRECT_SIGMATOL:
375 return NLOPT_XTOL_REACHED;
376 case DIRECT_OUT_OF_MEMORY:
377 return NLOPT_OUT_OF_MEMORY;
382 case NLOPT_GD_STOGO_RAND:
384 if (!stogo_minimize(n, f, f_data, x, minf, lb, ub, &stop,
385 algorithm == NLOPT_GD_STOGO
387 return NLOPT_FAILURE;
390 return NLOPT_FAILURE;
394 /* lacking a free/open-source license, we no longer use
395 Rowan's code, and instead use by "sbplx" re-implementation */
396 case NLOPT_LN_SUBPLEX: {
398 double *scale = (double *) malloc(sizeof(double) * n);
399 if (!scale) return NLOPT_OUT_OF_MEMORY;
400 for (i = 0; i < n; ++i)
401 scale[i] = initial_step(1, lb+i, ub+i, x+i);
402 iret = nlopt_subplex(f_bound, minf, x, n, &d, &stop, scale);
405 case -2: return NLOPT_INVALID_ARGS;
406 case -10: return NLOPT_MAXTIME_REACHED;
407 case -1: return NLOPT_MAXEVAL_REACHED;
408 case 0: return NLOPT_XTOL_REACHED;
409 case 1: return NLOPT_SUCCESS;
410 case 2: return NLOPT_MINF_MAX_REACHED;
411 case 20: return NLOPT_FTOL_REACHED;
412 case -200: return NLOPT_OUT_OF_MEMORY;
413 default: return NLOPT_FAILURE; /* unknown return code */
419 case NLOPT_LN_PRAXIS:
420 return praxis_(0.0, DBL_EPSILON,
421 initial_step(n, lb, ub, x), n, x, f_bound, &d,
425 case NLOPT_LD_LBFGS_NOCEDAL: {
426 int iret, *nbd = (int *) malloc(sizeof(int) * n);
427 if (!nbd) return NLOPT_OUT_OF_MEMORY;
428 for (i = 0; i < n; ++i) {
429 int linf = nlopt_isinf(lb[i]) && lb[i] < 0;
430 int uinf = nlopt_isinf(ub[i]) && ub[i] > 0;
431 nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
433 iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub,
434 MIN(n, 5), 0.0, ftol_rel,
435 xtol_abs ? *xtol_abs : xtol_rel,
440 case -1: return NLOPT_INVALID_ARGS;
441 case -2: default: return NLOPT_FAILURE;
445 *minf = f(n, x, NULL, f_data);
447 case 5: return NLOPT_MAXEVAL_REACHED;
448 case 2: return NLOPT_XTOL_REACHED;
449 case 1: return NLOPT_FTOL_REACHED;
450 default: return NLOPT_SUCCESS;
458 return luksan_plis(n, f, f_data, lb, ub, x, minf, &stop);
462 return luksan_plip(n, f, f_data, lb, ub, x, minf, &stop,
463 algorithm == NLOPT_LD_VAR1 ? 1 : 2);
465 case NLOPT_LD_TNEWTON:
466 case NLOPT_LD_TNEWTON_RESTART:
467 case NLOPT_LD_TNEWTON_PRECOND:
468 case NLOPT_LD_TNEWTON_PRECOND_RESTART:
469 return luksan_pnet(n, f, f_data, lb, ub, x, minf, &stop,
470 1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
471 1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
473 case NLOPT_GN_CRS2_LM:
474 return crs_minimize(n, f, f_data, lb, ub, x, minf, &stop, 0);
478 case NLOPT_GN_MLSL_LDS:
479 case NLOPT_GD_MLSL_LDS:
480 for (i = 0; i < n && stop.xtol_abs[i] <= 0; ++i) ;
481 if (stop.ftol_rel <= 0 && stop.ftol_abs <= 0 &&
482 stop.xtol_rel <= 0 && i == n) {
483 /* it is not sensible to call MLSL without *some*
484 nonzero tolerance for the local search */
485 stop.ftol_rel = 1e-15;
486 stop.xtol_rel = 1e-7;
488 return mlsl_minimize(n, f, f_data, lb, ub, x, minf, &stop,
489 (algorithm == NLOPT_GN_MLSL ||
490 algorithm == NLOPT_GN_MLSL_LDS)
491 ? local_search_alg_nonderiv
492 : local_search_alg_deriv,
493 local_search_maxeval,
494 algorithm >= NLOPT_GN_MLSL_LDS);
497 return mma_minimize(n, f, f_data, m, fc, fc_data, fc_datum_size,
498 lb, ub, x, minf, &stop,
499 local_search_alg_deriv, 1e-12, 100000);
501 case NLOPT_LN_COBYLA:
502 return cobyla_minimize(n, f, f_data,
503 m, fc, fc_data, fc_datum_size,
504 lb, ub, x, minf, &stop,
505 initial_step(n, lb, ub, x));
507 case NLOPT_LN_NEWUOA:
508 return newuoa(n, 2*n+1, x, 0, 0, initial_step(n, lb, ub, x),
509 &stop, minf, f_noderiv, &d);
511 case NLOPT_LN_NEWUOA_BOUND:
512 return newuoa(n, 2*n+1, x, lb, ub, initial_step(n, lb, ub, x),
513 &stop, minf, f_noderiv, &d);
515 case NLOPT_LN_NELDERMEAD:
519 double *scale = (double *) malloc(sizeof(double) * n);
520 if (!scale) return NLOPT_OUT_OF_MEMORY;
521 for (i = 0; i < n; ++i)
522 scale[i] = initial_step(1, lb+i, ub+i, x+i);
523 if (algorithm == NLOPT_LN_NELDERMEAD)
524 ret = nldrmd_minimize(n,f,f_data,lb,ub,x,minf,scale,&stop);
526 ret = sbplx_minimize(n,f,f_data,lb,ub,x,minf,scale,&stop);
531 case NLOPT_LN_AUGLAG:
532 case NLOPT_LN_AUGLAG_EQ:
533 return auglag_minimize(n, f, f_data,
534 m, fc, fc_data, fc_datum_size,
535 p, h, h_data, h_datum_size,
536 lb, ub, x, minf, &stop,
537 local_search_alg_nonderiv,
538 algorithm == NLOPT_LN_AUGLAG_EQ);
540 case NLOPT_LD_AUGLAG:
541 case NLOPT_LD_AUGLAG_EQ:
542 return auglag_minimize(n, f, f_data,
543 m, fc, fc_data, fc_datum_size,
544 p, h, h_data, h_datum_size,
545 lb, ub, x, minf, &stop,
546 local_search_alg_deriv,
547 algorithm == NLOPT_LD_AUGLAG_EQ);
550 return NLOPT_INVALID_ARGS;
553 return NLOPT_SUCCESS;
556 nlopt_result nlopt_minimize_econstrained(
557 nlopt_algorithm algorithm,
558 int n, nlopt_func f, void *f_data,
559 int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size,
560 int p, nlopt_func h, void *h_data, ptrdiff_t h_datum_size,
561 const double *lb, const double *ub, /* bounds */
562 double *x, /* in: initial guess, out: minimizer */
563 double *minf, /* out: minimum */
564 double minf_max, double ftol_rel, double ftol_abs,
565 double xtol_rel, const double *xtol_abs,
566 double htol_rel, double htol_abs,
567 int maxeval, double maxtime)
571 ret = nlopt_minimize_(algorithm, n, f, f_data,
572 m, fc, fc_data, fc_datum_size,
573 p, h, h_data, h_datum_size,
575 x, minf, minf_max, ftol_rel, ftol_abs,
576 xtol_rel, xtol_abs, htol_rel, htol_abs,
580 double *xtol = (double *) malloc(sizeof(double) * n);
581 if (!xtol) return NLOPT_OUT_OF_MEMORY;
582 for (i = 0; i < n; ++i) xtol[i] = -1;
583 ret = nlopt_minimize_(algorithm, n, f, f_data,
584 m, fc, fc_data, fc_datum_size,
585 p, h, h_data, h_datum_size,
587 x, minf, minf_max, ftol_rel, ftol_abs,
588 xtol_rel, xtol, htol_rel, htol_abs,
595 nlopt_result nlopt_minimize_constrained(
596 nlopt_algorithm algorithm,
597 int n, nlopt_func f, void *f_data,
598 int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size,
599 const double *lb, const double *ub, /* bounds */
600 double *x, /* in: initial guess, out: minimizer */
601 double *minf, /* out: minimum */
602 double minf_max, double ftol_rel, double ftol_abs,
603 double xtol_rel, const double *xtol_abs,
604 int maxeval, double maxtime)
606 return nlopt_minimize_econstrained(
607 algorithm, n, f, f_data,
608 m, fc, fc_data, fc_datum_size, 0, NULL, NULL, 0,
609 lb, ub, x, minf, minf_max, ftol_rel, ftol_abs,
610 xtol_rel, xtol_abs, ftol_rel, ftol_abs, maxeval, maxtime);
613 nlopt_result nlopt_minimize(
614 nlopt_algorithm algorithm,
615 int n, nlopt_func f, void *f_data,
616 const double *lb, const double *ub, /* bounds */
617 double *x, /* in: initial guess, out: minimizer */
618 double *minf, /* out: minimum */
619 double minf_max, double ftol_rel, double ftol_abs,
620 double xtol_rel, const double *xtol_abs,
621 int maxeval, double maxtime)
623 return nlopt_minimize_constrained(
624 algorithm, n, f, f_data, 0, NULL, NULL, 0,
625 lb, ub, x, minf, minf_max, ftol_rel, ftol_abs,
626 xtol_rel, xtol_abs, maxeval, maxtime);