1 /* Copyright (c) 2007-2010 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.
27 #include "nlopt-internal.h"
29 /*********************************************************************/
32 static int my_isnan(double x) { return x != x; }
33 # define isnan my_isnan
36 /*********************************************************************/
48 # include "l-bfgs-b.h"
59 #include "neldermead.h"
64 /*********************************************************************/
66 static double f_bound(int n, const double *x, void *data_)
69 nlopt_opt data = (nlopt_opt) data_;
72 /* some methods do not support bound constraints, but support
73 discontinuous objectives so we can just return Inf for invalid x */
74 for (i = 0; i < n; ++i)
75 if (x[i] < data->lb[i] || x[i] > data->ub[i])
78 f = data->f((unsigned) n, x, NULL, data->f_data);
79 return (isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f);
82 static double f_noderiv(int n, const double *x, void *data_)
84 nlopt_opt data = (nlopt_opt) data_;
85 return data->f((unsigned) n, x, NULL, data->f_data);
88 static double f_direct(int n, const double *x, int *undefined, void *data_)
90 nlopt_opt data = (nlopt_opt) data_;
93 f = data->f((unsigned) n, x, NULL, data->f_data);
94 *undefined = isnan(f) || nlopt_isinf(f);
95 for (i = 0; i < data->m && !*undefined; ++i)
96 if (data->fc[i].f((unsigned) n, x, NULL, data->fc[i].f_data) > 0)
101 /*********************************************************************/
103 /* get min(dx) for algorithms requiring a scalar initial step size */
104 static nlopt_result initial_step(nlopt_opt opt, const double *x, double *step)
106 unsigned freedx = 0, i;
110 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
111 return NLOPT_OUT_OF_MEMORY;
115 for (i = 0; i < opt->n; ++i)
116 if (*step > fabs(opt->dx[i]))
117 *step = fabs(opt->dx[i]);
119 if (freedx) { free(opt->dx); opt->dx = NULL; }
120 return NLOPT_SUCCESS;
123 /*********************************************************************/
125 /* return true if [lb,ub] is finite in every dimension (n dimensions) */
126 static int finite_domain(unsigned n, const double *lb, const double *ub)
129 for (i = 0; i < n; ++i)
130 if (nlopt_isinf(ub[i] - lb[i])) return 0;
134 /*********************************************************************/
136 #define POP(defaultpop) (opt->stochastic_population > 0 ? \
137 opt->stochastic_population : \
138 (nlopt_stochastic_population > 0 ? \
139 nlopt_stochastic_population : (defaultpop)))
141 /* unlike nlopt_optimize() below, only handles minimization case */
142 static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf)
144 const double *lb, *ub;
145 nlopt_algorithm algorithm;
146 nlopt_func f; void *f_data;
151 if (!opt || !x || !minf || !opt->f
152 || opt->maximize) return NLOPT_INVALID_ARGS;
154 /* copy a few params to local vars for convenience */
156 ni = (int) n; /* most of the subroutines take "int" arg */
157 lb = opt->lb; ub = opt->ub;
158 algorithm = opt->algorithm;
159 f = opt->f; f_data = opt->f_data;
161 if (n == 0) { /* trivial case: no degrees of freedom */
162 *minf = opt->f(n, x, NULL, opt->f_data);
163 return NLOPT_SUCCESS;
168 /* make sure rand generator is inited */
169 if (!nlopt_srand_called)
170 nlopt_srand_time(); /* default is non-deterministic */
172 /* check bound constraints */
173 for (i = 0; i < n; ++i)
174 if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
175 return NLOPT_INVALID_ARGS;
178 stop.minf_max = opt->stopval;
179 stop.ftol_rel = opt->ftol_rel;
180 stop.ftol_abs = opt->ftol_abs;
181 stop.xtol_rel = opt->xtol_rel;
182 stop.xtol_abs = opt->xtol_abs;
184 stop.maxeval = opt->maxeval;
185 stop.maxtime = opt->maxtime;
186 stop.start = nlopt_seconds();
189 case NLOPT_GN_DIRECT:
190 case NLOPT_GN_DIRECT_L:
191 case NLOPT_GN_DIRECT_L_RAND:
192 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
193 return cdirect(ni, f, f_data,
194 lb, ub, x, minf, &stop, 0.0,
195 (algorithm != NLOPT_GN_DIRECT)
196 + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND
197 ? 2 : (algorithm != NLOPT_GN_DIRECT))
198 + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND
199 ? 1 : (algorithm != NLOPT_GN_DIRECT)));
201 case NLOPT_GN_DIRECT_NOSCAL:
202 case NLOPT_GN_DIRECT_L_NOSCAL:
203 case NLOPT_GN_DIRECT_L_RAND_NOSCAL:
204 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
205 return cdirect_unscaled(ni, f, f_data, lb, ub, x, minf,
207 (algorithm != NLOPT_GN_DIRECT)
208 + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
209 + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
211 case NLOPT_GN_ORIG_DIRECT:
212 case NLOPT_GN_ORIG_DIRECT_L:
213 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
214 switch (direct_optimize(f_direct, opt, ni, lb, ub, x, minf,
215 stop.maxeval, -1, 0.0, 0.0,
216 pow(stop.xtol_rel, (double) n), -1.0,
219 algorithm == NLOPT_GN_ORIG_DIRECT
221 : DIRECT_GABLONSKY)) {
222 case DIRECT_INVALID_BOUNDS:
223 case DIRECT_MAXFEVAL_TOOBIG:
224 case DIRECT_INVALID_ARGS:
225 return NLOPT_INVALID_ARGS;
226 case DIRECT_INIT_FAILED:
227 case DIRECT_SAMPLEPOINTS_FAILED:
228 case DIRECT_SAMPLE_FAILED:
229 return NLOPT_FAILURE;
230 case DIRECT_MAXFEVAL_EXCEEDED:
231 case DIRECT_MAXITER_EXCEEDED:
232 return NLOPT_MAXEVAL_REACHED;
233 case DIRECT_GLOBAL_FOUND:
234 return NLOPT_MINF_MAX_REACHED;
236 case DIRECT_SIGMATOL:
237 return NLOPT_XTOL_REACHED;
238 case DIRECT_OUT_OF_MEMORY:
239 return NLOPT_OUT_OF_MEMORY;
244 case NLOPT_GD_STOGO_RAND:
246 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
247 if (!stogo_minimize(ni, f, f_data, x, minf, lb, ub, &stop,
248 algorithm == NLOPT_GD_STOGO
249 ? 0 : (int) POP(2*n)))
250 return NLOPT_FAILURE;
253 return NLOPT_FAILURE;
257 /* lacking a free/open-source license, we no longer use
258 Rowan's code, and instead use by "sbplx" re-implementation */
259 case NLOPT_LN_SUBPLEX: {
260 int iret, freedx = 0;
263 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
264 return NLOPT_OUT_OF_MEMORY;
266 iret = nlopt_subplex(f_bound, minf, x, n, opt, &stop, opt->dx);
267 if (freedx) { free(opt->dx); opt->dx = NULL; }
269 case -2: return NLOPT_INVALID_ARGS;
270 case -10: return NLOPT_MAXTIME_REACHED;
271 case -1: return NLOPT_MAXEVAL_REACHED;
272 case 0: return NLOPT_XTOL_REACHED;
273 case 1: return NLOPT_SUCCESS;
274 case 2: return NLOPT_MINF_MAX_REACHED;
275 case 20: return NLOPT_FTOL_REACHED;
276 case -200: return NLOPT_OUT_OF_MEMORY;
277 default: return NLOPT_FAILURE; /* unknown return code */
283 case NLOPT_LN_PRAXIS: {
285 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
286 return NLOPT_OUT_OF_MEMORY;
287 return praxis_(0.0, DBL_EPSILON,
288 step, ni, x, f_bound, opt, &stop, minf);
292 case NLOPT_LD_LBFGS_NOCEDAL: {
293 int iret, *nbd = (int *) malloc(sizeof(int) * n);
294 if (!nbd) return NLOPT_OUT_OF_MEMORY;
295 for (i = 0; i < n; ++i) {
296 int linf = nlopt_isinf(lb[i]) && lb[i] < 0;
297 int uinf = nlopt_isinf(ub[i]) && ub[i] > 0;
298 nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
300 iret = lbfgsb_minimize(ni, f, f_data, x, nbd, lb, ub,
301 MIN(ni, 5), 0.0, stop.ftol_rel,
302 stop.xtol_abs[0] > 0 ? stop.xtol_abs[0]
308 case -1: return NLOPT_INVALID_ARGS;
309 case -2: default: return NLOPT_FAILURE;
313 *minf = f(n, x, NULL, f_data);
315 case 5: return NLOPT_MAXEVAL_REACHED;
316 case 2: return NLOPT_XTOL_REACHED;
317 case 1: return NLOPT_FTOL_REACHED;
318 default: return NLOPT_SUCCESS;
326 return luksan_plis(ni, f, f_data, lb, ub, x, minf, &stop);
330 return luksan_plip(ni, f, f_data, lb, ub, x, minf, &stop,
331 algorithm == NLOPT_LD_VAR1 ? 1 : 2);
333 case NLOPT_LD_TNEWTON:
334 case NLOPT_LD_TNEWTON_RESTART:
335 case NLOPT_LD_TNEWTON_PRECOND:
336 case NLOPT_LD_TNEWTON_PRECOND_RESTART:
337 return luksan_pnet(ni, f, f_data, lb, ub, x, minf, &stop,
338 1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
339 1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
341 case NLOPT_GN_CRS2_LM:
342 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
343 return crs_minimize(ni, f, f_data, lb, ub, x, minf, &stop,
348 case NLOPT_GN_MLSL_LDS:
349 case NLOPT_GD_MLSL_LDS: {
350 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
351 nlopt_opt local_opt = opt->local_opt;
353 if (!local_opt) { /* default */
354 nlopt_algorithm local_alg = (algorithm == NLOPT_GN_MLSL ||
355 algorithm == NLOPT_GN_MLSL_LDS)
356 ? nlopt_local_search_alg_nonderiv
357 : nlopt_local_search_alg_deriv;
358 /* don't call MLSL recursively! */
359 if (local_alg >= NLOPT_GN_MLSL
360 && local_alg <= NLOPT_GD_MLSL_LDS)
361 local_alg = (algorithm == NLOPT_GN_MLSL ||
362 algorithm == NLOPT_GN_MLSL_LDS)
363 ? NLOPT_LN_COBYLA : NLOPT_LD_MMA;
364 local_opt = nlopt_create(local_alg, n);
365 if (!local_opt) return NLOPT_FAILURE;
366 nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
367 nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
368 nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
369 nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
370 nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
371 nlopt_set_initial_step(local_opt, opt->dx);
373 for (i = 0; i < n && stop.xtol_abs[i] > 0; ++i) ;
374 if (local_opt->ftol_rel <= 0 && local_opt->ftol_abs <= 0 &&
375 local_opt->xtol_rel <= 0 && i < n) {
376 /* it is not sensible to call MLSL without *some*
377 nonzero tolerance for the local search */
378 nlopt_set_ftol_rel(local_opt, 1e-15);
379 nlopt_set_xtol_rel(local_opt, 1e-7);
381 ret = mlsl_minimize(ni, f, f_data, lb, ub, x, minf, &stop,
382 local_opt, (int) POP(0),
383 algorithm >= NLOPT_GN_MLSL_LDS);
384 if (!opt->local_opt) nlopt_destroy(local_opt);
391 #define LO(param, def) (opt->local_opt ? opt->local_opt->param : (def))
392 dual_opt = nlopt_create(LO(algorithm,
393 nlopt_local_search_alg_deriv),
395 if (!dual_opt) return NLOPT_FAILURE;
396 nlopt_set_ftol_rel(dual_opt, LO(ftol_rel, 1e-12));
397 nlopt_set_ftol_abs(dual_opt, LO(ftol_abs, 0.0));
398 nlopt_set_maxeval(dual_opt, LO(maxeval, 100000));
401 ret = mma_minimize(ni, f, f_data, (int) (opt->m), opt->fc,
402 lb, ub, x, minf, &stop, dual_opt);
403 nlopt_destroy(dual_opt);
407 case NLOPT_LN_COBYLA: {
409 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
410 return NLOPT_OUT_OF_MEMORY;
411 return cobyla_minimize(ni, f, f_data,
413 lb, ub, x, minf, &stop,
417 case NLOPT_LN_NEWUOA: {
419 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
420 return NLOPT_OUT_OF_MEMORY;
421 return newuoa(ni, 2*n+1, x, 0, 0, step,
422 &stop, minf, f_noderiv, opt);
425 case NLOPT_LN_NEWUOA_BOUND: {
427 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
428 return NLOPT_OUT_OF_MEMORY;
429 return newuoa(ni, 2*n+1, x, lb, ub, step,
430 &stop, minf, f_noderiv, opt);
433 case NLOPT_LN_BOBYQA: {
435 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
436 return NLOPT_OUT_OF_MEMORY;
437 return bobyqa(ni, 2*n+1, x, lb, ub, step,
438 &stop, minf, f_noderiv, opt);
441 case NLOPT_LN_NELDERMEAD:
448 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
449 return NLOPT_OUT_OF_MEMORY;
451 if (algorithm == NLOPT_LN_NELDERMEAD)
452 ret= nldrmd_minimize(ni,f,f_data,lb,ub,x,minf,opt->dx,&stop);
454 ret= sbplx_minimize(ni,f,f_data,lb,ub,x,minf,opt->dx,&stop);
455 if (freedx) { free(opt->dx); opt->dx = NULL; }
459 case NLOPT_LN_AUGLAG:
460 case NLOPT_LN_AUGLAG_EQ:
461 case NLOPT_LD_AUGLAG:
462 case NLOPT_LD_AUGLAG_EQ: {
463 nlopt_opt local_opt = opt->local_opt;
465 if (!local_opt) { /* default */
466 local_opt = nlopt_create(
467 algorithm == NLOPT_LN_AUGLAG ||
468 algorithm == NLOPT_LN_AUGLAG_EQ
469 ? nlopt_local_search_alg_nonderiv
470 : nlopt_local_search_alg_deriv, n);
471 if (!local_opt) return NLOPT_FAILURE;
472 nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
473 nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
474 nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
475 nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
476 nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
477 nlopt_set_initial_step(local_opt, opt->dx);
479 ret = auglag_minimize(ni, f, f_data,
482 lb, ub, x, minf, &stop,
484 algorithm == NLOPT_LN_AUGLAG_EQ
485 || algorithm == NLOPT_LD_AUGLAG_EQ);
486 if (!opt->local_opt) nlopt_destroy(local_opt);
491 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
492 return isres_minimize(ni, f, f_data,
493 (int) (opt->m), opt->fc,
494 (int) (opt->p), opt->h,
495 lb, ub, x, minf, &stop,
499 return NLOPT_INVALID_ARGS;
502 return NLOPT_SUCCESS; /* never reached */
505 /*********************************************************************/
512 /* wrapper for maximizing: just flip the sign of f and grad */
513 static double f_max(unsigned n, const double *x, double *grad, void *data)
515 f_max_data *d = (f_max_data *) data;
516 double val = d->f(n, x, grad, d->f_data);
519 for (i = 0; i < n; ++i)
525 nlopt_result nlopt_optimize(nlopt_opt opt, double *x, double *opt_f)
527 nlopt_func f; void *f_data;
532 if (!opt || !opt_f || !opt->f) return NLOPT_INVALID_ARGS;
533 f = opt->f; f_data = opt->f_data;
535 /* for maximizing, just minimize the f_max wrapper, which
536 flips the sign of everything */
537 if ((maximize = opt->maximize)) {
538 fmd.f = f; fmd.f_data = f_data;
539 opt->f = f_max; opt->f_data = &fmd;
540 opt->stopval = -opt->stopval;
544 ret = nlopt_optimize_(opt, x, opt_f);
546 if (maximize) { /* restore original signs */
547 opt->maximize = maximize;
548 opt->stopval = -opt->stopval;
549 opt->f = f; opt->f_data = f_data;
556 /*********************************************************************/
558 nlopt_result nlopt_optimize_limited(nlopt_opt opt, double *x, double *minf,
559 int maxeval, double maxtime)
565 if (!opt) return NLOPT_INVALID_ARGS;
567 save_maxeval = nlopt_get_maxeval(opt);
568 save_maxtime = nlopt_get_maxtime(opt);
570 /* override opt limits if maxeval and/or maxtime are more stringent */
571 if (save_maxeval <= 0 || (maxeval > 0 && maxeval < save_maxeval))
572 nlopt_set_maxeval(opt, maxeval);
573 if (save_maxtime <= 0 || (maxtime > 0 && maxtime < save_maxtime))
574 nlopt_set_maxtime(opt, maxtime);
576 ret = nlopt_optimize(opt, x, minf);
578 nlopt_set_maxeval(opt, save_maxeval);
579 nlopt_set_maxtime(opt, save_maxtime);
584 /*********************************************************************/