Lanczos smoothing differentiators.

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 <boost/math/differentiation/lanczos_smoothing.hpp>
+
+namespace boost { namespace math { namespace differentiation {
+
+    template <class RandomAccessContainer>
+    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<double> 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<double> v(500);
+    // fill v with noisy data.
+    auto lanczos = lanczos_derivative(v, spacing);
+    // use lanczos with v . . .
+    std::vector<double> 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]

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]

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 <vector>
+#include <boost/assert.hpp>
+
+namespace boost {
+namespace math {
+namespace differentiation {
+
+namespace detail {
+template <typename Real>
+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 <class Real>
+std::vector<Real> 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<Real> f(n + 1, 0);
+    auto dlp = discrete_legendre<Real>(n);
+    std::vector<Real> coeffs(p+1, std::numeric_limits<Real>::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 <class Real>
+std::vector<Real> boundary_filter(int n, int p, int s) {
+    std::vector<Real> f(2 * n + 1, 0);
+    auto dlp = discrete_legendre<Real>(n);
+    Real sn = Real(s) / Real(n);
+    std::vector<Real> coeffs(p+1, std::numeric_limits<Real>::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 RandomAccessContainer>
+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<Real>(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<Real>(
+                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<Real>::quiet_NaN();
+    }
+
+private:
+    const RandomAccessContainer &m_v;
+    std::vector<Real> m_f;
+    std::vector<std::vector<Real>> 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

test/Jamfile.v2

@@ -902,6 +902,7 @@ test-suite misc :
    [ run test_constant_generate.cpp : : : release <define>USE_CPP_FLOAT=1 <exception-handling>off:<build>no  ]
    [ run test_cubic_b_spline.cpp ../../test/build//boost_unit_test_framework : : :  [ requires cxx11_smart_ptr cxx11_defaulted_functions ] <debug-symbols>off <toolset>msvc:<cxxflags>/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  ]

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 <random>
+#include <boost/accumulators/statistics/sum_kahan.hpp>
+#include <boost/test/included/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+#include <boost/math/differentiation/lanczos_smoothing.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+
+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<class Real>
+void test_dlp_norms()
+{
+    std::cout << "Testing Discrete Legendre Polynomial norms on type " << typeid(Real).name() << "\n";
+    Real tol = std::numeric_limits<Real>::epsilon();
+    auto dlp = discrete_legendre<Real>(1);
+    BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(0), 3, tol);
+    BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(1), 2, tol);
+    dlp = discrete_legendre<Real>(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<Real>(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<class Real>
+void test_dlp_evaluation()
+{
+    std::cout << "Testing evaluation of Discrete Legendre polynomials on type " << typeid(Real).name() << "\n";
+    Real tol = std::numeric_limits<Real>::epsilon();
+    size_t n = 25;
+    auto dlp = discrete_legendre<Real>(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<Real>(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<class Real>
+void test_dlp_next()
+{
+    std::cout << "Testing Discrete Legendre polynomial 'next' function on type " << typeid(Real).name() << "\n";
+    Real tol = std::numeric_limits<Real>::epsilon();
+
+    for(size_t n = 2; n < 20; ++n)
+    {
+        auto dlp = discrete_legendre<Real>(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<class Real>
+void test_dlp_derivatives()
+{
+    std::cout << "Testing Discrete Legendre polynomial derivatives on type " << typeid(Real).name() << "\n";
+    Real tol = 10*std::numeric_limits<Real>::epsilon();
+    int n = 25;
+    auto dlp = discrete_legendre<Real>(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<class Real>
+void test_interior_filter()
+{
+    std::cout << "Testing interior filter on type " << typeid(Real).name() << "\n";
+    Real tol = std::numeric_limits<Real>::epsilon();
+    for(int n = 1; n < 10; ++n)
+    {
+        for (int p = 1; p < n; p += 2)
+        {
+            auto f = interior_filter<Real>(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<class Real>
+void test_interior_lanczos()
+{
+    std::cout << "Testing interior Lanczos on type " << typeid(Real).name() << "\n";
+    Real tol = std::numeric_limits<Real>::epsilon();
+    std::vector<Real> 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<class Real>
+void test_boundary_filters()
+{
+    std::cout << "Testing boundary filters on type " << typeid(Real).name() << "\n";
+    Real tol = std::numeric_limits<Real>::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<Real>(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<class Real>
+void test_boundary_lanczos()
+{
+    std::cout << "Testing Lanczos boundary on type " << typeid(Real).name() << "\n";
+    Real tol = std::numeric_limits<Real>::epsilon();
+    std::vector<Real> 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<double>();
+    test_dlp_evaluation<double>();
+    test_dlp_derivatives<double>();
+    test_dlp_next<double>();
+    test_dlp_norms<cpp_bin_float_50>();
+    test_boundary_filters<double>();
+    test_boundary_filters<long double>();
+    test_boundary_filters<cpp_bin_float_50>();
+    test_boundary_lanczos<double>();
+    test_boundary_lanczos<long double>();
+    test_boundary_lanczos<cpp_bin_float_50>();
+
+    test_interior_filter<double>();
+    test_interior_filter<long double>();
+    test_interior_lanczos<double>();
+}