123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445 |
- [section:test_data Graphing, Profiling, and Generating Test Data for Special Functions]
- The class `test_data` and associated helper functions are designed so that in just
- a few lines of code you should be able to:
- * Profile a continued fraction, or infinite series for convergence and accuracy.
- * Generate csv data from a special function that can be imported into your favorite
- graphing program (or spreadsheet) for further analysis.
- * Generate high precision test data.
- [h4 Synopsis]
- #include <boost/math/tools/test_data.hpp>
- [important
- This is a non-core Boost.Math header that is predominantly used for internal
- maintenance of the library: as a result the library is located under
- `libs/math/include_private` and you will need to add that directory to
- your include path in order to use this feature.
- ]
- namespace boost{ namespace math{ namespace tools{
- enum parameter_type
- {
- random_in_range = 0,
- periodic_in_range = 1,
- power_series = 2,
- dummy_param = 0x80,
- };
- template <class T>
- struct parameter_info;
- template <class T>
- parameter_info<T> make_random_param(T start_range, T end_range, int n_points);
- template <class T>
- parameter_info<T> make_periodic_param(T start_range, T end_range, int n_points);
- template <class T>
- parameter_info<T> make_power_param(T basis, int start_exponent, int end_exponent);
- template <class T>
- bool get_user_parameter_info(parameter_info<T>& info, const char* param_name);
- template <class T>
- class test_data
- {
- public:
- typedef std::vector<T> row_type;
- typedef row_type value_type;
- private:
- typedef std::set<row_type> container_type;
- public:
- typedef typename container_type::reference reference;
- typedef typename container_type::const_reference const_reference;
- typedef typename container_type::iterator iterator;
- typedef typename container_type::const_iterator const_iterator;
- typedef typename container_type::difference_type difference_type;
- typedef typename container_type::size_type size_type;
- // creation:
- test_data(){}
- template <class F>
- test_data(F func, const parameter_info<T>& arg1);
- // insertion:
- template <class F>
- test_data& insert(F func, const parameter_info<T>& arg1);
- template <class F>
- test_data& insert(F func, const parameter_info<T>& arg1,
- const parameter_info<T>& arg2);
- template <class F>
- test_data& insert(F func, const parameter_info<T>& arg1,
- const parameter_info<T>& arg2,
- const parameter_info<T>& arg3);
- void clear();
- // access:
- iterator begin();
- iterator end();
- const_iterator begin()const;
- const_iterator end()const;
- bool operator==(const test_data& d)const;
- bool operator!=(const test_data& d)const;
- void swap(test_data& other);
- size_type size()const;
- size_type max_size()const;
- bool empty()const;
- bool operator < (const test_data& dat)const;
- bool operator <= (const test_data& dat)const;
- bool operator > (const test_data& dat)const;
- bool operator >= (const test_data& dat)const;
- };
- template <class charT, class traits, class T>
- std::basic_ostream<charT, traits>& write_csv(
- std::basic_ostream<charT, traits>& os,
- const test_data<T>& data);
- template <class charT, class traits, class T>
- std::basic_ostream<charT, traits>& write_csv(
- std::basic_ostream<charT, traits>& os,
- const test_data<T>& data,
- const charT* separator);
- template <class T>
- std::ostream& write_code(std::ostream& os,
- const test_data<T>& data,
- const char* name);
-
- }}} // namespaces
-
- [h4 Description]
- This tool is best illustrated with the following series of examples.
- The functionality of test_data is split into the following parts:
- * A functor that implements the function for which data is being generated:
- this is the bit you have to write.
- * One of more parameters that are to be passed to the functor, these are
- described in fairly abstract terms: give me N points distributed like /this/ etc.
- * The class test_data, that takes the functor and descriptions of the parameters
- and computes how ever many output points have been requested, these are stored
- in a sorted container.
- * Routines to iterate over the test_data container and output the data in either
- csv format, or as C++ source code (as a table using Boost.Array).
- [h5 Example 1: Output Data for Graph Plotting]
- For example, lets say we want to graph the lgamma function between -3 and 100,
- one could do this like so:
- #include <boost/math/tools/test_data.hpp>
- #include <boost/math/special_functions/gamma.hpp>
-
- int main()
- {
- using namespace boost::math::tools;
-
- // create an object to hold the data:
- test_data<double> data;
-
- // insert 500 points at uniform intervals between just after -3 and 100:
- double (*pf)(double) = boost::math::lgamma;
- data.insert(pf, make_periodic_param(-3.0 + 0.00001, 100.0, 500));
-
- // print out in csv format:
- write_csv(std::cout, data, ", ");
- return 0;
- }
-
- Which, when plotted, results in:
- [graph lgamma]
- [h5 Example 2: Creating Test Data]
- As a second example, let's suppose we want to create highly accurate test
- data for a special function. Since many special functions have two or
- more independent parameters, it's very hard to effectively cover all of
- the possible parameter space without generating gigabytes of data at
- great computational expense. A second best approach is to provide the tools
- by which a user (or the library maintainer) can quickly generate more data
- on demand to probe the function over a particular domain of interest.
- In this example we'll generate test data for the beta function using
- [@http://shoup.net/ntl/doc/RR.txt NTL::RR] at 1000 bit precision.
- Rather than call our generic
- version of the beta function, we'll implement a deliberately naive version
- of the beta function using lgamma, and rely on the high precision of the
- data type used to get results accurate to at least 128-bit precision. In this
- way our test data is independent of whatever clever tricks we may wish to
- use inside the our beta function.
- To start with then, here's the function object that creates the test data:
- #include <boost/math/tools/ntl.hpp>
- #include <boost/math/special_functions/gamma.hpp>
- #include <boost/math/tools/test_data.hpp>
- #include <fstream>
- #include <boost/math/tools/test_data.hpp>
- using namespace boost::math::tools;
- struct beta_data_generator
- {
- NTL::RR operator()(NTL::RR a, NTL::RR b)
- {
- //
- // If we throw a domain error then test_data will
- // ignore this input point. We'll use this to filter
- // out all cases where a < b since the beta function
- // is symmetrical in a and b:
- //
- if(a < b)
- throw std::domain_error("");
-
- // very naively calculate spots with lgamma:
- NTL::RR g1, g2, g3;
- int s1, s2, s3;
- g1 = boost::math::lgamma(a, &s1);
- g2 = boost::math::lgamma(b, &s2);
- g3 = boost::math::lgamma(a+b, &s3);
- g1 += g2 - g3;
- g1 = exp(g1);
- g1 *= s1 * s2 * s3;
- return g1;
- }
- };
- To create the data, we'll need to input the domains for a and b
- for which the function will be tested: the function `get_user_parameter_info`
- is designed for just that purpose. The start of main will look something like:
- // Set the precision on RR:
- NTL::RR::SetPrecision(1000); // bits.
- NTL::RR::SetOutputPrecision(40); // decimal digits.
- parameter_info<NTL::RR> arg1, arg2;
- test_data<NTL::RR> data;
- std::cout << "Welcome.\n"
- "This program will generate spot tests for the beta function:\n"
- " beta(a, b)\n\n";
- bool cont;
- std::string line;
- do{
- // prompt the user for the domain of a and b to test:
- get_user_parameter_info(arg1, "a");
- get_user_parameter_info(arg2, "b");
-
- // create the data:
- data.insert(beta_data_generator(), arg1, arg2);
-
- // see if the user want's any more domains tested:
- std::cout << "Any more data [y/n]?";
- std::getline(std::cin, line);
- boost::algorithm::trim(line);
- cont = (line == "y");
- }while(cont);
- [caution At this point one potential stumbling block should be mentioned:
- test_data<>::insert will create a matrix of test data when there are two
- or more parameters, so if we have two parameters and we're asked for
- a thousand points on each, that's a ['million test points in total].
- Don't say you weren't warned!]
- There's just one final step now, and that's to write the test data to file:
- std::cout << "Enter name of test data file [default=beta_data.ipp]";
- std::getline(std::cin, line);
- boost::algorithm::trim(line);
- if(line == "")
- line = "beta_data.ipp";
- std::ofstream ofs(line.c_str());
- write_code(ofs, data, "beta_data");
- The format of the test data looks something like:
- #define SC_(x) static_cast<T>(BOOST_JOIN(x, L))
- static const boost::array<boost::array<T, 3>, 1830>
- beta_med_data = {
- SC_(0.4883005917072296142578125),
- SC_(0.4883005917072296142578125),
- SC_(3.245912809500479157065104747353807392371),
- SC_(3.5808107852935791015625),
- SC_(0.4883005917072296142578125),
- SC_(1.007653173802923954909901438393379243537),
- /* ... lots of rows skipped */
- };
- The first two values in each row are the input parameters that were passed
- to our functor and the last value is the return value from the functor.
- Had our functor returned a __tuple rather than a value, then we would have had
- one entry for each element in the tuple in addition to the input parameters.
- The first #define serves two purposes:
- * It reduces the file sizes considerably: all those `static_cast`'s add up to a lot
- of bytes otherwise (they are needed to suppress compiler warnings when `T` is
- narrower than a `long double`).
- * It provides a useful customisation point: for example if we were testing
- a user-defined type that has more precision than a `long double` we could change
- it to:
- [^#define SC_(x) lexical_cast<T>(BOOST_STRINGIZE(x))]
-
- in order to ensure that no truncation of the values occurs prior to conversion
- to `T`. Note that this isn't used by default as it's rather hard on the compiler
- when the table is large.
- [h5 Example 3: Profiling a Continued Fraction for Convergence and Accuracy]
- Alternatively, lets say we want to profile a continued fraction for
- convergence and error. As an example, we'll use the continued fraction
- for the upper incomplete gamma function, the following function object
- returns the next a[sub N ] and b[sub N ] of the continued fraction
- each time it's invoked:
- template <class T>
- struct upper_incomplete_gamma_fract
- {
- private:
- T z, a;
- int k;
- public:
- typedef std::pair<T,T> result_type;
- upper_incomplete_gamma_fract(T a1, T z1)
- : z(z1-a1+1), a(a1), k(0)
- {
- }
- result_type operator()()
- {
- ++k;
- z += 2;
- return result_type(k * (a - k), z);
- }
- };
- We want to measure both the relative error, and the rate of convergence
- of this fraction, so we'll write a functor that returns both as a __tuple:
- class test_data will unpack the tuple for us, and create one column of data
- for each element in the tuple (in addition to the input parameters):
- #include <boost/math/tools/test_data.hpp>
- #include <boost/math/tools/test.hpp>
- #include <boost/math/special_functions/gamma.hpp>
- #include <boost/math/tools/ntl.hpp>
- #include <boost/math/tools/tuple.hpp>
- template <class T>
- struct profile_gamma_fraction
- {
- typedef ``__tuple``<T, T> result_type;
- result_type operator()(T val)
- {
- using namespace boost::math::tools;
- // estimate the true value, using arbitary precision
- // arithmetic and NTL::RR:
- NTL::RR rval(val);
- upper_incomplete_gamma_fract<NTL::RR> f1(rval, rval);
- NTL::RR true_val = continued_fraction_a(f1, 1000);
- //
- // Now get the aproximation at double precision, along with the number of
- // iterations required:
- boost::uintmax_t iters = std::numeric_limits<boost::uintmax_t>::max();
- upper_incomplete_gamma_fract<T> f2(val, val);
- T found_val = continued_fraction_a(f2, std::numeric_limits<T>::digits, iters);
- //
- // Work out the relative error, as measured in units of epsilon:
- T err = real_cast<T>(relative_error(true_val, NTL::RR(found_val)) / std::numeric_limits<T>::epsilon());
- //
- // now just return the results as a tuple:
- return boost::math::make_tuple(err, iters);
- }
- };
- Feeding that functor into test_data allows rapid output of csv data,
- for whatever type `T` we may be interested in:
-
- int main()
- {
- using namespace boost::math::tools;
- // create an object to hold the data:
- test_data<double> data;
- // insert 500 points at uniform intervals between just after 0 and 100:
- data.insert(profile_gamma_fraction<double>(), make_periodic_param(0.01, 100.0, 100));
- // print out in csv format:
- write_csv(std::cout, data, ", ");
- return 0;
- }
- This time there's no need to plot a graph, the first few rows are:
- a and z, Error/epsilon, Iterations required
-
- 0.01, 9723.14, 4726
- 1.0099, 9.54818, 87
- 2.0098, 3.84777, 40
- 3.0097, 0.728358, 25
- 4.0096, 2.39712, 21
- 5.0095, 0.233263, 16
- So it's pretty clear that this fraction shouldn't be used for small values of /a/ and /z/.
- [h4 reference]
- Most of this tool has been described already in the examples above, we'll
- just add the following notes on the non-member functions:
- template <class T>
- parameter_info<T> make_random_param(T start_range, T end_range, int n_points);
-
- Tells class test_data to test /n_points/ random values in the range
- [start_range,end_range].
- template <class T>
- parameter_info<T> make_periodic_param(T start_range, T end_range, int n_points);
-
- Tells class test_data to test /n_points/ evenly spaced values in the range
- [start_range,end_range].
- template <class T>
- parameter_info<T> make_power_param(T basis, int start_exponent, int end_exponent);
- Tells class test_data to test points of the form ['basis + R * 2[super expon]] for each
- /expon/ in the range [start_exponent, end_exponent], and /R/ a random number in \[0.5, 1\].
- template <class T>
- bool get_user_parameter_info(parameter_info<T>& info, const char* param_name);
-
- Prompts the user for the parameter range and form to use.
- Finally, if we don't want the parameter to be included in the output, we can tell
- test_data by setting it a "dummy parameter":
- parameter_info<double> p = make_random_param(2.0, 5.0, 10);
- p.type |= dummy_param;
-
- This is useful when the functor used transforms the parameter in some way
- before passing it to the function under test, usually the functor will then
- return both the transformed input and the result in a tuple, so there's no
- need for the original pseudo-parameter to be included in program output.
- [endsect] [/section:test_data Graphing, Profiling, and Generating Test Data for Special Functions]
- [/
- Copyright 2006 John Maddock and Paul A. Bristow.
- Distributed under 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).
- ]
|