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