]>
Commit | Line | Data |
---|---|---|
f67539c2 TL |
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> | |
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 | ||
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 | } |