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 <cmath>
14 #include <boost/core/demangle.hpp>
15 #include <boost/hana/for_each.hpp>
16 #include <boost/hana/ext/std/integer_sequence.hpp>
17 #include <boost/math/quadrature/wavelet_transforms.hpp>
18 #include <boost/math/tools/minima.hpp>
19 #include <boost/math/quadrature/trapezoidal.hpp>
20
21 #ifdef BOOST_HAS_FLOAT128
22 #include <boost/multiprecision/float128.hpp>
23 using boost::multiprecision::float128;
24 #endif
25
26
27 using boost::math::constants::pi;
28 using boost::math::constants::root_two;
29 using boost::math::quadrature::daubechies_wavelet_transform;
30 using boost::math::quadrature::trapezoidal;
31
32 template<typename Real, int p>
test_wavelet_transform()33 void test_wavelet_transform()
34 {
35 std::cout << "Testing wavelet transform of " << p << " vanishing moment Daubechies wavelet on type " << boost::core::demangle(typeid(Real).name()) << "\n";
36 auto psi = boost::math::daubechies_wavelet<Real, p>();
37
38 auto abs_psi = [&psi](Real x) {
39 return abs(psi(x));
40 };
41 auto [a, b] = psi.support();
42 auto psil1 = trapezoidal(abs_psi, a, b);
43 Real psi_sup_norm = 0;
44 for (double x = a; x < b; x += 0.0001)
45 {
46 Real y = psi(x);
47 if (std::abs(y) > psi_sup_norm)
48 {
49 psi_sup_norm = std::abs(y);
50 }
51 }
52 std::cout << "psi sup norm = " << psi_sup_norm << "\n";
53 // An even function:
54 auto f = [](Real x) {
55 return std::exp(-abs(x));
56 };
57 Real fmax = 1;
58 Real fl2 = 1;
59 Real fl1 = 2;
60 std::cout << "||f||_1 = " << fl1 << "\n";
61 std::cout << "||f||_2 = " << fl2 << "\n";
62
63 auto Wf = daubechies_wavelet_transform(f, psi);
64 for (double s = 0; s < 10; s += 0.01)
65 {
66 Real w1 = Wf(s, 0.0);
67 Real w2 = Wf(-s, 0.0);
68 // Since f is an even function, we get w1 = w2:
69 CHECK_ULP_CLOSE(w1, w2, 12);
70 }
71
72 // The wavelet transform with respect to Daubechies wavelets
73 for (double s = -10; s < 10; s += 0.1)
74 {
75 for (double t = -10; t < 10; t+= 0.1)
76 {
77 Real w = Wf(s, t);
78 // Integral inequality:
79 Real r1 = sqrt(abs(s))*fmax*psil1;
80 if(abs(w) > r1)
81 {
82 std::cerr << " Integral inequality | int fg| <= ||f||_infty ||psi||_1 is violated.\n";
83 }
84 if (s != 0)
85 {
86 Real r2 = fl1*psi_sup_norm/sqrt(abs(s));
87 if(abs(w) > r2)
88 {
89 std::cerr << " Integral inequality | int fg| <= ||f||_1 ||psi||_infty/sqrt(|s|) is violated.\n";
90 std::cerr << " Violation: " << abs(w) << " !<= " << r2 << "\n";
91 }
92 Real r3 = fmax*psil1/sqrt(abs(s));
93 if(abs(w) > r3)
94 {
95 std::cerr << " Integral inequality | int fg| <= ||f||_infty ||psi||_1/sqrt(|s|) is violated.\n";
96 std::cerr << " Computed = " << abs(w) << ", expected " << r3 << "\n";
97 }
98 }
99 if (abs(w) > fl2)
100 {
101 std::cerr << " Integral inequality |f psi_s,t| <= ||f||_2 ||psi||2 violated.\n";
102 }
103 Real r4 = sqrt(abs(s))*fl1*psi_sup_norm;
104 if (abs(w) > r4)
105 {
106 std::cerr << " Integral inequality |W[f](s,t)| <= sqrt(|s|)||f||_1 ||psi||_infty is violated.\n";
107 }
108 Real r5 = sqrt(abs(s))*fmax*psil1;
109 if (abs(w) > r5)
110 {
111 std::cerr << " Integral inequality |W[f](s,t)| <= sqrt(|s|)||f||_infty ||psi||_1 is violated.\n";
112 }
113
114 }
115 }
116 }
117
main()118 int main()
119 {
120 boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i) {
121 test_wavelet_transform<double, i+3>();
122 });
123 return boost::math::test::report_errors();
124 }
125