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