From 09213a0835f94a7945b0afcabffe2a6d3008fdfd Mon Sep 17 00:00:00 2001 From: Nick Date: Sat, 8 Feb 2020 12:46:28 -0500 Subject: [PATCH] Add septic Hermite interpolation. This still needs some tlc to get to 1ULP evaluation [CI SKIP] --- example/find_best_daubechies_interpolator.cpp | 103 +++++------ .../detail/septic_hermite_detail.hpp | 160 ++++++++++++++++++ .../math/interpolators/septic_hermite.hpp | 38 +++++ .../special_functions/daubechies_scaling.hpp | 50 +++++- 4 files changed, 282 insertions(+), 69 deletions(-) create mode 100644 include/boost/math/interpolators/detail/septic_hermite_detail.hpp create mode 100644 include/boost/math/interpolators/septic_hermite.hpp diff --git a/example/find_best_daubechies_interpolator.cpp b/example/find_best_daubechies_interpolator.cpp index 9ca360dd5..d9277f698 100644 --- a/example/find_best_daubechies_interpolator.cpp +++ b/example/find_best_daubechies_interpolator.cpp @@ -5,8 +5,9 @@ #include #include #include -#include #include +#include +#include #include #include #include @@ -22,64 +23,6 @@ using boost::multiprecision::float128; -/*template -void do_ulp() -{ - std::cout << "Creating ULP plot on type " << boost::core::demangle(typeid(Real).name()) << " and " << p << " vanishing moments.\n"; - using std::floor; - using std::ceil; - using std::abs; - int rmax = 14; - std::cout << "Computing phi_dense\n"; - std::future> f1 = std::async(std::launch::async, [&]{ return boost::math::detail::dyadic_grid(rmax); }); - std::future> f2 = std::async(std::launch::async, [&]{ return boost::math::detail::dyadic_grid(rmax); }); - - std::cout << "Computing phi and phi_prime\n"; - std::future> f3 = std::async(std::launch::async, [&]{ return boost::math::detail::dyadic_grid(rmax-2); }); - std::future> f4 = std::async(std::launch::async, [&]{ return boost::math::detail::dyadic_grid(rmax-2); }); - auto phi_dense = f1.get(); - auto phi_dense_prime = f2.get(); - auto phi_accurate = f3.get(); - auto phi_prime_accurate = f4.get(); - - Real dx_dense = (2*p-1)/static_cast(phi_dense.size()-1); - std::cout << "Done precomputing grids; downcasting now.\n"; - std::vector phi(phi_accurate.size()); - for (size_t i = 0; i < phi_accurate.size(); ++i) { - phi[i] = Real(phi_accurate[i]); - } - std::vector phi_prime(phi_accurate.size()); - for (size_t i = 0; i < phi_prime_accurate.size(); ++i) { - phi_prime[i] = Real(phi_prime_accurate[i]); - } - phi_accurate.resize(0); - phi_prime_accurate.resize(0); - - std::vector x(phi.size()); - Real dx = (2*p-1)/static_cast(x.size()-1); - std::cout << "dx = " << dx << "\n"; - for (size_t i = 0; i < x.size(); ++i) { - x[i] = i*dx; - } - - auto ch = boost::math::interpolators::cubic_hermite(std::move(x), std::move(phi), std::move(phi_prime)); - - std::vector x_acc(phi_dense.size()); - - for (size_t i = 0; i < x_acc.size(); ++i) { - x_acc[i] = i*dx_dense; - } - - auto acc = boost::math::interpolators::cubic_hermite(std::move(x_acc), std::move(phi_dense), std::move(phi_dense_prime)); - std::cout << "Writing ulp plot\n"; - std::string title = "daub" + std::to_string(p) + "_" + std::to_string(rmax-2) + "_" + boost::core::demangle(typeid(Real).name()) + ".svg"; - quicksvg::ulp_plot(ch, acc, Real(0), Real(2*p-1), - "ULP plot of Daubechies with " + std::to_string(rmax-2) + " refinements on type " + boost::core::demangle(typeid(Real).name()), - title, 15000, 1100, 10); - std::cout << "Done writing ulp plot\n"; - -}*/ - template void choose_refinement() @@ -486,11 +429,29 @@ void find_best_interpolator() m.insert({sup, "Third-order Taylor"}); } + { + auto phi_copy = phi; + auto x_copy = x; + auto phi_prime_copy = phi_prime; + auto phi_dbl_prime_copy = phi_dbl_prime; + auto phi_triple_prime_copy = phi_triple_prime; + auto sh = boost::math::interpolators::septic_hermite(std::move(x_copy), std::move(phi_copy), std::move(phi_prime_copy), std::move(phi_dbl_prime_copy), std::move(phi_triple_prime_copy)); + Real sup = 0; + for (size_t i = 0; i < phi_dense.size(); ++i) { + Real x = i*dx_dense; + Real diff = abs(phi_dense[i] - sh(x)); + if (diff > sup) { + sup = diff; + } + } + m.insert({sup, "septic_hermite_spline"}); + } + } std::string best = "none"; Real best_sup = 1000000000; - std::cout << std::setprecision(20) << std::fixed; + std::cout << std::setprecision(std::numeric_limits::digits10 + 3) << std::fixed; for (auto & e : m) { std::cout << "\t" << e.first << " is error of " << e.second << "\n"; @@ -505,12 +466,10 @@ void find_best_interpolator() int main() { - //do_ulp(); - //choose_refinement(); //choose_refinement(); // Says linear interpolation is the best: - find_best_interpolator(); + /*find_best_interpolator(); // Says linear interpolation is the best: find_best_interpolator(); @@ -528,5 +487,21 @@ int main() { find_best_interpolator(); // Says quintic_hermite_spline is best: - find_best_interpolator(); + find_best_interpolator();*/ + // Says quintic_hermite_spline is best: + find_best_interpolator(); + // Says septic_hermite_spline is best: + find_best_interpolator(); + /* + // Says septic_hermite_spline is best: + find_best_interpolator(); + // Says septic_hermite_spline is best: + find_best_interpolator(); + // Says septic_hermite_spline is best: + find_best_interpolator(); + // Says septic_hermite_spline is best: + find_best_interpolator();*/ + + // Says septic_hermite_spline is best: + find_best_interpolator(); } diff --git a/include/boost/math/interpolators/detail/septic_hermite_detail.hpp b/include/boost/math/interpolators/detail/septic_hermite_detail.hpp new file mode 100644 index 000000000..c3e4b8e0b --- /dev/null +++ b/include/boost/math/interpolators/detail/septic_hermite_detail.hpp @@ -0,0 +1,160 @@ +#ifndef BOOST_MATH_INTERPOLATORS_DETAIL_SEPTIC_HERMITE_DETAIL_HPP +#define BOOST_MATH_INTERPOLATORS_DETAIL_SEPTIC_HERMITE_DETAIL_HPP +#include +#include +#include +#include + +namespace boost::math::interpolators::detail { + +template +class septic_hermite_detail { +public: + using Real = typename RandomAccessContainer::value_type; + septic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3) + : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}, d2ydx2_{std::move(d2ydx2)}, d3ydx3_{std::move(d3ydx3)} + { + if (x_.size() != y_.size()) { + throw std::domain_error("Number of abscissas must = number of ordinates."); + } + if (x_.size() != dydx_.size()) { + throw std::domain_error("Numbers of derivatives must = number of abscissas."); + } + if (x_.size() != d2ydx2_.size()) { + throw std::domain_error("Number of second derivatives must equal number of abscissas."); + } + if (x_.size() != d3ydx3_.size()) { + throw std::domain_error("Number of third derivatives must equal number of abscissas."); + } + + if (x_.size() < 2) { + throw std::domain_error("At least 2 abscissas are required."); + } + Real x0 = x_[0]; + for (decltype(x_.size()) i = 1; i < x_.size(); ++i) { + Real x1 = x_[i]; + if (x1 <= x0) + { + throw std::domain_error("Abscissas must be sorted in strictly increasing order x0 < x1 < ... < x_{n-1}"); + } + x0 = x1; + } + } + + void push_back(Real x, Real y, Real dydx, Real d2ydx2, Real d3ydx3) { + using std::abs; + using std::isnan; + if (x <= x_.back()) { + throw std::domain_error("Calling push_back must preserve the monotonicity of the x's"); + } + x_.push_back(x); + y_.push_back(y); + dydx_.push_back(dydx); + d2ydx2_.push_back(d2ydx2); + d3ydx3_.push_back(d3ydx3); + } + + Real operator()(Real x) const { + if (x < x_[0] || x > x_.back()) { + std::ostringstream oss; + oss.precision(std::numeric_limits::digits10+3); + oss << "Requested abscissa x = " << x << ", which is outside of allowed range [" + << x_[0] << ", " << x_.back() << "]"; + throw std::domain_error(oss.str()); + } + // We need t := (x-x_k)/(x_{k+1}-x_k) \in [0,1) for this to work. + // Sadly this neccessitates this loathesome check, otherwise we get t = 1 at x = xf. + if (x == x_.back()) { + return y_.back(); + } + + auto it = std::upper_bound(x_.begin(), x_.end(), x); + auto i = std::distance(x_.begin(), it) -1; + Real x0 = *(it-1); + Real x1 = *it; + Real y0 = y_[i]; + Real y1 = y_[i+1]; + // Velocity: + Real v0 = dydx_[i]; + Real v1 = dydx_[i+1]; + // Acceleration: + Real a0 = d2ydx2_[i]; + Real a1 = d2ydx2_[i+1]; + // Jerk: + Real j0 = d3ydx3_[i]; + Real j1 = d3ydx3_[i+1]; + + Real dx = (x1-x0); + Real t = (x-x0)/dx; + + // See: + // http://seisweb.usask.ca/classes/GEOL481/2017/Labs/interpolation_utilities_matlab/shermite.m + Real t2 = t*t; + Real t3 = t2*t; + Real t4 = t3*t; + Real t5 = t4*t; + Real t6 = t5*t; + Real t7 = t6*t; + Real dx2 = dx*dx; + Real dx3 = dx2*dx; + Real z0 = 20*t7 - 70*t6 + 84*t5 - 35*t4 + 1; + Real z1 = 10*t7 - 36*t6 + 45*t5 - 20*t4 + t; + Real z2 = 2*t7 - 15*t6/2 + 10*t5 - 5*t4 + t2/2; + Real z3 = t7/6 - 2*t6/3 + t5 - 2*t4/3 + t3/6; + Real z4 = 1 - z0; + Real z5 = 10*t7 - 34*t6 + 39*t5 - 15*t4; + Real z6 = -2*t7 + 13*t6/2 - 7*t5 + 5*t4/2; + Real z7 = t7/6 - t6/2 + t5/2 - t4/6; + + Real y = z0*y0 + z1*dx*v0 + z2*dx2*a0 + z3*dx3*j0 + z4*y1 + z5*dx*v1 + z6*dx2*a1 + z7*dx3*j1; + return y; + } + + Real prime(Real x) const { + if (x < x_[0] || x > x_.back()) { + std::ostringstream oss; + oss.precision(std::numeric_limits::digits10+3); + oss << "Requested abscissa x = " << x << ", which is outside of allowed range [" + << x_[0] << ", " << x_.back() << "]"; + throw std::domain_error(oss.str()); + } + if (x == x_.back()) { + return dydx_.back(); + } + + auto it = std::upper_bound(x_.begin(), x_.end(), x); + auto i = std::distance(x_.begin(), it) -1; + Real x0 = *(it-1); + Real x1 = *it; + Real s0 = dydx_[i]; + Real s1 = dydx_[i+1]; + + // Ridiculous linear interpolation. Fine for now: + Real numerator = s0*(x1-x) + s1*(x-x0); + Real denominator = x1 - x0; + return numerator/denominator; + } + + + friend std::ostream& operator<<(std::ostream & os, const septic_hermite_detail & m) + { + os << "(x,y,y') = {"; + for (size_t i = 0; i < m.x_.size() - 1; ++i) { + os << "(" << m.x_[i] << ", " << m.y_[i] << ", " << m.dydx_[i] << ", " << m.d2ydx2_[i] << ", " << m.d3ydx3_[i] << "), "; + } + auto n = m.x_.size()-1; + os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ", " << m.d2ydx2_[n] << m.d3ydx3_[n] << ")}"; + return os; + } + + + +private: + RandomAccessContainer x_; + RandomAccessContainer y_; + RandomAccessContainer dydx_; + RandomAccessContainer d2ydx2_; + RandomAccessContainer d3ydx3_; +}; +} +#endif \ No newline at end of file diff --git a/include/boost/math/interpolators/septic_hermite.hpp b/include/boost/math/interpolators/septic_hermite.hpp new file mode 100644 index 000000000..c544d1c63 --- /dev/null +++ b/include/boost/math/interpolators/septic_hermite.hpp @@ -0,0 +1,38 @@ +#ifndef BOOST_MATH_INTERPOLATORS_SEPTIC_HERMITE_HPP +#define BOOST_MATH_INTERPOLATORS_SEPTIC_HERMITE_HPP +#include +#include +#include +#include + +namespace boost::math::interpolators { + +template +class septic_hermite { +public: + using Real = typename RandomAccessContainer::value_type; + septic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, + RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3) + : impl_(std::make_shared>(std::move(x), + std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3))) + {} + + Real operator()(Real x) const { + return impl_->operator()(x); + } + + Real prime(Real x) const { + return impl_->prime(x); + } + + friend std::ostream& operator<<(std::ostream & os, const septic_hermite & m) + { + os << *m.impl_; + return os; + } + +private: + std::shared_ptr> impl_; +}; +} +#endif \ No newline at end of file diff --git a/include/boost/math/special_functions/daubechies_scaling.hpp b/include/boost/math/special_functions/daubechies_scaling.hpp index c0e026080..591030adb 100644 --- a/include/boost/math/special_functions/daubechies_scaling.hpp +++ b/include/boost/math/special_functions/daubechies_scaling.hpp @@ -19,6 +19,7 @@ #include #include #include +#include namespace boost::math { @@ -226,6 +227,7 @@ public: // if necessary, compute the second derivative: std::vector d2ydx2; + std::vector d3ydx3; if constexpr (p >= 6) { std::future> t3 = std::async(std::launch::async, [&grid_refinements]() { if constexpr (std::is_same_v) { @@ -245,10 +247,35 @@ public: return w; } - return detail::dyadic_grid(grid_refinements); }); + return detail::dyadic_grid(grid_refinements); + }); + + if constexpr (p >= 10) { + std::future> t4 = std::async(std::launch::async, [&grid_refinements]() { + if constexpr (std::is_same_v) { + auto v = detail::dyadic_grid(grid_refinements); + std::vector w(v.size()); + for (size_t i = 0; i < v.size(); ++i) { + w[i] = v[i]; + } + return w; + } + else if constexpr (std::is_same_v) { + auto v = detail::dyadic_grid(grid_refinements); + std::vector w(v.size()); + for (size_t i = 0; i < v.size(); ++i) { + w[i] = v[i]; + } + return w; + } + + return detail::dyadic_grid(grid_refinements); }); + d3ydx3 = t4.get(); + } d2ydx2 = t3.get(); } + auto y = t0.get(); auto dydx = t1.get(); @@ -270,9 +297,13 @@ public: if constexpr (p == 4 || p == 5) { m_cbh = std::make_shared>>(std::move(x), std::move(y), std::move(dydx)); } - else if constexpr (p >= 6) { + if constexpr (p >= 6 && p <= 9) { m_qh = std::make_shared>>(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); } + if constexpr (p >= 10) { + m_sh = std::make_shared>>(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3)); + } + } @@ -286,9 +317,12 @@ public: if constexpr (p==4 || p ==5) { return m_cbh->operator()(x); } - else if constexpr (p >= 6) { + if constexpr (p >= 6 && p <= 9) { return m_qh->operator()(x); } + if constexpr (p >= 10) { + return m_sh->operator()(x); + } } Real prime(Real x) const { @@ -298,9 +332,12 @@ public: if constexpr (p==4 || p ==5) { return m_cbh->prime(x); } - else if constexpr (p >= 6) { + if constexpr (p >= 6 && p <= 9) { return m_qh->prime(x); } + if constexpr (p >= 10) { + return m_sh->prime(x); + } } std::pair support() const { @@ -315,9 +352,12 @@ private: std::shared_ptr>> m_lin; // need this for p = 4,5: std::shared_ptr>> m_cbh; - // need this for p >= 6: + // need this for p = 6,7,8,9: std::shared_ptr>> m_qh; + // need this for p >= 10: + std::shared_ptr>> m_sh; + /*Real constant_interpolation(Real x) const { if (x <= 0 || x >= 2*p-1) { return Real(0);