From b2f3054e2f7b52a42873c53697e1ebd24ef779c7 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Tue, 22 Jan 2019 15:46:19 -0700 Subject: [PATCH] Compute filters in higher precision and cast back to low precision so that high-p filters are accurate. [CI SKIP] --- doc/differentiation/lanczos_smoothing.qbk | 4 +- .../differentiation/lanczos_smoothing.hpp | 56 +++++-- test/lanczos_smoothing_test.cpp | 144 ++++-------------- 3 files changed, 74 insertions(+), 130 deletions(-) diff --git a/doc/differentiation/lanczos_smoothing.qbk b/doc/differentiation/lanczos_smoothing.qbk index fb3ba10eb..b9ec5a011 100644 --- a/doc/differentiation/lanczos_smoothing.qbk +++ b/doc/differentiation/lanczos_smoothing.qbk @@ -27,7 +27,9 @@ namespace boost::math::differentiation { template RandomAccessContainer operator()(RandomAccessContainer const & v) const; - void reset_spacing(Real spacing); + void set_spacing(Real spacing); + + Real get_spacing() const; }; } // namespaces diff --git a/include/boost/math/differentiation/lanczos_smoothing.hpp b/include/boost/math/differentiation/lanczos_smoothing.hpp index 0d39ae1b0..05eaf21ec 100644 --- a/include/boost/math/differentiation/lanczos_smoothing.hpp +++ b/include/boost/math/differentiation/lanczos_smoothing.hpp @@ -1,4 +1,4 @@ -// (C) Copyright Nick Thompson 2018. +// (C) 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) @@ -7,7 +7,6 @@ #define BOOST_MATH_DIFFERENTIATION_LANCZOS_SMOOTHING_HPP #include #include -#include namespace boost::math::differentiation { @@ -276,23 +275,53 @@ public: "The approximation order must be <= 2n"); BOOST_ASSERT_MSG(approximation_order >= 2, "The approximation order must be >= 2"); - m_f = detail::interior_filter(n, approximation_order); + + if constexpr (std::is_same_v || std::is_same_v) + { + auto interior = detail::interior_filter(n, approximation_order); + m_f.resize(interior.size()); + for (size_t j = 0; j < interior.size(); ++j) + { + m_f[j] = static_cast(interior[j]); + } + } + else + { + m_f = detail::interior_filter(n, approximation_order); + } m_boundary_filters.resize(n); + // This for loop is a natural candidate for parallelization. + // But does it matter? Probably not. for (size_t i = 0; i < n; ++i) { - // s = i - n; - int64_t s = static_cast(i) - static_cast(n); - m_boundary_filters[i] = detail::boundary_filter(n, approximation_order, s); + if constexpr (std::is_same_v || std::is_same_v) + { + int64_t s = static_cast(i) - static_cast(n); + auto bf = detail::boundary_filter(n, approximation_order, s); + m_boundary_filters[i].resize(bf.size()); + for (size_t j = 0; j < bf.size(); ++j) + { + m_boundary_filters[i][j] = static_cast(bf[j]); + } + } + else + { + int64_t s = static_cast(i) - static_cast(n); + m_boundary_filters[i] = detail::boundary_filter(n, approximation_order, s); + } } } else if constexpr (order == 2) { + // High precision isn't warranted for small p; only for large p. + // (The computation appears stable for large n.) + // But given that the filters are reusable for many vectors, + // it's better to do a high precision computation and then cast back, + // since the resulting cost is a factor of 2, and the cost of the filters not working is hours of debugging. if constexpr (std::is_same_v || std::is_same_v) { - std::cout << "We're using high precision!\n"; - using boost::multiprecision::cpp_bin_float_50; - auto f = detail::acceleration_boundary_filter(n, approximation_order, 0); + auto f = detail::acceleration_boundary_filter(n, approximation_order, 0); m_f.resize(n+1); for (size_t i = 0; i < m_f.size(); ++i) { @@ -302,7 +331,7 @@ public: for (size_t i = 0; i < n; ++i) { int64_t s = static_cast(i) - static_cast(n); - auto bf = detail::acceleration_boundary_filter(n, approximation_order, s); + auto bf = detail::acceleration_boundary_filter(n, approximation_order, s); m_boundary_filters[i].resize(bf.size()); for (size_t j = 0; j < bf.size(); ++j) { @@ -312,7 +341,8 @@ public: } else { - std::cout << "No high precision.\n"; + // Given that the purpose is denoising, for higher precision calculations, + // the default precision should be fine. auto f = detail::acceleration_boundary_filter(n, approximation_order, 0); m_f.resize(n+1); for (size_t i = 0; i < m_f.size(); ++i) @@ -333,13 +363,13 @@ public: } } - void reset_spacing(Real const & spacing) + void set_spacing(Real const & spacing) { BOOST_ASSERT_MSG(spacing > 0, "Spacing between samples must be > 0."); m_dt = spacing; } - Real spacing() const + Real get_spacing() const { return m_dt; } diff --git a/test/lanczos_smoothing_test.cpp b/test/lanczos_smoothing_test.cpp index d806a11a5..cf41c8016 100644 --- a/test/lanczos_smoothing_test.cpp +++ b/test/lanczos_smoothing_test.cpp @@ -1,5 +1,5 @@ /* - * Copyright Nick Thompson, 2017 + * 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) @@ -27,54 +27,6 @@ using boost::math::differentiation::detail::discrete_legendre; using boost::math::differentiation::detail::interior_filter; using boost::math::differentiation::detail::boundary_filter; -template -size_t l1_ulp_error(RandomAccessContainer1 const & v1, RandomAccessContainer2 const & v2) -{ - using std::abs; - using Real1 = typename RandomAccessContainer1::value_type; - using Real2 = typename RandomAccessContainer2::value_type; - BOOST_ASSERT_MSG(sizeof(Real1) <= sizeof(Real2), - "Second container must contain more accurate value than the first container."); - BOOST_ASSERT_MSG(v1.size() == v2.size(), - "Containers must have same size."); - size_t error = 0; - for (size_t i = 0; i < v1.size(); ++i) - { - auto err = abs(boost::math::float_distance(v1[i], static_cast(v2[i]))); - if (abs(v2[i]) > std::numeric_limits::epsilon()) - { - error += static_cast(err); - } - } - return error; -} - -template -std::tuple sup_ulp_error(RandomAccessContainer1 const & v1, RandomAccessContainer2 const & v2) -{ - using std::abs; - using Real1 = typename RandomAccessContainer1::value_type; - using Real2 = typename RandomAccessContainer2::value_type; - BOOST_ASSERT_MSG(sizeof(Real1) <= sizeof(Real2), - "Second container must contain more accurate value than the first container."); - BOOST_ASSERT_MSG(v1.size() == v2.size(), - "Containers must have same size."); - size_t error = 0; - Real1 worst_val1 = std::numeric_limits::quiet_NaN(); - Real2 worst_val2 = std::numeric_limits::quiet_NaN(); - for (size_t i = 0; i < v1.size(); ++i) - { - auto err = abs(boost::math::float_distance(v1[i], static_cast(v2[i]))); - if (err > error) - { - error = err; - worst_val1 = v1[i]; - worst_val2 = v2[i]; - } - } - return {error, worst_val1, worst_val2}; -} - template void test_dlp_norms() { @@ -526,40 +478,18 @@ template void test_acceleration_filters() { Real eps = std::numeric_limits::epsilon(); - std::cout << std::setprecision(std::numeric_limits::digits10); - std::cout << std::scientific; - for (size_t n = 1; n < 100; ++n) - { - for(size_t p = 3; p <= 2*n; ++p) - { - for(int64_t s = -int64_t(n); s <= 0 /*int64_t(n)*/; ++s) - { - auto v1 = boost::math::differentiation::detail::acceleration_boundary_filter(n,p,s); - auto v2 = boost::math::differentiation::detail::acceleration_boundary_filter(n,p,s); - size_t dist1 = l1_ulp_error(v1, v2); - auto [dist2, worst_val1, worst_val2] = sup_ulp_error(v1, v2); - std::cout << "(n,p,s) = (" << n << ", " << p << ", " << s << ") = " - << dist1 << ", sup = " << dist2 << ", worst = " << worst_val1 << ", actual = " << worst_val2 << "\n"; - } - } - } - - /*Real eps = std::numeric_limits::epsilon(); - for (size_t n = 1; n < 15; ++n) + for (size_t n = 1; n < 5; ++n) { for(size_t p = 3; p <= 2*n; ++p) { for(int64_t s = -int64_t(n); s <= 0; ++s) { - auto f = boost::math::differentiation::detail::acceleration_boundary_filter(n,p,s); - if (std::is_same_v) + auto g = boost::math::differentiation::detail::acceleration_boundary_filter(n,p,s); + + std::vector f(g.size()); + for (size_t i = 0; i < g.size(); ++i) { - //std::cout << "We're using high precision!\n"; - auto f1 = boost::math::differentiation::detail::acceleration_boundary_filter(n,p,s); - for (size_t i = 0; i < f1.size(); ++i) - { - f[i] = static_cast(f1[i]); - } + f[i] = static_cast(g[i]); } Real sum = 0; @@ -590,6 +520,9 @@ void test_acceleration_filters() BOOST_CHECK_CLOSE_FRACTION(sum, 2, 4*sqrt(eps)); // See unlabelled equation in McDevitt, 2012, just after equation 26: // It appears that there is an off-by-one error in that equation, since p + 1 moments don't vanish, only p. + // This test is itself suspect; the condition number of the moment sum is infinite. + // So the *slightest* error in the filter gets amplified by the test; in terms of the + // behavior of the actual filter, it's not a big deal. for (size_t l = 3; l <= p; ++l) { sum = 0; @@ -603,11 +536,11 @@ void test_acceleration_filters() c = (t - sum) - y; sum = t; } - BOOST_CHECK_SMALL(abs(sum), 500*l*sqrt(eps)); + BOOST_CHECK_SMALL(abs(sum), 500*sqrt(eps)); } } } - }*/ + } } template @@ -636,56 +569,32 @@ void test_lanczos_acceleration() } for (size_t i = 0; i < v.size(); ++i) { - BOOST_CHECK_CLOSE_FRACTION(lanczos(v, i), 14, 1000*eps); + BOOST_CHECK_CLOSE_FRACTION(lanczos(v, i), 14, 1500*eps); } - v.resize(2048); - Real step = two_pi()/v.size(); - for(size_t i = 0; i < v.size(); ++i) - { - Real x = i*step; - v[i] = sin(x); - } - - std::random_device rd{}; - auto seed = rd(); - std::cout << "seed = " << seed << "\n"; + // Now add noise, and kick up the smoothing of the Lanzcos derivative (increase n): + //std::random_device rd{}; + //auto seed = rd(); + //std::cout << "seed = " << seed << "\n"; + size_t seed = 2507134629; std::mt19937 gen(seed); - std::normal_distribution<> dis{0, 0.00001}; + Real std_dev = 0.1; + std::normal_distribution dis{0, std_dev}; for (size_t i = 0; i < v.size(); ++i) { v[i] += dis(gen); } - - - size_t n = 100; - lanczos = discrete_lanczos_derivative(step, n, 3); + lanczos = discrete_lanczos_derivative(Real(1), 18, 3); auto w = lanczos(v); - BOOST_TEST(w.size() == v.size()); - BOOST_CHECK_SMALL(w[0], 0.01); - for(size_t i = 1; i < n; ++i) + for (size_t i = 0; i < v.size(); ++i) { - BOOST_CHECK_CLOSE_FRACTION(w[i], -v[i], 0.01); + BOOST_CHECK_CLOSE_FRACTION(w[i], 14, std_dev/200); } - - for(size_t i = n; i < v.size() -n; ++i) - { - BOOST_CHECK_CLOSE_FRACTION(w[i], -v[i], 0.01); - } - - /*for(size_t i = v.size() - n; i < v.size(); ++i) - { - BOOST_CHECK_CLOSE_FRACTION(w[i], -v[i], 0.01); - }*/ - } BOOST_AUTO_TEST_CASE(lanczos_smoothing_test) { - - //test_acceleration_filters(); - //test_acceleration_filters(); - /*test_dlp_second_derivative(); + test_dlp_second_derivative(); test_dlp_norms(); test_dlp_evaluation(); test_dlp_derivatives(); @@ -701,7 +610,10 @@ BOOST_AUTO_TEST_CASE(lanczos_smoothing_test) test_interior_filter(); test_interior_filter(); - test_interior_lanczos();*/ + test_interior_lanczos(); + test_acceleration_filters(); + + test_lanczos_acceleration(); test_lanczos_acceleration(); }