Skip to content

Fix Daubechies derivative scaling and complete coverage #1455

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
NAThompson merged 1 commit into from
Sep 5, 2026
Merged

Fix Daubechies derivative scaling and complete coverage #1455

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 8 additions & 4 deletions example/daubechies_wavelets/find_best_daubechies_interpolator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)

#include <array>
#include <iostream>
#include <unordered_map>
#include <string>
Expand Down Expand Up @@ -153,11 +154,14 @@ void find_best_interpolator()
}

{
auto phi_copy = phi;
auto phi_prime_copy = phi_prime;
auto mh = boost::math::detail::matched_holder(std::move(phi_copy), std::move(phi_prime_copy), r, Real(0));
std::vector<std::array<Real, 2>> data(phi.size());
for (size_t i = 0; i < phi.size(); ++i)
{
data[i] = {phi[i], phi_prime[i]};
}
auto mh = boost::math::detail::matched_holder_aos(std::move(data), r, Real(0));
Real sup = 0;
// call to matched_holder is unchecked, so only go to phi_dense.size() -1.
// call to matched_holder_aos is unchecked, so only go to phi_dense.size() -1.
for (size_t i = 0; i < phi_dense.size() - 1; ++i)
{
Real x = i*dx_dense;
Expand Down
4 changes: 2 additions & 2 deletions include/boost/math/interpolators/detail/quintic_hermite_detail.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -378,7 +378,7 @@ class cardinal_quintic_hermite_detail {
Real d2y1 = d2y_[i+1];

Real d2ydx2 = 60*t*(1 - 3*t + 2*t*t)*(y1 - y0)*inv_dx_*inv_dx_;
d2ydx2 += (12*t)*((-3 + 8*t - 5*t*t)*dy0 - (2 - 7*t + 5*t*t)*dy1);
d2ydx2 += (12*t)*((-3 + 8*t - 5*t*t)*dy0 - (2 - 7*t + 5*t*t)*dy1)*inv_dx_*inv_dx_;
d2ydx2 += (1 - 9*t + 18*t*t - 10*t*t*t)*d2y0*(2*inv_dx_*inv_dx_) + t*(3 - 12*t + 10*t*t)*d2y1*(2*inv_dx_*inv_dx_);
return d2ydx2;
}
Expand Down Expand Up @@ -559,7 +559,7 @@ class cardinal_quintic_hermite_detail_aos {
Real d2y1 = data_[i+1][2];

Real d2ydx2 = 60*t*(1 - 3*t + 2*t*t)*(y1 - y0)*inv_dx_*inv_dx_;
d2ydx2 += (12*t)*((-3 + 8*t - 5*t*t)*dy0 - (2 - 7*t + 5*t*t)*dy1);
d2ydx2 += (12*t)*((-3 + 8*t - 5*t*t)*dy0 - (2 - 7*t + 5*t*t)*dy1)*inv_dx_*inv_dx_;
d2ydx2 += (1 - 9*t + 18*t*t - 10*t*t*t)*d2y0*(2*inv_dx_*inv_dx_) + t*(3 - 12*t + 10*t*t)*d2y1*(2*inv_dx_*inv_dx_);
return d2ydx2;
}
Expand Down
12 changes: 6 additions & 6 deletions include/boost/math/interpolators/detail/septic_hermite_detail.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -322,7 +322,7 @@ class cardinal_septic_hermite_detail {
}
if (x == xf)
{
return dy_.back()/inv_dx_;
return dy_.back()*inv_dx_;
}

return this->unchecked_prime(x);
Expand All @@ -337,7 +337,7 @@ class cardinal_septic_hermite_detail {
Real t = s3 - ii;
if (t==0)
{
return dy_[i]/inv_dx_;
return dy_[i]*inv_dx_;
}

Real y0 = y_[i];
Expand All @@ -361,7 +361,7 @@ class cardinal_septic_hermite_detail {
Real dydx = z0*(y1-y0)*inv_dx_;
dydx += (z1*dy0 + z2*dy1)*inv_dx_;
dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
dydx += t*t*(z5*j0 + z6*j1);
dydx += t*t*(z5*j0 + z6*j1)*inv_dx_;
return dydx;
}

Expand Down Expand Up @@ -417,7 +417,7 @@ class cardinal_septic_hermite_detail {
Real d2ydx2 = z0*(y1-y0)*inv_dx_*inv_dx_;
d2ydx2 += (z1*dy0 + z2*dy1)*inv_dx_*inv_dx_;
d2ydx2 += (z3*a0 + z4*a1)*2*inv_dx_*inv_dx_;
d2ydx2 += 6*(z5*j0 + z6*j1)/(inv_dx_*inv_dx_);
d2ydx2 += 6*(z5*j0 + z6*j1)*inv_dx_*inv_dx_;

return d2ydx2;
}
Expand Down Expand Up @@ -572,7 +572,7 @@ class cardinal_septic_hermite_detail_aos {
Real dydx = z0*(y1-y0)*inv_dx_;
dydx += (z1*dy0 + z2*dy1)*inv_dx_;
dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
dydx += t*t*(z5*j0 + z6*j1);
dydx += t*t*(z5*j0 + z6*j1)*inv_dx_;
return dydx;
}

Expand Down Expand Up @@ -627,7 +627,7 @@ class cardinal_septic_hermite_detail_aos {
Real d2ydx2 = z0*(y1-y0)*inv_dx_*inv_dx_;
d2ydx2 += (z1*dy0 + z2*dy1)*inv_dx_*inv_dx_;
d2ydx2 += (z3*a0 + z4*a1)*2*inv_dx_*inv_dx_;
d2ydx2 += 6*(z5*j0 + z6*j1)/(inv_dx_*inv_dx_);
d2ydx2 += 6*(z5*j0 + z6*j1)*inv_dx_*inv_dx_;

return d2ydx2;
}
Expand Down
88 changes: 0 additions & 88 deletions include/boost/math/special_functions/daubechies_scaling.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,50 +102,6 @@ std::vector<Real> daubechies_scaling_dyadic_grid(int64_t j_max)

namespace detail {

template<class RandomAccessContainer>
class matched_holder {
public:
using Real = typename RandomAccessContainer::value_type;

matched_holder(RandomAccessContainer && y, RandomAccessContainer && dydx, int grid_refinements, Real x0) : x0_{x0}, y_{std::move(y)}, dy_{std::move(dydx)}
{
inv_h_ = (1 << grid_refinements);
Real h = 1/inv_h_;
for (auto & dy : dy_)
{
dy *= h;
}
}

inline Real operator()(Real x) const
{
using std::floor;
using std::sqrt;
// This is the exact Holder exponent, but it's pessimistic almost everywhere!
// It's only exactly right at dyadic rationals.
//Real const alpha = 2 - log(1+sqrt(Real(3)))/log(Real(2));
// We're gonna use alpha = 1/2, rather than 0.5500...
Real s = (x-x0_)*inv_h_;
Real ii = floor(s);
auto i = static_cast<decltype(y_.size())>(ii);
Real t = s - ii;
Real dphi = dy_[i+1];
Real diff = y_[i+1] - y_[i];
return y_[i] + (2*dphi - diff)*t + 2*sqrt(t)*(diff-dphi);
}

int64_t bytes() const
{
return 2*y_.size()*sizeof(Real) + sizeof(*this);
}

private:
Real x0_;
Real inv_h_;
RandomAccessContainer y_;
RandomAccessContainer dy_;
};

template<class RandomAccessContainer>
class matched_holder_aos {
public:
Expand Down Expand Up @@ -189,50 +145,6 @@ class matched_holder_aos {
};


template<class RandomAccessContainer>
class linear_interpolation {
public:
using Real = typename RandomAccessContainer::value_type;

linear_interpolation(RandomAccessContainer && y, RandomAccessContainer && dydx, int grid_refinements) : y_{std::move(y)}, dydx_{std::move(dydx)}
{
s_ = (1 << grid_refinements);
}

inline Real operator()(Real x) const
{
using std::floor;
Real y = x*s_;
Real k = floor(y);

int64_t kk = static_cast<int64_t>(k);
Real t = y - k;
return (1-t)*y_[kk] + t*y_[kk+1];
}

inline Real prime(Real x) const
{
using std::floor;

Real y = x*s_;
Real k = floor(y);

int64_t kk = static_cast<int64_t>(k);
Real t = y - k;
return static_cast<Real>((Real(1)-t)*dydx_[kk] + t*dydx_[kk+1]);
}

int64_t bytes() const
{
return (1 + y_.size() + dydx_.size())*sizeof(Real) + sizeof(y_) + sizeof(dydx_);
}

private:
Real s_;
RandomAccessContainer y_;
RandomAccessContainer dydx_;
};

template<class RandomAccessContainer>
class linear_interpolation_aos {
public:
Expand Down
65 changes: 55 additions & 10 deletions test/daubechies_scaling_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -459,17 +459,55 @@ void test_quadratures()
}
}

void test_absolute_error_refinement_derivatives()
// Exercise both derivative-capable interpolators at points between grid nodes.
// Differentiating prime() checks double_prime() independently of the grid generator.
// Use coarse grids so the finite-difference stencil fits inside an interpolation interval.
template<class Real, int p>
void test_refinement_derivatives(int refinements)
{
using boost::math::differentiation::finite_difference_derivative;
auto phi = boost::math::daubechies_scaling<float, 19>(-2);
auto value = [&](float x) { return phi(x); };
auto prime = [&](float x) { return phi.prime(x); };
float x = 9.03125f;
float finite_difference_prime = finite_difference_derivative<decltype(value), float, 2>(value, x);
float finite_difference_double_prime = finite_difference_derivative<decltype(prime), float, 2>(prime, x);
CHECK_MOLLIFIED_CLOSE(finite_difference_prime, phi.prime(x), 0.005f);
CHECK_MOLLIFIED_CLOSE(finite_difference_double_prime, phi.double_prime(x), 0.005f);
auto f = boost::math::daubechies_scaling<Real, p>(refinements);
auto value = [&](Real x) { return f(x); };
auto prime = [&](Real x) { return f.prime(x); };
const Real tolerance = std::is_same_v<Real, float> ? Real(0.005) : Real(0.000001);
for (int i = 0; i < 8; ++i)
{
Real x = Real(p/3) + Real(i)/8 + Real(1)/32;
auto d1 = finite_difference_derivative<decltype(value), Real, 2>(value, x);
auto d2 = finite_difference_derivative<decltype(prime), Real, 2>(prime, x);
if constexpr (std::is_same_v<Real, float>)
{
// The automatic float step is too large for the curvature of these
// interpolants. This dyadic step stays well inside the grid cell.
Real h = Real(1)/1024;
d1 = (value(x + h) - value(x - h))/(2*h);
d2 = (prime(x + h) - prime(x - h))/(2*h);
}
CHECK_MOLLIFIED_CLOSE(d1, f.prime(x), tolerance);
CHECK_MOLLIFIED_CLOSE(d2, f.double_prime(x), tolerance);
}
}

// At shared dyadic nodes, changing the refinement must preserve function values.
template<class Real, int p>
void test_explicit_refinement()
{
auto coarse = boost::math::daubechies_scaling<Real, p>(3);
auto fine = boost::math::daubechies_scaling<Real, p>(4);
CHECK_EQUAL(true, fine.bytes() > coarse.bytes());
auto [a, b] = coarse.support();
for (Real x = a; x <= b; x += Real(1)/8)
{
CHECK_MOLLIFIED_CLOSE(coarse(x), fine(x), 32*std::numeric_limits<Real>::epsilon());
if constexpr (p > 2)
{
CHECK_MOLLIFIED_CLOSE(coarse.prime(x), fine.prime(x), 32*std::numeric_limits<Real>::epsilon());
}
if constexpr (p >= 6)
{
CHECK_MOLLIFIED_CLOSE(coarse.double_prime(x), fine.double_prime(x), 32*std::numeric_limits<Real>::epsilon());
}
}
}

int main()
Expand All @@ -479,7 +517,14 @@ int main()
test_quadratures<float, i+2>();
test_quadratures<double, i+2>();
});
test_absolute_error_refinement_derivatives();
test_refinement_derivatives<float, 6>(4);
test_refinement_derivatives<float, 19>(-2);
test_refinement_derivatives<double, 6>(4);
test_refinement_derivatives<double, 19>(4);
boost::hana::for_each(std::make_index_sequence<18>(), [&](auto i) {
test_explicit_refinement<float, i + 2>();
test_explicit_refinement<double, i + 2>();
});

test_agreement_with_ten_lectures();

Expand Down
87 changes: 77 additions & 10 deletions test/daubechies_wavelet_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,31 +119,98 @@ void test_quadratures()
}
}

void test_absolute_error_refinement_derivatives()
// Exercise both derivative-capable interpolators at points between grid nodes.
// Differentiating prime() checks double_prime() independently of the grid generator.
// Use coarse grids so the finite-difference stencil fits inside an interpolation interval.
template<class Real, int p>
void test_refinement_derivatives(int refinements)
{
using boost::math::differentiation::finite_difference_derivative;
auto psi = boost::math::daubechies_wavelet<float, 19>(-2);
auto value = [&](float x) { return psi(x); };
auto prime = [&](float x) { return psi.prime(x); };
float x = 0.1875f;
float finite_difference_prime = finite_difference_derivative<decltype(value), float, 2>(value, x);
float finite_difference_double_prime = finite_difference_derivative<decltype(prime), float, 2>(prime, x);
CHECK_MOLLIFIED_CLOSE(finite_difference_prime, psi.prime(x), 0.005f);
CHECK_MOLLIFIED_CLOSE(finite_difference_double_prime, psi.double_prime(x), 0.005f);
auto f = boost::math::daubechies_wavelet<Real, p>(refinements);
auto value = [&](Real x) { return f(x); };
auto prime = [&](Real x) { return f.prime(x); };
const Real tolerance = std::is_same_v<Real, float> ? Real(0.005) : Real(0.000001);
for (int i = 0; i < 8; ++i)
{
Real x = Real(0) + Real(i)/8 + Real(1)/32;
auto d1 = finite_difference_derivative<decltype(value), Real, 2>(value, x);
auto d2 = finite_difference_derivative<decltype(prime), Real, 2>(prime, x);
if constexpr (std::is_same_v<Real, float>)
{
// The automatic float step is too large for the curvature of these
// interpolants. This dyadic step stays well inside the grid cell.
Real h = Real(1)/1024;
d1 = (value(x + h) - value(x - h))/(2*h);
d2 = (prime(x + h) - prime(x - h))/(2*h);
}
CHECK_MOLLIFIED_CLOSE(d1, f.prime(x), tolerance);
CHECK_MOLLIFIED_CLOSE(d2, f.double_prime(x), tolerance);
}
}

// At shared dyadic nodes, changing the refinement must preserve function values.
template<class Real, int p>
void test_explicit_refinement()
{
auto coarse = boost::math::daubechies_wavelet<Real, p>(3);
auto fine = boost::math::daubechies_wavelet<Real, p>(4);
CHECK_EQUAL(true, fine.bytes() > coarse.bytes());
auto [a, b] = coarse.support();
for (Real x = a; x <= b; x += Real(1)/8)
{
CHECK_MOLLIFIED_CLOSE(coarse(x), fine(x), 32*std::numeric_limits<Real>::epsilon());
if constexpr (p > 2)
{
CHECK_MOLLIFIED_CLOSE(coarse.prime(x), fine.prime(x), 32*std::numeric_limits<Real>::epsilon());
}
if constexpr (p >= 6)
{
CHECK_MOLLIFIED_CLOSE(coarse.double_prime(x), fine.double_prime(x), 32*std::numeric_limits<Real>::epsilon());
}
}
}

template<class Real>
void test_haar_wavelet()
{
for (int refinements : {-2, -1, 3})
{
auto psi = boost::math::daubechies_wavelet<Real, 1>(refinements);
CHECK_EQUAL(Real(0), psi.support().first);
CHECK_EQUAL(Real(1), psi.support().second);
for (Real x : {Real(-1), Real(0), Real(0.5), Real(1), Real(2)})
{
CHECK_EQUAL(Real(0), psi(x));
}
CHECK_EQUAL(Real(1), psi(Real(0.25)));
CHECK_EQUAL(Real(-1), psi(Real(0.75)));
CHECK_EQUAL(Real(1), psi(std::nextafter(Real(0.5), Real(0))));
CHECK_EQUAL(Real(-1), psi(std::nextafter(Real(0.5), Real(1))));
}
CHECK_THROW((boost::math::daubechies_wavelet<Real, 1>(0)), std::domain_error);
}

int main()
{
#ifndef __MINGW32__
try
{
test_haar_wavelet<float>();
test_haar_wavelet<double>();
test_exact_value<double>();

boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i) {
test_quadratures<float, i + 3>();
test_quadratures<double, i + 3>();
});
test_absolute_error_refinement_derivatives();
test_refinement_derivatives<float, 6>(4);
test_refinement_derivatives<float, 19>(-2);
test_refinement_derivatives<double, 6>(4);
test_refinement_derivatives<double, 19>(4);
boost::hana::for_each(std::make_index_sequence<18>(), [&](auto i) {
test_explicit_refinement<float, i + 2>();
test_explicit_refinement<double, i + 2>();
});
}
catch (std::bad_alloc)
{
Expand Down
Loading
Loading