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"
65 /*********************************************************************/
67 static double f_bound(int n, const double *x, void *data_)
70 nlopt_opt data = (nlopt_opt) data_;
73 /* some methods do not support bound constraints, but support
74 discontinuous objectives so we can just return Inf for invalid x */
75 for (i = 0; i < n; ++i)
76 if (x[i] < data->lb[i] || x[i] > data->ub[i])
79 f = data->f((unsigned) n, x, NULL, data->f_data);
80 return (isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f);
83 static double f_noderiv(int n, const double *x, void *data_)
85 nlopt_opt data = (nlopt_opt) data_;
86 return data->f((unsigned) n, x, NULL, data->f_data);
89 static double f_direct(int n, const double *x, int *undefined, void *data_)
91 nlopt_opt data = (nlopt_opt) data_;
94 f = data->f((unsigned) n, x, NULL, data->f_data);
95 *undefined = isnan(f) || nlopt_isinf(f);
96 if (nlopt_get_force_stop(data)) return f;
97 for (i = 0; i < data->m && !*undefined; ++i) {
98 nlopt_eval_constraint(data->work, NULL, data->fc+i, (unsigned) n, x);
99 if (nlopt_get_force_stop(data)) return f;
100 for (j = 0; j < data->fc[i].m; ++j)
101 if (data->work[j] > 0)
107 /*********************************************************************/
109 /* get min(dx) for algorithms requiring a scalar initial step size */
110 static nlopt_result initial_step(nlopt_opt opt, const double *x, double *step)
112 unsigned freedx = 0, i;
116 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
117 return NLOPT_OUT_OF_MEMORY;
121 for (i = 0; i < opt->n; ++i)
122 if (*step > fabs(opt->dx[i]))
123 *step = fabs(opt->dx[i]);
125 if (freedx) { free(opt->dx); opt->dx = NULL; }
126 return NLOPT_SUCCESS;
129 /*********************************************************************/
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)
135 for (i = 0; i < n; ++i)
136 if (nlopt_isinf(ub[i] - lb[i])) return 0;
140 /*********************************************************************/
142 #define POP(defaultpop) (opt->stochastic_population > 0 ? \
143 opt->stochastic_population : \
144 (nlopt_stochastic_population > 0 ? \
145 nlopt_stochastic_population : (defaultpop)))
147 /* unlike nlopt_optimize() below, only handles minimization case */
148 static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf)
150 const double *lb, *ub;
151 nlopt_algorithm algorithm;
152 nlopt_func f; void *f_data;
157 if (!opt || !x || !minf || !opt->f
158 || opt->maximize) return NLOPT_INVALID_ARGS;
160 /* reset stopping flag */
161 nlopt_set_force_stop(opt, 0);
162 opt->force_stop_child = NULL;
164 /* copy a few params to local vars for convenience */
166 ni = (int) n; /* most of the subroutines take "int" arg */
167 lb = opt->lb; ub = opt->ub;
168 algorithm = opt->algorithm;
169 f = opt->f; f_data = opt->f_data;
171 if (n == 0) { /* trivial case: no degrees of freedom */
172 *minf = opt->f(n, x, NULL, opt->f_data);
173 return NLOPT_SUCCESS;
178 /* make sure rand generator is inited */
179 nlopt_srand_time_default(); /* default is non-deterministic */
181 /* check bound constraints */
182 for (i = 0; i < n; ++i)
183 if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
184 return NLOPT_INVALID_ARGS;
187 stop.minf_max = opt->stopval;
188 stop.ftol_rel = opt->ftol_rel;
189 stop.ftol_abs = opt->ftol_abs;
190 stop.xtol_rel = opt->xtol_rel;
191 stop.xtol_abs = opt->xtol_abs;
193 stop.maxeval = opt->maxeval;
194 stop.maxtime = opt->maxtime;
195 stop.start = nlopt_seconds();
196 stop.force_stop = &(opt->force_stop);
199 case NLOPT_GN_DIRECT:
200 case NLOPT_GN_DIRECT_L:
201 case NLOPT_GN_DIRECT_L_RAND:
202 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
203 return cdirect(ni, f, f_data,
204 lb, ub, x, minf, &stop, 0.0,
205 (algorithm != NLOPT_GN_DIRECT)
206 + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND
207 ? 2 : (algorithm != NLOPT_GN_DIRECT))
208 + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND
209 ? 1 : (algorithm != NLOPT_GN_DIRECT)));
211 case NLOPT_GN_DIRECT_NOSCAL:
212 case NLOPT_GN_DIRECT_L_NOSCAL:
213 case NLOPT_GN_DIRECT_L_RAND_NOSCAL:
214 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
215 return cdirect_unscaled(ni, f, f_data, lb, ub, x, minf,
217 (algorithm != NLOPT_GN_DIRECT)
218 + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
219 + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
221 case NLOPT_GN_ORIG_DIRECT:
222 case NLOPT_GN_ORIG_DIRECT_L: {
223 direct_return_code dret;
224 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
225 opt->work = (double*) malloc(sizeof(double) *
226 nlopt_max_constraint_dim(opt->m,
228 if (!opt->work) return NLOPT_OUT_OF_MEMORY;
229 dret = direct_optimize(f_direct, opt, ni, lb, ub, x, minf,
231 stop.start, stop.maxtime,
233 pow(stop.xtol_rel, (double) n), -1.0,
237 algorithm == NLOPT_GN_ORIG_DIRECT
240 free(opt->work); opt->work = NULL;
242 case DIRECT_INVALID_BOUNDS:
243 case DIRECT_MAXFEVAL_TOOBIG:
244 case DIRECT_INVALID_ARGS:
245 return NLOPT_INVALID_ARGS;
246 case DIRECT_INIT_FAILED:
247 case DIRECT_SAMPLEPOINTS_FAILED:
248 case DIRECT_SAMPLE_FAILED:
249 return NLOPT_FAILURE;
250 case DIRECT_MAXFEVAL_EXCEEDED:
251 case DIRECT_MAXITER_EXCEEDED:
252 return NLOPT_MAXEVAL_REACHED;
253 case DIRECT_MAXTIME_EXCEEDED:
254 return NLOPT_MAXTIME_REACHED;
255 case DIRECT_GLOBAL_FOUND:
256 return NLOPT_MINF_MAX_REACHED;
258 case DIRECT_SIGMATOL:
259 return NLOPT_XTOL_REACHED;
260 case DIRECT_OUT_OF_MEMORY:
261 return NLOPT_OUT_OF_MEMORY;
262 case DIRECT_FORCED_STOP:
263 return NLOPT_FORCED_STOP;
269 case NLOPT_GD_STOGO_RAND:
271 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
272 if (!stogo_minimize(ni, f, f_data, x, minf, lb, ub, &stop,
273 algorithm == NLOPT_GD_STOGO
274 ? 0 : (int) POP(2*n)))
275 return NLOPT_FAILURE;
278 return NLOPT_INVALID_ARGS;
282 /* lacking a free/open-source license, we no longer use
283 Rowan's code, and instead use by "sbplx" re-implementation */
284 case NLOPT_LN_SUBPLEX: {
285 int iret, freedx = 0;
288 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
289 return NLOPT_OUT_OF_MEMORY;
291 iret = nlopt_subplex(f_bound, minf, x, n, opt, &stop, opt->dx);
292 if (freedx) { free(opt->dx); opt->dx = NULL; }
294 case -2: return NLOPT_INVALID_ARGS;
295 case -20: return NLOPT_FORCED_STOP;
296 case -10: return NLOPT_MAXTIME_REACHED;
297 case -1: return NLOPT_MAXEVAL_REACHED;
298 case 0: return NLOPT_XTOL_REACHED;
299 case 1: return NLOPT_SUCCESS;
300 case 2: return NLOPT_MINF_MAX_REACHED;
301 case 20: return NLOPT_FTOL_REACHED;
302 case -200: return NLOPT_OUT_OF_MEMORY;
303 default: return NLOPT_FAILURE; /* unknown return code */
309 case NLOPT_LN_PRAXIS: {
311 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
312 return NLOPT_OUT_OF_MEMORY;
313 return praxis_(0.0, DBL_EPSILON,
314 step, ni, x, f_bound, opt, &stop, minf);
318 case NLOPT_LD_LBFGS_NOCEDAL: {
319 int iret, *nbd = (int *) malloc(sizeof(int) * n);
320 if (!nbd) return NLOPT_OUT_OF_MEMORY;
321 for (i = 0; i < n; ++i) {
322 int linf = nlopt_isinf(lb[i]) && lb[i] < 0;
323 int uinf = nlopt_isinf(ub[i]) && ub[i] > 0;
324 nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
326 iret = lbfgsb_minimize(ni, f, f_data, x, nbd, lb, ub,
327 ni < 5 ? ni : 5, 0.0, stop.ftol_rel,
328 stop.xtol_abs[0] > 0 ? stop.xtol_abs[0]
334 case -1: return NLOPT_INVALID_ARGS;
335 case -2: default: return NLOPT_FAILURE;
339 *minf = f(n, x, NULL, f_data);
341 case 5: return NLOPT_MAXEVAL_REACHED;
342 case 2: return NLOPT_XTOL_REACHED;
343 case 1: return NLOPT_FTOL_REACHED;
344 default: return NLOPT_SUCCESS;
352 return luksan_plis(ni, f, f_data, lb, ub, x, minf,
353 &stop, opt->subspace_dim);
357 return luksan_plip(ni, f, f_data, lb, ub, x, minf,
358 &stop, opt->subspace_dim,
359 algorithm == NLOPT_LD_VAR1 ? 1 : 2);
361 case NLOPT_LD_TNEWTON:
362 case NLOPT_LD_TNEWTON_RESTART:
363 case NLOPT_LD_TNEWTON_PRECOND:
364 case NLOPT_LD_TNEWTON_PRECOND_RESTART:
365 return luksan_pnet(ni, f, f_data, lb, ub, x, minf,
366 &stop, opt->subspace_dim,
367 1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
368 1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
370 case NLOPT_GN_CRS2_LM:
371 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
372 return crs_minimize(ni, f, f_data, lb, ub, x, minf, &stop,
376 case NLOPT_G_MLSL_LDS:
379 case NLOPT_GN_MLSL_LDS:
380 case NLOPT_GD_MLSL_LDS: {
381 nlopt_opt local_opt = opt->local_opt;
383 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
384 if (!local_opt && (algorithm == NLOPT_G_MLSL
385 || algorithm == NLOPT_G_MLSL_LDS))
386 return NLOPT_INVALID_ARGS;
387 if (!local_opt) { /* default */
388 nlopt_algorithm local_alg = (algorithm == NLOPT_GN_MLSL ||
389 algorithm == NLOPT_GN_MLSL_LDS)
390 ? nlopt_local_search_alg_nonderiv
391 : nlopt_local_search_alg_deriv;
392 /* don't call MLSL recursively! */
393 if (local_alg >= NLOPT_GN_MLSL
394 && local_alg <= NLOPT_GD_MLSL_LDS)
395 local_alg = (algorithm == NLOPT_GN_MLSL ||
396 algorithm == NLOPT_GN_MLSL_LDS)
397 ? NLOPT_LN_COBYLA : NLOPT_LD_MMA;
398 local_opt = nlopt_create(local_alg, n);
399 if (!local_opt) return NLOPT_FAILURE;
400 nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
401 nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
402 nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
403 nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
404 nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
406 if (opt->dx) nlopt_set_initial_step(local_opt, opt->dx);
407 for (i = 0; i < n && stop.xtol_abs[i] > 0; ++i) ;
408 if (local_opt->ftol_rel <= 0 && local_opt->ftol_abs <= 0 &&
409 local_opt->xtol_rel <= 0 && i < n) {
410 /* it is not sensible to call MLSL without *some*
411 nonzero tolerance for the local search */
412 nlopt_set_ftol_rel(local_opt, 1e-15);
413 nlopt_set_xtol_rel(local_opt, 1e-7);
415 opt->force_stop_child = local_opt;
416 ret = mlsl_minimize(ni, f, f_data, lb, ub, x, minf, &stop,
417 local_opt, (int) POP(0),
418 algorithm >= NLOPT_GN_MLSL_LDS &&
419 algorithm != NLOPT_G_MLSL);
420 opt->force_stop_child = NULL;
421 if (!opt->local_opt) nlopt_destroy(local_opt);
428 #define LO(param, def) (opt->local_opt ? opt->local_opt->param : (def))
429 dual_opt = nlopt_create(LO(algorithm,
430 nlopt_local_search_alg_deriv),
431 nlopt_count_constraints(opt->m,
433 if (!dual_opt) return NLOPT_FAILURE;
434 nlopt_set_ftol_rel(dual_opt, LO(ftol_rel, 1e-12));
435 nlopt_set_ftol_abs(dual_opt, LO(ftol_abs, 0.0));
436 nlopt_set_maxeval(dual_opt, LO(maxeval, 100000));
439 ret = mma_minimize(n, f, f_data, opt->m, opt->fc,
440 lb, ub, x, minf, &stop, dual_opt);
441 nlopt_destroy(dual_opt);
445 case NLOPT_LN_COBYLA: {
450 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
451 return NLOPT_OUT_OF_MEMORY;
453 return cobyla_minimize(n, f, f_data,
456 lb, ub, x, minf, &stop,
458 if (freedx) { free(opt->dx); opt->dx = NULL; }
462 case NLOPT_LN_NEWUOA: {
464 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
465 return NLOPT_OUT_OF_MEMORY;
466 return newuoa(ni, 2*n+1, x, 0, 0, step,
467 &stop, minf, f_noderiv, opt);
470 case NLOPT_LN_NEWUOA_BOUND: {
472 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
473 return NLOPT_OUT_OF_MEMORY;
474 return newuoa(ni, 2*n+1, x, lb, ub, step,
475 &stop, minf, f_noderiv, opt);
478 case NLOPT_LN_BOBYQA: {
483 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
484 return NLOPT_OUT_OF_MEMORY;
486 ret = bobyqa(ni, 2*n+1, x, lb, ub, opt->dx,
487 &stop, minf, opt->f, opt->f_data);
488 if (freedx) { free(opt->dx); opt->dx = NULL; }
492 case NLOPT_LN_NELDERMEAD:
499 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
500 return NLOPT_OUT_OF_MEMORY;
502 if (algorithm == NLOPT_LN_NELDERMEAD)
503 ret= nldrmd_minimize(ni,f,f_data,lb,ub,x,minf,opt->dx,&stop);
505 ret= sbplx_minimize(ni,f,f_data,lb,ub,x,minf,opt->dx,&stop);
506 if (freedx) { free(opt->dx); opt->dx = NULL; }
511 case NLOPT_AUGLAG_EQ:
512 case NLOPT_LN_AUGLAG:
513 case NLOPT_LN_AUGLAG_EQ:
514 case NLOPT_LD_AUGLAG:
515 case NLOPT_LD_AUGLAG_EQ: {
516 nlopt_opt local_opt = opt->local_opt;
518 if ((algorithm == NLOPT_AUGLAG || algorithm == NLOPT_AUGLAG_EQ)
520 return NLOPT_INVALID_ARGS;
521 if (!local_opt) { /* default */
522 local_opt = nlopt_create(
523 algorithm == NLOPT_LN_AUGLAG ||
524 algorithm == NLOPT_LN_AUGLAG_EQ
525 ? nlopt_local_search_alg_nonderiv
526 : nlopt_local_search_alg_deriv, n);
527 if (!local_opt) return NLOPT_FAILURE;
528 nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
529 nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
530 nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
531 nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
532 nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
534 if (opt->dx) nlopt_set_initial_step(local_opt, opt->dx);
535 opt->force_stop_child = local_opt;
536 ret = auglag_minimize(ni, f, f_data,
539 lb, ub, x, minf, &stop,
541 algorithm == NLOPT_AUGLAG_EQ
542 || algorithm == NLOPT_LN_AUGLAG_EQ
543 || algorithm == NLOPT_LD_AUGLAG_EQ);
544 opt->force_stop_child = NULL;
545 if (!opt->local_opt) nlopt_destroy(local_opt);
550 if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
551 return isres_minimize(ni, f, f_data,
552 (int) (opt->m), opt->fc,
553 (int) (opt->p), opt->h,
554 lb, ub, x, minf, &stop,
558 return nlopt_slsqp(n, f, f_data,
561 lb, ub, x, minf, &stop);
564 return NLOPT_INVALID_ARGS;
567 return NLOPT_SUCCESS; /* never reached */
570 /*********************************************************************/
577 /* wrapper for maximizing: just flip the sign of f and grad */
578 static double f_max(unsigned n, const double *x, double *grad, void *data)
580 f_max_data *d = (f_max_data *) data;
581 double val = d->f(n, x, grad, d->f_data);
584 for (i = 0; i < n; ++i)
591 NLOPT_STDCALL nlopt_optimize(nlopt_opt opt, double *x, double *opt_f)
593 nlopt_func f; void *f_data;
598 if (!opt || !opt_f || !opt->f) return NLOPT_INVALID_ARGS;
599 f = opt->f; f_data = opt->f_data;
601 /* for maximizing, just minimize the f_max wrapper, which
602 flips the sign of everything */
603 if ((maximize = opt->maximize)) {
604 fmd.f = f; fmd.f_data = f_data;
605 opt->f = f_max; opt->f_data = &fmd;
606 opt->stopval = -opt->stopval;
610 ret = nlopt_optimize_(opt, x, opt_f);
612 if (maximize) { /* restore original signs */
613 opt->maximize = maximize;
614 opt->stopval = -opt->stopval;
615 opt->f = f; opt->f_data = f_data;
622 /*********************************************************************/
624 nlopt_result nlopt_optimize_limited(nlopt_opt opt, double *x, double *minf,
625 int maxeval, double maxtime)
631 if (!opt) return NLOPT_INVALID_ARGS;
633 save_maxeval = nlopt_get_maxeval(opt);
634 save_maxtime = nlopt_get_maxtime(opt);
636 /* override opt limits if maxeval and/or maxtime are more stringent */
637 if (save_maxeval <= 0 || (maxeval > 0 && maxeval < save_maxeval))
638 nlopt_set_maxeval(opt, maxeval);
639 if (save_maxtime <= 0 || (maxtime > 0 && maxtime < save_maxtime))
640 nlopt_set_maxtime(opt, maxtime);
642 ret = nlopt_optimize(opt, x, minf);
644 nlopt_set_maxeval(opt, save_maxeval);
645 nlopt_set_maxtime(opt, save_maxtime);
650 /*********************************************************************/