// (C) Copyright John Maddock 2006. | |
// Use, modification and distribution are subject to the | |
// Boost Software License, Version 1.0. (See accompanying file | |
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
#ifndef BOOST_MATH_TOOLS_TEST_HPP | |
#define BOOST_MATH_TOOLS_TEST_HPP | |
#ifdef _MSC_VER | |
#pragma once | |
#endif | |
#include <boost/math/tools/config.hpp> | |
#include <boost/math/tools/stats.hpp> | |
#include <boost/math/special_functions/fpclassify.hpp> | |
#include <boost/test/test_tools.hpp> | |
#include <stdexcept> | |
namespace boost{ namespace math{ namespace tools{ | |
template <class T> | |
struct test_result | |
{ | |
private: | |
boost::math::tools::stats<T> stat; // Statistics for the test. | |
unsigned worst_case; // Index of the worst case test. | |
public: | |
test_result() { worst_case = 0; } | |
void set_worst(int i){ worst_case = i; } | |
void add(const T& point){ stat.add(point); } | |
// accessors: | |
unsigned worst()const{ return worst_case; } | |
T min BOOST_PREVENT_MACRO_SUBSTITUTION()const{ return (stat.min)(); } | |
T max BOOST_PREVENT_MACRO_SUBSTITUTION()const{ return (stat.max)(); } | |
T total()const{ return stat.total(); } | |
T mean()const{ return stat.mean(); } | |
boost::uintmax_t count()const{ return stat.count(); } | |
T variance()const{ return stat.variance(); } | |
T variance1()const{ return stat.variance1(); } | |
T rms()const{ return stat.rms(); } | |
test_result& operator+=(const test_result& t) | |
{ | |
if((t.stat.max)() > (stat.max)()) | |
worst_case = t.worst_case; | |
stat += t.stat; | |
return *this; | |
} | |
}; | |
template <class T> | |
struct calculate_result_type | |
{ | |
typedef typename T::value_type row_type; | |
typedef typename row_type::value_type value_type; | |
}; | |
template <class T> | |
T relative_error(T a, T b) | |
{ | |
BOOST_MATH_STD_USING | |
#ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS | |
// | |
// If math.h has no long double support we can't rely | |
// on the math functions generating exponents outside | |
// the range of a double: | |
// | |
T min_val = (std::max)( | |
tools::min_value<T>(), | |
static_cast<T>((std::numeric_limits<double>::min)())); | |
T max_val = (std::min)( | |
tools::max_value<T>(), | |
static_cast<T>((std::numeric_limits<double>::max)())); | |
#else | |
T min_val = tools::min_value<T>(); | |
T max_val = tools::max_value<T>(); | |
#endif | |
if((a != 0) && (b != 0)) | |
{ | |
// TODO: use isfinite: | |
if(fabs(b) >= max_val) | |
{ | |
if(fabs(a) >= max_val) | |
return 0; // one infinity is as good as another! | |
} | |
// If the result is denormalised, treat all denorms as equivalent: | |
if((a < min_val) && (a > 0)) | |
a = min_val; | |
else if((a > -min_val) && (a < 0)) | |
a = -min_val; | |
if((b < min_val) && (b > 0)) | |
b = min_val; | |
else if((b > -min_val) && (b < 0)) | |
b = -min_val; | |
return (std::max)(fabs((a-b)/a), fabs((a-b)/b)); | |
} | |
// Handle special case where one or both are zero: | |
if(min_val == 0) | |
return fabs(a-b); | |
if(fabs(a) < min_val) | |
a = min_val; | |
if(fabs(b) < min_val) | |
b = min_val; | |
return (std::max)(fabs((a-b)/a), fabs((a-b)/b)); | |
} | |
#if defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__) | |
template <> | |
inline double relative_error<double>(double a, double b) | |
{ | |
BOOST_MATH_STD_USING | |
// | |
// On Mac OS X we evaluate "double" functions at "long double" precision, | |
// but "long double" actually has a very slightly narrower range than "double"! | |
// Therefore use the range of "long double" as our limits since results outside | |
// that range may have been truncated to 0 or INF: | |
// | |
double min_val = (std::max)((double)tools::min_value<long double>(), tools::min_value<double>()); | |
double max_val = (std::min)((double)tools::max_value<long double>(), tools::max_value<double>()); | |
if((a != 0) && (b != 0)) | |
{ | |
// TODO: use isfinite: | |
if(b > max_val) | |
{ | |
if(a > max_val) | |
return 0; // one infinity is as good as another! | |
} | |
// If the result is denormalised, treat all denorms as equivalent: | |
if((a < min_val) && (a > 0)) | |
a = min_val; | |
else if((a > -min_val) && (a < 0)) | |
a = -min_val; | |
if((b < min_val) && (b > 0)) | |
b = min_val; | |
else if((b > -min_val) && (b < 0)) | |
b = -min_val; | |
return (std::max)(fabs((a-b)/a), fabs((a-b)/b)); | |
} | |
// Handle special case where one or both are zero: | |
if(min_val == 0) | |
return fabs(a-b); | |
if(fabs(a) < min_val) | |
a = min_val; | |
if(fabs(b) < min_val) | |
b = min_val; | |
return (std::max)(fabs((a-b)/a), fabs((a-b)/b)); | |
} | |
#endif | |
template <class T> | |
void set_output_precision(T) | |
{ | |
if(std::numeric_limits<T>::digits10) | |
{ | |
std::cout << std::setprecision(std::numeric_limits<T>::digits10 + 2); | |
} | |
} | |
template <class Seq> | |
void print_row(const Seq& row) | |
{ | |
set_output_precision(row[0]); | |
for(unsigned i = 0; i < row.size(); ++i) | |
{ | |
if(i) | |
std::cout << ", "; | |
std::cout << row[i]; | |
} | |
std::cout << std::endl; | |
} | |
// | |
// Function test accepts an matrix of input values (probably a 2D boost::array) | |
// and calls two functors for each row in the array - one calculates a value | |
// to test, and one extracts the expected value from the array (or possibly | |
// calculates it at high precision). The two functors are usually simple lambda | |
// expressions. | |
// | |
template <class A, class F1, class F2> | |
test_result<typename calculate_result_type<A>::value_type> test(const A& a, F1 test_func, F2 expect_func) | |
{ | |
typedef typename A::value_type row_type; | |
typedef typename row_type::value_type value_type; | |
test_result<value_type> result; | |
for(unsigned i = 0; i < a.size(); ++i) | |
{ | |
const row_type& row = a[i]; | |
value_type point; | |
try | |
{ | |
point = test_func(row); | |
} | |
catch(const std::underflow_error&) | |
{ | |
point = 0; | |
} | |
catch(const std::overflow_error&) | |
{ | |
point = std::numeric_limits<value_type>::has_infinity ? | |
std::numeric_limits<value_type>::infinity() | |
: tools::max_value<value_type>(); | |
} | |
catch(const std::exception& e) | |
{ | |
std::cerr << e.what() << std::endl; | |
print_row(row); | |
BOOST_ERROR("Unexpected exception."); | |
// so we don't get further errors: | |
point = expect_func(row); | |
} | |
value_type expected = expect_func(row); | |
value_type err = relative_error(point, expected); | |
#ifdef BOOST_INSTRUMENT | |
if(err != 0) | |
{ | |
std::cout << row[0] << " " << err; | |
if(std::numeric_limits<value_type>::is_specialized) | |
{ | |
std::cout << " (" << err / std::numeric_limits<value_type>::epsilon() << "eps)"; | |
} | |
std::cout << std::endl; | |
} | |
#endif | |
if(!(boost::math::isfinite)(point) && (boost::math::isfinite)(expected)) | |
{ | |
std::cout << "CAUTION: Found non-finite result, when a finite value was expected at entry " << i << "\n"; | |
std::cout << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl; | |
print_row(row); | |
BOOST_ERROR("Unexpected non-finite result"); | |
} | |
if(err > 0.5) | |
{ | |
std::cout << "CAUTION: Gross error found at entry " << i << ".\n"; | |
std::cout << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl; | |
print_row(row); | |
BOOST_ERROR("Gross error"); | |
} | |
result.add(err); | |
if((result.max)() == err) | |
result.set_worst(i); | |
} | |
return result; | |
} | |
} // namespace tools | |
} // namespace math | |
} // namespace boost | |
#endif | |