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)
8 #include "math_unit_test.hpp"
14 #include <boost/assert.hpp>
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_scaling.hpp>
20 #include <boost/math/filters/daubechies.hpp>
21 #include <boost/math/special_functions/detail/daubechies_scaling_integer_grid.hpp>
22 #include <boost/math/quadrature/trapezoidal.hpp>
23 #include <boost/math/special_functions/next.hpp>
25 #ifdef BOOST_HAS_FLOAT128
26 #include <boost/multiprecision/float128.hpp>
27 using boost::multiprecision::float128
;
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()
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
;
41 auto h
= daubechies_scaling_filter
<Real
, p
>();
42 auto g
= daubechies_wavelet_filter
<Real
, p
>();
44 auto inner
= std::inner_product(h
.begin(), h
.end(), g
.begin(), Real(0));
45 CHECK_MOLLIFIED_CLOSE(0, inner
, tol
);
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.
50 for (size_t j
= 0; j
< h
.size(); ++j
)
54 CHECK_MOLLIFIED_CLOSE(sqrt(static_cast<Real
>(2)), H0
, tol
);
56 // This is implied if we choose the scaling function to be an orthonormal basis of V0.
58 for (size_t j
= 0; j
< h
.size(); ++j
) {
61 CHECK_MOLLIFIED_CLOSE(1, scaling
, tol
);
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
)
71 for (size_t n
= 0; n
< h
.size(); ++n
)
73 Real t
= static_cast<Real
>(pow(n
, k
)*h
[n
]);
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
))
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";
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!
100 // Test that the filters agree with Daubechies, Ten Lenctures on Wavelets, Table 6.1:
101 void test_agreement_with_ten_lectures()
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
)
108 CHECK_ULP_CLOSE(h2
[i
], h2_
[i
], 3);
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
)
115 CHECK_ULP_CLOSE(h3
[i
], h3_
[i
], 5);
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
)
122 if(!CHECK_ULP_CLOSE(h4
[i
], h4_
[i
], 18))
124 std::cerr
<< " Index " << i
<< " incorrect.\n";
130 template<class Real1
, class Real2
, size_t p
>
131 void test_filter_ulp_distance()
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
>();
140 for (size_t i
= 0; i
< h1
.size(); ++i
)
142 if(!CHECK_ULP_CLOSE(h1
[i
], h2
[i
], 0))
144 std::cerr
<< " Index " << i
<< " at order " << p
<< " failed tolerance check\n";
150 template<class Real
, unsigned p
, unsigned order
>
151 void test_integer_grid()
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
>();
159 if constexpr (order
== 0)
161 auto cond
= summation_condition_number
<Real
>(0);
162 for (auto & x
: grid
)
166 CHECK_MOLLIFIED_CLOSE(1, cond
.sum(), 6*cond
.l1_norm()*unit_roundoff
);
169 if constexpr (order
== 1)
171 auto cond
= summation_condition_number
<Real
>(0);
172 for (size_t i
= 0; i
< grid
.size(); ++i
) {
175 CHECK_MOLLIFIED_CLOSE(Real(-1), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
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
) {
182 CHECK_MOLLIFIED_CLOSE(Real(0), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
185 if constexpr (order
== 2)
187 auto cond
= summation_condition_number
<Real
>(0);
188 for (size_t i
= 0; i
< grid
.size(); ++i
)
192 CHECK_MOLLIFIED_CLOSE(Real(2), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
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
)
200 CHECK_MOLLIFIED_CLOSE(Real(0), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
203 if constexpr (order
== 3)
205 auto cond
= summation_condition_number
<Real
>(0);
206 for (size_t i
= 0; i
< grid
.size(); ++i
)
208 cond
+= i
*i
*i
*grid
[i
];
210 CHECK_MOLLIFIED_CLOSE(Real(-6), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
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
)
218 CHECK_MOLLIFIED_CLOSE(Real(0), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
221 if constexpr (order
== 4)
223 auto cond
= summation_condition_number
<Real
>(0);
224 for (size_t i
= 0; i
< grid
.size(); ++i
)
226 cond
+= i
*i
*i
*i
*grid
[i
];
228 CHECK_MOLLIFIED_CLOSE(24, cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
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
)
236 CHECK_MOLLIFIED_CLOSE(Real(0), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
241 void test_dyadic_grid()
243 std::cout
<< "Testing dyadic grid on type " << boost::core::demangle(typeid(Real
).name()) << "\n";
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>();
248 BOOST_ASSERT(phik
.size() == phijk
.size());
250 for (size_t k
= 0; k
< phik
.size(); ++k
)
252 CHECK_ULP_CLOSE(phik
[k
], phijk
[k
], 0);
255 for (uint64_t j
= 1; j
< 10; ++j
)
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
)
261 CHECK_ULP_CLOSE(phik
[l
], phijk
[l
*(uint64_t(1)<<j
)], 0);
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
)
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
)
271 uint64_t idx
= l
*(uint64_t(1)<<j
) + k
;
272 if (idx
< phijk
.size())
277 CHECK_MOLLIFIED_CLOSE(Real(1), cond
.sum(), 10*cond()*std::numeric_limits
<Real
>::epsilon());
282 boost::hana::for_each(std::make_index_sequence
<18>(), f
);
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()
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
)
295 if(!CHECK_ULP_CLOSE(lin_3
[i
], phi1_3
[i
], 0))
297 std::cerr
<< " Index " << i
<< " is incorrect\n";
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};
305 for (size_t i
= 0; i
< lin_4
.size(); ++i
)
307 if(!CHECK_ULP_CLOSE(lin_4
[i
], phi1_4
[i
], 0))
309 std::cerr
<< " Index " << i
<< " is incorrect\n";
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
)
318 if(!CHECK_ULP_CLOSE(lin_5
[i
], phi1_5
[i
], 0))
320 std::cerr
<< " Index " << i
<< " is incorrect\n";
325 template<typename Real
, int p
>
326 void test_quadratures()
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)
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,
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
;
345 CHECK_ULP_CLOSE(Real(1), Q
, 32);
347 auto [a
, b
] = phi
.support();
348 // Now hit the boundary. Much can go wrong here; this just tests for segfaults:
352 for (int i
= 0; i
< samples
; ++i
)
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());
357 xhi
= std::nextafter(xhi
, (std::numeric_limits
<Real
>::max
)());
362 for (int i
= 0; i
< samples
; ++i
) {
363 BOOST_ASSERT(abs(phi(xlo
)) <= 5);
364 BOOST_ASSERT(abs(phi(xhi
)) <= 5);
365 xlo
= std::nextafter(xlo
, (std::numeric_limits
<Real
>::max
)());
366 xhi
= std::nextafter(xhi
, std::numeric_limits
<Real
>::lowest());
371 else if constexpr (p
> 2)
373 auto phi
= boost::math::daubechies_scaling
<Real
, p
>();
374 std::cout
<< "Scaling functor size is " << phi
.bytes() << " bytes" << std::endl
;
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)))
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";
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
)))
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";
395 std::random_device rd
;
396 Real t
= static_cast<Real
>(rd())/static_cast<Real
>((rd
.max
)());
398 Real dS
= phi
.prime(t
);
406 if(!CHECK_ULP_CLOSE(Real(1), S
, 64))
408 std::cerr
<< " Normalizing sum for " << p
<< " vanishing moment scaling function is incorrect.\n";
411 // The p = 3, 4 convergence rate is very slow, making this produce false positives:
414 if(!CHECK_MOLLIFIED_CLOSE(Real(0), dS
, 100*std::sqrt(std::numeric_limits
<Real
>::epsilon())))
416 std::cerr
<< " Derivative of normalizing sum for " << p
<< " vanishing moment scaling function doesn't vanish.\n";
420 // Test boundary for segfaults:
424 for (int i
= 0; i
< samples
; ++i
)
426 CHECK_ULP_CLOSE(Real(0), phi(xlo
), 0);
427 CHECK_ULP_CLOSE(Real(0), phi(xhi
), 0);
428 if constexpr (p
> 2) {
429 BOOST_ASSERT(abs(phi
.prime(xlo
)) <= 5);
430 BOOST_ASSERT(abs(phi
.prime(xhi
)) <= 5);
431 if constexpr (p
> 5) {
432 BOOST_ASSERT(abs(phi
.double_prime(xlo
)) <= 5);
433 BOOST_ASSERT(abs(phi
.double_prime(xhi
)) <= 5);
436 xlo
= std::nextafter(xlo
, std::numeric_limits
<Real
>::lowest());
437 xhi
= std::nextafter(xhi
, (std::numeric_limits
<Real
>::max
)());
442 for (int i
= 0; i
< samples
; ++i
) {
443 BOOST_ASSERT(abs(phi(xlo
)) <= 5);
444 BOOST_ASSERT(abs(phi(xhi
)) <= 5);
445 xlo
= std::nextafter(xlo
, (std::numeric_limits
<Real
>::max
)());
446 xhi
= std::nextafter(xhi
, std::numeric_limits
<Real
>::lowest());
454 boost::hana::for_each(std::make_index_sequence
<18>(), [&](auto i
){
455 test_quadratures
<float, i
+2>();
456 test_quadratures
<double, i
+2>();
459 test_agreement_with_ten_lectures();
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>();
467 test_first_derivative();
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>();
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>();
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>();
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>();
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
>();
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>();
529 boost::hana::for_each(std::make_index_sequence
<19>(), [&](auto i
){
530 test_daubechies_filters
<float128
, i
+1>();
533 #endif // compiler guard for CI
534 return boost::math::test::report_errors();