boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-27 19:12:08 +00:00
Code Issues Projects Releases Wiki Activity
Files
4a2f1febf8fc9255cd424f28f79a2dce41569e2a
math/test/daubechies_scaling_test.cpp
NAThompson b5188cfaba Continuous Daubechies wavelets: Fight bitrot.
2020-01-26 11:27:28 +08:00

482 lines
19 KiB
C++
Raw Blame History

/*
* Copyright Nick Thompson, 2019
* Use, modification and distribution are subject to the
* Boost Software License, Version 1.0. (See accompanying file
* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#include "math_unit_test.hpp"
#include <numeric>
#include <utility>
#include <boost/core/demangle.hpp>
#include <boost/hana/for_each.hpp>
#include <boost/hana/ext/std/integer_sequence.hpp>
#include <boost/math/tools/condition_numbers.hpp>
#include <boost/math/differentiation/finite_difference.hpp>
#include <boost/math/special_functions/daubechies_scaling.hpp>
#include <boost/math/filters/daubechies.hpp>
#include <boost/math/special_functions/detail/daubechies_scaling_integer_grid.hpp>
#include <boost/math/constants/constants.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;
#endif
using boost::math::constants::pi;
using boost::math::constants::root_two;
// Mallat, Theorem 7.4, characterization number 3:
// A conjugate mirror filter has p vanishing moments iff h^{(n)}(pi) = 0 for 0 <= n < p.
template<class Real, unsigned p>
void test_daubechies_filters()
{
std::cout << "Testing Daubechies filters with " << p << " vanishing moments on type " << boost::core::demangle(typeid(Real).name()) << "\n";
Real tol = 3*std::numeric_limits<Real>::epsilon();
using boost::math::filters::daubechies_scaling_filter;
using boost::math::filters::daubechies_wavelet_filter;
auto h = daubechies_scaling_filter<Real, p>();
auto g = daubechies_wavelet_filter<Real, p>();
auto inner = std::inner_product(h.begin(), h.end(), g.begin(), Real(0));
CHECK_MOLLIFIED_CLOSE(0, inner, tol);
// This is implied by Fourier transform of the two-scale dilatation equation;
// If this doesn't hold, the infinite product for m_0 diverges.
Real H0 = 0;
for (size_t j = 0; j < h.size(); ++j)
{
H0 += h[j];
}
CHECK_MOLLIFIED_CLOSE(root_two<Real>(), H0, tol);
// This is implied if we choose the scaling function to be an orthonormal basis of V0.
Real scaling = 0;
for (size_t j = 0; j < h.size(); ++j) {
scaling += h[j]*h[j];
}
CHECK_MOLLIFIED_CLOSE(1, scaling, tol);
using std::pow;
// Daubechies wavelet of order p has p vanishing moments.
// Unfortunately, the condition number of the sum is infinite.
// Hence we must scale the tolerance by the summation condition number to ensure that we don't get spurious test failures.
for (size_t k = 1; k < p; ++k)
{
Real hk = 0;
Real abs_hk = 0;
for (size_t n = 0; n < h.size(); ++n)
{
Real t = pow(n, k)*h[n];
if (n & 1)
{
hk -= t;
}
else
{
hk += t;
}
abs_hk += abs(t);
}
// Multiply the tolerance by the condition number:
Real cond = abs(hk) > 0 ? abs_hk/abs(hk) : 1/std::numeric_limits<Real>::epsilon();
if (!CHECK_MOLLIFIED_CLOSE(0, hk, cond*tol))
{
std::cerr << " The " << k << "th moment of the p = " << p << " filter did not vanish\n";
std::cerr << " Condition number = " << abs_hk/abs(hk) << "\n";
}
}
// For the scaling function to be orthonormal to its integer translates,
// sum h_k h_{k-2l} = \delta_{0,l}.
// See Theoretical Numerical Analysis, Atkinson, Exercise 4.5.2.
// This is the last condition we could test to ensure that the filters are correct,
// but I'm not gonna bother because it's painful!
}
// Test that the filters agree with Daubechies, Ten Lenctures on Wavelets, Table 6.1:
void test_agreement_with_ten_lectures()
{
std::cout << "Testing agreement with Ten Lectures\n";
std::array<double, 4> h2 = {0.4829629131445341, 0.8365163037378077, 0.2241438680420134, -0.1294095225512603};
auto h2_ = boost::math::filters::daubechies_scaling_filter<double, 2>();
for (size_t i = 0; i < h2.size(); ++i)
{
CHECK_ULP_CLOSE(h2[i], h2_[i], 3);
}
std::array<double, 6> h3 = {0.3326705529500825, 0.8068915093110924, 0.4598775021184914, -0.1350110200102546, -0.0854412738820267, 0.0352262918857095};
auto h3_ = boost::math::filters::daubechies_scaling_filter<double, 3>();
for (size_t i = 0; i < h3.size(); ++i)
{
CHECK_ULP_CLOSE(h3[i], h3_[i], 5);
}
std::array<double, 8> h4 = {0.2303778133088964, 0.7148465705529154, 0.6308807679298587, -0.0279837694168599, -0.1870348117190931, 0.0308413818355607, 0.0328830116668852 , -0.010597401785069};
auto h4_ = boost::math::filters::daubechies_scaling_filter<double, 4>();
for (size_t i = 0; i < h4.size(); ++i)
{
if(!CHECK_ULP_CLOSE(h4[i], h4_[i], 18)) {
std::cerr << " Index " << i << " incorrect.\n";
}
}
}
template<class Real1, class Real2, size_t p>
void test_filter_ulp_distance()
{
std::cout << "Testing filters ULP distance between types "
<< boost::core::demangle(typeid(Real1).name()) << "and"
<< boost::core::demangle(typeid(Real2).name()) << "\n";
using boost::math::filters::daubechies_scaling_filter;
auto h1 = daubechies_scaling_filter<Real1, p>();
auto h2 = daubechies_scaling_filter<Real2, p>();
for (size_t i = 0; i < h1.size(); ++i)
{
if(!CHECK_ULP_CLOSE(h1[i], h2[i], 0))
{
std::cerr << " Index " << i << " at order " << p << " failed tolerance check\n";
}
}
}
template<class Real, unsigned p, unsigned order>
void test_integer_grid()
{
std::cout << "Testing integer grid with " << p << " vanishing moments and " << order << " derivative on type " << boost::core::demangle(typeid(Real).name()) << "\n";
using boost::math::detail::daubechies_scaling_integer_grid;
using boost::math::tools::summation_condition_number;
Real unit_roundoff = std::numeric_limits<Real>::epsilon()/2;
auto grid = daubechies_scaling_integer_grid<Real, p, order>();
if constexpr (order == 0) {
auto cond = summation_condition_number<Real>(0);
for (auto & x : grid) {
cond += x;
}
CHECK_MOLLIFIED_CLOSE(1, cond.sum(), 6*cond.l1_norm()*unit_roundoff);
}
if constexpr (order == 1) {
auto cond = summation_condition_number<Real>(0);
for (size_t i = 0; i < grid.size(); ++i) {
cond += i*grid[i];
}
CHECK_MOLLIFIED_CLOSE(-1, cond.sum(), 2*cond.l1_norm()*unit_roundoff);
// Differentiate \sum_{k} \phi(x-k) = 1 to get this:
cond = summation_condition_number<Real>(0);
for (size_t i = 0; i < grid.size(); ++i) {
cond += grid[i];
}
CHECK_MOLLIFIED_CLOSE(0, cond.sum(), 2*cond.l1_norm()*unit_roundoff);
}
if constexpr (order == 2) {
auto cond = summation_condition_number<Real>(0);
for (size_t i = 0; i < grid.size(); ++i) {
cond += i*i*grid[i];
}
CHECK_MOLLIFIED_CLOSE(2, cond.sum(), 2*cond.l1_norm()*unit_roundoff);
// Differentiate \sum_{k} \phi(x-k) = 1 to get this:
cond = summation_condition_number<Real>(0);
for (size_t i = 0; i < grid.size(); ++i) {
cond += grid[i];
}
CHECK_MOLLIFIED_CLOSE(0, cond.sum(), 2*cond.l1_norm()*unit_roundoff);
}
if constexpr (order == 3) {
auto cond = summation_condition_number<Real>(0);
for (size_t i = 0; i < grid.size(); ++i) {
cond += i*i*i*grid[i];
}
CHECK_MOLLIFIED_CLOSE(-6, cond.sum(), 2*cond.l1_norm()*unit_roundoff);
// Differentiate \sum_{k} \phi(x-k) = 1 to get this:
cond = summation_condition_number<Real>(0);
for (size_t i = 0; i < grid.size(); ++i) {
cond += grid[i];
}
CHECK_MOLLIFIED_CLOSE(0, cond.sum(), 2*cond.l1_norm()*unit_roundoff);
}
if constexpr (order == 4) {
auto cond = summation_condition_number<Real>(0);
for (size_t i = 0; i < grid.size(); ++i) {
cond += i*i*i*i*grid[i];
}
CHECK_MOLLIFIED_CLOSE(24, cond.sum(), 2*cond.l1_norm()*unit_roundoff);
// Differentiate \sum_{k} \phi(x-k) = 1 to get this:
cond = summation_condition_number<Real>(0);
for (size_t i = 0; i < grid.size(); ++i) {
cond += grid[i];
}
CHECK_MOLLIFIED_CLOSE(0, cond.sum(), 2*cond.l1_norm()*unit_roundoff);
}
}
template<class Real>
void test_dyadic_grid()
{
std::cout << "Testing dyadic grid on type " << boost::core::demangle(typeid(Real).name()) << "\n";
boost::hana::for_each(std::make_index_sequence<13>(), [&](auto i){
auto phijk = boost::math::detail::dyadic_grid<Real, i+2, 0>(0);
auto phik = boost::math::detail::daubechies_scaling_integer_grid<Real, i+2, 0>();
assert(phik.size() == phijk.size());
for (size_t k = 0; k < phik.size(); ++k)
{
CHECK_ULP_CLOSE(phik[k], phijk[k], 0);
}
for (size_t j = 1; j < 10; ++j)
{
auto phijk = boost::math::detail::dyadic_grid<Real, i+2, 0>(j);
auto phik = boost::math::detail::daubechies_scaling_integer_grid<Real, i+2, 0>();
for (size_t i = 0; i < phik.size(); ++i)
{
CHECK_ULP_CLOSE(phik[i], phijk[i*(1<<j)], 0);
}
// This test is from Daubechies, Ten Lectures on Wavelets, Ch 7 "More About Compactly Supported Wavelets",
// page 245: \forall y \in \mathbb{R}, \sum_{n \in \mathbb{Z}} \phi(y+n) = 1
for (size_t k = 1; k < j; ++k)
{
auto cond = boost::math::tools::summation_condition_number<Real>(0);
for (size_t i = 0; i < phik.size(); ++i)
{
size_t idx = i*(1<<j) + k;
if (idx < phijk.size())
{
cond += phijk[idx];
}
}
CHECK_MOLLIFIED_CLOSE(1, cond.sum(), 10*cond()*std::numeric_limits<Real>::epsilon());
}
}
});
}
template<class Real>
void test_interpolation()
{
std::cout << "Testing constant interpolation on type " << boost::core::demangle(typeid(Real).name()) << "\n";
boost::hana::for_each(std::make_index_sequence<13>(), [&](auto i){
auto phik = boost::math::detail::daubechies_scaling_integer_grid<Real, i+2, 0>();
for (size_t j = 0; j < 5; ++j) {
auto phi = boost::math::daubechies_scaling<Real, i+2>(j);
assert(phik.size()==phi.support().second + 1);
for (size_t k = 1; k < phik.size(); ++k) {
auto expected = phik[k];
auto computed = phi.constant_interpolation(k);
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " Constant interpolation wrong at x = " << k << ", j_max = " << j << ", p = " << i+2 << "\n";
}
computed = phi.single_crank_linear(k);
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " Single crank linear interpolation wrong at x = " << k << ", j_max = " << j << ", p = " << i+2 << "\n";
}
computed = phi.first_order_taylor(k);
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " First order Taylor expansion wrong at x = " << k << ", j_max = " << j << ", p = " << i+2 << "\n";
}
computed = phi.constant_interpolation(k*(1+2*std::numeric_limits<Real>::epsilon()));
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " Constant interpolation wrong at x = " << k << ", j_max = " << j << ", p = " << i+2 << "\n";
}
computed = phi.constant_interpolation(k*(1-2*std::numeric_limits<Real>::epsilon()));
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " Constant interpolation wrong at x = " << k << ", j_max = " << j << ", p = " << i+2 << "\n";
}
computed = phi.linear_interpolation(k);
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " Linear interpolation wrong at x = " << k << ", j_max = " << j << ", p = " << i+2 << "\n";
}
}
for (size_t i = 0; i < phi.size() -1; ++i) {
Real x = phi.index_to_abscissa(i);
Real expected = phi[i];
Real computed = phi.constant_interpolation(x);
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " Constant interpolation wrong at x = " << x << ", j_max = " << j << ", p = " << i+2 << "\n";
}
computed = phi.linear_interpolation(x);
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " Linear interpolation wrong at x = " << x << ", j_max = " << j << ", p = " << i+2 << "\n";
}
computed = phi.single_crank_linear(x);
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " Single crank linear interpolation wrong at x = " << x << ", j_max = " << j << ", p = " << i+2 << "\n";
}
computed = phi.first_order_taylor(x);
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " First order Taylor expansion wrong at x = " << x << ", j_max = " << j << ", p = " << i+2 << "\n";
}
x += phi.spacing()/2;
computed = phi.linear_interpolation(x);
expected = phi[i]/2 + phi[i+1]/2;
if (!CHECK_ULP_CLOSE(expected, computed, 1)) {
std::cerr << " Linear interpolation wrong at x = " << x << ", j_max = " << j << ", p = " << i+2 << "\n";
}
x *= (1+std::numeric_limits<Real>::epsilon());
computed = phi.constant_interpolation(x);
expected = phi[i+1];
if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
std::cerr << " Linear interpolation wrong at x = " << x << ", j_max = " << j << ", p = " << i+2 << "\n";
}
}
}
});
}
// Taken from Lin, 2005, doi:10.1016/j.amc.2004.12.038,
// "Direct algorithm for computation of derivatives of the Daubechies basis functions"
void test_first_derivative()
{
auto phi1_3 = boost::math::detail::daubechies_scaling_integer_grid<long double, 3, 1>();
std::array<long double, 6> lin_3{0.0L, 1.638452340884085725014976L, -2.232758190463137395017742L, 0.5501593582740176149905562L, 0.04414649130503405501220997L, 0.0L};
for (size_t i = 0; i < lin_3.size(); ++i)
{
if(!CHECK_ULP_CLOSE(lin_3[i], phi1_3[i], 0))
{
std::cerr << " Index " << i << " is incorrect\n";
}
}
auto phi1_4 = boost::math::detail::daubechies_scaling_integer_grid<long double, 4, 1>();
std::array<long double, 8> lin_4 = {0.0L, 1.776072007522184640093776L, -2.785349397229543142492785L, 1.192452536632278174347632L, -0.1313745151846729587935189L, -0.05357102822023923595359996L,0.001770396479992522798495351L, 0.0L};
for (size_t i = 0; i < lin_4.size(); ++i)
{
if(!CHECK_ULP_CLOSE(lin_4[i], phi1_4[i], 0))
{
std::cerr << " Index " << i << " is incorrect\n";
}
}
std::array<long double, 10> lin_5 = {0.0L,1.558326313047001366564379L,-2.436012783189551921436896L,1.235905129801454293947039L,-0.3674377136938866359947561L,-0.02178035117564654658884556L,0.03234719350814368885815854L,-0.001335619912770701035229331L,-0.00001216838474354431384970525L,0.0L};
auto phi1_5 = boost::math::detail::daubechies_scaling_integer_grid<long double, 5, 1>();
for (size_t i = 0; i < lin_5.size(); ++i)
{
if(!CHECK_ULP_CLOSE(lin_5[i], phi1_5[i], 0))
{
std::cerr << " Index " << i << " is incorrect\n";
}
}
}
int main()
{
test_agreement_with_ten_lectures();
test_first_derivative();
test_dyadic_grid<float>();
test_dyadic_grid<double>();
test_dyadic_grid<long double>();
#ifdef BOOST_HAS_FLOAT128
test_dyadic_grid<float128>();
#endif
/*
test_interpolation<float>();
test_interpolation<double>();
test_interpolation<long double>();
#if BOOST_HAS_FLOAT128
test_interpolation<float128>();
#endif*/
// All scaling functions have a first derivative.
boost::hana::for_each(std::make_index_sequence<13>(), [&](auto idx){
test_integer_grid<float, idx+2, 0>();
test_integer_grid<float, idx+2, 1>();
test_integer_grid<double, idx+2, 0>();
test_integer_grid<double, idx+2, 1>();
test_integer_grid<long double, idx+2, 0>();
test_integer_grid<long double, idx+2, 1>();
#ifdef BOOST_HAS_FLOAT128
test_integer_grid<float128, idx+2, 0>();
test_integer_grid<float128, idx+2, 1>();
#endif
});
// 4-tap (2 vanishing moment) scaling function does not have a second derivative;
// all other scaling functions do.
boost::hana::for_each(std::make_index_sequence<13>(), [&](auto idx){
test_integer_grid<float, idx+3, 2>();
test_integer_grid<double, idx+3, 2>();
test_integer_grid<long double, idx+3, 2>();
#ifdef BOOST_HAS_FLOAT128
test_integer_grid<boost::multiprecision::float128, idx+3, 2>();
#endif
});
// 8-tap filter (4 vanishing moments) is the first to have a third derivative.
boost::hana::for_each(std::make_index_sequence<12>(), [&](auto idx){
test_integer_grid<float, idx+4, 3>();
test_integer_grid<double, idx+4, 3>();
test_integer_grid<long double, idx+4, 3>();
#ifdef BOOST_HAS_FLOAT128
test_integer_grid<boost::multiprecision::float128, idx+4, 3>();
#endif
});
// 10-tap filter (5 vanishing moments) is the first to have a fourth derivative.
boost::hana::for_each(std::make_index_sequence<11>(), [&](auto idx){
test_integer_grid<float, idx+5, 4>();
test_integer_grid<double, idx+5, 4>();
test_integer_grid<long double, idx+5, 4>();
#ifdef BOOST_HAS_FLOAT128
test_integer_grid<boost::multiprecision::float128, idx+5, 4>();
#endif
});
boost::hana::for_each(std::make_index_sequence<8>(), [&](auto i){
test_daubechies_filters<float, i+1>();
});
boost::hana::for_each(std::make_index_sequence<12>(), [&](auto i){
test_daubechies_filters<double, i+1>();
});
boost::hana::for_each(std::make_index_sequence<11>(), [&](auto i){
test_daubechies_filters<long double, i+1>();
});
#ifdef BOOST_HAS_FLOAT128
boost::hana::for_each(std::make_index_sequence<23>(), [&](auto i){
test_filter_ulp_distance<float128, long double, i+1>();
test_filter_ulp_distance<float128, double, i+1>();
test_filter_ulp_distance<float128, float, i+1>();
});
boost::hana::for_each(std::make_index_sequence<12>(), [&](auto i){
test_daubechies_filters<float128, i+1>();
});
#endif
return boost::math::test::report_errors();
}
Reference in New Issue View Git Blame Copy Permalink