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