libs/math/tools/hypergeometric_1F1_error_plot.cpp

   1 //  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 #define BOOST_ENABLE_ASSERT_HANDLER
   7 #define BOOST_MATH_MAX_SERIES_ITERATION_POLICY INT_MAX
   8 // for consistent behaviour across compilers/platforms:
   9 #define BOOST_MATH_PROMOTE_DOUBLE_POLICY false
  10 // overflow to infinity is OK, we treat these as zero error as long as the sign is correct!
  11 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
  12
  13 #include <iostream>
  14 #include <ctime>
  15 #include <boost/multiprecision/mpfr.hpp>
  16 #include <boost/math/special_functions/hypergeometric_1F1.hpp>
  17 #include <boost/math/special_functions/hypergeometric_pFq.hpp>
  18 #include <boost/math/special_functions/relative_difference.hpp>
  19
  20 #include <boost/random.hpp>
  21 #include <set>
  22 #include <fstream>
  23 #include <boost/iostreams/tee.hpp>
  24 #include <boost/iostreams/stream.hpp>
  25
  26 typedef double test_type;
  27 using boost::multiprecision::mpfr_float;
  28
  29 namespace boost {
  30    //
  31    // We convert assertions into exceptions, so we can log them and continue:
  32    //
  33    void assertion_failed(char const * expr, char const *, char const * file, long line)
  34    {
  35       std::ostringstream oss;
  36       oss << file << ":" << line << " Assertion failed: " << expr;
  37       throw std::runtime_error(oss.str());
  38    }
  39
  40 }
  41
  42 void print_value(double x, std::ostream& os = std::cout)
  43 {
  44    int e;
  45    double m = std::frexp(x, &e);
  46    m = std::ldexp(m, 54);
  47    e -= 54;
  48    boost::int64_t val = (boost::int64_t)m;
  49    BOOST_ASSERT(std::ldexp((double)val, e) == x);
  50    os << "std::ldexp((double)" << val << ", " << e << ")";
  51 }
  52
  53
  54 void print_row(double a, double b, double z, mpfr_float result, std::ostream& os = std::cout)
  55 {
  56    os << "  {{ ";
  57    print_value(a, os);
  58    os << ", ";
  59    print_value(b, os);
  60    os << ", ";
  61    print_value(z, os);
  62    os << ", SC_(" << std::setprecision(45) << result << ") }}" << std::endl;
  63 }
  64
  65 struct error_data
  66 {
  67    error_data(double a, double b, double z, boost::intmax_t e)
  68       : a(a), b(b), z(z), error(e) {}
  69    double a, b, z;
  70    boost::intmax_t error;
  71    bool operator<(const error_data& other)const
  72    {
  73       return error < other.error;
  74    }
  75 };
  76
  77 int main()
  78 {
  79    try {
  80       test_type max_a, max_b, max_z, min_a, min_b, min_z;
  81
  82       unsigned number_of_samples;
  83
  84       std::ofstream log_stream, incalculable_stream, unevaluated_stream, bins_stream;
  85       std::string basename;
  86
  87       std::cout << "Enter range for a: ";
  88       std::cin >> min_a >> max_a;
  89       std::cout << "Enter range for b: ";
  90       std::cin >> min_b >> max_b;
  91       std::cout << "Enter range for z: ";
  92       std::cin >> min_z >> max_z;
  93       std::cout << "Enter number of samples: ";
  94       std::cin >> number_of_samples;
  95       std::cout << "Enter basename for log files: ";
  96       std::cin >> basename;
  97
  98       typedef boost::iostreams::tee_device<std::ostream, std::ostream> tee_sink;
  99       typedef boost::iostreams::stream<tee_sink> tee_stream;
 100
 101       log_stream.open((basename + ".log").c_str());
 102       tee_stream tee_log(tee_sink(std::cout, log_stream));
 103       incalculable_stream.open((basename + "_incalculable.log").c_str());
 104       unevaluated_stream.open((basename + "_unevaluated.log").c_str());
 105       bins_stream.open((basename + "_bins.csv").c_str());
 106       tee_stream tee_bins(tee_sink(std::cout, bins_stream));
 107
 108       boost::random::mt19937 gen(std::time(0));
 109       boost::random::uniform_real_distribution<test_type> a_dist(min_a, max_a);
 110       boost::random::uniform_real_distribution<test_type> b_dist(min_b, max_b);
 111       boost::random::uniform_real_distribution<test_type> z_dist(min_z, max_z);
 112
 113       std::multiset<error_data> errors;
 114       std::map<std::pair<int, int>, int> bins;
 115
 116       unsigned incalculable = 0;
 117       unsigned evaluation_errors = 0;
 118       test_type max_error = 0;
 119
 120       do
 121       {
 122          test_type a = a_dist(gen);
 123          test_type b = b_dist(gen);
 124          test_type z = z_dist(gen);
 125          test_type found, expected;
 126          mpfr_float mp_expected;
 127
 128          try {
 129             mp_expected = boost::math::hypergeometric_pFq_precision({ mpfr_float(a) }, { mpfr_float(b) }, mpfr_float(z), 25, 200.0);
 130             expected = (test_type)mp_expected;
 131          }
 132          catch (const std::exception&)
 133          {
 134             // Unable to compute reference value:
 135             ++incalculable;
 136             tee_log << "Unable to compute reference value in reasonable time: " << std::endl;
 137             print_row(a, b, z, mpfr_float(0), tee_log);
 138             incalculable_stream << std::setprecision(6) << std::scientific << a << "," << b << "," << z << "\n";
 139             continue;
 140          }
 141          try
 142          {
 143             found = boost::math::hypergeometric_1F1(a, b, z);
 144          }
 145          catch (const std::exception&)
 146          {
 147             ++evaluation_errors;
 148             --number_of_samples;
 149             log_stream << "Unexpected exception calculating value: " << std::endl;
 150             print_row(a, b, z, mp_expected, log_stream);
 151             unevaluated_stream << std::setprecision(6) << std::scientific << a << "," << b << "," << z << "\n";
 152             continue;
 153          }
 154          test_type err = boost::math::epsilon_difference(found, expected);
 155          if (err > max_error)
 156          {
 157             tee_log << "New maximum error is: " << err << std::endl;
 158             print_row(a, b, z, mp_expected, tee_log);
 159             max_error = err;
 160          }
 161          try {
 162             errors.insert(error_data(a, b, z, boost::math::lltrunc(err)));
 163          }
 164          catch (...)
 165          {
 166             errors.insert(error_data(a, b, z, INT_MAX));
 167          }
 168          --number_of_samples;
 169          if (number_of_samples % 500 == 0)
 170             std::cout << number_of_samples << " samples to go" << std::endl;
 171       } while (number_of_samples);
 172
 173       tee_log << "Max error found was: " << max_error << std::endl;
 174
 175       unsigned current_bin = 0;
 176       unsigned lim = 1;
 177       unsigned old_lim = 0;
 178
 179       while (errors.size())
 180       {
 181          old_lim = lim;
 182          lim *= 2;
 183          //std::cout << "Enter upper limit for bin " << current_bin << ": ";
 184          //std::cin >> lim;
 185          auto p = errors.upper_bound(error_data(0, 0, 0, lim));
 186          int bin_count = std::distance(errors.begin(), p);
 187          if (bin_count)
 188          {
 189             std::ofstream os((basename + "_errors_" + std::to_string(current_bin + 1) + ".csv").c_str());
 190             os << "a,b,z,error\n";
 191             bins[std::make_pair(old_lim, lim)] = bin_count;
 192             for (auto pos = errors.begin(); pos != p; ++pos)
 193             {
 194                os << pos->a << "," << pos->b << "," << pos->z << "," << pos->error << "\n";
 195             }
 196             errors.erase(errors.begin(), p);
 197          }
 198          ++current_bin;
 199       }
 200
 201       tee_bins << "Results:\n\n";
 202       tee_bins << "#bin,Range,2^N,Count\n";
 203       int hash = 0;
 204       for (auto p = bins.begin(); p != bins.end(); ++p, ++hash)
 205       {
 206          tee_bins << hash << "," << p->first.first << "-" << p->first.second << "," << hash+1 << "," << p->second << std::endl;
 207       }
 208       if (evaluation_errors)
 209       {
 210          tee_bins << ",Failed,," << evaluation_errors << std::endl;
 211       }
 212       if (incalculable)
 213       {
 214          tee_bins << ",Incalculable,," << incalculable << std::endl;
 215       }
 216    }
 217    catch (const std::exception& e)
 218    {
 219       std::cout << "Terminating with unhandled exception: " << e.what() << std::endl;
 220    }
 221
 222    return 0;
 223 }
 224