boost/math/tools/test.hpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_TOOLS_TEST_HPP
   7 #define BOOST_MATH_TOOLS_TEST_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/tools/config.hpp>
  14 #include <boost/math/tools/stats.hpp>
  15 #include <boost/math/special_functions/fpclassify.hpp>
  16 #include <boost/test/test_tools.hpp>
  17 #include <stdexcept>
  18
  19 namespace boost{ namespace math{ namespace tools{
  20
  21 template <class T>
  22 struct test_result
  23 {
  24 private:
  25    boost::math::tools::stats<T> stat;   // Statistics for the test.
  26    unsigned worst_case;                 // Index of the worst case test.
  27 public:
  28    test_result() { worst_case = 0; }
  29    void set_worst(int i){ worst_case = i; }
  30    void add(const T& point){ stat.add(point); }
  31    // accessors:
  32    unsigned worst()const{ return worst_case; }
  33    T min BOOST_PREVENT_MACRO_SUBSTITUTION()const{ return (stat.min)(); }
  34    T max BOOST_PREVENT_MACRO_SUBSTITUTION()const{ return (stat.max)(); }
  35    T total()const{ return stat.total(); }
  36    T mean()const{ return stat.mean(); }
  37    boost::uintmax_t count()const{ return stat.count(); }
  38    T variance()const{ return stat.variance(); }
  39    T variance1()const{ return stat.variance1(); }
  40    T rms()const{ return stat.rms(); }
  41
  42    test_result& operator+=(const test_result& t)
  43    {
  44       if((t.stat.max)() > (stat.max)())
  45          worst_case = t.worst_case;
  46       stat += t.stat;
  47       return *this;
  48    }
  49 };
  50
  51 template <class T>
  52 struct calculate_result_type
  53 {
  54    typedef typename T::value_type row_type;
  55    typedef typename row_type::value_type value_type;
  56 };
  57
  58 template <class T>
  59 T relative_error(T a, T b)
  60 {
  61    BOOST_MATH_STD_USING
  62 #ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  63    //
  64    // If math.h has no long double support we can't rely
  65    // on the math functions generating exponents outside
  66    // the range of a double:
  67    //
  68    T min_val = (std::max)(
  69       tools::min_value<T>(),
  70       static_cast<T>((std::numeric_limits<double>::min)()));
  71    T max_val = (std::min)(
  72       tools::max_value<T>(),
  73       static_cast<T>((std::numeric_limits<double>::max)()));
  74 #else
  75    T min_val = tools::min_value<T>();
  76    T max_val = tools::max_value<T>();
  77 #endif
  78
  79    if((a != 0) && (b != 0))
  80    {
  81       // TODO: use isfinite:
  82       if(fabs(b) >= max_val)
  83       {
  84          if(fabs(a) >= max_val)
  85             return 0;  // one infinity is as good as another!
  86       }
  87       // If the result is denormalised, treat all denorms as equivalent:
  88       if((a < min_val) && (a > 0))
  89          a = min_val;
  90       else if((a > -min_val) && (a < 0))
  91          a = -min_val;
  92       if((b < min_val) && (b > 0))
  93          b = min_val;
  94       else if((b > -min_val) && (b < 0))
  95          b = -min_val;
  96       return (std::max)(fabs((a-b)/a), fabs((a-b)/b));
  97    }
  98
  99    // Handle special case where one or both are zero:
 100    if(min_val == 0)
 101       return fabs(a-b);
 102    if(fabs(a) < min_val)
 103       a = min_val;
 104    if(fabs(b) < min_val)
 105       b = min_val;
 106    return (std::max)(fabs((a-b)/a), fabs((a-b)/b));
 107 }
 108
 109 #if defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__)
 110 template <>
 111 inline double relative_error<double>(double a, double b)
 112 {
 113    BOOST_MATH_STD_USING
 114    //
 115    // On Mac OS X we evaluate "double" functions at "long double" precision,
 116    // but "long double" actually has a very slightly narrower range than "double"!
 117    // Therefore use the range of "long double" as our limits since results outside
 118    // that range may have been truncated to 0 or INF:
 119    //
 120    double min_val = (std::max)((double)tools::min_value<long double>(), tools::min_value<double>());
 121    double max_val = (std::min)((double)tools::max_value<long double>(), tools::max_value<double>());
 122
 123    if((a != 0) && (b != 0))
 124    {
 125       // TODO: use isfinite:
 126       if(b > max_val)
 127       {
 128          if(a > max_val)
 129             return 0;  // one infinity is as good as another!
 130       }
 131       // If the result is denormalised, treat all denorms as equivalent:
 132       if((a < min_val) && (a > 0))
 133          a = min_val;
 134       else if((a > -min_val) && (a < 0))
 135          a = -min_val;
 136       if((b < min_val) && (b > 0))
 137          b = min_val;
 138       else if((b > -min_val) && (b < 0))
 139          b = -min_val;
 140       return (std::max)(fabs((a-b)/a), fabs((a-b)/b));
 141    }
 142
 143    // Handle special case where one or both are zero:
 144    if(min_val == 0)
 145       return fabs(a-b);
 146    if(fabs(a) < min_val)
 147       a = min_val;
 148    if(fabs(b) < min_val)
 149       b = min_val;
 150    return (std::max)(fabs((a-b)/a), fabs((a-b)/b));
 151 }
 152 #endif
 153
 154 template <class T>
 155 void set_output_precision(T)
 156 {
 157 #ifdef BOOST_MSVC
 158 #pragma warning(push)
 159 #pragma warning(disable:4127)
 160 #endif
 161    if(std::numeric_limits<T>::digits10)
 162    {
 163       std::cout << std::setprecision(std::numeric_limits<T>::digits10 + 2);
 164    }
 165 #ifdef BOOST_MSVC
 166 #pragma warning(pop)
 167 #endif
 168 }
 169
 170 template <class Seq>
 171 void print_row(const Seq& row)
 172 {
 173    set_output_precision(row[0]);
 174    for(unsigned i = 0; i < row.size(); ++i)
 175    {
 176       if(i)
 177          std::cout << ", ";
 178       std::cout << row[i];
 179    }
 180    std::cout << std::endl;
 181 }
 182
 183 //
 184 // Function test accepts an matrix of input values (probably a 2D boost::array)
 185 // and calls two functors for each row in the array - one calculates a value
 186 // to test, and one extracts the expected value from the array (or possibly
 187 // calculates it at high precision).  The two functors are usually simple lambda
 188 // expressions.
 189 //
 190 template <class A, class F1, class F2>
 191 test_result<typename calculate_result_type<A>::value_type> test(const A& a, F1 test_func, F2 expect_func)
 192 {
 193    typedef typename A::value_type         row_type;
 194    typedef typename row_type::value_type  value_type;
 195
 196    test_result<value_type> result;
 197
 198    for(unsigned i = 0; i < a.size(); ++i)
 199    {
 200       const row_type& row = a[i];
 201       value_type point;
 202       try
 203       {
 204          point = test_func(row);
 205       }
 206       catch(const std::underflow_error&)
 207       {
 208          point = 0;
 209       }
 210       catch(const std::overflow_error&)
 211       {
 212          point = std::numeric_limits<value_type>::has_infinity ?
 213             std::numeric_limits<value_type>::infinity()
 214             : tools::max_value<value_type>();
 215       }
 216       catch(const std::exception& e)
 217       {
 218          std::cerr << e.what() << std::endl;
 219          print_row(row);
 220          BOOST_ERROR("Unexpected exception.");
 221          // so we don't get further errors:
 222          point = expect_func(row);
 223       }
 224       value_type expected = expect_func(row);
 225       value_type err = relative_error(point, expected);
 226 #ifdef BOOST_INSTRUMENT
 227       if(err != 0)
 228       {
 229          std::cout << row[0] << " " << err;
 230          if(std::numeric_limits<value_type>::is_specialized)
 231          {
 232             std::cout << " (" << err / std::numeric_limits<value_type>::epsilon() << "eps)";
 233          }
 234          std::cout << std::endl;
 235       }
 236 #endif
 237       if(!(boost::math::isfinite)(point) && (boost::math::isfinite)(expected))
 238       {
 239          std::cout << "CAUTION: Found non-finite result, when a finite value was expected at entry " << i << "\n";
 240          std::cout << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl;
 241          print_row(row);
 242          BOOST_ERROR("Unexpected non-finite result");
 243       }
 244       if(err > 0.5)
 245       {
 246          std::cout << "CAUTION: Gross error found at entry " << i << ".\n";
 247          std::cout << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl;
 248          print_row(row);
 249          BOOST_ERROR("Gross error");
 250       }
 251       result.add(err);
 252       if((result.max)() == err)
 253          result.set_worst(i);
 254    }
 255    return result;
 256 }
 257
 258 template <class Real, class A, class F1, class F2>
 259 test_result<Real> test_hetero(const A& a, F1 test_func, F2 expect_func)
 260 {
 261    typedef typename A::value_type         row_type;
 262    typedef Real                          value_type;
 263
 264    test_result<value_type> result;
 265
 266    for(unsigned i = 0; i < a.size(); ++i)
 267    {
 268       const row_type& row = a[i];
 269       value_type point;
 270       try
 271       {
 272          point = test_func(row);
 273       }
 274       catch(const std::underflow_error&)
 275       {
 276          point = 0;
 277       }
 278       catch(const std::overflow_error&)
 279       {
 280          point = std::numeric_limits<value_type>::has_infinity ?
 281             std::numeric_limits<value_type>::infinity()
 282             : tools::max_value<value_type>();
 283       }
 284       catch(const std::exception& e)
 285       {
 286          std::cerr << e.what() << std::endl;
 287          print_row(row);
 288          BOOST_ERROR("Unexpected exception.");
 289          // so we don't get further errors:
 290          point = expect_func(row);
 291       }
 292       value_type expected = expect_func(row);
 293       value_type err = relative_error(point, expected);
 294 #ifdef BOOST_INSTRUMENT
 295       if(err != 0)
 296       {
 297          std::cout << row[0] << " " << err;
 298          if(std::numeric_limits<value_type>::is_specialized)
 299          {
 300             std::cout << " (" << err / std::numeric_limits<value_type>::epsilon() << "eps)";
 301          }
 302          std::cout << std::endl;
 303       }
 304 #endif
 305       if(!(boost::math::isfinite)(point) && (boost::math::isfinite)(expected))
 306       {
 307          std::cout << "CAUTION: Found non-finite result, when a finite value was expected at entry " << i << "\n";
 308          std::cout << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl;
 309          print_row(row);
 310          BOOST_ERROR("Unexpected non-finite result");
 311       }
 312       if(err > 0.5)
 313       {
 314          std::cout << "CAUTION: Gross error found at entry " << i << ".\n";
 315          std::cout << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl;
 316          print_row(row);
 317          BOOST_ERROR("Gross error");
 318       }
 319       result.add(err);
 320       if((result.max)() == err)
 321          result.set_worst(i);
 322    }
 323    return result;
 324 }
 325
 326 } // namespace tools
 327 } // namespace math
 328 } // namespace boost
 329
 330 #endif
 331
 332