math/test/exp_sinh_quadrature_test.cpp

// Copyright Nick Thompson, 2017
// 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_TEST_MODULE exp_sinh_quadrature_test

#include <boost/math/concepts/real_concept.hpp>
#include <boost/test/included/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/quadrature/exp_sinh.hpp>
#include <boost/math/special_functions/sinc.hpp>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/math/special_functions/next.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/sinc.hpp>
#include <boost/type_traits/is_class.hpp>

using std::exp;
using std::cos;
using std::tan;
using std::log;
using std::sqrt;
using std::abs;
using std::sinh;
using std::cosh;
using std::pow;
using std::atan;
using boost::multiprecision::cpp_bin_float_50;
using boost::multiprecision::cpp_bin_float_100;
using boost::multiprecision::cpp_bin_float_quad;
using boost::math::constants::pi;
using boost::math::constants::half_pi;
using boost::math::constants::two_div_pi;
using boost::math::constants::half;
using boost::math::constants::third;
using boost::math::constants::half;
using boost::math::constants::third;
using boost::math::constants::catalan;
using boost::math::constants::ln_two;
using boost::math::constants::root_two;
using boost::math::constants::root_two_pi;
using boost::math::constants::root_pi;
using boost::math::quadrature::exp_sinh;

#if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) && !defined(TEST5) && !defined(TEST6) && !defined(TEST7)
#  define TEST1
#  define TEST2
#  define TEST3
#  define TEST4
#  define TEST5
#  define TEST6
#  define TEST7
#endif

#ifdef BOOST_MSVC
#pragma warning (disable:4127)
#endif

//
// Coefficient generation code:
//
template <class T>
void print_levels(const T& v, const char* suffix)
{
   std::cout << "{\n";
   for (unsigned i = 0; i < v.size(); ++i)
   {
      std::cout << "      { ";
      for (unsigned j = 0; j < v[i].size(); ++j)
      {
         std::cout << v[i][j] << suffix << ", ";
      }
      std::cout << "},\n";
   }
   std::cout << "   };\n";
}

template <class T>
void print_levels(const std::pair<T, T>& p, const char* suffix = "")
{
   std::cout << "   static const std::vector<std::vector<Real> > abscissa = ";
   print_levels(p.first, suffix);
   std::cout << "   static const std::vector<std::vector<Real> > weights = ";
   print_levels(p.second, suffix);
}

template <class Real, class TargetType>
std::pair<std::vector<std::vector<Real>>, std::vector<std::vector<Real>> > generate_constants(unsigned max_rows)
{
   using boost::math::constants::half_pi;
   using boost::math::constants::two_div_pi;
   using boost::math::constants::pi;
   auto g = [](Real t) { return exp(half_pi<Real>()*sinh(t)); };
   auto w = [](Real t) { return cosh(t)*half_pi<Real>()*exp(half_pi<Real>()*sinh(t)); };

   std::vector<std::vector<Real>> abscissa, weights;

   std::vector<Real> temp;

   Real tmp = (Real(boost::math::tools::log_min_value<TargetType>()) + log(Real(boost::math::tools::epsilon<TargetType>())))*0.5f;
   Real t_min = asinh(two_div_pi<Real>()*tmp);
   // truncate t_min to an exact binary value:
   t_min = floor(t_min * 128) / 128;

   std::cout << "m_t_min = " << t_min << ";\n";

   // t_max is chosen to make g'(t_max) ~ sqrt(max) (g' grows faster than g).
   // This will allow some flexibility on the users part; they can at least square a number function without overflow.
   // But there is no unique choice; the further out we can evaluate the function, the better we can do on slowly decaying integrands.
   const Real t_max = log(2 * two_div_pi<Real>()*log(2 * two_div_pi<Real>()*sqrt(Real(boost::math::tools::max_value<TargetType>()))));

   Real h = 1;
   for (Real t = t_min; t < t_max; t += h)
   {
      temp.push_back(g(t));
   }
   abscissa.push_back(temp);
   temp.clear();

   for (Real t = t_min; t < t_max; t += h)
   {
      temp.push_back(w(t * h));
   }
   weights.push_back(temp);
   temp.clear();

   for (unsigned row = 1; row < max_rows; ++row)
   {
      h /= 2;
      for (Real t = t_min + h; t < t_max; t += 2 * h)
         temp.push_back(g(t));
      abscissa.push_back(temp);
      temp.clear();
   }
   h = 1;
   for (unsigned row = 1; row < max_rows; ++row)
   {
      h /= 2;
      for (Real t = t_min + h; t < t_max; t += 2 * h)
         temp.push_back(w(t));
      weights.push_back(temp);
      temp.clear();
   }

   return std::make_pair(abscissa, weights);
}


template <class Real>
const exp_sinh<Real>& get_integrator()
{
   static const exp_sinh<Real> integrator(14);
   return integrator;
}

template <class Real>
Real get_convergence_tolerance()
{
   return boost::math::tools::root_epsilon<Real>();
}

template<class Real>
void test_right_limit_infinite()
{
    std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = 10 * boost::math::tools::epsilon<Real>();
    Real Q;
    Real Q_expected;
    Real error;
    Real L1;
    auto integrator = get_integrator<Real>();

    // Example 12
    const auto f2 = [](const Real& t) { return exp(-t)/sqrt(t); };
    Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = root_pi<Real>();
    Real tol_mult = 1;
    // Multiprecision type have higher error rates, probably evaluation of f() is less accurate:
    if (std::numeric_limits<Real>::digits10 > std::numeric_limits<long double>::digits10)
       tol_mult = 12;
    else if (std::numeric_limits<Real>::digits10 > std::numeric_limits<double>::digits10)
       tol_mult = 5;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * tol_mult);
    // The integrand is strictly positive, so it coincides with the value of the integral:
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * tol_mult);

    auto f3 = [](Real t)->Real { Real z = exp(-t); if (z == 0) { return z; } return z*cos(t); };
    Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = half<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    Q = integrator.integrate(f3, 10, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = boost::lexical_cast<Real>("-6.6976341310426674140007086979326069121526743314567805278252392932e-6");
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol);
    // Integrating through zero risks precision loss:
    Q = integrator.integrate(f3, -10, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = boost::lexical_cast<Real>("-15232.3213626280525704332288302799653087046646639974940243044623285817777006");
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, std::numeric_limits<Real>::digits10 > 30 ? 1000 * tol : tol);

    auto f4 = [](Real t)->Real { return 1/(1+t*t); };
    Q = integrator.integrate(f4, 1, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = pi<Real>()/4;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
    Q = integrator.integrate(f4, 20, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = boost::lexical_cast<Real>("0.0499583957219427614100062870348448814912770804235071744108534548299835954767");
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
    Q = integrator.integrate(f4, 500, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = boost::lexical_cast<Real>("0.0019999973333397333150476759363217553199063513829126652556286269630");
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
}

template<class Real>
void test_left_limit_infinite()
{
    std::cout << "Testing left limit infinite for 1/(1+t^2) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = 10 * boost::math::tools::epsilon<Real>();
    Real Q;
    Real Q_expected;
    Real error;
    Real L1;
    auto integrator = get_integrator<Real>();

    // Example 11:
    auto f1 = [](const Real& t) { return 1/(1+t*t);};
    Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), 0, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = half_pi<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
    Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), -20, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = boost::lexical_cast<Real>("0.0499583957219427614100062870348448814912770804235071744108534548299835954767");
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
    Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), -500, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = boost::lexical_cast<Real>("0.0019999973333397333150476759363217553199063513829126652556286269630");
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
}


// Some examples of tough integrals from NR, section 4.5.4:
template<class Real>
void test_nr_examples()
{
    using std::sin;
    using std::cos;
    using std::pow;
    using std::exp;
    using std::sqrt;
    std::cout << "Testing type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = 10 * boost::math::tools::epsilon<Real>();
    std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
    Real Q;
    Real Q_expected;
    Real L1;
    Real error;
    auto integrator = get_integrator<Real>();

    auto f0 = [] (Real)->Real { return (Real) 0; };
    Q = integrator.integrate(f0, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = 0;
    BOOST_CHECK_CLOSE_FRACTION(Q, 0.0f, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, 0.0f, tol);

    auto f = [](const Real& x) { return 1/(1+x*x); };
    Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = half_pi<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);

    auto f1 = [](Real x)->Real {
        Real z1 = exp(-x);
        if (z1 == 0)
        {
            return (Real) 0;
        }
        Real z2 = pow(x, -3*half<Real>())*z1;
        if (z2 == 0)
        {
            return (Real) 0;
        }
        return sin(x*half<Real>())*z2;
    };

    Q = integrator.integrate(f1, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = sqrt(pi<Real>()*(sqrt((Real) 5) - 2));

    // The integrand is oscillatory; the accuracy is low.
    Real tol_mul = 1;
    if (std::numeric_limits<Real>::digits10 > 40)
       tol_mul = 500000;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol);

    auto f2 = [](Real x)->Real { return pow(x, -(Real) 2/(Real) 7)*exp(-x*x); };
    Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = half<Real>()*boost::math::tgamma((Real) 5/ (Real) 14);
    tol_mul = 1;
    if (std::numeric_limits<Real>::is_specialized == false)
       tol_mul = 3;
    if (std::numeric_limits<Real>::digits10 > 40)
       tol_mul = 100;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol_mul * tol);

    auto f3 = [](Real x)->Real { return (Real) 1/ (sqrt(x)*(1+x)); };
    Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = pi<Real>();

    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10*boost::math::tools::epsilon<Real>());
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, 10*boost::math::tools::epsilon<Real>());

    auto f4 = [](const Real& t) { return exp(-t*t*half<Real>()); };
    Q = integrator.integrate(f4, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = root_two_pi<Real>()/2;
    tol_mul = 1;
    if (std::numeric_limits<Real>::digits10 > 40)
       tol_mul = 5000;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol_mul * tol);

    auto f5 = [](const Real& t) { return 1/cosh(t);};
    Q = integrator.integrate(f5, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = half_pi<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * 12);   // Fails at float precision without higher error rate
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * 12);
}

// Definite integrals found in the CRC Handbook of Mathematical Formulas
template<class Real>
void test_crc()
{
    using std::sin;
    using std::pow;
    using std::exp;
    using std::sqrt;
    using std::log;
    using std::cos;
    std::cout << "Testing integral from CRC handbook on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = 10 * boost::math::tools::epsilon<Real>();
    std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
    Real Q;
    Real Q_expected;
    Real L1;
    Real error;
    auto integrator = get_integrator<Real>();

    auto f0 = [](const Real& x) { return log(x)*exp(-x); };
    Q = integrator.integrate(f0, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = -boost::math::constants::euler<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);

    // Test the integral representation of the gamma function:
    auto f1 = [](Real t)->Real { Real x = exp(-t);
        if(x == 0)
        {
            return (Real) 0;
        }
        return pow(t, (Real) 12 - 1)*x;
    };

    Q = integrator.integrate(f1, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = boost::math::tgamma(12.0f);
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);

    // Integral representation of the modified bessel function:
    // K_5(12)
    auto f2 = [](Real t)->Real {
        Real x = exp(-12*cosh(t));
        if (x == 0)
        {
            return (Real) 0;
        }
        return x*cosh(5*t);
    };
    Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = boost::math::cyl_bessel_k<int, Real>(5, 12);
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    // Laplace transform of cos(at)
    Real a = 20;
    Real s = 1;
    auto f3 = [&](Real t)->Real {
        Real x = exp(-s*t);
        if (x == 0)
        {
            return (Real) 0;
        }
        return cos(a*t)*x;
    };

    // For high oscillation frequency, the quadrature sum is ill-conditioned.
    Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = s/(a*a+s*s);
    // Since the integrand is oscillatory, we increase the tolerance:
    Real tol_mult = 10;
    // Multiprecision type have higher error rates, probably evaluation of f() is less accurate:
    if (!boost::is_class<Real>::value)
    {
       if (std::numeric_limits<Real>::digits10 > std::numeric_limits<double>::digits10)
          tol_mult = 5000; // we should really investigate this more??
       BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mult*tol);
    }

    //
    // This one doesn't pass for real_concept..
    //
    if (std::numeric_limits<Real>::is_specialized)
    {
       // Laplace transform of J_0(t):
       auto f4 = [&](Real t)->Real {
          Real x = exp(-s*t);
          if (x == 0)
          {
             return (Real)0;
          }
          return boost::math::cyl_bessel_j(0, t)*x;
       };

       Q = integrator.integrate(f4, get_convergence_tolerance<Real>(), &error, &L1);
       Q_expected = 1 / sqrt(1 + s*s);
       tol_mult = 3;
       // Multiprecision type have higher error rates, probably evaluation of f() is less accurate:
       if (std::numeric_limits<Real>::digits10 > std::numeric_limits<long double>::digits10)
          tol_mult = 750;
       BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mult * tol);
    }
    auto f6 = [](const Real& t) { return exp(-t*t)*log(t);};
    Q = integrator.integrate(f6, get_convergence_tolerance<Real>(), &error, &L1);
    Q_expected = -boost::math::constants::root_pi<Real>()*(boost::math::constants::euler<Real>() + 2*ln_two<Real>())/4;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);

    // CRC Section 5.5, integral 591
    // The parameter p allows us to control the strength of the singularity.
    // Rapid convergence is not guaranteed for this function, as the branch cut makes it non-analytic on a disk.
    // This converges only when our test type has an extended exponent range as all the area of the integral
    // occurs so close to 0 (or 1) that we need abscissa values exceptionally small to find it.
    // "There's a lot of room at the bottom".
    // This version is transformed via argument substitution (exp(-x) for x) so that the integral is spread
    // over (0, INF).
    tol *= boost::math::tools::digits<Real>() > 100 ? 100000 : 75;
    for (Real pn = 99; pn > 0; pn -= 10) {
       Real p = pn / 100;
       auto f = [&](Real x)->Real
       {
          return exp(-x * (1 - p) + p * log(-boost::math::expm1(-x)));
       };
       Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
       Q_expected = 1 / boost::math::sinc_pi(p*pi<Real>());
       /*
       std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << p << std::endl;
       std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << Q << std::endl;
       std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << Q_expected << std::endl;
       std::cout << fabs((Q - Q_expected) / Q_expected) << std::endl;
       */
       BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    }
    // and for p < 1:
    for (Real p = -0.99; p < 0; p += 0.1) {
       auto f = [&](Real x)->Real
       {
          return exp(-p * log(-boost::math::expm1(-x)) - (1 + p) * x);
       };
       Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
       Q_expected = 1 / boost::math::sinc_pi(p*pi<Real>());
       BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    }
}


BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test)
{
   //
   // Uncomment to generate the coefficients:
   //

   /*
   std::cout << std::scientific << std::setprecision(8);
   print_levels(generate_constants<cpp_bin_float_100, float>(8), "f");
   std::cout << std::setprecision(18);
   print_levels(generate_constants<cpp_bin_float_100, double>(8), "");
   std::cout << std::setprecision(35);
   print_levels(generate_constants<cpp_bin_float_100, cpp_bin_float_quad>(8), "L");
   */


#ifdef TEST1
    test_left_limit_infinite<float>();
    test_right_limit_infinite<float>();
    test_nr_examples<float>();
    test_crc<float>();
#endif
#ifdef TEST2
    test_left_limit_infinite<double>();
    test_right_limit_infinite<double>();
    test_nr_examples<double>();
    test_crc<double>();
#endif
#ifdef TEST3
    test_left_limit_infinite<long double>();
    test_right_limit_infinite<long double>();
    test_nr_examples<long double>();
    test_crc<long double>();
#endif
#ifdef TEST4
    test_left_limit_infinite<cpp_bin_float_quad>();
    test_right_limit_infinite<cpp_bin_float_quad>();
    test_nr_examples<cpp_bin_float_quad>();
    test_crc<cpp_bin_float_quad>();
#endif

#ifdef TEST5
    test_left_limit_infinite<boost::math::concepts::real_concept>();
    test_right_limit_infinite<boost::math::concepts::real_concept>();
    test_nr_examples<boost::math::concepts::real_concept>();
    test_crc<boost::math::concepts::real_concept>();
#endif

#ifdef TEST6
    test_left_limit_infinite<boost::multiprecision::cpp_bin_float_50>();
    test_right_limit_infinite<boost::multiprecision::cpp_bin_float_50>();
    test_nr_examples<boost::multiprecision::cpp_bin_float_50>();
    test_crc<boost::multiprecision::cpp_bin_float_50>();
#endif
#ifdef TEST7
    test_left_limit_infinite<boost::multiprecision::cpp_dec_float_50>();
    test_right_limit_infinite<boost::multiprecision::cpp_dec_float_50>();
    test_nr_examples<boost::multiprecision::cpp_dec_float_50>();
    test_crc<boost::multiprecision::cpp_dec_float_50>();
#endif
}