ceph/src/boost/boost/math/special_functions/daubechies_scaling.hpp

   1 /*
   2  * Copyright Nick Thompson, John Maddock 2020
   3  * Use, modification and distribution are subject to the
   4  * Boost Software License, Version 1.0. (See accompanying file
   5  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7
   8 #ifndef BOOST_MATH_SPECIAL_DAUBECHIES_SCALING_HPP
   9 #define BOOST_MATH_SPECIAL_DAUBECHIES_SCALING_HPP
  10 #include <vector>
  11 #include <array>
  12 #include <cmath>
  13 #include <thread>
  14 #include <future>
  15 #include <iostream>
  16 #include <boost/math/special_functions/detail/daubechies_scaling_integer_grid.hpp>
  17 #include <boost/math/filters/daubechies.hpp>
  18 #include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
  19 #include <boost/math/interpolators/detail/quintic_hermite_detail.hpp>
  20 #include <boost/math/interpolators/detail/septic_hermite_detail.hpp>
  21
  22 namespace boost::math {
  23
  24 template<class Real, int p, int order>
  25 std::vector<Real> daubechies_scaling_dyadic_grid(int64_t j_max)
  26 {
  27     using std::isnan;
  28     using std::sqrt;
  29     auto c = boost::math::filters::daubechies_scaling_filter<Real, p>();
  30     Real scale = sqrt(static_cast<Real>(2))*(1 << order);
  31     for (auto & x : c)
  32     {
  33         x *= scale;
  34     }
  35
  36     auto phik = detail::daubechies_scaling_integer_grid<Real, p, order>();
  37
  38     // Maximum sensible j for 32 bit floats is j_max = 22:
  39     if (std::is_same_v<Real, float>)
  40     {
  41         if (j_max > 23)
  42         {
  43             throw std::logic_error("Requested dyadic grid more dense than number of representables on the interval.");
  44         }
  45     }
  46     std::vector<Real> v(2*p + (2*p-1)*((1<<j_max) -1), std::numeric_limits<Real>::quiet_NaN());
  47     v[0] = 0;
  48     v[v.size()-1] = 0;
  49     for (int64_t i = 0; i < (int64_t) phik.size(); ++i) {
  50         v[i*(1uLL<<j_max)] = phik[i];
  51     }
  52
  53     for (int64_t j = 1; j <= j_max; ++j)
  54     {
  55         int64_t k_max = v.size()/(int64_t(1) << (j_max-j));
  56         for (int64_t k = 1; k < k_max;  k += 2)
  57         {
  58             // Where this value will go:
  59             int64_t delivery_idx = k*(1uLL << (j_max-j));
  60             // This is a nice check, but we've tested this exhaustively, and it's an expensive check:
  61             //if (delivery_idx >= (int64_t) v.size()) {
  62             //    std::cerr << "Delivery index out of range!\n";
  63             //    continue;
  64             //}
  65             Real term = 0;
  66             for (int64_t l = 0; l < (int64_t) c.size(); ++l)
  67             {
  68                 int64_t idx = k*(int64_t(1) << (j_max - j + 1)) - l*(int64_t(1) << j_max);
  69                 if (idx < 0)
  70                 {
  71                     break;
  72                 }
  73                 if (idx < (int64_t) v.size())
  74                 {
  75                     term += c[l]*v[idx];
  76                 }
  77             }
  78             // Again, another nice check:
  79             //if (!isnan(v[delivery_idx])) {
  80             //    std::cerr << "Delivery index already populated!, = " << v[delivery_idx] << "\n";
  81             //    std::cerr << "would overwrite with " << term << "\n";
  82             //}
  83             v[delivery_idx] = term;
  84         }
  85     }
  86     return v;
  87 }
  88
  89 namespace detail {
  90
  91 template<class RandomAccessContainer>
  92 class matched_holder {
  93 public:
  94     using Real = typename RandomAccessContainer::value_type;
  95
  96     matched_holder(RandomAccessContainer && y, RandomAccessContainer && dydx, int grid_refinements, Real x0) : x0_{x0}, y_{std::move(y)}, dy_{std::move(dydx)}
  97     {
  98         inv_h_ = (1 << grid_refinements);
  99         Real h = 1/inv_h_;
 100         for (auto & dy : dy_)
 101         {
 102             dy *= h;
 103         }
 104     }
 105
 106     inline Real operator()(Real x) const
 107     {
 108         using std::floor;
 109         using std::sqrt;
 110         // This is the exact Holder exponent, but it's pessimistic almost everywhere!
 111         // It's only exactly right at dyadic rationals.
 112         //Real const alpha = 2 - log(1+sqrt(Real(3)))/log(Real(2));
 113         // We're gonna use alpha = 1/2, rather than 0.5500...
 114         Real s = (x-x0_)*inv_h_;
 115         Real ii = floor(s);
 116         auto i = static_cast<decltype(y_.size())>(ii);
 117         Real t = s - ii;
 118         Real dphi = dy_[i+1];
 119         Real diff = y_[i+1] - y_[i];
 120         return y_[i] + (2*dphi - diff)*t + 2*sqrt(t)*(diff-dphi);
 121     }
 122
 123     int64_t bytes() const
 124     {
 125         return 2*y_.size()*sizeof(Real) + sizeof(this);
 126     }
 127
 128 private:
 129     Real x0_;
 130     Real inv_h_;
 131     RandomAccessContainer y_;
 132     RandomAccessContainer dy_;
 133 };
 134
 135 template<class RandomAccessContainer>
 136 class matched_holder_aos {
 137 public:
 138     using Point = typename RandomAccessContainer::value_type;
 139     using Real = typename Point::value_type;
 140
 141     matched_holder_aos(RandomAccessContainer && data, int grid_refinements, Real x0) : x0_{x0}, data_{std::move(data)}
 142     {
 143         inv_h_ = Real(1uLL << grid_refinements);
 144         Real h = 1/inv_h_;
 145         for (auto & datum : data_)
 146         {
 147             datum[1] *= h;
 148         }
 149     }
 150
 151     inline Real operator()(Real x) const
 152     {
 153         using std::floor;
 154         using std::sqrt;
 155         Real s = (x-x0_)*inv_h_;
 156         Real ii = floor(s);
 157         auto i = static_cast<decltype(data_.size())>(ii);
 158         Real t = s - ii;
 159         Real y0 = data_[i][0];
 160         Real y1 = data_[i+1][0];
 161         Real dphi = data_[i+1][1];
 162         Real diff = y1 - y0;
 163         return y0 + (2*dphi - diff)*t + 2*sqrt(t)*(diff-dphi);
 164     }
 165
 166     int64_t bytes() const
 167     {
 168         return data_.size()*data_[0].size()*sizeof(Real) + sizeof(this);
 169     }
 170
 171 private:
 172     Real x0_;
 173     Real inv_h_;
 174     RandomAccessContainer data_;
 175 };
 176
 177
 178 template<class RandomAccessContainer>
 179 class linear_interpolation {
 180 public:
 181     using Real = typename RandomAccessContainer::value_type;
 182
 183     linear_interpolation(RandomAccessContainer && y, RandomAccessContainer && dydx, int grid_refinements) : y_{std::move(y)}, dydx_{std::move(dydx)}
 184     {
 185         s_ = (1 << grid_refinements);
 186     }
 187
 188     inline Real operator()(Real x) const
 189     {
 190         using std::floor;
 191         Real y = x*s_;
 192         Real k = floor(y);
 193
 194         int64_t kk = static_cast<int64_t>(k);
 195         Real t = y - k;
 196         return (1-t)*y_[kk] + t*y_[kk+1];
 197     }
 198
 199     inline Real prime(Real x) const
 200     {
 201         using std::floor;
 202         Real y = x*s_;
 203         Real k = floor(y);
 204
 205         int64_t kk = static_cast<int64_t>(k);
 206         Real t = y - k;
 207         return (1-t)*dydx_[kk] + t*dydx_[kk+1];
 208     }
 209
 210     int64_t bytes() const
 211     {
 212         return (1 + y_.size() + dydx_.size())*sizeof(Real) + sizeof(y_) + sizeof(dydx_);
 213     }
 214
 215 private:
 216     Real s_;
 217     RandomAccessContainer y_;
 218     RandomAccessContainer dydx_;
 219 };
 220
 221 template<class RandomAccessContainer>
 222 class linear_interpolation_aos {
 223 public:
 224     using Point = typename RandomAccessContainer::value_type;
 225     using Real = typename Point::value_type;
 226
 227     linear_interpolation_aos(RandomAccessContainer && data, int grid_refinements, Real x0) : x0_{x0}, data_{std::move(data)}
 228     {
 229         s_ = Real(1uLL << grid_refinements);
 230     }
 231
 232     inline Real operator()(Real x) const
 233     {
 234         using std::floor;
 235         Real y = (x-x0_)*s_;
 236         Real k = floor(y);
 237
 238         int64_t kk = static_cast<int64_t>(k);
 239         Real t = y - k;
 240         return (t != 0) ? (1-t)*data_[kk][0] + t*data_[kk+1][0] : data_[kk][0];
 241     }
 242
 243     inline Real prime(Real x) const
 244     {
 245         using std::floor;
 246         Real y = (x-x0_)*s_;
 247         Real k = floor(y);
 248
 249         int64_t kk = static_cast<int64_t>(k);
 250         Real t = y - k;
 251         return t != 0 ? (1-t)*data_[kk][1] + t*data_[kk+1][1] : data_[kk][1];
 252     }
 253
 254     int64_t bytes() const
 255     {
 256         return sizeof(this) + data_.size()*data_[0].size()*sizeof(Real);
 257     }
 258
 259 private:
 260     Real x0_;
 261     Real s_;
 262     RandomAccessContainer data_;
 263 };
 264
 265
 266 template <class T>
 267 struct daubechies_eval_type
 268 {
 269    typedef T type;
 270
 271    static const std::vector<T>& vector_cast(const std::vector<T>& v) { return v; }
 272
 273 };
 274 template <>
 275 struct daubechies_eval_type<float>
 276 {
 277    typedef double type;
 278
 279    inline static std::vector<float> vector_cast(const std::vector<double>& v)
 280    {
 281       std::vector<float> result(v.size());
 282       for (unsigned i = 0; i < v.size(); ++i)
 283          result[i] = static_cast<float>(v[i]);
 284       return result;
 285    }
 286 };
 287 template <>
 288 struct daubechies_eval_type<double>
 289 {
 290    typedef long double type;
 291
 292    inline static std::vector<double> vector_cast(const std::vector<long double>& v)
 293    {
 294       std::vector<double> result(v.size());
 295       for (unsigned i = 0; i < v.size(); ++i)
 296          result[i] = static_cast<double>(v[i]);
 297       return result;
 298    }
 299 };
 300
 301 struct null_interpolator
 302 {
 303    template <class T>
 304    T operator()(const T&)
 305    {
 306       return 1;
 307    }
 308 };
 309
 310 } // namespace detail
 311
 312 template<class Real, int p>
 313 class daubechies_scaling {
 314    //
 315    // Some type manipulation so we know the type of the interpolator, and the vector type it requires:
 316    //
 317    typedef std::vector<std::array<Real, p < 6 ? 2 : p < 10 ? 3 : 4>> vector_type;
 318    //
 319    // List our interpolators:
 320    //
 321    typedef std::tuple<
 322       detail::null_interpolator, detail::matched_holder_aos<vector_type>, detail::linear_interpolation_aos<vector_type>,
 323       interpolators::detail::cardinal_cubic_hermite_detail_aos<vector_type>, interpolators::detail::cardinal_quintic_hermite_detail_aos<vector_type>,
 324       interpolators::detail::cardinal_septic_hermite_detail_aos<vector_type> > interpolator_list;
 325    //
 326    // Select the one we need:
 327    //
 328    typedef std::tuple_element_t<
 329       p == 1 ? 0 :
 330       p == 2 ? 1 :
 331       p == 3 ? 2 :
 332       p <= 5 ? 3 :
 333       p <= 9 ? 4 : 5, interpolator_list> interpolator_type;
 334
 335 public:
 336    daubechies_scaling(int grid_refinements = -1)
 337    {
 338       static_assert(p < 20, "Daubechies scaling functions are only implemented for p < 20.");
 339       static_assert(p > 0, "Daubechies scaling functions must have at least 1 vanishing moment.");
 340       if constexpr (p == 1)
 341       {
 342          return;
 343       }
 344       else {
 345          if (grid_refinements < 0)
 346          {
 347             if (std::is_same_v<Real, float>)
 348             {
 349                if (grid_refinements == -2)
 350                {
 351                   // Control absolute error:
 352                   //                          p= 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
 353                   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 };
 354                   grid_refinements = r[p];
 355                }
 356                else
 357                {
 358                   // Control relative error:
 359                   //                          p= 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
 360                   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 };
 361                   grid_refinements = r[p];
 362                }
 363             }
 364             else if (std::is_same_v<Real, double>)
 365             {
 366                //                          p= 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
 367                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 };
 368                grid_refinements = r[p];
 369             }
 370             else
 371             {
 372                grid_refinements = 21;
 373             }
 374          }
 375
 376          // Compute the refined grid:
 377          // 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.
 378          std::future<std::vector<Real>> t0 = std::async(std::launch::async, [&grid_refinements]() {
 379             // Computing in higher precision and downcasting is essential for 1ULP evaluation in float precision:
 380             auto v = daubechies_scaling_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 0>(grid_refinements);
 381             return detail::daubechies_eval_type<Real>::vector_cast(v);
 382             });
 383          // Compute the derivative of the refined grid:
 384          std::future<std::vector<Real>> t1 = std::async(std::launch::async, [&grid_refinements]() {
 385             auto v = daubechies_scaling_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 1>(grid_refinements);
 386             return detail::daubechies_eval_type<Real>::vector_cast(v);
 387             });
 388
 389          // if necessary, compute the second and third derivative:
 390          std::vector<Real> d2ydx2;
 391          std::vector<Real> d3ydx3;
 392          if constexpr (p >= 6) {
 393             std::future<std::vector<Real>> t3 = std::async(std::launch::async, [&grid_refinements]() {
 394                auto v = daubechies_scaling_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 2>(grid_refinements);
 395                return detail::daubechies_eval_type<Real>::vector_cast(v);
 396                });
 397
 398             if constexpr (p >= 10) {
 399                std::future<std::vector<Real>> t4 = std::async(std::launch::async, [&grid_refinements]() {
 400                   auto v = daubechies_scaling_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 3>(grid_refinements);
 401                   return detail::daubechies_eval_type<Real>::vector_cast(v);
 402                   });
 403                d3ydx3 = t4.get();
 404             }
 405             d2ydx2 = t3.get();
 406          }
 407
 408
 409          auto y = t0.get();
 410          auto dydx = t1.get();
 411
 412          if constexpr (p >= 2)
 413          {
 414             vector_type data(y.size());
 415             for (size_t i = 0; i < y.size(); ++i)
 416             {
 417                data[i][0] = y[i];
 418                data[i][1] = dydx[i];
 419                if constexpr (p >= 6)
 420                   data[i][2] = d2ydx2[i];
 421                if constexpr (p >= 10)
 422                   data[i][3] = d3ydx3[i];
 423             }
 424             if constexpr (p <= 3)
 425                m_interpolator = std::make_shared<interpolator_type>(std::move(data), grid_refinements, Real(0));
 426             else
 427                m_interpolator = std::make_shared<interpolator_type>(std::move(data), Real(0), Real(1) / (1 << grid_refinements));
 428          }
 429          else
 430             m_interpolator = std::make_shared<detail::null_interpolator>();
 431       }
 432    }
 433
 434     inline Real operator()(Real x) const
 435     {
 436         if (x <= 0 || x >= 2*p-1)
 437         {
 438             return 0;
 439         }
 440         return (*m_interpolator)(x);
 441     }
 442
 443     inline Real prime(Real x) const
 444     {
 445         static_assert(p > 2, "The 3-vanishing moment Daubechies scaling function is the first which is continuously differentiable.");
 446         if (x <= 0 || x >= 2*p-1)
 447         {
 448             return 0;
 449         }
 450         return m_interpolator->prime(x);
 451     }
 452
 453     inline Real double_prime(Real x) const
 454     {
 455         static_assert(p >= 6, "Second derivatives require at least 6 vanishing moments.");
 456         if (x <= 0 || x >= 2*p - 1)
 457         {
 458             return Real(0);
 459         }
 460         return m_interpolator->double_prime(x);
 461     }
 462
 463     std::pair<Real, Real> support() const
 464     {
 465         return {Real(0), Real(2*p-1)};
 466     }
 467
 468     int64_t bytes() const
 469     {
 470        return m_interpolator->bytes() + sizeof(m_interpolator);
 471     }
 472
 473 private:
 474    std::shared_ptr<interpolator_type> m_interpolator;
 475 };
 476
 477 }
 478 #endif