math/differentiation/lanczos_smoothing.hpp

0001 //  (C) Copyright Nick Thompson 2019.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_DIFFERENTIATION_LANCZOS_SMOOTHING_HPP
0007 #define BOOST_MATH_DIFFERENTIATION_LANCZOS_SMOOTHING_HPP
0008 #include <cmath> // for std::abs
0009 #include <cstddef>
0010 #include <limits> // to nan initialize
0011 #include <vector>
0012 #include <string>
0013 #include <cstdint>
0014 #include <stdexcept>
0015 #include <type_traits>
0016 #include <boost/math/tools/assert.hpp>
0017
0018 #include <boost/math/tools/is_standalone.hpp>
0019 #ifndef BOOST_MATH_STANDALONE
0020 #include <boost/config.hpp>
0021 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
0022 #error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
0023 #endif
0024 #endif
0025
0026 namespace boost::math::differentiation {
0027
0028 namespace detail {
0029 template <typename Real>
0030 class discrete_legendre {
0031   public:
0032     explicit discrete_legendre(std::size_t n, Real x) : m_n{n}, m_r{2}, m_x{x},
0033                                                         m_qrm2{1}, m_qrm1{x},
0034                                                         m_qrm2p{0}, m_qrm1p{1},
0035                                                         m_qrm2pp{0}, m_qrm1pp{0}
0036     {
0037         using std::abs;
0038         BOOST_MATH_ASSERT_MSG(abs(m_x) <= 1, "Three term recurrence is stable only for |x| <=1.");
0039         // The integer n indexes a family of discrete Legendre polynomials indexed by k <= 2*n
0040     }
0041
0042     Real norm_sq(int r) const
0043     {
0044         Real prod = Real(2) / Real(2 * r + 1);
0045         for (int k = -r; k <= r; ++k) {
0046             prod *= Real(2 * m_n + 1 + k) / Real(2 * m_n);
0047         }
0048         return prod;
0049     }
0050
0051     Real next()
0052     {
0053         Real N = 2 * m_n + 1;
0054         Real num = (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) * m_qrm2;
0055         Real tmp = (2 * m_r - 1) * m_x * m_qrm1 - num / Real(4 * m_n * m_n);
0056         m_qrm2 = m_qrm1;
0057         m_qrm1 = tmp / m_r;
0058         ++m_r;
0059         return m_qrm1;
0060     }
0061
0062     Real next_prime()
0063     {
0064         Real N = 2 * m_n + 1;
0065         Real s = (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) / Real(4 * m_n * m_n);
0066         Real tmp1 = ((2 * m_r - 1) * m_x * m_qrm1 - s * m_qrm2) / m_r;
0067         Real tmp2 = ((2 * m_r - 1) * (m_qrm1 + m_x * m_qrm1p) - s * m_qrm2p) / m_r;
0068         m_qrm2 = m_qrm1;
0069         m_qrm1 = tmp1;
0070         m_qrm2p = m_qrm1p;
0071         m_qrm1p = tmp2;
0072         ++m_r;
0073         return m_qrm1p;
0074     }
0075
0076     Real next_dbl_prime()
0077     {
0078         Real N = 2*m_n + 1;
0079         Real trm1 = 2*m_r - 1;
0080         Real s = (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) / Real(4 * m_n * m_n);
0081         Real rqrpp = 2*trm1*m_qrm1p + trm1*m_x*m_qrm1pp - s*m_qrm2pp;
0082         Real tmp1 = ((2 * m_r - 1) * m_x * m_qrm1 - s * m_qrm2) / m_r;
0083         Real tmp2 = ((2 * m_r - 1) * (m_qrm1 + m_x * m_qrm1p) - s * m_qrm2p) / m_r;
0084         m_qrm2 = m_qrm1;
0085         m_qrm1 = tmp1;
0086         m_qrm2p = m_qrm1p;
0087         m_qrm1p = tmp2;
0088         m_qrm2pp = m_qrm1pp;
0089         m_qrm1pp = rqrpp/m_r;
0090         ++m_r;
0091         return m_qrm1pp;
0092     }
0093
0094     Real operator()(Real x, std::size_t k)
0095     {
0096         BOOST_MATH_ASSERT_MSG(k <= 2 * m_n, "r <= 2n is required.");
0097         if (k == 0)
0098         {
0099             return 1;
0100         }
0101         if (k == 1)
0102         {
0103             return x;
0104         }
0105         Real qrm2 = 1;
0106         Real qrm1 = x;
0107         Real N = 2 * m_n + 1;
0108         for (std::size_t r = 2; r <= k; ++r) {
0109             Real num = (r - 1) * (N * N - (r - 1) * (r - 1)) * qrm2;
0110             Real tmp = (2 * r - 1) * x * qrm1 - num / Real(4 * m_n * m_n);
0111             qrm2 = qrm1;
0112             qrm1 = tmp / r;
0113         }
0114         return qrm1;
0115     }
0116
0117     Real prime(Real x, std::size_t k) {
0118         BOOST_MATH_ASSERT_MSG(k <= 2 * m_n, "r <= 2n is required.");
0119         if (k == 0) {
0120             return 0;
0121         }
0122         if (k == 1) {
0123             return 1;
0124         }
0125         Real qrm2 = 1;
0126         Real qrm1 = x;
0127         Real qrm2p = 0;
0128         Real qrm1p = 1;
0129         Real N = 2 * m_n + 1;
0130         for (std::size_t r = 2; r <= k; ++r) {
0131             Real s =
0132                 (r - 1) * (N * N - (r - 1) * (r - 1)) / Real(4 * m_n * m_n);
0133             Real tmp1 = ((2 * r - 1) * x * qrm1 - s * qrm2) / r;
0134             Real tmp2 = ((2 * r - 1) * (qrm1 + x * qrm1p) - s * qrm2p) / r;
0135             qrm2 = qrm1;
0136             qrm1 = tmp1;
0137             qrm2p = qrm1p;
0138             qrm1p = tmp2;
0139         }
0140         return qrm1p;
0141     }
0142
0143   private:
0144     std::size_t m_n;
0145     std::size_t m_r;
0146     Real m_x;
0147     Real m_qrm2;
0148     Real m_qrm1;
0149     Real m_qrm2p;
0150     Real m_qrm1p;
0151     Real m_qrm2pp;
0152     Real m_qrm1pp;
0153 };
0154
0155 template <class Real>
0156 std::vector<Real> interior_velocity_filter(std::size_t n, std::size_t p) {
0157     auto dlp = discrete_legendre<Real>(n, 0);
0158     std::vector<Real> coeffs(p+1);
0159     coeffs[1] = 1/dlp.norm_sq(1);
0160     for (std::size_t l = 3; l < p + 1; l += 2)
0161     {
0162         dlp.next_prime();
0163         coeffs[l] = dlp.next_prime()/ dlp.norm_sq(l);
0164     }
0165
0166     // We could make the filter length n, as f[0] = 0,
0167     // but that'd make the indexing awkward when applying the filter.
0168     std::vector<Real> f(n + 1);
0169     // This value should never be read, but this is the correct value *if it is read*.
0170     // Hmm, should it be a nan then? I'm not gonna agonize.
0171     f[0] = 0;
0172     for (std::size_t j = 1; j < f.size(); ++j)
0173     {
0174         Real arg = Real(j) / Real(n);
0175         dlp = discrete_legendre<Real>(n, arg);
0176         f[j] = coeffs[1]*arg;
0177         for (std::size_t l = 3; l <= p; l += 2)
0178         {
0179             dlp.next();
0180             f[j] += coeffs[l]*dlp.next();
0181         }
0182         f[j] /= (n * n);
0183     }
0184     return f;
0185 }
0186
0187 template <class Real>
0188 std::vector<Real> boundary_velocity_filter(std::size_t n, std::size_t p, int64_t s)
0189 {
0190     std::vector<Real> coeffs(p+1, std::numeric_limits<Real>::quiet_NaN());
0191     Real sn = Real(s) / Real(n);
0192     auto dlp = discrete_legendre<Real>(n, sn);
0193     coeffs[0] = 0;
0194     coeffs[1] = 1/dlp.norm_sq(1);
0195     for (std::size_t l = 2; l < p + 1; ++l)
0196     {
0197         // Calculation of the norms is common to all filters,
0198         // so it seems like an obvious optimization target.
0199         // I tried this: The spent in computing the norms time is not negligible,
0200         // but still a small fraction of the total compute time.
0201         // Hence I'm not refactoring out these norm calculations.
0202         coeffs[l] = dlp.next_prime()/ dlp.norm_sq(l);
0203     }
0204
0205     std::vector<Real> f(2*n + 1);
0206     for (std::size_t k = 0; k < f.size(); ++k)
0207     {
0208         Real j = Real(k) - Real(n);
0209         Real arg = j/Real(n);
0210         dlp = discrete_legendre<Real>(n, arg);
0211         f[k] = coeffs[1]*arg;
0212         for (std::size_t l = 2; l <= p; ++l)
0213         {
0214             f[k] += coeffs[l]*dlp.next();
0215         }
0216         f[k] /= (n * n);
0217     }
0218     return f;
0219 }
0220
0221 template <class Real>
0222 std::vector<Real> acceleration_filter(std::size_t n, std::size_t p, int64_t s)
0223 {
0224     BOOST_MATH_ASSERT_MSG(p <= 2*n, "Approximation order must be <= 2*n");
0225     BOOST_MATH_ASSERT_MSG(p > 2, "Approximation order must be > 2");
0226
0227     std::vector<Real> coeffs(p+1, std::numeric_limits<Real>::quiet_NaN());
0228     Real sn = Real(s) / Real(n);
0229     auto dlp = discrete_legendre<Real>(n, sn);
0230     coeffs[0] = 0;
0231     coeffs[1] = 0;
0232     for (std::size_t l = 2; l < p + 1; ++l)
0233     {
0234         coeffs[l] = dlp.next_dbl_prime()/ dlp.norm_sq(l);
0235     }
0236
0237     std::vector<Real> f(2*n + 1, 0);
0238     for (std::size_t k = 0; k < f.size(); ++k)
0239     {
0240         Real j = Real(k) - Real(n);
0241         Real arg = j/Real(n);
0242         dlp = discrete_legendre<Real>(n, arg);
0243         for (std::size_t l = 2; l <= p; ++l)
0244         {
0245             f[k] += coeffs[l]*dlp.next();
0246         }
0247         f[k] /= (n * n * n);
0248     }
0249     return f;
0250 }
0251
0252
0253 } // namespace detail
0254
0255 template <typename Real, std::size_t order = 1>
0256 class discrete_lanczos_derivative {
0257 public:
0258     discrete_lanczos_derivative(Real const & spacing,
0259                                 std::size_t n = 18,
0260                                 std::size_t approximation_order = 3)
0261         : m_dt{spacing}
0262     {
0263         static_assert(!std::is_integral_v<Real>,
0264                       "Spacing must be a floating point type.");
0265         BOOST_MATH_ASSERT_MSG(spacing > 0,
0266                          "Spacing between samples must be > 0.");
0267
0268         if constexpr (order == 1)
0269         {
0270             BOOST_MATH_ASSERT_MSG(approximation_order <= 2 * n,
0271                              "The approximation order must be <= 2n");
0272             BOOST_MATH_ASSERT_MSG(approximation_order >= 2,
0273                              "The approximation order must be >= 2");
0274
0275             if constexpr (std::is_same_v<Real, float> || std::is_same_v<Real, double>)
0276             {
0277                 auto interior = detail::interior_velocity_filter<long double>(n, approximation_order);
0278                 m_f.resize(interior.size());
0279                 for (std::size_t j = 0; j < interior.size(); ++j)
0280                 {
0281                     m_f[j] = static_cast<Real>(interior[j])/m_dt;
0282                 }
0283             }
0284             else
0285             {
0286                 m_f = detail::interior_velocity_filter<Real>(n, approximation_order);
0287                 for (auto & x : m_f)
0288                 {
0289                     x /= m_dt;
0290                 }
0291             }
0292
0293             m_boundary_filters.resize(n);
0294             // This for loop is a natural candidate for parallelization.
0295             // But does it matter? Probably not.
0296             for (std::size_t i = 0; i < n; ++i)
0297             {
0298                 if constexpr (std::is_same_v<Real, float> || std::is_same_v<Real, double>)
0299                 {
0300                     int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);
0301                     auto bf = detail::boundary_velocity_filter<long double>(n, approximation_order, s);
0302                     m_boundary_filters[i].resize(bf.size());
0303                     for (std::size_t j = 0; j < bf.size(); ++j)
0304                     {
0305                         m_boundary_filters[i][j] = static_cast<Real>(bf[j])/m_dt;
0306                     }
0307                 }
0308                 else
0309                 {
0310                     int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);
0311                     m_boundary_filters[i] = detail::boundary_velocity_filter<Real>(n, approximation_order, s);
0312                     for (auto & bf : m_boundary_filters[i])
0313                     {
0314                         bf /= m_dt;
0315                     }
0316                 }
0317             }
0318         }
0319         else if constexpr (order == 2)
0320         {
0321             // High precision isn't warranted for small p; only for large p.
0322             // (The computation appears stable for large n.)
0323             // But given that the filters are reusable for many vectors,
0324             // it's better to do a high precision computation and then cast back,
0325             // since the resulting cost is a factor of 2, and the cost of the filters not working is hours of debugging.
0326             if constexpr (std::is_same_v<Real, double> || std::is_same_v<Real, float>)
0327             {
0328                 auto f = detail::acceleration_filter<long double>(n, approximation_order, 0);
0329                 m_f.resize(n+1);
0330                 for (std::size_t i = 0; i < m_f.size(); ++i)
0331                 {
0332                     m_f[i] = static_cast<Real>(f[i+n])/(m_dt*m_dt);
0333                 }
0334                 m_boundary_filters.resize(n);
0335                 for (std::size_t i = 0; i < n; ++i)
0336                 {
0337                     int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);
0338                     auto bf = detail::acceleration_filter<long double>(n, approximation_order, s);
0339                     m_boundary_filters[i].resize(bf.size());
0340                     for (std::size_t j = 0; j < bf.size(); ++j)
0341                     {
0342                         m_boundary_filters[i][j] = static_cast<Real>(bf[j])/(m_dt*m_dt);
0343                     }
0344                 }
0345             }
0346             else
0347             {
0348                 // Given that the purpose is denoising, for higher precision calculations,
0349                 // the default precision should be fine.
0350                 auto f = detail::acceleration_filter<Real>(n, approximation_order, 0);
0351                 m_f.resize(n+1);
0352                 for (std::size_t i = 0; i < m_f.size(); ++i)
0353                 {
0354                     m_f[i] = f[i+n]/(m_dt*m_dt);
0355                 }
0356                 m_boundary_filters.resize(n);
0357                 for (std::size_t i = 0; i < n; ++i)
0358                 {
0359                     int64_t s = static_cast<int64_t>(i) - static_cast<int64_t>(n);
0360                     m_boundary_filters[i] = detail::acceleration_filter<Real>(n, approximation_order, s);
0361                     for (auto & bf : m_boundary_filters[i])
0362                     {
0363                         bf /= (m_dt*m_dt);
0364                     }
0365                 }
0366             }
0367         }
0368         else
0369         {
0370             BOOST_MATH_ASSERT_MSG(false, "Derivatives of order 3 and higher are not implemented.");
0371         }
0372     }
0373
0374     Real get_spacing() const
0375     {
0376         return m_dt;
0377     }
0378
0379     template<class RandomAccessContainer>
0380     Real operator()(RandomAccessContainer const & v, std::size_t i) const
0381     {
0382         static_assert(std::is_same_v<typename RandomAccessContainer::value_type, Real>,
0383                       "The type of the values in the vector provided does not match the type in the filters.");
0384
0385         BOOST_MATH_ASSERT_MSG(std::size(v) >= m_boundary_filters[0].size(),
0386             "Vector must be at least as long as the filter length");
0387
0388         if constexpr (order==1)
0389         {
0390             if (i >= m_f.size() - 1 && i <= std::size(v) - m_f.size())
0391             {
0392                 // The filter has length >= 1:
0393                 Real dvdt = m_f[1] * (v[i + 1] - v[i - 1]);
0394                 for (std::size_t j = 2; j < m_f.size(); ++j)
0395                 {
0396                     dvdt += m_f[j] * (v[i + j] - v[i - j]);
0397                 }
0398                 return dvdt;
0399             }
0400
0401             // m_f.size() = N+1
0402             if (i < m_f.size() - 1)
0403             {
0404                 auto &bf = m_boundary_filters[i];
0405                 Real dvdt = bf[0]*v[0];
0406                 for (std::size_t j = 1; j < bf.size(); ++j)
0407                 {
0408                     dvdt += bf[j] * v[j];
0409                 }
0410                 return dvdt;
0411             }
0412
0413             if (i > std::size(v) - m_f.size() && i < std::size(v))
0414             {
0415                 int k = std::size(v) - 1 - i;
0416                 auto &bf = m_boundary_filters[k];
0417                 Real dvdt = bf[0]*v[std::size(v)-1];
0418                 for (std::size_t j = 1; j < bf.size(); ++j)
0419                 {
0420                     dvdt += bf[j] * v[std::size(v) - 1 - j];
0421                 }
0422                 return -dvdt;
0423             }
0424         }
0425         else if constexpr (order==2)
0426         {
0427             if (i >= m_f.size() - 1 && i <= std::size(v) - m_f.size())
0428             {
0429                 Real d2vdt2 = m_f[0]*v[i];
0430                 for (std::size_t j = 1; j < m_f.size(); ++j)
0431                 {
0432                     d2vdt2 += m_f[j] * (v[i + j] + v[i - j]);
0433                 }
0434                 return d2vdt2;
0435             }
0436
0437             // m_f.size() = N+1
0438             if (i < m_f.size() - 1)
0439             {
0440                 auto &bf = m_boundary_filters[i];
0441                 Real d2vdt2 = bf[0]*v[0];
0442                 for (std::size_t j = 1; j < bf.size(); ++j)
0443                 {
0444                     d2vdt2 += bf[j] * v[j];
0445                 }
0446                 return d2vdt2;
0447             }
0448
0449             if (i > std::size(v) - m_f.size() && i < std::size(v))
0450             {
0451                 int k = std::size(v) - 1 - i;
0452                 auto &bf = m_boundary_filters[k];
0453                 Real d2vdt2 = bf[0] * v[std::size(v) - 1];
0454                 for (std::size_t j = 1; j < bf.size(); ++j)
0455                 {
0456                     d2vdt2 += bf[j] * v[std::size(v) - 1 - j];
0457                 }
0458                 return d2vdt2;
0459             }
0460         }
0461
0462         // OOB access:
0463         std::string msg = "Out of bounds access in Lanczos derivative.";
0464         msg += "Input vector has length " + std::to_string(std::size(v)) + ", but user requested access at index " + std::to_string(i) + ".";
0465         throw std::out_of_range(msg);
0466         return std::numeric_limits<Real>::quiet_NaN();
0467     }
0468
0469     template<class RandomAccessContainer>
0470     void operator()(RandomAccessContainer const & v, RandomAccessContainer & w) const
0471     {
0472         static_assert(std::is_same_v<typename RandomAccessContainer::value_type, Real>,
0473                       "The type of the values in the vector provided does not match the type in the filters.");
0474         if (&w[0] == &v[0])
0475         {
0476             throw std::logic_error("This transform cannot be performed in-place.");
0477         }
0478
0479         if (std::size(v) < m_boundary_filters[0].size())
0480         {
0481             std::string msg = "The input vector must be at least as long as the filter length. ";
0482             msg += "The input vector has length = " + std::to_string(std::size(v)) + ", the filter has length " + std::to_string(m_boundary_filters[0].size());
0483             throw std::length_error(msg);
0484         }
0485
0486         if (std::size(w) < std::size(v))
0487         {
0488             std::string msg = "The output vector (containing the derivative) must be at least as long as the input vector.";
0489             msg += "The output vector has length = " + std::to_string(std::size(w)) + ", the input vector has length " + std::to_string(std::size(v));
0490             throw std::length_error(msg);
0491         }
0492
0493         if constexpr (order==1)
0494         {
0495             for (std::size_t i = 0; i < m_f.size() - 1; ++i)
0496             {
0497                 auto &bf = m_boundary_filters[i];
0498                 Real dvdt = bf[0] * v[0];
0499                 for (std::size_t j = 1; j < bf.size(); ++j)
0500                 {
0501                     dvdt += bf[j] * v[j];
0502                 }
0503                 w[i] = dvdt;
0504             }
0505
0506             for(std::size_t i = m_f.size() - 1; i <= std::size(v) - m_f.size(); ++i)
0507             {
0508                 Real dvdt = m_f[1] * (v[i + 1] - v[i - 1]);
0509                 for (std::size_t j = 2; j < m_f.size(); ++j)
0510                 {
0511                     dvdt += m_f[j] *(v[i + j] - v[i - j]);
0512                 }
0513                 w[i] = dvdt;
0514             }
0515
0516
0517             for(std::size_t i = std::size(v) - m_f.size() + 1; i < std::size(v); ++i)
0518             {
0519                 int k = std::size(v) - 1 - i;
0520                 auto &f = m_boundary_filters[k];
0521                 Real dvdt = f[0] * v[std::size(v) - 1];;
0522                 for (std::size_t j = 1; j < f.size(); ++j)
0523                 {
0524                     dvdt += f[j] * v[std::size(v) - 1 - j];
0525                 }
0526                 w[i] = -dvdt;
0527             }
0528         }
0529         else if constexpr (order==2)
0530         {
0531             // m_f.size() = N+1
0532             for (std::size_t i = 0; i < m_f.size() - 1; ++i)
0533             {
0534                 auto &bf = m_boundary_filters[i];
0535                 Real d2vdt2 = 0;
0536                 for (std::size_t j = 0; j < bf.size(); ++j)
0537                 {
0538                     d2vdt2 += bf[j] * v[j];
0539                 }
0540                 w[i] = d2vdt2;
0541             }
0542
0543             for (std::size_t i = m_f.size() - 1; i <= std::size(v) - m_f.size(); ++i)
0544             {
0545                 Real d2vdt2 = m_f[0]*v[i];
0546                 for (std::size_t j = 1; j < m_f.size(); ++j)
0547                 {
0548                     d2vdt2 += m_f[j] * (v[i + j] + v[i - j]);
0549                 }
0550                 w[i] = d2vdt2;
0551             }
0552
0553             for (std::size_t i = std::size(v) - m_f.size() + 1; i < std::size(v); ++i)
0554             {
0555                 int k = std::size(v) - 1 - i;
0556                 auto &bf = m_boundary_filters[k];
0557                 Real d2vdt2 = bf[0] * v[std::size(v) - 1];
0558                 for (std::size_t j = 1; j < bf.size(); ++j)
0559                 {
0560                     d2vdt2 += bf[j] * v[std::size(v) - 1 - j];
0561                 }
0562                 w[i] = d2vdt2;
0563             }
0564         }
0565     }
0566
0567     template<class RandomAccessContainer>
0568     RandomAccessContainer operator()(RandomAccessContainer const & v) const
0569     {
0570         RandomAccessContainer w(std::size(v));
0571         this->operator()(v, w);
0572         return w;
0573     }
0574
0575
0576     // Don't copy; too big.
0577     discrete_lanczos_derivative( const discrete_lanczos_derivative & ) = delete;
0578     discrete_lanczos_derivative& operator=(const discrete_lanczos_derivative&) = delete;
0579
0580     // Allow moves:
0581     discrete_lanczos_derivative(discrete_lanczos_derivative&&) noexcept = default;
0582     discrete_lanczos_derivative& operator=(discrete_lanczos_derivative&&) noexcept = default;
0583
0584 private:
0585     std::vector<Real> m_f;
0586     std::vector<std::vector<Real>> m_boundary_filters;
0587     Real m_dt;
0588 };
0589
0590 } // namespaces
0591 #endif