]>
Commit | Line | Data |
---|---|---|
f67539c2 TL |
1 | /* |
2 | * Copyright Nick Thompson, John Maddock 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> | |
1e59de90 | 14 | #include <boost/assert.hpp> |
f67539c2 TL |
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> | |
f67539c2 TL |
19 | #include <boost/math/special_functions/daubechies_scaling.hpp> |
20 | #include <boost/math/filters/daubechies.hpp> | |
21 | #include <boost/math/special_functions/detail/daubechies_scaling_integer_grid.hpp> | |
f67539c2 TL |
22 | #include <boost/math/quadrature/trapezoidal.hpp> |
23 | #include <boost/math/special_functions/next.hpp> | |
24 | ||
25 | #ifdef BOOST_HAS_FLOAT128 | |
26 | #include <boost/multiprecision/float128.hpp> | |
27 | using boost::multiprecision::float128; | |
28 | #endif | |
29 | ||
20effc67 | 30 | using std::sqrt; |
f67539c2 TL |
31 | // Mallat, Theorem 7.4, characterization number 3: |
32 | // A conjugate mirror filter has p vanishing moments iff h^{(n)}(pi) = 0 for 0 <= n < p. | |
33 | template<class Real, unsigned p> | |
34 | void test_daubechies_filters() | |
35 | { | |
36 | std::cout << "Testing Daubechies filters with " << p << " vanishing moments on type " << boost::core::demangle(typeid(Real).name()) << "\n"; | |
37 | Real tol = 3*std::numeric_limits<Real>::epsilon(); | |
38 | using boost::math::filters::daubechies_scaling_filter; | |
39 | using boost::math::filters::daubechies_wavelet_filter; | |
40 | ||
41 | auto h = daubechies_scaling_filter<Real, p>(); | |
42 | auto g = daubechies_wavelet_filter<Real, p>(); | |
43 | ||
44 | auto inner = std::inner_product(h.begin(), h.end(), g.begin(), Real(0)); | |
45 | CHECK_MOLLIFIED_CLOSE(0, inner, tol); | |
46 | ||
47 | // This is implied by Fourier transform of the two-scale dilatation equation; | |
48 | // If this doesn't hold, the infinite product for m_0 diverges. | |
49 | Real H0 = 0; | |
50 | for (size_t j = 0; j < h.size(); ++j) | |
51 | { | |
52 | H0 += h[j]; | |
53 | } | |
20effc67 | 54 | CHECK_MOLLIFIED_CLOSE(sqrt(static_cast<Real>(2)), H0, tol); |
f67539c2 TL |
55 | |
56 | // This is implied if we choose the scaling function to be an orthonormal basis of V0. | |
57 | Real scaling = 0; | |
58 | for (size_t j = 0; j < h.size(); ++j) { | |
59 | scaling += h[j]*h[j]; | |
60 | } | |
61 | CHECK_MOLLIFIED_CLOSE(1, scaling, tol); | |
62 | ||
63 | using std::pow; | |
64 | // Daubechies wavelet of order p has p vanishing moments. | |
65 | // Unfortunately, the condition number of the sum is infinite. | |
66 | // Hence we must scale the tolerance by the summation condition number to ensure that we don't get spurious test failures. | |
67 | for (size_t k = 1; k < p && k < 9; ++k) | |
68 | { | |
69 | Real hk = 0; | |
70 | Real abs_hk = 0; | |
71 | for (size_t n = 0; n < h.size(); ++n) | |
72 | { | |
73 | Real t = static_cast<Real>(pow(n, k)*h[n]); | |
74 | if (n & 1) | |
75 | { | |
76 | hk -= t; | |
77 | } | |
78 | else | |
79 | { | |
80 | hk += t; | |
81 | } | |
82 | abs_hk += abs(t); | |
83 | } | |
84 | // Multiply the tolerance by the condition number: | |
85 | Real cond = abs(hk) > 0 ? abs_hk/abs(hk) : 1/std::numeric_limits<Real>::epsilon(); | |
86 | if (!CHECK_MOLLIFIED_CLOSE(0, hk, 2*cond*tol)) | |
87 | { | |
88 | std::cerr << " The " << k << "th moment of the p = " << p << " filter did not vanish\n"; | |
89 | std::cerr << " Condition number = " << abs_hk/abs(hk) << "\n"; | |
90 | } | |
91 | } | |
92 | ||
93 | // For the scaling function to be orthonormal to its integer translates, | |
94 | // sum h_k h_{k-2l} = \delta_{0,l}. | |
95 | // See Theoretical Numerical Analysis, Atkinson, Exercise 4.5.2. | |
96 | // This is the last condition we could test to ensure that the filters are correct, | |
97 | // but I'm not gonna bother because it's painful! | |
98 | } | |
99 | ||
100 | // Test that the filters agree with Daubechies, Ten Lenctures on Wavelets, Table 6.1: | |
101 | void test_agreement_with_ten_lectures() | |
102 | { | |
103 | std::cout << "Testing agreement with Ten Lectures\n"; | |
104 | std::array<double, 4> h2 = {0.4829629131445341, 0.8365163037378077, 0.2241438680420134, -0.1294095225512603}; | |
105 | auto h2_ = boost::math::filters::daubechies_scaling_filter<double, 2>(); | |
106 | for (size_t i = 0; i < h2.size(); ++i) | |
107 | { | |
108 | CHECK_ULP_CLOSE(h2[i], h2_[i], 3); | |
109 | } | |
110 | ||
111 | std::array<double, 6> h3 = {0.3326705529500825, 0.8068915093110924, 0.4598775021184914, -0.1350110200102546, -0.0854412738820267, 0.0352262918857095}; | |
112 | auto h3_ = boost::math::filters::daubechies_scaling_filter<double, 3>(); | |
113 | for (size_t i = 0; i < h3.size(); ++i) | |
114 | { | |
115 | CHECK_ULP_CLOSE(h3[i], h3_[i], 5); | |
116 | } | |
117 | ||
118 | std::array<double, 8> h4 = {0.2303778133088964, 0.7148465705529154, 0.6308807679298587, -0.0279837694168599, -0.1870348117190931, 0.0308413818355607, 0.0328830116668852 , -0.010597401785069}; | |
119 | auto h4_ = boost::math::filters::daubechies_scaling_filter<double, 4>(); | |
120 | for (size_t i = 0; i < h4.size(); ++i) | |
121 | { | |
122 | if(!CHECK_ULP_CLOSE(h4[i], h4_[i], 18)) | |
123 | { | |
124 | std::cerr << " Index " << i << " incorrect.\n"; | |
125 | } | |
126 | } | |
127 | ||
128 | } | |
129 | ||
130 | template<class Real1, class Real2, size_t p> | |
131 | void test_filter_ulp_distance() | |
132 | { | |
133 | std::cout << "Testing filters ULP distance between types " | |
134 | << boost::core::demangle(typeid(Real1).name()) << "and" | |
135 | << boost::core::demangle(typeid(Real2).name()) << "\n"; | |
136 | using boost::math::filters::daubechies_scaling_filter; | |
137 | auto h1 = daubechies_scaling_filter<Real1, p>(); | |
138 | auto h2 = daubechies_scaling_filter<Real2, p>(); | |
139 | ||
140 | for (size_t i = 0; i < h1.size(); ++i) | |
141 | { | |
142 | if(!CHECK_ULP_CLOSE(h1[i], h2[i], 0)) | |
143 | { | |
144 | std::cerr << " Index " << i << " at order " << p << " failed tolerance check\n"; | |
145 | } | |
146 | } | |
147 | } | |
148 | ||
149 | ||
150 | template<class Real, unsigned p, unsigned order> | |
151 | void test_integer_grid() | |
152 | { | |
153 | std::cout << "Testing integer grid with " << p << " vanishing moments and " << order << " derivative on type " << boost::core::demangle(typeid(Real).name()) << "\n"; | |
154 | using boost::math::detail::daubechies_scaling_integer_grid; | |
155 | using boost::math::tools::summation_condition_number; | |
156 | Real unit_roundoff = std::numeric_limits<Real>::epsilon()/2; | |
157 | auto grid = daubechies_scaling_integer_grid<Real, p, order>(); | |
158 | ||
159 | if constexpr (order == 0) | |
160 | { | |
161 | auto cond = summation_condition_number<Real>(0); | |
162 | for (auto & x : grid) | |
163 | { | |
164 | cond += x; | |
165 | } | |
166 | CHECK_MOLLIFIED_CLOSE(1, cond.sum(), 6*cond.l1_norm()*unit_roundoff); | |
167 | } | |
168 | ||
169 | if constexpr (order == 1) | |
170 | { | |
171 | auto cond = summation_condition_number<Real>(0); | |
172 | for (size_t i = 0; i < grid.size(); ++i) { | |
173 | cond += i*grid[i]; | |
174 | } | |
175 | CHECK_MOLLIFIED_CLOSE(Real(-1), cond.sum(), 2*cond.l1_norm()*unit_roundoff); | |
176 | ||
177 | // Differentiate \sum_{k} \phi(x-k) = 1 to get this: | |
178 | cond = summation_condition_number<Real>(0); | |
179 | for (size_t i = 0; i < grid.size(); ++i) { | |
180 | cond += grid[i]; | |
181 | } | |
182 | CHECK_MOLLIFIED_CLOSE(Real(0), cond.sum(), 2*cond.l1_norm()*unit_roundoff); | |
183 | } | |
184 | ||
185 | if constexpr (order == 2) | |
186 | { | |
187 | auto cond = summation_condition_number<Real>(0); | |
188 | for (size_t i = 0; i < grid.size(); ++i) | |
189 | { | |
190 | cond += i*i*grid[i]; | |
191 | } | |
192 | CHECK_MOLLIFIED_CLOSE(Real(2), cond.sum(), 2*cond.l1_norm()*unit_roundoff); | |
193 | ||
194 | // Differentiate \sum_{k} \phi(x-k) = 1 to get this: | |
195 | cond = summation_condition_number<Real>(0); | |
196 | for (size_t i = 0; i < grid.size(); ++i) | |
197 | { | |
198 | cond += grid[i]; | |
199 | } | |
200 | CHECK_MOLLIFIED_CLOSE(Real(0), cond.sum(), 2*cond.l1_norm()*unit_roundoff); | |
201 | } | |
202 | ||
203 | if constexpr (order == 3) | |
204 | { | |
205 | auto cond = summation_condition_number<Real>(0); | |
206 | for (size_t i = 0; i < grid.size(); ++i) | |
207 | { | |
208 | cond += i*i*i*grid[i]; | |
209 | } | |
210 | CHECK_MOLLIFIED_CLOSE(Real(-6), cond.sum(), 2*cond.l1_norm()*unit_roundoff); | |
211 | ||
212 | // Differentiate \sum_{k} \phi(x-k) = 1 to get this: | |
213 | cond = summation_condition_number<Real>(0); | |
214 | for (size_t i = 0; i < grid.size(); ++i) | |
215 | { | |
216 | cond += grid[i]; | |
217 | } | |
218 | CHECK_MOLLIFIED_CLOSE(Real(0), cond.sum(), 2*cond.l1_norm()*unit_roundoff); | |
219 | } | |
220 | ||
221 | if constexpr (order == 4) | |
222 | { | |
223 | auto cond = summation_condition_number<Real>(0); | |
224 | for (size_t i = 0; i < grid.size(); ++i) | |
225 | { | |
226 | cond += i*i*i*i*grid[i]; | |
227 | } | |
228 | CHECK_MOLLIFIED_CLOSE(24, cond.sum(), 2*cond.l1_norm()*unit_roundoff); | |
229 | ||
230 | // Differentiate \sum_{k} \phi(x-k) = 1 to get this: | |
231 | cond = summation_condition_number<Real>(0); | |
232 | for (size_t i = 0; i < grid.size(); ++i) | |
233 | { | |
234 | cond += grid[i]; | |
235 | } | |
236 | CHECK_MOLLIFIED_CLOSE(Real(0), cond.sum(), 2*cond.l1_norm()*unit_roundoff); | |
237 | } | |
238 | } | |
239 | ||
240 | template<class Real> | |
241 | void test_dyadic_grid() | |
242 | { | |
243 | std::cout << "Testing dyadic grid on type " << boost::core::demangle(typeid(Real).name()) << "\n"; | |
244 | auto f = [&](auto i) | |
245 | { | |
246 | auto phijk = boost::math::daubechies_scaling_dyadic_grid<Real, i+2, 0>(0); | |
247 | auto phik = boost::math::detail::daubechies_scaling_integer_grid<Real, i+2, 0>(); | |
1e59de90 | 248 | BOOST_ASSERT(phik.size() == phijk.size()); |
f67539c2 TL |
249 | |
250 | for (size_t k = 0; k < phik.size(); ++k) | |
251 | { | |
252 | CHECK_ULP_CLOSE(phik[k], phijk[k], 0); | |
253 | } | |
254 | ||
255 | for (uint64_t j = 1; j < 10; ++j) | |
256 | { | |
257 | phijk = boost::math::daubechies_scaling_dyadic_grid<Real, i+2, 0>(j); | |
258 | phik = boost::math::detail::daubechies_scaling_integer_grid<Real, i+2, 0>(); | |
259 | for (uint64_t l = 0; l < static_cast<uint64_t>(phik.size()); ++l) | |
260 | { | |
261 | CHECK_ULP_CLOSE(phik[l], phijk[l*(uint64_t(1)<<j)], 0); | |
262 | } | |
263 | ||
264 | // This test is from Daubechies, Ten Lectures on Wavelets, Ch 7 "More About Compactly Supported Wavelets", | |
265 | // page 245: \forall y \in \mathbb{R}, \sum_{n \in \mathbb{Z}} \phi(y+n) = 1 | |
266 | for (size_t k = 1; k < j; ++k) | |
267 | { | |
268 | auto cond = boost::math::tools::summation_condition_number<Real>(0); | |
269 | for (uint64_t l = 0; l < static_cast<uint64_t>(phik.size()); ++l) | |
270 | { | |
271 | uint64_t idx = l*(uint64_t(1)<<j) + k; | |
272 | if (idx < phijk.size()) | |
273 | { | |
274 | cond += phijk[idx]; | |
275 | } | |
276 | } | |
277 | CHECK_MOLLIFIED_CLOSE(Real(1), cond.sum(), 10*cond()*std::numeric_limits<Real>::epsilon()); | |
278 | } | |
279 | } | |
280 | }; | |
281 | ||
282 | boost::hana::for_each(std::make_index_sequence<18>(), f); | |
283 | } | |
284 | ||
285 | ||
286 | // Taken from Lin, 2005, doi:10.1016/j.amc.2004.12.038, | |
287 | // "Direct algorithm for computation of derivatives of the Daubechies basis functions" | |
288 | void test_first_derivative() | |
289 | { | |
290 | auto phi1_3 = boost::math::detail::daubechies_scaling_integer_grid<long double, 3, 1>(); | |
291 | std::array<long double, 6> lin_3{0.0L, 1.638452340884085725014976L, -2.232758190463137395017742L, | |
292 | 0.5501593582740176149905562L, 0.04414649130503405501220997L, 0.0L}; | |
293 | for (size_t i = 0; i < lin_3.size(); ++i) | |
294 | { | |
295 | if(!CHECK_ULP_CLOSE(lin_3[i], phi1_3[i], 0)) | |
296 | { | |
297 | std::cerr << " Index " << i << " is incorrect\n"; | |
298 | } | |
299 | } | |
300 | ||
301 | auto phi1_4 = boost::math::detail::daubechies_scaling_integer_grid<long double, 4, 1>(); | |
302 | std::array<long double, 8> lin_4 = {0.0L, 1.776072007522184640093776L, -2.785349397229543142492785L, 1.192452536632278174347632L, | |
303 | -0.1313745151846729587935189L, -0.05357102822023923595359996L,0.001770396479992522798495351L, 0.0L}; | |
304 | ||
305 | for (size_t i = 0; i < lin_4.size(); ++i) | |
306 | { | |
307 | if(!CHECK_ULP_CLOSE(lin_4[i], phi1_4[i], 0)) | |
308 | { | |
309 | std::cerr << " Index " << i << " is incorrect\n"; | |
310 | } | |
311 | } | |
312 | ||
313 | std::array<long double, 10> lin_5 = {0.0L, 1.558326313047001366564379L, -2.436012783189551921436896L, 1.235905129801454293947039L, -0.3674377136938866359947561L, | |
314 | -0.02178035117564654658884556L,0.03234719350814368885815854L,-0.001335619912770701035229331L,-0.00001216838474354431384970525L,0.0L}; | |
315 | auto phi1_5 = boost::math::detail::daubechies_scaling_integer_grid<long double, 5, 1>(); | |
316 | for (size_t i = 0; i < lin_5.size(); ++i) | |
317 | { | |
318 | if(!CHECK_ULP_CLOSE(lin_5[i], phi1_5[i], 0)) | |
319 | { | |
320 | std::cerr << " Index " << i << " is incorrect\n"; | |
321 | } | |
322 | } | |
323 | } | |
324 | ||
325 | template<typename Real, int p> | |
326 | void test_quadratures() | |
327 | { | |
328 | std::cout << "Testing " << p << " vanishing moment scaling function quadratures on type " << boost::core::demangle(typeid(Real).name()) << "\n"; | |
329 | using boost::math::quadrature::trapezoidal; | |
330 | if constexpr (p == 2) | |
331 | { | |
332 | // 2phi is truly bizarre, because two successive trapezoidal estimates are always bitwise equal, | |
333 | // whereas the third is way different. I don' t think that's a reasonable thing to optimize for, | |
334 | // so one-off it is. | |
335 | Real h = Real(1)/Real(256); | |
336 | auto phi = boost::math::daubechies_scaling<Real, p>(); | |
337 | std::cout << "Scaling functor size is " << phi.bytes() << " bytes" << std::endl; | |
338 | Real t = 0; | |
339 | Real Q = 0; | |
340 | while (t < 3) { | |
341 | Q += phi(t); | |
342 | t += h; | |
343 | } | |
344 | Q *= h; | |
345 | CHECK_ULP_CLOSE(Real(1), Q, 32); | |
346 | ||
347 | auto [a, b] = phi.support(); | |
348 | // Now hit the boundary. Much can go wrong here; this just tests for segfaults: | |
349 | int samples = 500; | |
350 | Real xlo = a; | |
351 | Real xhi = b; | |
352 | for (int i = 0; i < samples; ++i) | |
353 | { | |
354 | CHECK_ULP_CLOSE(Real(0), phi(xlo), 0); | |
355 | CHECK_ULP_CLOSE(Real(0), phi(xhi), 0); | |
356 | xlo = std::nextafter(xlo, std::numeric_limits<Real>::lowest()); | |
1e59de90 | 357 | xhi = std::nextafter(xhi, (std::numeric_limits<Real>::max)()); |
f67539c2 TL |
358 | } |
359 | ||
360 | xlo = a; | |
361 | xhi = b; | |
362 | for (int i = 0; i < samples; ++i) { | |
1e59de90 TL |
363 | BOOST_ASSERT(abs(phi(xlo)) <= 5); |
364 | BOOST_ASSERT(abs(phi(xhi)) <= 5); | |
365 | xlo = std::nextafter(xlo, (std::numeric_limits<Real>::max)()); | |
f67539c2 TL |
366 | xhi = std::nextafter(xhi, std::numeric_limits<Real>::lowest()); |
367 | } | |
368 | ||
369 | return; | |
370 | } | |
371 | else if constexpr (p > 2) | |
372 | { | |
373 | auto phi = boost::math::daubechies_scaling<Real, p>(); | |
374 | std::cout << "Scaling functor size is " << phi.bytes() << " bytes" << std::endl; | |
375 | ||
376 | Real tol = std::numeric_limits<Real>::epsilon(); | |
377 | Real error_estimate = std::numeric_limits<Real>::quiet_NaN(); | |
378 | Real L1 = std::numeric_limits<Real>::quiet_NaN(); | |
379 | auto [a, b] = phi.support(); | |
380 | Real Q = trapezoidal(phi, a, b, tol, 15, &error_estimate, &L1); | |
381 | if (!CHECK_MOLLIFIED_CLOSE(Real(1), Q, Real(0.0001))) | |
382 | { | |
383 | std::cerr << " Quadrature of " << p << " vanishing moment scaling function is not equal 1.\n"; | |
384 | std::cerr << " Error estimate is " << error_estimate << ", L1 norm is " << L1 << "\n"; | |
385 | } | |
386 | ||
387 | auto phi_sq = [phi](Real x) { Real t = phi(x); return t*t; }; | |
388 | Q = trapezoidal(phi, a, b, tol, 15, &error_estimate, &L1); | |
389 | if (!CHECK_MOLLIFIED_CLOSE(Real(1), Q, 20*std::sqrt(std::numeric_limits<Real>::epsilon())/(p*p))) | |
390 | { | |
391 | std::cerr << " L2 norm of " << p << " vanishing moment scaling function is not equal 1.\n"; | |
392 | std::cerr << " Error estimate is " << error_estimate << ", L1 norm is " << L1 << "\n"; | |
393 | } | |
394 | ||
395 | std::random_device rd; | |
1e59de90 | 396 | Real t = static_cast<Real>(rd())/static_cast<Real>((rd.max)()); |
f67539c2 TL |
397 | Real S = phi(t); |
398 | Real dS = phi.prime(t); | |
399 | while (t < b) | |
400 | { | |
401 | t += 1; | |
402 | S += phi(t); | |
403 | dS += phi.prime(t); | |
404 | } | |
405 | ||
406 | if(!CHECK_ULP_CLOSE(Real(1), S, 64)) | |
407 | { | |
408 | std::cerr << " Normalizing sum for " << p << " vanishing moment scaling function is incorrect.\n"; | |
409 | } | |
410 | ||
411 | // The p = 3, 4 convergence rate is very slow, making this produce false positives: | |
412 | if constexpr(p > 4) | |
413 | { | |
414 | if(!CHECK_MOLLIFIED_CLOSE(Real(0), dS, 100*std::sqrt(std::numeric_limits<Real>::epsilon()))) | |
415 | { | |
416 | std::cerr << " Derivative of normalizing sum for " << p << " vanishing moment scaling function doesn't vanish.\n"; | |
417 | } | |
418 | } | |
419 | ||
420 | // Test boundary for segfaults: | |
421 | int samples = 500; | |
422 | Real xlo = a; | |
423 | Real xhi = b; | |
424 | for (int i = 0; i < samples; ++i) | |
425 | { | |
426 | CHECK_ULP_CLOSE(Real(0), phi(xlo), 0); | |
427 | CHECK_ULP_CLOSE(Real(0), phi(xhi), 0); | |
428 | if constexpr (p > 2) { | |
1e59de90 TL |
429 | BOOST_ASSERT(abs(phi.prime(xlo)) <= 5); |
430 | BOOST_ASSERT(abs(phi.prime(xhi)) <= 5); | |
f67539c2 | 431 | if constexpr (p > 5) { |
1e59de90 TL |
432 | BOOST_ASSERT(abs(phi.double_prime(xlo)) <= 5); |
433 | BOOST_ASSERT(abs(phi.double_prime(xhi)) <= 5); | |
f67539c2 TL |
434 | } |
435 | } | |
436 | xlo = std::nextafter(xlo, std::numeric_limits<Real>::lowest()); | |
1e59de90 | 437 | xhi = std::nextafter(xhi, (std::numeric_limits<Real>::max)()); |
f67539c2 TL |
438 | } |
439 | ||
440 | xlo = a; | |
441 | xhi = b; | |
442 | for (int i = 0; i < samples; ++i) { | |
1e59de90 TL |
443 | BOOST_ASSERT(abs(phi(xlo)) <= 5); |
444 | BOOST_ASSERT(abs(phi(xhi)) <= 5); | |
445 | xlo = std::nextafter(xlo, (std::numeric_limits<Real>::max)()); | |
f67539c2 TL |
446 | xhi = std::nextafter(xhi, std::numeric_limits<Real>::lowest()); |
447 | } | |
448 | } | |
449 | } | |
450 | ||
451 | int main() | |
452 | { | |
1e59de90 | 453 | #ifndef __MINGW32__ |
f67539c2 TL |
454 | boost::hana::for_each(std::make_index_sequence<18>(), [&](auto i){ |
455 | test_quadratures<float, i+2>(); | |
456 | test_quadratures<double, i+2>(); | |
457 | }); | |
458 | ||
459 | test_agreement_with_ten_lectures(); | |
460 | ||
461 | boost::hana::for_each(std::make_index_sequence<19>(), [&](auto i){ | |
462 | test_daubechies_filters<float, i+1>(); | |
463 | test_daubechies_filters<double, i+1>(); | |
464 | test_daubechies_filters<long double, i+1>(); | |
465 | }); | |
466 | ||
467 | test_first_derivative(); | |
468 | ||
469 | // All scaling functions have a first derivative. | |
470 | boost::hana::for_each(std::make_index_sequence<18>(), [&](auto idx){ | |
471 | test_integer_grid<float, idx+2, 0>(); | |
472 | test_integer_grid<float, idx+2, 1>(); | |
473 | test_integer_grid<double, idx+2, 0>(); | |
474 | test_integer_grid<double, idx+2, 1>(); | |
475 | test_integer_grid<long double, idx+2, 0>(); | |
476 | test_integer_grid<long double, idx+2, 1>(); | |
477 | #ifdef BOOST_HAS_FLOAT128 | |
478 | test_integer_grid<float128, idx+2, 0>(); | |
479 | test_integer_grid<float128, idx+2, 1>(); | |
480 | #endif | |
481 | }); | |
482 | ||
483 | // 4-tap (2 vanishing moment) scaling function does not have a second derivative; | |
484 | // all other scaling functions do. | |
485 | boost::hana::for_each(std::make_index_sequence<17>(), [&](auto idx){ | |
486 | test_integer_grid<float, idx+3, 2>(); | |
487 | test_integer_grid<double, idx+3, 2>(); | |
488 | test_integer_grid<long double, idx+3, 2>(); | |
489 | #ifdef BOOST_HAS_FLOAT128 | |
490 | test_integer_grid<boost::multiprecision::float128, idx+3, 2>(); | |
491 | #endif | |
492 | }); | |
493 | ||
494 | // 8-tap filter (4 vanishing moments) is the first to have a third derivative. | |
495 | boost::hana::for_each(std::make_index_sequence<16>(), [&](auto idx){ | |
496 | test_integer_grid<float, idx+4, 3>(); | |
497 | test_integer_grid<double, idx+4, 3>(); | |
498 | test_integer_grid<long double, idx+4, 3>(); | |
499 | #ifdef BOOST_HAS_FLOAT128 | |
500 | test_integer_grid<boost::multiprecision::float128, idx+4, 3>(); | |
501 | #endif | |
502 | }); | |
503 | ||
504 | // 10-tap filter (5 vanishing moments) is the first to have a fourth derivative. | |
505 | boost::hana::for_each(std::make_index_sequence<15>(), [&](auto idx){ | |
506 | test_integer_grid<float, idx+5, 4>(); | |
507 | test_integer_grid<double, idx+5, 4>(); | |
508 | test_integer_grid<long double, idx+5, 4>(); | |
509 | #ifdef BOOST_HAS_FLOAT128 | |
510 | test_integer_grid<boost::multiprecision::float128, idx+5, 4>(); | |
511 | #endif | |
512 | }); | |
513 | ||
514 | test_dyadic_grid<float>(); | |
515 | test_dyadic_grid<double>(); | |
516 | test_dyadic_grid<long double>(); | |
517 | #ifdef BOOST_HAS_FLOAT128 | |
518 | test_dyadic_grid<float128>(); | |
519 | #endif | |
520 | ||
521 | ||
522 | #ifdef BOOST_HAS_FLOAT128 | |
523 | boost::hana::for_each(std::make_index_sequence<19>(), [&](auto i){ | |
524 | test_filter_ulp_distance<float128, long double, i+1>(); | |
525 | test_filter_ulp_distance<float128, double, i+1>(); | |
526 | test_filter_ulp_distance<float128, float, i+1>(); | |
527 | }); | |
528 | ||
529 | boost::hana::for_each(std::make_index_sequence<19>(), [&](auto i){ | |
530 | test_daubechies_filters<float128, i+1>(); | |
531 | }); | |
532 | #endif | |
1e59de90 | 533 | #endif // compiler guard for CI |
f67539c2 TL |
534 | return boost::math::test::report_errors(); |
535 | } |