1 /* Copyright (c) 2007-2009 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 /*********************************************************************/
69 const double *lb, *ub;
74 static double f_bound(int n, const double *x, void *data_)
77 nlopt_data *data = (nlopt_data *) data_;
80 /* some methods do not support bound constraints, but support
81 discontinuous objectives so we can just return Inf for invalid x */
82 for (i = 0; i < n; ++i)
83 if (x[i] < data->lb[i] || x[i] > data->ub[i])
86 f = data->f(n, x, NULL, data->f_data);
87 return (isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f);
90 static double f_noderiv(int n, const double *x, void *data_)
92 nlopt_data *data = (nlopt_data *) data_;
93 return data->f(n, x, NULL, data->f_data);
96 static double f_direct(int n, const double *x, int *undefined, void *data_)
98 nlopt_data *data = (nlopt_data *) data_;
100 f = data->f(n, x, NULL, data->f_data);
101 *undefined = isnan(f) || nlopt_isinf(f);
105 /*********************************************************************/
107 /* get min(dx) for algorithms requiring a scalar initial step size */
108 static nlopt_result initial_step(nlopt_opt opt, const double *x, double *step)
114 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
115 return NLOPT_OUT_OF_MEMORY;
119 for (i = 0; i < opt->n; ++i)
120 if (*step > fabs(opt->dx[i]))
121 *step = fabs(opt->dx[i]);
123 if (freedx) { free(opt->dx); opt->dx = NULL; }
124 return NLOPT_SUCCESS;
127 /*********************************************************************/
129 #define POP(defaultpop) (opt->stochastic_population > 0 ? \
130 opt->stochastic_population : \
131 (nlopt_stochastic_population > 0 ? \
132 nlopt_stochastic_population : (defaultpop)))
134 nlopt_result nlopt_optimize(nlopt_opt opt, double *x, double *minf)
136 const double *lb, *ub;
137 nlopt_algorithm algorithm;
138 nlopt_func f; void *f_data;
143 if (!opt || !x || !minf || !opt->f) return NLOPT_INVALID_ARGS;
145 /* copy a few params to local vars for convenience */
147 lb = opt->lb; ub = opt->ub;
148 algorithm = opt->algorithm;
149 f = opt->f; f_data = opt->f_data;
151 if (n == 0) { /* trivial case: no degrees of freedom */
152 *minf = f(n, x, NULL, f_data);
153 return NLOPT_SUCCESS;
163 /* make sure rand generator is inited */
164 if (!nlopt_srand_called)
165 nlopt_srand_time(); /* default is non-deterministic */
167 /* check bound constraints */
168 for (i = 0; i < n; ++i)
169 if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
170 return NLOPT_INVALID_ARGS;
173 stop.minf_max = opt->minf_max;
174 stop.ftol_rel = opt->ftol_rel;
175 stop.ftol_abs = opt->ftol_abs;
176 stop.xtol_rel = opt->xtol_rel;
177 stop.xtol_abs = opt->xtol_abs;
179 stop.maxeval = opt->maxeval;
180 stop.maxtime = opt->maxtime;
181 stop.start = nlopt_seconds();
184 case NLOPT_GN_DIRECT:
185 case NLOPT_GN_DIRECT_L:
186 case NLOPT_GN_DIRECT_L_RAND:
187 return cdirect(n, f, f_data,
188 lb, ub, x, minf, &stop, 0.0,
189 (algorithm != NLOPT_GN_DIRECT)
190 + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND
191 ? 2 : (algorithm != NLOPT_GN_DIRECT))
192 + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND
193 ? 1 : (algorithm != NLOPT_GN_DIRECT)));
195 case NLOPT_GN_DIRECT_NOSCAL:
196 case NLOPT_GN_DIRECT_L_NOSCAL:
197 case NLOPT_GN_DIRECT_L_RAND_NOSCAL:
198 return cdirect_unscaled(n, f, f_data, lb, ub, x, minf,
200 (algorithm != NLOPT_GN_DIRECT)
201 + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
202 + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
204 case NLOPT_GN_ORIG_DIRECT:
205 case NLOPT_GN_ORIG_DIRECT_L:
206 switch (direct_optimize(f_direct, &d, n, lb, ub, x, minf,
207 stop.maxeval, -1, 0.0, 0.0,
208 pow(stop.xtol_rel, (double) n), -1.0,
211 algorithm == NLOPT_GN_ORIG_DIRECT
213 : DIRECT_GABLONSKY)) {
214 case DIRECT_INVALID_BOUNDS:
215 case DIRECT_MAXFEVAL_TOOBIG:
216 case DIRECT_INVALID_ARGS:
217 return NLOPT_INVALID_ARGS;
218 case DIRECT_INIT_FAILED:
219 case DIRECT_SAMPLEPOINTS_FAILED:
220 case DIRECT_SAMPLE_FAILED:
221 return NLOPT_FAILURE;
222 case DIRECT_MAXFEVAL_EXCEEDED:
223 case DIRECT_MAXITER_EXCEEDED:
224 return NLOPT_MAXEVAL_REACHED;
225 case DIRECT_GLOBAL_FOUND:
226 return NLOPT_MINF_MAX_REACHED;
228 case DIRECT_SIGMATOL:
229 return NLOPT_XTOL_REACHED;
230 case DIRECT_OUT_OF_MEMORY:
231 return NLOPT_OUT_OF_MEMORY;
236 case NLOPT_GD_STOGO_RAND:
238 if (!stogo_minimize(n, f, f_data, x, minf, lb, ub, &stop,
239 algorithm == NLOPT_GD_STOGO
241 return NLOPT_FAILURE;
244 return NLOPT_FAILURE;
248 /* lacking a free/open-source license, we no longer use
249 Rowan's code, and instead use by "sbplx" re-implementation */
250 case NLOPT_LN_SUBPLEX: {
251 int iret, freedx = 0;
254 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
255 return NLOPT_OUT_OF_MEMORY;
257 iret = nlopt_subplex(f_bound, minf, x, n, &d, &stop, opt->dx);
258 if (freedx) { free(opt->dx); opt->dx = NULL; }
260 case -2: return NLOPT_INVALID_ARGS;
261 case -10: return NLOPT_MAXTIME_REACHED;
262 case -1: return NLOPT_MAXEVAL_REACHED;
263 case 0: return NLOPT_XTOL_REACHED;
264 case 1: return NLOPT_SUCCESS;
265 case 2: return NLOPT_MINF_MAX_REACHED;
266 case 20: return NLOPT_FTOL_REACHED;
267 case -200: return NLOPT_OUT_OF_MEMORY;
268 default: return NLOPT_FAILURE; /* unknown return code */
274 case NLOPT_LN_PRAXIS: {
276 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
277 return NLOPT_OUT_OF_MEMORY;
278 return praxis_(0.0, DBL_EPSILON,
279 step, n, x, f_bound, &d, &stop, minf);
283 case NLOPT_LD_LBFGS_NOCEDAL: {
284 int iret, *nbd = (int *) malloc(sizeof(int) * n);
285 if (!nbd) return NLOPT_OUT_OF_MEMORY;
286 for (i = 0; i < n; ++i) {
287 int linf = nlopt_isinf(lb[i]) && lb[i] < 0;
288 int uinf = nlopt_isinf(ub[i]) && ub[i] > 0;
289 nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
291 iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub,
292 MIN(n, 5), 0.0, stop.ftol_rel,
293 stop.xtol_abs[0] > 0 ? stop.xtol_abs[0]
299 case -1: return NLOPT_INVALID_ARGS;
300 case -2: default: return NLOPT_FAILURE;
304 *minf = f(n, x, NULL, f_data);
306 case 5: return NLOPT_MAXEVAL_REACHED;
307 case 2: return NLOPT_XTOL_REACHED;
308 case 1: return NLOPT_FTOL_REACHED;
309 default: return NLOPT_SUCCESS;
317 return luksan_plis(n, f, f_data, lb, ub, x, minf, &stop);
321 return luksan_plip(n, f, f_data, lb, ub, x, minf, &stop,
322 algorithm == NLOPT_LD_VAR1 ? 1 : 2);
324 case NLOPT_LD_TNEWTON:
325 case NLOPT_LD_TNEWTON_RESTART:
326 case NLOPT_LD_TNEWTON_PRECOND:
327 case NLOPT_LD_TNEWTON_PRECOND_RESTART:
328 return luksan_pnet(n, f, f_data, lb, ub, x, minf, &stop,
329 1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
330 1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
332 case NLOPT_GN_CRS2_LM:
333 return crs_minimize(n, f, f_data, lb, ub, x, minf, &stop,
338 case NLOPT_GN_MLSL_LDS:
339 case NLOPT_GD_MLSL_LDS: {
340 nlopt_opt local_opt = opt->local_opt;
342 if (!local_opt) { /* default */
343 local_opt = nlopt_create((algorithm == NLOPT_GN_MLSL ||
344 algorithm == NLOPT_GN_MLSL_LDS)
345 ? nlopt_local_search_alg_nonderiv
346 : nlopt_local_search_alg_deriv, n);
347 if (!local_opt) return NLOPT_FAILURE;
348 nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
349 nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
350 nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
351 nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
352 nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
353 nlopt_set_initial_step(local_opt, opt->dx);
355 for (i = 0; i < n && stop.xtol_abs[i] > 0; ++i) ;
356 if (local_opt->ftol_rel <= 0 && local_opt->ftol_abs <= 0 &&
357 local_opt->xtol_rel <= 0 && i < n) {
358 /* it is not sensible to call MLSL without *some*
359 nonzero tolerance for the local search */
360 nlopt_set_ftol_rel(local_opt, 1e-15);
361 nlopt_set_xtol_rel(local_opt, 1e-7);
363 ret = mlsl_minimize(n, f, f_data, lb, ub, x, minf, &stop,
365 algorithm >= NLOPT_GN_MLSL_LDS);
366 if (!opt->local_opt) nlopt_destroy(local_opt);
373 #define LO(param, def) (opt->local_opt ? opt->local_opt->param : (def))
374 dual_opt = nlopt_create(LO(algorithm,
375 nlopt_local_search_alg_deriv),
377 if (!dual_opt) return NLOPT_FAILURE;
378 nlopt_set_ftol_rel(dual_opt, LO(ftol_rel, 1e-12));
379 nlopt_set_ftol_abs(dual_opt, LO(ftol_abs, 0.0));
380 nlopt_set_maxeval(dual_opt, LO(maxeval, 100000));
383 ret = mma_minimize(n, f, f_data, opt->m, opt->fc,
384 lb, ub, x, minf, &stop, dual_opt);
385 nlopt_destroy(dual_opt);
389 case NLOPT_LN_COBYLA: {
391 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
392 return NLOPT_OUT_OF_MEMORY;
393 return cobyla_minimize(n, f, f_data,
395 lb, ub, x, minf, &stop,
399 case NLOPT_LN_NEWUOA: {
401 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
402 return NLOPT_OUT_OF_MEMORY;
403 return newuoa(n, 2*n+1, x, 0, 0, step,
404 &stop, minf, f_noderiv, &d);
407 case NLOPT_LN_NEWUOA_BOUND: {
409 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
410 return NLOPT_OUT_OF_MEMORY;
411 return newuoa(n, 2*n+1, x, lb, ub, step,
412 &stop, minf, f_noderiv, &d);
415 case NLOPT_LN_BOBYQA: {
417 if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
418 return NLOPT_OUT_OF_MEMORY;
419 return bobyqa(n, 2*n+1, x, lb, ub, step,
420 &stop, minf, f_noderiv, &d);
423 case NLOPT_LN_NELDERMEAD:
430 if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
431 return NLOPT_OUT_OF_MEMORY;
433 if (algorithm == NLOPT_LN_NELDERMEAD)
434 ret= nldrmd_minimize(n,f,f_data,lb,ub,x,minf,opt->dx,&stop);
436 ret= sbplx_minimize(n,f,f_data,lb,ub,x,minf,opt->dx,&stop);
437 if (freedx) { free(opt->dx); opt->dx = NULL; }
441 case NLOPT_LN_AUGLAG:
442 case NLOPT_LN_AUGLAG_EQ:
443 case NLOPT_LD_AUGLAG:
444 case NLOPT_LD_AUGLAG_EQ: {
445 nlopt_opt local_opt = opt->local_opt;
447 if (!local_opt) { /* default */
448 local_opt = nlopt_create(
449 algorithm == NLOPT_LN_AUGLAG ||
450 algorithm == NLOPT_LN_AUGLAG_EQ
451 ? nlopt_local_search_alg_nonderiv
452 : nlopt_local_search_alg_deriv, n);
453 if (!local_opt) return NLOPT_FAILURE;
454 nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
455 nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
456 nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
457 nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
458 nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
459 nlopt_set_initial_step(local_opt, opt->dx);
461 ret = auglag_minimize(n, f, f_data,
464 lb, ub, x, minf, &stop,
466 algorithm == NLOPT_LN_AUGLAG_EQ
467 || algorithm == NLOPT_LD_AUGLAG_EQ);
468 if (!opt->local_opt) nlopt_destroy(local_opt);
473 return isres_minimize(n, f, f_data,
476 lb, ub, x, minf, &stop,
480 return NLOPT_INVALID_ARGS;
483 return NLOPT_SUCCESS; /* never reached */
486 /*********************************************************************/
488 nlopt_result nlopt_optimize_limited(nlopt_opt opt, double *x, double *minf,
489 int maxeval, double maxtime)
495 if (!opt) return NLOPT_INVALID_ARGS;
497 save_maxeval = nlopt_get_maxeval(opt);
498 save_maxtime = nlopt_get_maxtime(opt);
500 /* override opt limits if maxeval and/or maxtime are more stringent */
501 if (save_maxeval <= 0 || (maxeval > 0 && maxeval < save_maxeval))
502 nlopt_set_maxeval(opt, maxeval);
503 if (save_maxtime <= 0 || (maxtime > 0 && maxtime < save_maxtime))
504 nlopt_set_maxtime(opt, maxtime);
506 ret = nlopt_optimize(opt, x, minf);
508 nlopt_set_maxeval(opt, save_maxeval);
509 nlopt_set_maxtime(opt, save_maxtime);
514 /*********************************************************************/