mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
1F1: begin to document error rates and testing methodology.
[CI SKIP]
This commit is contained in:
jzmaddock
217
tools/hypergeometric_1F1_error_plot.cpp
Normal file
217
tools/hypergeometric_1F1_error_plot.cpp
Normal file
@@ -0,0 +1,217 @@
|
||||
|
||||
#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/math/special_functions/hypergeometric_1F1.hpp>
|
||||
#include <boost/math/special_functions/hypergeometric_pFq.hpp>
|
||||
#include <boost/math/special_functions/relative_difference.hpp>
|
||||
|
||||
#include <boost/multiprecision/mpfr.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, int e)
|
||||
: a(a), b(b), z(z), error(e) {}
|
||||
double a, b, z;
|
||||
int 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, 3.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;
|
||||
tee_log << "Unexpected exception calculating value: " << std::endl;
|
||||
print_row(a, b, z, mp_expected, tee_log);
|
||||
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::itrunc(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())
|
||||
{
|
||||
std::ofstream os((basename + "_errors_" + std::to_string(current_bin + 1) + ".csv").c_str());
|
||||
os << "a,b,z,error\n";
|
||||
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);
|
||||
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;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user