libs/math/test/ooura_fourier_integral_test.cpp

   1 // Copyright Nick Thompson, 2019
   2 // Use, modification and distribution are subject to the
   3 // Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6 #define BOOST_TEST_MODULE test_ooura_fourier_transform
   7
   8 #include <cmath>
   9 #include <iostream>
  10 #include <boost/type_index.hpp>
  11 #include <boost/test/included/unit_test.hpp>
  12 #include <boost/test/tools/floating_point_comparison.hpp>
  13 #include <boost/math/quadrature/ooura_fourier_integrals.hpp>
  14 #include <boost/multiprecision/cpp_bin_float.hpp>
  15
  16 using boost::math::quadrature::ooura_fourier_sin;
  17 using boost::math::quadrature::ooura_fourier_cos;
  18 using boost::math::constants::pi;
  19
  20
  21 float float_tol = 10*std::numeric_limits<float>::epsilon();
  22 ooura_fourier_sin<float> float_sin_integrator(float_tol);
  23
  24 double double_tol = 10*std::numeric_limits<double>::epsilon();
  25 ooura_fourier_sin<double> double_sin_integrator(double_tol);
  26
  27 long double long_double_tol = 10*std::numeric_limits<long double>::epsilon();
  28 ooura_fourier_sin<long double> long_double_sin_integrator(long_double_tol);
  29
  30 template<class Real>
  31 auto get_sin_integrator() {
  32     if constexpr (std::is_same_v<Real, float>) {
  33         return float_sin_integrator;
  34     }
  35     if constexpr (std::is_same_v<Real, double>) {
  36         return double_sin_integrator;
  37     }
  38     if constexpr (std::is_same_v<Real, long double>) {
  39         return long_double_sin_integrator;
  40     }
  41 }
  42
  43 ooura_fourier_cos<float> float_cos_integrator(float_tol);
  44 ooura_fourier_cos<double> double_cos_integrator(double_tol);
  45 ooura_fourier_cos<long double> long_double_cos_integrator(long_double_tol);
  46
  47 template<class Real>
  48 auto get_cos_integrator() {
  49     if constexpr (std::is_same_v<Real, float>) {
  50         return float_cos_integrator;
  51     }
  52     if constexpr (std::is_same_v<Real, double>) {
  53         return double_cos_integrator;
  54     }
  55     if constexpr (std::is_same_v<Real, long double>) {
  56         return long_double_cos_integrator;
  57     }
  58 }
  59
  60
  61 template<class Real>
  62 void test_ooura_eta()
  63 {
  64     using boost::math::quadrature::detail::ooura_eta;
  65     std::cout << "Testing eta function on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  66     {
  67         Real x = 0;
  68         Real alpha = 7;
  69         auto [eta, eta_prime] = ooura_eta(x, alpha);
  70         BOOST_CHECK_SMALL(eta, (std::numeric_limits<Real>::min)());
  71         BOOST_CHECK_CLOSE_FRACTION(eta_prime, 2 + alpha + Real(1)/Real(4), 10*std::numeric_limits<Real>::epsilon());
  72     }
  73
  74     {
  75         Real alpha = 4;
  76         for (Real z = 0.125; z < 500; z += 0.125) {
  77             Real x = std::log(z);
  78             auto [eta, eta_prime] = ooura_eta(x, alpha);
  79             BOOST_CHECK_CLOSE_FRACTION(eta, 2*x + alpha*(1-1/z) + (z-1)/4, 10*std::numeric_limits<Real>::epsilon());
  80             BOOST_CHECK_CLOSE_FRACTION(eta_prime, 2 + alpha/z + z/4, 10*std::numeric_limits<Real>::epsilon());
  81         }
  82     }
  83 }
  84
  85 template<class Real>
  86 void test_ooura_sin_nodes_and_weights()
  87 {
  88     using boost::math::quadrature::detail::ooura_sin_node_and_weight;
  89     using boost::math::quadrature::detail::ooura_eta;
  90     std::cout << "Testing nodes and weights on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  91     {
  92         long n = 1;
  93         Real alpha = 1;
  94         Real h = 1;
  95         auto [node, weight] = ooura_sin_node_and_weight(n, h, alpha);
  96         Real expected_node = pi<Real>()/(1-exp(-ooura_eta(n*h, alpha).first));
  97         BOOST_CHECK_CLOSE_FRACTION(node,  expected_node,10*std::numeric_limits<Real>::epsilon());
  98     }
  99 }
 100
 101 template<class Real>
 102 void test_ooura_alpha() {
 103     std::cout << "Testing Ooura alpha on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 104     using std::sqrt;
 105     using std::log1p;
 106     using boost::math::quadrature::detail::calculate_ooura_alpha;
 107     Real alpha = calculate_ooura_alpha(Real(1));
 108     Real expected = 1/sqrt(16 + 4*log1p(pi<Real>()));
 109     BOOST_CHECK_CLOSE_FRACTION(alpha, expected, 10*std::numeric_limits<Real>::epsilon());
 110 }
 111
 112 void test_node_weight_precision_agreement()
 113 {
 114     using std::abs;
 115     using boost::math::quadrature::detail::ooura_sin_node_and_weight;
 116     using boost::math::quadrature::detail::ooura_eta;
 117     using boost::multiprecision::cpp_bin_float_quad;
 118     std::cout << "Testing agreement in two different precisions of nodes and weights\n";
 119     cpp_bin_float_quad alpha_quad = 1;
 120     long int_max = 128;
 121     cpp_bin_float_quad h_quad = 1/cpp_bin_float_quad(int_max);
 122     double alpha_dbl = 1;
 123     double h_dbl = static_cast<double>(h_quad);
 124     std::cout << std::fixed;
 125     for (long n = -1; n > -6*int_max; --n) {
 126         auto [node_dbl, weight_dbl] = ooura_sin_node_and_weight(n, h_dbl, alpha_dbl);
 127         auto p = ooura_sin_node_and_weight(n, h_quad, alpha_quad);
 128         double node_quad = static_cast<double>(p.first);
 129         double weight_quad = static_cast<double>(p.second);
 130         auto node_dist = abs(boost::math::float_distance(node_quad, node_dbl));
 131         if ( (weight_quad < 0 && weight_dbl > 0) || (weight_dbl < 0 && weight_quad > 0) ){
 132             std::cout << "Weights at different precisions have different signs!\n";
 133         } else {
 134             auto weight_dist = abs(boost::math::float_distance(weight_quad, weight_dbl));
 135             if (weight_dist > 100) {
 136                 std::cout << std::fixed;
 137                 std::cout <<"n =" << n << ", x = " << n*h_dbl << ", node distance = " << node_dist << ", weight distance = " << weight_dist << "\n";
 138                 std::cout << std::scientific;
 139                 std::cout << "computed weight = " << weight_dbl << ", actual weight = " << weight_quad << "\n";
 140             }
 141         }
 142     }
 143
 144 }
 145
 146 template<class Real>
 147 void test_sinc()
 148 {
 149     std::cout << "Testing sinc integral on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 150     using std::numeric_limits;
 151     Real tol = 50*numeric_limits<Real>::epsilon();
 152     auto integrator = get_sin_integrator<Real>();
 153     auto f = [](Real x)->Real { return 1/x; };
 154     Real omega = 1;
 155     while (omega < 10)
 156     {
 157         auto [Is, err] = integrator.integrate(f, omega);
 158         BOOST_CHECK_CLOSE_FRACTION(Is, pi<Real>()/2, tol);
 159
 160         auto [Isn, errn] = integrator.integrate(f, -omega);
 161         BOOST_CHECK_CLOSE_FRACTION(Isn, -pi<Real>()/2, tol);
 162         omega += 1;
 163     }
 164 }
 165
 166
 167 template<class Real>
 168 void test_exp()
 169 {
 170     std::cout << "Testing exponential integral on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 171     using std::exp;
 172     using std::numeric_limits;
 173     Real tol = 50*numeric_limits<Real>::epsilon();
 174     auto integrator = get_sin_integrator<Real>();
 175     auto f = [](Real x)->Real {return exp(-x);};
 176     Real omega = 1;
 177     while (omega < 5)
 178     {
 179         auto [Is, err] = integrator.integrate(f, omega);
 180         Real exact = omega/(1+omega*omega);
 181         BOOST_CHECK_CLOSE_FRACTION(Is, exact, tol);
 182         omega += 1;
 183     }
 184 }
 185
 186
 187 template<class Real>
 188 void test_root()
 189 {
 190     std::cout << "Testing integral of sin(kx)/sqrt(x) on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 191     using std::sqrt;
 192     using std::numeric_limits;
 193     Real tol = 10*numeric_limits<Real>::epsilon();
 194     auto integrator = get_sin_integrator<Real>();
 195     auto f = [](Real x)->Real { return 1/sqrt(x);};
 196     Real omega = 1;
 197     while (omega < 5) {
 198         auto [Is, err] = integrator.integrate(f, omega);
 199         Real exact = sqrt(pi<Real>()/(2*omega));
 200         BOOST_CHECK_CLOSE_FRACTION(Is, exact, 10*tol);
 201         omega += 1;
 202     }
 203 }
 204
 205 // See: https://scicomp.stackexchange.com/questions/32790/numerical-evaluation-of-highly-oscillatory-integral/32799#32799
 206 template<class Real>
 207 Real asymptotic(Real lambda) {
 208     using std::sin;
 209     using std::cos;
 210     using boost::math::constants::pi;
 211     Real I1 = cos(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/lambda);
 212     Real I2 = sin(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/(lambda*lambda*lambda))/8;
 213     return I1 + I2;
 214 }
 215
 216 template<class Real>
 217 void test_double_osc()
 218 {
 219     std::cout << "Testing double oscillation on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 220     using std::sqrt;
 221     using std::numeric_limits;
 222     auto integrator = get_sin_integrator<Real>();
 223     Real lambda = 7;
 224     auto f = [&lambda](Real x)->Real { return cos(lambda*cos(x))/x; };
 225     Real omega = 1;
 226     auto [Is, err] = integrator.integrate(f, omega);
 227     Real exact = asymptotic(lambda);
 228     BOOST_CHECK_CLOSE_FRACTION(2*Is, exact, 0.05);
 229 }
 230
 231 template<class Real>
 232 void test_zero_integrand()
 233 {
 234     // Make sure relative error tolerance doesn't break on zero integrand:
 235     std::cout << "Testing zero integrand on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 236     using std::sqrt;
 237     using std::numeric_limits;
 238     auto integrator = get_sin_integrator<Real>();
 239     auto f = [](Real /* x */)->Real { return Real(0); };
 240     Real omega = 1;
 241     auto [Is, err] = integrator.integrate(f, omega);
 242     Real exact = 0;
 243     BOOST_CHECK_EQUAL(Is, exact);
 244 }
 245
 246
 247 // This works, but doesn't recover the precision you want in a unit test:
 248 // template<class Real>
 249 // void test_log()
 250 // {
 251 //     std::cout << "Testing integral of log(x)sin(x) on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 252 //     using std::log;
 253 //     using std::exp;
 254 //     using std::numeric_limits;
 255 //     using boost::math::constants::euler;
 256 //     Real tol = 1000*numeric_limits<Real>::epsilon();
 257 //     auto f = [](Real x)->Real { return exp(-100*numeric_limits<Real>::epsilon()*x)*log(x);};
 258 //     Real omega = 1;
 259 //     Real Is = ooura_fourier_sin<decltype(f), Real>(f, omega, sqrt(numeric_limits<Real>::epsilon())/100);
 260 //     BOOST_CHECK_CLOSE_FRACTION(Is, -euler<Real>(), tol);
 261 // }
 262
 263
 264 template<class Real>
 265 void test_cos_integral1()
 266 {
 267     std::cout << "Testing integral of cos(x)/(x*x+1) on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 268     using std::exp;
 269     using boost::math::constants::half_pi;
 270     using boost::math::constants::e;
 271     using std::numeric_limits;
 272     Real tol = 10*numeric_limits<Real>::epsilon();
 273
 274     auto integrator = get_cos_integrator<Real>();
 275     auto f = [](Real x)->Real { return 1/(x*x+1);};
 276     Real omega = 1;
 277     auto [Is, err] = integrator.integrate(f, omega);
 278     Real exact = half_pi<Real>()/e<Real>();
 279     BOOST_CHECK_CLOSE_FRACTION(Is, exact, tol);
 280 }
 281
 282 template<class Real>
 283 void test_cos_integral2()
 284 {
 285     std::cout << "Testing integral of exp(-a*x) on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 286     using std::exp;
 287     using boost::math::constants::half_pi;
 288     using boost::math::constants::e;
 289     using std::numeric_limits;
 290     Real tol = 10*numeric_limits<Real>::epsilon();
 291
 292     auto integrator = get_cos_integrator<Real>();
 293     for (Real a = 1; a < 5; ++a) {
 294         auto f = [&a](Real x)->Real { return exp(-a*x);};
 295         for(Real omega = 1; omega < 5; ++omega) {
 296             auto [Is, err] = integrator.integrate(f, omega);
 297             Real exact = a/(a*a+omega*omega);
 298             BOOST_CHECK_CLOSE_FRACTION(Is, exact, 10*tol);
 299         }
 300     }
 301 }
 302
 303 template<class Real>
 304 void test_nodes()
 305 {
 306     std::cout << "Testing nodes and weights on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 307     auto sin_integrator = get_sin_integrator<Real>();
 308
 309     auto const & big_nodes = sin_integrator.big_nodes();
 310     for (auto & node_row : big_nodes) {
 311         Real t0 = node_row[0];
 312         for (size_t i = 1; i < node_row.size(); ++i) {
 313             Real t1 = node_row[i];
 314             BOOST_CHECK(t1 > t0);
 315             t0 = t1;
 316         }
 317     }
 318
 319     auto const & little_nodes = sin_integrator.little_nodes();
 320     for (auto & node_row : little_nodes) {
 321         Real t0 = node_row[0];
 322         for (size_t i = 1; i < node_row.size(); ++i) {
 323             Real t1 = node_row[i];
 324             BOOST_CHECK(t1 < t0);
 325             t0 = t1;
 326         }
 327     }
 328 }
 329
 330
 331 BOOST_AUTO_TEST_CASE(ooura_fourier_transform_test)
 332 {
 333     test_cos_integral1<float>();
 334     test_cos_integral1<double>();
 335     test_cos_integral1<long double>();
 336
 337     test_cos_integral2<float>();
 338     test_cos_integral2<double>();
 339     test_cos_integral2<long double>();
 340
 341     //test_node_weight_precision_agreement();
 342     test_zero_integrand<float>();
 343     test_zero_integrand<double>();
 344
 345     test_ooura_eta<float>();
 346     test_ooura_eta<double>();
 347     test_ooura_eta<long double>();
 348
 349     test_ooura_sin_nodes_and_weights<float>();
 350     test_ooura_sin_nodes_and_weights<double>();
 351     test_ooura_sin_nodes_and_weights<long double>();
 352
 353     test_ooura_alpha<float>();
 354     test_ooura_alpha<double>();
 355     test_ooura_alpha<long double>();
 356
 357     test_sinc<float>();
 358     test_sinc<double>();
 359     test_sinc<long double>();
 360
 361     test_exp<float>();
 362     test_exp<double>();
 363     test_exp<long double>();
 364
 365     test_root<float>();
 366     test_root<double>();
 367
 368     test_double_osc<float>();
 369     test_double_osc<double>();
 370     // Takes too long!
 371     //test_double_osc<long double>();
 372
 373     // This test should be last:
 374     test_nodes<float>();
 375     test_nodes<double>();
 376     test_nodes<long double>();
 377 }