boost/libs/math/tools/hypergeometric_1F1_error_plot.cpp

//  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)

#define BOOST_ENABLE_ASSERT_HANDLER
#define BOOST_MATH_MAX_SERIES_ITERATION_POLICY INT_MAX
// for consistent behaviour across compilers/platforms:
#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false
// overflow to infinity is OK, we treat these as zero error as long as the sign is correct!
#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error

#include <iostream>
#include <ctime>
#include <boost/multiprecision/mpfr.hpp>
#include <boost/math/special_functions/hypergeometric_1F1.hpp>
#include <boost/math/special_functions/hypergeometric_pFq.hpp>
#include <boost/math/special_functions/relative_difference.hpp>

#include <boost/random.hpp>
#include <set>
#include <fstream>
#include <boost/iostreams/tee.hpp>
#include <boost/iostreams/stream.hpp>

typedef double test_type;
using boost::multiprecision::mpfr_float;

namespace boost {
   //
   // We convert assertions into exceptions, so we can log them and continue:
   //
   void assertion_failed(char const * expr, char const *, char const * file, long line)
   {
      std::ostringstream oss;
      oss << file << ":" << line << " Assertion failed: " << expr;
      throw std::runtime_error(oss.str());
   }

}

void print_value(double x, std::ostream& os = std::cout)
{
   int e;
   double m = std::frexp(x, &e);
   m = std::ldexp(m, 54);
   e -= 54;
   boost::int64_t val = (boost::int64_t)m;
   BOOST_ASSERT(std::ldexp((double)val, e) == x);
   os << "std::ldexp((double)" << val << ", " << e << ")";
}


void print_row(double a, double b, double z, mpfr_float result, std::ostream& os = std::cout)
{
   os << "  {{ ";
   print_value(a, os);
   os << ", ";
   print_value(b, os);
   os << ", ";
   print_value(z, os);
   os << ", SC_(" << std::setprecision(45) << result << ") }}" << std::endl;
}

struct error_data
{
   error_data(double a, double b, double z, boost::intmax_t e)
      : a(a), b(b), z(z), error(e) {}
   double a, b, z;
   boost::intmax_t error;
   bool operator<(const error_data& other)const
   {
      return error < other.error;
   }
};

int main()
{
   try {
      test_type max_a, max_b, max_z, min_a, min_b, min_z;

      unsigned number_of_samples;

      std::ofstream log_stream, incalculable_stream, unevaluated_stream, bins_stream;
      std::string basename;

      std::cout << "Enter range for a: ";
      std::cin >> min_a >> max_a;
      std::cout << "Enter range for b: ";
      std::cin >> min_b >> max_b;
      std::cout << "Enter range for z: ";
      std::cin >> min_z >> max_z;
      std::cout << "Enter number of samples: ";
      std::cin >> number_of_samples;
      std::cout << "Enter basename for log files: ";
      std::cin >> basename;

      typedef boost::iostreams::tee_device<std::ostream, std::ostream> tee_sink;
      typedef boost::iostreams::stream<tee_sink> tee_stream;

      log_stream.open((basename + ".log").c_str());
      tee_stream tee_log(tee_sink(std::cout, log_stream));
      incalculable_stream.open((basename + "_incalculable.log").c_str());
      unevaluated_stream.open((basename + "_unevaluated.log").c_str());
      bins_stream.open((basename + "_bins.csv").c_str());
      tee_stream tee_bins(tee_sink(std::cout, bins_stream));

      boost::random::mt19937 gen(std::time(0));
      boost::random::uniform_real_distribution<test_type> a_dist(min_a, max_a);
      boost::random::uniform_real_distribution<test_type> b_dist(min_b, max_b);
      boost::random::uniform_real_distribution<test_type> z_dist(min_z, max_z);

      std::multiset<error_data> errors;
      std::map<std::pair<int, int>, int> bins;

      unsigned incalculable = 0;
      unsigned evaluation_errors = 0;
      test_type max_error = 0;

      do
      {
         test_type a = a_dist(gen);
         test_type b = b_dist(gen);
         test_type z = z_dist(gen);
         test_type found, expected;
         mpfr_float mp_expected;

         try {
            mp_expected = boost::math::hypergeometric_pFq_precision({ mpfr_float(a) }, { mpfr_float(b) }, mpfr_float(z), 25, 200.0);
            expected = (test_type)mp_expected;
         }
         catch (const std::exception&)
         {
            // Unable to compute reference value:
            ++incalculable;
            tee_log << "Unable to compute reference value in reasonable time: " << std::endl;
            print_row(a, b, z, mpfr_float(0), tee_log);
            incalculable_stream << std::setprecision(6) << std::scientific << a << "," << b << "," << z << "\n";
            continue;
         }
         try
         {
            found = boost::math::hypergeometric_1F1(a, b, z);
         }
         catch (const std::exception&)
         {
            ++evaluation_errors;
            --number_of_samples;
            log_stream << "Unexpected exception calculating value: " << std::endl;
            print_row(a, b, z, mp_expected, log_stream);
            unevaluated_stream << std::setprecision(6) << std::scientific << a << "," << b << "," << z << "\n";
            continue;
         }
         test_type err = boost::math::epsilon_difference(found, expected);
         if (err > max_error)
         {
            tee_log << "New maximum error is: " << err << std::endl;
            print_row(a, b, z, mp_expected, tee_log);
            max_error = err;
         }
         try {
            errors.insert(error_data(a, b, z, boost::math::lltrunc(err)));
         }
         catch (...)
         {
            errors.insert(error_data(a, b, z, INT_MAX));
         }
         --number_of_samples;
         if (number_of_samples % 500 == 0)
            std::cout << number_of_samples << " samples to go" << std::endl;
      } while (number_of_samples);

      tee_log << "Max error found was: " << max_error << std::endl;

      unsigned current_bin = 0;
      unsigned lim = 1;
      unsigned old_lim = 0;

      while (errors.size())
      {
         old_lim = lim;
         lim *= 2;
         //std::cout << "Enter upper limit for bin " << current_bin << ": ";
         //std::cin >> lim;
         auto p = errors.upper_bound(error_data(0, 0, 0, lim));
         int bin_count = std::distance(errors.begin(), p);
         if (bin_count)
         {
            std::ofstream os((basename + "_errors_" + std::to_string(current_bin + 1) + ".csv").c_str());
            os << "a,b,z,error\n";
            bins[std::make_pair(old_lim, lim)] = bin_count;
            for (auto pos = errors.begin(); pos != p; ++pos)
            {
               os << pos->a << "," << pos->b << "," << pos->z << "," << pos->error << "\n";
            }
            errors.erase(errors.begin(), p);
         }
         ++current_bin;
      }

      tee_bins << "Results:\n\n";
      tee_bins << "#bin,Range,2^N,Count\n";
      int hash = 0;
      for (auto p = bins.begin(); p != bins.end(); ++p, ++hash)
      {
         tee_bins << hash << "," << p->first.first << "-" << p->first.second << "," << hash+1 << "," << p->second << std::endl;
      }
      if (evaluation_errors)
      {
         tee_bins << ",Failed,," << evaluation_errors << std::endl;
      }
      if (incalculable)
      {
         tee_bins << ",Incalculable,," << incalculable << std::endl;
      }
   }
   catch (const std::exception& e)
   {
      std::cout << "Terminating with unhandled exception: " << e.what() << std::endl;
   }

   return 0;
}