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))
  35         return 0;
  36     return (fabs(vnew - vold) < abstol || fabs(vnew - vold) < reltol * (fabs(vnew) + fabs(vold)) * 0.5 || (reltol > 0 && vnew == vold));        /* catch vnew == vold == 0 */
  37 }
  38
  39 int nlopt_stop_ftol(const nlopt_stopping * s, double f, double oldf)
  40 {
  41     return (relstop(oldf, f, s->ftol_rel, s->ftol_abs));
  42 }
  43
  44 int nlopt_stop_f(const nlopt_stopping * s, double f, double oldf)
  45 {
  46     return (f <= s->minf_max || nlopt_stop_ftol(s, f, oldf));
  47 }
  48
  49 int nlopt_stop_x(const nlopt_stopping * s, const double *x, const double *oldx)
  50 {
  51     unsigned i;
  52     for (i = 0; i < s->n; ++i)
  53         if (!relstop(oldx[i], x[i], s->xtol_rel, s->xtol_abs[i]))
  54             return 0;
  55     return 1;
  56 }
  57
  58 int nlopt_stop_dx(const nlopt_stopping * s, const double *x, const double *dx)
  59 {
  60     unsigned i;
  61     for (i = 0; i < s->n; ++i)
  62         if (!relstop(x[i] - dx[i], x[i], s->xtol_rel, s->xtol_abs[i]))
  63             return 0;
  64     return 1;
  65 }
  66
  67 static double sc(double x, double smin, double smax)
  68 {
  69     return smin + x * (smax - smin);
  70 }
  71
  72 /* some of the algorithms rescale x to a unit hypercube, so we need to
  73    scale back before we can compare to the tolerances */
  74 int nlopt_stop_xs(const nlopt_stopping * s, const double *xs, const double *oldxs, const double *scale_min, const double *scale_max)
  75 {
  76     unsigned i;
  77     for (i = 0; i < s->n; ++i)
  78         if (relstop(sc(oldxs[i], scale_min[i], scale_max[i]), sc(xs[i], scale_min[i], scale_max[i]), s->xtol_rel, s->xtol_abs[i]))
  79             return 1;
  80     return 0;
  81 }
  82
  83 int nlopt_stop_evals(const nlopt_stopping * s)
  84 {
  85     return (s->maxeval > 0 && *(s->nevals_p) >= s->maxeval);
  86 }
  87
  88 int nlopt_stop_time_(double start, double maxtime)
  89 {
  90     return (maxtime > 0 && nlopt_seconds() - start >= maxtime);
  91 }
  92
  93 int nlopt_stop_time(const nlopt_stopping * s)
  94 {
  95     return nlopt_stop_time_(s->start, s->maxtime);
  96 }
  97
  98 int nlopt_stop_evalstime(const nlopt_stopping * stop)
  99 {
 100     return nlopt_stop_evals(stop) || nlopt_stop_time(stop);
 101 }
 102
 103 int nlopt_stop_forced(const nlopt_stopping * stop)
 104 {
 105     return stop->force_stop && *(stop->force_stop);
 106 }
 107
 108 unsigned nlopt_count_constraints(unsigned p, const nlopt_constraint * c)
 109 {
 110     unsigned i, count = 0;
 111     for (i = 0; i < p; ++i)
 112         count += c[i].m;
 113     return count;
 114 }
 115
 116 unsigned nlopt_max_constraint_dim(unsigned p, const nlopt_constraint * c)
 117 {
 118     unsigned i, max_dim = 0;
 119     for (i = 0; i < p; ++i)
 120         if (c[i].m > max_dim)
 121             max_dim = c[i].m;
 122     return max_dim;
 123 }
 124
 125 void nlopt_eval_constraint(double *result, double *grad, const nlopt_constraint * c, unsigned n, const double *x)
 126 {
 127     if (c->f)
 128         result[0] = c->f(n, x, grad, c->f_data);
 129     else
 130         c->mf(c->m, result, n, x, grad, c->f_data);
 131 }
 132
 133 char *nlopt_vsprintf(char *p, const char *format, va_list ap)
 134 {
 135     size_t len = strlen(format) + 128;
 136     int ret;
 137
 138     p = (char *) realloc(p, len);
 139     if (!p)
 140         abort();
 141
 142     /* TODO: check HAVE_VSNPRINTF, and fallback to vsprintf otherwise */
 143     while ((ret = vsnprintf(p, len, format, ap)) < 0 || (size_t) ret >= len) {
 144         /* C99 vsnprintf returns the required number of bytes (excluding \0)
 145            if the buffer is too small; older versions (e.g. MS) return -1 */
 146         len = ret >= 0 ? (size_t) (ret + 1) : (len * 3) >> 1;
 147         p = (char *) realloc(p, len);
 148         if (!p)
 149             abort();
 150     }
 151     return p;
 152 }
 153
 154 void nlopt_stop_msg(const nlopt_stopping * s, const char *format, ...)
 155 {
 156     va_list ap;
 157     if (s->stop_msg) {
 158         va_start(ap, format);
 159         *(s->stop_msg) = nlopt_vsprintf(*(s->stop_msg), format, ap);
 160         va_end(ap);
 161     }
 162 }
 163
 164 /*************************************************************************/
 165
 166 int nlopt_isinf(double x)
 167 {
 168     return (fabs(x) >= HUGE_VAL * 0.99)
 169 #if defined(HAVE_ISINF)
 170         || isinf(x)
 171 #else
 172         || (!nlopt_isnan(x) && nlopt_isnan(x - x))
 173 #endif
 174         ;
 175 }
 176
 177 int nlopt_isfinite(double x)
 178 {
 179     return (fabs(x) <= DBL_MAX)
 180 #if defined(HAVE_ISFINITE)
 181         || isfinite(x)
 182 #elif defined(_WIN32)
 183         || _finite(x)
 184 #endif
 185         ;
 186 }
 187
 188 int nlopt_istiny(double x)
 189 {
 190     if (x == 0.0)
 191         return 1;
 192     else {
 193 #if defined(HAVE_FPCLASSIFY)
 194         return fpclassify(x) == FP_SUBNORMAL;
 195 #elif defined(_WIN32)
 196         int c = _fpclass(x);
 197         return c == _FPCLASS_ND || c == _FPCLASS_PD;
 198 #else
 199         return fabs(x) < 2.2250738585072014e-308;       /* assume IEEE 754 double */
 200 #endif
 201     }
 202 }
 203
 204 int nlopt_isnan(double x)
 205 {
 206 #if defined(HAVE_ISNAN)
 207     return isnan(x);
 208 #elif defined(_WIN32)
 209     return _isnan(x);
 210 #else
 211     return (x != x);            /* might fail with aggressive optimization */
 212 #endif
 213 }