src/util/stop.c

   1 /* Copyright (c) 2007-2014 Massachusetts Institute of Technology
   2  *
   3  * Permission is hereby granted, free of charge, to any person obtaining
   4  * a copy of this software and associated documentation files (the
   5  * "Software"), to deal in the Software without restriction, including
   6  * without limitation the rights to use, copy, modify, merge, publish,
   7  * distribute, sublicense, and/or sell copies of the Software, and to
   8  * permit persons to whom the Software is furnished to do so, subject to
   9  * the following conditions:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  16  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  17  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  18  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  19  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  20  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  21  */
  22
  23 #include <math.h>
  24 #include <float.h>
  25 #include <string.h>
  26 #include <stdio.h>
  27 #include <stdarg.h>
  28 #include "nlopt-util.h"
  29
  30 /* utility routines to implement the various stopping criteria */
  31
  32 static int relstop(double vold, double vnew, double reltol, double abstol)
  33 {
  34      if (nlopt_isinf(vold)) return 0;
  35      return(fabs(vnew - vold) < abstol
  36             || fabs(vnew - vold) < reltol * (fabs(vnew) + fabs(vold)) * 0.5
  37             || (reltol > 0 && vnew == vold)); /* catch vnew == vold == 0 */
  38 }
  39
  40 int nlopt_stop_ftol(const nlopt_stopping *s, double f, double oldf)
  41 {
  42      return (relstop(oldf, f, s->ftol_rel, s->ftol_abs));
  43 }
  44
  45 int nlopt_stop_f(const nlopt_stopping *s, double f, double oldf)
  46 {
  47      return (f <= s->minf_max || nlopt_stop_ftol(s, f, oldf));
  48 }
  49
  50 int nlopt_stop_x(const nlopt_stopping *s, const double *x, const double *oldx)
  51 {
  52      unsigned i;
  53      for (i = 0; i < s->n; ++i)
  54           if (!relstop(oldx[i], x[i], s->xtol_rel, s->xtol_abs[i]))
  55                return 0;
  56      return 1;
  57 }
  58
  59 int nlopt_stop_dx(const nlopt_stopping *s, const double *x, const double *dx)
  60 {
  61      unsigned i;
  62      for (i = 0; i < s->n; ++i)
  63           if (!relstop(x[i] - dx[i], x[i], s->xtol_rel, s->xtol_abs[i]))
  64                return 0;
  65      return 1;
  66 }
  67
  68 static double sc(double x, double smin, double smax)
  69 {
  70      return smin + x * (smax - smin);
  71 }
  72
  73 /* some of the algorithms rescale x to a unit hypercube, so we need to
  74    scale back before we can compare to the tolerances */
  75 int nlopt_stop_xs(const nlopt_stopping *s,
  76                   const double *xs, const double *oldxs,
  77                   const double *scale_min, const double *scale_max)
  78 {
  79      unsigned i;
  80      for (i = 0; i < s->n; ++i)
  81           if (relstop(sc(oldxs[i], scale_min[i], scale_max[i]),
  82                       sc(xs[i], scale_min[i], scale_max[i]),
  83                       s->xtol_rel, s->xtol_abs[i]))
  84                return 1;
  85      return 0;
  86 }
  87
  88 int nlopt_stop_evals(const nlopt_stopping *s)
  89 {
  90      return (s->maxeval > 0 && *(s->nevals_p) >= s->maxeval);
  91 }
  92
  93 int nlopt_stop_time_(double start, double maxtime)
  94 {
  95      return (maxtime > 0 && nlopt_seconds() - start >= maxtime);
  96 }
  97
  98 int nlopt_stop_time(const nlopt_stopping *s)
  99 {
 100      return nlopt_stop_time_(s->start, s->maxtime);
 101 }
 102
 103 int nlopt_stop_evalstime(const nlopt_stopping *stop)
 104 {
 105      return nlopt_stop_evals(stop) || nlopt_stop_time(stop);
 106 }
 107
 108 int nlopt_stop_forced(const nlopt_stopping *stop)
 109 {
 110      return stop->force_stop && *(stop->force_stop);
 111 }
 112
 113 unsigned nlopt_count_constraints(unsigned p, const nlopt_constraint *c)
 114 {
 115      unsigned i, count = 0;
 116      for (i = 0; i < p; ++i)
 117           count += c[i].m;
 118      return count;
 119 }
 120
 121 unsigned nlopt_max_constraint_dim(unsigned p, const nlopt_constraint *c)
 122 {
 123      unsigned i, max_dim = 0;
 124      for (i = 0; i < p; ++i)
 125           if (c[i].m > max_dim)
 126                max_dim = c[i].m;
 127      return max_dim;
 128 }
 129
 130 void nlopt_eval_constraint(double *result, double *grad,
 131                            const nlopt_constraint *c,
 132                            unsigned n, const double *x)
 133 {
 134      if (c->f)
 135           result[0] = c->f(n, x, grad, c->f_data);
 136      else
 137           c->mf(c->m, result, n, x, grad, c->f_data);
 138 }
 139
 140 char *nlopt_vsprintf(char *p, const char *format, va_list ap)
 141 {
 142     size_t len = strlen(format) + 128;
 143     int ret;
 144
 145     p = (char *) realloc(p, len);
 146     if (!p) abort();
 147
 148     /* TODO: check HAVE_VSNPRINTF, and fallback to vsprintf otherwise */
 149     while ((ret = vsnprintf(p, len, format, ap)) < 0 || (size_t)ret >= len) {
 150         /* C99 vsnprintf returns the required number of bytes (excluding \0)
 151            if the buffer is too small; older versions (e.g. MS) return -1 */
 152         len = ret >= 0 ? (size_t)(ret + 1) : (len*3)>>1;
 153         p = (char *) realloc(p, len);
 154         if (!p) abort();
 155     }
 156     return p;
 157 }
 158
 159 void nlopt_stop_msg(const nlopt_stopping *s, const char *format, ...)
 160 {
 161     va_list ap;
 162     if (s->stop_msg) {
 163         va_start(ap, format);
 164         *(s->stop_msg) = nlopt_vsprintf(*(s->stop_msg), format, ap);
 165         va_end(ap);
 166     }
 167 }
 168
 169 /*************************************************************************/
 170
 171 int nlopt_isinf(double x)
 172 {
 173     return (fabs(x) >= HUGE_VAL * 0.99)
 174 #if defined(HAVE_ISINF)
 175         || isinf(x)
 176 #else
 177         || (!nlopt_isnan(x) && nlopt_isnan(x - x))
 178 #endif
 179     ;
 180 }
 181
 182 int nlopt_isfinite(double x)
 183 {
 184     return (fabs(x) <= DBL_MAX)
 185 #if defined(HAVE_ISFINITE)
 186         || isfinite(x)
 187 #elif defined(_WIN32)
 188         || _finite(x)
 189 #endif
 190     ;
 191 }
 192
 193 int nlopt_istiny(double x)
 194 {
 195     if (x == 0.0)
 196         return 1;
 197     else {
 198 #if defined(HAVE_FPCLASSIFY)
 199         return fpclassify(x) == FP_SUBNORMAL;
 200 #elif defined(_WIN32)
 201         int c = _fpclass(x);
 202         return c == _FPCLASS_ND || c == _FPCLASS_PD;
 203 #else
 204         return fabs(x) < 2.2250738585072014e-308; /* assume IEEE 754 double */
 205 #endif
 206     }
 207 }
 208
 209 int nlopt_isnan(double x)
 210 {
 211 #if defined(HAVE_ISNAN)
 212     return isnan(x);
 213 #elif defined(_WIN32)
 214     return _isnan(x);
 215 #else
 216     return (x != x); /* might fail with aggressive optimization */
 217 #endif
 218 }