From a250691ed190cac0bf261eee92a10fed43fec1f4 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Thu, 27 Dec 2018 18:59:44 -0700 Subject: [PATCH 01/10] Begin rearrangement. --- .../numerical_differentiation.qbk | 10 +- .../differentiation/finite_difference.hpp | 266 +++++++++++++++++ .../math/tools/numerical_differentiation.hpp | 267 +----------------- test/test_numerical_differentiation.cpp | 6 +- 4 files changed, 284 insertions(+), 265 deletions(-) create mode 100644 include/boost/math/differentiation/finite_difference.hpp diff --git a/doc/differentiation/numerical_differentiation.qbk b/doc/differentiation/numerical_differentiation.qbk index 857627bd4..6f7956b8d 100644 --- a/doc/differentiation/numerical_differentiation.qbk +++ b/doc/differentiation/numerical_differentiation.qbk @@ -10,10 +10,10 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) [heading Synopsis] `` -#include +#include -namespace boost::math::tools { +namespace boost { namespace math { namespace differentiation { template Real complex_step_derivative(const F f, Real x); @@ -21,19 +21,19 @@ namespace boost::math::tools { template Real finite_difference_derivative(const F f, Real x, Real* error = nullptr); -} // namespaces +}}} // namespaces `` [heading Description] -The function `finite_difference_derivative` calculates a finite-difference approximation to the derivative of of a function /f/ at point /x/. +The function `finite_difference_derivative` calculates a finite-difference approximation to the derivative of a function /f/ at point /x/. A basic usage is auto f = [](double x) { return std::exp(x); }; double x = 1.7; double dfdx = finite_difference_derivative(f, x); -Finite differencing is complicated, as finite-difference approximations to the derivative are /infinitely/ ill-conditioned. +Finite differencing is complicated, as differentiation is /infinitely/ ill-conditioned. In addition, for any function implemented in finite-precision arithmetic, the "true" derivative is /zero/ almost everywhere, and undefined at representables. However, some tricks allow for reasonable results to be obtained in many cases. diff --git a/include/boost/math/differentiation/finite_difference.hpp b/include/boost/math/differentiation/finite_difference.hpp new file mode 100644 index 000000000..c7615b6b2 --- /dev/null +++ b/include/boost/math/differentiation/finite_difference.hpp @@ -0,0 +1,266 @@ +// (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_FINITE_DIFFERENCE_HPP +#define BOOST_MATH_DIFFERENTIATION_FINITE_DIFFERENCE_HPP + +/* + * Performs numerical differentiation by finite-differences. + * + * All numerical differentiation using finite-differences are ill-conditioned, and these routines are no exception. + * A simple argument demonstrates that the error is unbounded as h->0. + * Take the one sides finite difference formula f'(x) = (f(x+h)-f(x))/h. + * The evaluation of f induces an error as well as the error from the finite-difference approximation, giving + * |f'(x) - (f(x+h) -f(x))/h| < h|f''(x)|/2 + (|f(x)|+|f(x+h)|)eps/h =: g(h), where eps is the unit roundoff for the type. + * It is reasonable to choose h in a way that minimizes the maximum error bound g(h). + * The value of h that minimizes g is h = sqrt(2eps(|f(x)| + |f(x+h)|)/|f''(x)|), and for this value of h the error bound is + * sqrt(2eps(|f(x+h) +f(x)||f''(x)|)). + * In fact it is not necessary to compute the ratio (|f(x+h)| + |f(x)|)/|f''(x)|; the error bound of ~\sqrt{\epsilon} still holds if we set it to one. + * + * + * For more details on this method of analysis, see + * + * http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf + * http://web.archive.org/web/20150420195907/http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf + * + * + * It can be shown on general grounds that when choosing the optimal h, the maximum error in f'(x) is ~(|f(x)|eps)^k/k+1|f^(k-1)(x)|^1/k+1. + * From this we can see that full precision can be recovered in the limit k->infinity. + * + * References: + * + * 1) Fornberg, Bengt. "Generation of finite difference formulas on arbitrarily spaced grids." Mathematics of computation 51.184 (1988): 699-706. + * + * + * The second algorithm, the complex step derivative, is not ill-conditioned. + * However, it requires that your function can be evaluated at complex arguments. + * The idea is that f(x+ih) = f(x) +ihf'(x) - h^2f''(x) + ... so f'(x) \approx Im[f(x+ih)]/h. + * No subtractive cancellation occurs. The error is ~ eps|f'(x)| + eps^2|f'''(x)|/6; hard to beat that. + * + * References: + * + * 1) Squire, William, and George Trapp. "Using complex variables to estimate derivatives of real functions." Siam Review 40.1 (1998): 110-112. + */ + +#include +#include + +namespace boost{ namespace math{ namespace differentiation { + +namespace detail { + template + Real make_xph_representable(Real x, Real h) + { + using std::numeric_limits; + // Redefine h so that x + h is representable. Not using this trick leads to large error. + // The compiler flag -ffast-math evaporates these operations . . . + Real temp = x + h; + h = temp - x; + // Handle the case x + h == x: + if (h == 0) + { + h = boost::math::nextafter(x, (numeric_limits::max)()) - x; + } + return h; + } +} + +template +Real complex_step_derivative(const F f, Real x) +{ + // Is it really this easy? Yes. + // Note that some authors recommend taking the stepsize h to be smaller than epsilon(), some recommending use of the min(). + // This idea was tested over a few billion test cases and found the make the error *much* worse. + // Even 2eps and eps/2 made the error worse, which was surprising. + using std::complex; + using std::numeric_limits; + constexpr const Real step = (numeric_limits::epsilon)(); + constexpr const Real inv_step = 1/(numeric_limits::epsilon)(); + return f(complex(x, step)).imag()*inv_step; +} + +namespace detail { + + template + struct fd_tag {}; + + template + Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<1>&) + { + using std::sqrt; + using std::pow; + using std::abs; + using std::numeric_limits; + + const Real eps = (numeric_limits::epsilon)(); + // Error bound ~eps^1/2 + // Note that this estimate of h differs from the best estimate by a factor of sqrt((|f(x)| + |f(x+h)|)/|f''(x)|). + // Since this factor is invariant under the scaling f -> kf, then we are somewhat justified in approximating it by 1. + // This approximation will get better as we move to higher orders of accuracy. + Real h = 2 * sqrt(eps); + h = detail::make_xph_representable(x, h); + + Real yh = f(x + h); + Real y0 = f(x); + Real diff = yh - y0; + if (error) + { + Real ym = f(x - h); + Real ypph = abs(yh - 2 * y0 + ym) / h; + // h*|f''(x)|*0.5 + (|f(x+h)+|f(x)|)*eps/h + *error = ypph / 2 + (abs(yh) + abs(y0))*eps / h; + } + return diff / h; + } + + template + Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<2>&) + { + using std::sqrt; + using std::pow; + using std::abs; + using std::numeric_limits; + + const Real eps = (numeric_limits::epsilon)(); + // Error bound ~eps^2/3 + // See the previous discussion to understand determination of h and the error bound. + // Series[(f[x+h] - f[x-h])/(2*h), {h, 0, 4}] + Real h = pow(3 * eps, static_cast(1) / static_cast(3)); + h = detail::make_xph_representable(x, h); + + Real yh = f(x + h); + Real ymh = f(x - h); + Real diff = yh - ymh; + if (error) + { + Real yth = f(x + 2 * h); + Real ymth = f(x - 2 * h); + *error = eps * (abs(yh) + abs(ymh)) / (2 * h) + abs((yth - ymth) / 2 - diff) / (6 * h); + } + + return diff / (2 * h); + } + + template + Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<4>&) + { + using std::sqrt; + using std::pow; + using std::abs; + using std::numeric_limits; + + const Real eps = (numeric_limits::epsilon)(); + // Error bound ~eps^4/5 + Real h = pow(11.25*eps, (Real)1 / (Real)5); + h = detail::make_xph_representable(x, h); + Real ymth = f(x - 2 * h); + Real yth = f(x + 2 * h); + Real yh = f(x + h); + Real ymh = f(x - h); + Real y2 = ymth - yth; + Real y1 = yh - ymh; + if (error) + { + // Mathematica code to extract the remainder: + // Series[(f[x-2*h]+ 8*f[x+h] - 8*f[x-h] - f[x+2*h])/(12*h), {h, 0, 7}] + Real y_three_h = f(x + 3 * h); + Real y_m_three_h = f(x - 3 * h); + // Error from fifth derivative: + *error = abs((y_three_h - y_m_three_h) / 2 + 2 * (ymth - yth) + 5 * (yh - ymh) / 2) / (30 * h); + // Error from function evaluation: + *error += eps * (abs(yth) + abs(ymth) + 8 * (abs(ymh) + abs(yh))) / (12 * h); + } + return (y2 + 8 * y1) / (12 * h); + } + + template + Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<6>&) + { + using std::sqrt; + using std::pow; + using std::abs; + using std::numeric_limits; + + const Real eps = (numeric_limits::epsilon)(); + // Error bound ~eps^6/7 + // Error: h^6f^(7)(x)/140 + 5|f(x)|eps/h + Real h = pow(eps / 168, (Real)1 / (Real)7); + h = detail::make_xph_representable(x, h); + + Real yh = f(x + h); + Real ymh = f(x - h); + Real y1 = yh - ymh; + Real y2 = f(x - 2 * h) - f(x + 2 * h); + Real y3 = f(x + 3 * h) - f(x - 3 * h); + + if (error) + { + // Mathematica code to generate fd scheme for 7th derivative: + // Sum[(-1)^i*Binomial[7, i]*(f[x+(3-i)*h] + f[x+(4-i)*h])/2, {i, 0, 7}] + // Mathematica to demonstrate that this is a finite difference formula for 7th derivative: + // Series[(f[x+4*h]-f[x-4*h] + 6*(f[x-3*h] - f[x+3*h]) + 14*(f[x-h] - f[x+h] + f[x+2*h] - f[x-2*h]))/2, {h, 0, 15}] + Real y7 = (f(x + 4 * h) - f(x - 4 * h) - 6 * y3 - 14 * y1 - 14 * y2) / 2; + *error = abs(y7) / (140 * h) + 5 * (abs(yh) + abs(ymh))*eps / h; + } + return (y3 + 9 * y2 + 45 * y1) / (60 * h); + } + + template + Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<8>&) + { + using std::sqrt; + using std::pow; + using std::abs; + using std::numeric_limits; + + const Real eps = (numeric_limits::epsilon)(); + // Error bound ~eps^8/9. + // In double precision, we only expect to lose two digits of precision while using this formula, at the cost of 8 function evaluations. + // Error: h^8|f^(9)(x)|/630 + 7|f(x)|eps/h assuming 7 unstabilized additions. + // Mathematica code to get the error: + // Series[(f[x+h]-f[x-h])*(4/5) + (1/5)*(f[x-2*h] - f[x+2*h]) + (4/105)*(f[x+3*h] - f[x-3*h]) + (1/280)*(f[x-4*h] - f[x+4*h]), {h, 0, 9}] + // If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h. + Real h = pow(551.25*eps, (Real)1 / (Real)9); + h = detail::make_xph_representable(x, h); + + Real yh = f(x + h); + Real ymh = f(x - h); + Real y1 = yh - ymh; + Real y2 = f(x - 2 * h) - f(x + 2 * h); + Real y3 = f(x + 3 * h) - f(x - 3 * h); + Real y4 = f(x - 4 * h) - f(x + 4 * h); + + Real tmp1 = 3 * y4 / 8 + 4 * y3; + Real tmp2 = 21 * y2 + 84 * y1; + + if (error) + { + // Mathematica code to generate fd scheme for 7th derivative: + // Sum[(-1)^i*Binomial[9, i]*(f[x+(4-i)*h] + f[x+(5-i)*h])/2, {i, 0, 9}] + // Mathematica to demonstrate that this is a finite difference formula for 7th derivative: + // Series[(f[x+5*h]-f[x- 5*h])/2 + 4*(f[x-4*h] - f[x+4*h]) + 27*(f[x+3*h] - f[x-3*h])/2 + 24*(f[x-2*h] - f[x+2*h]) + 21*(f[x+h] - f[x-h]), {h, 0, 15}] + Real f9 = (f(x + 5 * h) - f(x - 5 * h)) / 2 + 4 * y4 + 27 * y3 / 2 + 24 * y2 + 21 * y1; + *error = abs(f9) / (630 * h) + 7 * (abs(yh) + abs(ymh))*eps / h; + } + return (tmp1 + tmp2) / (105 * h); + } + + template + Real finite_difference_derivative(const F f, Real x, Real* error, const tag&) + { + // Always fails, but condition is template-arg-dependent so only evaluated if we get instantiated. + BOOST_STATIC_ASSERT_MSG(sizeof(Real) == 0, "Finite difference not implemented for this order: try 1, 2, 4, 6 or 8"); + } + +} + +template +inline Real finite_difference_derivative(const F f, Real x, Real* error = nullptr) +{ + return detail::finite_difference_derivative(f, x, error, detail::fd_tag()); +} + +}}} // namespaces +#endif diff --git a/include/boost/math/tools/numerical_differentiation.hpp b/include/boost/math/tools/numerical_differentiation.hpp index 34fef0db8..eb4bdcb44 100644 --- a/include/boost/math/tools/numerical_differentiation.hpp +++ b/include/boost/math/tools/numerical_differentiation.hpp @@ -2,266 +2,19 @@ // 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_TOOLS_NUMERICAL_DIFFERENTIATION_HPP #define BOOST_MATH_TOOLS_NUMERICAL_DIFFERENTIATION_HPP +#include +#include -/* - * Performs numerical differentiation by finite-differences. - * - * All numerical differentiation using finite-differences are ill-conditioned, and these routines are no exception. - * A simple argument demonstrates that the error is unbounded as h->0. - * Take the one sides finite difference formula f'(x) = (f(x+h)-f(x))/h. - * The evaluation of f induces an error as well as the error from the finite-difference approximation, giving - * |f'(x) - (f(x+h) -f(x))/h| < h|f''(x)|/2 + (|f(x)|+|f(x+h)|)eps/h =: g(h), where eps is the unit roundoff for the type. - * It is reasonable to choose h in a way that minimizes the maximum error bound g(h). - * The value of h that minimizes g is h = sqrt(2eps(|f(x)| + |f(x+h)|)/|f''(x)|), and for this value of h the error bound is - * sqrt(2eps(|f(x+h) +f(x)||f''(x)|)). - * In fact it is not necessary to compute the ratio (|f(x+h)| + |f(x)|)/|f''(x)|; the error bound of ~\sqrt{\epsilon} still holds if we set it to one. - * - * - * For more details on this method of analysis, see - * - * http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf - * http://web.archive.org/web/20150420195907/http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf - * - * - * It can be shown on general grounds that when choosing the optimal h, the maximum error in f'(x) is ~(|f(x)|eps)^k/k+1|f^(k-1)(x)|^1/k+1. - * From this we can see that full precision can be recovered in the limit k->infinity. - * - * References: - * - * 1) Fornberg, Bengt. "Generation of finite difference formulas on arbitrarily spaced grids." Mathematics of computation 51.184 (1988): 699-706. - * - * - * The second algorithm, the complex step derivative, is not ill-conditioned. - * However, it requires that your function can be evaluated at complex arguments. - * The idea is that f(x+ih) = f(x) +ihf'(x) - h^2f''(x) + ... so f'(x) \approx Im[f(x+ih)]/h. - * No subtractive cancellation occurs. The error is ~ eps|f'(x)| + eps^2|f'''(x)|/6; hard to beat that. - * - * References: - * - * 1) Squire, William, and George Trapp. "Using complex variables to estimate derivatives of real functions." Siam Review 40.1 (1998): 110-112. - */ +BOOST_HEADER_DEPRECATED(""); -#include -#include - -// Make sure everyone is informed that C++17 is required: -namespace boost{ namespace math{ namespace tools { - -namespace detail { - template - Real make_xph_representable(Real x, Real h) - { - using std::numeric_limits; - // Redefine h so that x + h is representable. Not using this trick leads to large error. - // The compiler flag -ffast-math evaporates these operations . . . - Real temp = x + h; - h = temp - x; - // Handle the case x + h == x: - if (h == 0) - { - h = boost::math::nextafter(x, (numeric_limits::max)()) - x; - } - return h; - } +namespace boost { + namespace math { + namespace differentiation { + using finite_difference_derivative; + using complex_step_derivative; + } + } } - -template -Real complex_step_derivative(const F f, Real x) -{ - // Is it really this easy? Yes. - // Note that some authors recommend taking the stepsize h to be smaller than epsilon(), some recommending use of the min(). - // This idea was tested over a few billion test cases and found the make the error *much* worse. - // Even 2eps and eps/2 made the error worse, which was surprising. - using std::complex; - using std::numeric_limits; - constexpr const Real step = (numeric_limits::epsilon)(); - constexpr const Real inv_step = 1/(numeric_limits::epsilon)(); - return f(complex(x, step)).imag()*inv_step; -} - -namespace detail { - - template - struct fd_tag {}; - - template - Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<1>&) - { - using std::sqrt; - using std::pow; - using std::abs; - using std::numeric_limits; - - const Real eps = (numeric_limits::epsilon)(); - // Error bound ~eps^1/2 - // Note that this estimate of h differs from the best estimate by a factor of sqrt((|f(x)| + |f(x+h)|)/|f''(x)|). - // Since this factor is invariant under the scaling f -> kf, then we are somewhat justified in approximating it by 1. - // This approximation will get better as we move to higher orders of accuracy. - Real h = 2 * sqrt(eps); - h = detail::make_xph_representable(x, h); - - Real yh = f(x + h); - Real y0 = f(x); - Real diff = yh - y0; - if (error) - { - Real ym = f(x - h); - Real ypph = abs(yh - 2 * y0 + ym) / h; - // h*|f''(x)|*0.5 + (|f(x+h)+|f(x)|)*eps/h - *error = ypph / 2 + (abs(yh) + abs(y0))*eps / h; - } - return diff / h; - } - - template - Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<2>&) - { - using std::sqrt; - using std::pow; - using std::abs; - using std::numeric_limits; - - const Real eps = (numeric_limits::epsilon)(); - // Error bound ~eps^2/3 - // See the previous discussion to understand determination of h and the error bound. - // Series[(f[x+h] - f[x-h])/(2*h), {h, 0, 4}] - Real h = pow(3 * eps, static_cast(1) / static_cast(3)); - h = detail::make_xph_representable(x, h); - - Real yh = f(x + h); - Real ymh = f(x - h); - Real diff = yh - ymh; - if (error) - { - Real yth = f(x + 2 * h); - Real ymth = f(x - 2 * h); - *error = eps * (abs(yh) + abs(ymh)) / (2 * h) + abs((yth - ymth) / 2 - diff) / (6 * h); - } - - return diff / (2 * h); - } - - template - Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<4>&) - { - using std::sqrt; - using std::pow; - using std::abs; - using std::numeric_limits; - - const Real eps = (numeric_limits::epsilon)(); - // Error bound ~eps^4/5 - Real h = pow(11.25*eps, (Real)1 / (Real)5); - h = detail::make_xph_representable(x, h); - Real ymth = f(x - 2 * h); - Real yth = f(x + 2 * h); - Real yh = f(x + h); - Real ymh = f(x - h); - Real y2 = ymth - yth; - Real y1 = yh - ymh; - if (error) - { - // Mathematica code to extract the remainder: - // Series[(f[x-2*h]+ 8*f[x+h] - 8*f[x-h] - f[x+2*h])/(12*h), {h, 0, 7}] - Real y_three_h = f(x + 3 * h); - Real y_m_three_h = f(x - 3 * h); - // Error from fifth derivative: - *error = abs((y_three_h - y_m_three_h) / 2 + 2 * (ymth - yth) + 5 * (yh - ymh) / 2) / (30 * h); - // Error from function evaluation: - *error += eps * (abs(yth) + abs(ymth) + 8 * (abs(ymh) + abs(yh))) / (12 * h); - } - return (y2 + 8 * y1) / (12 * h); - } - - template - Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<6>&) - { - using std::sqrt; - using std::pow; - using std::abs; - using std::numeric_limits; - - const Real eps = (numeric_limits::epsilon)(); - // Error bound ~eps^6/7 - // Error: h^6f^(7)(x)/140 + 5|f(x)|eps/h - Real h = pow(eps / 168, (Real)1 / (Real)7); - h = detail::make_xph_representable(x, h); - - Real yh = f(x + h); - Real ymh = f(x - h); - Real y1 = yh - ymh; - Real y2 = f(x - 2 * h) - f(x + 2 * h); - Real y3 = f(x + 3 * h) - f(x - 3 * h); - - if (error) - { - // Mathematica code to generate fd scheme for 7th derivative: - // Sum[(-1)^i*Binomial[7, i]*(f[x+(3-i)*h] + f[x+(4-i)*h])/2, {i, 0, 7}] - // Mathematica to demonstrate that this is a finite difference formula for 7th derivative: - // Series[(f[x+4*h]-f[x-4*h] + 6*(f[x-3*h] - f[x+3*h]) + 14*(f[x-h] - f[x+h] + f[x+2*h] - f[x-2*h]))/2, {h, 0, 15}] - Real y7 = (f(x + 4 * h) - f(x - 4 * h) - 6 * y3 - 14 * y1 - 14 * y2) / 2; - *error = abs(y7) / (140 * h) + 5 * (abs(yh) + abs(ymh))*eps / h; - } - return (y3 + 9 * y2 + 45 * y1) / (60 * h); - } - - template - Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<8>&) - { - using std::sqrt; - using std::pow; - using std::abs; - using std::numeric_limits; - - const Real eps = (numeric_limits::epsilon)(); - // Error bound ~eps^8/9. - // In double precision, we only expect to lose two digits of precision while using this formula, at the cost of 8 function evaluations. - // Error: h^8|f^(9)(x)|/630 + 7|f(x)|eps/h assuming 7 unstabilized additions. - // Mathematica code to get the error: - // Series[(f[x+h]-f[x-h])*(4/5) + (1/5)*(f[x-2*h] - f[x+2*h]) + (4/105)*(f[x+3*h] - f[x-3*h]) + (1/280)*(f[x-4*h] - f[x+4*h]), {h, 0, 9}] - // If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h. - Real h = pow(551.25*eps, (Real)1 / (Real)9); - h = detail::make_xph_representable(x, h); - - Real yh = f(x + h); - Real ymh = f(x - h); - Real y1 = yh - ymh; - Real y2 = f(x - 2 * h) - f(x + 2 * h); - Real y3 = f(x + 3 * h) - f(x - 3 * h); - Real y4 = f(x - 4 * h) - f(x + 4 * h); - - Real tmp1 = 3 * y4 / 8 + 4 * y3; - Real tmp2 = 21 * y2 + 84 * y1; - - if (error) - { - // Mathematica code to generate fd scheme for 7th derivative: - // Sum[(-1)^i*Binomial[9, i]*(f[x+(4-i)*h] + f[x+(5-i)*h])/2, {i, 0, 9}] - // Mathematica to demonstrate that this is a finite difference formula for 7th derivative: - // Series[(f[x+5*h]-f[x- 5*h])/2 + 4*(f[x-4*h] - f[x+4*h]) + 27*(f[x+3*h] - f[x-3*h])/2 + 24*(f[x-2*h] - f[x+2*h]) + 21*(f[x+h] - f[x-h]), {h, 0, 15}] - Real f9 = (f(x + 5 * h) - f(x - 5 * h)) / 2 + 4 * y4 + 27 * y3 / 2 + 24 * y2 + 21 * y1; - *error = abs(f9) / (630 * h) + 7 * (abs(yh) + abs(ymh))*eps / h; - } - return (tmp1 + tmp2) / (105 * h); - } - - template - Real finite_difference_derivative(const F f, Real x, Real* error, const tag&) - { - // Always fails, but condition is template-arg-dependent so only evaluated if we get instantiated. - BOOST_STATIC_ASSERT_MSG(sizeof(Real) == 0, "Finite difference not implemented for this order: try 1, 2, 4, 6 or 8"); - } - -} - -template -inline Real finite_difference_derivative(const F f, Real x, Real* error = nullptr) -{ - return detail::finite_difference_derivative(f, x, error, detail::fd_tag()); -} - -}}} // namespaces #endif diff --git a/test/test_numerical_differentiation.cpp b/test/test_numerical_differentiation.cpp index 01f42e16e..d49a75daa 100644 --- a/test/test_numerical_differentiation.cpp +++ b/test/test_numerical_differentiation.cpp @@ -16,12 +16,12 @@ #include #include #include -#include +#include using std::abs; using std::pow; -using boost::math::tools::finite_difference_derivative; -using boost::math::tools::complex_step_derivative; +using boost::math::differentiation::finite_difference_derivative; +using boost::math::differentiation::complex_step_derivative; using boost::math::cyl_bessel_j; using boost::math::cyl_bessel_j_prime; using boost::math::constants::half; From a27c4bc96cfd757bfdb6371ce2be7c1c8483fcaa Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Thu, 27 Dec 2018 20:38:27 -0700 Subject: [PATCH 02/10] Also change concept and include test [CI SKIP] --- include/boost/math/tools/numerical_differentiation.hpp | 8 -------- .../numerical_differentiation_concept_test.cpp | 4 ++-- test/compile_test/numerical_differentiation_incl_test.cpp | 6 +++--- 3 files changed, 5 insertions(+), 13 deletions(-) diff --git a/include/boost/math/tools/numerical_differentiation.hpp b/include/boost/math/tools/numerical_differentiation.hpp index eb4bdcb44..52d278889 100644 --- a/include/boost/math/tools/numerical_differentiation.hpp +++ b/include/boost/math/tools/numerical_differentiation.hpp @@ -9,12 +9,4 @@ BOOST_HEADER_DEPRECATED(""); -namespace boost { - namespace math { - namespace differentiation { - using finite_difference_derivative; - using complex_step_derivative; - } - } -} #endif diff --git a/test/compile_test/numerical_differentiation_concept_test.cpp b/test/compile_test/numerical_differentiation_concept_test.cpp index edf1dfa3d..7fcb898b8 100644 --- a/test/compile_test/numerical_differentiation_concept_test.cpp +++ b/test/compile_test/numerical_differentiation_concept_test.cpp @@ -4,12 +4,12 @@ // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) // #include -#include +#include void compile_and_link_test() { boost::math::concepts::std_real_concept x = 0; auto f = [](boost::math::concepts::std_real_concept x) { return x; }; - boost::math::tools::finite_difference_derivative(f, x); + boost::math::differentiation::finite_difference_derivative(f, x); } diff --git a/test/compile_test/numerical_differentiation_incl_test.cpp b/test/compile_test/numerical_differentiation_incl_test.cpp index 56d2ea607..ab7c32cc6 100644 --- a/test/compile_test/numerical_differentiation_incl_test.cpp +++ b/test/compile_test/numerical_differentiation_incl_test.cpp @@ -6,7 +6,7 @@ // Basic sanity check that header // #includes all the files that it needs to. // -#include +#include // // Note this header includes no other headers, this is // important if this test is to be meaningful: @@ -17,8 +17,8 @@ void compile_and_link_test() { auto f = [](double x) { return x; }; double x = 0; - check_result(boost::math::tools::finite_difference_derivative(f, x)); + check_result(boost::math::differentiation::finite_difference_derivative(f, x)); auto g = [](std::complex x){ return x;}; - check_result(boost::math::tools::complex_step_derivative(g, x)); + check_result(boost::math::differentiation::complex_step_derivative(g, x)); } From bf7b29f13c3fb7b90b258941d9213ac50d2e267c Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 29 Dec 2018 18:41:38 -0700 Subject: [PATCH 03/10] Fix typo. --- test/compile_test/numerical_differentiation_incl_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/compile_test/numerical_differentiation_incl_test.cpp b/test/compile_test/numerical_differentiation_incl_test.cpp index ab7c32cc6..7286949a1 100644 --- a/test/compile_test/numerical_differentiation_incl_test.cpp +++ b/test/compile_test/numerical_differentiation_incl_test.cpp @@ -6,7 +6,7 @@ // Basic sanity check that header // #includes all the files that it needs to. // -#include +#include // // Note this header includes no other headers, this is // important if this test is to be meaningful: From 7bbf05e8ba3dafc9a433347c8fef59af4093eb37 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Mon, 31 Dec 2018 20:11:25 -0700 Subject: [PATCH 04/10] 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(); +} From 1cc2ec907daa8e2d9af51ddb72b239fa92a3acfa Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Tue, 1 Jan 2019 21:14:50 -0700 Subject: [PATCH 05/10] Add example of differentiating the LIGO data [CI SKIP] --- doc/differentiation/lanczos_smoothing.qbk | 12 +- doc/graphs/ligo_derivative.svg | 2146 +++++++++++++++++++++ 2 files changed, 2157 insertions(+), 1 deletion(-) create mode 100644 doc/graphs/ligo_derivative.svg diff --git a/doc/differentiation/lanczos_smoothing.qbk b/doc/differentiation/lanczos_smoothing.qbk index 213fc964e..38162d508 100644 --- a/doc/differentiation/lanczos_smoothing.qbk +++ b/doc/differentiation/lanczos_smoothing.qbk @@ -17,7 +17,7 @@ namespace boost { namespace math { namespace differentiation { template class lanczos_derivative { public: - lanczos_smoothing_derivative(RandomAccessContainer const & v, + lanczos_derivative(RandomAccessContainer const & v, typename RandomAccessContainer::value_type spacing = 1, size_t n = 18, size_t approximation_order = 3); @@ -96,6 +96,16 @@ This is far from negligible, and as such we provide API calls which allow the fi 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 Example] + +We have extracted some data from the [@https://www.gw-openscience.org/data/ LIGO signal] and differentiated it +using the (/n/, /p/) = (60, 4) Lanczos smoothing derivative, as well as using the (/n/, /p/) = (4, 8) (nonsmoothing) derivative. + +[graph ligo_derivative] + +The original data is in orange, the smoothing derivative in blue, and the non-smoothing standard finite difference formula is in gray. +(Each time series has been rescaled to fit in the same graph.) +We can see that the smoothing derivative tracking the increase and decrease in the trend well, whereas the standard finite difference formula produces nonsense. [heading References] diff --git a/doc/graphs/ligo_derivative.svg b/doc/graphs/ligo_derivative.svg new file mode 100644 index 000000000..fd1c68201 --- /dev/null +++ b/doc/graphs/ligo_derivative.svg @@ -0,0 +1,2146 @@ + + + + + + + + +-1.843e-19 + +-1.238e-19 + +-6.319e-20 + +-2.614e-21 + +5.796e-20 + +1.185e-19 + +1.791e-19 + +2.397e-19 + +0.01707 + +0.03413 + +0.0512 + +0.06826 + +0.08533 + +0.1024 + +0.1195 + +0.1365 + +0.1536 + +0.1707 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From b2c0f9eac27f75dcb53f69a45c931da4270105dd Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Tue, 1 Jan 2019 21:31:56 -0700 Subject: [PATCH 06/10] Remove grammar errors and reduce point radius. [CI SKIP] --- doc/differentiation/lanczos_smoothing.qbk | 2 +- doc/graphs/ligo_derivative.svg | 4200 ++++++++++----------- 2 files changed, 2101 insertions(+), 2101 deletions(-) diff --git a/doc/differentiation/lanczos_smoothing.qbk b/doc/differentiation/lanczos_smoothing.qbk index 38162d508..af3437676 100644 --- a/doc/differentiation/lanczos_smoothing.qbk +++ b/doc/differentiation/lanczos_smoothing.qbk @@ -105,7 +105,7 @@ using the (/n/, /p/) = (60, 4) Lanczos smoothing derivative, as well as using th The original data is in orange, the smoothing derivative in blue, and the non-smoothing standard finite difference formula is in gray. (Each time series has been rescaled to fit in the same graph.) -We can see that the smoothing derivative tracking the increase and decrease in the trend well, whereas the standard finite difference formula produces nonsense. +We can see that the smoothing derivative tracks the increase and decrease in the trend well, whereas the standard finite difference formula produces nonsense and amplifies noise. [heading References] diff --git a/doc/graphs/ligo_derivative.svg b/doc/graphs/ligo_derivative.svg index fd1c68201..0fd967eca 100644 --- a/doc/graphs/ligo_derivative.svg +++ b/doc/graphs/ligo_derivative.svg @@ -42,2105 +42,2105 @@ 0.1536 0.1707 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From e070ed17e79e850e4100c9e8dcf65520662c1e4d Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Wed, 2 Jan 2019 12:38:58 -0700 Subject: [PATCH 07/10] Remove sign-compare warnings. Take advice of cppcheck. Grammar in documentation [CI SKIP] --- doc/differentiation/lanczos_smoothing.qbk | 36 +++++++------- .../differentiation/lanczos_smoothing.hpp | 37 ++++++++------- test/lanczos_smoothing_test.cpp | 47 +++++++++---------- 3 files changed, 60 insertions(+), 60 deletions(-) diff --git a/doc/differentiation/lanczos_smoothing.qbk b/doc/differentiation/lanczos_smoothing.qbk index af3437676..95218e82b 100644 --- a/doc/differentiation/lanczos_smoothing.qbk +++ b/doc/differentiation/lanczos_smoothing.qbk @@ -15,37 +15,39 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) namespace boost { namespace math { namespace differentiation { template - class lanczos_derivative { + class discrete_lanczos_derivative { public: - lanczos_derivative(RandomAccessContainer const & v, - typename RandomAccessContainer::value_type spacing = 1, - size_t n = 18, - size_t approximation_order = 3); + discrete_lanczos_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/. +The `discrete_lanczos_derivative` class 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); + using boost::math::differentiation::discrete_lanczos_derivative; + auto lanczos = discrete_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/. +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/: @@ -54,11 +56,12 @@ You can play around with the approximation order /p/ and the filter length /n/: auto lanczos = lanczos_derivative(v, spacing, n, p); double dvdt = lanczos[24]; -If /p=2n/, then the Lanczos smoothing derivative is not smoothing: +If /p=2n/, then the discrete Lanczos 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. +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. @@ -68,15 +71,16 @@ 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 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. +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 +the complexity of the constructor call is O(/n/[super 2]/p/). 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: diff --git a/include/boost/math/differentiation/lanczos_smoothing.hpp b/include/boost/math/differentiation/lanczos_smoothing.hpp index 5696fc454..11f346c41 100644 --- a/include/boost/math/differentiation/lanczos_smoothing.hpp +++ b/include/boost/math/differentiation/lanczos_smoothing.hpp @@ -16,7 +16,7 @@ namespace detail { template class discrete_legendre { public: - discrete_legendre(int n) : m_n{n} + explicit discrete_legendre(int n) : m_n{n} { // The integer n indexes a family of discrete Legendre polynomials indexed by k <= 2*n } @@ -128,7 +128,7 @@ class discrete_legendre { }; template -std::vector interior_filter(int n, int p) { +std::vector interior_filter(size_t n, size_t 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); @@ -136,7 +136,7 @@ std::vector interior_filter(int n, int p) { 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) + for (size_t l = 3; l < p + 1; l += 2) { dlp.next_prime(); coeffs[l] = dlp.next_prime()/ dlp.norm_sq(l); @@ -147,7 +147,7 @@ std::vector interior_filter(int n, int p) { Real arg = Real(j) / Real(n); dlp.initialize_recursion(arg); f[j] = coeffs[1]*arg; - for (int l = 3; l <= p; l += 2) + for (size_t l = 3; l <= p; l += 2) { dlp.next(); f[j] += coeffs[l]*dlp.next(); @@ -158,7 +158,8 @@ std::vector interior_filter(int n, int p) { } template -std::vector boundary_filter(int n, int p, int s) { +std::vector boundary_filter(size_t n, size_t p, int64_t s) +{ std::vector f(2 * n + 1, 0); auto dlp = discrete_legendre(n); Real sn = Real(s) / Real(n); @@ -166,7 +167,7 @@ std::vector boundary_filter(int n, int p, int s) { dlp.initialize_recursion(sn); coeffs[0] = 0; coeffs[1] = 1/dlp.norm_sq(1); - for (int l = 2; l < p + 1; ++l) + for (size_t l = 2; l < p + 1; ++l) { // Calculation of the norms is common to all filters, // so it seems like an obvious optimization target. @@ -176,14 +177,14 @@ std::vector boundary_filter(int n, int p, int s) { coeffs[l] = dlp.next_prime()/ dlp.norm_sq(l); } - for (int k = 0; k < f.size(); ++k) + for (size_t k = 0; k < f.size(); ++k) { - int j = k - n; + Real j = Real(k) - Real(n); f[k] = 0; - Real arg = Real(j) / Real(n); + Real arg = j/Real(n); dlp.initialize_recursion(arg); f[k] = coeffs[1]*arg; - for (int l = 2; l <= p; ++l) + for (size_t l = 2; l <= p; ++l) { f[k] += coeffs[l]*dlp.next(); } @@ -195,13 +196,13 @@ std::vector boundary_filter(int n, int p, int s) { } // namespace detail template -class lanczos_derivative { +class discrete_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) + discrete_lanczos_derivative(RandomAccessContainer const &v, + Real const & spacing = 1, + size_t filter_length = 18, + size_t approximation_order = 3) : m_v{v}, dt{spacing} { BOOST_ASSERT_MSG(approximation_order <= 2 * filter_length, @@ -218,8 +219,9 @@ public: for (size_t i = 0; i < filter_length; ++i) { // s = i - n; + int64_t s = static_cast(i) - static_cast(filter_length); boundary_filters[i] = detail::boundary_filter( - filter_length, approximation_order, i - filter_length); + filter_length, approximation_order, s); } } @@ -230,7 +232,7 @@ public: m_v = v; } - void reset_spacing(Real spacing) + void reset_spacing(Real const & spacing) { BOOST_ASSERT_MSG(spacing > 0, "Spacing between samples must be > 0."); dt = spacing; @@ -289,7 +291,6 @@ private: 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 diff --git a/test/lanczos_smoothing_test.cpp b/test/lanczos_smoothing_test.cpp index fd7315e88..aef9d94b2 100644 --- a/test/lanczos_smoothing_test.cpp +++ b/test/lanczos_smoothing_test.cpp @@ -18,7 +18,7 @@ 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::discrete_lanczos_derivative; using boost::math::differentiation::detail::discrete_legendre; using boost::math::differentiation::detail::interior_filter; using boost::math::differentiation::detail::boundary_filter; @@ -156,7 +156,7 @@ void test_interior_filter() // 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) + for (size_t j = 0; j < f.size(); ++j) { sum += j*f[j]; } @@ -165,7 +165,7 @@ void test_interior_filter() for (int l = 3; l <= p; l += 2) { sum = 0; - for (int j = 0; j < f.size(); ++j) + for (size_t j = 0; j < f.size(); ++j) { // The condition number of this sum is infinite! // No need to get to worked up about the tolerance. @@ -189,15 +189,13 @@ 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; - } + std::fill(v.begin(), v.end(), 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); + auto lsd = discrete_lanczos_derivative(v, 0.1, n, p); for (size_t m = n; m < v.size() - n; ++m) { Real dvdt = lsd[m]; @@ -216,7 +214,7 @@ void test_interior_lanczos() { for (size_t p = 2; p < 2*n; p += 2) { - auto lsd = lanczos_derivative(v, Real(1), n, p); + auto lsd = discrete_lanczos_derivative(v, Real(1), n, p); for (size_t m = n; m < v.size() - n; ++m) { Real dvdt = lsd[m]; @@ -240,7 +238,7 @@ void test_interior_lanczos() { for (size_t p = 2; p < 2*n; p += 2) { - auto lsd = lanczos_derivative(v, Real(1), n, p); + auto lsd = discrete_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)); @@ -258,7 +256,7 @@ void test_interior_lanczos() { for (size_t p = 2; p < 2*n; p += 2) { - auto lsd = lanczos_derivative(v, Real(1), n, p); + auto lsd = discrete_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)); @@ -277,7 +275,7 @@ void test_interior_lanczos() { for (size_t p = 3; p < 100 && p < n/2; p += 2) { - auto lsd = lanczos_derivative(v, Real(1), n, p); + auto lsd = discrete_lanczos_derivative(v, Real(1), n, p); for (size_t m = n; m < v.size() - n && m < n + 10; ++m) { @@ -313,9 +311,9 @@ void test_boundary_filters() sum = 0; c = 0; - for (int k = 0; k < f.size(); ++k) + for (size_t k = 0; k < f.size(); ++k) { - int j = k - n; + Real j = Real(k) - Real(n); // note the shifted index here: Real x = (j-s)*f[k]; Real y = x - c; @@ -330,12 +328,12 @@ void test_boundary_filters() { sum = 0; c = 0; - for (int k = 0; k < f.size(); ++k) + for (size_t k = 0; k < f.size(); ++k) { - int j = k - n; + Real j = Real(k) - Real(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 x = pow(j-s, l)*f[k]; Real y = x - c; Real t = sum + y; c = (t-sum) -y; @@ -360,16 +358,13 @@ 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; - } + std::vector v(500, 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); + auto lsd = discrete_lanczos_derivative(v, 0.0125, n, p); for (size_t m = 0; m < n; ++m) { Real dvdt = lsd[m]; @@ -392,7 +387,7 @@ void test_boundary_lanczos() { for (size_t p = 2; p < 2*n; ++p) { - auto lsd = lanczos_derivative(v, Real(1), n, p); + auto lsd = discrete_lanczos_derivative(v, Real(1), n, p); for (size_t m = 0; m < n; ++m) { Real dvdt = lsd[m]; @@ -416,7 +411,7 @@ void test_boundary_lanczos() { for (size_t p = 2; p < 2*n; ++p) { - auto lsd = lanczos_derivative(v, Real(1), n, p); + auto lsd = discrete_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)); @@ -439,7 +434,7 @@ void test_boundary_lanczos() { for (size_t p = 2; p < n; ++p) { - auto lsd = lanczos_derivative(v, Real(1), n, p); + auto lsd = discrete_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); From d76d49533ab0d917579a471f1d542c02612ed421 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Wed, 2 Jan 2019 12:53:04 -0700 Subject: [PATCH 08/10] Cleanup [CI SKIP] --- doc/differentiation/lanczos_smoothing.qbk | 13 ++++---- .../differentiation/lanczos_smoothing.hpp | 30 +++++++++---------- 2 files changed, 22 insertions(+), 21 deletions(-) diff --git a/doc/differentiation/lanczos_smoothing.qbk b/doc/differentiation/lanczos_smoothing.qbk index 95218e82b..f75be257d 100644 --- a/doc/differentiation/lanczos_smoothing.qbk +++ b/doc/differentiation/lanczos_smoothing.qbk @@ -12,24 +12,25 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) `` #include -namespace boost { namespace math { namespace differentiation { +namespace boost::math::differentiation { template class discrete_lanczos_derivative { public: + using Real = typename RandomAccessContainer::value_type; discrete_lanczos_derivative(RandomAccessContainer const & v, - typename RandomAccessContainer::value_type spacing = 1, - size_t n = 18, - size_t approximation_order = 3); + Real spacing = 1, + size_t n = 18, + size_t approximation_order = 3); - typename RandomAccessContainer::value_type operator[](size_t i) const; + Real operator[](size_t i) const; void reset_data(RandomAccessContainer const &v); void reset_spacing(Real spacing); }; -}}} // namespaces +} // namespaces `` [heading Description] diff --git a/include/boost/math/differentiation/lanczos_smoothing.hpp b/include/boost/math/differentiation/lanczos_smoothing.hpp index 11f346c41..1a57f2020 100644 --- a/include/boost/math/differentiation/lanczos_smoothing.hpp +++ b/include/boost/math/differentiation/lanczos_smoothing.hpp @@ -8,15 +8,18 @@ #include #include -namespace boost { -namespace math { -namespace differentiation { +namespace boost::math::differentiation { namespace detail { template class discrete_legendre { public: - explicit discrete_legendre(int n) : m_n{n} + explicit discrete_legendre(size_t n) : m_n{n}, m_r{2}, + m_x{std::numeric_limits::quiet_NaN()}, + m_qrm2{std::numeric_limits::quiet_NaN()}, + m_qrm1{std::numeric_limits::quiet_NaN()}, + m_qrm2p{std::numeric_limits::quiet_NaN()}, + m_qrm1p{std::numeric_limits::quiet_NaN()} { // The integer n indexes a family of discrete Legendre polynomials indexed by k <= 2*n } @@ -68,7 +71,7 @@ class discrete_legendre { } - Real operator()(Real x, int k) + Real operator()(Real x, size_t k) { BOOST_ASSERT_MSG(k <= 2 * m_n, "r <= 2n is required."); if (k == 0) @@ -82,7 +85,7 @@ class discrete_legendre { Real qrm2 = 1; Real qrm1 = x; Real N = 2 * m_n + 1; - for (int r = 2; r <= k; ++r) { + for (size_t 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; @@ -91,7 +94,7 @@ class discrete_legendre { return qrm1; } - Real prime(Real x, int k) { + Real prime(Real x, size_t k) { BOOST_ASSERT_MSG(k <= 2 * m_n, "r <= 2n is required."); if (k == 0) { return 0; @@ -104,7 +107,7 @@ class discrete_legendre { Real qrm2p = 0; Real qrm1p = 1; Real N = 2 * m_n + 1; - for (int r = 2; r <= k; ++r) { + for (size_t 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; @@ -118,13 +121,13 @@ class discrete_legendre { } private: - int m_n; + size_t m_n; + size_t m_r; + Real m_x; Real m_qrm2; Real m_qrm1; Real m_qrm2p; Real m_qrm1p; - int m_r; - Real m_x; }; template @@ -291,8 +294,5 @@ private: Real dt; }; - -} // namespace differentiation -} // namespace math -} // namespace boost +} // namespaces #endif From 95f993c9bcffa05cd1e7ef741589b3b525c92b33 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Thu, 3 Jan 2019 11:55:29 -0700 Subject: [PATCH 09/10] Add denoising second derivative. --- doc/differentiation/lanczos_smoothing.qbk | 11 +- .../differentiation/lanczos_smoothing.hpp | 207 ++++++++++++++---- test/lanczos_smoothing_test.cpp | 89 ++++++++ 3 files changed, 263 insertions(+), 44 deletions(-) diff --git a/doc/differentiation/lanczos_smoothing.qbk b/doc/differentiation/lanczos_smoothing.qbk index f75be257d..169836f9f 100644 --- a/doc/differentiation/lanczos_smoothing.qbk +++ b/doc/differentiation/lanczos_smoothing.qbk @@ -14,7 +14,7 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) namespace boost::math::differentiation { - template + template class discrete_lanczos_derivative { public: using Real = typename RandomAccessContainer::value_type; @@ -45,6 +45,15 @@ A basic usage is auto lanczos = discrete_lanczos_derivative(v, spacing); double dvdt = lanczos[30]; +Noise-suppressing second derivatives can also be computed. +The syntax is as follows: + + std::vector v(500); + // fill v with noisy data. + auto lanczos = lanczos_derivative(v, spacing); + // evaluate: + double dvdt = lanczos[25]; + 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, diff --git a/include/boost/math/differentiation/lanczos_smoothing.hpp b/include/boost/math/differentiation/lanczos_smoothing.hpp index 1a57f2020..c3b78f0a4 100644 --- a/include/boost/math/differentiation/lanczos_smoothing.hpp +++ b/include/boost/math/differentiation/lanczos_smoothing.hpp @@ -19,7 +19,9 @@ class discrete_legendre { m_qrm2{std::numeric_limits::quiet_NaN()}, m_qrm1{std::numeric_limits::quiet_NaN()}, m_qrm2p{std::numeric_limits::quiet_NaN()}, - m_qrm1p{std::numeric_limits::quiet_NaN()} + m_qrm1p{std::numeric_limits::quiet_NaN()}, + m_qrm2pp{std::numeric_limits::quiet_NaN()}, + m_qrm1pp{std::numeric_limits::quiet_NaN()} { // The integer n indexes a family of discrete Legendre polynomials indexed by k <= 2*n } @@ -35,10 +37,16 @@ class discrete_legendre { void initialize_recursion(Real x) { + using std::abs; + BOOST_ASSERT_MSG(abs(x) <= 1, "Three term recurrence is stable only for |x| <=1."); m_qrm2 = 1; m_qrm1 = x; + // Derivatives: m_qrm2p = 0; m_qrm1p = 1; + // Second derivatives: + m_qrm2pp = 0; + m_qrm1pp = 0; m_r = 2; m_x = x; @@ -58,8 +66,7 @@ class discrete_legendre { 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 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; @@ -70,6 +77,23 @@ class discrete_legendre { return m_qrm1p; } + Real next_dbl_prime() + { + Real N = 2*m_n + 1; + Real trm1 = 2*m_r - 1; + Real s = (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) / Real(4 * m_n * m_n); + Real rqrpp = 2*trm1*m_qrm1p + trm1*m_x*m_qrm1pp - s*m_qrm2pp; + 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_qrm2pp = m_qrm1pp; + m_qrm1pp = rqrpp/m_r; + ++m_r; + return m_qrm1pp; + } Real operator()(Real x, size_t k) { @@ -128,6 +152,8 @@ class discrete_legendre { Real m_qrm1; Real m_qrm2p; Real m_qrm1p; + Real m_qrm2pp; + Real m_qrm1pp; }; template @@ -196,35 +222,91 @@ std::vector boundary_filter(size_t n, size_t p, int64_t s) return f; } +template +std::vector acceleration_boundary_filter(size_t n, size_t p, int64_t s) +{ + BOOST_ASSERT_MSG(p <= 2*n, "Approximation order must be <= 2*n"); + BOOST_ASSERT_MSG(p > 2, "Approximation order must be > 2"); + std::vector f(2 * n + 1, 0); + auto dlp = discrete_legendre(n); + Real sn = Real(s) / Real(n); + std::vector coeffs(p+2, std::numeric_limits::quiet_NaN()); + dlp.initialize_recursion(sn); + coeffs[0] = 0; + coeffs[1] = 0; + for (size_t l = 2; l < p + 2; ++l) + { + coeffs[l] = dlp.next_dbl_prime()/ dlp.norm_sq(l); + } + + for (size_t k = 0; k < f.size(); ++k) + { + Real j = Real(k) - Real(n); + f[k] = 0; + Real arg = j/Real(n); + dlp.initialize_recursion(arg); + f[k] = coeffs[1]*arg; + for (size_t l = 2; l <= p; ++l) + { + f[k] += coeffs[l]*dlp.next(); + } + f[k] /= (n * n * n); + } + return f; +} + + } // namespace detail -template +template class discrete_lanczos_derivative { public: using Real = typename RandomAccessContainer::value_type; discrete_lanczos_derivative(RandomAccessContainer const &v, Real const & spacing = 1, - size_t filter_length = 18, + size_t n = 18, size_t 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, + BOOST_ASSERT_MSG(size(v) >= 2*n+1, "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) + if constexpr (order == 1) { - // s = i - n; - int64_t s = static_cast(i) - static_cast(filter_length); - boundary_filters[i] = detail::boundary_filter( - filter_length, approximation_order, s); + BOOST_ASSERT_MSG(approximation_order <= 2 * n, + "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); + + boundary_filters.resize(n); + for (size_t i = 0; i < n; ++i) + { + // s = i - n; + int64_t s = static_cast(i) - static_cast(n); + boundary_filters[i] = detail::boundary_filter(n, approximation_order, s); + } + } + else if constexpr (order == 2) + { + 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) + { + m_f[i] = f[i+n]; + } + boundary_filters.resize(n); + for (size_t i = 0; i < n; ++i) + { + int64_t s = static_cast(i) - static_cast(n); + boundary_filters[i] = detail::acceleration_boundary_filter(n, approximation_order, s); + } + } + else + { + BOOST_ASSERT_MSG(false, "Derivatives of order 3 and higher are not implemented."); } } @@ -248,38 +330,77 @@ public: Real operator[](size_t i) const { - using std::size; - if (i >= m_f.size() - 1 && i <= size(m_v) - m_f.size()) + if constexpr (order==1) { - Real dv = 0; - for (size_t j = 1; j < m_f.size(); ++j) + using std::size; + if (i >= m_f.size() - 1 && i <= size(m_v) - m_f.size()) { - dv += m_f[j] * (m_v[i + j] - m_v[i - j]); + 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; } - 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) + // m_f.size() = N+1 + if (i < m_f.size() - 1) { - dv += f[j] * m_v[m_v.size() - 1 - j]; + 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; + } + } + else if constexpr (order==2) + { + using std::size; + if (i >= m_f.size() - 1 && i <= size(m_v) - m_f.size()) + { + Real d2v = m_f[0]*m_v[i]; + for (size_t j = 1; j < m_f.size(); ++j) + { + d2v += m_f[j] * (m_v[i + j] + m_v[i - j]); + } + return d2v / (dt*dt); + } + + // m_f.size() = N+1 + if (i < m_f.size() - 1) + { + auto &f = boundary_filters[i]; + Real d2v = 0; + for (size_t j = 0; j < f.size(); ++j) { + d2v += f[j] * m_v[j]; + } + return d2v / (dt*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 d2v = 0; + for (size_t j = 0; j < f.size(); ++j) + { + d2v += f[j] * m_v[m_v.size() - 1 - j]; + } + return d2v / dt; } - return -dv / dt; } // OOB access: diff --git a/test/lanczos_smoothing_test.cpp b/test/lanczos_smoothing_test.cpp index aef9d94b2..bebc2d562 100644 --- a/test/lanczos_smoothing_test.cpp +++ b/test/lanczos_smoothing_test.cpp @@ -143,6 +143,19 @@ void test_dlp_derivatives() BOOST_CHECK_CLOSE_FRACTION(q2p, expected, tol); } +template +void test_dlp_second_derivative() +{ + std::cout << "Testing Discrete Legendre polynomial derivatives on type " << typeid(Real).name() << "\n"; + int n = 25; + auto dlp = discrete_legendre(n); + Real x = Real(1)/Real(3); + dlp.initialize_recursion(x); + Real q2pp = dlp.next_dbl_prime(); + BOOST_TEST(q2pp == 3); +} + + template void test_interior_filter() { @@ -444,8 +457,82 @@ void test_boundary_lanczos() } +template +void test_acceleration_filters() +{ + Real eps = std::numeric_limits::epsilon(); + for (size_t n = 1; n < 8; ++n) + { + for(size_t p = 3; p <= 2*n; ++p) + { + for(int64_t s = -int64_t(n); s <= 0 /*int64_t(n)*/; ++s) + { + auto f = boost::math::differentiation::detail::acceleration_boundary_filter(n,p,s); + Real sum = 0; + for (auto & x : f) + { + sum += x; + } + BOOST_CHECK_SMALL(abs(sum), sqrt(eps)); + + sum = 0; + for (size_t k = 0; k < f.size(); ++k) + { + Real j = Real(k) - Real(n); + sum += (j-s)*f[k]; + } + BOOST_CHECK_SMALL(sum, sqrt(eps)); + + sum = 0; + for (size_t k = 0; k < f.size(); ++k) + { + Real j = Real(k) - Real(n); + sum += (j-s)*(j-s)*f[k]; + } + BOOST_CHECK_CLOSE_FRACTION(sum, 2, 4*sqrt(eps)); + // More moments vanish, but the moment sums become increasingly poorly conditioned. + // We can check them later if we wish. + } + } + } +} + +template +void test_lanczos_acceleration() +{ + Real eps = std::numeric_limits::epsilon(); + std::vector v(100, 7); + auto lanczos = discrete_lanczos_derivative(v, Real(1), 4, 3); + for (size_t i = 0; i < v.size(); ++i) + { + BOOST_CHECK_SMALL(lanczos[i], eps); + } + + for(size_t i = 0; i < v.size(); ++i) + { + v[i] = 7*i + 6; + } + for (size_t i = 0; i < v.size(); ++i) + { + BOOST_CHECK_SMALL(lanczos[i], 200*eps); + } + + for(size_t i = 0; i < v.size(); ++i) + { + v[i] = 7*i*i + 9*i + 6; + } + for (size_t i = 0; i < v.size(); ++i) + { + BOOST_CHECK_CLOSE_FRACTION(lanczos[i], 14, 1000*eps); + } + + +} + BOOST_AUTO_TEST_CASE(lanczos_smoothing_test) { + test_acceleration_filters(); + test_dlp_second_derivative(); test_dlp_norms(); test_dlp_evaluation(); test_dlp_derivatives(); @@ -461,4 +548,6 @@ BOOST_AUTO_TEST_CASE(lanczos_smoothing_test) test_interior_filter(); test_interior_filter(); test_interior_lanczos(); + + test_lanczos_acceleration(); } From bfa7619954c317b0c37da1c9ef096e33b684a536 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Fri, 4 Jan 2019 12:50:58 -0700 Subject: [PATCH 10/10] Refactor so as to not store a reference member, make call threadsafe, compute entire vector in one go. --- doc/differentiation/lanczos_smoothing.qbk | 52 ++--- .../differentiation/lanczos_smoothing.hpp | 194 +++++++++++++----- test/lanczos_smoothing_test.cpp | 65 +++--- 3 files changed, 205 insertions(+), 106 deletions(-) diff --git a/doc/differentiation/lanczos_smoothing.qbk b/doc/differentiation/lanczos_smoothing.qbk index 169836f9f..fb3ba10eb 100644 --- a/doc/differentiation/lanczos_smoothing.qbk +++ b/doc/differentiation/lanczos_smoothing.qbk @@ -14,18 +14,18 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) namespace boost::math::differentiation { - template + template class discrete_lanczos_derivative { public: - using Real = typename RandomAccessContainer::value_type; - discrete_lanczos_derivative(RandomAccessContainer const & v, - Real spacing = 1, + discrete_lanczos_derivative(Real spacing, size_t n = 18, size_t approximation_order = 3); - Real operator[](size_t i) const; + template + Real operator()(RandomAccessContainer const & v, size_t i) const; - void reset_data(RandomAccessContainer const &v); + template + RandomAccessContainer operator()(RandomAccessContainer const & v) const; void reset_spacing(Real spacing); }; @@ -35,24 +35,29 @@ namespace boost::math::differentiation { [heading Description] -The `discrete_lanczos_derivative` class calculates a finite-difference approximation to the derivative of a noisy sequence of equally-spaced values /v/ at an index /i/. +The `discrete_lanczos_derivative` class calculates a finite-difference approximation to the derivative of a noisy sequence of equally-spaced values /v/. A basic usage is std::vector v(500); // fill v with noisy data. double spacing = 0.001; using boost::math::differentiation::discrete_lanczos_derivative; - auto lanczos = discrete_lanczos_derivative(v, spacing); - double dvdt = lanczos[30]; + auto lanczos = discrete_lanczos_derivative(spacing); + // Compute dvdt at index 30: + double dvdt30 = lanczos(v, 30); + // Compute derivative of entire vector: + std::vector dvdt = lanczos(v); Noise-suppressing second derivatives can also be computed. The syntax is as follows: std::vector v(500); // fill v with noisy data. - auto lanczos = lanczos_derivative(v, spacing); - // evaluate: - double dvdt = lanczos[25]; + auto lanczos = lanczos_derivative(spacing); + // evaluate second derivative at a point: + double d2vdt2 = lanczos(v, 18); + // evaluate second derivative of entire vector: + std::vector d2vdt2 = lanczos(v); 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], @@ -63,8 +68,8 @@ 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]; + auto lanczos = lanczos_derivative(spacing, n, p); + double dvdt = lanczos(v, i); If /p=2n/, then the discrete Lanczos derivative is not smoothing: It reduces to the standard /2n+1/-point finite-difference formula. @@ -79,7 +84,7 @@ 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. +though it may be sensible to take /p/ less than this. For a sinusoidal signal contaminated with AWGN, we ran a few tests showing that for SNR = 1, p = n/8 gave the best results, @@ -93,18 +98,17 @@ Since each filter has length /2n+1/ and there are /n/ filters, whose element eac the complexity of the constructor call is O(/n/[super 2]/p/). 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: +This is far from negligible, and as such the filters may be used with multiple data, and even shared between threads: - 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 . . . + std::vector v1(500); + std::vector v2(500); + // fill v1, v2 with noisy data. + auto lanczos = lanczos_derivative(spacing); + std::vector dv1dt = lanczos(v1); // threadsafe + std::vector dv2dt = lanczos(v2); // threadsafe // need to use a different spacing? - lanczos.reset_spacing(0.02); + lanczos.reset_spacing(0.02); // not threadsafe The implementation follows [@https://doi.org/10.1080/00207160.2012.666348 McDevitt, 2012], diff --git a/include/boost/math/differentiation/lanczos_smoothing.hpp b/include/boost/math/differentiation/lanczos_smoothing.hpp index c3b78f0a4..6c1349cdf 100644 --- a/include/boost/math/differentiation/lanczos_smoothing.hpp +++ b/include/boost/math/differentiation/lanczos_smoothing.hpp @@ -258,20 +258,16 @@ std::vector acceleration_boundary_filter(size_t n, size_t p, int64_t s) } // namespace detail -template +template class discrete_lanczos_derivative { public: - using Real = typename RandomAccessContainer::value_type; - discrete_lanczos_derivative(RandomAccessContainer const &v, - Real const & spacing = 1, - size_t n = 18, - size_t approximation_order = 3) - : m_v{v}, dt{spacing} + discrete_lanczos_derivative(Real const & spacing, + size_t n = 18, + size_t approximation_order = 3) + : m_dt{spacing} { + static_assert(!std::is_integral_v, "Spacing must be a floating point type."); BOOST_ASSERT_MSG(spacing > 0, "Spacing between samples must be > 0."); - using std::size; - BOOST_ASSERT_MSG(size(v) >= 2*n+1, - "Vector must be at least as long as the filter length"); if constexpr (order == 1) { @@ -281,12 +277,12 @@ public: "The approximation order must be >= 2"); m_f = detail::interior_filter(n, approximation_order); - boundary_filters.resize(n); + m_boundary_filters.resize(n); for (size_t i = 0; i < n; ++i) { // s = i - n; int64_t s = static_cast(i) - static_cast(n); - boundary_filters[i] = detail::boundary_filter(n, approximation_order, s); + m_boundary_filters[i] = detail::boundary_filter(n, approximation_order, s); } } else if constexpr (order == 2) @@ -297,11 +293,11 @@ public: { m_f[i] = f[i+n]; } - boundary_filters.resize(n); + m_boundary_filters.resize(n); for (size_t i = 0; i < n; ++i) { int64_t s = static_cast(i) - static_cast(n); - boundary_filters[i] = detail::acceleration_boundary_filter(n, approximation_order, s); + m_boundary_filters[i] = detail::acceleration_boundary_filter(n, approximation_order, s); } } else @@ -310,96 +306,96 @@ public: } } - 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 const & spacing) { BOOST_ASSERT_MSG(spacing > 0, "Spacing between samples must be > 0."); - dt = spacing; + m_dt = spacing; } Real spacing() const { - return dt; + return m_dt; } - Real operator[](size_t i) const + template + Real operator()(RandomAccessContainer const & v, size_t i) const { + static_assert(std::is_same_v, + "The type of the values in the vector provided does not match the type in the filters."); + using std::size; + BOOST_ASSERT_MSG(size(v) >= m_boundary_filters[0].size(), + "Vector must be at least as long as the filter length"); + if constexpr (order==1) { - using std::size; - if (i >= m_f.size() - 1 && i <= size(m_v) - m_f.size()) + if (i >= m_f.size() - 1 && i <= size(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]); + dv += m_f[j] * (v[i + j] - v[i - j]); } - return dv / dt; + return dv / m_dt; } // m_f.size() = N+1 if (i < m_f.size() - 1) { - auto &f = boundary_filters[i]; + auto &bf = m_boundary_filters[i]; Real dv = 0; - for (size_t j = 0; j < f.size(); ++j) { - dv += f[j] * m_v[j]; + for (size_t j = 0; j < bf.size(); ++j) + { + dv += bf[j] * v[j]; } - return dv / dt; + return dv / m_dt; } - if (i > size(m_v) - m_f.size() && i < size(m_v)) + if (i > size(v) - m_f.size() && i < size(v)) { - int k = size(m_v) - 1 - i; - auto &f = boundary_filters[k]; + int k = size(v) - 1 - i; + auto &bf = m_boundary_filters[k]; Real dv = 0; - for (size_t j = 0; j < f.size(); ++j) + for (size_t j = 0; j < bf.size(); ++j) { - dv += f[j] * m_v[m_v.size() - 1 - j]; + dv += bf[j] * v[size(v) - 1 - j]; } - return -dv / dt; + return -dv / m_dt; } } else if constexpr (order==2) { - using std::size; - if (i >= m_f.size() - 1 && i <= size(m_v) - m_f.size()) + if (i >= m_f.size() - 1 && i <= size(v) - m_f.size()) { - Real d2v = m_f[0]*m_v[i]; + Real d2v = m_f[0]*v[i]; for (size_t j = 1; j < m_f.size(); ++j) { - d2v += m_f[j] * (m_v[i + j] + m_v[i - j]); + d2v += m_f[j] * (v[i + j] + v[i - j]); } - return d2v / (dt*dt); + return d2v / (m_dt*m_dt); } // m_f.size() = N+1 if (i < m_f.size() - 1) { - auto &f = boundary_filters[i]; + auto &bf = m_boundary_filters[i]; Real d2v = 0; - for (size_t j = 0; j < f.size(); ++j) { - d2v += f[j] * m_v[j]; + for (size_t j = 0; j < bf.size(); ++j) + { + d2v += bf[j] * v[j]; } - return d2v / (dt*dt); + return d2v / (m_dt*m_dt); } - if (i > size(m_v) - m_f.size() && i < size(m_v)) + if (i > size(v) - m_f.size() && i < size(v)) { - int k = size(m_v) - 1 - i; - auto &f = boundary_filters[k]; + int k = size(v) - 1 - i; + auto &bf = m_boundary_filters[k]; Real d2v = 0; - for (size_t j = 0; j < f.size(); ++j) + for (size_t j = 0; j < bf.size(); ++j) { - d2v += f[j] * m_v[m_v.size() - 1 - j]; + d2v += bf[j] * v[size(v) - 1 - j]; } - return d2v / dt; + return d2v / (m_dt*m_dt); } } @@ -408,11 +404,97 @@ public: return std::numeric_limits::quiet_NaN(); } + template + RandomAccessContainer operator()(RandomAccessContainer const & v) const + { + static_assert(std::is_same_v, + "The type of the values in the vector provided does not match the type in the filters."); + using std::size; + BOOST_ASSERT_MSG(size(v) >= m_boundary_filters[0].size(), + "Vector must be at least as long as the filter length"); + + RandomAccessContainer w(size(v)); + if constexpr (order==1) + { + for (size_t i = 0; i < m_f.size() - 1; ++i) + { + auto &bf = m_boundary_filters[i]; + Real dv = 0; + for (size_t j = 0; j < bf.size(); ++j) + { + dv += bf[j] * v[j]; + } + w[i] = dv / m_dt; + } + + for(size_t i = m_f.size() - 1; i <= size(v) - m_f.size(); ++i) + { + Real dv = 0; + for (size_t j = 1; j < m_f.size(); ++j) + { + dv += m_f[j] * (v[i + j] - v[i - j]); + } + w[i] = dv / m_dt; + } + + + for(size_t i = size(v) - m_f.size() + 1; i < size(v); ++i) + { + int k = size(v) - 1 - i; + auto &f = m_boundary_filters[k]; + Real dv = 0; + for (size_t j = 0; j < f.size(); ++j) + { + dv += f[j] * v[size(v) - 1 - j]; + } + w[i] = -dv / m_dt; + } + } + else if constexpr (order==2) + { + // m_f.size() = N+1 + for (size_t i = 0; i < m_f.size() - 1; ++i) + { + auto &bf = m_boundary_filters[i]; + Real d2v = 0; + for (size_t j = 0; j < bf.size(); ++j) + { + d2v += bf[j] * v[j]; + } + w[i] = d2v / (m_dt*m_dt); + } + + for (size_t i = m_f.size() - 1; i <= size(v) - m_f.size(); ++i) + { + Real d2v = m_f[0]*v[i]; + for (size_t j = 1; j < m_f.size(); ++j) + { + d2v += m_f[j] * (v[i + j] + v[i - j]); + } + w[i] = d2v / (m_dt*m_dt); + } + + + for (size_t i = size(v) - m_f.size() + 1; i < size(v); ++i) + { + int k = size(v) - 1 - i; + auto &bf = m_boundary_filters[k]; + Real d2v = 0; + for (size_t j = 0; j < bf.size(); ++j) + { + d2v += bf[j] * v[size(v) - 1 - j]; + } + w[i] = d2v / (m_dt*m_dt); + } + } + + return w; + } + private: - const RandomAccessContainer &m_v; std::vector m_f; - std::vector> boundary_filters; - Real dt; + std::vector> m_boundary_filters; + Real m_dt; }; } // namespaces diff --git a/test/lanczos_smoothing_test.cpp b/test/lanczos_smoothing_test.cpp index bebc2d562..cf2d958af 100644 --- a/test/lanczos_smoothing_test.cpp +++ b/test/lanczos_smoothing_test.cpp @@ -208,12 +208,17 @@ void test_interior_lanczos() { for (size_t p = 2; p < 2*n; p += 2) { - auto lsd = discrete_lanczos_derivative(v, 0.1, n, p); + auto dld = discrete_lanczos_derivative(Real(0.1), n, p); for (size_t m = n; m < v.size() - n; ++m) { - Real dvdt = lsd[m]; + Real dvdt = dld(v, m); BOOST_CHECK_SMALL(dvdt, tol); } + auto dvdt = dld(v); + for (size_t m = n; m < v.size() - n; ++m) + { + BOOST_CHECK_SMALL(dvdt[m], tol); + } } } @@ -227,12 +232,17 @@ void test_interior_lanczos() { for (size_t p = 2; p < 2*n; p += 2) { - auto lsd = discrete_lanczos_derivative(v, Real(1), n, p); + auto dld = discrete_lanczos_derivative(Real(1), n, p); for (size_t m = n; m < v.size() - n; ++m) { - Real dvdt = lsd[m]; + Real dvdt = dld(v, m); BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, 2000*tol); } + auto dvdt = dld(v); + for (size_t m = n; m < v.size() - n; ++m) + { + BOOST_CHECK_CLOSE_FRACTION(dvdt[m], 7, 2000*tol); + } } } @@ -241,7 +251,6 @@ void test_interior_lanczos() //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); @@ -251,10 +260,10 @@ void test_interior_lanczos() { for (size_t p = 2; p < 2*n; p += 2) { - auto lsd = discrete_lanczos_derivative(v, Real(1), n, p); + auto dld = discrete_lanczos_derivative(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)); + BOOST_CHECK_CLOSE_FRACTION(dld(v, m), Real(7), Real(0.0042)); } } } @@ -269,10 +278,10 @@ void test_interior_lanczos() { for (size_t p = 2; p < 2*n; p += 2) { - auto lsd = discrete_lanczos_derivative(v, Real(1), n, p); + auto dld = discrete_lanczos_derivative(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)); + BOOST_CHECK_CLOSE_FRACTION(dld(v,m), Real(30*m + 7), Real(0.00008)); } } } @@ -288,11 +297,11 @@ void test_interior_lanczos() { for (size_t p = 3; p < 100 && p < n/2; p += 2) { - auto lsd = discrete_lanczos_derivative(v, Real(1), n, p); + auto dld = discrete_lanczos_derivative(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)); + BOOST_CHECK_CLOSE_FRACTION(dld(v,m), omega*cos(omega*m), Real(0.03)); } } } @@ -377,15 +386,15 @@ void test_boundary_lanczos() { for (size_t p = 2; p < 2*n; ++p) { - auto lsd = discrete_lanczos_derivative(v, 0.0125, n, p); + auto lsd = discrete_lanczos_derivative(Real(0.0125), n, p); for (size_t m = 0; m < n; ++m) { - Real dvdt = lsd[m]; + Real dvdt = lsd(v,m); BOOST_CHECK_SMALL(dvdt, 4*sqrt(tol)); } for (size_t m = v.size() - n; m < v.size(); ++m) { - Real dvdt = lsd[m]; + Real dvdt = lsd(v,m); BOOST_CHECK_SMALL(dvdt, 4*sqrt(tol)); } } @@ -400,16 +409,16 @@ void test_boundary_lanczos() { for (size_t p = 2; p < 2*n; ++p) { - auto lsd = discrete_lanczos_derivative(v, Real(1), n, p); + auto lsd = discrete_lanczos_derivative(Real(1), n, p); for (size_t m = 0; m < n; ++m) { - Real dvdt = lsd[m]; + Real dvdt = lsd(v,m); BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, sqrt(tol)); } for (size_t m = v.size() - n; m < v.size(); ++m) { - Real dvdt = lsd[m]; + Real dvdt = lsd(v,m); BOOST_CHECK_CLOSE_FRACTION(dvdt, 7, 4*sqrt(tol)); } } @@ -424,10 +433,10 @@ void test_boundary_lanczos() { for (size_t p = 2; p < 2*n; ++p) { - auto lsd = discrete_lanczos_derivative(v, Real(1), n, p); + auto lsd = discrete_lanczos_derivative(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)); + BOOST_CHECK_CLOSE_FRACTION(lsd(v,m), 30*m+7, 30*sqrt(tol)); } } } @@ -447,14 +456,18 @@ void test_boundary_lanczos() { for (size_t p = 2; p < n; ++p) { - auto lsd = discrete_lanczos_derivative(v, Real(1), n, p); + auto lsd = discrete_lanczos_derivative(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_CHECK_CLOSE_FRACTION(lsd(v,m), 30*m+7, 0.005); + } + auto dvdt = lsd(v); + for (size_t m = 0; m < v.size(); ++m) + { + BOOST_CHECK_CLOSE_FRACTION(dvdt[m], 30*m+7, 0.005); } } } - } template @@ -502,10 +515,10 @@ void test_lanczos_acceleration() { Real eps = std::numeric_limits::epsilon(); std::vector v(100, 7); - auto lanczos = discrete_lanczos_derivative(v, Real(1), 4, 3); + auto lanczos = discrete_lanczos_derivative(Real(1), 4, 3); for (size_t i = 0; i < v.size(); ++i) { - BOOST_CHECK_SMALL(lanczos[i], eps); + BOOST_CHECK_SMALL(lanczos(v, i), eps); } for(size_t i = 0; i < v.size(); ++i) @@ -514,7 +527,7 @@ void test_lanczos_acceleration() } for (size_t i = 0; i < v.size(); ++i) { - BOOST_CHECK_SMALL(lanczos[i], 200*eps); + BOOST_CHECK_SMALL(lanczos(v,i), 200*eps); } for(size_t i = 0; i < v.size(); ++i) @@ -523,7 +536,7 @@ void test_lanczos_acceleration() } for (size_t i = 0; i < v.size(); ++i) { - BOOST_CHECK_CLOSE_FRACTION(lanczos[i], 14, 1000*eps); + BOOST_CHECK_CLOSE_FRACTION(lanczos(v, i), 14, 1000*eps); }