1[/ 2Copyright (c) 2020 Nick Thompson 3Use, modification and distribution are subject to the 4Boost Software License, Version 1.0. (See accompanying file 5LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 6] 7 8[section:quintic_hermite Quintic Hermite interpolation] 9 10[heading Synopsis] 11`` 12#include <boost/math/interpolators/quintic_hermite.hpp> 13 14namespace boost::math::interpolators { 15 16template<class RandomAccessContainer> 17class quintic_hermite { 18public: 19 using Real = typename RandomAccessContainer::value_type; 20 quintic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2) 21 22 inline Real operator()(Real x) const; 23 24 inline Real prime(Real x) const; 25 26 inline Real double_prime(Real x) const; 27 28 std::pair<Real, Real> domain() const; 29 30 friend std::ostream& operator<<(std::ostream & os, const quintic_hermite & m); 31 32 void push_back(Real x, Real y, Real dydx, Real d2ydx2); 33}; 34 35template<class RandomAccessContainer> 36class cardinal_quintic_hermite { 37public: 38 using Real = typename RandomAccessContainer::value_type; 39 cardinal_quintic_hermite(RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, Real x0, Real dx); 40 41 inline Real operator()(Real x) const; 42 43 inline Real prime(Real x) const; 44 45 inline Real double_prime(Real x) const; 46 47 std::pair<Real, Real> domain() const; 48}; 49 50template<class RandomAccessContainer> 51class cardinal_quintic_hermite_aos { 52public: 53 using Point = typename RandomAccessContainer::value_type; 54 using Real = typename Point::value_type; 55 cardinal_quintic_hermite_aos(RandomAccessContainer && data, Real x0, Real dx) 56 57 inline Real operator()(Real x) const; 58 59 inline Real prime(Real x) const; 60 61 inline Real double_prime(Real x) const; 62 63 std::pair<Real, Real> domain() const; 64 65} 66`` 67 68 69[heading Quintic Hermite Interpolation] 70 71The quintic Hermite interpolator takes a list of possibly non-uniformly spaced abscissas, ordinates, and their velocities and accelerations which are used to construct a quintic interpolating polynomial between segments. 72This is useful for taking solution skeletons from ODE steppers and turning them into a continuous function, provided that the right-hand side /f/(/x/, /y/) is differentiable along the solution path. 73The interpolant is /C/[super 2] and its evaluation has [bigo](log(/N/)) complexity. 74An example usage is as follows: 75 76 std::vector<double> x{1,2,3, 4, 5, 6}; 77 std::vector<double> y{7,8,9,10,11,12}; 78 std::vector<double> dydx{1,1,1,1,1,1}; 79 std::vector<double> d2ydx2{0,0,0,0,0,0}; 80 81 using boost::math::interpolators::quintic_hermite; 82 auto spline = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); 83 // evaluate at a point: 84 double z = spline(3.4); 85 // evaluate derivative at a point: 86 double zprime = spline.prime(3.4); 87 88Periodically, it is helpful to see what data the interpolator has. 89This can be achieved via 90 91 std::cout << spline << "\n"; 92 93Note that the interpolator is pimpl'd, so that copying the class is cheap, and hence it can be shared between threads. 94(The call operator and `.prime()` are threadsafe.) 95 96The interpolator can be updated in constant time. 97Hence we can use `boost::circular_buffer` to do real-time interpolation. 98 99 100 #include <boost/circular_buffer.hpp> 101 ... 102 boost::circular_buffer<double> initial_x{1,2,3,4}; 103 boost::circular_buffer<double> initial_y{4,5,6,7}; 104 boost::circular_buffer<double> initial_dydx{1,1,1,1}; 105 boost::circular_buffer<double> initial_d2ydx2{0,0,0,0}; 106 auto circular_akima = quintic_hermite(std::move(initial_x), std::move(initial_y), std::move(initial_dydx), std::move(initial_d2ydx2)); 107 // interpolate via call operation: 108 double y = circular_akima(3.5); 109 // add new data: 110 circular_akima.push_back(5, 8, 1, 0); 111 // interpolate at 4.5: 112 y = circular_akima(4.5); 113 114[$../graphs/quintic_sine_approximation.svg] 115 116For equispaced data, we can use `cardinal_quintic_hermite` or `cardinal_quintic_hermite_aos` to get constant-time evaluation. 117This is useful in memory-constrained or performance critical applications where data is equispaced. 118 119[heading Complexity and Performance] 120 121The following google benchmark demonstrates the cost of the call operator for this interpolator: 122 123``` 124Run on (16 X 4300 MHz CPU s) 125CPU Caches: 126 L1 Data 32K (x8) 127 L1 Instruction 32K (x8) 128 L2 Unified 1024K (x8) 129 L3 Unified 11264K (x1) 130Load Average: 0.92, 0.64, 0.35 131-------------------------------------------------- 132Benchmark Time 133-------------------------------------------------- 134QuinticHermite<double>/8 8.14 ns 135QuinticHermite<double>/16 8.69 ns 136QuinticHermite<double>/32 9.42 ns 137QuinticHermite<double>/64 9.90 ns 138QuinticHermite<double>/128 10.4 ns 139QuinticHermite<double>/256 10.9 ns 140QuinticHermite<double>/512 11.6 ns 141QuinticHermite<double>/1024 12.3 ns 142QuinticHermite<double>/2048 12.8 ns 143QuinticHermite<double>/4096 13.6 ns 144QuinticHermite<double>/8192 14.6 ns 145QuinticHermite<double>/16384 15.5 ns 146QuinticHermite<double>/32768 17.4 ns 147QuinticHermite<double>/65536 18.5 ns 148QuinticHermite<double>/131072 18.8 ns 149QuinticHermite<double>/262144 19.8 ns 150QuinticHermite<double>/524288 20.5 ns 151QuinticHermite<double>/1048576 21.6 ns 152QuinticHermite<double>/2097152 22.5 ns 153QuinticHermite<double>/4194304 24.2 ns 154QuinticHermite<double>/8388608 26.6 ns 155QuinticHermite<double>/16777216 26.7 ns 156QuinticHermite<double>_BigO 1.14 lgN 157CardinalQuinticHermite<double>/256 5.22 ns 158CardinalQuinticHermite<double>/512 5.21 ns 159CardinalQuinticHermite<double>/1024 5.21 ns 160CardinalQuinticHermite<double>/2048 5.32 ns 161CardinalQuinticHermite<double>/4096 5.33 ns 162CardinalQuinticHermite<double>/8192 5.50 ns 163CardinalQuinticHermite<double>/16384 5.74 ns 164CardinalQuinticHermite<double>/32768 7.74 ns 165CardinalQuinticHermite<double>/65536 10.6 ns 166CardinalQuinticHermite<double>/131072 10.7 ns 167CardinalQuinticHermite<double>/262144 10.6 ns 168CardinalQuinticHermite<double>/524288 10.5 ns 169CardinalQuinticHermite<double>/1048576 10.6 ns 170CardinalQuinticHermite<double>_BigO 7.57 (1) 171CardinalQuinticHermiteAOS<double>/256 5.27 ns 172CardinalQuinticHermiteAOS<double>/512 5.26 ns 173CardinalQuinticHermiteAOS<double>/1024 5.26 ns 174CardinalQuinticHermiteAOS<double>/2048 5.28 ns 175CardinalQuinticHermiteAOS<double>/4096 5.30 ns 176CardinalQuinticHermiteAOS<double>/8192 5.41 ns 177CardinalQuinticHermiteAOS<double>/16384 5.89 ns 178CardinalQuinticHermiteAOS<double>/32768 5.97 ns 179CardinalQuinticHermiteAOS<double>/65536 5.96 ns 180CardinalQuinticHermiteAOS<double>/131072 5.92 ns 181CardinalQuinticHermiteAOS<double>/262144 5.94 ns 182CardinalQuinticHermiteAOS<double>/524288 5.96 ns 183CardinalQuinticHermiteAOS<double>/1048576 5.93 ns 184CardinalQuinticHermiteAOS<double>_BigO 5.64 (1) 185``` 186 187 188[endsect] 189[/section:quintic_hermite] 190
