1 [section:test_data Graphing, Profiling, and Generating Test Data for Special Functions]
3 The class `test_data` and associated helper functions are designed so that in just
4 a few lines of code you should be able to:
6 * Profile a continued fraction, or infinite series for convergence and accuracy.
7 * Generate csv data from a special function that can be imported into your favorite
8 graphing program (or spreadsheet) for further analysis.
9 * Generate high precision test data.
13 #include <boost/math/tools/test_data.hpp>
16 This is a non-core Boost.Math header that is predominantly used for internal
17 maintenance of the library: as a result the library is located under
18 `libs/math/include_private` and you will need to add that directory to
19 your include path in order to use this feature.
22 namespace boost{ namespace math{ namespace tools{
27 periodic_in_range = 1,
33 struct parameter_info;
36 parameter_info<T> make_random_param(T start_range, T end_range, int n_points);
39 parameter_info<T> make_periodic_param(T start_range, T end_range, int n_points);
42 parameter_info<T> make_power_param(T basis, int start_exponent, int end_exponent);
45 bool get_user_parameter_info(parameter_info<T>& info, const char* param_name);
51 typedef std::vector<T> row_type;
52 typedef row_type value_type;
54 typedef std::set<row_type> container_type;
56 typedef typename container_type::reference reference;
57 typedef typename container_type::const_reference const_reference;
58 typedef typename container_type::iterator iterator;
59 typedef typename container_type::const_iterator const_iterator;
60 typedef typename container_type::difference_type difference_type;
61 typedef typename container_type::size_type size_type;
66 test_data(F func, const parameter_info<T>& arg1);
70 test_data& insert(F func, const parameter_info<T>& arg1);
73 test_data& insert(F func, const parameter_info<T>& arg1,
74 const parameter_info<T>& arg2);
77 test_data& insert(F func, const parameter_info<T>& arg1,
78 const parameter_info<T>& arg2,
79 const parameter_info<T>& arg3);
86 const_iterator begin()const;
87 const_iterator end()const;
88 bool operator==(const test_data& d)const;
89 bool operator!=(const test_data& d)const;
90 void swap(test_data& other);
91 size_type size()const;
92 size_type max_size()const;
95 bool operator < (const test_data& dat)const;
96 bool operator <= (const test_data& dat)const;
97 bool operator > (const test_data& dat)const;
98 bool operator >= (const test_data& dat)const;
101 template <class charT, class traits, class T>
102 std::basic_ostream<charT, traits>& write_csv(
103 std::basic_ostream<charT, traits>& os,
104 const test_data<T>& data);
106 template <class charT, class traits, class T>
107 std::basic_ostream<charT, traits>& write_csv(
108 std::basic_ostream<charT, traits>& os,
109 const test_data<T>& data,
110 const charT* separator);
113 std::ostream& write_code(std::ostream& os,
114 const test_data<T>& data,
121 This tool is best illustrated with the following series of examples.
123 The functionality of test_data is split into the following parts:
125 * A functor that implements the function for which data is being generated:
126 this is the bit you have to write.
127 * One of more parameters that are to be passed to the functor, these are
128 described in fairly abstract terms: give me N points distributed like /this/ etc.
129 * The class test_data, that takes the functor and descriptions of the parameters
130 and computes how ever many output points have been requested, these are stored
131 in a sorted container.
132 * Routines to iterate over the test_data container and output the data in either
133 csv format, or as C++ source code (as a table using Boost.Array).
135 [h5 Example 1: Output Data for Graph Plotting]
137 For example, lets say we want to graph the lgamma function between -3 and 100,
138 one could do this like so:
140 #include <boost/math/tools/test_data.hpp>
141 #include <boost/math/special_functions/gamma.hpp>
145 using namespace boost::math::tools;
147 // create an object to hold the data:
148 test_data<double> data;
150 // insert 500 points at uniform intervals between just after -3 and 100:
151 double (*pf)(double) = boost::math::lgamma;
152 data.insert(pf, make_periodic_param(-3.0 + 0.00001, 100.0, 500));
154 // print out in csv format:
155 write_csv(std::cout, data, ", ");
159 Which, when plotted, results in:
163 [h5 Example 2: Creating Test Data]
165 As a second example, let's suppose we want to create highly accurate test
166 data for a special function. Since many special functions have two or
167 more independent parameters, it's very hard to effectively cover all of
168 the possible parameter space without generating gigabytes of data at
169 great computational expense. A second best approach is to provide the tools
170 by which a user (or the library maintainer) can quickly generate more data
171 on demand to probe the function over a particular domain of interest.
173 In this example we'll generate test data for the beta function using
174 [@http://shoup.net/ntl/doc/RR.txt NTL::RR] at 1000 bit precision.
175 Rather than call our generic
176 version of the beta function, we'll implement a deliberately naive version
177 of the beta function using lgamma, and rely on the high precision of the
178 data type used to get results accurate to at least 128-bit precision. In this
179 way our test data is independent of whatever clever tricks we may wish to
180 use inside the our beta function.
182 To start with then, here's the function object that creates the test data:
184 #include <boost/math/tools/ntl.hpp>
185 #include <boost/math/special_functions/gamma.hpp>
186 #include <boost/math/tools/test_data.hpp>
189 #include <boost/math/tools/test_data.hpp>
191 using namespace boost::math::tools;
193 struct beta_data_generator
195 NTL::RR operator()(NTL::RR a, NTL::RR b)
198 // If we throw a domain error then test_data will
199 // ignore this input point. We'll use this to filter
200 // out all cases where a < b since the beta function
201 // is symmetrical in a and b:
204 throw std::domain_error("");
206 // very naively calculate spots with lgamma:
209 g1 = boost::math::lgamma(a, &s1);
210 g2 = boost::math::lgamma(b, &s2);
211 g3 = boost::math::lgamma(a+b, &s3);
219 To create the data, we'll need to input the domains for a and b
220 for which the function will be tested: the function `get_user_parameter_info`
221 is designed for just that purpose. The start of main will look something like:
223 // Set the precision on RR:
224 NTL::RR::SetPrecision(1000); // bits.
225 NTL::RR::SetOutputPrecision(40); // decimal digits.
227 parameter_info<NTL::RR> arg1, arg2;
228 test_data<NTL::RR> data;
230 std::cout << "Welcome.\n"
231 "This program will generate spot tests for the beta function:\n"
238 // prompt the user for the domain of a and b to test:
239 get_user_parameter_info(arg1, "a");
240 get_user_parameter_info(arg2, "b");
243 data.insert(beta_data_generator(), arg1, arg2);
245 // see if the user want's any more domains tested:
246 std::cout << "Any more data [y/n]?";
247 std::getline(std::cin, line);
248 boost::algorithm::trim(line);
249 cont = (line == "y");
252 [caution At this point one potential stumbling block should be mentioned:
253 test_data<>::insert will create a matrix of test data when there are two
254 or more parameters, so if we have two parameters and we're asked for
255 a thousand points on each, that's a ['million test points in total].
256 Don't say you weren't warned!]
258 There's just one final step now, and that's to write the test data to file:
260 std::cout << "Enter name of test data file [default=beta_data.ipp]";
261 std::getline(std::cin, line);
262 boost::algorithm::trim(line);
264 line = "beta_data.ipp";
265 std::ofstream ofs(line.c_str());
266 write_code(ofs, data, "beta_data");
268 The format of the test data looks something like:
270 #define SC_(x) static_cast<T>(BOOST_JOIN(x, L))
271 static const boost::array<boost::array<T, 3>, 1830>
273 SC_(0.4883005917072296142578125),
274 SC_(0.4883005917072296142578125),
275 SC_(3.245912809500479157065104747353807392371),
276 SC_(3.5808107852935791015625),
277 SC_(0.4883005917072296142578125),
278 SC_(1.007653173802923954909901438393379243537),
279 /* ... lots of rows skipped */
282 The first two values in each row are the input parameters that were passed
283 to our functor and the last value is the return value from the functor.
284 Had our functor returned a __tuple rather than a value, then we would have had
285 one entry for each element in the tuple in addition to the input parameters.
287 The first #define serves two purposes:
289 * It reduces the file sizes considerably: all those `static_cast`'s add up to a lot
290 of bytes otherwise (they are needed to suppress compiler warnings when `T` is
291 narrower than a `long double`).
292 * It provides a useful customisation point: for example if we were testing
293 a user-defined type that has more precision than a `long double` we could change
296 [^#define SC_(x) lexical_cast<T>(BOOST_STRINGIZE(x))]
298 in order to ensure that no truncation of the values occurs prior to conversion
299 to `T`. Note that this isn't used by default as it's rather hard on the compiler
300 when the table is large.
302 [h5 Example 3: Profiling a Continued Fraction for Convergence and Accuracy]
304 Alternatively, lets say we want to profile a continued fraction for
305 convergence and error. As an example, we'll use the continued fraction
306 for the upper incomplete gamma function, the following function object
307 returns the next a[sub N ] and b[sub N ] of the continued fraction
308 each time it's invoked:
311 struct upper_incomplete_gamma_fract
317 typedef std::pair<T,T> result_type;
319 upper_incomplete_gamma_fract(T a1, T z1)
320 : z(z1-a1+1), a(a1), k(0)
324 result_type operator()()
328 return result_type(k * (a - k), z);
332 We want to measure both the relative error, and the rate of convergence
333 of this fraction, so we'll write a functor that returns both as a __tuple:
334 class test_data will unpack the tuple for us, and create one column of data
335 for each element in the tuple (in addition to the input parameters):
337 #include <boost/math/tools/test_data.hpp>
338 #include <boost/math/tools/test.hpp>
339 #include <boost/math/special_functions/gamma.hpp>
340 #include <boost/math/tools/ntl.hpp>
341 #include <boost/math/tools/tuple.hpp>
344 struct profile_gamma_fraction
346 typedef ``__tuple``<T, T> result_type;
348 result_type operator()(T val)
350 using namespace boost::math::tools;
351 // estimate the true value, using arbitary precision
352 // arithmetic and NTL::RR:
354 upper_incomplete_gamma_fract<NTL::RR> f1(rval, rval);
355 NTL::RR true_val = continued_fraction_a(f1, 1000);
357 // Now get the aproximation at double precision, along with the number of
358 // iterations required:
359 boost::uintmax_t iters = std::numeric_limits<boost::uintmax_t>::max();
360 upper_incomplete_gamma_fract<T> f2(val, val);
361 T found_val = continued_fraction_a(f2, std::numeric_limits<T>::digits, iters);
363 // Work out the relative error, as measured in units of epsilon:
364 T err = real_cast<T>(relative_error(true_val, NTL::RR(found_val)) / std::numeric_limits<T>::epsilon());
366 // now just return the results as a tuple:
367 return boost::math::make_tuple(err, iters);
371 Feeding that functor into test_data allows rapid output of csv data,
372 for whatever type `T` we may be interested in:
376 using namespace boost::math::tools;
377 // create an object to hold the data:
378 test_data<double> data;
379 // insert 500 points at uniform intervals between just after 0 and 100:
380 data.insert(profile_gamma_fraction<double>(), make_periodic_param(0.01, 100.0, 100));
381 // print out in csv format:
382 write_csv(std::cout, data, ", ");
386 This time there's no need to plot a graph, the first few rows are:
388 a and z, Error/epsilon, Iterations required
397 So it's pretty clear that this fraction shouldn't be used for small values
402 Most of this tool has been described already in the examples above, we'll
403 just add the following notes on the non-member functions:
406 parameter_info<T> make_random_param(T start_range, T end_range, int n_points);
408 Tells class test_data to test /n_points/ random values in the range
409 [start_range,end_range].
412 parameter_info<T> make_periodic_param(T start_range, T end_range, int n_points);
414 Tells class test_data to test /n_points/ evenly spaced values in the range
415 [start_range,end_range].
418 parameter_info<T> make_power_param(T basis, int start_exponent, int end_exponent);
420 Tells class test_data to test points of the form ['basis + R * 2[super expon]] for each
421 /expon/ in the range [start_exponent, end_exponent], and /R/ a random number in \[0.5, 1\].
424 bool get_user_parameter_info(parameter_info<T>& info, const char* param_name);
426 Prompts the user for the parameter range and form to use.
428 Finally, if we don't want the parameter to be included in the output, we can tell
429 test_data by setting it a "dummy parameter":
431 parameter_info<double> p = make_random_param(2.0, 5.0, 10);
432 p.type |= dummy_param;
434 This is useful when the functor used transforms the parameter in some way
435 before passing it to the function under test, usually the functor will then
436 return both the transformed input and the result in a tuple, so there's no
437 need for the original pseudo-parameter to be included in program output.
439 [endsect][/section:test_data Graphing, Profiling, and Generating Test Data for Special Functions]
442 Copyright 2006 John Maddock and Paul A. Bristow.
443 Distributed under the Boost Software License, Version 1.0.
444 (See accompanying file LICENSE_1_0.txt or copy at
445 http://www.boost.org/LICENSE_1_0.txt).