1 /*
2 * Copyright Nick Thompson, 2020
3 * Use, modification and distribution are subject to the
4 * Boost Software License, Version 1.0. (See accompanying file
5 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 */
7
8 #include "math_unit_test.hpp"
9 #include <numeric>
10 #include <utility>
11 #include <iomanip>
12 #include <iostream>
13 #include <random>
14 #include <cmath>
15 #include <boost/core/demangle.hpp>
16 #include <boost/hana/for_each.hpp>
17 #include <boost/hana/ext/std/integer_sequence.hpp>
18 #include <boost/math/tools/condition_numbers.hpp>
19 #include <boost/math/special_functions/daubechies_wavelet.hpp>
20 #include <boost/math/special_functions/next.hpp>
21 #include <boost/math/quadrature/trapezoidal.hpp>
22
23 #ifdef BOOST_HAS_FLOAT128
24 #include <boost/multiprecision/float128.hpp>
25 using boost::multiprecision::float128;
26 #endif
27
28
29 using boost::math::constants::pi;
30 using boost::math::constants::root_two;
31
32 template<typename Real>
test_exact_value()33 void test_exact_value()
34 {
35 // The global phase of the wavelet is not constrained by anything other than convention.
36 // Make sure that our conventions match the rest of the world:
37 auto psi = boost::math::daubechies_wavelet<Real, 2>(2);
38 Real computed = psi(1);
39 Real expected = -1.366025403784439;
40 CHECK_MOLLIFIED_CLOSE(expected, computed, 0.0001);
41 }
42
43 template<typename Real, int p>
test_quadratures()44 void test_quadratures()
45 {
46 std::cout << "Testing quadratures of " << p << " vanishing moment Daubechies wavelet on type " << boost::core::demangle(typeid(Real).name()) << "\n";
47 using boost::math::quadrature::trapezoidal;
48 auto psi = boost::math::daubechies_wavelet<Real, p>();
49 std::cout << "Wavelet functor size is " << psi.bytes() << " bytes" << std::endl;
50
51 Real tol = std::numeric_limits<Real>::epsilon();
52 Real error_estimate = std::numeric_limits<Real>::quiet_NaN();
53 Real L1 = std::numeric_limits<Real>::quiet_NaN();
54 auto [a, b] = psi.support();
55 CHECK_ULP_CLOSE(Real(-p+1), a, 0);
56 CHECK_ULP_CLOSE(Real(p), b, 0);
57 // A wavelet is a function of zero average; ensure the quadrature over its support is zero.
58 Real Q = trapezoidal(psi, a, b, tol, 15, &error_estimate, &L1);
59 if (!CHECK_MOLLIFIED_CLOSE(Real(0), Q, Real(0.0001)))
60 {
61 std::cerr << " Quadrature of " << p << " vanishing moment wavelet does not vanish.\n";
62 std::cerr << " Error estimate: " << error_estimate << ", L1 norm: " << L1 << "\n";
63 }
64 auto psi_sq = [psi](Real x) {
65 Real t = psi(x);
66 return t*t;
67 };
68 Q = trapezoidal(psi_sq, a, b, tol, 15, &error_estimate, &L1);
69 Real quad_tol = 2000*std::sqrt(std::numeric_limits<Real>::epsilon())/(p*p*p);
70 if (!CHECK_MOLLIFIED_CLOSE(Real(1), Q, quad_tol))
71 {
72 std::cerr << " L2 norm of " << p << " vanishing moment wavelet does not vanish.\n";
73 std::cerr << " Error estimate: " << error_estimate << ", L1 norm: " << L1 << "\n";
74 }
75 // psi is orthogonal to its integer translates: \int \psi(x-k) \psi(x) \, \mathrm{d}x = 0
76 // g_n = 1/sqrt(2) <psi(t/2), phi(t-n)> (Mallat, 7.55)
77
78 // Now hit the boundary. Much can go wrong here; this just tests for segfaults:
79 int samples = 500;
80 Real xlo = a;
81 Real xhi = b;
82 for (int i = 0; i < samples; ++i)
83 {
84 CHECK_ULP_CLOSE(Real(0), psi(xlo), 0);
85 CHECK_ULP_CLOSE(Real(0), psi(xhi), 0);
86 if constexpr (p > 2)
87 {
88 CHECK_ULP_CLOSE(Real(0), psi.prime(xlo), 0);
89 CHECK_ULP_CLOSE(Real(0), psi.prime(xhi), 0);
90 if constexpr (p >= 6) {
91 CHECK_ULP_CLOSE(Real(0), psi.double_prime(xlo), 0);
92 CHECK_ULP_CLOSE(Real(0), psi.double_prime(xhi), 0);
93 }
94 }
95 xlo = std::nextafter(xlo, std::numeric_limits<Real>::lowest());
96 xhi = std::nextafter(xhi, std::numeric_limits<Real>::max());
97 }
98
99 xlo = a;
100 xhi = b;
101 for (int i = 0; i < samples; ++i) {
102 std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10);
103 assert(abs(psi(xlo)) <= 5);
104 assert(abs(psi(xhi)) <= 5);
105 if constexpr (p > 2)
106 {
107 assert(abs(psi.prime(xlo)) <= 5);
108 assert(abs(psi.prime(xhi)) <= 5);
109 if constexpr (p >= 6)
110 {
111 assert(abs(psi.double_prime(xlo)) <= 5);
112 assert(abs(psi.double_prime(xhi)) <= 5);
113 }
114 }
115 xlo = std::nextafter(xlo, std::numeric_limits<Real>::max());
116 xhi = std::nextafter(xhi, std::numeric_limits<Real>::lowest());
117 }
118 }
119
main()120 int main()
121 {
122 test_exact_value<double>();
123
124 boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i){
125 test_quadratures<float, i+3>();
126 test_quadratures<double, i+3>();
127 });
128
129 return boost::math::test::report_errors();
130 }
131