From 7bbf05e8ba3dafc9a433347c8fef59af4093eb37 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Mon, 31 Dec 2018 20:11:25 -0700 Subject: [PATCH] Lanczos smoothing differentiators. --- doc/differentiation/lanczos_smoothing.qbk | 110 ++++ doc/math.qbk | 1 + .../differentiation/lanczos_smoothing.hpp | 297 +++++++++++ test/Jamfile.v2 | 1 + test/lanczos_smoothing_test.cpp | 469 ++++++++++++++++++ 5 files changed, 878 insertions(+) create mode 100644 doc/differentiation/lanczos_smoothing.qbk create mode 100644 include/boost/math/differentiation/lanczos_smoothing.hpp create mode 100644 test/lanczos_smoothing_test.cpp diff --git a/doc/differentiation/lanczos_smoothing.qbk b/doc/differentiation/lanczos_smoothing.qbk new file mode 100644 index 000000000..213fc964e --- /dev/null +++ b/doc/differentiation/lanczos_smoothing.qbk @@ -0,0 +1,110 @@ +[/ +Copyright (c) 2019 Nick Thompson +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) +] + +[section:diff Lanczos Smoothing Derivatives] + +[heading Synopsis] + +`` +#include + +namespace boost { namespace math { namespace differentiation { + + template + class lanczos_derivative { + public: + lanczos_smoothing_derivative(RandomAccessContainer const & v, + typename RandomAccessContainer::value_type spacing = 1, + size_t n = 18, + size_t approximation_order = 3); + + typename RandomAccessContainer::value_type operator[](size_t i) const; + + void reset_data(RandomAccessContainer const &v); + + void reset_spacing(Real spacing); + } + +}}} // namespaces +`` + +[heading Description] + +The function `lanczos_derivative` calculates a finite-difference approximation to the derivative of a noisy sequence of equally-spaced values /v/ at an index /i/. +A basic usage is + + std::vector v(500); + // fill v with noisy data. + double spacing = 0.001; + using boost::math::differentiation::lanczos_derivative; + auto lanczos = lanczos_derivative(v, spacing); + double dvdt = lanczos[30]; + +If the data has variance \u03C3[super 2], then the variance of the computed derivative is roughly \u03C3[super 2]/p/[super 3] /n/[super -3] \u0394 /t/[super -2], +i.e., it increases cubically with the approximation order /p/, linearly with the data variance, and decreases at the cube of the filter length /n/. +In addition, we must not forget the discretization error which is /O/(\u0394 /t/[super /p/]). +You can play around with the approximation order /p/ and the filter length /n/: + + size_t n = 12; + size_t p = 2; + auto lanczos = lanczos_derivative(v, spacing, n, p); + double dvdt = lanczos[24]; + +If /p=2n/, then the Lanczos smoothing derivative is not smoothing: +It reduces to the standard /2n+1/-point finite-difference formula. +For /p>2n/, an assertion is hit as the filter is undefined. + +In our tests with AWGN, we have found the error decreases monotonically with /n/, as is expected from the theory discussed above. +So the choice of /n/ is simple: +As high as possible given your speed requirements (larger /n/ implies a longer filter and hence more compute), +balanced against the danger of overfitting and averaging over non-stationarity. + +The choice of approximation order /p/ for a given /n/ is more difficult. +If your signal is believed to be a polynomial, +it does not make sense to set /p/ to larger than the polynomial degree- +though it may be sensible to take /p/ less than the polynomial degree. + +For a sinusoidal signal contaminated with AWGN, we ran a few tests showing that for SNR = 1, p = n/8 gave the best results, +for SNR = 10, p = n/7 was the best, and for SNR = 100, p = n/6 was the most reasonable choice. +For SNR = 0.1, the method appears to be useless. +The user is urged to use these results with caution-they have no theoretical backing and are extrapolated from a single case. + +The filters are (regrettably) computed at runtime-the vast number of combinations of approximation order and filter length makes the number of filters that must be stored excessive for compile-time data. +Hence the constructor call computes the filters. +Since each filter has length /2n+1/ and there are /n/ filters, whose element each consist of /p/ summands, +the complexity of the constructor call is O(/n/[super 2]/p/).e +This is not cheap-though for most cases small /p/ and /n/ not too large (< 20) is desired. +However, for concreteness, on the author's 2.7GHz Intel laptop CPU, the /n=16/, /p=3/ filter takes 9 microseconds to compute. +This is far from negligible, and as such we provide API calls which allow the filters to be used with multiple data: + + + std::vector v(500); + // fill v with noisy data. + auto lanczos = lanczos_derivative(v, spacing); + // use lanczos with v . . . + std::vector w(500); + lanczos.reset_data(w); + // use lanczos with w . . . + // need to use a different spacing? + lanczos.reset_spacing(0.02); + + +The implementation follows [@https://doi.org/10.1080/00207160.2012.666348 McDevitt, 2012], +who vastly expanded the ideas of Lanczos to create a very general framework for numerically differentiating noisy equispaced data. + + + +[heading References] + +* Corless, Robert M., and Nicolas Fillion. ['A graduate introduction to numerical methods.] AMC 10 (2013): 12. + +* Lanczos, Cornelius. ['Applied analysis.] Courier Corporation, 1988. + +* Timothy J. McDevitt (2012): ['Discrete Lanczos derivatives of noisy data], International Journal of Computer Mathematics, 89:7, 916-931 + + +[endsect] diff --git a/doc/math.qbk b/doc/math.qbk index 711e88f8c..0090e4ea0 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -661,6 +661,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include quadrature/double_exponential.qbk] [include quadrature/naive_monte_carlo.qbk] [include differentiation/numerical_differentiation.qbk] +[include differentiation/lanczos_smoothing.qbk] [endmathpart] [include complex/complex-tr1.qbk] diff --git a/include/boost/math/differentiation/lanczos_smoothing.hpp b/include/boost/math/differentiation/lanczos_smoothing.hpp new file mode 100644 index 000000000..5696fc454 --- /dev/null +++ b/include/boost/math/differentiation/lanczos_smoothing.hpp @@ -0,0 +1,297 @@ +// (C) Copyright Nick Thompson 2018. +// 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) + +#ifndef BOOST_MATH_DIFFERENTIATION_LANCZOS_SMOOTHING_HPP +#define BOOST_MATH_DIFFERENTIATION_LANCZOS_SMOOTHING_HPP +#include +#include + +namespace boost { +namespace math { +namespace differentiation { + +namespace detail { +template +class discrete_legendre { + public: + discrete_legendre(int n) : m_n{n} + { + // The integer n indexes a family of discrete Legendre polynomials indexed by k <= 2*n + } + + Real norm_sq(int r) + { + Real prod = Real(2) / Real(2 * r + 1); + for (int k = -r; k <= r; ++k) { + prod *= Real(2 * m_n + 1 + k) / Real(2 * m_n); + } + return prod; + } + + void initialize_recursion(Real x) + { + m_qrm2 = 1; + m_qrm1 = x; + m_qrm2p = 0; + m_qrm1p = 1; + + m_r = 2; + m_x = x; + } + + Real next() + { + Real N = 2 * m_n + 1; + Real num = (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) * m_qrm2; + Real tmp = (2 * m_r - 1) * m_x * m_qrm1 - num / Real(4 * m_n * m_n); + m_qrm2 = m_qrm1; + m_qrm1 = tmp / m_r; + ++m_r; + return m_qrm1; + } + + Real next_prime() + { + Real N = 2 * m_n + 1; + Real s = + (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) / Real(4 * m_n * m_n); + Real tmp1 = ((2 * m_r - 1) * m_x * m_qrm1 - s * m_qrm2) / m_r; + Real tmp2 = ((2 * m_r - 1) * (m_qrm1 + m_x * m_qrm1p) - s * m_qrm2p) / m_r; + m_qrm2 = m_qrm1; + m_qrm1 = tmp1; + m_qrm2p = m_qrm1p; + m_qrm1p = tmp2; + ++m_r; + return m_qrm1p; + } + + + Real operator()(Real x, int k) + { + BOOST_ASSERT_MSG(k <= 2 * m_n, "r <= 2n is required."); + if (k == 0) + { + return 1; + } + if (k == 1) + { + return x; + } + Real qrm2 = 1; + Real qrm1 = x; + Real N = 2 * m_n + 1; + for (int r = 2; r <= k; ++r) { + Real num = (r - 1) * (N * N - (r - 1) * (r - 1)) * qrm2; + Real tmp = (2 * r - 1) * x * qrm1 - num / Real(4 * m_n * m_n); + qrm2 = qrm1; + qrm1 = tmp / r; + } + return qrm1; + } + + Real prime(Real x, int k) { + BOOST_ASSERT_MSG(k <= 2 * m_n, "r <= 2n is required."); + if (k == 0) { + return 0; + } + if (k == 1) { + return 1; + } + Real qrm2 = 1; + Real qrm1 = x; + Real qrm2p = 0; + Real qrm1p = 1; + Real N = 2 * m_n + 1; + for (int r = 2; r <= k; ++r) { + Real s = + (r - 1) * (N * N - (r - 1) * (r - 1)) / Real(4 * m_n * m_n); + Real tmp1 = ((2 * r - 1) * x * qrm1 - s * qrm2) / r; + Real tmp2 = ((2 * r - 1) * (qrm1 + x * qrm1p) - s * qrm2p) / r; + qrm2 = qrm1; + qrm1 = tmp1; + qrm2p = qrm1p; + qrm1p = tmp2; + } + return qrm1p; + } + + private: + int m_n; + Real m_qrm2; + Real m_qrm1; + Real m_qrm2p; + Real m_qrm1p; + int m_r; + Real m_x; +}; + +template +std::vector interior_filter(int n, int p) { + // We could make the filter length n, as f[0] = 0, + // but that'd make the indexing awkward when applying the filter. + std::vector f(n + 1, 0); + auto dlp = discrete_legendre(n); + std::vector coeffs(p+1, std::numeric_limits::quiet_NaN()); + dlp.initialize_recursion(0); + coeffs[1] = 1/dlp.norm_sq(1); + for (int l = 3; l < p + 1; l += 2) + { + dlp.next_prime(); + coeffs[l] = dlp.next_prime()/ dlp.norm_sq(l); + } + + for (size_t j = 1; j < f.size(); ++j) + { + Real arg = Real(j) / Real(n); + dlp.initialize_recursion(arg); + f[j] = coeffs[1]*arg; + for (int l = 3; l <= p; l += 2) + { + dlp.next(); + f[j] += coeffs[l]*dlp.next(); + } + f[j] /= (n * n); + } + return f; +} + +template +std::vector boundary_filter(int n, int p, int s) { + std::vector f(2 * n + 1, 0); + auto dlp = discrete_legendre(n); + Real sn = Real(s) / Real(n); + std::vector coeffs(p+1, std::numeric_limits::quiet_NaN()); + dlp.initialize_recursion(sn); + coeffs[0] = 0; + coeffs[1] = 1/dlp.norm_sq(1); + for (int l = 2; l < p + 1; ++l) + { + // Calculation of the norms is common to all filters, + // so it seems like an obvious optimization target. + // I tried this: The spent in computing the norms time is not negligible, + // but still a small fraction of the total compute time. + // Hence I'm not refactoring out these norm calculations. + coeffs[l] = dlp.next_prime()/ dlp.norm_sq(l); + } + + for (int k = 0; k < f.size(); ++k) + { + int j = k - n; + f[k] = 0; + Real arg = Real(j) / Real(n); + dlp.initialize_recursion(arg); + f[k] = coeffs[1]*arg; + for (int l = 2; l <= p; ++l) + { + f[k] += coeffs[l]*dlp.next(); + } + f[k] /= (n * n); + } + return f; +} + +} // namespace detail + +template +class lanczos_derivative { +public: + using Real = typename RandomAccessContainer::value_type; + lanczos_derivative(RandomAccessContainer const &v, + Real spacing = 1, + int filter_length = 18, + int approximation_order = 3) + : m_v{v}, dt{spacing} + { + BOOST_ASSERT_MSG(approximation_order <= 2 * filter_length, + "The approximation order must be <= 2n"); + BOOST_ASSERT_MSG(approximation_order >= 2, + "The approximation order must be >= 2"); + BOOST_ASSERT_MSG(spacing > 0, "Spacing between samples must be > 0."); + using std::size; + BOOST_ASSERT_MSG(size(v) >= filter_length, + "Vector must be at least as long as the filter length"); + m_f = detail::interior_filter(filter_length, approximation_order); + + boundary_filters.resize(filter_length); + for (size_t i = 0; i < filter_length; ++i) + { + // s = i - n; + boundary_filters[i] = detail::boundary_filter( + filter_length, approximation_order, i - filter_length); + } + } + + void reset_data(RandomAccessContainer const &v) + { + using std::size; + BOOST_ASSERT_MSG(size(v) >= m_f.size(), "Vector must be at least as long as the filter length"); + m_v = v; + } + + void reset_spacing(Real spacing) + { + BOOST_ASSERT_MSG(spacing > 0, "Spacing between samples must be > 0."); + dt = spacing; + } + + Real spacing() const + { + return dt; + } + + Real operator[](size_t i) const + { + using std::size; + if (i >= m_f.size() - 1 && i <= size(m_v) - m_f.size()) + { + Real dv = 0; + for (size_t j = 1; j < m_f.size(); ++j) + { + dv += m_f[j] * (m_v[i + j] - m_v[i - j]); + } + return dv / dt; + } + + // m_f.size() = N+1 + if (i < m_f.size() - 1) + { + auto &f = boundary_filters[i]; + Real dv = 0; + for (size_t j = 0; j < f.size(); ++j) { + dv += f[j] * m_v[j]; + } + return dv / dt; + } + + if (i > size(m_v) - m_f.size() && i < size(m_v)) + { + int k = size(m_v) - 1 - i; + auto &f = boundary_filters[k]; + Real dv = 0; + for (size_t j = 0; j < f.size(); ++j) + { + dv += f[j] * m_v[m_v.size() - 1 - j]; + } + return -dv / dt; + } + + // OOB access: + BOOST_ASSERT_MSG(false, "Out of bounds access in Lanczos derivative"); + return std::numeric_limits::quiet_NaN(); + } + +private: + const RandomAccessContainer &m_v; + std::vector m_f; + std::vector> boundary_filters; + Real dt; +}; + +// We can also implement lanczos_acceleration, but let's get the API for lanczos_derivative nailed down before doing so. + +} // namespace differentiation +} // namespace math +} // namespace boost +#endif diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 911572d66..ac7926265 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -902,6 +902,7 @@ test-suite misc : [ run test_constant_generate.cpp : : : release USE_CPP_FLOAT=1 off:no ] [ run test_cubic_b_spline.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_smart_ptr cxx11_defaulted_functions ] off msvc:/bigobj release ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] # does not in fact require C++17 constexpr; requires C++17 std::size. + [ run lanczos_smoothing_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] ] [ run test_real_concept.cpp ../../test/build//boost_unit_test_framework ] [ run test_remez.cpp pch ../../test/build//boost_unit_test_framework ] [ run test_roots.cpp pch ../../test/build//boost_unit_test_framework ] diff --git a/test/lanczos_smoothing_test.cpp b/test/lanczos_smoothing_test.cpp new file mode 100644 index 000000000..fd7315e88 --- /dev/null +++ b/test/lanczos_smoothing_test.cpp @@ -0,0 +1,469 @@ +/* + * Copyright Nick Thompson, 2017 + * 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) + */ +#define BOOST_TEST_MODULE lanczos_smoothing_test + +#include +#include +#include +#include +#include +#include + +using std::abs; +using std::pow; +using std::sqrt; +using boost::multiprecision::cpp_bin_float_50; +using boost::multiprecision::cpp_bin_float_100; +using boost::math::differentiation::lanczos_derivative; +using boost::math::differentiation::detail::discrete_legendre; +using boost::math::differentiation::detail::interior_filter; +using boost::math::differentiation::detail::boundary_filter; + +template +void test_dlp_norms() +{ + std::cout << "Testing Discrete Legendre Polynomial norms on type " << typeid(Real).name() << "\n"; + Real tol = std::numeric_limits::epsilon(); + auto dlp = discrete_legendre(1); + BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(0), 3, tol); + BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(1), 2, tol); + dlp = discrete_legendre(2); + BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(0), Real(5)/Real(2), tol); + BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(1), Real(5)/Real(4), tol); + BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(2), Real(3*3*7)/Real(pow(2,6)), 2*tol); + dlp = discrete_legendre(200); + for(size_t r = 0; r < 10; ++r) + { + Real calc = dlp.norm_sq(r); + Real expected = Real(2)/Real(2*r+1); + // As long as r << n, ||q_r||^2 -> 2/(2r+1) as n->infty + BOOST_CHECK_CLOSE_FRACTION(calc, expected, 0.05); + } + +} + +template +void test_dlp_evaluation() +{ + std::cout << "Testing evaluation of Discrete Legendre polynomials on type " << typeid(Real).name() << "\n"; + Real tol = std::numeric_limits::epsilon(); + size_t n = 25; + auto dlp = discrete_legendre(n); + Real x = 0.72; + Real q0 = dlp(x, 0); + BOOST_TEST(q0 == 1); + Real q1 = dlp(x, 1); + BOOST_TEST(q1 == x); + Real q2 = dlp(x, 2); + int N = 2*n+1; + Real expected = 0.5*(3*x*x - Real(N*N - 1)/Real(4*n*n)); + BOOST_CHECK_CLOSE_FRACTION(q2, expected, tol); + Real q3 = dlp(x, 3); + expected = (x/3)*(5*expected - (Real(N*N - 4))/(2*n*n)); + BOOST_CHECK_CLOSE_FRACTION(q3, expected, 2*tol); + + // q_r(x) is even for even r, and odd for odd r: + for (size_t n = 8; n < 22; ++n) + { + dlp = discrete_legendre(n); + for(size_t r = 2; r <= n; ++r) + { + if (r & 1) + { + Real q1 = dlp(x, r); + Real q2 = -dlp(-x, r); + BOOST_CHECK_CLOSE_FRACTION(q1, q2, tol); + } + else + { + Real q1 = dlp(x, r); + Real q2 = dlp(-x, r); + BOOST_CHECK_CLOSE_FRACTION(q1, q2, tol); + } + + Real l2_sq = 0; + for (int j = -(int)n; j <= (int) n; ++j) + { + Real y = Real(j)/Real(n); + Real term = dlp(y, r); + l2_sq += term*term; + } + l2_sq /= n; + Real l2_sq_expected = dlp.norm_sq(r); + BOOST_CHECK_CLOSE_FRACTION(l2_sq, l2_sq_expected, 20*tol); + } + } +} + +template +void test_dlp_next() +{ + std::cout << "Testing Discrete Legendre polynomial 'next' function on type " << typeid(Real).name() << "\n"; + Real tol = std::numeric_limits::epsilon(); + + for(size_t n = 2; n < 20; ++n) + { + auto dlp = discrete_legendre(n); + for(Real x = -1; x <= 1; x += 0.1) + { + dlp.initialize_recursion(x); + for (size_t k = 2; k < n; ++k) + { + BOOST_CHECK_CLOSE(dlp.next(), dlp(x, k), tol); + } + + dlp.initialize_recursion(x); + for (size_t k = 2; k < n; ++k) + { + BOOST_CHECK_CLOSE(dlp.next_prime(), dlp.prime(x, k), tol); + } + } + } +} + + +template +void test_dlp_derivatives() +{ + std::cout << "Testing Discrete Legendre polynomial derivatives on type " << typeid(Real).name() << "\n"; + Real tol = 10*std::numeric_limits::epsilon(); + int n = 25; + auto dlp = discrete_legendre(n); + Real x = 0.72; + Real q0p = dlp.prime(x, 0); + BOOST_TEST(q0p == 0); + Real q1p = dlp.prime(x, 1); + BOOST_TEST(q1p == 1); + Real q2p = dlp.prime(x, 2); + Real expected = 3*x; + BOOST_CHECK_CLOSE_FRACTION(q2p, expected, tol); +} + +template +void test_interior_filter() +{ + std::cout << "Testing interior filter on type " << typeid(Real).name() << "\n"; + Real tol = std::numeric_limits::epsilon(); + for(int n = 1; n < 10; ++n) + { + for (int p = 1; p < n; p += 2) + { + auto f = interior_filter(n,p); + // Since we only store half the filter coefficients, + // we need to reindex the moment sums: + Real sum = 0; + for (int j = 0; j < f.size(); ++j) + { + sum += j*f[j]; + } + BOOST_CHECK_CLOSE_FRACTION(2*sum, 1, 1000*tol); + + for (int l = 3; l <= p; l += 2) + { + sum = 0; + for (int j = 0; j < f.size(); ++j) + { + // The condition number of this sum is infinite! + // No need to get to worked up about the tolerance. + sum += pow(Real(j), l)*f[j]; + } + BOOST_CHECK_SMALL(sum, sqrt(tol)/100); + } + //std::cout << "(n,p) = (" << n << "," << p << ") = {"; + //for (auto & x : f) + //{ + // std::cout << x << ", "; + //} + //std::cout << "}\n"; + } + } +} + +template +void test_interior_lanczos() +{ + std::cout << "Testing interior Lanczos on type " << typeid(Real).name() << "\n"; + Real tol = std::numeric_limits::epsilon(); + std::vector v(500); + for (auto & x : v) + { + x = 7; + } + for (size_t n = 1; n < 10; ++n) + { + for (size_t p = 2; p < 2*n; p += 2) + { + auto lsd = lanczos_derivative(v, 0.1, n, p); + for (size_t m = n; m < v.size() - n; ++m) + { + Real dvdt = lsd[m]; + BOOST_CHECK_SMALL(dvdt, tol); + } + } + } + + + for(size_t i = 0; i < v.size(); ++i) + { + v[i] = 7*i+8; + } + + for (size_t n = 1; n < 10; ++n) + { + for (size_t p = 2; p < 2*n; p += 2) + { + auto lsd = lanczos_derivative(v, Real(1), n, p); + for (size_t m = n; m < v.size() - n; ++m) + { + Real dvdt = lsd[m]; + BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, 2000*tol); + } + } + } + + //std::random_device rd{}; + //auto seed = rd(); + //std::cout << "Seed = " << seed << "\n"; + std::mt19937 gen(4172378669); + std::normal_distribution<> dis{0, 0.01}; + std::cout << std::fixed; + for (size_t i = 0; i < v.size(); ++i) + { + v[i] = 7*i+8 + dis(gen); + } + + for (size_t n = 1; n < 10; ++n) + { + for (size_t p = 2; p < 2*n; p += 2) + { + auto lsd = lanczos_derivative(v, Real(1), n, p); + for (size_t m = n; m < v.size() - n; ++m) + { + BOOST_CHECK_CLOSE_FRACTION(lsd[m], Real(7), Real(0.0042)); + } + } + } + + + for (size_t i = 0; i < v.size(); ++i) + { + v[i] = 15*i*i + 7*i+8 + dis(gen); + } + + for (size_t n = 1; n < 10; ++n) + { + for (size_t p = 2; p < 2*n; p += 2) + { + auto lsd = lanczos_derivative(v, Real(1), n, p); + for (size_t m = n; m < v.size() - n; ++m) + { + BOOST_CHECK_CLOSE_FRACTION(lsd[m], Real(30*m + 7), Real(0.00008)); + } + } + } + + std::normal_distribution<> dis1{0, 0.0001}; + Real omega = Real(1)/Real(16); + for (size_t i = 0; i < v.size(); ++i) + { + v[i] = sin(i*omega) + dis1(gen); + } + + for (size_t n = 10; n < 20; ++n) + { + for (size_t p = 3; p < 100 && p < n/2; p += 2) + { + auto lsd = lanczos_derivative(v, Real(1), n, p); + + for (size_t m = n; m < v.size() - n && m < n + 10; ++m) + { + BOOST_CHECK_CLOSE_FRACTION(lsd[m], omega*cos(omega*m), Real(0.03)); + } + } + } +} + +template +void test_boundary_filters() +{ + std::cout << "Testing boundary filters on type " << typeid(Real).name() << "\n"; + Real tol = std::numeric_limits::epsilon(); + for(int n = 1; n < 5; ++n) + { + for (int p = 1; p < 2*n+1; ++p) + { + for (int s = -n; s <= n; ++s) + { + auto f = boundary_filter(n, p, s); + // Sum is zero: + Real sum = 0; + Real c = 0; + for (auto & x : f) + { + Real y = x - c; + Real t = sum + y; + c = (t-sum) -y; + sum = t; + } + BOOST_CHECK_SMALL(sum, 200*tol); + + sum = 0; + c = 0; + for (int k = 0; k < f.size(); ++k) + { + int j = k - n; + // note the shifted index here: + Real x = (j-s)*f[k]; + Real y = x - c; + Real t = sum + y; + c = (t-sum) -y; + sum = t; + } + BOOST_CHECK_CLOSE_FRACTION(sum, 1, 350*tol); + + + for (int l = 2; l <= p; ++l) + { + sum = 0; + c = 0; + for (int k = 0; k < f.size(); ++k) + { + int j = k - n; + // The condition number of this sum is infinite! + // No need to get to worked up about the tolerance. + Real x = pow(Real(j-s), l)*f[k]; + Real y = x - c; + Real t = sum + y; + c = (t-sum) -y; + sum = t; + } + BOOST_CHECK_SMALL(sum, sqrt(tol)/10); + } + + //std::cout << "(n,p,s) = ("<< n << ", " << p << "," << s << ") = {"; + //for (auto & x : f) + //{ + // std::cout << x << ", "; + //} + //std::cout << "}\n";*/ + } + } + } +} + +template +void test_boundary_lanczos() +{ + std::cout << "Testing Lanczos boundary on type " << typeid(Real).name() << "\n"; + Real tol = std::numeric_limits::epsilon(); + std::vector v(500); + for (auto & x : v) + { + x = 7; + } + for (size_t n = 1; n < 10; ++n) + { + for (size_t p = 2; p < 2*n; ++p) + { + auto lsd = lanczos_derivative(v, 0.0125, n, p); + for (size_t m = 0; m < n; ++m) + { + Real dvdt = lsd[m]; + BOOST_CHECK_SMALL(dvdt, 4*sqrt(tol)); + } + for (size_t m = v.size() - n; m < v.size(); ++m) + { + Real dvdt = lsd[m]; + BOOST_CHECK_SMALL(dvdt, 4*sqrt(tol)); + } + } + } + + for(size_t i = 0; i < v.size(); ++i) + { + v[i] = 7*i+8; + } + + for (size_t n = 3; n < 10; ++n) + { + for (size_t p = 2; p < 2*n; ++p) + { + auto lsd = lanczos_derivative(v, Real(1), n, p); + for (size_t m = 0; m < n; ++m) + { + Real dvdt = lsd[m]; + BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, sqrt(tol)); + } + + for (size_t m = v.size() - n; m < v.size(); ++m) + { + Real dvdt = lsd[m]; + BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, 4*sqrt(tol)); + } + } + } + + for (size_t i = 0; i < v.size(); ++i) + { + v[i] = 15*i*i + 7*i+8; + } + + for (size_t n = 1; n < 10; ++n) + { + for (size_t p = 2; p < 2*n; ++p) + { + auto lsd = lanczos_derivative(v, Real(1), n, p); + for (size_t m = 0; m < v.size(); ++m) + { + BOOST_CHECK_CLOSE_FRACTION(lsd[m], 30*m+7, 30*sqrt(tol)); + } + } + } + + // Demonstrate that the boundary filters are also denoising: + //std::random_device rd{}; + //auto seed = rd(); + //std::cout << "seed = " << seed << "\n"; + std::mt19937 gen(311354333); + std::normal_distribution<> dis{0, 0.01}; + for (size_t i = 0; i < v.size(); ++i) + { + v[i] += dis(gen); + } + + for (size_t n = 1; n < 10; ++n) + { + for (size_t p = 2; p < n; ++p) + { + auto lsd = lanczos_derivative(v, Real(1), n, p); + for (size_t m = 0; m < v.size(); ++m) + { + BOOST_CHECK_CLOSE_FRACTION(lsd[m], 30*m+7, 0.005); + } + } + } + +} + +BOOST_AUTO_TEST_CASE(lanczos_smoothing_test) +{ + test_dlp_norms(); + test_dlp_evaluation(); + test_dlp_derivatives(); + test_dlp_next(); + test_dlp_norms(); + test_boundary_filters(); + test_boundary_filters(); + test_boundary_filters(); + test_boundary_lanczos(); + test_boundary_lanczos(); + test_boundary_lanczos(); + + test_interior_filter(); + test_interior_filter(); + test_interior_lanczos(); +}