123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377 |
- // Copyright Nick Thompson, 2019
- // 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 test_ooura_fourier_transform
- #include <cmath>
- #include <iostream>
- #include <boost/type_index.hpp>
- #include <boost/test/included/unit_test.hpp>
- #include <boost/test/tools/floating_point_comparison.hpp>
- #include <boost/math/quadrature/ooura_fourier_integrals.hpp>
- #include <boost/multiprecision/cpp_bin_float.hpp>
- using boost::math::quadrature::ooura_fourier_sin;
- using boost::math::quadrature::ooura_fourier_cos;
- using boost::math::constants::pi;
- float float_tol = 10*std::numeric_limits<float>::epsilon();
- ooura_fourier_sin<float> float_sin_integrator(float_tol);
- double double_tol = 10*std::numeric_limits<double>::epsilon();
- ooura_fourier_sin<double> double_sin_integrator(double_tol);
- long double long_double_tol = 10*std::numeric_limits<long double>::epsilon();
- ooura_fourier_sin<long double> long_double_sin_integrator(long_double_tol);
- template<class Real>
- auto get_sin_integrator() {
- if constexpr (std::is_same_v<Real, float>) {
- return float_sin_integrator;
- }
- if constexpr (std::is_same_v<Real, double>) {
- return double_sin_integrator;
- }
- if constexpr (std::is_same_v<Real, long double>) {
- return long_double_sin_integrator;
- }
- }
- ooura_fourier_cos<float> float_cos_integrator(float_tol);
- ooura_fourier_cos<double> double_cos_integrator(double_tol);
- ooura_fourier_cos<long double> long_double_cos_integrator(long_double_tol);
- template<class Real>
- auto get_cos_integrator() {
- if constexpr (std::is_same_v<Real, float>) {
- return float_cos_integrator;
- }
- if constexpr (std::is_same_v<Real, double>) {
- return double_cos_integrator;
- }
- if constexpr (std::is_same_v<Real, long double>) {
- return long_double_cos_integrator;
- }
- }
- template<class Real>
- void test_ooura_eta()
- {
- using boost::math::quadrature::detail::ooura_eta;
- std::cout << "Testing eta function on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- {
- Real x = 0;
- Real alpha = 7;
- auto [eta, eta_prime] = ooura_eta(x, alpha);
- BOOST_CHECK_SMALL(eta, (std::numeric_limits<Real>::min)());
- BOOST_CHECK_CLOSE_FRACTION(eta_prime, 2 + alpha + Real(1)/Real(4), 10*std::numeric_limits<Real>::epsilon());
- }
- {
- Real alpha = 4;
- for (Real z = 0.125; z < 500; z += 0.125) {
- Real x = std::log(z);
- auto [eta, eta_prime] = ooura_eta(x, alpha);
- BOOST_CHECK_CLOSE_FRACTION(eta, 2*x + alpha*(1-1/z) + (z-1)/4, 10*std::numeric_limits<Real>::epsilon());
- BOOST_CHECK_CLOSE_FRACTION(eta_prime, 2 + alpha/z + z/4, 10*std::numeric_limits<Real>::epsilon());
- }
- }
- }
- template<class Real>
- void test_ooura_sin_nodes_and_weights()
- {
- using boost::math::quadrature::detail::ooura_sin_node_and_weight;
- using boost::math::quadrature::detail::ooura_eta;
- std::cout << "Testing nodes and weights on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- {
- long n = 1;
- Real alpha = 1;
- Real h = 1;
- auto [node, weight] = ooura_sin_node_and_weight(n, h, alpha);
- Real expected_node = pi<Real>()/(1-exp(-ooura_eta(n*h, alpha).first));
- BOOST_CHECK_CLOSE_FRACTION(node, expected_node,10*std::numeric_limits<Real>::epsilon());
- }
- }
- template<class Real>
- void test_ooura_alpha() {
- std::cout << "Testing Ooura alpha on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- using std::sqrt;
- using std::log1p;
- using boost::math::quadrature::detail::calculate_ooura_alpha;
- Real alpha = calculate_ooura_alpha(Real(1));
- Real expected = 1/sqrt(16 + 4*log1p(pi<Real>()));
- BOOST_CHECK_CLOSE_FRACTION(alpha, expected, 10*std::numeric_limits<Real>::epsilon());
- }
- void test_node_weight_precision_agreement()
- {
- using std::abs;
- using boost::math::quadrature::detail::ooura_sin_node_and_weight;
- using boost::math::quadrature::detail::ooura_eta;
- using boost::multiprecision::cpp_bin_float_quad;
- std::cout << "Testing agreement in two different precisions of nodes and weights\n";
- cpp_bin_float_quad alpha_quad = 1;
- long int_max = 128;
- cpp_bin_float_quad h_quad = 1/cpp_bin_float_quad(int_max);
- double alpha_dbl = 1;
- double h_dbl = static_cast<double>(h_quad);
- std::cout << std::fixed;
- for (long n = -1; n > -6*int_max; --n) {
- auto [node_dbl, weight_dbl] = ooura_sin_node_and_weight(n, h_dbl, alpha_dbl);
- auto p = ooura_sin_node_and_weight(n, h_quad, alpha_quad);
- double node_quad = static_cast<double>(p.first);
- double weight_quad = static_cast<double>(p.second);
- auto node_dist = abs(boost::math::float_distance(node_quad, node_dbl));
- if ( (weight_quad < 0 && weight_dbl > 0) || (weight_dbl < 0 && weight_quad > 0) ){
- std::cout << "Weights at different precisions have different signs!\n";
- } else {
- auto weight_dist = abs(boost::math::float_distance(weight_quad, weight_dbl));
- if (weight_dist > 100) {
- std::cout << std::fixed;
- std::cout <<"n =" << n << ", x = " << n*h_dbl << ", node distance = " << node_dist << ", weight distance = " << weight_dist << "\n";
- std::cout << std::scientific;
- std::cout << "computed weight = " << weight_dbl << ", actual weight = " << weight_quad << "\n";
- }
- }
- }
- }
- template<class Real>
- void test_sinc()
- {
- std::cout << "Testing sinc integral on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- using std::numeric_limits;
- Real tol = 50*numeric_limits<Real>::epsilon();
- auto integrator = get_sin_integrator<Real>();
- auto f = [](Real x)->Real { return 1/x; };
- Real omega = 1;
- while (omega < 10)
- {
- auto [Is, err] = integrator.integrate(f, omega);
- BOOST_CHECK_CLOSE_FRACTION(Is, pi<Real>()/2, tol);
- auto [Isn, errn] = integrator.integrate(f, -omega);
- BOOST_CHECK_CLOSE_FRACTION(Isn, -pi<Real>()/2, tol);
- omega += 1;
- }
- }
- template<class Real>
- void test_exp()
- {
- std::cout << "Testing exponential integral on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- using std::exp;
- using std::numeric_limits;
- Real tol = 50*numeric_limits<Real>::epsilon();
- auto integrator = get_sin_integrator<Real>();
- auto f = [](Real x)->Real {return exp(-x);};
- Real omega = 1;
- while (omega < 5)
- {
- auto [Is, err] = integrator.integrate(f, omega);
- Real exact = omega/(1+omega*omega);
- BOOST_CHECK_CLOSE_FRACTION(Is, exact, tol);
- omega += 1;
- }
- }
- template<class Real>
- void test_root()
- {
- std::cout << "Testing integral of sin(kx)/sqrt(x) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- using std::sqrt;
- using std::numeric_limits;
- Real tol = 10*numeric_limits<Real>::epsilon();
- auto integrator = get_sin_integrator<Real>();
- auto f = [](Real x)->Real { return 1/sqrt(x);};
- Real omega = 1;
- while (omega < 5) {
- auto [Is, err] = integrator.integrate(f, omega);
- Real exact = sqrt(pi<Real>()/(2*omega));
- BOOST_CHECK_CLOSE_FRACTION(Is, exact, 10*tol);
- omega += 1;
- }
- }
- // See: https://scicomp.stackexchange.com/questions/32790/numerical-evaluation-of-highly-oscillatory-integral/32799#32799
- template<class Real>
- Real asymptotic(Real lambda) {
- using std::sin;
- using std::cos;
- using boost::math::constants::pi;
- Real I1 = cos(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/lambda);
- Real I2 = sin(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/(lambda*lambda*lambda))/8;
- return I1 + I2;
- }
- template<class Real>
- void test_double_osc()
- {
- std::cout << "Testing double oscillation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- using std::sqrt;
- using std::numeric_limits;
- auto integrator = get_sin_integrator<Real>();
- Real lambda = 7;
- auto f = [&lambda](Real x)->Real { return cos(lambda*cos(x))/x; };
- Real omega = 1;
- auto [Is, err] = integrator.integrate(f, omega);
- Real exact = asymptotic(lambda);
- BOOST_CHECK_CLOSE_FRACTION(2*Is, exact, 0.05);
- }
- template<class Real>
- void test_zero_integrand()
- {
- // Make sure relative error tolerance doesn't break on zero integrand:
- std::cout << "Testing zero integrand on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- using std::sqrt;
- using std::numeric_limits;
- auto integrator = get_sin_integrator<Real>();
- auto f = [](Real /* x */)->Real { return Real(0); };
- Real omega = 1;
- auto [Is, err] = integrator.integrate(f, omega);
- Real exact = 0;
- BOOST_CHECK_EQUAL(Is, exact);
- }
- // This works, but doesn't recover the precision you want in a unit test:
- // template<class Real>
- // void test_log()
- // {
- // std::cout << "Testing integral of log(x)sin(x) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- // using std::log;
- // using std::exp;
- // using std::numeric_limits;
- // using boost::math::constants::euler;
- // Real tol = 1000*numeric_limits<Real>::epsilon();
- // auto f = [](Real x)->Real { return exp(-100*numeric_limits<Real>::epsilon()*x)*log(x);};
- // Real omega = 1;
- // Real Is = ooura_fourier_sin<decltype(f), Real>(f, omega, sqrt(numeric_limits<Real>::epsilon())/100);
- // BOOST_CHECK_CLOSE_FRACTION(Is, -euler<Real>(), tol);
- // }
- template<class Real>
- void test_cos_integral1()
- {
- std::cout << "Testing integral of cos(x)/(x*x+1) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- using std::exp;
- using boost::math::constants::half_pi;
- using boost::math::constants::e;
- using std::numeric_limits;
- Real tol = 10*numeric_limits<Real>::epsilon();
- auto integrator = get_cos_integrator<Real>();
- auto f = [](Real x)->Real { return 1/(x*x+1);};
- Real omega = 1;
- auto [Is, err] = integrator.integrate(f, omega);
- Real exact = half_pi<Real>()/e<Real>();
- BOOST_CHECK_CLOSE_FRACTION(Is, exact, tol);
- }
- template<class Real>
- void test_cos_integral2()
- {
- std::cout << "Testing integral of exp(-a*x) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- using std::exp;
- using boost::math::constants::half_pi;
- using boost::math::constants::e;
- using std::numeric_limits;
- Real tol = 10*numeric_limits<Real>::epsilon();
- auto integrator = get_cos_integrator<Real>();
- for (Real a = 1; a < 5; ++a) {
- auto f = [&a](Real x)->Real { return exp(-a*x);};
- for(Real omega = 1; omega < 5; ++omega) {
- auto [Is, err] = integrator.integrate(f, omega);
- Real exact = a/(a*a+omega*omega);
- BOOST_CHECK_CLOSE_FRACTION(Is, exact, 10*tol);
- }
- }
- }
- template<class Real>
- void test_nodes()
- {
- std::cout << "Testing nodes and weights on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- auto sin_integrator = get_sin_integrator<Real>();
- auto const & big_nodes = sin_integrator.big_nodes();
- for (auto & node_row : big_nodes) {
- Real t0 = node_row[0];
- for (size_t i = 1; i < node_row.size(); ++i) {
- Real t1 = node_row[i];
- BOOST_CHECK(t1 > t0);
- t0 = t1;
- }
- }
- auto const & little_nodes = sin_integrator.little_nodes();
- for (auto & node_row : little_nodes) {
- Real t0 = node_row[0];
- for (size_t i = 1; i < node_row.size(); ++i) {
- Real t1 = node_row[i];
- BOOST_CHECK(t1 < t0);
- t0 = t1;
- }
- }
- }
- BOOST_AUTO_TEST_CASE(ooura_fourier_transform_test)
- {
- test_cos_integral1<float>();
- test_cos_integral1<double>();
- test_cos_integral1<long double>();
- test_cos_integral2<float>();
- test_cos_integral2<double>();
- test_cos_integral2<long double>();
- //test_node_weight_precision_agreement();
- test_zero_integrand<float>();
- test_zero_integrand<double>();
- test_ooura_eta<float>();
- test_ooura_eta<double>();
- test_ooura_eta<long double>();
- test_ooura_sin_nodes_and_weights<float>();
- test_ooura_sin_nodes_and_weights<double>();
- test_ooura_sin_nodes_and_weights<long double>();
- test_ooura_alpha<float>();
- test_ooura_alpha<double>();
- test_ooura_alpha<long double>();
- test_sinc<float>();
- test_sinc<double>();
- test_sinc<long double>();
- test_exp<float>();
- test_exp<double>();
- test_exp<long double>();
- test_root<float>();
- test_root<double>();
- test_double_osc<float>();
- test_double_osc<double>();
- // Takes too long!
- //test_double_osc<long double>();
- // This test should be last:
- test_nodes<float>();
- test_nodes<double>();
- test_nodes<long double>();
- }
|