592 lines · plain
1// (C) Copyright Nick Thompson 2019.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_DIFFERENTIATION_LANCZOS_SMOOTHING_HPP7#define BOOST_MATH_DIFFERENTIATION_LANCZOS_SMOOTHING_HPP8#include <cmath> // for std::abs9#include <cstddef>10#include <limits> // to nan initialize11#include <vector>12#include <string>13#include <cstdint>14#include <stdexcept>15#include <type_traits>16#include <boost/math/tools/assert.hpp>17 18#include <boost/math/tools/is_standalone.hpp>19#ifndef BOOST_MATH_STANDALONE20#include <boost/config.hpp>21#ifdef BOOST_MATH_NO_CXX17_IF_CONSTEXPR22#error "The header <boost/math/norms.hpp> can only be used in C++17 and later."23#endif24#endif25 26namespace boost::math::differentiation {27 28namespace detail {29template <typename Real>30class discrete_legendre {31 public:32 explicit discrete_legendre(std::size_t n, Real x) : m_n{n}, m_r{2}, m_x{x},33 m_qrm2{1}, m_qrm1{x},34 m_qrm2p{0}, m_qrm1p{1},35 m_qrm2pp{0}, m_qrm1pp{0}36 {37 using std::abs;38 BOOST_MATH_ASSERT_MSG(abs(m_x) <= 1, "Three term recurrence is stable only for |x| <=1.");39 // The integer n indexes a family of discrete Legendre polynomials indexed by k <= 2*n40 }41 42 Real norm_sq(int r) const43 {44 Real prod = Real(2) / Real(2 * r + 1);45 for (int k = -r; k <= r; ++k) {46 prod *= Real(2 * m_n + 1 + k) / Real(2 * m_n);47 }48 return prod;49 }50 51 Real next()52 {53 Real N = 2 * m_n + 1;54 Real num = (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) * m_qrm2;55 Real tmp = (2 * m_r - 1) * m_x * m_qrm1 - num / Real(4 * m_n * m_n);56 m_qrm2 = m_qrm1;57 m_qrm1 = tmp / m_r;58 ++m_r;59 return m_qrm1;60 }61 62 Real next_prime()63 {64 Real N = 2 * m_n + 1;65 Real s = (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) / Real(4 * m_n * m_n);66 Real tmp1 = ((2 * m_r - 1) * m_x * m_qrm1 - s * m_qrm2) / m_r;67 Real tmp2 = ((2 * m_r - 1) * (m_qrm1 + m_x * m_qrm1p) - s * m_qrm2p) / m_r;68 m_qrm2 = m_qrm1;69 m_qrm1 = tmp1;70 m_qrm2p = m_qrm1p;71 m_qrm1p = tmp2;72 ++m_r;73 return m_qrm1p;74 }75 76 Real next_dbl_prime()77 {78 Real N = 2*m_n + 1;79 Real trm1 = 2*m_r - 1;80 Real s = (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) / Real(4 * m_n * m_n);81 Real rqrpp = 2*trm1*m_qrm1p + trm1*m_x*m_qrm1pp - s*m_qrm2pp;82 Real tmp1 = ((2 * m_r - 1) * m_x * m_qrm1 - s * m_qrm2) / m_r;83 Real tmp2 = ((2 * m_r - 1) * (m_qrm1 + m_x * m_qrm1p) - s * m_qrm2p) / m_r;84 m_qrm2 = m_qrm1;85 m_qrm1 = tmp1;86 m_qrm2p = m_qrm1p;87 m_qrm1p = tmp2;88 m_qrm2pp = m_qrm1pp;89 m_qrm1pp = rqrpp/m_r;90 ++m_r;91 return m_qrm1pp;92 }93 94 Real operator()(Real x, std::size_t k)95 {96 BOOST_MATH_ASSERT_MSG(k <= 2 * m_n, "r <= 2n is required.");97 if (k == 0)98 {99 return 1;100 }101 if (k == 1)102 {103 return x;104 }105 Real qrm2 = 1;106 Real qrm1 = x;107 Real N = 2 * m_n + 1;108 for (std::size_t r = 2; r <= k; ++r) {109 Real num = (r - 1) * (N * N - (r - 1) * (r - 1)) * qrm2;110 Real tmp = (2 * r - 1) * x * qrm1 - num / Real(4 * m_n * m_n);111 qrm2 = qrm1;112 qrm1 = tmp / r;113 }114 return qrm1;115 }116 117 Real prime(Real x, std::size_t k) {118 BOOST_MATH_ASSERT_MSG(k <= 2 * m_n, "r <= 2n is required.");119 if (k == 0) {120 return 0;121 }122 if (k == 1) {123 return 1;124 }125 Real qrm2 = 1;126 Real qrm1 = x;127 Real qrm2p = 0;128 Real qrm1p = 1;129 Real N = 2 * m_n + 1;130 for (std::size_t r = 2; r <= k; ++r) {131 Real s =132 (r - 1) * (N * N - (r - 1) * (r - 1)) / Real(4 * m_n * m_n);133 Real tmp1 = ((2 * r - 1) * x * qrm1 - s * qrm2) / r;134 Real tmp2 = ((2 * r - 1) * (qrm1 + x * qrm1p) - s * qrm2p) / r;135 qrm2 = qrm1;136 qrm1 = tmp1;137 qrm2p = qrm1p;138 qrm1p = tmp2;139 }140 return qrm1p;141 }142 143 private:144 std::size_t m_n;145 std::size_t m_r;146 Real m_x;147 Real m_qrm2;148 Real m_qrm1;149 Real m_qrm2p;150 Real m_qrm1p;151 Real m_qrm2pp;152 Real m_qrm1pp;153};154 155template <class Real>156std::vector<Real> interior_velocity_filter(std::size_t n, std::size_t p) {157 auto dlp = discrete_legendre<Real>(n, 0);158 std::vector<Real> coeffs(p+1);159 coeffs[1] = 1/dlp.norm_sq(1);160 for (std::size_t l = 3; l < p + 1; l += 2)161 {162 dlp.next_prime();163 coeffs[l] = dlp.next_prime()/ dlp.norm_sq(l);164 }165 166 // We could make the filter length n, as f[0] = 0,167 // but that'd make the indexing awkward when applying the filter.168 std::vector<Real> f(n + 1);169 // This value should never be read, but this is the correct value *if it is read*.170 // Hmm, should it be a nan then? I'm not gonna agonize.171 f[0] = 0;172 for (std::size_t j = 1; j < f.size(); ++j)173 {174 Real arg = Real(j) / Real(n);175 dlp = discrete_legendre<Real>(n, arg);176 f[j] = coeffs[1]*arg;177 for (std::size_t l = 3; l <= p; l += 2)178 {179 dlp.next();180 f[j] += coeffs[l]*dlp.next();181 }182 f[j] /= (n * n);183 }184 return f;185}186 187template <class Real>188std::vector<Real> boundary_velocity_filter(std::size_t n, std::size_t p, int64_t s)189{190 std::vector<Real> coeffs(p+1, std::numeric_limits<Real>::quiet_NaN());191 Real sn = Real(s) / Real(n);192 auto dlp = discrete_legendre<Real>(n, sn);193 coeffs[0] = 0;194 coeffs[1] = 1/dlp.norm_sq(1);195 for (std::size_t l = 2; l < p + 1; ++l)196 {197 // Calculation of the norms is common to all filters,198 // so it seems like an obvious optimization target.199 // I tried this: The spent in computing the norms time is not negligible,200 // but still a small fraction of the total compute time.201 // Hence I'm not refactoring out these norm calculations.202 coeffs[l] = dlp.next_prime()/ dlp.norm_sq(l);203 }204 205 std::vector<Real> f(2*n + 1);206 for (std::size_t k = 0; k < f.size(); ++k)207 {208 Real j = Real(k) - Real(n);209 Real arg = j/Real(n);210 dlp = discrete_legendre<Real>(n, arg);211 f[k] = coeffs[1]*arg;212 for (std::size_t l = 2; l <= p; ++l)213 {214 f[k] += coeffs[l]*dlp.next();215 }216 f[k] /= (n * n);217 }218 return f;219}220 221template <class Real>222std::vector<Real> acceleration_filter(std::size_t n, std::size_t p, int64_t s)223{224 BOOST_MATH_ASSERT_MSG(p <= 2*n, "Approximation order must be <= 2*n");225 BOOST_MATH_ASSERT_MSG(p > 2, "Approximation order must be > 2");226 227 std::vector<Real> coeffs(p+1, std::numeric_limits<Real>::quiet_NaN());228 Real sn = Real(s) / Real(n);229 auto dlp = discrete_legendre<Real>(n, sn);230 coeffs[0] = 0;231 coeffs[1] = 0;232 for (std::size_t l = 2; l < p + 1; ++l)233 {234 coeffs[l] = dlp.next_dbl_prime()/ dlp.norm_sq(l);235 }236 237 std::vector<Real> f(2*n + 1, 0);238 for (std::size_t k = 0; k < f.size(); ++k)239 {240 Real j = Real(k) - Real(n);241 Real arg = j/Real(n);242 dlp = discrete_legendre<Real>(n, arg);243 for (std::size_t l = 2; l <= p; ++l)244 {245 f[k] += coeffs[l]*dlp.next();246 }247 f[k] /= (n * n * n);248 }249 return f;250}251 252 253} // namespace detail254 255template <typename Real, std::size_t order = 1>256class discrete_lanczos_derivative {257public:258 discrete_lanczos_derivative(Real const & spacing,259 std::size_t n = 18,260 std::size_t approximation_order = 3)261 : m_dt{spacing}262 {263 static_assert(!std::is_integral_v<Real>,264 "Spacing must be a floating point type.");265 BOOST_MATH_ASSERT_MSG(spacing > 0,266 "Spacing between samples must be > 0.");267 268 if constexpr (order == 1)269 {270 BOOST_MATH_ASSERT_MSG(approximation_order <= 2 * n,271 "The approximation order must be <= 2n");272 BOOST_MATH_ASSERT_MSG(approximation_order >= 2,273 "The approximation order must be >= 2");274 275 if constexpr (std::is_same_v<Real, float> || std::is_same_v<Real, double>)276 {277 auto interior = detail::interior_velocity_filter<long double>(n, approximation_order);278 m_f.resize(interior.size());279 for (std::size_t j = 0; j < interior.size(); ++j)280 {281 m_f[j] = static_cast<Real>(interior[j])/m_dt;282 }283 }284 else285 {286 m_f = detail::interior_velocity_filter<Real>(n, approximation_order);287 for (auto & x : m_f)288 {289 x /= m_dt;290 }291 }292 293 m_boundary_filters.resize(n);294 // This for loop is a natural candidate for parallelization.295 // But does it matter? Probably not.296 for (std::size_t i = 0; i < n; ++i)297 {298 if constexpr (std::is_same_v<Real, float> || std::is_same_v<Real, double>)299 {300 int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);301 auto bf = detail::boundary_velocity_filter<long double>(n, approximation_order, s);302 m_boundary_filters[i].resize(bf.size());303 for (std::size_t j = 0; j < bf.size(); ++j)304 {305 m_boundary_filters[i][j] = static_cast<Real>(bf[j])/m_dt;306 }307 }308 else309 {310 int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);311 m_boundary_filters[i] = detail::boundary_velocity_filter<Real>(n, approximation_order, s);312 for (auto & bf : m_boundary_filters[i])313 {314 bf /= m_dt;315 }316 }317 }318 }319 else if constexpr (order == 2)320 {321 // High precision isn't warranted for small p; only for large p.322 // (The computation appears stable for large n.)323 // But given that the filters are reusable for many vectors,324 // it's better to do a high precision computation and then cast back,325 // since the resulting cost is a factor of 2, and the cost of the filters not working is hours of debugging.326 if constexpr (std::is_same_v<Real, double> || std::is_same_v<Real, float>)327 {328 auto f = detail::acceleration_filter<long double>(n, approximation_order, 0);329 m_f.resize(n+1);330 for (std::size_t i = 0; i < m_f.size(); ++i)331 {332 m_f[i] = static_cast<Real>(f[i+n])/(m_dt*m_dt);333 }334 m_boundary_filters.resize(n);335 for (std::size_t i = 0; i < n; ++i)336 {337 int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);338 auto bf = detail::acceleration_filter<long double>(n, approximation_order, s);339 m_boundary_filters[i].resize(bf.size());340 for (std::size_t j = 0; j < bf.size(); ++j)341 {342 m_boundary_filters[i][j] = static_cast<Real>(bf[j])/(m_dt*m_dt);343 }344 }345 }346 else347 {348 // Given that the purpose is denoising, for higher precision calculations,349 // the default precision should be fine.350 auto f = detail::acceleration_filter<Real>(n, approximation_order, 0);351 m_f.resize(n+1);352 for (std::size_t i = 0; i < m_f.size(); ++i)353 {354 m_f[i] = f[i+n]/(m_dt*m_dt);355 }356 m_boundary_filters.resize(n);357 for (std::size_t i = 0; i < n; ++i)358 {359 int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);360 m_boundary_filters[i] = detail::acceleration_filter<Real>(n, approximation_order, s);361 for (auto & bf : m_boundary_filters[i])362 {363 bf /= (m_dt*m_dt);364 }365 }366 }367 }368 else369 {370 BOOST_MATH_ASSERT_MSG(false, "Derivatives of order 3 and higher are not implemented.");371 }372 }373 374 Real get_spacing() const375 {376 return m_dt;377 }378 379 template<class RandomAccessContainer>380 Real operator()(RandomAccessContainer const & v, std::size_t i) const381 {382 static_assert(std::is_same_v<typename RandomAccessContainer::value_type, Real>,383 "The type of the values in the vector provided does not match the type in the filters.");384 385 BOOST_MATH_ASSERT_MSG(std::size(v) >= m_boundary_filters[0].size(),386 "Vector must be at least as long as the filter length");387 388 if constexpr (order==1)389 {390 if (i >= m_f.size() - 1 && i <= std::size(v) - m_f.size())391 {392 // The filter has length >= 1:393 Real dvdt = m_f[1] * (v[i + 1] - v[i - 1]);394 for (std::size_t j = 2; j < m_f.size(); ++j)395 {396 dvdt += m_f[j] * (v[i + j] - v[i - j]);397 }398 return dvdt;399 }400 401 // m_f.size() = N+1402 if (i < m_f.size() - 1)403 {404 auto &bf = m_boundary_filters[i];405 Real dvdt = bf[0]*v[0];406 for (std::size_t j = 1; j < bf.size(); ++j)407 {408 dvdt += bf[j] * v[j];409 }410 return dvdt;411 }412 413 if (i > std::size(v) - m_f.size() && i < std::size(v))414 {415 int k = std::size(v) - 1 - i;416 auto &bf = m_boundary_filters[k];417 Real dvdt = bf[0]*v[std::size(v)-1];418 for (std::size_t j = 1; j < bf.size(); ++j)419 {420 dvdt += bf[j] * v[std::size(v) - 1 - j];421 }422 return -dvdt;423 }424 }425 else if constexpr (order==2)426 {427 if (i >= m_f.size() - 1 && i <= std::size(v) - m_f.size())428 {429 Real d2vdt2 = m_f[0]*v[i];430 for (std::size_t j = 1; j < m_f.size(); ++j)431 {432 d2vdt2 += m_f[j] * (v[i + j] + v[i - j]);433 }434 return d2vdt2;435 }436 437 // m_f.size() = N+1438 if (i < m_f.size() - 1)439 {440 auto &bf = m_boundary_filters[i];441 Real d2vdt2 = bf[0]*v[0];442 for (std::size_t j = 1; j < bf.size(); ++j)443 {444 d2vdt2 += bf[j] * v[j];445 }446 return d2vdt2;447 }448 449 if (i > std::size(v) - m_f.size() && i < std::size(v))450 {451 int k = std::size(v) - 1 - i;452 auto &bf = m_boundary_filters[k];453 Real d2vdt2 = bf[0] * v[std::size(v) - 1];454 for (std::size_t j = 1; j < bf.size(); ++j)455 {456 d2vdt2 += bf[j] * v[std::size(v) - 1 - j];457 }458 return d2vdt2;459 }460 }461 462 // OOB access:463 std::string msg = "Out of bounds access in Lanczos derivative.";464 msg += "Input vector has length " + std::to_string(std::size(v)) + ", but user requested access at index " + std::to_string(i) + ".";465 throw std::out_of_range(msg);466 return std::numeric_limits<Real>::quiet_NaN();467 }468 469 template<class RandomAccessContainer>470 void operator()(RandomAccessContainer const & v, RandomAccessContainer & w) const471 {472 static_assert(std::is_same_v<typename RandomAccessContainer::value_type, Real>,473 "The type of the values in the vector provided does not match the type in the filters.");474 if (&w[0] == &v[0])475 {476 throw std::logic_error("This transform cannot be performed in-place.");477 }478 479 if (std::size(v) < m_boundary_filters[0].size())480 {481 std::string msg = "The input vector must be at least as long as the filter length. ";482 msg += "The input vector has length = " + std::to_string(std::size(v)) + ", the filter has length " + std::to_string(m_boundary_filters[0].size());483 throw std::length_error(msg);484 }485 486 if (std::size(w) < std::size(v))487 {488 std::string msg = "The output vector (containing the derivative) must be at least as long as the input vector.";489 msg += "The output vector has length = " + std::to_string(std::size(w)) + ", the input vector has length " + std::to_string(std::size(v));490 throw std::length_error(msg);491 }492 493 if constexpr (order==1)494 {495 for (std::size_t i = 0; i < m_f.size() - 1; ++i)496 {497 auto &bf = m_boundary_filters[i];498 Real dvdt = bf[0] * v[0];499 for (std::size_t j = 1; j < bf.size(); ++j)500 {501 dvdt += bf[j] * v[j];502 }503 w[i] = dvdt;504 }505 506 for(std::size_t i = m_f.size() - 1; i <= std::size(v) - m_f.size(); ++i)507 {508 Real dvdt = m_f[1] * (v[i + 1] - v[i - 1]);509 for (std::size_t j = 2; j < m_f.size(); ++j)510 {511 dvdt += m_f[j] *(v[i + j] - v[i - j]);512 }513 w[i] = dvdt;514 }515 516 517 for(std::size_t i = std::size(v) - m_f.size() + 1; i < std::size(v); ++i)518 {519 int k = std::size(v) - 1 - i;520 auto &f = m_boundary_filters[k];521 Real dvdt = f[0] * v[std::size(v) - 1];;522 for (std::size_t j = 1; j < f.size(); ++j)523 {524 dvdt += f[j] * v[std::size(v) - 1 - j];525 }526 w[i] = -dvdt;527 }528 }529 else if constexpr (order==2)530 {531 // m_f.size() = N+1532 for (std::size_t i = 0; i < m_f.size() - 1; ++i)533 {534 auto &bf = m_boundary_filters[i];535 Real d2vdt2 = 0;536 for (std::size_t j = 0; j < bf.size(); ++j)537 {538 d2vdt2 += bf[j] * v[j];539 }540 w[i] = d2vdt2;541 }542 543 for (std::size_t i = m_f.size() - 1; i <= std::size(v) - m_f.size(); ++i)544 {545 Real d2vdt2 = m_f[0]*v[i];546 for (std::size_t j = 1; j < m_f.size(); ++j)547 {548 d2vdt2 += m_f[j] * (v[i + j] + v[i - j]);549 }550 w[i] = d2vdt2;551 }552 553 for (std::size_t i = std::size(v) - m_f.size() + 1; i < std::size(v); ++i)554 {555 int k = std::size(v) - 1 - i;556 auto &bf = m_boundary_filters[k];557 Real d2vdt2 = bf[0] * v[std::size(v) - 1];558 for (std::size_t j = 1; j < bf.size(); ++j)559 {560 d2vdt2 += bf[j] * v[std::size(v) - 1 - j];561 }562 w[i] = d2vdt2;563 }564 }565 }566 567 template<class RandomAccessContainer>568 RandomAccessContainer operator()(RandomAccessContainer const & v) const569 {570 RandomAccessContainer w(std::size(v));571 this->operator()(v, w);572 return w;573 }574 575 576 // Don't copy; too big.577 discrete_lanczos_derivative( const discrete_lanczos_derivative & ) = delete;578 discrete_lanczos_derivative& operator=(const discrete_lanczos_derivative&) = delete;579 580 // Allow moves:581 discrete_lanczos_derivative(discrete_lanczos_derivative&&) noexcept = default;582 discrete_lanczos_derivative& operator=(discrete_lanczos_derivative&&) noexcept = default;583 584private:585 std::vector<Real> m_f;586 std::vector<std::vector<Real>> m_boundary_filters;587 Real m_dt;588};589 590} // namespaces591#endif592