more robustness to underflow etc. (fix #36)

diff --git a/api/general.c b/api/general.c
index dfe20102de0d8d9a2a1d8ddb46bd75948b5fee41..dc2fb3818323b232f58ca2358ff4d085fe66abf7 100644 (file)
--- 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 0d679dc0ba6c1dd3678ed48fb3b259cc3d2248aa..289dcdb782d19d04c3a098cebd345d9da429d2bb 100644 (file)
--- 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 69c7693cb2e1597e8af459d6668ec5dca2b2524d..b2ce5febb55fc73a7214c67a2f530c07047d11f9 100644 (file)
--- 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 9945b607583eb4b7d716636e30957a3f834bf065..031778922a48d10dc700cca1bb400dd688a18e05 100644 (file)
--- 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 bb3066b8f3c478340c59fed31ff9f872fb2332eb..ce394d8bd189b9a15ffbef81e2513bfbec8d51f1 100644 (file)
--- 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 cb8a4f9fa1d75bd9ac02b0fe27d75d9a33204288..57dbaccac46399849af98ec6975e90e4020ec9fb 100644 (file)
--- 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 <math.h>
+], [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 <math.h>
 ], [double x = copysign(3.14159, -2.7183);], ok=yes, ok=no)

diff --git a/mma/mma.c b/mma/mma.c
index d7e9725ad5c21e01e4c6380a6c915c6747fc52de..afaa024f3ddbf65dc34c6d6fe22bfc72cf715335 100644 (file)
--- 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 a2f5fd61b428de12d6919b072ba64d91f393fc7c..5df932eb0a911b5c47f63b921cf594b6f0c5f58d 100644 (file)
--- 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,