]>
Commit | Line | Data |
---|---|---|
f67539c2 TL |
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> | |
f67539c2 TL |
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; | |
20effc67 | 28 | using std::sqrt; |
f67539c2 | 29 | auto c = boost::math::filters::daubechies_scaling_filter<Real, p>(); |
20effc67 | 30 | Real scale = sqrt(static_cast<Real>(2))*(1 << order); |
f67539c2 TL |
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 |