math/special_functions/daubechies_scaling.hpp

0001 /*
0002  * Copyright Nick Thompson, John Maddock 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_SCALING_HPP
0009 #define BOOST_MATH_SPECIAL_DAUBECHIES_SCALING_HPP
0010
0011 #include <cstdint>
0012 #include <cstring>
0013 #include <cmath>
0014 #include <vector>
0015 #include <array>
0016 #include <thread>
0017 #include <future>
0018 #include <iostream>
0019 #include <memory>
0020 #include <boost/math/special_functions/detail/daubechies_scaling_integer_grid.hpp>
0021 #include <boost/math/filters/daubechies.hpp>
0022 #include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
0023 #include <boost/math/interpolators/detail/quintic_hermite_detail.hpp>
0024 #include <boost/math/interpolators/detail/septic_hermite_detail.hpp>
0025
0026 #include <boost/math/tools/is_standalone.hpp>
0027 #ifndef BOOST_MATH_STANDALONE
0028 #  include <boost/config.hpp>
0029 #  ifdef BOOST_NO_CXX17_IF_CONSTEXPR
0030 #    error "The header <boost/special_functions/daubechies_scaling.hpp> can only be used in C++17 and later."
0031 #  endif
0032 #endif
0033
0034 namespace boost::math {
0035
0036 template<class Real, int p, int order>
0037 std::vector<Real> daubechies_scaling_dyadic_grid(int64_t j_max)
0038 {
0039     using std::isnan;
0040     using std::sqrt;
0041     auto c = boost::math::filters::daubechies_scaling_filter<Real, p>();
0042     Real scale = sqrt(static_cast<Real>(2))*(1 << order);
0043     for (auto & x : c)
0044     {
0045         x *= scale;
0046     }
0047
0048     auto phik = detail::daubechies_scaling_integer_grid<Real, p, order>();
0049
0050     // Maximum sensible j for 32 bit floats is j_max = 22:
0051     if (std::is_same_v<Real, float>)
0052     {
0053         if (j_max > 23)
0054         {
0055             throw std::logic_error("Requested dyadic grid more dense than number of representables on the interval.");
0056         }
0057     }
0058     std::vector<Real> v(2*p + (2*p-1)*((1<<j_max) -1), std::numeric_limits<Real>::quiet_NaN());
0059     v[0] = 0;
0060     v[v.size()-1] = 0;
0061     for (int64_t i = 0; i < static_cast<int64_t>(phik.size()); ++i) {
0062         v[i*(1uLL<<j_max)] = phik[i];
0063     }
0064
0065     for (int64_t j = 1; j <= j_max; ++j)
0066     {
0067         int64_t k_max = v.size()/(int64_t(1) << (j_max-j));
0068         for (int64_t k = 1; k < k_max;  k += 2)
0069         {
0070             // Where this value will go:
0071             int64_t delivery_idx = k*(1uLL << (j_max-j));
0072             // This is a nice check, but we've tested this exhaustively, and it's an expensive check:
0073             //if (delivery_idx >= static_cast<int64_t>(v.size())) {
0074             //    std::cerr << "Delivery index out of range!\n";
0075             //    continue;
0076             //}
0077             Real term = 0;
0078             for (int64_t l = 0; l < static_cast<int64_t>(c.size()); ++l)
0079             {
0080                 int64_t idx = k*(int64_t(1) << (j_max - j + 1)) - l*(int64_t(1) << j_max);
0081                 if (idx < 0)
0082                 {
0083                     break;
0084                 }
0085                 if (idx < static_cast<int64_t>(v.size()))
0086                 {
0087                     term += c[l]*v[idx];
0088                 }
0089             }
0090             // Again, another nice check:
0091             //if (!isnan(v[delivery_idx])) {
0092             //    std::cerr << "Delivery index already populated!, = " << v[delivery_idx] << "\n";
0093             //    std::cerr << "would overwrite with " << term << "\n";
0094             //}
0095             v[delivery_idx] = term;
0096         }
0097     }
0098     return v;
0099 }
0100
0101 namespace detail {
0102
0103 template<class RandomAccessContainer>
0104 class matched_holder {
0105 public:
0106     using Real = typename RandomAccessContainer::value_type;
0107
0108     matched_holder(RandomAccessContainer && y, RandomAccessContainer && dydx, int grid_refinements, Real x0) : x0_{x0}, y_{std::move(y)}, dy_{std::move(dydx)}
0109     {
0110         inv_h_ = (1 << grid_refinements);
0111         Real h = 1/inv_h_;
0112         for (auto & dy : dy_)
0113         {
0114             dy *= h;
0115         }
0116     }
0117
0118     inline Real operator()(Real x) const
0119     {
0120         using std::floor;
0121         using std::sqrt;
0122         // This is the exact Holder exponent, but it's pessimistic almost everywhere!
0123         // It's only exactly right at dyadic rationals.
0124         //Real const alpha = 2 - log(1+sqrt(Real(3)))/log(Real(2));
0125         // We're gonna use alpha = 1/2, rather than 0.5500...
0126         Real s = (x-x0_)*inv_h_;
0127         Real ii = floor(s);
0128         auto i = static_cast<decltype(y_.size())>(ii);
0129         Real t = s - ii;
0130         Real dphi = dy_[i+1];
0131         Real diff = y_[i+1] - y_[i];
0132         return y_[i] + (2*dphi - diff)*t + 2*sqrt(t)*(diff-dphi);
0133     }
0134
0135     int64_t bytes() const
0136     {
0137         return 2*y_.size()*sizeof(Real) + sizeof(*this);
0138     }
0139
0140 private:
0141     Real x0_;
0142     Real inv_h_;
0143     RandomAccessContainer y_;
0144     RandomAccessContainer dy_;
0145 };
0146
0147 template<class RandomAccessContainer>
0148 class matched_holder_aos {
0149 public:
0150     using Point = typename RandomAccessContainer::value_type;
0151     using Real = typename Point::value_type;
0152
0153     matched_holder_aos(RandomAccessContainer && data, int grid_refinements, Real x0) : x0_{x0}, data_{std::move(data)}
0154     {
0155         inv_h_ = Real(1uLL << grid_refinements);
0156         Real h = 1/inv_h_;
0157         for (auto & datum : data_)
0158         {
0159             datum[1] *= h;
0160         }
0161     }
0162
0163     inline Real operator()(Real x) const
0164     {
0165         using std::floor;
0166         using std::sqrt;
0167         Real s = (x-x0_)*inv_h_;
0168         Real ii = floor(s);
0169         auto i = static_cast<decltype(data_.size())>(ii);
0170         Real t = s - ii;
0171         Real y0 = data_[i][0];
0172         Real y1 = data_[i+1][0];
0173         Real dphi = data_[i+1][1];
0174         Real diff = y1 - y0;
0175         return y0 + (2*dphi - diff)*t + 2*sqrt(t)*(diff-dphi);
0176     }
0177
0178     int64_t bytes() const
0179     {
0180         return data_.size()*data_[0].size()*sizeof(Real) + sizeof(*this);
0181     }
0182
0183 private:
0184     Real x0_;
0185     Real inv_h_;
0186     RandomAccessContainer data_;
0187 };
0188
0189
0190 template<class RandomAccessContainer>
0191 class linear_interpolation {
0192 public:
0193     using Real = typename RandomAccessContainer::value_type;
0194
0195     linear_interpolation(RandomAccessContainer && y, RandomAccessContainer && dydx, int grid_refinements) : y_{std::move(y)}, dydx_{std::move(dydx)}
0196     {
0197         s_ = (1 << grid_refinements);
0198     }
0199
0200     inline Real operator()(Real x) const
0201     {
0202         using std::floor;
0203         Real y = x*s_;
0204         Real k = floor(y);
0205
0206         int64_t kk = static_cast<int64_t>(k);
0207         Real t = y - k;
0208         return (1-t)*y_[kk] + t*y_[kk+1];
0209     }
0210
0211     inline Real prime(Real x) const
0212     {
0213         using std::floor;
0214
0215         Real y = x*s_;
0216         Real k = floor(y);
0217
0218         int64_t kk = static_cast<int64_t>(k);
0219         Real t = y - k;
0220         return static_cast<Real>((Real(1)-t)*dydx_[kk] + t*dydx_[kk+1]);
0221     }
0222
0223     int64_t bytes() const
0224     {
0225         return (1 + y_.size() + dydx_.size())*sizeof(Real) + sizeof(y_) + sizeof(dydx_);
0226     }
0227
0228 private:
0229     Real s_;
0230     RandomAccessContainer y_;
0231     RandomAccessContainer dydx_;
0232 };
0233
0234 template<class RandomAccessContainer>
0235 class linear_interpolation_aos {
0236 public:
0237     using Point = typename RandomAccessContainer::value_type;
0238     using Real = typename Point::value_type;
0239
0240     linear_interpolation_aos(RandomAccessContainer && data, int grid_refinements, Real x0) : x0_{x0}, data_{std::move(data)}
0241     {
0242         s_ = Real(1uLL << grid_refinements);
0243     }
0244
0245     inline Real operator()(Real x) const
0246     {
0247         using std::floor;
0248         Real y = (x-x0_)*s_;
0249         Real k = floor(y);
0250
0251         int64_t kk = static_cast<int64_t>(k);
0252         Real t = y - k;
0253         return (t != 0) ? (1-t)*data_[kk][0] + t*data_[kk+1][0] : data_[kk][0];
0254     }
0255
0256     inline Real prime(Real x) const
0257     {
0258         using std::floor;
0259         Real y = (x-x0_)*s_;
0260         Real k = floor(y);
0261
0262         int64_t kk = static_cast<int64_t>(k);
0263         Real t = y - k;
0264         return t != 0 ? (1-t)*data_[kk][1] + t*data_[kk+1][1] : data_[kk][1];
0265     }
0266
0267     int64_t bytes() const
0268     {
0269         return sizeof(*this) + data_.size()*data_[0].size()*sizeof(Real);
0270     }
0271
0272 private:
0273     Real x0_;
0274     Real s_;
0275     RandomAccessContainer data_;
0276 };
0277
0278
0279 template <class T>
0280 struct daubechies_eval_type
0281 {
0282     using type = T;
0283
0284     static const std::vector<T>& vector_cast(const std::vector<T>& v) { return v; }
0285
0286 };
0287 template <>
0288 struct daubechies_eval_type<float>
0289 {
0290     using type = double;
0291
0292     inline static std::vector<float> vector_cast(const std::vector<double>& v)
0293     {
0294         std::vector<float> result(v.size());
0295         for (unsigned i = 0; i < v.size(); ++i)
0296             result[i] = static_cast<float>(v[i]);
0297         return result;
0298     }
0299 };
0300 template <>
0301 struct daubechies_eval_type<double>
0302 {
0303     using type = long double;
0304
0305     inline static std::vector<double> vector_cast(const std::vector<long double>& v)
0306     {
0307         std::vector<double> result(v.size());
0308         for (unsigned i = 0; i < v.size(); ++i)
0309             result[i] = static_cast<double>(v[i]);
0310         return result;
0311     }
0312 };
0313
0314 struct null_interpolator
0315 {
0316     template <class T>
0317     T operator()(const T&)
0318     {
0319         return 1;
0320     }
0321 };
0322
0323 } // namespace detail
0324
0325 template<class Real, int p>
0326 class daubechies_scaling {
0327     //
0328     // Some type manipulation so we know the type of the interpolator, and the vector type it requires:
0329     //
0330     using vector_type = std::vector<std::array<Real, p < 6 ? 2 : p < 10 ? 3 : 4>>;
0331     //
0332     // List our interpolators:
0333     //
0334     using interpolator_list = std::tuple<
0335         detail::null_interpolator, detail::matched_holder_aos<vector_type>, detail::linear_interpolation_aos<vector_type>,
0336         interpolators::detail::cardinal_cubic_hermite_detail_aos<vector_type>, interpolators::detail::cardinal_quintic_hermite_detail_aos<vector_type>,
0337         interpolators::detail::cardinal_septic_hermite_detail_aos<vector_type> >;
0338     //
0339     // Select the one we need:
0340     //
0341     using interpolator_type = std::tuple_element_t<
0342         p == 1 ? 0 :
0343         p == 2 ? 1 :
0344         p == 3 ? 2 :
0345         p <= 5 ? 3 :
0346         p <= 9 ? 4 : 5, interpolator_list>;
0347
0348 public:
0349     daubechies_scaling(int grid_refinements = -1)
0350     {
0351         static_assert(p < 20, "Daubechies scaling functions are only implemented for p < 20.");
0352         static_assert(p > 0, "Daubechies scaling functions must have at least 1 vanishing moment.");
0353         if constexpr (p == 1)
0354         {
0355             return;
0356         }
0357         else {
0358             if (grid_refinements < 0)
0359             {
0360                 if (std::is_same_v<Real, float>)
0361                 {
0362                 if (grid_refinements == -2)
0363                 {
0364                     // Control absolute error:
0365                     //                          p= 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
0366                     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 };
0367                     grid_refinements = r[p];
0368                 }
0369                 else
0370                 {
0371                     // Control relative error:
0372                     //                          p= 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
0373                     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 };
0374                     grid_refinements = r[p];
0375                 }
0376                 }
0377                 else if (std::is_same_v<Real, double>)
0378                 {
0379                     //                          p= 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
0380                     std::array<int, 20> r{ -1, -1, 21, 21, 21, 21, 21, 21, 21, 21, 20, 20, 19, 19, 18, 18, 18, 18, 18, 18 };
0381                     grid_refinements = r[p];
0382                 }
0383                 else
0384                 {
0385                     grid_refinements = 21;
0386                 }
0387             }
0388
0389             // Compute the refined grid:
0390             // 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.
0391             std::future<std::vector<Real>> t0 = std::async(std::launch::async, [&grid_refinements]() {
0392                 // Computing in higher precision and downcasting is essential for 1ULP evaluation in float precision:
0393                 auto v = daubechies_scaling_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 0>(grid_refinements);
0394                 return detail::daubechies_eval_type<Real>::vector_cast(v);
0395                 });
0396             // Compute the derivative of the refined grid:
0397             std::future<std::vector<Real>> t1 = std::async(std::launch::async, [&grid_refinements]() {
0398                 auto v = daubechies_scaling_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 1>(grid_refinements);
0399                 return detail::daubechies_eval_type<Real>::vector_cast(v);
0400                 });
0401
0402             // if necessary, compute the second and third derivative:
0403             std::vector<Real> d2ydx2;
0404             std::vector<Real> d3ydx3;
0405             if constexpr (p >= 6) {
0406                 std::future<std::vector<Real>> t3 = std::async(std::launch::async, [&grid_refinements]() {
0407                 auto v = daubechies_scaling_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 2>(grid_refinements);
0408                 return detail::daubechies_eval_type<Real>::vector_cast(v);
0409                 });
0410
0411                 if constexpr (p >= 10) {
0412                 std::future<std::vector<Real>> t4 = std::async(std::launch::async, [&grid_refinements]() {
0413                     auto v = daubechies_scaling_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 3>(grid_refinements);
0414                     return detail::daubechies_eval_type<Real>::vector_cast(v);
0415                     });
0416                 d3ydx3 = t4.get();
0417                 }
0418                 d2ydx2 = t3.get();
0419             }
0420
0421
0422             auto y = t0.get();
0423             auto dydx = t1.get();
0424
0425             if constexpr (p >= 2)
0426             {
0427                 vector_type data(y.size());
0428                 for (size_t i = 0; i < y.size(); ++i)
0429                 {
0430                     data[i][0] = y[i];
0431                     data[i][1] = dydx[i];
0432                     if constexpr (p >= 6)
0433                         data[i][2] = d2ydx2[i];
0434                     if constexpr (p >= 10)
0435                         data[i][3] = d3ydx3[i];
0436                 }
0437                 if constexpr (p <= 3)
0438                     m_interpolator = std::make_shared<interpolator_type>(std::move(data), grid_refinements, Real(0));
0439                 else
0440                     m_interpolator = std::make_shared<interpolator_type>(std::move(data), Real(0), Real(1) / (1 << grid_refinements));
0441             }
0442             else
0443                 m_interpolator = std::make_shared<detail::null_interpolator>();
0444         }
0445     }
0446
0447     inline Real operator()(Real x) const
0448     {
0449         if (x <= 0 || x >= 2*p-1)
0450         {
0451             return 0;
0452         }
0453         return (*m_interpolator)(x);
0454     }
0455
0456     inline Real prime(Real x) const
0457     {
0458         static_assert(p > 2, "The 3-vanishing moment Daubechies scaling function is the first which is continuously differentiable.");
0459         if (x <= Real(0) || x >= 2*p-1)
0460         {
0461             return 0;
0462         }
0463         return m_interpolator->prime(x);
0464     }
0465
0466     inline Real double_prime(Real x) const
0467     {
0468         static_assert(p >= 6, "Second derivatives require at least 6 vanishing moments.");
0469         if (x <= 0 || x >= 2*p - 1)
0470         {
0471             return Real(0);
0472         }
0473         return m_interpolator->double_prime(x);
0474     }
0475
0476     std::pair<Real, Real> support() const
0477     {
0478         return {Real(0), Real(2*p-1)};
0479     }
0480
0481     int64_t bytes() const
0482     {
0483         return m_interpolator->bytes() + sizeof(m_interpolator);
0484     }
0485
0486 private:
0487    std::shared_ptr<interpolator_type> m_interpolator;
0488 };
0489
0490 }
0491 #endif