From 9a4365882ce430ae1a4a0ad43542b1fb8870ffb1 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 5 Sep 2026 10:41:28 -0700 Subject: [PATCH] Fix Daubechies derivative scaling and complete coverage Exercise explicit refinements, absolute-error mode, interior derivatives, and the Haar wavelet. Remove obsolete interpolators and migrate the remaining example to matched_holder_aos. Correct missing grid-spacing factors in cardinal quintic and septic derivatives, including septic node evaluation. Add cubic-polynomial regressions for both data layouts and non-unit grid spacing. Fixes boostorg/math#1101 --- .../find_best_daubechies_interpolator.cpp | 12 ++- .../detail/quintic_hermite_detail.hpp | 4 +- .../detail/septic_hermite_detail.hpp | 12 +-- .../special_functions/daubechies_scaling.hpp | 88 ------------------- test/daubechies_scaling_test.cpp | 65 +++++++++++--- test/daubechies_wavelet_test.cpp | 87 +++++++++++++++--- test/quintic_hermite_test.cpp | 31 +++++++ test/septic_hermite_test.cpp | 32 +++++++ 8 files changed, 211 insertions(+), 120 deletions(-) diff --git a/example/daubechies_wavelets/find_best_daubechies_interpolator.cpp b/example/daubechies_wavelets/find_best_daubechies_interpolator.cpp index 56c5ed3cba..381ec85218 100644 --- a/example/daubechies_wavelets/find_best_daubechies_interpolator.cpp +++ b/example/daubechies_wavelets/find_best_daubechies_interpolator.cpp @@ -4,6 +4,7 @@ // (See accompanying file LICENSE_1_0.txt // or copy at http://www.boost.org/LICENSE_1_0.txt) +#include #include #include #include @@ -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> 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; diff --git a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp index 85b1458fcd..baf55843af 100644 --- a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp +++ b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp @@ -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; } @@ -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; } diff --git a/include/boost/math/interpolators/detail/septic_hermite_detail.hpp b/include/boost/math/interpolators/detail/septic_hermite_detail.hpp index e4bcdf1f1f..299ea6f82f 100644 --- a/include/boost/math/interpolators/detail/septic_hermite_detail.hpp +++ b/include/boost/math/interpolators/detail/septic_hermite_detail.hpp @@ -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); @@ -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]; @@ -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; } @@ -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; } @@ -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; } @@ -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; } diff --git a/include/boost/math/special_functions/daubechies_scaling.hpp b/include/boost/math/special_functions/daubechies_scaling.hpp index db07f47e0d..81d1975675 100644 --- a/include/boost/math/special_functions/daubechies_scaling.hpp +++ b/include/boost/math/special_functions/daubechies_scaling.hpp @@ -102,50 +102,6 @@ std::vector daubechies_scaling_dyadic_grid(int64_t j_max) namespace detail { -template -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(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 matched_holder_aos { public: @@ -189,50 +145,6 @@ class matched_holder_aos { }; -template -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(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(k); - Real t = y - k; - return static_cast((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 linear_interpolation_aos { public: diff --git a/test/daubechies_scaling_test.cpp b/test/daubechies_scaling_test.cpp index fba8c3c8eb..a587e4bb31 100644 --- a/test/daubechies_scaling_test.cpp +++ b/test/daubechies_scaling_test.cpp @@ -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 +void test_refinement_derivatives(int refinements) { using boost::math::differentiation::finite_difference_derivative; - auto phi = boost::math::daubechies_scaling(-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(value, x); - float finite_difference_double_prime = finite_difference_derivative(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(refinements); + auto value = [&](Real x) { return f(x); }; + auto prime = [&](Real x) { return f.prime(x); }; + const Real tolerance = std::is_same_v ? 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(value, x); + auto d2 = finite_difference_derivative(prime, x); + if constexpr (std::is_same_v) + { + // 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 +void test_explicit_refinement() +{ + auto coarse = boost::math::daubechies_scaling(3); + auto fine = boost::math::daubechies_scaling(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::epsilon()); + if constexpr (p > 2) + { + CHECK_MOLLIFIED_CLOSE(coarse.prime(x), fine.prime(x), 32*std::numeric_limits::epsilon()); + } + if constexpr (p >= 6) + { + CHECK_MOLLIFIED_CLOSE(coarse.double_prime(x), fine.double_prime(x), 32*std::numeric_limits::epsilon()); + } + } } int main() @@ -479,7 +517,14 @@ int main() test_quadratures(); test_quadratures(); }); - test_absolute_error_refinement_derivatives(); + test_refinement_derivatives(4); + test_refinement_derivatives(-2); + test_refinement_derivatives(4); + test_refinement_derivatives(4); + boost::hana::for_each(std::make_index_sequence<18>(), [&](auto i) { + test_explicit_refinement(); + test_explicit_refinement(); + }); test_agreement_with_ten_lectures(); diff --git a/test/daubechies_wavelet_test.cpp b/test/daubechies_wavelet_test.cpp index 5afaf4d825..7fbc7f0f97 100644 --- a/test/daubechies_wavelet_test.cpp +++ b/test/daubechies_wavelet_test.cpp @@ -119,17 +119,75 @@ 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 +void test_refinement_derivatives(int refinements) { using boost::math::differentiation::finite_difference_derivative; - auto psi = boost::math::daubechies_wavelet(-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(value, x); - float finite_difference_double_prime = finite_difference_derivative(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(refinements); + auto value = [&](Real x) { return f(x); }; + auto prime = [&](Real x) { return f.prime(x); }; + const Real tolerance = std::is_same_v ? 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(value, x); + auto d2 = finite_difference_derivative(prime, x); + if constexpr (std::is_same_v) + { + // 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 +void test_explicit_refinement() +{ + auto coarse = boost::math::daubechies_wavelet(3); + auto fine = boost::math::daubechies_wavelet(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::epsilon()); + if constexpr (p > 2) + { + CHECK_MOLLIFIED_CLOSE(coarse.prime(x), fine.prime(x), 32*std::numeric_limits::epsilon()); + } + if constexpr (p >= 6) + { + CHECK_MOLLIFIED_CLOSE(coarse.double_prime(x), fine.double_prime(x), 32*std::numeric_limits::epsilon()); + } + } +} + +template +void test_haar_wavelet() +{ + for (int refinements : {-2, -1, 3}) + { + auto psi = boost::math::daubechies_wavelet(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(0)), std::domain_error); } int main() @@ -137,13 +195,22 @@ int main() #ifndef __MINGW32__ try { + test_haar_wavelet(); + test_haar_wavelet(); test_exact_value(); boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i) { test_quadratures(); test_quadratures(); }); - test_absolute_error_refinement_derivatives(); + test_refinement_derivatives(4); + test_refinement_derivatives(-2); + test_refinement_derivatives(4); + test_refinement_derivatives(4); + boost::hana::for_each(std::make_index_sequence<18>(), [&](auto i) { + test_explicit_refinement(); + test_explicit_refinement(); + }); } catch (std::bad_alloc) { diff --git a/test/quintic_hermite_test.cpp b/test/quintic_hermite_test.cpp index d2deec7709..c5007ff3c9 100644 --- a/test/quintic_hermite_test.cpp +++ b/test/quintic_hermite_test.cpp @@ -524,8 +524,39 @@ void test_cardinal_quartic() } +// Non-unit spacing exposes missing chain-rule factors in derivative evaluation. +void test_cardinal_derivative_scaling() +{ + for (double dx : {0.125, 2.0}) + { + std::vector> data(4); + std::vector y(4), dy(4), d2y(4); + for (size_t i = 0; i < data.size(); ++i) + { + double x = -1 + i*dx; + data[i] = {x*x*x, 3*x*x, 6*x}; + y[i] = x*x*x; + dy[i] = 3*x*x; + d2y[i] = 6*x; + } + auto f = boost::math::interpolators::detail::cardinal_quintic_hermite_detail_aos(std::move(data), -1.0, dx); + auto g = boost::math::interpolators::detail::cardinal_quintic_hermite_detail(std::move(y), std::move(dy), std::move(d2y), -1.0, dx); + for (int i = 0; i <= 24; ++i) + { + double x = -1 + i*dx/8; + CHECK_MOLLIFIED_CLOSE(x*x*x, f(x), 1e-12); + CHECK_MOLLIFIED_CLOSE(3*x*x, f.prime(x), 1e-12); + CHECK_MOLLIFIED_CLOSE(6*x, f.double_prime(x), 1e-12); + CHECK_MOLLIFIED_CLOSE(x*x*x, g(x), 1e-12); + CHECK_MOLLIFIED_CLOSE(3*x*x, g.prime(x), 1e-12); + CHECK_MOLLIFIED_CLOSE(6*x, g.double_prime(x), 1e-12); + } + } +} + int main() { + test_cardinal_derivative_scaling(); #ifdef __STDCPP_FLOAT32_T__ test_constant(); test_linear(); diff --git a/test/septic_hermite_test.cpp b/test/septic_hermite_test.cpp index 119e626191..7b3920f6e3 100644 --- a/test/septic_hermite_test.cpp +++ b/test/septic_hermite_test.cpp @@ -497,8 +497,40 @@ void test_interpolation_condition() } +// Non-unit spacing exposes missing chain-rule factors in derivative evaluation. +void test_cardinal_derivative_scaling() +{ + for (double dx : {0.125, 2.0}) + { + std::vector> data(4); + std::vector y(4), dy(4), d2y(4), d3y(4); + for (size_t i = 0; i < data.size(); ++i) + { + double x = -1 + i*dx; + data[i] = {x*x*x, 3*x*x, 6*x, 6}; + y[i] = x*x*x; + dy[i] = 3*x*x; + d2y[i] = 6*x; + d3y[i] = 6; + } + auto f = boost::math::interpolators::detail::cardinal_septic_hermite_detail_aos(std::move(data), -1.0, dx); + auto g = boost::math::interpolators::detail::cardinal_septic_hermite_detail(std::move(y), std::move(dy), std::move(d2y), std::move(d3y), -1.0, dx); + for (int i = 0; i <= 24; ++i) + { + double x = -1 + i*dx/8; + CHECK_MOLLIFIED_CLOSE(x*x*x, f(x), 1e-12); + CHECK_MOLLIFIED_CLOSE(3*x*x, f.prime(x), 1e-12); + CHECK_MOLLIFIED_CLOSE(6*x, f.double_prime(x), 1e-12); + CHECK_MOLLIFIED_CLOSE(x*x*x, g(x), 1e-12); + CHECK_MOLLIFIED_CLOSE(3*x*x, g.prime(x), 1e-12); + CHECK_MOLLIFIED_CLOSE(6*x, g.double_prime(x), 1e-12); + } + } +} + int main() { + test_cardinal_derivative_scaling(); #ifdef __STDCPP_FLOAT32_T__ test_constant(); test_linear();