267 lines · plain
1/*2 * Copyright Nick Thompson, 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_WAVELET_HPP9#define BOOST_MATH_SPECIAL_DAUBECHIES_WAVELET_HPP10#include <vector>11#include <array>12#include <cmath>13#include <thread>14#include <future>15#include <iostream>16#include <boost/math/constants/constants.hpp>17#include <boost/math/special_functions/detail/daubechies_scaling_integer_grid.hpp>18#include <boost/math/special_functions/daubechies_scaling.hpp>19#include <boost/math/filters/daubechies.hpp>20#include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>21#include <boost/math/interpolators/detail/quintic_hermite_detail.hpp>22#include <boost/math/interpolators/detail/septic_hermite_detail.hpp>23 24#include <boost/math/tools/is_standalone.hpp>25#ifndef BOOST_MATH_STANDALONE26#include <boost/config.hpp>27#ifdef BOOST_MATH_NO_CXX17_IF_CONSTEXPR28#error "The header <boost/math/norms.hpp> can only be used in C++17 and later."29#endif30#endif31 32namespace boost::math {33 34 template<class Real, int p, int order>35 std::vector<Real> daubechies_wavelet_dyadic_grid(int64_t j_max)36 {37 if (j_max == 0)38 {39 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.");40 }41 auto phijk = daubechies_scaling_dyadic_grid<Real, p, order>(j_max - 1);42 //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})43 //For derivatives just map c_k -> 2^order c_k.44 auto d = boost::math::filters::daubechies_scaling_filter<Real, p>();45 Real scale = boost::math::constants::root_two<Real>() * (1 << order);46 for (size_t i = 0; i < d.size(); ++i)47 {48 d[i] *= scale;49 if (!(i & 1))50 {51 d[i] = -d[i];52 }53 }54 55 std::vector<Real> v(2 * p + (2 * p - 1) * ((int64_t(1) << j_max) - 1), std::numeric_limits<Real>::quiet_NaN());56 v[0] = 0;57 v[v.size() - 1] = 0;58 59 for (int64_t l = 1; l < static_cast<int64_t>(v.size() - 1); ++l)60 {61 Real term = 0;62 for (int64_t k = 0; k < static_cast<int64_t>(d.size()); ++k)63 {64 int64_t idx = (int64_t(1) << (j_max - 1)) * (1 - 2 * p + k) + l;65 if (idx < 0 || idx >= static_cast<int64_t>(phijk.size()))66 {67 continue;68 }69 term += d[k] * phijk[idx];70 }71 v[l] = term;72 }73 74 return v;75 }76 77 78 template<class Real, int p>79 class daubechies_wavelet {80 //81 // Some type manipulation so we know the type of the interpolator, and the vector type it requires:82 //83 using vector_type = std::vector < std::array < Real, p < 6 ? 2 : p < 10 ? 3 : 4>>;84 //85 // List our interpolators:86 //87 using interpolator_list = std::tuple<88 detail::null_interpolator, detail::matched_holder_aos<vector_type>, detail::linear_interpolation_aos<vector_type>,89 interpolators::detail::cardinal_cubic_hermite_detail_aos<vector_type>, interpolators::detail::cardinal_quintic_hermite_detail_aos<vector_type>,90 interpolators::detail::cardinal_septic_hermite_detail_aos<vector_type> > ;91 //92 // Select the one we need:93 //94 using interpolator_type = std::tuple_element_t<95 p == 1 ? 0 :96 p == 2 ? 1 :97 p == 3 ? 2 :98 p <= 5 ? 3 :99 p <= 9 ? 4 : 5, interpolator_list>;100 public:101 explicit daubechies_wavelet(int grid_refinements = -1)102 {103 static_assert(p < 20, "Daubechies wavelets are only implemented for p < 20.");104 static_assert(p > 0, "Daubechies wavelets must have at least 1 vanishing moment.");105 if (grid_refinements == 0)106 {107 throw std::domain_error("The wavelet requires at least 1 grid refinement.");108 }109 if constexpr (p == 1)110 {111 return;112 }113 else114 {115 if (grid_refinements < 0)116 {117 if constexpr (std::is_same_v<Real, float>)118 {119 if (grid_refinements == -2)120 {121 // Control absolute error:122 // p= 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19123 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 };124 grid_refinements = r[p];125 }126 else127 {128 // Control relative error:129 // p= 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19130 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 };131 grid_refinements = r[p];132 }133 }134 else if constexpr (std::is_same_v<Real, double>)135 {136 // p= 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19137 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 };138 grid_refinements = r[p];139 }140 else141 {142 grid_refinements = 21;143 }144 }145 146 // Compute the refined grid:147 // 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.148 std::future<std::vector<Real>> t0 = std::async(std::launch::async, [&grid_refinements]() {149 // Computing in higher precision and downcasting is essential for 1ULP evaluation in float precision:150 auto v = daubechies_wavelet_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 0>(grid_refinements);151 return detail::daubechies_eval_type<Real>::vector_cast(v);152 });153 // Compute the derivative of the refined grid:154 std::future<std::vector<Real>> t1 = std::async(std::launch::async, [&grid_refinements]() {155 auto v = daubechies_wavelet_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 1>(grid_refinements);156 return detail::daubechies_eval_type<Real>::vector_cast(v);157 });158 159 // if necessary, compute the second and third derivative:160 std::vector<Real> d2ydx2;161 std::vector<Real> d3ydx3;162 if constexpr (p >= 6) {163 std::future<std::vector<Real>> t3 = std::async(std::launch::async, [&grid_refinements]() {164 auto v = daubechies_wavelet_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 2>(grid_refinements);165 return detail::daubechies_eval_type<Real>::vector_cast(v);166 });167 168 if constexpr (p >= 10) {169 std::future<std::vector<Real>> t4 = std::async(std::launch::async, [&grid_refinements]() {170 auto v = daubechies_wavelet_dyadic_grid<typename detail::daubechies_eval_type<Real>::type, p, 3>(grid_refinements);171 return detail::daubechies_eval_type<Real>::vector_cast(v);172 });173 d3ydx3 = t4.get();174 }175 d2ydx2 = t3.get();176 }177 178 179 auto y = t0.get();180 auto dydx = t1.get();181 182 if constexpr (p >= 2)183 {184 vector_type data(y.size());185 for (size_t i = 0; i < y.size(); ++i)186 {187 data[i][0] = y[i];188 data[i][1] = dydx[i];189 if constexpr (p >= 6)190 data[i][2] = d2ydx2[i];191 if constexpr (p >= 10)192 data[i][3] = d3ydx3[i];193 }194 if constexpr (p <= 3)195 m_interpolator = std::make_shared<interpolator_type>(std::move(data), grid_refinements, Real(-p + 1));196 else197 m_interpolator = std::make_shared<interpolator_type>(std::move(data), Real(-p + 1), Real(1) / (1 << grid_refinements));198 }199 else200 m_interpolator = std::make_shared<detail::null_interpolator>();201 }202 }203 204 205 inline Real operator()(Real x) const206 {207 if (x <= -p + 1 || x >= p)208 {209 return 0;210 }211 212 if constexpr (p == 1)213 {214 if (x < Real(1) / Real(2))215 {216 return 1;217 }218 else if (x == Real(1) / Real(2))219 {220 return 0;221 }222 return -1;223 }224 else225 {226 return (*m_interpolator)(x);227 }228 }229 230 inline Real prime(Real x) const231 {232 static_assert(p > 2, "The 3-vanishing moment Daubechies wavelet is the first which is continuously differentiable.");233 if (x <= -p + 1 || x >= p)234 {235 return 0;236 }237 return m_interpolator->prime(x);238 }239 240 inline Real double_prime(Real x) const241 {242 static_assert(p >= 6, "Second derivatives of Daubechies wavelets require at least 6 vanishing moments.");243 if (x <= -p + 1 || x >= p)244 {245 return Real(0);246 }247 return m_interpolator->double_prime(x);248 }249 250 std::pair<Real, Real> support() const251 {252 return std::make_pair(Real(-p + 1), Real(p));253 }254 255 int64_t bytes() const256 {257 return m_interpolator->bytes() + sizeof(*this);258 }259 260 private:261 std::shared_ptr<interpolator_type> m_interpolator;262 };263 264}265 266#endif // BOOST_MATH_SPECIAL_DAUBECHIES_WAVELET_HPP267