mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
118 lines
3.9 KiB
C++
118 lines
3.9 KiB
C++
/*
|
|
* Copyright Nick Thompson, 2020
|
|
* 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 <iomanip>
|
|
#include <iostream>
|
|
#include <random>
|
|
#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/special_functions/daubechies_wavelet.hpp>
|
|
#include <boost/math/quadrature/trapezoidal.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;
|
|
|
|
|
|
template<class Real>
|
|
void test_wavelet_dyadic_grid()
|
|
{
|
|
std::cout << "Testing wavelet dyadic grid on type " << boost::core::demangle(typeid(Real).name()) << "\n";
|
|
auto f = [&](auto i)
|
|
{
|
|
auto phijk = boost::math::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 (int64_t j = 1; j < 10; ++j)
|
|
{
|
|
phijk = boost::math::dyadic_grid<Real, i+2, 0>(j);
|
|
phik = boost::math::detail::daubechies_scaling_integer_grid<Real, i+2, 0>();
|
|
for (int64_t l = 0; l < static_cast<int64_t>(phik.size()); ++l)
|
|
{
|
|
CHECK_ULP_CLOSE(phik[l], phijk[l*(int64_t(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 (int64_t l = 0; l < static_cast<int64_t>(phik.size()); ++l)
|
|
{
|
|
int64_t idx = l*(int64_t(1)<<j) + k;
|
|
if (idx < phijk.size())
|
|
{
|
|
cond += phijk[idx];
|
|
}
|
|
}
|
|
CHECK_MOLLIFIED_CLOSE(Real(1), cond.sum(), 10*cond()*std::numeric_limits<Real>::epsilon());
|
|
}
|
|
}
|
|
};
|
|
|
|
boost::hana::for_each(std::make_index_sequence<18>(), f);
|
|
}
|
|
|
|
|
|
template<typename Real, int p>
|
|
void test_quadratures()
|
|
{
|
|
using boost::math::quadrature::trapezoidal;
|
|
auto psi = boost::math::daubechies_wavelet<Real, p>();
|
|
|
|
Real tol = std::numeric_limits<Real>::epsilon();
|
|
Real error_estimate = std::numeric_limits<Real>::quiet_NaN();
|
|
Real L1 = std::numeric_limits<Real>::quiet_NaN();
|
|
auto [a, b] = psi.support();
|
|
CHECK_ULP_CLOSE(Real(-p+1), a, 0);
|
|
CHECK_ULP_CLOSE(Real(p), b, 0);
|
|
// A wavelet is a function of zero average; ensure the quadrature over its support is zero.
|
|
/*Real Q = trapezoidal(psi, a, b, tol, 15, &error_estimate, &L1);
|
|
if (!CHECK_MOLLIFIED_CLOSE(Real(0), Q, Real(0.0001)))
|
|
{
|
|
std::cerr << " Quadrature of " << p << " vanishing moment wavelet does not vanish.\n";
|
|
}*/
|
|
// psi is orthogonal to its integer translates: \int \psi(x-k) \psi(x) \, \mathrm{d}x = 0
|
|
// psi has L2 norm 1:
|
|
|
|
// g_n = 1/sqrt(2) <psi(t/2), phi(t-n)> (Mallat, 7.55)
|
|
}
|
|
|
|
int main()
|
|
{
|
|
boost::hana::for_each(std::make_index_sequence<18>(), [&](auto i){
|
|
test_quadratures<float, i+2>();
|
|
test_quadratures<double, i+2>();
|
|
});
|
|
|
|
|
|
test_wavelet_dyadic_grid<float>();
|
|
test_wavelet_dyadic_grid<double>();
|
|
test_wavelet_dyadic_grid<long double>();
|
|
#ifdef BOOST_HAS_FLOAT128
|
|
test_wavelet_dyadic_grid<float128>();
|
|
#endif
|
|
|
|
return boost::math::test::report_errors();
|
|
}
|