math/special_functions/daubechies_wavelet.hpp

0001 /*
0002  * Copyright Nick Thompson, 2020
0003  * Use, modification and distribution are subject to the
0004  * Boost Software License, Version 1.0. (See accompanying file
0005  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006  */
0007
0008 #ifndef BOOST_MATH_SPECIAL_DAUBECHIES_WAVELET_HPP
0009 #define BOOST_MATH_SPECIAL_DAUBECHIES_WAVELET_HPP
0010 #include <vector>
0011 #include <array>
0012 #include <cmath>
0013 #include <thread>
0014 #include <future>
0015 #include <iostream>
0016 #include <boost/math/constants/constants.hpp>
0017 #include <boost/math/special_functions/detail/daubechies_scaling_integer_grid.hpp>
0018 #include <boost/math/special_functions/daubechies_scaling.hpp>
0019 #include <boost/math/filters/daubechies.hpp>
0020 #include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
0021 #include <boost/math/interpolators/detail/quintic_hermite_detail.hpp>
0022 #include <boost/math/interpolators/detail/septic_hermite_detail.hpp>
0023
0024 #include <boost/math/tools/is_standalone.hpp>
0025 #ifndef BOOST_MATH_STANDALONE
0026 #include <boost/config.hpp>
0027 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
0028 #error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
0029 #endif
0030 #endif
0031
0032 namespace boost::math {
0033
0034    template<class Real, int p, int order>
0035    std::vector<Real> daubechies_wavelet_dyadic_grid(int64_t j_max)
0036    {
0037       if (j_max == 0)
0038       {
0039          throw std::domain_error("The wavelet dyadic grid is refined from the scaling integer grid, so its minimum amount of data is half integer widths.");
0040       }
0041       auto phijk = daubechies_scaling_dyadic_grid<Real, p, order>(j_max - 1);
0042       //psi_j[l] = psi(-p+1 + l/2^j) = \sum_{k=0}^{2p-1} (-1)^k c_k \phi(1-2p+k + l/2^{j-1})
0043       //For derivatives just map c_k -> 2^order c_k.
0044       auto d = boost::math::filters::daubechies_scaling_filter<Real, p>();
0045       Real scale = boost::math::constants::root_two<Real>() * (1 << order);
0046       for (size_t i = 0; i < d.size(); ++i)
0047       {
0048          d[i] *= scale;
0049          if (!(i & 1))
0050          {
0051             d[i] = -d[i];
0052          }
0053       }
0054
0055       std::vector<Real> v(2 * p + (2 * p - 1) * ((int64_t(1) << j_max) - 1), std::numeric_limits<Real>::quiet_NaN());
0056       v[0] = 0;
0057       v[v.size() - 1] = 0;
0058
0059       for (int64_t l = 1; l < static_cast<int64_t>(v.size() - 1); ++l)
0060       {
0061          Real term = 0;
0062          for (int64_t k = 0; k < static_cast<int64_t>(d.size()); ++k)
0063          {
0064             int64_t idx = (int64_t(1) << (j_max - 1)) * (1 - 2 * p + k) + l;
0065             if (idx < 0 || idx >= static_cast<int64_t>(phijk.size()))
0066             {
0067                continue;
0068             }
0069             term += d[k] * phijk[idx];
0070          }
0071          v[l] = term;
0072       }
0073
0074       return v;
0075    }
0076
0077
0078    template<class Real, int p>
0079    class daubechies_wavelet {
0080       //
0081       // Some type manipulation so we know the type of the interpolator, and the vector type it requires:
0082       //
0083       using vector_type = std::vector < std::array < Real, p < 6 ? 2 : p < 10 ? 3 : 4>>;
0084       //
0085       // List our interpolators:
0086       //
0087       using interpolator_list = std::tuple<
0088          detail::null_interpolator, detail::matched_holder_aos<vector_type>, detail::linear_interpolation_aos<vector_type>,
0089          interpolators::detail::cardinal_cubic_hermite_detail_aos<vector_type>, interpolators::detail::cardinal_quintic_hermite_detail_aos<vector_type>,
0090          interpolators::detail::cardinal_septic_hermite_detail_aos<vector_type> > ;
0091       //
0092       // Select the one we need:
0093       //
0094       using interpolator_type = std::tuple_element_t<
0095          p == 1 ? 0 :
0096          p == 2 ? 1 :
0097          p == 3 ? 2 :
0098          p <= 5 ? 3 :
0099          p <= 9 ? 4 : 5, interpolator_list>;
0100    public:
0101       explicit daubechies_wavelet(int grid_refinements = -1)
0102       {
0103          static_assert(p < 20, "Daubechies wavelets are only implemented for p < 20.");
0104          static_assert(p > 0, "Daubechies wavelets must have at least 1 vanishing moment.");
0105          if (grid_refinements == 0)
0106          {
0107             throw std::domain_error("The wavelet requires at least 1 grid refinement.");
0108          }
0109          if constexpr (p == 1)
0110          {
0111             return;
0112          }
0113          else
0114          {
0115             if (grid_refinements < 0)
0116             {
0117                if constexpr (std::is_same_v<Real, float>)
0118                {
0119                   if (grid_refinements == -2)
0120                   {
0121                      // Control absolute error:
0122                      //                          p= 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
0123                      std::array<int, 20> r{ -1, -1, 18, 19, 16, 11,  8,  7,  7,  7,  5,  5,  4,  4,  4,  4,  3,  3,  3,  3 };
0124                      grid_refinements = r[p];
0125                   }
0126                   else
0127                   {
0128                      // Control relative error:
0129                      //                          p= 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
0130                      std::array<int, 20> r{ -1, -1, 21, 21, 21, 17, 16, 15, 14, 13, 12, 11, 11, 11, 11, 11, 11, 11, 11, 11 };
0131                      grid_refinements = r[p];
0132                   }
0133                }
0134                else if constexpr (std::is_same_v<Real, double>)
0135                {
0136                   //                          p= 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
0137                   std::array<int, 20> r{ -1, -1, 21, 21, 21, 21, 21, 21, 21, 21, 20, 20, 19, 18, 18, 18, 18, 18, 18, 18 };
0138                   grid_refinements = r[p];
0139                }
0140                else
0141                {
0142                   grid_refinements = 21;
0143                }
0144             }
0145
0146             // Compute the refined grid:
0147             // In fact for float precision I know the grid must be computed in double precision and then cast back down, or else parts of the support are systematically inaccurate.
0148             std::future<std::vector<Real>> t0 = std::async(std::launch::async, [&grid_refinements]() {
0149                // Computing in higher precision and downcasting is essential for 1ULP evaluation in float precision:
0150                auto v = daubechies_wavelet_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 0>(grid_refinements);
0151                return detail::daubechies_eval_type<Real>::vector_cast(v);
0152                });
0153             // Compute the derivative of the refined grid:
0154             std::future<std::vector<Real>> t1 = std::async(std::launch::async, [&grid_refinements]() {
0155                auto v = daubechies_wavelet_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 1>(grid_refinements);
0156                return detail::daubechies_eval_type<Real>::vector_cast(v);
0157                });
0158
0159             // if necessary, compute the second and third derivative:
0160             std::vector<Real> d2ydx2;
0161             std::vector<Real> d3ydx3;
0162             if constexpr (p >= 6) {
0163                std::future<std::vector<Real>> t3 = std::async(std::launch::async, [&grid_refinements]() {
0164                   auto v = daubechies_wavelet_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 2>(grid_refinements);
0165                   return detail::daubechies_eval_type<Real>::vector_cast(v);
0166                   });
0167
0168                if constexpr (p >= 10) {
0169                   std::future<std::vector<Real>> t4 = std::async(std::launch::async, [&grid_refinements]() {
0170                      auto v = daubechies_wavelet_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 3>(grid_refinements);
0171                      return detail::daubechies_eval_type<Real>::vector_cast(v);
0172                      });
0173                   d3ydx3 = t4.get();
0174                }
0175                d2ydx2 = t3.get();
0176             }
0177
0178
0179             auto y = t0.get();
0180             auto dydx = t1.get();
0181
0182             if constexpr (p >= 2)
0183             {
0184                vector_type data(y.size());
0185                for (size_t i = 0; i < y.size(); ++i)
0186                {
0187                   data[i][0] = y[i];
0188                   data[i][1] = dydx[i];
0189                   if constexpr (p >= 6)
0190                      data[i][2] = d2ydx2[i];
0191                   if constexpr (p >= 10)
0192                      data[i][3] = d3ydx3[i];
0193                }
0194                if constexpr (p <= 3)
0195                   m_interpolator = std::make_shared<interpolator_type>(std::move(data), grid_refinements, Real(-p + 1));
0196                else
0197                   m_interpolator = std::make_shared<interpolator_type>(std::move(data), Real(-p + 1), Real(1) / (1 << grid_refinements));
0198             }
0199             else
0200                m_interpolator = std::make_shared<detail::null_interpolator>();
0201          }
0202       }
0203
0204
0205       inline Real operator()(Real x) const
0206       {
0207          if (x <= -p + 1 || x >= p)
0208          {
0209             return 0;
0210          }
0211
0212          if constexpr (p == 1)
0213          {
0214             if (x < Real(1) / Real(2))
0215             {
0216                return 1;
0217             }
0218             else if (x == Real(1) / Real(2))
0219             {
0220                return 0;
0221             }
0222             return -1;
0223          }
0224          else
0225          {
0226             return (*m_interpolator)(x);
0227          }
0228       }
0229
0230       inline Real prime(Real x) const
0231       {
0232          static_assert(p > 2, "The 3-vanishing moment Daubechies wavelet is the first which is continuously differentiable.");
0233          if (x <= -p + 1 || x >= p)
0234          {
0235             return 0;
0236          }
0237          return m_interpolator->prime(x);
0238       }
0239
0240       inline Real double_prime(Real x) const
0241       {
0242          static_assert(p >= 6, "Second derivatives of Daubechies wavelets require at least 6 vanishing moments.");
0243          if (x <= -p + 1 || x >= p)
0244          {
0245             return Real(0);
0246          }
0247          return m_interpolator->double_prime(x);
0248       }
0249
0250       std::pair<Real, Real> support() const
0251       {
0252          return std::make_pair(Real(-p + 1), Real(p));
0253       }
0254
0255       int64_t bytes() const
0256       {
0257          return m_interpolator->bytes() + sizeof(*this);
0258       }
0259
0260    private:
0261       std::shared_ptr<interpolator_type> m_interpolator;
0262    };
0263
0264 }
0265
0266 #endif // BOOST_MATH_SPECIAL_DAUBECHIES_WAVELET_HPP