diff --git a/doc/differentiation/lanczos_smoothing.qbk b/doc/differentiation/lanczos_smoothing.qbk new file mode 100644 index 000000000..7a7876e31 --- /dev/null +++ b/doc/differentiation/lanczos_smoothing.qbk @@ -0,0 +1,163 @@ +[/ +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::math::differentiation { + + template + class discrete_lanczos_derivative { + public: + discrete_lanczos_derivative(Real spacing, + size_t n = 18, + size_t approximation_order = 3); + + template + Real operator()(RandomAccessContainer const & v, size_t i) const; + + template + void operator()(RandomAccessContainer const & v, RandomAccessContainer & dvdt) const; + + template + RandomAccessContainer operator()(RandomAccessContainer const & v) const; + + Real get_spacing() const; + }; + +} // namespaces +`` + +[heading Description] + +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(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(spacing); + // evaluate second derivative at a point: + double d2vdt2 = lanczos(v, 18); + // evaluate second derivative of entire vector: + std::vector d2vdt2 = lanczos(v); + +For memory conscious programmers, you can reuse the memory space for the derivative as follows: + + std::vector v(500); + std::vector dvdt(500); + // . . . define spacing, create and instance of discrete_lanczos_derivative + + // populate dvdt, perhaps in a loop: + lanczos(v, dvdt); + + +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(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. +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 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, +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. +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/). +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 the filters may be used with multiple data, and even shared between threads: + + + 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); // not threadsafe + + +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 tracks the increase and decrease in the trend well, whereas the standard finite difference formula produces nonsense and amplifies noise. + +[heading Caveats] + +The computation of the filters is ill-conditioned for large /p/. +In order to mitigate this problem, we have computed the filters in higher precision and cast the results to the desired type. +However, this simply pushes the problem to larger /p/. +In practice, this is not a problem, as large /p/ corresponds to less powerful denoising, but keep it in mind. + +In addition, the `-ffast-math` flag has a very large effect on the speed of these functions. +In our benchmarks, they were 50% faster with the flag enabled, which is much larger than the usual performance increases we see by turning on this flag. +Hence, if the default performance is not sufficient, this flag is available, though it comes with its own problems. + +This requires C++17 `if constexpr`. + +[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/fp_utilities/condition_numbers.qbk b/doc/fp_utilities/condition_numbers.qbk new file mode 100644 index 000000000..093acf3e8 --- /dev/null +++ b/doc/fp_utilities/condition_numbers.qbk @@ -0,0 +1,118 @@ +[/ +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:cond Condition Numbers] + +[heading Synopsis] + +`` +#include + + +namespace boost::math::tools { + +template +class summation_condition_number { +public: + summation_condition_number(Real const x = 0); + + void operator+=(Real const & x); + + inline void operator-=(Real const & x); + + [[nodiscard]] Real operator()() const; + + [[nodiscard]] Real sum() const; + + [[nodiscard]] Real l1_norm() const; +}; + +template +auto evaluation_condition_number(F const & f, Real const & x); + +} // namespaces +`` + +[heading Summation Condition Number] + +Here we compute the condition number of the alternating harmonic sum: + + using boost::math::tools::summation_condition_number; + auto cond = summation_condition_number(); + float max_n = 10000000; + for (float n = 1; n < max_n; n += 2) + { + cond += 1/n; + cond -= 1/(n+1); + } + std::cout << std::setprecision(std::numeric_limits::digits10); + std::cout << "ln(2) = " << boost::math::constants::ln_two() << "\n"; + std::cout << "Sum = " << cond.sum() << "\n"; + std::cout << "Condition number = " << cond() << "\n"; + +Output: + + ln(2) = 0.693147 + Sum = 0.693137 + Condition number = 22.22282 + +We see that we have lost roughly two digits of accuracy, +consistent with the heuristic that if the condition number is 10[super /k/], +then we lose /k/ significant digits in the sum. + +Our guess it that if you're worried about whether your sum is ill-conditioned, +the /last/ thing you want is for your condition number estimate to be inaccurate as well. +Since the condition number estimate relies on computing the (perhaps ill-conditioned) sum, +we have defaulted the accumulation to use Kahan summation: + + auto cond = boost::math::tools::summation_condition_number(); // will use Kahan summation. + // ... + +Output: + + ln(2) = 0.693147 + Kahan sum = 0.693147 + Condition number = 22.2228 + +If you are interested, the L[sub 1] norm is also generated by this computation, so you may query it if you like: + + float l1 = cond.l1_norm(); + // l1 = 15.4 + +[heading Condition Number of Function Evaluation] + +The [@https://en.wikipedia.org/wiki/Condition_number condition number] of function evaluation is defined as the absolute value of /xf/'(/x/)/f/(/x/)[super -1]. +It is computed as follows: + + using boost::math::tools::evaluation_condition_number; + auto f = [](double x)->double { return std::log(x); }; + double x = 1.125; + double cond = evaluation_condition_number(f, 1.125); + // cond = 1/log(x) + +[heading Caveats] + +The condition number of function evaluation gives us a measure of how sensitive our function is to roundoff error. +Unfortunately, evaluating the condition number requires evaluating the function and its derivative, +and this calculation is itself inaccurate whenever the condition number of function evaluation is large. +Sadly, this is also the regime when you are most interested in evaluating a condition number! + +This seems to be a fundamental problem. +However, it should not be necessary to evaluate the condition number to high accuracy, +valuable insights can be obtained simply by looking at the change in condition number as the function +evolves over its domain. + + +[heading References] + +* Gautschi, Walter. ['Orthogonal polynomials: computation and approximation] Oxford University Press on Demand, 2004. + +* Higham, Nicholas J. ['The accuracy of floating point summation.] SIAM Journal on Scientific Computing 14.4 (1993): 783-799. + +* Higham, Nicholas J. ['Accuracy and stability of numerical algorithms.] Vol. 80. Siam, 2002. + +[endsect] diff --git a/doc/graphs/ligo_derivative.svg b/doc/graphs/ligo_derivative.svg new file mode 100644 index 000000000..0fd967eca --- /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 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/math.qbk b/doc/math.qbk index ae927707a..8e2d23476 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -539,6 +539,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include fp_utilities/fp_facets.qbk] [include fp_utilities/float_next.qbk] [include fp_utilities/float_comparison.qbk] +[include fp_utilities/condition_numbers.qbk] [endmathpart] [mathpart cstdfloat..Specified-width floating-point typedefs] @@ -669,6 +670,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include quadrature/naive_monte_carlo.qbk] [include differentiation/numerical_differentiation.qbk] [include differentiation/autodiff.qbk] +[include differentiation/lanczos_smoothing.qbk] [endmathpart] [include complex/complex-tr1.qbk] diff --git a/doc/overview/issues.qbk b/doc/overview/issues.qbk index 34ac6b22d..8c96b5641 100644 --- a/doc/overview/issues.qbk +++ b/doc/overview/issues.qbk @@ -65,8 +65,6 @@ as a better approximation for very large degrees of freedom? [h4 Feature Requests] -We have a request for the Lambert W function, see [@https://svn.boost.org/trac/boost/ticket/11027 #11027]. - The following table lists distributions that are found in other packages but which are not yet present here, the more frequently the distribution is found, the higher the priority for implementing it: diff --git a/doc/overview/roadmap.qbk b/doc/overview/roadmap.qbk index 6f2a9fd6b..f4fcd2b6d 100644 --- a/doc/overview/roadmap.qbk +++ b/doc/overview/roadmap.qbk @@ -7,6 +7,19 @@ All bug reports including closed ones can be viewed [@https://svn.boost.org/trac/boost/query?status=assigned&status=closed&status=new&status=reopened&component=math&col=id&col=summary&col=status&col=type&col=milestone&col=component&order=priority here] and [@https://github.com/boostorg/math/issues?utf8=%E2%9C%93&q=is%3Aissue here]. +[h4 Math-2.8.1 (Boost-1.70)] + +* Add Lanczos smoothing derivatives +* Move `numerical_differentiation.hpp` from `boost/math/tools/` to `boost/math/differentiation/finite_difference.hpp`. +* Add mean, variance, skewness, kurtosis, median, Gini coefficient, and median absolute deviation to `tools/univariate_statistics.hpp`. +* Add correlation coefficients and covariance to `tools/bivariate_statistics.hpp` +* Add absolute Gini coefficient, Hoyer sparsity, oracle SNR, and the /M/[sub 2]/M/[sub 4] SNR estimator to `tools/signal_statistics.hpp`. +* Add total variation, l0, l1, l2, and sup norms, as well as corresponding distance functions to `tools/norms.hpp`. +* Add move constructors for polynomials, support complex coefficients, add `.prime()` and `.integrate()` methods. +* Add `quadratic_roots` to `tools/roots.hpp`. +* Add support for complex-valued functions to Newton's method in `roots.hpp`. +* Add Catmull-Rom interpolator. + [h4 Math-2.8.0 (Boost-1.69)] * Add LambertW functions. diff --git a/include/boost/math/differentiation/finite_difference.hpp b/include/boost/math/differentiation/finite_difference.hpp index c7615b6b2..215f5b965 100644 --- a/include/boost/math/differentiation/finite_difference.hpp +++ b/include/boost/math/differentiation/finite_difference.hpp @@ -248,7 +248,7 @@ namespace detail { } template - Real finite_difference_derivative(const F f, Real x, Real* error, const tag&) + Real finite_difference_derivative(const F, Real, Real*, 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"); diff --git a/include/boost/math/differentiation/lanczos_smoothing.hpp b/include/boost/math/differentiation/lanczos_smoothing.hpp new file mode 100644 index 000000000..e022ca5fe --- /dev/null +++ b/include/boost/math/differentiation/lanczos_smoothing.hpp @@ -0,0 +1,580 @@ +// (C) Copyright Nick Thompson 2019. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_DIFFERENTIATION_LANCZOS_SMOOTHING_HPP +#define BOOST_MATH_DIFFERENTIATION_LANCZOS_SMOOTHING_HPP +#include // for std::abs +#include // to nan initialize +#include +#include +#include +#include + +namespace boost::math::differentiation { + +namespace detail { +template +class discrete_legendre { + public: + explicit discrete_legendre(std::size_t n, Real x) : m_n{n}, m_r{2}, m_x{x}, + m_qrm2{1}, m_qrm1{x}, + m_qrm2p{0}, m_qrm1p{1}, + m_qrm2pp{0}, m_qrm1pp{0} + { + using std::abs; + BOOST_ASSERT_MSG(abs(m_x) <= 1, "Three term recurrence is stable only for |x| <=1."); + // The integer n indexes a family of discrete Legendre polynomials indexed by k <= 2*n + } + + Real norm_sq(int r) const + { + 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; + } + + 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 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, std::size_t 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 (std::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; + qrm1 = tmp / r; + } + return qrm1; + } + + Real prime(Real x, std::size_t 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 (std::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; + Real tmp2 = ((2 * r - 1) * (qrm1 + x * qrm1p) - s * qrm2p) / r; + qrm2 = qrm1; + qrm1 = tmp1; + qrm2p = qrm1p; + qrm1p = tmp2; + } + return qrm1p; + } + + private: + std::size_t m_n; + std::size_t m_r; + Real m_x; + Real m_qrm2; + Real m_qrm1; + Real m_qrm2p; + Real m_qrm1p; + Real m_qrm2pp; + Real m_qrm1pp; +}; + +template +std::vector interior_velocity_filter(std::size_t n, std::size_t p) { + auto dlp = discrete_legendre(n, 0); + std::vector coeffs(p+1); + coeffs[1] = 1/dlp.norm_sq(1); + for (std::size_t l = 3; l < p + 1; l += 2) + { + dlp.next_prime(); + coeffs[l] = dlp.next_prime()/ dlp.norm_sq(l); + } + + // 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); + // This value should never be read, but this is the correct value *if it is read*. + // Hmm, should it be a nan then? I'm not gonna agonize. + f[0] = 0; + for (std::size_t j = 1; j < f.size(); ++j) + { + Real arg = Real(j) / Real(n); + dlp = discrete_legendre(n, arg); + f[j] = coeffs[1]*arg; + for (std::size_t l = 3; l <= p; l += 2) + { + dlp.next(); + f[j] += coeffs[l]*dlp.next(); + } + f[j] /= (n * n); + } + return f; +} + +template +std::vector boundary_velocity_filter(std::size_t n, std::size_t p, int64_t s) +{ + std::vector coeffs(p+1, std::numeric_limits::quiet_NaN()); + Real sn = Real(s) / Real(n); + auto dlp = discrete_legendre(n, sn); + coeffs[0] = 0; + coeffs[1] = 1/dlp.norm_sq(1); + for (std::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. + // 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); + } + + std::vector f(2*n + 1); + for (std::size_t k = 0; k < f.size(); ++k) + { + Real j = Real(k) - Real(n); + Real arg = j/Real(n); + dlp = discrete_legendre(n, arg); + f[k] = coeffs[1]*arg; + for (std::size_t l = 2; l <= p; ++l) + { + f[k] += coeffs[l]*dlp.next(); + } + f[k] /= (n * n); + } + return f; +} + +template +std::vector acceleration_filter(std::size_t n, std::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 coeffs(p+1, std::numeric_limits::quiet_NaN()); + Real sn = Real(s) / Real(n); + auto dlp = discrete_legendre(n, sn); + coeffs[0] = 0; + coeffs[1] = 0; + for (std::size_t l = 2; l < p + 1; ++l) + { + coeffs[l] = dlp.next_dbl_prime()/ dlp.norm_sq(l); + } + + std::vector f(2*n + 1, 0); + for (std::size_t k = 0; k < f.size(); ++k) + { + Real j = Real(k) - Real(n); + Real arg = j/Real(n); + dlp = discrete_legendre(n, arg); + for (std::size_t l = 2; l <= p; ++l) + { + f[k] += coeffs[l]*dlp.next(); + } + f[k] /= (n * n * n); + } + return f; +} + + +} // namespace detail + +template +class discrete_lanczos_derivative { +public: + discrete_lanczos_derivative(Real const & spacing, + std::size_t n = 18, + std::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."); + + if constexpr (order == 1) + { + 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"); + + if constexpr (std::is_same_v || std::is_same_v) + { + auto interior = detail::interior_velocity_filter(n, approximation_order); + m_f.resize(interior.size()); + for (std::size_t j = 0; j < interior.size(); ++j) + { + m_f[j] = static_cast(interior[j])/m_dt; + } + } + else + { + m_f = detail::interior_velocity_filter(n, approximation_order); + for (auto & x : m_f) + { + x /= m_dt; + } + } + + m_boundary_filters.resize(n); + // This for loop is a natural candidate for parallelization. + // But does it matter? Probably not. + for (std::size_t i = 0; i < n; ++i) + { + if constexpr (std::is_same_v || std::is_same_v) + { + int64_t s = static_cast(i) - static_cast(n); + auto bf = detail::boundary_velocity_filter(n, approximation_order, s); + m_boundary_filters[i].resize(bf.size()); + for (std::size_t j = 0; j < bf.size(); ++j) + { + m_boundary_filters[i][j] = static_cast(bf[j])/m_dt; + } + } + else + { + int64_t s = static_cast(i) - static_cast(n); + m_boundary_filters[i] = detail::boundary_velocity_filter(n, approximation_order, s); + for (auto & bf : m_boundary_filters[i]) + { + bf /= m_dt; + } + } + } + } + else if constexpr (order == 2) + { + // High precision isn't warranted for small p; only for large p. + // (The computation appears stable for large n.) + // But given that the filters are reusable for many vectors, + // it's better to do a high precision computation and then cast back, + // since the resulting cost is a factor of 2, and the cost of the filters not working is hours of debugging. + if constexpr (std::is_same_v || std::is_same_v) + { + auto f = detail::acceleration_filter(n, approximation_order, 0); + m_f.resize(n+1); + for (std::size_t i = 0; i < m_f.size(); ++i) + { + m_f[i] = static_cast(f[i+n])/(m_dt*m_dt); + } + m_boundary_filters.resize(n); + for (std::size_t i = 0; i < n; ++i) + { + int64_t s = static_cast(i) - static_cast(n); + auto bf = detail::acceleration_filter(n, approximation_order, s); + m_boundary_filters[i].resize(bf.size()); + for (std::size_t j = 0; j < bf.size(); ++j) + { + m_boundary_filters[i][j] = static_cast(bf[j])/(m_dt*m_dt); + } + } + } + else + { + // Given that the purpose is denoising, for higher precision calculations, + // the default precision should be fine. + auto f = detail::acceleration_filter(n, approximation_order, 0); + m_f.resize(n+1); + for (std::size_t i = 0; i < m_f.size(); ++i) + { + m_f[i] = f[i+n]/(m_dt*m_dt); + } + m_boundary_filters.resize(n); + for (std::size_t i = 0; i < n; ++i) + { + int64_t s = static_cast(i) - static_cast(n); + m_boundary_filters[i] = detail::acceleration_filter(n, approximation_order, s); + for (auto & bf : m_boundary_filters[i]) + { + bf /= (m_dt*m_dt); + } + } + } + } + else + { + BOOST_ASSERT_MSG(false, "Derivatives of order 3 and higher are not implemented."); + } + } + + Real get_spacing() const + { + return m_dt; + } + + template + Real operator()(RandomAccessContainer const & v, std::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."); + + BOOST_ASSERT_MSG(std::size(v) >= m_boundary_filters[0].size(), + "Vector must be at least as long as the filter length"); + + if constexpr (order==1) + { + if (i >= m_f.size() - 1 && i <= std::size(v) - m_f.size()) + { + // The filter has length >= 1: + Real dvdt = m_f[1] * (v[i + 1] - v[i - 1]); + for (std::size_t j = 2; j < m_f.size(); ++j) + { + dvdt += m_f[j] * (v[i + j] - v[i - j]); + } + return dvdt; + } + + // m_f.size() = N+1 + if (i < m_f.size() - 1) + { + auto &bf = m_boundary_filters[i]; + Real dvdt = bf[0]*v[0]; + for (std::size_t j = 1; j < bf.size(); ++j) + { + dvdt += bf[j] * v[j]; + } + return dvdt; + } + + if (i > std::size(v) - m_f.size() && i < std::size(v)) + { + int k = std::size(v) - 1 - i; + auto &bf = m_boundary_filters[k]; + Real dvdt = bf[0]*v[std::size(v)-1]; + for (std::size_t j = 1; j < bf.size(); ++j) + { + dvdt += bf[j] * v[std::size(v) - 1 - j]; + } + return -dvdt; + } + } + else if constexpr (order==2) + { + if (i >= m_f.size() - 1 && i <= std::size(v) - m_f.size()) + { + Real d2vdt2 = m_f[0]*v[i]; + for (std::size_t j = 1; j < m_f.size(); ++j) + { + d2vdt2 += m_f[j] * (v[i + j] + v[i - j]); + } + return d2vdt2; + } + + // m_f.size() = N+1 + if (i < m_f.size() - 1) + { + auto &bf = m_boundary_filters[i]; + Real d2vdt2 = bf[0]*v[0]; + for (std::size_t j = 1; j < bf.size(); ++j) + { + d2vdt2 += bf[j] * v[j]; + } + return d2vdt2; + } + + if (i > std::size(v) - m_f.size() && i < std::size(v)) + { + int k = std::size(v) - 1 - i; + auto &bf = m_boundary_filters[k]; + Real d2vdt2 = bf[0] * v[std::size(v) - 1]; + for (std::size_t j = 1; j < bf.size(); ++j) + { + d2vdt2 += bf[j] * v[std::size(v) - 1 - j]; + } + return d2vdt2; + } + } + + // OOB access: + std::string msg = "Out of bounds access in Lanczos derivative."; + msg += "Input vector has length " + std::to_string(std::size(v)) + ", but user requested access at index " + std::to_string(i) + "."; + throw std::out_of_range(msg); + return std::numeric_limits::quiet_NaN(); + } + + template + void operator()(RandomAccessContainer const & v, RandomAccessContainer & w) const + { + static_assert(std::is_same_v, + "The type of the values in the vector provided does not match the type in the filters."); + if (&w[0] == &v[0]) + { + throw std::logic_error("This transform cannot be performed in-place."); + } + + if (std::size(v) < m_boundary_filters[0].size()) + { + std::string msg = "The input vector must be at least as long as the filter length. "; + msg += "The input vector has length = " + std::to_string(std::size(v)) + ", the filter has length " + std::to_string(m_boundary_filters[0].size()); + throw std::length_error(msg); + } + + if (std::size(w) < std::size(v)) + { + std::string msg = "The output vector (containing the derivative) must be at least as long as the input vector."; + msg += "The output vector has length = " + std::to_string(std::size(w)) + ", the input vector has length " + std::to_string(std::size(v)); + throw std::length_error(msg); + } + + if constexpr (order==1) + { + for (std::size_t i = 0; i < m_f.size() - 1; ++i) + { + auto &bf = m_boundary_filters[i]; + Real dvdt = bf[0] * v[0]; + for (std::size_t j = 1; j < bf.size(); ++j) + { + dvdt += bf[j] * v[j]; + } + w[i] = dvdt; + } + + for(std::size_t i = m_f.size() - 1; i <= std::size(v) - m_f.size(); ++i) + { + Real dvdt = m_f[1] * (v[i + 1] - v[i - 1]); + for (std::size_t j = 2; j < m_f.size(); ++j) + { + dvdt += m_f[j] *(v[i + j] - v[i - j]); + } + w[i] = dvdt; + } + + + for(std::size_t i = std::size(v) - m_f.size() + 1; i < std::size(v); ++i) + { + int k = std::size(v) - 1 - i; + auto &f = m_boundary_filters[k]; + Real dvdt = f[0] * v[std::size(v) - 1];; + for (std::size_t j = 1; j < f.size(); ++j) + { + dvdt += f[j] * v[std::size(v) - 1 - j]; + } + w[i] = -dvdt; + } + } + else if constexpr (order==2) + { + // m_f.size() = N+1 + for (std::size_t i = 0; i < m_f.size() - 1; ++i) + { + auto &bf = m_boundary_filters[i]; + Real d2vdt2 = 0; + for (std::size_t j = 0; j < bf.size(); ++j) + { + d2vdt2 += bf[j] * v[j]; + } + w[i] = d2vdt2; + } + + for (std::size_t i = m_f.size() - 1; i <= std::size(v) - m_f.size(); ++i) + { + Real d2vdt2 = m_f[0]*v[i]; + for (std::size_t j = 1; j < m_f.size(); ++j) + { + d2vdt2 += m_f[j] * (v[i + j] + v[i - j]); + } + w[i] = d2vdt2; + } + + for (std::size_t i = std::size(v) - m_f.size() + 1; i < std::size(v); ++i) + { + int k = std::size(v) - 1 - i; + auto &bf = m_boundary_filters[k]; + Real d2vdt2 = bf[0] * v[std::size(v) - 1]; + for (std::size_t j = 1; j < bf.size(); ++j) + { + d2vdt2 += bf[j] * v[std::size(v) - 1 - j]; + } + w[i] = d2vdt2; + } + } + } + + template + RandomAccessContainer operator()(RandomAccessContainer const & v) const + { + RandomAccessContainer w(std::size(v)); + this->operator()(v, w); + return w; + } + + + // Don't copy; too big. + discrete_lanczos_derivative( const discrete_lanczos_derivative & ) = delete; + discrete_lanczos_derivative& operator=(const discrete_lanczos_derivative&) = delete; + + // Allow moves: + discrete_lanczos_derivative(discrete_lanczos_derivative&&) = default; + discrete_lanczos_derivative& operator=(discrete_lanczos_derivative&&) = default; + +private: + std::vector m_f; + std::vector> m_boundary_filters; + Real m_dt; +}; + +} // namespaces +#endif diff --git a/include/boost/math/tools/condition_numbers.hpp b/include/boost/math/tools/condition_numbers.hpp new file mode 100644 index 000000000..afbea75cb --- /dev/null +++ b/include/boost/math/tools/condition_numbers.hpp @@ -0,0 +1,139 @@ +// (C) Copyright Nick Thompson 2019. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP +#define BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP +#include +#include + +namespace boost::math::tools { + +template +class summation_condition_number { +public: + summation_condition_number(Real const x = 0) + { + using std::abs; + m_l1 = abs(x); + m_sum = x; + m_c = 0; + } + + void operator+=(Real const & x) + { + using std::abs; + // No need to Kahan the l1 calc; it's well conditioned: + m_l1 += abs(x); + if constexpr(kahan) + { + Real y = x - m_c; + Real t = m_sum + y; + m_c = (t-m_sum) -y; + m_sum = t; + } + else + { + m_sum += x; + } + } + + inline void operator-=(Real const & x) + { + this->operator+=(-x); + } + + // Is operator*= relevant? Presumably everything gets rescaled, + // (m_sum -> k*m_sum, m_l1->k*m_l1, m_c->k*m_c), + // but is this sensible? More important is it useful? + // In addition, it might change the condition number. + + [[nodiscard]] Real operator()() const + { + using std::abs; + if (m_sum == Real(0) && m_l1 != Real(0)) + { + return std::numeric_limits::infinity(); + } + return m_l1/abs(m_sum); + } + + [[nodiscard]] Real sum() const + { + // Higham, 1993, "The Accuracy of Floating Point Summation": + // "In [17] and [18], Kahan describes a variation of compensated summation in which the final sum is also corrected + // thus s=s+e is appended to the algorithm above)." + return m_sum + m_c; + } + + [[nodiscard]] Real l1_norm() const + { + return m_l1; + } + +private: + Real m_l1; + Real m_sum; + Real m_c; +}; + +template +auto evaluation_condition_number(F const & f, Real const & x) +{ + using std::abs; + using std::isnan; + using std::sqrt; + using boost::math::differentiation::finite_difference_derivative; + + Real fx = f(x); + if (isnan(fx)) + { + return std::numeric_limits::quiet_NaN(); + } + bool caught_exception = false; + Real fp; + try + { + fp = finite_difference_derivative(f, x); + } + catch(...) + { + caught_exception = true; + } + + if (isnan(fp) || caught_exception) + { + // Check if the right derivative exists: + fp = finite_difference_derivative(f, x); + if (isnan(fp)) + { + // Check if a left derivative exists: + const Real eps = (std::numeric_limits::epsilon)(); + Real h = - 2 * sqrt(eps); + h = boost::math::differentiation::detail::make_xph_representable(x, h); + Real yh = f(x + h); + Real y0 = f(x); + Real diff = yh - y0; + fp = diff / h; + if (isnan(fp)) + { + return std::numeric_limits::quiet_NaN(); + } + } + } + + if (fx == 0) + { + if (x==0 || fp==0) + { + return std::numeric_limits::quiet_NaN(); + } + return std::numeric_limits::infinity(); + } + + return abs(x*fp/fx); +} + +} +#endif diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index f3e0f6f55..c49dcccdc 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -906,6 +906,8 @@ test-suite misc : [ run norms_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run signal_statistics_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run bivariate_statistics_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] + [ run lanczos_smoothing_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] + [ run condition_number_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ 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/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index 78f6f5365..55c36ee81 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -175,7 +175,7 @@ void test_affine_invariance() std::vector> v(100); std::vector> u(100); std::mt19937_64 gen(438232); - Real inv_denom = (Real) 100/( (Real) gen.max() + (Real) 2); + Real inv_denom = (Real) 100/( (Real) (gen.max)() + (Real) 2); for(size_t j = 0; j < dimension; ++j) { v[0][j] = gen()*inv_denom; diff --git a/test/condition_number_test.cpp b/test/condition_number_test.cpp new file mode 100644 index 000000000..2adb9be27 --- /dev/null +++ b/test/condition_number_test.cpp @@ -0,0 +1,122 @@ +// (C) Copyright Nick Thompson, 2019 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_TEST_MODULE condition_number_test + +#include +#include +#include +#include +#include +#include +#include +#include + +using std::abs; +using boost::math::constants::half; +using boost::math::constants::ln_two; +using boost::multiprecision::cpp_bin_float_50; +using boost::math::tools::summation_condition_number; +using boost::math::tools::evaluation_condition_number; + +template +void test_summation_condition_number() +{ + Real tol = 1000*std::numeric_limits::epsilon(); + auto cond = summation_condition_number(); + // I've checked that the condition number increases with max_n, + // and that the computed sum gets more accurate with increasing max_n. + // But the CI system would die with more terms. + Real max_n = 10000; + for (Real n = 1; n < max_n; n += 2) + { + cond += 1/n; + cond -= 1/(n+1); + } + + BOOST_CHECK_CLOSE_FRACTION(cond.sum(), ln_two(), tol); + BOOST_TEST(cond() > 14); +} + +template +void test_exponential_sum() +{ + using std::exp; + using std::abs; + Real eps = std::numeric_limits::epsilon(); + for (Real x = -20; x <= -1; x += 0.5) + { + auto cond = summation_condition_number(1); + size_t n = 1; + Real term = x; + while(n++ < 1000) + { + cond += term; + term *= (x/n); + } + BOOST_CHECK_CLOSE_FRACTION(exp(x), cond.sum(), eps*cond()); + BOOST_CHECK_CLOSE_FRACTION(exp(2*abs(x)), cond(), eps*cond()); + } +} + + + +template +void test_evaluation_condition_number() +{ + using std::abs; + using std::log; + using std::sqrt; + using std::exp; + using std::sin; + using std::tan; + Real tol = sqrt(std::numeric_limits::epsilon()); + + auto f1 = [](auto x) { return log(x); }; + for (Real x = 1.125; x < 8; x += 0.125) + { + Real cond = evaluation_condition_number(f1, x); + BOOST_CHECK_CLOSE_FRACTION(cond, 1/log(x), tol); + } + + auto f2 = [](auto x) { return exp(x); }; + for (Real x = 1.125; x < 8; x += 0.125) + { + Real cond = evaluation_condition_number(f2, x); + BOOST_CHECK_CLOSE_FRACTION(cond, x, tol); + } + + auto f3 = [](auto x) { return sin(x); }; + for (Real x = 1.125; x < 8; x += 0.125) + { + Real cond = evaluation_condition_number(f3, x); + BOOST_CHECK_CLOSE_FRACTION(cond, abs(x/tan(x)), tol); + } + + // Test a function which right differentiable: + using boost::math::constants::e; + auto f4 = [](Real x) { return boost::math::lambert_w0(x); }; + Real cond = evaluation_condition_number(f4, -1/e()); + if (std::is_same_v) + { + BOOST_CHECK_GE(cond, 30); + } + else + { + BOOST_CHECK_GE(cond, 4900); + } +} + + +BOOST_AUTO_TEST_CASE(numerical_differentiation_test) +{ + test_summation_condition_number(); + test_summation_condition_number(); + test_evaluation_condition_number(); + test_evaluation_condition_number(); + test_evaluation_condition_number(); + test_evaluation_condition_number(); + test_exponential_sum(); +} diff --git a/test/lanczos_smoothing_test.cpp b/test/lanczos_smoothing_test.cpp new file mode 100644 index 000000000..5e5169298 --- /dev/null +++ b/test/lanczos_smoothing_test.cpp @@ -0,0 +1,708 @@ +/* + * Copyright Nick Thompson, 2019 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ +#define BOOST_TEST_MODULE lanczos_smoothing_test + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // for float_distance +#include + +using std::abs; +using std::pow; +using std::sqrt; +using std::sin; +using boost::math::constants::two_pi; +using boost::multiprecision::cpp_bin_float_50; +using boost::multiprecision::cpp_bin_float_100; +using boost::math::differentiation::discrete_lanczos_derivative; +using boost::math::differentiation::detail::discrete_legendre; +using boost::math::differentiation::detail::interior_velocity_filter; +using boost::math::differentiation::detail::boundary_velocity_filter; +using boost::math::tools::summation_condition_number; + +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, Real(0)); + BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(0), 3, tol); + BOOST_CHECK_CLOSE_FRACTION(dlp.norm_sq(1), 2, tol); + dlp = discrete_legendre(2, Real(0)); + 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, Real(0)); + 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; + Real x = 0.72; + auto dlp = discrete_legendre(n, x); + 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, x); + 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) + { + for(Real x = -1; x <= 1; x += 0.1) + { + auto dlp = discrete_legendre(n, x); + for (size_t k = 2; k < n; ++k) + { + BOOST_CHECK_CLOSE(dlp.next(), dlp(x, k), tol); + } + + dlp = discrete_legendre(n, 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; + Real x = 0.72; + auto dlp = discrete_legendre(n, x); + 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_dlp_second_derivative() +{ + std::cout << "Testing Discrete Legendre polynomial derivatives on type " << typeid(Real).name() << "\n"; + int n = 25; + Real x = Real(1)/Real(3); + auto dlp = discrete_legendre(n, x); + Real q2pp = dlp.next_dbl_prime(); + BOOST_TEST(q2pp == 3); +} + + +template +void test_interior_velocity_filter() +{ + using boost::math::constants::half; + 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_velocity_filter(n,p); + // Since we only store half the filter coefficients, + // we need to reindex the moment sums: + auto cond = summation_condition_number(0); + for (size_t j = 0; j < f.size(); ++j) + { + cond += j*f[j]; + } + BOOST_CHECK_CLOSE_FRACTION(cond.sum(), half(), 2*cond()*tol); + + for (int l = 3; l <= p; l += 2) + { + cond = summation_condition_number(0); + for (size_t j = 0; j < f.size() - 1; ++j) + { + cond += pow(Real(j), l)*f[j]; + } + Real expected = -pow(Real(f.size() - 1), l)*f[f.size()-1]; + BOOST_CHECK_CLOSE_FRACTION(expected, cond.sum(), 7*cond()*tol); + } + //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); + 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 dld = discrete_lanczos_derivative(Real(0.1), n, p); + for (size_t m = n; m < v.size() - n; ++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); + } + } + } + + + 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 dld = discrete_lanczos_derivative(Real(1), n, p); + for (size_t m = n; m < v.size() - n; ++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); + } + } + } + + //std::random_device rd{}; + //auto seed = rd(); + //std::cout << "Seed = " << seed << "\n"; + std::mt19937 gen(4172378669); + std::normal_distribution<> dis{0, 0.01}; + 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 dld = discrete_lanczos_derivative(Real(1), n, p); + for (size_t m = n; m < v.size() - n; ++m) + { + BOOST_CHECK_CLOSE_FRACTION(dld(v, 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 dld = discrete_lanczos_derivative(Real(1), n, p); + for (size_t m = n; m < v.size() - n; ++m) + { + BOOST_CHECK_CLOSE_FRACTION(dld(v,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 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(dld(v,m), omega*cos(omega*m), Real(0.03)); + } + } + } +} + +template +void test_boundary_velocity_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_velocity_filter(n, p, s); + // Sum is zero: + auto cond = summation_condition_number(0); + for (size_t i = 0; i < f.size() - 1; ++i) + { + cond += f[i]; + } + + BOOST_CHECK_CLOSE_FRACTION(cond.sum(), -f[f.size()-1], 6*cond()*tol); + + cond = summation_condition_number(0); + for (size_t k = 0; k < f.size(); ++k) + { + Real j = Real(k) - Real(n); + // note the shifted index here: + cond += (j-s)*f[k]; + } + BOOST_CHECK_CLOSE_FRACTION(cond.sum(), 1, 6*cond()*tol); + + + for (int l = 2; l <= p; ++l) + { + cond = summation_condition_number(0); + for (size_t k = 0; k < f.size() - 1; ++k) + { + Real j = Real(k) - Real(n); + // The condition number of this sum is infinite! + // No need to get to worked up about the tolerance. + cond += pow(j-s, l)*f[k]; + } + + Real expected = -pow(Real(f.size()-1) - Real(n) - Real(s), l)*f[f.size()-1]; + if (expected == 0) + { + BOOST_CHECK_SMALL(cond.sum(), cond()*tol); + } + else + { + BOOST_CHECK_CLOSE_FRACTION(expected, cond.sum(), 200*cond()*tol); + } + } + + //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, 7); + + for (size_t n = 1; n < 10; ++n) + { + for (size_t p = 2; p < 2*n; ++p) + { + auto lsd = discrete_lanczos_derivative(Real(0.0125), n, p); + for (size_t m = 0; m < n; ++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(v,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 = discrete_lanczos_derivative(Real(1), n, p); + for (size_t m = 0; m < n; ++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(v,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 = discrete_lanczos_derivative(Real(1), n, p); + for (size_t m = 0; m < v.size(); ++m) + { + BOOST_CHECK_CLOSE_FRACTION(lsd(v,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 = discrete_lanczos_derivative(Real(1), n, p); + for (size_t m = 0; m < v.size(); ++m) + { + 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 +void test_acceleration_filters() +{ + Real eps = std::numeric_limits::epsilon(); + for (size_t n = 1; n < 5; ++n) + { + for(size_t p = 3; p <= 2*n; ++p) + { + for(int64_t s = -int64_t(n); s <= 0; ++s) + { + auto g = boost::math::differentiation::detail::acceleration_filter(n,p,s); + + std::vector f(g.size()); + for (size_t i = 0; i < g.size(); ++i) + { + f[i] = static_cast(g[i]); + } + + auto cond = summation_condition_number(0); + + for (size_t i = 0; i < f.size() - 1; ++i) + { + cond += f[i]; + } + BOOST_CHECK_CLOSE_FRACTION(cond.sum(), -f[f.size()-1], cond()*eps); + + + cond = summation_condition_number(0); + for (size_t k = 0; k < f.size() -1; ++k) + { + Real j = Real(k) - Real(n); + cond += (j-s)*f[k]; + } + Real expected = -(Real(f.size()-1)- Real(n) - s)*f[f.size()-1]; + BOOST_CHECK_CLOSE_FRACTION(cond.sum(), expected, cond()*eps); + + cond = summation_condition_number(0); + for (size_t k = 0; k < f.size(); ++k) + { + Real j = Real(k) - Real(n); + cond += (j-s)*(j-s)*f[k]; + } + BOOST_CHECK_CLOSE_FRACTION(cond.sum(), 2, cond()*eps); + // See unlabelled equation in McDevitt, 2012, just after equation 26: + // It appears that there is an off-by-one error in that equation, since p + 1 moments don't vanish, only p. + // This test is itself suspect; the condition number of the moment sum is infinite. + // So the *slightest* error in the filter gets amplified by the test; in terms of the + // behavior of the actual filter, it's not a big deal. + for (size_t l = 3; l <= p; ++l) + { + cond = summation_condition_number(0); + for (size_t k = 0; k < f.size() - 1; ++k) + { + Real j = Real(k) - Real(n); + cond += pow((j-s), l)*f[k]; + } + Real expected = -pow(Real(f.size()- 1 - n -s), l)*f[f.size()-1]; + BOOST_CHECK_CLOSE_FRACTION(cond.sum(), expected, cond()*eps); + } + } + } + } +} + +template +void test_lanczos_acceleration() +{ + Real eps = std::numeric_limits::epsilon(); + std::vector v(100, 7); + auto lanczos = discrete_lanczos_derivative(Real(1), 4, 3); + for (size_t i = 0; i < v.size(); ++i) + { + BOOST_CHECK_SMALL(lanczos(v, 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(v,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(v, i), 14, 1500*eps); + } + + // Now add noise, and kick up the smoothing of the Lanzcos derivative (increase n): + //std::random_device rd{}; + //auto seed = rd(); + //std::cout << "seed = " << seed << "\n"; + size_t seed = 2507134629; + std::mt19937 gen(seed); + Real std_dev = 0.1; + std::normal_distribution dis{0, std_dev}; + for (size_t i = 0; i < v.size(); ++i) + { + v[i] += dis(gen); + } + lanczos = discrete_lanczos_derivative(Real(1), 18, 3); + auto w = lanczos(v); + for (size_t i = 0; i < v.size(); ++i) + { + BOOST_CHECK_CLOSE_FRACTION(w[i], 14, std_dev/200); + } +} + +template +void test_rescaling() +{ + std::cout << "Test rescaling on type " << typeid(Real).name() << "\n"; + Real tol = std::numeric_limits::epsilon(); + std::vector v(500); + for(size_t i = 0; i < v.size(); ++i) + { + v[i] = 7*i*i + 9*i + 6; + } + std::vector dvdt1(500); + std::vector dvdt2(500); + auto lanczos1 = discrete_lanczos_derivative(Real(1)); + auto lanczos2 = discrete_lanczos_derivative(Real(2)); + + lanczos1(v, dvdt1); + lanczos2(v, dvdt2); + + for(size_t i = 0; i < v.size(); ++i) + { + BOOST_CHECK_CLOSE_FRACTION(dvdt1[i], 2*dvdt2[i], tol); + } + + auto lanczos3 = discrete_lanczos_derivative(Real(1)); + auto lanczos4 = discrete_lanczos_derivative(Real(2)); + + + std::vector dv2dt21(500); + std::vector dv2dt22(500); + + for(size_t i = 0; i < v.size(); ++i) + { + BOOST_CHECK_CLOSE_FRACTION(dv2dt21[i], 4*dv2dt22[i], tol); + } +} + +template +void test_data_representations() +{ + std::cout << "Test rescaling on type " << typeid(Real).name() << "\n"; + Real tol = 150*std::numeric_limits::epsilon(); + std::array v; + for(size_t i = 0; i < v.size(); ++i) + { + v[i] = 9*i + 6; + } + std::array dvdt; + auto lanczos = discrete_lanczos_derivative(Real(1)); + + lanczos(v, dvdt); + + for(size_t i = 0; i < v.size(); ++i) + { + BOOST_CHECK_CLOSE_FRACTION(dvdt[i], 9, tol); + } + + boost::numeric::ublas::vector w(500); + boost::numeric::ublas::vector dwdt(500); + for(size_t i = 0; i < w.size(); ++i) + { + w[i] = 9*i + 6; + } + + lanczos(w, dwdt); + + for(size_t i = 0; i < v.size(); ++i) + { + BOOST_CHECK_CLOSE_FRACTION(dwdt[i], 9, tol); + } + + auto v1 = boost::make_iterator_range(v.begin(), v.end()); + auto v2 = boost::make_iterator_range(dvdt.begin(), dvdt.end()); + lanczos(v1, v2); + + for(size_t i = 0; i < v2.size(); ++i) + { + BOOST_CHECK_CLOSE_FRACTION(v2[i], 9, tol); + } + + auto lanczos2 = discrete_lanczos_derivative(Real(1)); + + lanczos2(v1, v2); + + for(size_t i = 0; i < v2.size(); ++i) + { + BOOST_CHECK_SMALL(v2[i], tol); + } + +} + +BOOST_AUTO_TEST_CASE(lanczos_smoothing_test) +{ + test_dlp_second_derivative(); + test_dlp_norms(); + test_dlp_evaluation(); + test_dlp_derivatives(); + test_dlp_next(); + + // Takes too long! + //test_dlp_norms(); + test_boundary_velocity_filters(); + test_boundary_velocity_filters(); + // Takes too long! + //test_boundary_velocity_filters(); + test_boundary_lanczos(); + test_boundary_lanczos(); + // Takes too long! + //test_boundary_lanczos(); + + test_interior_velocity_filter(); + test_interior_velocity_filter(); + test_interior_lanczos(); + + test_acceleration_filters(); + + test_lanczos_acceleration(); + test_lanczos_acceleration(); + + test_rescaling(); + test_data_representations(); +} diff --git a/test/univariate_statistics_test.cpp b/test/univariate_statistics_test.cpp index d39b24eba..a2db405c7 100644 --- a/test/univariate_statistics_test.cpp +++ b/test/univariate_statistics_test.cpp @@ -748,7 +748,6 @@ int main() test_median_absolute_deviation(); test_median_absolute_deviation(); test_median_absolute_deviation(); - test_median_absolute_deviation(); test_gini_coefficient(); test_gini_coefficient();