From e537f38c63b2100c715157f698ffce724ab2bdf2 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 12 May 2015 16:08:23 -0400 Subject: [PATCH] more robustness to underflow etc. (fix #36) --- api/general.c | 27 ++++++++++++ api/optimize.c | 11 +---- api/options.c | 104 +++++++++++++++++++++++++++++++++++----------- bobyqa/bobyqa.c | 8 ++++ cobyla/cobyla.c | 5 +++ configure.ac | 8 ++++ mma/mma.c | 19 ++++----- util/nlopt-util.h | 2 + 8 files changed, 138 insertions(+), 46 deletions(-) diff --git a/api/general.c b/api/general.c index dfe2010..dc2fb38 100644 --- a/api/general.c +++ b/api/general.c @@ -39,6 +39,33 @@ int nlopt_isfinite(double x) { return fabs(x) <= DBL_MAX; } +int nlopt_istiny(double x) +{ +#if defined(HAVE_FPCLASSIFY) + return x == 0.0 || fpclassify(x) == FP_SUBNORMAL; +#elif defined(_WIN32) + if (x == 0.0) + return 1; + else { + int c = _fpclass(x); + return c == _FPCLASS_ND || c == _FPCLASS_PD; + } +#else + return fabs(x) < 2.2250738585072014e-308; /* assume IEEE 754 double */ +#endif +} + +int nlopt_isnan(double x) +{ +#if defined(HAVE_ISNAN) + return isnan(x); +#elif defined(_WIN32) + return _isnan(x); +#else + return x != x; /* might fail with aggressive optimization */ +#endif +} + /*************************************************************************/ void NLOPT_STDCALL nlopt_version(int *major, int *minor, int *bugfix) diff --git a/api/optimize.c b/api/optimize.c index 0d679dc..289dcdb 100644 --- a/api/optimize.c +++ b/api/optimize.c @@ -28,13 +28,6 @@ /*********************************************************************/ -#ifndef HAVE_ISNAN -static int my_isnan(double x) { return x != x; } -# define isnan my_isnan -#endif - -/*********************************************************************/ - #include "praxis.h" #include "direct.h" @@ -74,7 +67,7 @@ static double f_bound(int n, const double *x, void *data_) return HUGE_VAL; f = data->f((unsigned) n, x, NULL, data->f_data); - return (isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f); + return (nlopt_isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f); } static double f_noderiv(int n, const double *x, void *data_) @@ -90,7 +83,7 @@ static double f_direct(int n, const double *x, int *undefined, void *data_) double f; unsigned i, j; f = data->f((unsigned) n, x, NULL, data->f_data); - *undefined = isnan(f) || nlopt_isinf(f); + *undefined = nlopt_isnan(f) || nlopt_isinf(f); if (nlopt_get_force_stop(data)) return f; for (i = 0; i < data->m && !*undefined; ++i) { nlopt_eval_constraint(work, NULL, data->fc+i, (unsigned) n, x); diff --git a/api/options.c b/api/options.c index 69c7693..b2ce5fe 100644 --- a/api/options.c +++ b/api/options.c @@ -29,6 +29,8 @@ #include "nlopt-internal.h" +#define ERR(err, opt, msg) (nlopt_set_errmsg(opt, msg) ? err : err) + /*************************************************************************/ void NLOPT_STDCALL nlopt_destroy(nlopt_opt opt) @@ -222,6 +224,7 @@ nlopt_result NLOPT_STDCALL nlopt_set_precond_min_objective(nlopt_opt opt, void *f_data) { if (opt) { + nlopt_unset_errmsg(opt); if (opt->munge_on_destroy) opt->munge_on_destroy(opt->f_data); opt->f = f; opt->f_data = f_data; opt->pre = pre; opt->maximize = 0; @@ -244,6 +247,7 @@ nlopt_result NLOPT_STDCALL nlopt_set_precond_max_objective(nlopt_opt opt, void *f_data) { if (opt) { + nlopt_unset_errmsg(opt); if (opt->munge_on_destroy) opt->munge_on_destroy(opt->f_data); opt->f = f; opt->f_data = f_data; opt->pre = pre; opt->maximize = 1; @@ -265,8 +269,13 @@ nlopt_result NLOPT_STDCALL nlopt_set_max_objective(nlopt_opt opt, nlopt_result NLOPT_STDCALL nlopt_set_lower_bounds(nlopt_opt opt, const double *lb) { + nlopt_unset_errmsg(opt); if (opt && (opt->n == 0 || lb)) { + int i; memcpy(opt->lb, lb, sizeof(double) * (opt->n)); + for (i = 0; i < opt->n; ++i) + if (opt->lb[i] < opt->ub[i] && nlopt_istiny(opt->ub[i] - opt->lb[i])) + opt->lb[i] = opt->ub[i]; return NLOPT_SUCCESS; } return NLOPT_INVALID_ARGS; @@ -275,10 +284,14 @@ NLOPT_STDCALL nlopt_set_lower_bounds(nlopt_opt opt, const double *lb) nlopt_result NLOPT_STDCALL nlopt_set_lower_bounds1(nlopt_opt opt, double lb) { + nlopt_unset_errmsg(opt); if (opt) { unsigned i; for (i = 0; i < opt->n; ++i) - opt->lb[i] = lb; + opt->lb[i] = lb; { + if (opt->lb[i] < opt->ub[i] && nlopt_istiny(opt->ub[i] - opt->lb[i])) + opt->lb[i] = opt->ub[i]; + } return NLOPT_SUCCESS; } return NLOPT_INVALID_ARGS; @@ -287,6 +300,7 @@ NLOPT_STDCALL nlopt_set_lower_bounds1(nlopt_opt opt, double lb) nlopt_result NLOPT_STDCALL nlopt_get_lower_bounds(const nlopt_opt opt, double *lb) { + nlopt_unset_errmsg(opt); if (opt && (opt->n == 0 || lb)) { memcpy(lb, opt->lb, sizeof(double) * (opt->n)); return NLOPT_SUCCESS; @@ -297,8 +311,13 @@ NLOPT_STDCALL nlopt_get_lower_bounds(const nlopt_opt opt, double *lb) nlopt_result NLOPT_STDCALL nlopt_set_upper_bounds(nlopt_opt opt, const double *ub) { + nlopt_unset_errmsg(opt); if (opt && (opt->n == 0 || ub)) { + int i; memcpy(opt->ub, ub, sizeof(double) * (opt->n)); + for (i = 0; i < opt->n; ++i) + if (opt->lb[i] < opt->ub[i] && nlopt_istiny(opt->ub[i] - opt->lb[i])) + opt->ub[i] = opt->lb[i]; return NLOPT_SUCCESS; } return NLOPT_INVALID_ARGS; @@ -307,10 +326,14 @@ NLOPT_STDCALL nlopt_set_upper_bounds(nlopt_opt opt, const double *ub) nlopt_result NLOPT_STDCALL nlopt_set_upper_bounds1(nlopt_opt opt, double ub) { + nlopt_unset_errmsg(opt); if (opt) { unsigned i; - for (i = 0; i < opt->n; ++i) - opt->ub[i] = ub; + for (i = 0; i < opt->n; ++i) { + opt->ub[i] = ub; + if (opt->lb[i] < opt->ub[i] && nlopt_istiny(opt->ub[i] - opt->lb[i])) + opt->ub[i] = opt->lb[i]; + } return NLOPT_SUCCESS; } return NLOPT_INVALID_ARGS; @@ -319,6 +342,7 @@ NLOPT_STDCALL nlopt_set_upper_bounds1(nlopt_opt opt, double ub) nlopt_result NLOPT_STDCALL nlopt_get_upper_bounds(const nlopt_opt opt, double *ub) { + nlopt_unset_errmsg(opt); if (opt && (opt->n == 0 || ub)) { memcpy(ub, opt->ub, sizeof(double) * (opt->n)); return NLOPT_SUCCESS; @@ -339,6 +363,7 @@ nlopt_result NLOPT_STDCALL nlopt_remove_inequality_constraints(nlopt_opt opt) { unsigned i; + nlopt_unset_errmsg(opt); if (!opt) return NLOPT_INVALID_ARGS; if (opt->munge_on_destroy) { nlopt_munge munge = opt->munge_on_destroy; @@ -353,7 +378,8 @@ NLOPT_STDCALL nlopt_remove_inequality_constraints(nlopt_opt opt) return NLOPT_SUCCESS; } -static nlopt_result add_constraint(unsigned *m, unsigned *m_alloc, +static nlopt_result add_constraint(nlopt_opt opt, + unsigned *m, unsigned *m_alloc, nlopt_constraint **c, unsigned fm, nlopt_func fc, nlopt_mfunc mfc, nlopt_precond pre, @@ -366,8 +392,9 @@ static nlopt_result add_constraint(unsigned *m, unsigned *m_alloc, if ((fc && mfc) || (fc && fm != 1) || (!fc && !mfc)) return NLOPT_INVALID_ARGS; if (tol) - for (i = 0; i < fm; ++i) if (tol[i] < 0) return NLOPT_INVALID_ARGS; - + for (i = 0; i < fm; ++i) if (tol[i] < 0) + return ERR(NLOPT_INVALID_ARGS, opt, "negative constraint tolerance"); + tolcopy = (double *) malloc(sizeof(double) * fm); if (fm && !tolcopy) return NLOPT_OUT_OF_MEMORY; if (tol) @@ -416,12 +443,15 @@ NLOPT_STDCALL nlopt_add_inequality_mconstraint(nlopt_opt opt, unsigned m, const double *tol) { nlopt_result ret; + nlopt_unset_errmsg(opt); if (!m) { /* empty constraints are always ok */ if (opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data); return NLOPT_SUCCESS; } - if (!opt || !inequality_ok(opt->algorithm)) ret = NLOPT_INVALID_ARGS; - else ret = add_constraint(&opt->m, &opt->m_alloc, &opt->fc, + if (!opt) ret = NLOPT_INVALID_ARGS; + else if (!inequality_ok(opt->algorithm)) + ret = ERR(NLOPT_INVALID_ARGS, opt, "invalid algorithm for constraints"); + else ret = add_constraint(opt, &opt->m, &opt->m_alloc, &opt->fc, m, NULL, fc, NULL, fc_data, tol); if (ret < 0 && opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data); @@ -436,8 +466,11 @@ NLOPT_STDCALL nlopt_add_precond_inequality_constraint(nlopt_opt opt, double tol) { nlopt_result ret; - if (!opt || !inequality_ok(opt->algorithm)) ret = NLOPT_INVALID_ARGS; - else ret = add_constraint(&opt->m, &opt->m_alloc, &opt->fc, + nlopt_unset_errmsg(opt); + if (!opt) ret = NLOPT_INVALID_ARGS; + else if (!inequality_ok(opt->algorithm)) + ret = ERR(NLOPT_INVALID_ARGS, opt, "invalid algorithm for constraints"); + else ret = add_constraint(opt, &opt->m, &opt->m_alloc, &opt->fc, 1, fc, NULL, pre, fc_data, &tol); if (ret < 0 && opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data); @@ -457,6 +490,7 @@ nlopt_result NLOPT_STDCALL nlopt_remove_equality_constraints(nlopt_opt opt) { unsigned i; + nlopt_unset_errmsg(opt); if (!opt) return NLOPT_INVALID_ARGS; if (opt->munge_on_destroy) { nlopt_munge munge = opt->munge_on_destroy; @@ -485,14 +519,17 @@ NLOPT_STDCALL nlopt_add_equality_mconstraint(nlopt_opt opt, unsigned m, const double *tol) { nlopt_result ret; + nlopt_unset_errmsg(opt); if (!m) { /* empty constraints are always ok */ if (opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data); return NLOPT_SUCCESS; } - if (!opt || !equality_ok(opt->algorithm) - || nlopt_count_constraints(opt->p, opt->h) + m > opt->n) - ret = NLOPT_INVALID_ARGS; - else ret = add_constraint(&opt->p, &opt->p_alloc, &opt->h, + if (!opt) ret = NLOPT_INVALID_ARGS; + else if (!equality_ok(opt->algorithm)) + ret = ERR(NLOPT_INVALID_ARGS, opt, "invalid algorithm for constraints"); + else if (nlopt_count_constraints(opt->p, opt->h) + m > opt->n) + ret = ERR(NLOPT_INVALID_ARGS, opt, "too many equality constraints"); + else ret = add_constraint(opt, &opt->p, &opt->p_alloc, &opt->h, m, NULL, fc, NULL, fc_data, tol); if (ret < 0 && opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data); @@ -507,10 +544,13 @@ NLOPT_STDCALL nlopt_add_precond_equality_constraint(nlopt_opt opt, double tol) { nlopt_result ret; - if (!opt || !equality_ok(opt->algorithm) - || nlopt_count_constraints(opt->p, opt->h) + 1 > opt->n) - ret = NLOPT_INVALID_ARGS; - else ret = add_constraint(&opt->p, &opt->p_alloc, &opt->h, + nlopt_unset_errmsg(opt); + if (!opt) ret = NLOPT_INVALID_ARGS; + else if (!equality_ok(opt->algorithm)) + ret = ERR(NLOPT_INVALID_ARGS, opt, "invalid algorithm for constraints"); + else if (nlopt_count_constraints(opt->p, opt->h) + 1 > opt->n) + ret = ERR(NLOPT_INVALID_ARGS, opt, "too many equality constraints"); + else ret = add_constraint(opt, &opt->p, &opt->p_alloc, &opt->h, 1, fc, NULL, pre, fc_data, &tol); if (ret < 0 && opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data); @@ -531,6 +571,7 @@ NLOPT_STDCALL nlopt_add_equality_constraint(nlopt_opt opt, nlopt_result NLOPT_STDCALL nlopt_set_##param(nlopt_opt opt, T arg) \ { \ if (opt) { \ + nlopt_unset_errmsg(opt); \ opt->arg = arg; \ return NLOPT_SUCCESS; \ } \ @@ -555,6 +596,7 @@ nlopt_result NLOPT_STDCALL nlopt_set_xtol_abs(nlopt_opt opt, const double *xtol_abs) { if (opt) { + nlopt_unset_errmsg(opt); memcpy(opt->xtol_abs, xtol_abs, opt->n * sizeof(double)); return NLOPT_SUCCESS; } @@ -566,6 +608,7 @@ NLOPT_STDCALL nlopt_set_xtol_abs1(nlopt_opt opt, double xtol_abs) { if (opt) { unsigned i; + nlopt_unset_errmsg(opt); for (i = 0; i < opt->n; ++i) opt->xtol_abs[i] = xtol_abs; return NLOPT_SUCCESS; @@ -590,6 +633,7 @@ nlopt_result NLOPT_STDCALL nlopt_set_force_stop(nlopt_opt opt, int force_stop) { if (opt) { + nlopt_unset_errmsg(opt); opt->force_stop = force_stop; if (opt->force_stop_child) return nlopt_set_force_stop(opt->force_stop_child, force_stop); @@ -615,7 +659,9 @@ NLOPT_STDCALL 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_unset_errmsg(opt); + if (local_opt && local_opt->n != opt->n) + return ERR(NLOPT_INVALID_ARGS, opt, "dimension mismatch in local optimizer"); nlopt_destroy(opt->local_opt); opt->local_opt = nlopt_copy(local_opt); if (local_opt) { @@ -643,7 +689,9 @@ GETSET(vector_storage, unsigned, vector_storage) nlopt_result NLOPT_STDCALL nlopt_set_initial_step1(nlopt_opt opt, double dx) { unsigned i; - if (!opt || dx == 0) return NLOPT_INVALID_ARGS; + if (!opt) return NLOPT_INVALID_ARGS; + nlopt_unset_errmsg(opt); + if (dx == 0) return ERR(NLOPT_INVALID_ARGS, opt, "zero step size"); if (!opt->dx && opt->n > 0) { opt->dx = (double *) malloc(sizeof(double) * (opt->n)); if (!opt->dx) return NLOPT_OUT_OF_MEMORY; @@ -657,11 +705,13 @@ NLOPT_STDCALL nlopt_set_initial_step(nlopt_opt opt, const double *dx) { unsigned i; if (!opt) return NLOPT_INVALID_ARGS; + nlopt_unset_errmsg(opt); if (!dx) { free(opt->dx); opt->dx = NULL; return NLOPT_SUCCESS; } - for (i = 0; i < opt->n; ++i) if (dx[i] == 0) return NLOPT_INVALID_ARGS; + for (i = 0; i < opt->n; ++i) + if (dx[i] == 0) return ERR(NLOPT_INVALID_ARGS, opt, "zero step size"); if (!opt->dx && nlopt_set_initial_step1(opt, 1) == NLOPT_OUT_OF_MEMORY) return NLOPT_OUT_OF_MEMORY; memcpy(opt->dx, dx, sizeof(double) * (opt->n)); @@ -673,6 +723,7 @@ NLOPT_STDCALL nlopt_get_initial_step(const nlopt_opt opt, const double *x, double *dx) { if (!opt) return NLOPT_INVALID_ARGS; + nlopt_unset_errmsg(opt); if (!opt->n) return NLOPT_SUCCESS; if (!opt->dx) { nlopt_opt o = (nlopt_opt) opt; /* discard const temporarily */ @@ -692,6 +743,7 @@ NLOPT_STDCALL nlopt_set_default_initial_step(nlopt_opt opt, const double *x) const double *lb, *ub; unsigned i; + nlopt_unset_errmsg(opt); if (!opt || !x) return NLOPT_INVALID_ARGS; lb = opt->lb; ub = opt->ub; @@ -720,10 +772,10 @@ NLOPT_STDCALL nlopt_set_default_initial_step(nlopt_opt opt, const double *x) && fabs(x[i] - lb[i]) < fabs(step)) step = (x[i] - lb[i]) * 1.1; } - if (nlopt_isinf(step) || step == 0) { + if (nlopt_isinf(step) || nlopt_istiny(step)) { step = x[i]; } - if (nlopt_isinf(step) || step == 0) + if (nlopt_isinf(step) || step == 0.0) step = 1; opt->dx[i] = step; @@ -767,8 +819,10 @@ const char *nlopt_set_errmsg(nlopt_opt opt, const char *format, ...) void nlopt_unset_errmsg(nlopt_opt opt) { - free(opt->errmsg); - opt->errmsg = NULL; + if (opt) { + free(opt->errmsg); + opt->errmsg = NULL; + } } const char *nlopt_get_errmsg(nlopt_opt opt) diff --git a/bobyqa/bobyqa.c b/bobyqa/bobyqa.c index 9945b60..0317789 100644 --- a/bobyqa/bobyqa.c +++ b/bobyqa/bobyqa.c @@ -3096,6 +3096,11 @@ nlopt_result bobyqa(int n, int npt, double *x, equal in all directions */ s = nlopt_compute_rescaling(U(n), dx); if (!s) return NLOPT_OUT_OF_MEMORY; + for (j = 0; j < n; ++j) + if (s[j] == 0 || !nlopt_isfinite(s[j])) { + nlopt_stop_msg(stop, "invalid scaling %g of dimension %d: possible over/underflow?", s[j], j); + ret = NLOPT_INVALID_ARGS; goto done; + } /* this statement must go before goto done, so that --x occurs */ nlopt_rescale(U(n), s, x, x); --x; @@ -3171,6 +3176,7 @@ nlopt_result bobyqa(int n, int npt, double *x, if (npt < n + 2 || npt > (n + 2) * np / 2) { /* Return from BOBYQA because NPT is not in the required interval */ ret = NLOPT_INVALID_ARGS; + nlopt_stop_msg(stop, "invalid number of sampling points"); goto done; } @@ -3216,6 +3222,8 @@ nlopt_result bobyqa(int n, int npt, double *x, /* Return from BOBYQA because one of the differences XU(I)-XL(I)s is less than 2*RHOBEG. */ ret = NLOPT_INVALID_ARGS; + nlopt_stop_msg(stop, "insufficient space between the bounds: %g - %g < %g", + xu[j], xl[j], rhobeg+rhobeg); goto done; } jsl = isl + j - 1; diff --git a/cobyla/cobyla.c b/cobyla/cobyla.c index bb3066b..ce394d8 100644 --- a/cobyla/cobyla.c +++ b/cobyla/cobyla.c @@ -198,6 +198,11 @@ nlopt_result cobyla_minimize(unsigned n, nlopt_func f, void *f_data, s.scale = nlopt_compute_rescaling(n, dx); if (!s.scale) { ret = NLOPT_OUT_OF_MEMORY; goto done; } + for (j = 0; j < n; ++j) + if (s.scale[j] == 0 || !nlopt_isfinite(s.scale[j])) { + nlopt_stop_msg(stop, "invalid scaling %g of dimension %d: possible over/underflow?", s.scale[j], j); + ret = NLOPT_INVALID_ARGS; goto done; + } s.lb = nlopt_new_rescaled(n, s.scale, lb); if (!s.lb) { ret = NLOPT_OUT_OF_MEMORY; goto done; } diff --git a/configure.ac b/configure.ac index cb8a4f9..57dbacc 100644 --- a/configure.ac +++ b/configure.ac @@ -86,6 +86,14 @@ if test "$ok" = "yes"; then fi AC_MSG_RESULT(${ok}) +AC_MSG_CHECKING([for fpclassify]) +AC_TRY_LINK([#include +], [if (!fpclassify(3.14159)) fpclassify(2.7183);], ok=yes, ok=no) +if test "$ok" = "yes"; then + AC_DEFINE(HAVE_FPCLASSIFY,1,[Define if the fpclassify() function/macro is available.]) +fi +AC_MSG_RESULT(${ok}) + AC_MSG_CHECKING([for copysign]) AC_TRY_LINK([#include ], [double x = copysign(3.14159, -2.7183);], ok=yes, ok=no) diff --git a/mma/mma.c b/mma/mma.c index d7e9725..afaa024 100644 --- a/mma/mma.c +++ b/mma/mma.c @@ -33,11 +33,6 @@ unsigned mma_verbose = 0; /* > 0 for verbose output */ #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) -#ifndef HAVE_ISNAN -static int my_isnan(double x) { return x != x; } -# define isnan my_isnan -#endif - /* magic minimum value for rho in MMA ... the 2002 paper says it should be a "fixed, strictly positive `small' number, e.g. 1e-5" ... grrr, I hate these magic numbers, which seem like they @@ -80,7 +75,7 @@ static double dual_func(unsigned m, const double *y, double *grad, void *d_) val = d->gval = fval; d->wval = 0; for (i = 0; i < m; ++i) - val += y[i] * (gcval[i] = isnan(fcval[i]) ? 0 : fcval[i]); + val += y[i] * (gcval[i] = nlopt_isnan(fcval[i]) ? 0 : fcval[i]); for (j = 0; j < n; ++j) { double u, v, dx, denominv, c, sigma2, dx2; @@ -105,7 +100,7 @@ static double dual_func(unsigned m, const double *y, double *grad, void *d_) u = dfdx[j]; v = fabs(dfdx[j]) * sigma[j] + 0.5 * rho; - for (i = 0; i < m; ++i) if (!isnan(fcval[i])) { + for (i = 0; i < m; ++i) if (!nlopt_isnan(fcval[i])) { u += dfcdx[i*n + j] * y[i]; v += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i]; } @@ -128,7 +123,7 @@ static double dual_func(unsigned m, const double *y, double *grad, void *d_) d->gval += (dfdx[j] * c + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * dx2) * denominv; d->wval += 0.5 * dx2 * denominv; - for (i = 0; i < m; ++i) if (!isnan(fcval[i])) + for (i = 0; i < m; ++i) if (!nlopt_isnan(fcval[i])) gcval[i] += (dfcdx[i*n+j] * c + (fabs(dfcdx[i*n+j])*sigma[j] + 0.5*rhoc[i]) * dx2) * denominv; @@ -226,7 +221,7 @@ nlopt_result mma_minimize(unsigned n, nlopt_func f, void *f_data, if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; } } for (i = 0; i < m; ++i) { - feasible = feasible && (fcval[i] <= 0 || isnan(fcval[i])); + feasible = feasible && (fcval[i] <= 0 || nlopt_isnan(fcval[i])); if (fcval[i] > infeasibility) infeasibility = fcval[i]; } /* For non-feasible initial points, set a finite (large) @@ -308,10 +303,10 @@ nlopt_result mma_minimize(unsigned n, nlopt_func f, void *f_data, for (i = ifc = 0; ifc < mfc; ++ifc) { unsigned i0 = i, inext = i + fc[ifc].m; for (; i < inext; ++i) - if (!isnan(fcval_cur[i])) { + if (!nlopt_isnan(fcval_cur[i])) { feasible_cur = feasible_cur && (fcval_cur[i] <= fc[ifc].tol[i-i0]); - if (!isnan(fcval[i])) + if (!nlopt_isnan(fcval[i])) inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]); else if (fcval_cur[i] > 0) @@ -359,7 +354,7 @@ nlopt_result mma_minimize(unsigned n, nlopt_func f, void *f_data, if (fcur > dd.gval) rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval)); for (i = 0; i < m; ++i) - if (!isnan(fcval_cur[i]) && fcval_cur[i] > dd.gcval[i]) + if (!nlopt_isnan(fcval_cur[i]) && fcval_cur[i] > dd.gcval[i]) rhoc[i] = MIN(10*rhoc[i], 1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i]) diff --git a/util/nlopt-util.h b/util/nlopt-util.h index a2f5fd6..5df932e 100644 --- a/util/nlopt-util.h +++ b/util/nlopt-util.h @@ -48,6 +48,8 @@ extern "C" int nlopt_isinf(double x); int nlopt_isfinite(double x); +int nlopt_istiny(double x); +int nlopt_isnan(double x); /* re-entrant qsort */ extern void nlopt_qsort_r(void *base_, size_t nmemb, size_t size, void *thunk, -- 2.30.2