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/core/demangle.hpp>
15 #include <boost/hana/for_each.hpp>
16 #include <boost/hana/ext/std/integer_sequence.hpp>
17 #include <boost/math/tools/condition_numbers.hpp>
18 #include <boost/math/special_functions/daubechies_scaling.hpp>
19 #include <boost/math/filters/daubechies.hpp>
20 #include <boost/math/special_functions/detail/daubechies_scaling_integer_grid.hpp>
21 #include <boost/math/quadrature/trapezoidal.hpp>
22 #include <boost/math/special_functions/next.hpp>
24 #ifdef BOOST_HAS_FLOAT128
25 #include <boost/multiprecision/float128.hpp>
26 using boost::multiprecision::float128
;
30 // Mallat, Theorem 7.4, characterization number 3:
31 // A conjugate mirror filter has p vanishing moments iff h^{(n)}(pi) = 0 for 0 <= n < p.
32 template<class Real
, unsigned p
>
33 void test_daubechies_filters()
35 std::cout
<< "Testing Daubechies filters with " << p
<< " vanishing moments on type " << boost::core::demangle(typeid(Real
).name()) << "\n";
36 Real tol
= 3*std::numeric_limits
<Real
>::epsilon();
37 using boost::math::filters::daubechies_scaling_filter
;
38 using boost::math::filters::daubechies_wavelet_filter
;
40 auto h
= daubechies_scaling_filter
<Real
, p
>();
41 auto g
= daubechies_wavelet_filter
<Real
, p
>();
43 auto inner
= std::inner_product(h
.begin(), h
.end(), g
.begin(), Real(0));
44 CHECK_MOLLIFIED_CLOSE(0, inner
, tol
);
46 // This is implied by Fourier transform of the two-scale dilatation equation;
47 // If this doesn't hold, the infinite product for m_0 diverges.
49 for (size_t j
= 0; j
< h
.size(); ++j
)
53 CHECK_MOLLIFIED_CLOSE(sqrt(static_cast<Real
>(2)), H0
, tol
);
55 // This is implied if we choose the scaling function to be an orthonormal basis of V0.
57 for (size_t j
= 0; j
< h
.size(); ++j
) {
60 CHECK_MOLLIFIED_CLOSE(1, scaling
, tol
);
63 // Daubechies wavelet of order p has p vanishing moments.
64 // Unfortunately, the condition number of the sum is infinite.
65 // Hence we must scale the tolerance by the summation condition number to ensure that we don't get spurious test failures.
66 for (size_t k
= 1; k
< p
&& k
< 9; ++k
)
70 for (size_t n
= 0; n
< h
.size(); ++n
)
72 Real t
= static_cast<Real
>(pow(n
, k
)*h
[n
]);
83 // Multiply the tolerance by the condition number:
84 Real cond
= abs(hk
) > 0 ? abs_hk
/abs(hk
) : 1/std::numeric_limits
<Real
>::epsilon();
85 if (!CHECK_MOLLIFIED_CLOSE(0, hk
, 2*cond
*tol
))
87 std::cerr
<< " The " << k
<< "th moment of the p = " << p
<< " filter did not vanish\n";
88 std::cerr
<< " Condition number = " << abs_hk
/abs(hk
) << "\n";
92 // For the scaling function to be orthonormal to its integer translates,
93 // sum h_k h_{k-2l} = \delta_{0,l}.
94 // See Theoretical Numerical Analysis, Atkinson, Exercise 4.5.2.
95 // This is the last condition we could test to ensure that the filters are correct,
96 // but I'm not gonna bother because it's painful!
99 // Test that the filters agree with Daubechies, Ten Lenctures on Wavelets, Table 6.1:
100 void test_agreement_with_ten_lectures()
102 std::cout
<< "Testing agreement with Ten Lectures\n";
103 std::array
<double, 4> h2
= {0.4829629131445341, 0.8365163037378077, 0.2241438680420134, -0.1294095225512603};
104 auto h2_
= boost::math::filters::daubechies_scaling_filter
<double, 2>();
105 for (size_t i
= 0; i
< h2
.size(); ++i
)
107 CHECK_ULP_CLOSE(h2
[i
], h2_
[i
], 3);
110 std::array
<double, 6> h3
= {0.3326705529500825, 0.8068915093110924, 0.4598775021184914, -0.1350110200102546, -0.0854412738820267, 0.0352262918857095};
111 auto h3_
= boost::math::filters::daubechies_scaling_filter
<double, 3>();
112 for (size_t i
= 0; i
< h3
.size(); ++i
)
114 CHECK_ULP_CLOSE(h3
[i
], h3_
[i
], 5);
117 std::array
<double, 8> h4
= {0.2303778133088964, 0.7148465705529154, 0.6308807679298587, -0.0279837694168599, -0.1870348117190931, 0.0308413818355607, 0.0328830116668852 , -0.010597401785069};
118 auto h4_
= boost::math::filters::daubechies_scaling_filter
<double, 4>();
119 for (size_t i
= 0; i
< h4
.size(); ++i
)
121 if(!CHECK_ULP_CLOSE(h4
[i
], h4_
[i
], 18))
123 std::cerr
<< " Index " << i
<< " incorrect.\n";
129 template<class Real1
, class Real2
, size_t p
>
130 void test_filter_ulp_distance()
132 std::cout
<< "Testing filters ULP distance between types "
133 << boost::core::demangle(typeid(Real1
).name()) << "and"
134 << boost::core::demangle(typeid(Real2
).name()) << "\n";
135 using boost::math::filters::daubechies_scaling_filter
;
136 auto h1
= daubechies_scaling_filter
<Real1
, p
>();
137 auto h2
= daubechies_scaling_filter
<Real2
, p
>();
139 for (size_t i
= 0; i
< h1
.size(); ++i
)
141 if(!CHECK_ULP_CLOSE(h1
[i
], h2
[i
], 0))
143 std::cerr
<< " Index " << i
<< " at order " << p
<< " failed tolerance check\n";
149 template<class Real
, unsigned p
, unsigned order
>
150 void test_integer_grid()
152 std::cout
<< "Testing integer grid with " << p
<< " vanishing moments and " << order
<< " derivative on type " << boost::core::demangle(typeid(Real
).name()) << "\n";
153 using boost::math::detail::daubechies_scaling_integer_grid
;
154 using boost::math::tools::summation_condition_number
;
155 Real unit_roundoff
= std::numeric_limits
<Real
>::epsilon()/2;
156 auto grid
= daubechies_scaling_integer_grid
<Real
, p
, order
>();
158 if constexpr (order
== 0)
160 auto cond
= summation_condition_number
<Real
>(0);
161 for (auto & x
: grid
)
165 CHECK_MOLLIFIED_CLOSE(1, cond
.sum(), 6*cond
.l1_norm()*unit_roundoff
);
168 if constexpr (order
== 1)
170 auto cond
= summation_condition_number
<Real
>(0);
171 for (size_t i
= 0; i
< grid
.size(); ++i
) {
174 CHECK_MOLLIFIED_CLOSE(Real(-1), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
176 // Differentiate \sum_{k} \phi(x-k) = 1 to get this:
177 cond
= summation_condition_number
<Real
>(0);
178 for (size_t i
= 0; i
< grid
.size(); ++i
) {
181 CHECK_MOLLIFIED_CLOSE(Real(0), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
184 if constexpr (order
== 2)
186 auto cond
= summation_condition_number
<Real
>(0);
187 for (size_t i
= 0; i
< grid
.size(); ++i
)
191 CHECK_MOLLIFIED_CLOSE(Real(2), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
193 // Differentiate \sum_{k} \phi(x-k) = 1 to get this:
194 cond
= summation_condition_number
<Real
>(0);
195 for (size_t i
= 0; i
< grid
.size(); ++i
)
199 CHECK_MOLLIFIED_CLOSE(Real(0), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
202 if constexpr (order
== 3)
204 auto cond
= summation_condition_number
<Real
>(0);
205 for (size_t i
= 0; i
< grid
.size(); ++i
)
207 cond
+= i
*i
*i
*grid
[i
];
209 CHECK_MOLLIFIED_CLOSE(Real(-6), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
211 // Differentiate \sum_{k} \phi(x-k) = 1 to get this:
212 cond
= summation_condition_number
<Real
>(0);
213 for (size_t i
= 0; i
< grid
.size(); ++i
)
217 CHECK_MOLLIFIED_CLOSE(Real(0), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
220 if constexpr (order
== 4)
222 auto cond
= summation_condition_number
<Real
>(0);
223 for (size_t i
= 0; i
< grid
.size(); ++i
)
225 cond
+= i
*i
*i
*i
*grid
[i
];
227 CHECK_MOLLIFIED_CLOSE(24, cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
229 // Differentiate \sum_{k} \phi(x-k) = 1 to get this:
230 cond
= summation_condition_number
<Real
>(0);
231 for (size_t i
= 0; i
< grid
.size(); ++i
)
235 CHECK_MOLLIFIED_CLOSE(Real(0), cond
.sum(), 2*cond
.l1_norm()*unit_roundoff
);
240 void test_dyadic_grid()
242 std::cout
<< "Testing dyadic grid on type " << boost::core::demangle(typeid(Real
).name()) << "\n";
245 auto phijk
= boost::math::daubechies_scaling_dyadic_grid
<Real
, i
+2, 0>(0);
246 auto phik
= boost::math::detail::daubechies_scaling_integer_grid
<Real
, i
+2, 0>();
247 assert(phik
.size() == phijk
.size());
249 for (size_t k
= 0; k
< phik
.size(); ++k
)
251 CHECK_ULP_CLOSE(phik
[k
], phijk
[k
], 0);
254 for (uint64_t j
= 1; j
< 10; ++j
)
256 phijk
= boost::math::daubechies_scaling_dyadic_grid
<Real
, i
+2, 0>(j
);
257 phik
= boost::math::detail::daubechies_scaling_integer_grid
<Real
, i
+2, 0>();
258 for (uint64_t l
= 0; l
< static_cast<uint64_t>(phik
.size()); ++l
)
260 CHECK_ULP_CLOSE(phik
[l
], phijk
[l
*(uint64_t(1)<<j
)], 0);
263 // This test is from Daubechies, Ten Lectures on Wavelets, Ch 7 "More About Compactly Supported Wavelets",
264 // page 245: \forall y \in \mathbb{R}, \sum_{n \in \mathbb{Z}} \phi(y+n) = 1
265 for (size_t k
= 1; k
< j
; ++k
)
267 auto cond
= boost::math::tools::summation_condition_number
<Real
>(0);
268 for (uint64_t l
= 0; l
< static_cast<uint64_t>(phik
.size()); ++l
)
270 uint64_t idx
= l
*(uint64_t(1)<<j
) + k
;
271 if (idx
< phijk
.size())
276 CHECK_MOLLIFIED_CLOSE(Real(1), cond
.sum(), 10*cond()*std::numeric_limits
<Real
>::epsilon());
281 boost::hana::for_each(std::make_index_sequence
<18>(), f
);
285 // Taken from Lin, 2005, doi:10.1016/j.amc.2004.12.038,
286 // "Direct algorithm for computation of derivatives of the Daubechies basis functions"
287 void test_first_derivative()
289 auto phi1_3
= boost::math::detail::daubechies_scaling_integer_grid
<long double, 3, 1>();
290 std::array
<long double, 6> lin_3
{0.0L, 1.638452340884085725014976L, -2.232758190463137395017742L,
291 0.5501593582740176149905562L, 0.04414649130503405501220997L, 0.0L};
292 for (size_t i
= 0; i
< lin_3
.size(); ++i
)
294 if(!CHECK_ULP_CLOSE(lin_3
[i
], phi1_3
[i
], 0))
296 std::cerr
<< " Index " << i
<< " is incorrect\n";
300 auto phi1_4
= boost::math::detail::daubechies_scaling_integer_grid
<long double, 4, 1>();
301 std::array
<long double, 8> lin_4
= {0.0L, 1.776072007522184640093776L, -2.785349397229543142492785L, 1.192452536632278174347632L,
302 -0.1313745151846729587935189L, -0.05357102822023923595359996L,0.001770396479992522798495351L, 0.0L};
304 for (size_t i
= 0; i
< lin_4
.size(); ++i
)
306 if(!CHECK_ULP_CLOSE(lin_4
[i
], phi1_4
[i
], 0))
308 std::cerr
<< " Index " << i
<< " is incorrect\n";
312 std::array
<long double, 10> lin_5
= {0.0L, 1.558326313047001366564379L, -2.436012783189551921436896L, 1.235905129801454293947039L, -0.3674377136938866359947561L,
313 -0.02178035117564654658884556L,0.03234719350814368885815854L,-0.001335619912770701035229331L,-0.00001216838474354431384970525L,0.0L};
314 auto phi1_5
= boost::math::detail::daubechies_scaling_integer_grid
<long double, 5, 1>();
315 for (size_t i
= 0; i
< lin_5
.size(); ++i
)
317 if(!CHECK_ULP_CLOSE(lin_5
[i
], phi1_5
[i
], 0))
319 std::cerr
<< " Index " << i
<< " is incorrect\n";
324 template<typename Real
, int p
>
325 void test_quadratures()
327 std::cout
<< "Testing " << p
<< " vanishing moment scaling function quadratures on type " << boost::core::demangle(typeid(Real
).name()) << "\n";
328 using boost::math::quadrature::trapezoidal
;
329 if constexpr (p
== 2)
331 // 2phi is truly bizarre, because two successive trapezoidal estimates are always bitwise equal,
332 // whereas the third is way different. I don' t think that's a reasonable thing to optimize for,
334 Real h
= Real(1)/Real(256);
335 auto phi
= boost::math::daubechies_scaling
<Real
, p
>();
336 std::cout
<< "Scaling functor size is " << phi
.bytes() << " bytes" << std::endl
;
344 CHECK_ULP_CLOSE(Real(1), Q
, 32);
346 auto [a
, b
] = phi
.support();
347 // Now hit the boundary. Much can go wrong here; this just tests for segfaults:
351 for (int i
= 0; i
< samples
; ++i
)
353 CHECK_ULP_CLOSE(Real(0), phi(xlo
), 0);
354 CHECK_ULP_CLOSE(Real(0), phi(xhi
), 0);
355 xlo
= std::nextafter(xlo
, std::numeric_limits
<Real
>::lowest());
356 xhi
= std::nextafter(xhi
, std::numeric_limits
<Real
>::max());
361 for (int i
= 0; i
< samples
; ++i
) {
362 assert(abs(phi(xlo
)) <= 5);
363 assert(abs(phi(xhi
)) <= 5);
364 xlo
= std::nextafter(xlo
, std::numeric_limits
<Real
>::max());
365 xhi
= std::nextafter(xhi
, std::numeric_limits
<Real
>::lowest());
370 else if constexpr (p
> 2)
372 auto phi
= boost::math::daubechies_scaling
<Real
, p
>();
373 std::cout
<< "Scaling functor size is " << phi
.bytes() << " bytes" << std::endl
;
375 Real tol
= std::numeric_limits
<Real
>::epsilon();
376 Real error_estimate
= std::numeric_limits
<Real
>::quiet_NaN();
377 Real L1
= std::numeric_limits
<Real
>::quiet_NaN();
378 auto [a
, b
] = phi
.support();
379 Real Q
= trapezoidal(phi
, a
, b
, tol
, 15, &error_estimate
, &L1
);
380 if (!CHECK_MOLLIFIED_CLOSE(Real(1), Q
, Real(0.0001)))
382 std::cerr
<< " Quadrature of " << p
<< " vanishing moment scaling function is not equal 1.\n";
383 std::cerr
<< " Error estimate is " << error_estimate
<< ", L1 norm is " << L1
<< "\n";
386 auto phi_sq
= [phi
](Real x
) { Real t
= phi(x
); return t
*t
; };
387 Q
= trapezoidal(phi
, a
, b
, tol
, 15, &error_estimate
, &L1
);
388 if (!CHECK_MOLLIFIED_CLOSE(Real(1), Q
, 20*std::sqrt(std::numeric_limits
<Real
>::epsilon())/(p
*p
)))
390 std::cerr
<< " L2 norm of " << p
<< " vanishing moment scaling function is not equal 1.\n";
391 std::cerr
<< " Error estimate is " << error_estimate
<< ", L1 norm is " << L1
<< "\n";
394 std::random_device rd
;
395 Real t
= static_cast<Real
>(rd())/static_cast<Real
>(rd
.max());
397 Real dS
= phi
.prime(t
);
405 if(!CHECK_ULP_CLOSE(Real(1), S
, 64))
407 std::cerr
<< " Normalizing sum for " << p
<< " vanishing moment scaling function is incorrect.\n";
410 // The p = 3, 4 convergence rate is very slow, making this produce false positives:
413 if(!CHECK_MOLLIFIED_CLOSE(Real(0), dS
, 100*std::sqrt(std::numeric_limits
<Real
>::epsilon())))
415 std::cerr
<< " Derivative of normalizing sum for " << p
<< " vanishing moment scaling function doesn't vanish.\n";
419 // Test boundary for segfaults:
423 for (int i
= 0; i
< samples
; ++i
)
425 CHECK_ULP_CLOSE(Real(0), phi(xlo
), 0);
426 CHECK_ULP_CLOSE(Real(0), phi(xhi
), 0);
427 if constexpr (p
> 2) {
428 assert(abs(phi
.prime(xlo
)) <= 5);
429 assert(abs(phi
.prime(xhi
)) <= 5);
430 if constexpr (p
> 5) {
431 assert(abs(phi
.double_prime(xlo
)) <= 5);
432 assert(abs(phi
.double_prime(xhi
)) <= 5);
435 xlo
= std::nextafter(xlo
, std::numeric_limits
<Real
>::lowest());
436 xhi
= std::nextafter(xhi
, std::numeric_limits
<Real
>::max());
441 for (int i
= 0; i
< samples
; ++i
) {
442 assert(abs(phi(xlo
)) <= 5);
443 assert(abs(phi(xhi
)) <= 5);
444 xlo
= std::nextafter(xlo
, std::numeric_limits
<Real
>::max());
445 xhi
= std::nextafter(xhi
, std::numeric_limits
<Real
>::lowest());
452 boost::hana::for_each(std::make_index_sequence
<18>(), [&](auto i
){
453 test_quadratures
<float, i
+2>();
454 test_quadratures
<double, i
+2>();
457 test_agreement_with_ten_lectures();
459 boost::hana::for_each(std::make_index_sequence
<19>(), [&](auto i
){
460 test_daubechies_filters
<float, i
+1>();
461 test_daubechies_filters
<double, i
+1>();
462 test_daubechies_filters
<long double, i
+1>();
465 test_first_derivative();
467 // All scaling functions have a first derivative.
468 boost::hana::for_each(std::make_index_sequence
<18>(), [&](auto idx
){
469 test_integer_grid
<float, idx
+2, 0>();
470 test_integer_grid
<float, idx
+2, 1>();
471 test_integer_grid
<double, idx
+2, 0>();
472 test_integer_grid
<double, idx
+2, 1>();
473 test_integer_grid
<long double, idx
+2, 0>();
474 test_integer_grid
<long double, idx
+2, 1>();
475 #ifdef BOOST_HAS_FLOAT128
476 test_integer_grid
<float128
, idx
+2, 0>();
477 test_integer_grid
<float128
, idx
+2, 1>();
481 // 4-tap (2 vanishing moment) scaling function does not have a second derivative;
482 // all other scaling functions do.
483 boost::hana::for_each(std::make_index_sequence
<17>(), [&](auto idx
){
484 test_integer_grid
<float, idx
+3, 2>();
485 test_integer_grid
<double, idx
+3, 2>();
486 test_integer_grid
<long double, idx
+3, 2>();
487 #ifdef BOOST_HAS_FLOAT128
488 test_integer_grid
<boost::multiprecision::float128
, idx
+3, 2>();
492 // 8-tap filter (4 vanishing moments) is the first to have a third derivative.
493 boost::hana::for_each(std::make_index_sequence
<16>(), [&](auto idx
){
494 test_integer_grid
<float, idx
+4, 3>();
495 test_integer_grid
<double, idx
+4, 3>();
496 test_integer_grid
<long double, idx
+4, 3>();
497 #ifdef BOOST_HAS_FLOAT128
498 test_integer_grid
<boost::multiprecision::float128
, idx
+4, 3>();
502 // 10-tap filter (5 vanishing moments) is the first to have a fourth derivative.
503 boost::hana::for_each(std::make_index_sequence
<15>(), [&](auto idx
){
504 test_integer_grid
<float, idx
+5, 4>();
505 test_integer_grid
<double, idx
+5, 4>();
506 test_integer_grid
<long double, idx
+5, 4>();
507 #ifdef BOOST_HAS_FLOAT128
508 test_integer_grid
<boost::multiprecision::float128
, idx
+5, 4>();
512 test_dyadic_grid
<float>();
513 test_dyadic_grid
<double>();
514 test_dyadic_grid
<long double>();
515 #ifdef BOOST_HAS_FLOAT128
516 test_dyadic_grid
<float128
>();
520 #ifdef BOOST_HAS_FLOAT128
521 boost::hana::for_each(std::make_index_sequence
<19>(), [&](auto i
){
522 test_filter_ulp_distance
<float128
, long double, i
+1>();
523 test_filter_ulp_distance
<float128
, double, i
+1>();
524 test_filter_ulp_distance
<float128
, float, i
+1>();
527 boost::hana::for_each(std::make_index_sequence
<19>(), [&](auto i
){
528 test_daubechies_filters
<float128
, i
+1>();
532 return boost::math::test::report_errors();