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