From 2ba6e6f6866ff706fea1ca46177ffd09f2af306b Mon Sep 17 00:00:00 2001 From: stevenj Date: Sat, 3 Apr 2010 12:00:38 -0400 Subject: [PATCH] new, extensible "object-oriented" API, first draft darcs-hash:20100403160038-c8de0-232c66be6bea9251b847be53380e00a5952f695d.gz --- api/Makefile.am | 3 +- api/deprecated.c | 161 +++++++++++ api/general.c | 115 ++++++++ api/nlopt-internal.h | 82 ++++++ api/nlopt.c | 659 ------------------------------------------- api/nlopt.h | 126 ++++++++- api/optimize.c | 511 +++++++++++++++++++++++++++++++++ api/options.c | 453 +++++++++++++++++++++++++++++ auglag/auglag.c | 108 ++++--- auglag/auglag.h | 8 +- cobyla/cobyla.c | 11 +- cobyla/cobyla.h | 3 +- isres/isres.c | 31 +- isres/isres.h | 6 +- mlsl/mlsl.c | 26 +- mlsl/mlsl.h | 3 +- mma/mma.c | 44 ++- mma/mma.h | 6 +- util/Makefile.am | 2 + util/nlopt-util.h | 19 +- 20 files changed, 1580 insertions(+), 797 deletions(-) create mode 100644 api/deprecated.c create mode 100644 api/general.c create mode 100644 api/nlopt-internal.h delete mode 100644 api/nlopt.c create mode 100644 api/optimize.c create mode 100644 api/options.c diff --git a/api/Makefile.am b/api/Makefile.am index 3e2f5ce..7b9cfb4 100644 --- a/api/Makefile.am +++ b/api/Makefile.am @@ -4,7 +4,8 @@ include_HEADERS = nlopt.h nlopt.f noinst_LTLIBRARIES = libapi.la dist_man_MANS = nlopt_minimize.3 nlopt_minimize_constrained.3 -libapi_la_SOURCES = nlopt.c nlopt.h f77api.c f77funcs.h +libapi_la_SOURCES = general.c options.c optimize.c deprecated.c \ +nlopt-internal.h nlopt.h f77api.c f77funcs.h BUILT_SOURCES = nlopt.f diff --git a/api/deprecated.c b/api/deprecated.c new file mode 100644 index 0000000..11539bd --- /dev/null +++ b/api/deprecated.c @@ -0,0 +1,161 @@ +/* Copyright (c) 2007-2008 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include "nlopt.h" + +/*************************************************************************/ + +nlopt_algorithm nlopt_local_search_alg_deriv = NLOPT_LD_MMA; +nlopt_algorithm nlopt_local_search_alg_nonderiv = NLOPT_LN_COBYLA; +int nlopt_local_search_maxeval = -1; /* no maximum by default */ + +void nlopt_get_local_search_algorithm(nlopt_algorithm *deriv, + nlopt_algorithm *nonderiv, + int *maxeval) +{ + *deriv = nlopt_local_search_alg_deriv; + *nonderiv = nlopt_local_search_alg_nonderiv; + *maxeval = nlopt_local_search_maxeval; +} + +void nlopt_set_local_search_algorithm(nlopt_algorithm deriv, + nlopt_algorithm nonderiv, + int maxeval) +{ + nlopt_local_search_alg_deriv = deriv; + nlopt_local_search_alg_nonderiv = nonderiv; + nlopt_local_search_maxeval = maxeval; +} + +/*************************************************************************/ + +int nlopt_stochastic_population = 0; + +int nlopt_get_stochastic_population(void) { + return nlopt_stochastic_population; } +void nlopt_set_stochastic_population(int pop) { + nlopt_stochastic_population = pop; } + +/*************************************************************************/ + +nlopt_result nlopt_minimize_econstrained( + nlopt_algorithm algorithm, + int n, nlopt_func f, void *f_data, + int m, nlopt_func fc, void *fc_data_, ptrdiff_t fc_datum_size, + int p, nlopt_func h, void *h_data_, ptrdiff_t h_datum_size, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, /* out: minimum */ + double minf_max, double ftol_rel, double ftol_abs, + double xtol_rel, const double *xtol_abs, + double htol_rel, double htol_abs, + int maxeval, double maxtime) +{ + char *fc_data = (char *) fc_data_; + char *h_data = (char *) h_data_; + nlopt_opt opt = nlopt_create(algorithm, n); + nlopt_result ret; + int i; + + if (!opt) return NLOPT_INVALID_ARGS; + + ret = nlopt_set_min_objective(opt, f, f_data); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + + for (i = 0; i < m; ++i) { + ret = nlopt_add_inequality_constraint(opt, fc, + fc_data + i*fc_datum_size, + 0.0); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + } + + (void) htol_rel; /* unused */ + for (i = 0; i < p; ++i) { + ret = nlopt_add_equality_constraint(opt, h, + h_data + i*h_datum_size, + htol_abs); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + } + + ret = nlopt_set_lower_bounds(opt, lb); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + ret = nlopt_set_upper_bounds(opt, ub); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + + ret = nlopt_set_stopval(opt, minf_max); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + + ret = nlopt_set_ftol_rel(opt, ftol_rel); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + ret = nlopt_set_ftol_abs(opt, ftol_abs); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + + ret = nlopt_set_xtol_rel(opt, xtol_rel); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + ret = nlopt_set_xtol_abs(opt, xtol_abs); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + + ret = nlopt_set_maxeval(opt, maxeval); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + + ret = nlopt_set_maxtime(opt, maxtime); + if (ret != NLOPT_SUCCESS) { nlopt_destroy(opt); return ret; } + + ret = nlopt_optimize(opt, x, minf); + + nlopt_destroy(opt); + return ret; +} + +nlopt_result nlopt_minimize_constrained( + nlopt_algorithm algorithm, + int n, nlopt_func f, void *f_data, + int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, /* out: minimum */ + double minf_max, double ftol_rel, double ftol_abs, + double xtol_rel, const double *xtol_abs, + int maxeval, double maxtime) +{ + return nlopt_minimize_econstrained( + algorithm, n, f, f_data, + m, fc, fc_data, fc_datum_size, 0, NULL, NULL, 0, + lb, ub, x, minf, minf_max, ftol_rel, ftol_abs, + xtol_rel, xtol_abs, ftol_rel, ftol_abs, maxeval, maxtime); +} + +nlopt_result nlopt_minimize( + nlopt_algorithm algorithm, + int n, nlopt_func f, void *f_data, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, /* out: minimum */ + double minf_max, double ftol_rel, double ftol_abs, + double xtol_rel, const double *xtol_abs, + int maxeval, double maxtime) +{ + return nlopt_minimize_constrained( + algorithm, n, f, f_data, 0, NULL, NULL, 0, + lb, ub, x, minf, minf_max, ftol_rel, ftol_abs, + xtol_rel, xtol_abs, maxeval, maxtime); +} diff --git a/api/general.c b/api/general.c new file mode 100644 index 0000000..3cb66c0 --- /dev/null +++ b/api/general.c @@ -0,0 +1,115 @@ +/* Copyright (c) 2007-2008 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include + +#include "nlopt.h" +#include "nlopt-util.h" + +/*************************************************************************/ + +int nlopt_isinf(double x) { + return fabs(x) >= HUGE_VAL * 0.99 +#ifdef HAVE_ISINF + || isinf(x) +#endif + ; +} +/*************************************************************************/ + +void nlopt_version(int *major, int *minor, int *bugfix) +{ + *major = MAJOR_VERSION; + *minor = MINOR_VERSION; + *bugfix = BUGFIX_VERSION; +} + +/*************************************************************************/ + +static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = { + "DIRECT (global, no-derivative)", + "DIRECT-L (global, no-derivative)", + "Randomized DIRECT-L (global, no-derivative)", + "Unscaled DIRECT (global, no-derivative)", + "Unscaled DIRECT-L (global, no-derivative)", + "Unscaled Randomized DIRECT-L (global, no-derivative)", + "Original DIRECT version (global, no-derivative)", + "Original DIRECT-L version (global, no-derivative)", +#ifdef WITH_CXX + "StoGO (global, derivative-based)", + "StoGO with randomized search (global, derivative-based)", +#else + "StoGO (NOT COMPILED)", + "StoGO randomized (NOT COMPILED)", +#endif +#ifdef WITH_NOCEDAL_LBFGS + "original NON-FREE L-BFGS code by Nocedal et al. (local, deriv.-based)", +#else + "original NON-FREE L-BFGS code by Nocedal et al. (NOT COMPILED)", +#endif + "Limited-memory BFGS (L-BFGS) (local, derivative-based)", + "Principal-axis, praxis (local, no-derivative)", + "Limited-memory variable-metric, rank 1 (local, derivative-based)", + "Limited-memory variable-metric, rank 2 (local, derivative-based)", + "Truncated Newton (local, derivative-based)", + "Truncated Newton with restarting (local, derivative-based)", + "Preconditioned truncated Newton (local, derivative-based)", + "Preconditioned truncated Newton with restarting (local, derivative-based)", + "Controlled random search (CRS2) with local mutation (global, no-derivative)", + "Multi-level single-linkage (MLSL), random (global, no-derivative)", + "Multi-level single-linkage (MLSL), random (global, derivative)", + "Multi-level single-linkage (MLSL), quasi-random (global, no-derivative)", + "Multi-level single-linkage (MLSL), quasi-random (global, derivative)", + "Method of Moving Asymptotes (MMA) (local, derivative)", + "COBYLA (Constrained Optimization BY Linear Approximations) (local, no-derivative)", + "NEWUOA unconstrained optimization via quadratic models (local, no-derivative)", + "Bound-constrained optimization via NEWUOA-based quadratic models (local, no-derivative)", + "Nelder-Mead simplex algorithm (local, no-derivative)", + "Sbplx variant of Nelder-Mead (re-implementation of Rowan's Subplex) (local, no-derivative)", + "Augmented Lagrangian method (local, no-derivative)", + "Augmented Lagrangian method (local, derivative)", + "Augmented Lagrangian method for equality constraints (local, no-derivative)", + "Augmented Lagrangian method for equality constraints (local, derivative)", + "BOBYQA bound-constrained optimization via quadratic models (local, no-derivative)", + "ISRES evolutionary constrained optimization (global, no-derivative)", +}; + +const char *nlopt_algorithm_name(nlopt_algorithm a) +{ + if (a < 0 || a >= NLOPT_NUM_ALGORITHMS) return "UNKNOWN"; + return nlopt_algorithm_names[a]; +} + +/*************************************************************************/ + +int nlopt_srand_called = 0; +void nlopt_srand(unsigned long seed) { + nlopt_srand_called = 1; + nlopt_init_genrand(seed); +} + +void nlopt_srand_time() { + nlopt_srand(nlopt_time_seed()); +} + +/*************************************************************************/ diff --git a/api/nlopt-internal.h b/api/nlopt-internal.h new file mode 100644 index 0000000..e74ea54 --- /dev/null +++ b/api/nlopt-internal.h @@ -0,0 +1,82 @@ +/* Copyright (c) 2007-2009 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#ifndef NLOPT_INTERNAL_H +#define NLOPT_INTERNAL_H + +#include "nlopt.h" +#include "nlopt-util.h" + +#ifdef __cplusplus +extern "C" +{ +#endif /* __cplusplus */ + +/*********************************************************************/ + +struct nlopt_opt_s { + nlopt_algorithm algorithm; /* the optimization algorithm (immutable) */ + int n; /* the dimension of the problem (immutable) */ + + nlopt_func f; void *f_data; /* objective function to minimize */ + + double *lb, *ub; /* lower and upper bounds (length n) */ + + int m; /* number of inequality constraints */ + int m_alloc; /* number of inequality constraints allocated */ + nlopt_constraint *fc; /* inequality constraints, length m_alloc */ + + int p; /* number of equality constraints */ + int p_alloc; /* number of inequality constraints allocated */ + nlopt_constraint *h; /* equality constraints, length p_alloc */ + + /* stopping criteria */ + double minf_max; /* stop when f < minf_max */ + double ftol_rel, ftol_abs; /* relative/absolute f tolerances */ + double xtol_rel, *xtol_abs; /* rel/abs x tolerances */ + int maxeval; /* max # evaluations */ + double maxtime; /* max time (seconds) */ + + /* algorithm-specific parameters */ + nlopt_opt local_opt; /* local optimizer */ + int stochastic_population; /* population size for stochastic algs */ + double *dx; /* initial step sizes (length n) for nonderivative algs */ +}; + +/*********************************************************************/ +extern int nlopt_srand_called; /* whether the random seed is initialized */ + +/*********************************************************************/ +/* global defaults set by deprecated API: */ + +extern nlopt_algorithm nlopt_local_search_alg_deriv; +extern nlopt_algorithm nlopt_local_search_alg_nonderiv; +extern int nlopt_local_search_maxeval; +extern int nlopt_stochastic_population; + +/*********************************************************************/ + +#ifdef __cplusplus +} /* extern "C" */ +#endif /* __cplusplus */ + +#endif /* NLOPT_INTERNAL_H */ diff --git a/api/nlopt.c b/api/nlopt.c deleted file mode 100644 index 7696ca4..0000000 --- a/api/nlopt.c +++ /dev/null @@ -1,659 +0,0 @@ -/* Copyright (c) 2007-2008 Massachusetts Institute of Technology - * - * Permission is hereby granted, free of charge, to any person obtaining - * a copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be - * included in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE - * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION - * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION - * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -#include -#include -#include -#include - -#include "nlopt.h" -#include "nlopt-util.h" - -#define MIN(a,b) ((a) < (b) ? (a) : (b)) - -/*************************************************************************/ - -int nlopt_isinf(double x) { - return fabs(x) >= HUGE_VAL * 0.99 -#ifdef HAVE_ISINF - || isinf(x) -#endif - ; -} - -#ifndef HAVE_ISNAN -static int my_isnan(double x) { return x != x; } -# define isnan my_isnan -#endif - -/*************************************************************************/ - -void nlopt_version(int *major, int *minor, int *bugfix) -{ - *major = MAJOR_VERSION; - *minor = MINOR_VERSION; - *bugfix = BUGFIX_VERSION; -} - -/*************************************************************************/ - -static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = { - "DIRECT (global, no-derivative)", - "DIRECT-L (global, no-derivative)", - "Randomized DIRECT-L (global, no-derivative)", - "Unscaled DIRECT (global, no-derivative)", - "Unscaled DIRECT-L (global, no-derivative)", - "Unscaled Randomized DIRECT-L (global, no-derivative)", - "Original DIRECT version (global, no-derivative)", - "Original DIRECT-L version (global, no-derivative)", -#ifdef WITH_CXX - "StoGO (global, derivative-based)", - "StoGO with randomized search (global, derivative-based)", -#else - "StoGO (NOT COMPILED)", - "StoGO randomized (NOT COMPILED)", -#endif -#ifdef WITH_NOCEDAL_LBFGS - "original NON-FREE L-BFGS code by Nocedal et al. (local, deriv.-based)", -#else - "original NON-FREE L-BFGS code by Nocedal et al. (NOT COMPILED)", -#endif - "Limited-memory BFGS (L-BFGS) (local, derivative-based)", - "Principal-axis, praxis (local, no-derivative)", - "Limited-memory variable-metric, rank 1 (local, derivative-based)", - "Limited-memory variable-metric, rank 2 (local, derivative-based)", - "Truncated Newton (local, derivative-based)", - "Truncated Newton with restarting (local, derivative-based)", - "Preconditioned truncated Newton (local, derivative-based)", - "Preconditioned truncated Newton with restarting (local, derivative-based)", - "Controlled random search (CRS2) with local mutation (global, no-derivative)", - "Multi-level single-linkage (MLSL), random (global, no-derivative)", - "Multi-level single-linkage (MLSL), random (global, derivative)", - "Multi-level single-linkage (MLSL), quasi-random (global, no-derivative)", - "Multi-level single-linkage (MLSL), quasi-random (global, derivative)", - "Method of Moving Asymptotes (MMA) (local, derivative)", - "COBYLA (Constrained Optimization BY Linear Approximations) (local, no-derivative)", - "NEWUOA unconstrained optimization via quadratic models (local, no-derivative)", - "Bound-constrained optimization via NEWUOA-based quadratic models (local, no-derivative)", - "Nelder-Mead simplex algorithm (local, no-derivative)", - "Sbplx variant of Nelder-Mead (re-implementation of Rowan's Subplex) (local, no-derivative)", - "Augmented Lagrangian method (local, no-derivative)", - "Augmented Lagrangian method (local, derivative)", - "Augmented Lagrangian method for equality constraints (local, no-derivative)", - "Augmented Lagrangian method for equality constraints (local, derivative)", - "BOBYQA bound-constrained optimization via quadratic models (local, no-derivative)", - "ISRES evolutionary constrained optimization (global, no-derivative)", -}; - -const char *nlopt_algorithm_name(nlopt_algorithm a) -{ - if (a < 0 || a >= NLOPT_NUM_ALGORITHMS) return "UNKNOWN"; - return nlopt_algorithm_names[a]; -} - -/*************************************************************************/ - -static int nlopt_srand_called = 0; -void nlopt_srand(unsigned long seed) { - nlopt_srand_called = 1; - nlopt_init_genrand(seed); -} - -void nlopt_srand_time() { - nlopt_srand(nlopt_time_seed()); -} - -/*************************************************************************/ - -/* crude heuristics for initial step size of nonderivative algorithms */ -static double initial_step(int n, - const double *lb, const double *ub, const double *x) -{ - int i; - double step = HUGE_VAL; - for (i = 0; i < n; ++i) { - if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i]) - && (ub[i] - lb[i]) * 0.25 < step && ub[i] > lb[i]) - step = (ub[i] - lb[i]) * 0.25; - if (!nlopt_isinf(ub[i]) - && ub[i] - x[i] < step && ub[i] > x[i]) - step = ub[i] - x[i]; - if (!nlopt_isinf(lb[i]) - && x[i] - lb[i] < step && x[i] > lb[i]) - step = x[i] - lb[i]; - } - if (nlopt_isinf(step)) - for (i = 0; i < n; ++i) { - if (!nlopt_isinf(ub[i]) - && ub[i] - x[i] < step && ub[i] > x[i] + 1e-4) - step = ub[i] - x[i]; - if (!nlopt_isinf(lb[i]) - && x[i] - lb[i] < step && x[i] > lb[i] + 1e-4) - step = x[i] - lb[i]; - } - if (nlopt_isinf(step)) { - step = 0; - for (i = 0; i < n; ++i) - if (fabs(x[i]) * 0.25 > step) - step = fabs(x[i]) * 0.25; - if (step == 0) - step = 1; - } - return step; -} - -/*************************************************************************/ - -typedef struct { - nlopt_func f; - void *f_data; - const double *lb, *ub; -} nlopt_data; - -#include "praxis.h" - -static double f_bound(int n, const double *x, void *data_) -{ - int i; - nlopt_data *data = (nlopt_data *) data_; - double f; - - /* some methods do not support bound constraints, but support - discontinuous objectives so we can just return Inf for invalid x */ - for (i = 0; i < n; ++i) - if (x[i] < data->lb[i] || x[i] > data->ub[i]) - return HUGE_VAL; - - f = data->f(n, x, NULL, data->f_data); - return (isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f); -} - -static double f_noderiv(int n, const double *x, void *data_) -{ - nlopt_data *data = (nlopt_data *) data_; - return data->f(n, x, NULL, data->f_data); -} - -#include "direct.h" - -static double f_direct(int n, const double *x, int *undefined, void *data_) -{ - nlopt_data *data = (nlopt_data *) data_; - double f; - f = data->f(n, x, NULL, data->f_data); - *undefined = isnan(f) || nlopt_isinf(f); - return f; -} - -#ifdef WITH_CXX -# include "stogo.h" -#endif - -#include "cdirect.h" - -#ifdef WITH_NOCEDAL -# include "l-bfgs-b.h" -#endif - -#include "luksan.h" - -#include "crs.h" - -#include "mlsl.h" -#include "mma.h" -#include "cobyla.h" -#include "newuoa.h" -#include "neldermead.h" -#include "auglag.h" -#include "bobyqa.h" -#include "isres.h" - -#define AUGLAG_ALG(a) ((a) == NLOPT_LN_AUGLAG || \ - (a) == NLOPT_LN_AUGLAG_EQ || \ - (a) == NLOPT_LD_AUGLAG || \ - (a) == NLOPT_LD_AUGLAG_EQ) - -/*************************************************************************/ - -/* for "hybrid" algorithms that combine global search with some - local search algorithm, most of the time we anticipate that people - will stick with the default local search algorithm, so we - don't add this as a parameter to nlopt_minimize. Instead, the user - can call nlopt_{set/get}_hybrid_local_algorithm to get/set the defaults. */ - -/* default local-search algorithms */ -static nlopt_algorithm local_search_alg_deriv = NLOPT_LD_MMA; -static nlopt_algorithm local_search_alg_nonderiv = NLOPT_LN_COBYLA; - -static int local_search_maxeval = -1; /* no maximum by default */ - -void nlopt_get_local_search_algorithm(nlopt_algorithm *deriv, - nlopt_algorithm *nonderiv, - int *maxeval) -{ - *deriv = local_search_alg_deriv; - *nonderiv = local_search_alg_nonderiv; - *maxeval = local_search_maxeval; -} - -void nlopt_set_local_search_algorithm(nlopt_algorithm deriv, - nlopt_algorithm nonderiv, - int maxeval) -{ - local_search_alg_deriv = deriv; - local_search_alg_nonderiv = nonderiv; - local_search_maxeval = maxeval; -} - -/*************************************************************************/ - -/* For stochastic algorithms, there is often an initial "population" - size to seed the search. Like above, we let the user - call nlopt_{set/get}_stochastic_population in order to get/set the - defaults. The special stochastic population size of "0" means - that the optimization algorithm should pick its default population. */ - -static int stochastic_population = 0; - -int nlopt_get_stochastic_population(void) { return stochastic_population; } -void nlopt_set_stochastic_population(int pop) { stochastic_population = pop; } - -/*************************************************************************/ - -/* same as nlopt_minimize_econstrained, - but xtol_abs is required to be non-NULL */ -static nlopt_result nlopt_minimize_( - nlopt_algorithm algorithm, - int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size, - int p, nlopt_func h, void *h_data, ptrdiff_t h_datum_size, - const double *lb, const double *ub, /* bounds */ - double *x, /* in: initial guess, out: minimizer */ - double *minf, /* out: minimum */ - double minf_max, double ftol_rel, double ftol_abs, - double xtol_rel, const double *xtol_abs, - double htol_rel, double htol_abs, - int maxeval, double maxtime) -{ - int i; - nlopt_data d; - nlopt_stopping stop; - - /* some basic argument checks */ - if (!minf || !f || n < 0 || m < 0 || p < 0 || (p > 0 && !h) - || (m > 0 && !fc)) return NLOPT_INVALID_ARGS; - if (n == 0) { /* trivial case: no degrees of freedom */ - *minf = f(n, x, NULL, f_data); - return NLOPT_SUCCESS; - } - else if (n < 0 || !lb || !ub || !x) - return NLOPT_INVALID_ARGS; - - /* nonlinear constraints are only supported with some algorithms */ - if (m != 0 && algorithm != NLOPT_LD_MMA && algorithm != NLOPT_LN_COBYLA - && !AUGLAG_ALG(algorithm) && algorithm != NLOPT_GN_ISRES) - return NLOPT_INVALID_ARGS; - - /* nonlinear equality constraints (h(x) = 0) only via some algorithms */ - if (p != 0 && !AUGLAG_ALG(algorithm) && algorithm != NLOPT_GN_ISRES) - return NLOPT_INVALID_ARGS; - - *minf = HUGE_VAL; - - d.f = f; - d.f_data = f_data; - d.lb = lb; - d.ub = ub; - - /* make sure rand generator is inited */ - if (!nlopt_srand_called) - nlopt_srand_time(); /* default is non-deterministic */ - - /* check bound constraints */ - for (i = 0; i < n; ++i) - if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i]) - return NLOPT_INVALID_ARGS; - - stop.n = n; - stop.minf_max = (isnan(minf_max) - || (nlopt_isinf(minf_max) && minf_max < 0)) - ? -HUGE_VAL : minf_max; - stop.ftol_rel = ftol_rel; - stop.ftol_abs = ftol_abs; - stop.xtol_rel = xtol_rel; - stop.xtol_abs = xtol_abs; - stop.htol_rel = htol_rel; - stop.htol_abs = htol_abs; - stop.nevals = 0; - stop.maxeval = maxeval; - stop.maxtime = maxtime; - stop.start = nlopt_seconds(); - - switch (algorithm) { - case NLOPT_GN_DIRECT: - case NLOPT_GN_DIRECT_L: - case NLOPT_GN_DIRECT_L_RAND: - return cdirect(n, 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)) - + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT))); - - case NLOPT_GN_DIRECT_NOSCAL: - case NLOPT_GN_DIRECT_L_NOSCAL: - case NLOPT_GN_DIRECT_L_RAND_NOSCAL: - return cdirect_unscaled(n, 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)) - + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT))); - - case NLOPT_GN_ORIG_DIRECT: - case NLOPT_GN_ORIG_DIRECT_L: - switch (direct_optimize(f_direct, &d, n, lb, ub, x, minf, - maxeval, -1, 0.0, 0.0, - pow(xtol_rel, (double) n), -1.0, - stop.minf_max, 0.0, - NULL, - algorithm == NLOPT_GN_ORIG_DIRECT - ? DIRECT_ORIGINAL - : DIRECT_GABLONSKY)) { - case DIRECT_INVALID_BOUNDS: - case DIRECT_MAXFEVAL_TOOBIG: - case DIRECT_INVALID_ARGS: - return NLOPT_INVALID_ARGS; - case DIRECT_INIT_FAILED: - case DIRECT_SAMPLEPOINTS_FAILED: - case DIRECT_SAMPLE_FAILED: - return NLOPT_FAILURE; - case DIRECT_MAXFEVAL_EXCEEDED: - case DIRECT_MAXITER_EXCEEDED: - return NLOPT_MAXEVAL_REACHED; - case DIRECT_GLOBAL_FOUND: - return NLOPT_MINF_MAX_REACHED; - case DIRECT_VOLTOL: - case DIRECT_SIGMATOL: - return NLOPT_XTOL_REACHED; - case DIRECT_OUT_OF_MEMORY: - return NLOPT_OUT_OF_MEMORY; - break; - } - - case NLOPT_GD_STOGO: - case NLOPT_GD_STOGO_RAND: -#ifdef WITH_CXX - if (!stogo_minimize(n, f, f_data, x, minf, lb, ub, &stop, - algorithm == NLOPT_GD_STOGO - ? 0 : 2*n)) - return NLOPT_FAILURE; - break; -#else - return NLOPT_FAILURE; -#endif - -#if 0 - /* lacking a free/open-source license, we no longer use - Rowan's code, and instead use by "sbplx" re-implementation */ - case NLOPT_LN_SUBPLEX: { - int iret; - double *scale = (double *) malloc(sizeof(double) * n); - if (!scale) return NLOPT_OUT_OF_MEMORY; - for (i = 0; i < n; ++i) - scale[i] = initial_step(1, lb+i, ub+i, x+i); - iret = nlopt_subplex(f_bound, minf, x, n, &d, &stop, scale); - free(scale); - switch (iret) { - case -2: return NLOPT_INVALID_ARGS; - case -10: return NLOPT_MAXTIME_REACHED; - case -1: return NLOPT_MAXEVAL_REACHED; - case 0: return NLOPT_XTOL_REACHED; - case 1: return NLOPT_SUCCESS; - case 2: return NLOPT_MINF_MAX_REACHED; - case 20: return NLOPT_FTOL_REACHED; - case -200: return NLOPT_OUT_OF_MEMORY; - default: return NLOPT_FAILURE; /* unknown return code */ - } - break; - } -#endif - - case NLOPT_LN_PRAXIS: - return praxis_(0.0, DBL_EPSILON, - initial_step(n, lb, ub, x), n, x, f_bound, &d, - &stop, minf); - -#ifdef WITH_NOCEDAL - case NLOPT_LD_LBFGS_NOCEDAL: { - int iret, *nbd = (int *) malloc(sizeof(int) * n); - if (!nbd) return NLOPT_OUT_OF_MEMORY; - for (i = 0; i < n; ++i) { - int linf = nlopt_isinf(lb[i]) && lb[i] < 0; - int uinf = nlopt_isinf(ub[i]) && ub[i] > 0; - nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2)); - } - iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub, - MIN(n, 5), 0.0, ftol_rel, - xtol_abs ? *xtol_abs : xtol_rel, - maxeval); - free(nbd); - if (iret <= 0) { - switch (iret) { - case -1: return NLOPT_INVALID_ARGS; - case -2: default: return NLOPT_FAILURE; - } - } - else { - *minf = f(n, x, NULL, f_data); - switch (iret) { - case 5: return NLOPT_MAXEVAL_REACHED; - case 2: return NLOPT_XTOL_REACHED; - case 1: return NLOPT_FTOL_REACHED; - default: return NLOPT_SUCCESS; - } - } - break; - } -#endif - - case NLOPT_LD_LBFGS: - return luksan_plis(n, f, f_data, lb, ub, x, minf, &stop); - - case NLOPT_LD_VAR1: - case NLOPT_LD_VAR2: - return luksan_plip(n, f, f_data, lb, ub, x, minf, &stop, - algorithm == NLOPT_LD_VAR1 ? 1 : 2); - - case NLOPT_LD_TNEWTON: - case NLOPT_LD_TNEWTON_RESTART: - case NLOPT_LD_TNEWTON_PRECOND: - case NLOPT_LD_TNEWTON_PRECOND_RESTART: - return luksan_pnet(n, f, f_data, lb, ub, x, minf, &stop, - 1 + (algorithm - NLOPT_LD_TNEWTON) % 2, - 1 + (algorithm - NLOPT_LD_TNEWTON) / 2); - - case NLOPT_GN_CRS2_LM: - return crs_minimize(n, f, f_data, lb, ub, x, minf, &stop, - stochastic_population, 0); - - case NLOPT_GN_MLSL: - case NLOPT_GD_MLSL: - case NLOPT_GN_MLSL_LDS: - case NLOPT_GD_MLSL_LDS: - for (i = 0; i < n && stop.xtol_abs[i] <= 0; ++i) ; - if (stop.ftol_rel <= 0 && stop.ftol_abs <= 0 && - stop.xtol_rel <= 0 && i == n) { - /* it is not sensible to call MLSL without *some* - nonzero tolerance for the local search */ - stop.ftol_rel = 1e-15; - stop.xtol_rel = 1e-7; - } - return mlsl_minimize(n, f, f_data, lb, ub, x, minf, &stop, - (algorithm == NLOPT_GN_MLSL || - algorithm == NLOPT_GN_MLSL_LDS) - ? local_search_alg_nonderiv - : local_search_alg_deriv, - local_search_maxeval, - stochastic_population, - algorithm >= NLOPT_GN_MLSL_LDS); - - case NLOPT_LD_MMA: - return mma_minimize(n, f, f_data, m, fc, fc_data, fc_datum_size, - lb, ub, x, minf, &stop, - local_search_alg_deriv, 1e-12, 100000); - - case NLOPT_LN_COBYLA: - return cobyla_minimize(n, f, f_data, - m, fc, fc_data, fc_datum_size, - lb, ub, x, minf, &stop, - initial_step(n, lb, ub, x)); - - case NLOPT_LN_NEWUOA: - return newuoa(n, 2*n+1, x, 0, 0, initial_step(n, lb, ub, x), - &stop, minf, f_noderiv, &d); - - case NLOPT_LN_NEWUOA_BOUND: - return newuoa(n, 2*n+1, x, lb, ub, initial_step(n, lb, ub, x), - &stop, minf, f_noderiv, &d); - - case NLOPT_LN_BOBYQA: - return bobyqa(n, 2*n+1, x, lb, ub, initial_step(n, lb, ub, x), - &stop, minf, f_noderiv, &d); - - case NLOPT_LN_NELDERMEAD: - case NLOPT_LN_SBPLX: - { - nlopt_result ret; - double *scale = (double *) malloc(sizeof(double) * n); - if (!scale) return NLOPT_OUT_OF_MEMORY; - for (i = 0; i < n; ++i) - scale[i] = initial_step(1, lb+i, ub+i, x+i); - if (algorithm == NLOPT_LN_NELDERMEAD) - ret = nldrmd_minimize(n,f,f_data,lb,ub,x,minf,scale,&stop); - else - ret = sbplx_minimize(n,f,f_data,lb,ub,x,minf,scale,&stop); - free(scale); - return ret; - } - - case NLOPT_LN_AUGLAG: - case NLOPT_LN_AUGLAG_EQ: - return auglag_minimize(n, f, f_data, - m, fc, fc_data, fc_datum_size, - p, h, h_data, h_datum_size, - lb, ub, x, minf, &stop, - local_search_alg_nonderiv, - algorithm == NLOPT_LN_AUGLAG_EQ); - - case NLOPT_LD_AUGLAG: - case NLOPT_LD_AUGLAG_EQ: - return auglag_minimize(n, f, f_data, - m, fc, fc_data, fc_datum_size, - p, h, h_data, h_datum_size, - lb, ub, x, minf, &stop, - local_search_alg_deriv, - algorithm == NLOPT_LD_AUGLAG_EQ); - - case NLOPT_GN_ISRES: - return isres_minimize(n, f, f_data, - m, fc, fc_data, fc_datum_size, - p, h, h_data, h_datum_size, - lb, ub, x, minf, &stop, - stochastic_population); - - default: - return NLOPT_INVALID_ARGS; - } - - return NLOPT_SUCCESS; -} - -nlopt_result nlopt_minimize_econstrained( - nlopt_algorithm algorithm, - int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size, - int p, nlopt_func h, void *h_data, ptrdiff_t h_datum_size, - const double *lb, const double *ub, /* bounds */ - double *x, /* in: initial guess, out: minimizer */ - double *minf, /* out: minimum */ - double minf_max, double ftol_rel, double ftol_abs, - double xtol_rel, const double *xtol_abs, - double htol_rel, double htol_abs, - int maxeval, double maxtime) -{ - nlopt_result ret; - if (xtol_abs) - ret = nlopt_minimize_(algorithm, n, f, f_data, - m, fc, fc_data, fc_datum_size, - p, h, h_data, h_datum_size, - lb, ub, - x, minf, minf_max, ftol_rel, ftol_abs, - xtol_rel, xtol_abs, htol_rel, htol_abs, - maxeval, maxtime); - else { - int i; - double *xtol = (double *) malloc(sizeof(double) * n); - if (!xtol) return NLOPT_OUT_OF_MEMORY; - for (i = 0; i < n; ++i) xtol[i] = -1; - ret = nlopt_minimize_(algorithm, n, f, f_data, - m, fc, fc_data, fc_datum_size, - p, h, h_data, h_datum_size, - lb, ub, - x, minf, minf_max, ftol_rel, ftol_abs, - xtol_rel, xtol, htol_rel, htol_abs, - maxeval, maxtime); - free(xtol); - } - return ret; -} - -nlopt_result nlopt_minimize_constrained( - nlopt_algorithm algorithm, - int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size, - const double *lb, const double *ub, /* bounds */ - double *x, /* in: initial guess, out: minimizer */ - double *minf, /* out: minimum */ - double minf_max, double ftol_rel, double ftol_abs, - double xtol_rel, const double *xtol_abs, - int maxeval, double maxtime) -{ - return nlopt_minimize_econstrained( - algorithm, n, f, f_data, - m, fc, fc_data, fc_datum_size, 0, NULL, NULL, 0, - lb, ub, x, minf, minf_max, ftol_rel, ftol_abs, - xtol_rel, xtol_abs, ftol_rel, ftol_abs, maxeval, maxtime); -} - -nlopt_result nlopt_minimize( - nlopt_algorithm algorithm, - int n, nlopt_func f, void *f_data, - const double *lb, const double *ub, /* bounds */ - double *x, /* in: initial guess, out: minimizer */ - double *minf, /* out: minimum */ - double minf_max, double ftol_rel, double ftol_abs, - double xtol_rel, const double *xtol_abs, - int maxeval, double maxtime) -{ - return nlopt_minimize_constrained( - algorithm, n, f, f_data, 0, NULL, NULL, 0, - lb, ub, x, minf, minf_max, ftol_rel, ftol_abs, - xtol_rel, xtol_abs, maxeval, maxtime); -} diff --git a/api/nlopt.h b/api/nlopt.h index 61ae819..3a91177 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -144,6 +144,111 @@ typedef enum { NLOPT_MAXTIME_REACHED = 6 } nlopt_result; + +NLOPT_EXTERN void nlopt_srand(unsigned long seed); +NLOPT_EXTERN void nlopt_srand_time(void); + +NLOPT_EXTERN void nlopt_version(int *major, int *minor, int *bugfix); + +/*************************** OBJECT-ORIENTED API **************************/ +/* The style here is that we create an nlopt_opt "object" (an opaque pointer), + then set various optimization parameters, and then execute the + algorithm. In this way, we can add more and more optimization parameters + (including algorithm-specific ones) without breaking backwards + compatibility, having functions with zillions of parameters, or + relying non-reentrantly on global variables.*/ + +struct nlopt_opt_s; /* opaque structure, defined internally */ +typedef struct nlopt_opt_s *nlopt_opt; + +/* the only immutable parameters of an optimization are the algorithm and + the dimension n of the problem, since changing either of these could + have side-effects on lots of other parameters */ +NLOPT_EXTERN nlopt_opt nlopt_create(nlopt_algorithm algorithm, int n); +NLOPT_EXTERN void nlopt_destroy(nlopt_opt opt); +NLOPT_EXTERN nlopt_opt nlopt_copy(const nlopt_opt opt); + +NLOPT_EXTERN nlopt_result nlopt_optimize(nlopt_opt opt, double *x, + double *opt_f); + +NLOPT_EXTERN nlopt_result nlopt_set_min_objective(nlopt_opt opt, nlopt_func f, + void *f_data); + +/* constraints: */ + +NLOPT_EXTERN nlopt_result nlopt_set_lower_bounds(nlopt_opt opt, + const double *lb); +NLOPT_EXTERN nlopt_result nlopt_set_lower_bounds1(nlopt_opt opt, double lb); +NLOPT_EXTERN void nlopt_get_lower_bounds(const nlopt_opt opt, double *lb); +NLOPT_EXTERN nlopt_result nlopt_set_upper_bounds(nlopt_opt opt, + const double *ub); +NLOPT_EXTERN nlopt_result nlopt_set_upper_bounds1(nlopt_opt opt, double ub); +NLOPT_EXTERN void nlopt_get_upper_bounds(const nlopt_opt opt, double *ub); + +NLOPT_EXTERN nlopt_result nlopt_remove_inequality_constraints(nlopt_opt opt); +NLOPT_EXTERN nlopt_result nlopt_add_inequality_constraint(nlopt_opt opt, + nlopt_func fc, + void *fc_data, + double tol); + +NLOPT_EXTERN nlopt_result nlopt_remove_equality_constraints(nlopt_opt opt); +NLOPT_EXTERN nlopt_result nlopt_add_equality_constraint(nlopt_opt opt, + nlopt_func h, + void *h_data, + double tol); + +/* stopping criteria: */ + +NLOPT_EXTERN nlopt_result nlopt_set_stopval(nlopt_opt opt, double minf_max); +NLOPT_EXTERN double nlopt_get_stopval(const nlopt_opt opt); + +NLOPT_EXTERN nlopt_result nlopt_set_ftol_rel(nlopt_opt opt, double tol); +NLOPT_EXTERN double nlopt_get_ftol_rel(const nlopt_opt opt); +NLOPT_EXTERN nlopt_result nlopt_set_ftol_abs(nlopt_opt opt, double tol); +NLOPT_EXTERN double nlopt_get_ftol_abs(const nlopt_opt opt); + +NLOPT_EXTERN nlopt_result nlopt_set_xtol_rel(nlopt_opt opt, double tol); +NLOPT_EXTERN double nlopt_get_xtol_rel(const nlopt_opt opt); +NLOPT_EXTERN nlopt_result nlopt_set_xtol_abs1(nlopt_opt opt, double tol); +NLOPT_EXTERN nlopt_result nlopt_set_xtol_abs(nlopt_opt opt, const double *tol); +NLOPT_EXTERN void nlopt_get_xtol_abs(const nlopt_opt opt, double *tol); + +NLOPT_EXTERN nlopt_result nlopt_set_maxeval(nlopt_opt opt, int maxeval); +NLOPT_EXTERN int nlopt_get_maxeval(nlopt_opt opt); + +NLOPT_EXTERN nlopt_result nlopt_set_maxtime(nlopt_opt opt, double maxtime); +NLOPT_EXTERN double nlopt_get_maxtime(nlopt_opt opt); + +/* more algorithm-specific parameters */ + +NLOPT_EXTERN nlopt_result nlopt_set_local_optimizer(nlopt_opt opt, + const nlopt_opt local_opt); + +NLOPT_EXTERN nlopt_result nlopt_set_population(nlopt_opt opt, int pop); +NLOPT_EXTERN int nlopt_get_population(const nlopt_opt opt); + +NLOPT_EXTERN nlopt_result nlopt_set_default_initial_step(nlopt_opt opt, + const double *x); +NLOPT_EXTERN nlopt_result nlopt_set_initial_step(nlopt_opt opt, + const double *dx); +NLOPT_EXTERN nlopt_result nlopt_set_initial_step1(nlopt_opt opt, double dx); +NLOPT_EXTERN nlopt_result nlopt_get_initial_step(const nlopt_opt opt, + const double *x, double *dx); + +/*************************** DEPRECATED API **************************/ +/* The new "object-oriented" API is preferred, since it allows us to + gracefully add new features and algorithm-specific options in a + re-entrant way, and we can automatically assume reasonable defaults + for unspecified parameters. */ + +/* Where possible (e.g. for gcc >= 3.1), enable a compiler warning + for code that uses a deprecated function */ +#if defined(__GNUC__) && (__GNUC__ > 3 || (__GNUC__==3 && __GNUC_MINOR__ > 0)) +# define NLOPT_DEPRECATED __attribute__((deprecated)) +#else +# define NLOPT_DEPRECATED +#endif + NLOPT_EXTERN nlopt_result nlopt_minimize( nlopt_algorithm algorithm, int n, nlopt_func f, void *f_data, @@ -152,7 +257,7 @@ NLOPT_EXTERN nlopt_result nlopt_minimize( double *minf, /* out: minimum */ double minf_max, double ftol_rel, double ftol_abs, double xtol_rel, const double *xtol_abs, - int maxeval, double maxtime); + int maxeval, double maxtime) NLOPT_DEPRECATED; NLOPT_EXTERN nlopt_result nlopt_minimize_constrained( nlopt_algorithm algorithm, @@ -163,7 +268,7 @@ NLOPT_EXTERN nlopt_result nlopt_minimize_constrained( double *minf, /* out: minimum */ double minf_max, double ftol_rel, double ftol_abs, double xtol_rel, const double *xtol_abs, - int maxeval, double maxtime); + int maxeval, double maxtime) NLOPT_DEPRECATED; NLOPT_EXTERN nlopt_result nlopt_minimize_econstrained( nlopt_algorithm algorithm, @@ -176,22 +281,19 @@ NLOPT_EXTERN nlopt_result nlopt_minimize_econstrained( double minf_max, double ftol_rel, double ftol_abs, double xtol_rel, const double *xtol_abs, double htol_rel, double htol_abs, - int maxeval, double maxtime); - -NLOPT_EXTERN void nlopt_srand(unsigned long seed); -NLOPT_EXTERN void nlopt_srand_time(void); - -NLOPT_EXTERN void nlopt_version(int *major, int *minor, int *bugfix); + int maxeval, double maxtime) NLOPT_DEPRECATED; NLOPT_EXTERN void nlopt_get_local_search_algorithm(nlopt_algorithm *deriv, nlopt_algorithm *nonderiv, - int *maxeval); + int *maxeval) NLOPT_DEPRECATED; NLOPT_EXTERN void nlopt_set_local_search_algorithm(nlopt_algorithm deriv, nlopt_algorithm nonderiv, - int maxeval); + int maxeval) NLOPT_DEPRECATED; + +NLOPT_EXTERN int nlopt_get_stochastic_population(void) NLOPT_DEPRECATED; +NLOPT_EXTERN void nlopt_set_stochastic_population(int pop) NLOPT_DEPRECATED; -NLOPT_EXTERN int nlopt_get_stochastic_population(void); -NLOPT_EXTERN void nlopt_set_stochastic_population(int pop); +/*********************************************************************/ #ifdef __cplusplus } /* extern "C" */ diff --git a/api/optimize.c b/api/optimize.c new file mode 100644 index 0000000..8cb7e67 --- /dev/null +++ b/api/optimize.c @@ -0,0 +1,511 @@ +/* Copyright (c) 2007-2009 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include + +#include "nlopt-internal.h" + +/*********************************************************************/ + +#ifndef HAVE_ISNAN +static int my_isnan(double x) { return x != x; } +# define isnan my_isnan +#endif + +/*********************************************************************/ + +#include "praxis.h" +#include "direct.h" + +#ifdef WITH_CXX +# include "stogo.h" +#endif + +#include "cdirect.h" + +#ifdef WITH_NOCEDAL +# include "l-bfgs-b.h" +#endif + +#include "luksan.h" + +#include "crs.h" + +#include "mlsl.h" +#include "mma.h" +#include "cobyla.h" +#include "newuoa.h" +#include "neldermead.h" +#include "auglag.h" +#include "bobyqa.h" +#include "isres.h" + +/*********************************************************************/ + +typedef struct { + nlopt_func f; + void *f_data; + const double *lb, *ub; +} nlopt_data; + +#include "praxis.h" + +static double f_bound(int n, const double *x, void *data_) +{ + int i; + nlopt_data *data = (nlopt_data *) data_; + double f; + + /* some methods do not support bound constraints, but support + discontinuous objectives so we can just return Inf for invalid x */ + for (i = 0; i < n; ++i) + if (x[i] < data->lb[i] || x[i] > data->ub[i]) + return HUGE_VAL; + + f = data->f(n, x, NULL, data->f_data); + return (isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f); +} + +static double f_noderiv(int n, const double *x, void *data_) +{ + nlopt_data *data = (nlopt_data *) data_; + return data->f(n, x, NULL, data->f_data); +} + +static double f_direct(int n, const double *x, int *undefined, void *data_) +{ + nlopt_data *data = (nlopt_data *) data_; + double f; + f = data->f(n, x, NULL, data->f_data); + *undefined = isnan(f) || nlopt_isinf(f); + return f; +} + +/*********************************************************************/ + +/* get min(dx) for algorithms requiring a scalar initial step size */ +static nlopt_result initial_step(nlopt_opt opt, const double *x, double *step) +{ + int freedx = 0, i; + + if (!opt->dx) { + freedx = 1; + if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS) + return NLOPT_OUT_OF_MEMORY; + } + + *step = HUGE_VAL; + for (i = 0; i < opt->n; ++i) + if (*step > fabs(opt->dx[i])) + *step = fabs(opt->dx[i]); + + if (freedx) { free(opt->dx); opt->dx = NULL; } + return NLOPT_SUCCESS; +} + +/*********************************************************************/ + +#define POP(defaultpop) (opt->stochastic_population > 0 ? \ + opt->stochastic_population : \ + (nlopt_stochastic_population > 0 ? \ + nlopt_stochastic_population : (defaultpop))) + +nlopt_result nlopt_optimize(nlopt_opt opt, double *x, double *minf) +{ + const double *lb, *ub; + nlopt_algorithm algorithm; + nlopt_func f; void *f_data; + int n, i; + nlopt_data d; + nlopt_stopping stop; + + if (!opt || !x || !minf || !opt->f) return NLOPT_INVALID_ARGS; + + /* copy a few params to local vars for convenience */ + n = opt->n; + lb = opt->lb; ub = opt->ub; + algorithm = opt->algorithm; + f = opt->f; f_data = opt->f_data; + + if (n == 0) { /* trivial case: no degrees of freedom */ + *minf = f(n, x, NULL, f_data); + return NLOPT_SUCCESS; + } + + *minf = HUGE_VAL; + + d.f = f; + d.f_data = f_data; + d.lb = lb; + d.ub = ub; + + /* make sure rand generator is inited */ + if (!nlopt_srand_called) + nlopt_srand_time(); /* default is non-deterministic */ + + /* check bound constraints */ + for (i = 0; i < n; ++i) + if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i]) + return NLOPT_INVALID_ARGS; + + stop.n = n; + stop.minf_max = opt->minf_max; + stop.ftol_rel = opt->ftol_rel; + stop.ftol_abs = opt->ftol_abs; + stop.xtol_rel = opt->xtol_rel; + stop.xtol_abs = opt->xtol_abs; + stop.nevals = 0; + stop.maxeval = opt->maxeval; + stop.maxtime = opt->maxtime; + stop.start = nlopt_seconds(); + + switch (algorithm) { + case NLOPT_GN_DIRECT: + case NLOPT_GN_DIRECT_L: + case NLOPT_GN_DIRECT_L_RAND: + return cdirect(n, 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)) + + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND + ? 1 : (algorithm != NLOPT_GN_DIRECT))); + + case NLOPT_GN_DIRECT_NOSCAL: + case NLOPT_GN_DIRECT_L_NOSCAL: + case NLOPT_GN_DIRECT_L_RAND_NOSCAL: + return cdirect_unscaled(n, 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)) + + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT))); + + case NLOPT_GN_ORIG_DIRECT: + case NLOPT_GN_ORIG_DIRECT_L: + switch (direct_optimize(f_direct, &d, n, lb, ub, x, minf, + stop.maxeval, -1, 0.0, 0.0, + pow(stop.xtol_rel, (double) n), -1.0, + stop.minf_max, 0.0, + NULL, + algorithm == NLOPT_GN_ORIG_DIRECT + ? DIRECT_ORIGINAL + : DIRECT_GABLONSKY)) { + case DIRECT_INVALID_BOUNDS: + case DIRECT_MAXFEVAL_TOOBIG: + case DIRECT_INVALID_ARGS: + return NLOPT_INVALID_ARGS; + case DIRECT_INIT_FAILED: + case DIRECT_SAMPLEPOINTS_FAILED: + case DIRECT_SAMPLE_FAILED: + return NLOPT_FAILURE; + case DIRECT_MAXFEVAL_EXCEEDED: + case DIRECT_MAXITER_EXCEEDED: + return NLOPT_MAXEVAL_REACHED; + case DIRECT_GLOBAL_FOUND: + return NLOPT_MINF_MAX_REACHED; + case DIRECT_VOLTOL: + case DIRECT_SIGMATOL: + return NLOPT_XTOL_REACHED; + case DIRECT_OUT_OF_MEMORY: + return NLOPT_OUT_OF_MEMORY; + break; + } + + case NLOPT_GD_STOGO: + case NLOPT_GD_STOGO_RAND: +#ifdef WITH_CXX + if (!stogo_minimize(n, f, f_data, x, minf, lb, ub, &stop, + algorithm == NLOPT_GD_STOGO + ? 0 : POP(2*n))) + return NLOPT_FAILURE; + break; +#else + return NLOPT_FAILURE; +#endif + +#if 0 + /* lacking a free/open-source license, we no longer use + Rowan's code, and instead use by "sbplx" re-implementation */ + case NLOPT_LN_SUBPLEX: { + int iret, freedx = 0; + if (!opt->dx) { + freedx = 1; + if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS) + return NLOPT_OUT_OF_MEMORY; + } + iret = nlopt_subplex(f_bound, minf, x, n, &d, &stop, opt->dx); + if (freedx) { free(opt->dx); opt->dx = NULL; } + switch (iret) { + case -2: return NLOPT_INVALID_ARGS; + case -10: return NLOPT_MAXTIME_REACHED; + case -1: return NLOPT_MAXEVAL_REACHED; + case 0: return NLOPT_XTOL_REACHED; + case 1: return NLOPT_SUCCESS; + case 2: return NLOPT_MINF_MAX_REACHED; + case 20: return NLOPT_FTOL_REACHED; + case -200: return NLOPT_OUT_OF_MEMORY; + default: return NLOPT_FAILURE; /* unknown return code */ + } + break; + } +#endif + + case NLOPT_LN_PRAXIS: { + double step; + if (initial_step(opt, x, &step) != NLOPT_SUCCESS) + return NLOPT_OUT_OF_MEMORY; + return praxis_(0.0, DBL_EPSILON, + step, n, x, f_bound, &d, &stop, minf); + } + +#ifdef WITH_NOCEDAL + case NLOPT_LD_LBFGS_NOCEDAL: { + int iret, *nbd = (int *) malloc(sizeof(int) * n); + if (!nbd) return NLOPT_OUT_OF_MEMORY; + for (i = 0; i < n; ++i) { + int linf = nlopt_isinf(lb[i]) && lb[i] < 0; + int uinf = nlopt_isinf(ub[i]) && ub[i] > 0; + nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2)); + } + iret = lbfgsb_minimize(n, f, f_data, x, nbd, lb, ub, + MIN(n, 5), 0.0, stop.ftol_rel, + stop.xtol_abs[0] > 0 ? stop.xtol_abs[0] + : stop.xtol_rel, + stop.maxeval); + free(nbd); + if (iret <= 0) { + switch (iret) { + case -1: return NLOPT_INVALID_ARGS; + case -2: default: return NLOPT_FAILURE; + } + } + else { + *minf = f(n, x, NULL, f_data); + switch (iret) { + case 5: return NLOPT_MAXEVAL_REACHED; + case 2: return NLOPT_XTOL_REACHED; + case 1: return NLOPT_FTOL_REACHED; + default: return NLOPT_SUCCESS; + } + } + break; + } +#endif + + case NLOPT_LD_LBFGS: + return luksan_plis(n, f, f_data, lb, ub, x, minf, &stop); + + case NLOPT_LD_VAR1: + case NLOPT_LD_VAR2: + return luksan_plip(n, f, f_data, lb, ub, x, minf, &stop, + algorithm == NLOPT_LD_VAR1 ? 1 : 2); + + case NLOPT_LD_TNEWTON: + case NLOPT_LD_TNEWTON_RESTART: + case NLOPT_LD_TNEWTON_PRECOND: + case NLOPT_LD_TNEWTON_PRECOND_RESTART: + return luksan_pnet(n, f, f_data, lb, ub, x, minf, &stop, + 1 + (algorithm - NLOPT_LD_TNEWTON) % 2, + 1 + (algorithm - NLOPT_LD_TNEWTON) / 2); + + case NLOPT_GN_CRS2_LM: + return crs_minimize(n, f, f_data, lb, ub, x, minf, &stop, + POP(0), 0); + + case NLOPT_GN_MLSL: + case NLOPT_GD_MLSL: + case NLOPT_GN_MLSL_LDS: + case NLOPT_GD_MLSL_LDS: { + nlopt_opt local_opt = opt->local_opt; + nlopt_result ret; + if (!local_opt) { /* default */ + local_opt = nlopt_create((algorithm == NLOPT_GN_MLSL || + algorithm == NLOPT_GN_MLSL_LDS) + ? nlopt_local_search_alg_nonderiv + : nlopt_local_search_alg_deriv, n); + if (!local_opt) return NLOPT_FAILURE; + nlopt_set_ftol_rel(local_opt, opt->ftol_rel); + nlopt_set_ftol_abs(local_opt, opt->ftol_abs); + nlopt_set_xtol_rel(local_opt, opt->xtol_rel); + nlopt_set_xtol_abs(local_opt, opt->xtol_abs); + nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval); + } + for (i = 0; i < n && stop.xtol_abs[i] > 0; ++i) ; + if (local_opt->ftol_rel <= 0 && local_opt->ftol_abs <= 0 && + local_opt->xtol_rel <= 0 && i < n) { + /* it is not sensible to call MLSL without *some* + nonzero tolerance for the local search */ + nlopt_set_ftol_rel(local_opt, 1e-15); + nlopt_set_xtol_rel(local_opt, 1e-7); + } + ret = mlsl_minimize(n, f, f_data, lb, ub, x, minf, &stop, + local_opt, POP(0), + algorithm >= NLOPT_GN_MLSL_LDS); + if (!opt->local_opt) nlopt_destroy(local_opt); + return ret; + } + + case NLOPT_LD_MMA: { + nlopt_opt dual_opt; + nlopt_result ret; +#define LO(param, def) (opt->local_opt ? opt->local_opt->param : (def)) + dual_opt = nlopt_create(LO(algorithm, + nlopt_local_search_alg_deriv), + opt->m); + if (!dual_opt) return NLOPT_FAILURE; + nlopt_set_ftol_rel(dual_opt, LO(ftol_rel, 1e-12)); + nlopt_set_ftol_abs(dual_opt, LO(ftol_abs, 0.0)); + nlopt_set_maxeval(dual_opt, LO(maxeval, 100000)); +#undef LO + + ret = mma_minimize(n, f, f_data, opt->m, opt->fc, + lb, ub, x, minf, &stop, dual_opt); + nlopt_destroy(dual_opt); + return ret; + } + + case NLOPT_LN_COBYLA: { + double step; + if (initial_step(opt, x, &step) != NLOPT_SUCCESS) + return NLOPT_OUT_OF_MEMORY; + return cobyla_minimize(n, f, f_data, + opt->m, opt->fc, + lb, ub, x, minf, &stop, + step); + } + + case NLOPT_LN_NEWUOA: { + double step; + if (initial_step(opt, x, &step) != NLOPT_SUCCESS) + return NLOPT_OUT_OF_MEMORY; + return newuoa(n, 2*n+1, x, 0, 0, step, + &stop, minf, f_noderiv, &d); + } + + case NLOPT_LN_NEWUOA_BOUND: { + double step; + if (initial_step(opt, x, &step) != NLOPT_SUCCESS) + return NLOPT_OUT_OF_MEMORY; + return newuoa(n, 2*n+1, x, lb, ub, step, + &stop, minf, f_noderiv, &d); + } + + case NLOPT_LN_BOBYQA: { + double step; + if (initial_step(opt, x, &step) != NLOPT_SUCCESS) + return NLOPT_OUT_OF_MEMORY; + return bobyqa(n, 2*n+1, x, lb, ub, step, + &stop, minf, f_noderiv, &d); + } + + case NLOPT_LN_NELDERMEAD: + case NLOPT_LN_SBPLX: + { + nlopt_result ret; + int freedx = 0; + if (!opt->dx) { + freedx = 1; + if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS) + return NLOPT_OUT_OF_MEMORY; + } + if (algorithm == NLOPT_LN_NELDERMEAD) + ret= nldrmd_minimize(n,f,f_data,lb,ub,x,minf,opt->dx,&stop); + else + ret= sbplx_minimize(n,f,f_data,lb,ub,x,minf,opt->dx,&stop); + if (freedx) { free(opt->dx); opt->dx = NULL; } + return ret; + } + + case NLOPT_LN_AUGLAG: + case NLOPT_LN_AUGLAG_EQ: + case NLOPT_LD_AUGLAG: + case NLOPT_LD_AUGLAG_EQ: { + nlopt_opt local_opt = opt->local_opt; + nlopt_result ret; + if (!local_opt) { /* default */ + local_opt = nlopt_create( + algorithm == NLOPT_LN_AUGLAG || + algorithm == NLOPT_LN_AUGLAG_EQ + ? nlopt_local_search_alg_nonderiv + : nlopt_local_search_alg_deriv, n); + if (!local_opt) return NLOPT_FAILURE; + nlopt_set_ftol_rel(local_opt, opt->ftol_rel); + nlopt_set_ftol_abs(local_opt, opt->ftol_abs); + nlopt_set_xtol_rel(local_opt, opt->xtol_rel); + nlopt_set_xtol_abs(local_opt, opt->xtol_abs); + } + ret = auglag_minimize(n, f, f_data, + opt->m, opt->fc, + opt->p, opt->h, + lb, ub, x, minf, &stop, + local_opt, + algorithm == NLOPT_LN_AUGLAG_EQ + || algorithm == NLOPT_LD_AUGLAG_EQ); + if (!opt->local_opt) nlopt_destroy(local_opt); + return ret; + } + + case NLOPT_GN_ISRES: + return isres_minimize(n, f, f_data, + opt->m, opt->fc, + opt->p, opt->h, + lb, ub, x, minf, &stop, + POP(0)); + + default: + return NLOPT_INVALID_ARGS; + } + + return NLOPT_SUCCESS; /* never reached */ +} + +/*********************************************************************/ + +nlopt_result nlopt_optimize_limited(nlopt_opt opt, double *x, double *minf, + int maxeval, double maxtime) +{ + int save_maxeval; + double save_maxtime; + nlopt_result ret; + + if (!opt) return NLOPT_INVALID_ARGS; + + save_maxeval = nlopt_get_maxeval(opt); + save_maxtime = nlopt_get_maxtime(opt); + + /* override opt limits if maxeval and/or maxtime are more stringent */ + if (save_maxeval <= 0 || (maxeval > 0 && maxeval < save_maxeval)) + nlopt_set_maxeval(opt, maxeval); + if (save_maxtime <= 0 || (maxtime > 0 && maxtime < save_maxtime)) + nlopt_set_maxtime(opt, maxtime); + + ret = nlopt_optimize(opt, x, minf); + + nlopt_set_maxeval(opt, save_maxeval); + nlopt_set_maxtime(opt, save_maxtime); + + return ret; +} + +/*********************************************************************/ diff --git a/api/options.c b/api/options.c new file mode 100644 index 0000000..4b0403e --- /dev/null +++ b/api/options.c @@ -0,0 +1,453 @@ +/* Copyright (c) 2007-2009 Massachusetts Institute of Technology + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include +#include +#include + +#include "nlopt-internal.h" +#include "nlopt-util.h" + +#define U(n) ((unsigned int) (n)) + +/*************************************************************************/ + +void nlopt_destroy(nlopt_opt opt) +{ + if (opt) { + free(opt->lb); free(opt->ub); + free(opt->xtol_abs); + free(opt->fc); + free(opt->h); + nlopt_destroy(opt->local_opt); + free(opt->dx); + } +} + +nlopt_opt nlopt_create(nlopt_algorithm algorithm, int n) +{ + nlopt_opt opt; + + if (n < 0 || algorithm < 0 || algorithm >= NLOPT_NUM_ALGORITHMS) + return NULL; + + opt = (nlopt_opt) malloc(sizeof(struct nlopt_opt_s)); + if (opt) { + opt->algorithm = algorithm; + opt->n = n; + opt->f = NULL; opt->f_data = NULL; + + opt->lb = opt->ub = NULL; + opt->m = opt->m_alloc = 0; + opt->fc = NULL; + opt->p = opt->p_alloc = 0; + opt->h = NULL; + + opt->minf_max = -HUGE_VAL; + opt->ftol_rel = opt->ftol_abs = 0; + opt->xtol_rel = 0; opt->xtol_abs = NULL; + opt->maxeval = 0; + opt->maxtime = 0; + + opt->local_opt = NULL; + opt->stochastic_population = 0; + opt->dx = NULL; + + if (n > 0) { + opt->lb = (double *) malloc(sizeof(double) * U(n)); + if (!opt->lb) goto oom; + opt->ub = (double *) malloc(sizeof(double) * U(n)); + if (!opt->ub) goto oom; + opt->xtol_abs = (double *) malloc(sizeof(double) * U(n)); + if (!opt->xtol_abs) goto oom; + nlopt_set_lower_bounds1(opt, -HUGE_VAL); + nlopt_set_upper_bounds1(opt, +HUGE_VAL); + nlopt_set_xtol_abs1(opt, 0.0); + } + } + + return opt; + +oom: + nlopt_destroy(opt); + return NULL; +} + +nlopt_opt nlopt_copy(const nlopt_opt opt) +{ + nlopt_opt nopt = NULL; + if (opt) { + nopt = (nlopt_opt) malloc(sizeof(struct nlopt_opt_s)); + *nopt = *opt; + nopt->lb = nopt->ub = nopt->xtol_abs = NULL; + nopt->fc = nopt->h = NULL; + nopt->m_alloc = nopt->p_alloc = 0; + nopt->local_opt = NULL; + nopt->dx = NULL; + + if (opt->n > 0) { + nopt->lb = (double *) malloc(sizeof(double) * U(opt->n)); + if (!opt->lb) goto oom; + nopt->ub = (double *) malloc(sizeof(double) * U(opt->n)); + if (!opt->ub) goto oom; + nopt->xtol_abs = (double *) malloc(sizeof(double) * U(opt->n)); + if (!opt->xtol_abs) goto oom; + + memcpy(nopt->lb, opt->lb, sizeof(double) * U(opt->n)); + memcpy(nopt->ub, opt->ub, sizeof(double) * U(opt->n)); + memcpy(nopt->xtol_abs, opt->xtol_abs, sizeof(double) * U(opt->n)); + } + + if (opt->m) { + nopt->m_alloc = opt->m; + nopt->fc = (nlopt_constraint *) malloc(sizeof(nlopt_constraint) + * U(opt->m)); + if (!nopt->fc) goto oom; + memcpy(nopt->fc, opt->fc, sizeof(nlopt_constraint) * U(opt->m)); + } + + if (opt->p) { + nopt->p_alloc = opt->p; + nopt->h = (nlopt_constraint *) malloc(sizeof(nlopt_constraint) + * U(opt->p)); + if (!nopt->h) goto oom; + memcpy(nopt->h, opt->h, sizeof(nlopt_constraint) * U(opt->p)); + } + + if (opt->local_opt) { + nopt->local_opt = nlopt_copy(opt->local_opt); + if (!nopt->local_opt) goto oom; + } + + if (opt->dx) { + nopt->dx = (double *) malloc(sizeof(double) * U(opt->n)); + if (!nopt->dx) goto oom; + memcpy(nopt->dx, opt->dx, sizeof(double) * U(opt->n)); + } + } + return nopt; + +oom: + nlopt_destroy(nopt); + return NULL; +} + +/*************************************************************************/ + +nlopt_result nlopt_set_min_objective(nlopt_opt opt, nlopt_func f, void *f_data) +{ + if (opt && f) { + opt->f = f; opt->f_data = f_data; + return NLOPT_SUCCESS; + } + return NLOPT_INVALID_ARGS; +} + +/*************************************************************************/ + +nlopt_result nlopt_set_lower_bounds(nlopt_opt opt, const double *lb) +{ + if (opt) { + memcpy(opt->lb, lb, sizeof(double) * U(opt->n)); + return NLOPT_SUCCESS; + } + return NLOPT_INVALID_ARGS; +} + +nlopt_result nlopt_set_lower_bounds1(nlopt_opt opt, double lb) +{ + if (opt) { + int i; + for (i = 0; i < opt->n; ++i) + opt->lb[i] = lb; + return NLOPT_SUCCESS; + } + return NLOPT_INVALID_ARGS; +} + +nlopt_result nlopt_set_upper_bounds(nlopt_opt opt, const double *ub) +{ + if (opt) { + memcpy(opt->ub, ub, sizeof(double) * U(opt->n)); + return NLOPT_SUCCESS; + } + return NLOPT_INVALID_ARGS; +} + +nlopt_result nlopt_set_upper_bounds1(nlopt_opt opt, double ub) +{ + if (opt) { + int i; + for (i = 0; i < opt->n; ++i) + opt->ub[i] = ub; + return NLOPT_SUCCESS; + } + return NLOPT_INVALID_ARGS; +} + +/*************************************************************************/ + +#define AUGLAG_ALG(a) ((a) == NLOPT_LN_AUGLAG || \ + (a) == NLOPT_LN_AUGLAG_EQ || \ + (a) == NLOPT_LD_AUGLAG || \ + (a) == NLOPT_LD_AUGLAG_EQ) + +nlopt_result nlopt_remove_inequality_constraints(nlopt_opt opt) +{ + free(opt->fc); + opt->fc = NULL; + opt->m = opt->m_alloc = 0; + return NLOPT_SUCCESS; +} + +static nlopt_result add_constraint(int *m, int *m_alloc, + nlopt_constraint **c, + nlopt_func fc, void *fc_data, + double tol) +{ + *m += 1; + if (*m > *m_alloc) { + /* allocate by repeated doubling so that + we end up with O(log m) mallocs rather than O(m). */ + *m_alloc = 2 * (*m); + *c = (nlopt_constraint *) realloc(*c, + sizeof(nlopt_constraint) + * U(*m_alloc)); + if (!*c) { + *m_alloc = *m = 0; + return NLOPT_OUT_OF_MEMORY; + } + } + + (*c)[*m - 1].f = fc; + (*c)[*m - 1].f_data = fc_data; + (*c)[*m - 1].tol = tol; + return NLOPT_SUCCESS; +} + +nlopt_result nlopt_add_inequality_constraint(nlopt_opt opt, + nlopt_func fc, void *fc_data, + double tol) +{ + if (opt && fc && tol >= 0) { + /* nonlinear constraints are only supported with some algorithms */ + if (opt->algorithm != NLOPT_LD_MMA + && opt->algorithm != NLOPT_LN_COBYLA + && !AUGLAG_ALG(opt->algorithm) + && opt->algorithm != NLOPT_GN_ISRES) + return NLOPT_INVALID_ARGS; + + return add_constraint(&opt->m, &opt->m_alloc, &opt->fc, + fc, fc_data, tol); + } + return NLOPT_INVALID_ARGS; +} + +nlopt_result nlopt_remove_equality_constraints(nlopt_opt opt) +{ + free(opt->h); + opt->h = NULL; + opt->p = opt->p_alloc = 0; + return NLOPT_SUCCESS; +} + +nlopt_result nlopt_add_equality_constraint(nlopt_opt opt, + nlopt_func h, void *h_data, + double tol) +{ + if (opt && h && tol >= 0) { + /* equality constraints (h(x) = 0) only via some algorithms */ + if (!AUGLAG_ALG(opt->algorithm) && opt->algorithm != NLOPT_GN_ISRES) + return NLOPT_INVALID_ARGS; + + return add_constraint(&opt->p, &opt->p_alloc, &opt->h, + h, h_data, tol); + } + return NLOPT_INVALID_ARGS; +} + +/*************************************************************************/ + +#define SET(param, T, arg) \ + nlopt_result nlopt_set_##param(nlopt_opt opt, T arg) \ + { \ + if (opt) { \ + opt->arg = arg; \ + return NLOPT_SUCCESS; \ + } \ + return NLOPT_INVALID_ARGS; \ + } + + +#define GET(param, T, arg) T nlopt_get_##param(const nlopt_opt opt) { \ + return opt->arg; \ + } + +#define GETSET(param, T, arg) GET(param, T, arg) SET(param, T, arg) + +GETSET(stopval, double, minf_max) + +GETSET(ftol_rel, double, ftol_rel) +GETSET(ftol_abs, double, ftol_abs) +GETSET(xtol_rel, double, xtol_rel) + +nlopt_result nlopt_set_xtol_abs(nlopt_opt opt, const double *xtol_abs) +{ + if (opt) { + memcpy(opt->xtol_abs, xtol_abs, opt->n & sizeof(double)); + return NLOPT_SUCCESS; + } + return NLOPT_INVALID_ARGS; +} + +nlopt_result nlopt_set_xtol_abs1(nlopt_opt opt, const double xtol_abs) +{ + if (opt) { + int i; + for (i = 0; i < opt->n; ++i) + opt->xtol_abs[i] = xtol_abs; + return NLOPT_SUCCESS; + } + return NLOPT_INVALID_ARGS; +} + +void nlopt_get_xtol_abs(const nlopt_opt opt, double *xtol_abs) +{ + memcpy(xtol_abs, opt->xtol_abs, opt->n & sizeof(double)); +} + +GETSET(maxeval, int, maxeval) + +GETSET(maxtime, double, maxtime) + +/*************************************************************************/ + +nlopt_result nlopt_set_local_optimizer(nlopt_opt opt, + const nlopt_opt local_opt) +{ + if (opt) { + if (local_opt && local_opt->n != opt->n) return NLOPT_INVALID_ARGS; + nlopt_destroy(opt->local_opt); + opt->local_opt = nlopt_copy(local_opt); + if (local_opt) { + if (!opt->local_opt) return NLOPT_OUT_OF_MEMORY; + nlopt_set_lower_bounds(opt->local_opt, opt->lb); + nlopt_set_upper_bounds(opt->local_opt, opt->ub); + nlopt_remove_inequality_constraints(opt->local_opt); + nlopt_remove_equality_constraints(opt->local_opt); + } + return NLOPT_SUCCESS; + } + return NLOPT_INVALID_ARGS; +} + +/*************************************************************************/ + +GETSET(population, int, stochastic_population) + +/*************************************************************************/ + +nlopt_result nlopt_set_initial_step1(nlopt_opt opt, double dx) +{ + int i; + if (!opt || dx == 0) return NLOPT_INVALID_ARGS; + if (!opt->dx && opt->n > 0) { + opt->dx = (double *) malloc(sizeof(double) * U(opt->n)); + if (!opt->dx) return NLOPT_OUT_OF_MEMORY; + } + for (i = 0; i < opt->n; ++i) opt->dx[i] = dx; + return NLOPT_SUCCESS; +} + +nlopt_result nlopt_set_initial_step(nlopt_opt opt, const double *dx) +{ + int i; + if (!opt || !dx) return NLOPT_INVALID_ARGS; + for (i = 0; i < opt->n; ++i) if (dx[i] == 0) return NLOPT_INVALID_ARGS; + if (!opt->dx && nlopt_set_initial_step1(opt, 1) == NLOPT_OUT_OF_MEMORY) + return NLOPT_OUT_OF_MEMORY; + memcpy(opt->dx, dx, sizeof(double) * U(opt->n)); + return NLOPT_SUCCESS; +} + +nlopt_result nlopt_get_initial_step(nlopt_opt opt, const double *x, double *dx) +{ + if (!opt) return NLOPT_INVALID_ARGS; + if (!opt->n) return NLOPT_SUCCESS; + if (!opt->dx) { + nlopt_result ret = nlopt_set_default_initial_step(opt, x); + if (ret != NLOPT_SUCCESS) return ret; + memcpy(dx, opt->dx, sizeof(double) * U(opt->n)); + free(opt->dx); opt->dx = NULL; /* don't save, since x-dependent */ + } + else + memcpy(dx, opt->dx, sizeof(double) * U(opt->n)); + return NLOPT_SUCCESS; +} + +nlopt_result nlopt_set_default_initial_step(nlopt_opt opt, const double *x) +{ + const double *lb, *ub; + int i; + + if (!opt || !x) return NLOPT_INVALID_ARGS; + lb = opt->lb; ub = opt->ub; + + if (!opt->dx && nlopt_set_initial_step1(opt, 1) == NLOPT_OUT_OF_MEMORY) + return NLOPT_OUT_OF_MEMORY; + + /* crude heuristics for initial step size of nonderivative algorithms */ + for (i = 0; i < opt->n; ++i) { + double step = HUGE_VAL; + + if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i]) + && (ub[i] - lb[i]) * 0.25 < step && ub[i] > lb[i]) + step = (ub[i] - lb[i]) * 0.25; + if (!nlopt_isinf(ub[i]) + && ub[i] - x[i] < step && ub[i] > x[i]) + step = (ub[i] - x[i]) * 0.75; + if (!nlopt_isinf(lb[i]) + && x[i] - lb[i] < step && x[i] > lb[i]) + step = (x[i] - lb[i]) * 0.75; + + if (nlopt_isinf(step)) { + if (!nlopt_isinf(ub[i]) + && fabs(ub[i] - x[i]) < fabs(step)) + step = (ub[i] - x[i]) * 1.1; + if (!nlopt_isinf(lb[i]) + && fabs(x[i] - lb[i]) < fabs(step)) + step = (x[i] - lb[i]) * 1.1; + } + if (nlopt_isinf(step) || step == 0) { + step = fabs(x[i]) * 0.25; + } + if (nlopt_isinf(step) || step == 0) + step = 1; + + opt->dx[i] = step; + } + return NLOPT_SUCCESS; +} + +/*************************************************************************/ diff --git a/auglag/auglag.c b/auglag/auglag.c index ad47f98..04731bc 100644 --- a/auglag/auglag.c +++ b/auglag/auglag.c @@ -14,8 +14,8 @@ int auglag_verbose = 1; typedef struct { nlopt_func f; void *f_data; - int m; nlopt_func fc; char *fc_data; ptrdiff_t fc_datum_size; - int p; nlopt_func h; char *h_data; ptrdiff_t h_datum_size; + int m; nlopt_constraint *fc; + int p; nlopt_constraint *h; double rho, *lambda, *mu; double *gradtmp; nlopt_stopping *stop; @@ -35,16 +35,14 @@ static double auglag(int n, const double *x, double *grad, void *data) for (i = 0; i < d->p; ++i) { double h; - h = d->h(n, x, gradtmp, d->h_data + d->h_datum_size * i) - + lambda[i] / rho; + h = d->h[i].f(n, x, gradtmp, d->h[i].f_data) + lambda[i] / rho; L += 0.5 * rho * h*h; if (grad) for (j = 0; j < n; ++j) grad[j] += (rho * h) * gradtmp[j]; } for (i = 0; i < d->m; ++i) { double fc; - fc = d->fc(n, x, gradtmp, d->fc_data + d->fc_datum_size * i) - + mu[i] / rho; + fc = d->fc[i].f(n, x, gradtmp, d->fc[i].f_data) + mu[i] / rho; if (fc > 0) { L += 0.5 * rho * fc*fc; if (grad) for (j = 0; j < n; ++j) @@ -60,30 +58,27 @@ static double auglag(int n, const double *x, double *grad, void *data) /***************************************************************************/ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, - void *fc_data, ptrdiff_t fc_datum_size, - int p, nlopt_func h, - void *h_data, ptrdiff_t h_datum_size, + int m, nlopt_constraint *fc, + int p, nlopt_constraint *h, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, nlopt_stopping *stop, - nlopt_algorithm sub_alg, int sub_has_fc) + nlopt_opt sub_opt, int sub_has_fc) { auglag_data d; nlopt_result ret = NLOPT_SUCCESS; double ICM = HUGE_VAL; - int i; + double *xcur = NULL, fcur; + int i, feasible; /* magic parameters from Birgin & Martinez */ const double tau = 0.5, gam = 10; const double lam_min = -1e20, lam_max = 1e20, mu_max = 1e20; d.f = f; d.f_data = f_data; - d.m = m; d.fc = fc; d.fc_data = (char *) fc_data; - d.fc_datum_size = fc_datum_size; - d.p = p; d.h = h; d.h_data = (char *) h_data; - d.h_datum_size = h_datum_size; + d.m = m; d.fc = fc; + d.p = p; d.h = h; d.stop = stop; /* whether we handle inequality constraints via the augmented @@ -93,8 +88,22 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, else m = 0; - d.gradtmp = (double *) malloc(sizeof(double) * (n + d.p + d.m)); - if (!d.gradtmp) return NLOPT_OUT_OF_MEMORY; + ret = nlopt_set_min_objective(sub_opt, auglag, &d); if (ret<0) return ret; + ret = nlopt_set_lower_bounds(sub_opt, lb); if (ret<0) return ret; + ret = nlopt_set_upper_bounds(sub_opt, ub); if (ret<0) return ret; + ret = nlopt_remove_inequality_constraints(sub_opt); if (ret<0) return ret; + ret = nlopt_remove_equality_constraints(sub_opt); if (ret<0) return ret; + for (i = 0; i < m; ++i) { + ret = nlopt_add_inequality_constraint(sub_opt, fc[i].f, fc[i].f_data, + fc[i].tol); + if (ret < 0) return ret; + } + + xcur = (double *) malloc(sizeof(double) * (2*n + d.p + d.m)); + if (!xcur) return NLOPT_OUT_OF_MEMORY; + memcpy(xcur, x, sizeof(double) * n); + + d.gradtmp = xcur + n; memset(d.gradtmp, 0, sizeof(double) * (n + d.p + d.m)); d.lambda = d.gradtmp + n; d.mu = d.lambda + d.p; @@ -102,48 +111,52 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, /* starting rho suggested by B & M */ { double con2 = 0; - *minf = f(n, x, NULL, f_data); + d.stop->nevals++; + fcur = f(n, xcur, NULL, f_data); + feasible = 1; for (i = 0; i < d.p; ++i) { - double hi = h(n, x, NULL, d.h_data + i*h_datum_size); + double hi = h[i].f(n, xcur, NULL, d.h[i].f_data); + feasible = feasible && fabs(hi) <= h[i].tol; con2 += hi * hi; } for (i = 0; i < d.m; ++i) { - double fci = fc(n, x, NULL, d.fc_data + i*fc_datum_size); + double fci = fc[i].f(n, xcur, NULL, d.fc[i].f_data); + feasible = feasible && fci <= fc[i].tol; if (fci > 0) con2 += fci * fci; } d.rho = MAX(1e-6, MIN(10, 2 * fabs(*minf) / con2)); } - *minf = HUGE_VAL; + if (feasible) + *minf = fcur; + else + *minf = HUGE_VAL; do { double prev_ICM = ICM; - - ret = nlopt_minimize_constrained(sub_alg, n, auglag, &d, - m, fc, fc_data, fc_datum_size, - lb, ub, x, minf, - -HUGE_VAL, - stop->ftol_rel, stop->ftol_abs, - stop->xtol_rel, stop->xtol_abs, - stop->maxeval - stop->nevals, - stop->maxtime - (nlopt_seconds() - - stop->start)); + + ret = nlopt_optimize_limited(sub_opt, xcur, minf, + stop->maxeval - stop->nevals, + stop->maxtime - (nlopt_seconds() + - stop->start)); if (ret < 0) break; - if (nlopt_stop_evals(stop)) {ret = NLOPT_MAXEVAL_REACHED; break;} - if (nlopt_stop_time(stop)) {ret = NLOPT_MAXTIME_REACHED; break;} - *minf = f(n, x, NULL, f_data); + d.stop->nevals++; + fcur = f(n, xcur, NULL, f_data); ICM = 0; + feasible = 1; for (i = 0; i < d.p; ++i) { - double hi = h(n, x, NULL, d.h_data + i*h_datum_size); + double hi = h[i].f(n, xcur, NULL, d.h[i].f_data); double newlam = d.lambda[i] + d.rho * hi; + feasible = feasible && fabs(hi) <= h[i].tol; ICM = MAX(ICM, fabs(hi)); d.lambda[i] = MIN(MAX(lam_min, newlam), lam_max); } for (i = 0; i < d.m; ++i) { - double fci = fc(n, x, NULL, d.fc_data + i*fc_datum_size); + double fci = fc[i].f(n, xcur, NULL, d.fc[i].f_data); double newmu = d.mu[i] + d.rho * fci; + feasible = feasible && fci <= fc[i].tol; ICM = MAX(ICM, fabs(MAX(fci, -d.mu[i] / d.rho))); d.mu[i] = MIN(MAX(0.0, newmu), mu_max); } @@ -158,12 +171,25 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, printf("\n"); } - if (ICM <= stop->htol_abs || ICM <= stop->htol_rel * fabs(*minf)) { - ret = NLOPT_FTOL_REACHED; - break; + /* only check f & x convergence for feasible points... + for this to be effective on active constraints, the user + must set some nonzero tolerance for each constraint */ + if (feasible && fcur < *minf) { + ret = NLOPT_SUCCESS; + if (fcur < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED; + if (nlopt_stop_ftol(stop, fcur,*minf)) ret = NLOPT_FTOL_REACHED; + if (nlopt_stop_x(stop, xcur, x)) ret = NLOPT_XTOL_REACHED; + *minf = fcur; + memcpy(x, xcur, sizeof(double) * n); + if (ret != NLOPT_SUCCESS) break; } + + if (nlopt_stop_evals(stop)) {ret = NLOPT_MAXEVAL_REACHED; break;} + if (nlopt_stop_time(stop)) {ret = NLOPT_MAXTIME_REACHED; break;} + + /* TODO: use some stopping criterion on ICM? */ } while (1); - free(d.gradtmp); + free(xcur); return ret; } diff --git a/auglag/auglag.h b/auglag/auglag.h index 4bbbcfa..8afe4cc 100644 --- a/auglag/auglag.h +++ b/auglag/auglag.h @@ -34,15 +34,13 @@ extern "C" extern int auglag_verbose; nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, - void *fc_data, ptrdiff_t fc_datum_size, - int p, nlopt_func h, - void *h_data, ptrdiff_t h_datum_size, + int m, nlopt_constraint *fc, + int p, nlopt_constraint *h, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, nlopt_stopping *stop, - nlopt_algorithm sub_alg, int sub_has_fc); + nlopt_opt sub_opt, int sub_has_fc); #ifdef __cplusplus } /* extern "C" */ diff --git a/cobyla/cobyla.c b/cobyla/cobyla.c index 388a873..c45b626 100644 --- a/cobyla/cobyla.c +++ b/cobyla/cobyla.c @@ -65,9 +65,7 @@ typedef struct { nlopt_func f; void *f_data; int m_orig; - nlopt_func fc; - char *fc_data; - ptrdiff_t fc_datum_size; + nlopt_constraint *fc; double *xtmp; const double *lb, *ub; } func_wrap_state; @@ -93,7 +91,7 @@ static int func_wrap(int n, int m, double *x, double *f, double *con, *f = s->f(n, xtmp, NULL, s->f_data); for (i = 0; i < s->m_orig; ++i) - con[i] = -s->fc(n, xtmp, NULL, s->fc_data + i * s->fc_datum_size); + con[i] = -s->fc[i].f(n, xtmp, NULL, s->fc[i].f_data); for (j = 0; j < n; ++j) { if (!nlopt_isinf(lb[j])) con[i++] = x[j] - lb[j]; @@ -105,8 +103,7 @@ static int func_wrap(int n, int m, double *x, double *f, double *con, } nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, - void *fc_data_, ptrdiff_t fc_datum_size, + int m, nlopt_constraint *fc, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, @@ -119,7 +116,7 @@ nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data, s.f = f; s.f_data = f_data; s.m_orig = m; - s.fc = fc; s.fc_data = (char*) fc_data_; s.fc_datum_size = fc_datum_size; + s.fc = fc; s.lb = lb; s.ub = ub; s.xtmp = (double *) malloc(sizeof(double) * n); if (!s.xtmp) return NLOPT_OUT_OF_MEMORY; diff --git a/cobyla/cobyla.h b/cobyla/cobyla.h index 4f36e55..2a85e53 100644 --- a/cobyla/cobyla.h +++ b/cobyla/cobyla.h @@ -98,8 +98,7 @@ extern nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, /* NLopt-style interface function */ nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, - void *fc_data_, ptrdiff_t fc_datum_size, + int m, nlopt_constraint *fc, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, diff --git a/isres/isres.c b/isres/isres.c index c41a8ee..243a198 100644 --- a/isres/isres.c +++ b/isres/isres.c @@ -54,10 +54,8 @@ static int key_compare(void *keys_, const void *a_, const void *b_) } nlopt_result isres_minimize(int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, /* fc <= 0 constraints */ - void *fc_data_, ptrdiff_t fc_datum_size, - int p, nlopt_func h, /* h == 0 constraints */ - void *h_data_, ptrdiff_t h_datum_size, + int m, nlopt_constraint *fc, /* fc <= 0 */ + int p, nlopt_constraint *h, /* h == 0 */ const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, @@ -71,8 +69,6 @@ nlopt_result isres_minimize(int n, nlopt_func f, void *f_data, const double SURVIVOR = 1.0/7.0; /* survivor fraction, from paper */ int survivors; nlopt_result ret = NLOPT_SUCCESS; - char *fc_data = (char *) fc_data_; - char *h_data = (char *) h_data_; double *sigmas = 0, *xs; /* population-by-n arrays (row-major) */ double *fval; /* population array of function vals */ double *penalty; /* population array of penalty vals */ @@ -122,20 +118,21 @@ nlopt_result isres_minimize(int n, nlopt_func f, void *f_data, /* evaluate f and constraint violations for whole population */ for (k = 0; k < population; ++k) { + int feasible = 1; double gpenalty; stop->nevals++; fval[k] = f(n, xs + k*n, NULL, f_data); penalty[k] = 0; for (c = 0; c < m; ++c) { /* inequality constraints */ - double gval = fc(n, xs + k*n, NULL, - fc_data + c * fc_datum_size); + double gval = fc[c].f(n, xs + k*n, NULL, fc[c].f_data); + if (gval > fc[c].tol) feasible = 0; if (gval < 0) gval = 0; penalty[k] += gval*gval; } gpenalty = penalty[k]; for (c = m; c < mp; ++c) { /* equality constraints */ - double hval = h(n, xs + k*n, NULL, - h_data + (c-m) * h_datum_size); + double hval = h[c-m].f(n, xs + k*n, NULL, h[c-m].f_data); + if (fabs(hval) > h[c-m].tol) feasible = 0; penalty[k] += hval*hval; } if (penalty[k] > 0) all_feasible = 0; @@ -146,22 +143,24 @@ nlopt_result isres_minimize(int n, nlopt_func f, void *f_data, we decide which solution is the "best" so far? ... need some total order on the solutions? */ - if (penalty[k] <= minf_penalty + if ((penalty[k] <= minf_penalty || feasible) && (fval[k] <= *minf || minf_gpenalty > 0) - && (penalty[k] != minf_penalty || fval[k] != *minf)) { - if (fval[k] < stop->minf_max && penalty[k] == 0) + && ((feasible ? 0 : penalty[k]) != minf_penalty + || fval[k] != *minf)) { + if (fval[k] < stop->minf_max && feasible) ret = NLOPT_MINF_MAX_REACHED; else if (!nlopt_isinf(*minf)) { if (nlopt_stop_f(stop, fval[k], *minf) - && nlopt_stop_f(stop, penalty[k], minf_penalty)) + && nlopt_stop_f(stop, feasible ? 0 : penalty[k], + minf_penalty)) ret = NLOPT_FTOL_REACHED; else if (nlopt_stop_x(stop, xs+k*n, x)) ret = NLOPT_XTOL_REACHED; } memcpy(x, xs+k*n, sizeof(double)*n); *minf = fval[k]; - minf_penalty = penalty[k]; - minf_gpenalty = gpenalty; + minf_penalty = feasible ? 0 : penalty[k]; + minf_gpenalty = feasible ? 0 : gpenalty; if (ret != NLOPT_SUCCESS) goto done; } diff --git a/isres/isres.h b/isres/isres.h index bd506fa..a994aa9 100644 --- a/isres/isres.h +++ b/isres/isres.h @@ -32,10 +32,8 @@ extern "C" #endif /* __cplusplus */ nlopt_result isres_minimize(int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, /* fc <= 0 constraints */ - void *fc_data_, ptrdiff_t fc_datum_size, - int p, nlopt_func h, /* h == 0 constraints */ - void *h_data_, ptrdiff_t h_datum_size, + int m, nlopt_constraint *fc, /* fc <= 0 */ + int p, nlopt_constraint *h, /* h == 0 */ const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, diff --git a/mlsl/mlsl.c b/mlsl/mlsl.c index cf35730..103e9fb 100644 --- a/mlsl/mlsl.c +++ b/mlsl/mlsl.c @@ -94,9 +94,6 @@ typedef struct { nlopt_sobol s; /* sobol data for LDS point generation, or NULL to use pseudo-random numbers */ - nlopt_algorithm local_alg; /* local search algorithm */ - int local_maxeval; /* max # local iterations (0 if unlimited) */ - double R_prefactor, dlm, dbound, gamma; /* parameters of MLSL */ int N; /* number of pts to add per iteration */ } mlsl_data; @@ -281,8 +278,7 @@ nlopt_result mlsl_minimize(int n, nlopt_func f, void *f_data, double *x, /* in: initial guess, out: minimizer */ double *minf, nlopt_stopping *stop, - nlopt_algorithm local_alg, - int local_maxeval, + nlopt_opt local_opt, int Nsamples, /* #samples/iteration (0=default) */ int lds) /* random or low-discrepancy seq. (lds) */ { @@ -304,7 +300,11 @@ nlopt_result mlsl_minimize(int n, nlopt_func f, void *f_data, rb_tree_init(&d.pts, pt_compare); rb_tree_init(&d.lms, lm_compare); d.s = lds ? nlopt_sobol_create((unsigned) n) : NULL; - d.local_alg = local_alg; d.local_maxeval = local_maxeval; + + nlopt_set_min_objective(local_opt, fcount, &d); + nlopt_set_lower_bounds(local_opt, lb); + nlopt_set_upper_bounds(local_opt, ub); + nlopt_set_stopval(local_opt, stop->minf_max); d.gamma = MLSL_GAMMA; @@ -392,16 +392,10 @@ nlopt_result mlsl_minimize(int n, nlopt_func f, void *f_data, lm = (double *) malloc(sizeof(double) * (n+1)); if (!lm) { ret = NLOPT_OUT_OF_MEMORY; goto done; } memcpy(lm+1, p->x, sizeof(double) * n); - lret = nlopt_minimize(local_alg, n, fcount, &d, - lb, ub, lm+1, lm, - stop->minf_max, - stop->ftol_rel, stop->ftol_abs, - stop->xtol_rel, stop->xtol_abs, - local_maxeval > 0 ? - MIN(local_maxeval, - stop->maxeval - stop->nevals) - : stop->maxeval - stop->nevals, - stop->maxtime - (t - stop->start)); + lret = nlopt_optimize_limited(local_opt, lm+1, lm, + stop->maxeval - stop->nevals, + stop->maxtime - + (t - stop->start)); p->minimized = 1; if (lret < 0) { free(lm); ret = lret; goto done; } if (!rb_tree_insert(&d.lms, lm)) { diff --git a/mlsl/mlsl.h b/mlsl/mlsl.h index e8dcf19..c7d39cd 100644 --- a/mlsl/mlsl.h +++ b/mlsl/mlsl.h @@ -36,8 +36,7 @@ nlopt_result mlsl_minimize(int n, nlopt_func f, void *f_data, double *x, /* in: initial guess, out: minimizer */ double *minf, nlopt_stopping *stop, - nlopt_algorithm local_alg, - int local_maxeval, + nlopt_opt local_opt, int Nsamples, /* #samples/iteration (0=default) */ int lds); diff --git a/mma/mma.c b/mma/mma.c index ab12dc1..cd8cfed 100644 --- a/mma/mma.c +++ b/mma/mma.c @@ -146,21 +146,18 @@ static double dual_func(int m, const double *y, double *grad, void *d_) function returns NaN, that constraint becomes inactive. */ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, - void *fc_data_, ptrdiff_t fc_datum_size, + int m, nlopt_constraint *fc, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, nlopt_stopping *stop, - nlopt_algorithm dual_alg, - double dual_tolrel, int dual_maxeval) + nlopt_opt dual_opt) { nlopt_result ret = NLOPT_SUCCESS; double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur; double *dfcdx, *dfcdx_cur; double *fcval, *fcval_cur, *rhoc, *gcval, *y, *dual_lb, *dual_ub; int i, j, k = 0; - char *fc_data = (char *) fc_data_; dual_data dd; int feasible; double infeasibility; @@ -213,7 +210,7 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, feasible = 1; infeasibility = 0; for (i = 0; i < m; ++i) { - fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_datum_size * i); + fcval[i] = fc[i].f(n, x, dfcdx + i*n, fc[i].f_data); feasible = feasible && (fcval[i] <= 0 || isnan(fcval[i])); if (fcval[i] > infeasibility) infeasibility = fcval[i]; } @@ -231,6 +228,10 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, if (!feasible) for (i = 0; i < m; ++i) dual_ub[i] = 1e40; + nlopt_set_min_objective(dual_opt, dual_func, &dd); + nlopt_set_lower_bounds(dual_opt, dual_lb); + nlopt_set_upper_bounds(dual_opt, dual_ub); + while (1) { /* outer iterations */ double fprev = fcur; if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED; @@ -248,25 +249,13 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, /* solve dual problem */ dd.rho = rho; dd.count = 0; - dual_solution: save_verbose = mma_verbose; - mma_verbose = 0; - reti = nlopt_minimize( - dual_alg, m, dual_func, &dd, - dual_lb, dual_ub, y, &min_dual, - -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval, - stop->maxtime - (nlopt_seconds() - stop->start)); + mma_verbose = 0; /* no recursive verbosity */ + reti = nlopt_optimize_limited(dual_opt, y, &min_dual, + 0, + stop->maxtime - (nlopt_seconds() + - stop->start)); mma_verbose = save_verbose; - if (reti == NLOPT_FAILURE && dual_alg != NLOPT_LD_MMA) { - /* LBFGS etc. converge quickly but are sometimes - very finicky if there are any rounding errors in - the gradient, etcetera; if it fails, try again - with MMA called recursively for the dual */ - dual_alg = NLOPT_LD_MMA; - if (mma_verbose) - printf("MMA: switching to recursive MMA for dual\n"); - goto dual_solution; - } if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) { ret = reti; goto done; @@ -287,10 +276,11 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, new_infeasible_constraint = 0; inner_done = dd.gval >= fcur; for (i = 0; i < m; ++i) { - fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n, - fc_data + fc_datum_size * i); + fcval_cur[i] = fc[i].f(n, xcur, dfcdx_cur + i*n, + fc[i].f_data); if (!isnan(fcval_cur[i])) { - feasible_cur = feasible_cur && (fcval_cur[i] <= 0); + feasible_cur = feasible_cur + && (fcval_cur[i] <= fc[i].tol); if (!isnan(fcval[i])) inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]); @@ -317,7 +307,7 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, again (if inner_done), although the constraints may be violated slightly by rounding errors etc. so we must be a little careful about checking feasibility */ - if (feasible_cur) { + if (infeasibility_cur == 0) { if (!feasible) /* reset upper bounds to infinity */ for (i = 0; i < m; ++i) dual_ub[i] = HUGE_VAL; feasible = 1; diff --git a/mma/mma.h b/mma/mma.h index 8148d11..b626419 100644 --- a/mma/mma.h +++ b/mma/mma.h @@ -34,14 +34,12 @@ extern "C" extern int mma_verbose; nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, - int m, nlopt_func fc, - void *fc_data_, ptrdiff_t fc_datum_size, + int m, nlopt_constraint *fc, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, nlopt_stopping *stop, - nlopt_algorithm dual_alg, - double dual_tolrel, int dual_maxeval); + nlopt_opt dual_opt); #ifdef __cplusplus } /* extern "C" */ diff --git a/util/Makefile.am b/util/Makefile.am index 41b501f..04ec0c9 100644 --- a/util/Makefile.am +++ b/util/Makefile.am @@ -1,3 +1,5 @@ +AM_CPPFLAGS = -I$(top_srcdir)/api + noinst_LTLIBRARIES = libutil.la libutil_la_SOURCES = mt19937ar.c sobolseq.c soboldata.h timer.c stop.c nlopt-util.h redblack.c redblack.h qsort_r.c diff --git a/util/nlopt-util.h b/util/nlopt-util.h index a6b9744..e52710b 100644 --- a/util/nlopt-util.h +++ b/util/nlopt-util.h @@ -27,6 +27,8 @@ #include #include "config.h" +#include "nlopt.h" + /* workaround for Solaris + gcc 3.4.x bug (see configure.ac) */ #if defined(__GNUC__) && defined(REPLACEMENT_HUGE_VAL) # undef HUGE_VAL @@ -71,7 +73,6 @@ typedef struct { double ftol_abs; double xtol_rel; const double *xtol_abs; - double htol_rel, htol_abs; int nevals, maxeval; double maxtime, start; } nlopt_stopping; @@ -88,6 +89,22 @@ extern int nlopt_stop_evals(const nlopt_stopping *stop); extern int nlopt_stop_time(const nlopt_stopping *stop); extern int nlopt_stop_evalstime(const nlopt_stopping *stop); +/* for local optimizations, temporarily setting eval/time limits */ +extern nlopt_result nlopt_optimize_limited(nlopt_opt opt, + double *x, double *minf, + int maxevals, double maxtime); + +/* data structure for nonlinear inequality or equality constraint + (f <= 0 or f = 0, respectively). tol (>= 0) is a tolerance + that is used for stopping criteria -- the point is considered + "feasible" for purposes of stopping of the constraint is violated + by at most tol. */ +typedef struct { + nlopt_func f; + void *f_data; + double tol; +} nlopt_constraint; + #ifdef __cplusplus } /* extern "C" */ #endif /* __cplusplus */ -- 2.30.2