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:cubic_hermite Cubic Hermite interpolation] 9 10[heading Synopsis] 11`` 12 #include <boost/math/interpolators/cubic_hermite.hpp> 13`` 14 15 namespace boost::math::interpolators { 16 17 template <class RandomAccessContainer> 18 class cubic_hermite 19 { 20 public: 21 22 using Real = RandomAccessContainer::value_type; 23 24 cubic_hermite(RandomAccessContainer&& abscissas, RandomAccessContainer&& ordinates, RandomAccessContainer&& derivatives); 25 26 Real operator()(Real x) const; 27 28 Real prime(Real x) const; 29 30 void push_back(Real x, Real y, Real dydx); 31 32 std::pair<Real, Real> domain() const; 33 34 friend std::ostream& operator<<(std::ostream & os, const cubic_hermite & m); 35 }; 36 37 template<class RandomAccessContainer> 38 class cardinal_cubic_hermite { 39 public: 40 using Real = typename RandomAccessContainer::value_type; 41 42 cardinal_cubic_hermite(RandomAccessContainer && y, RandomAccessContainer && dydx, Real x0, Real dx) 43 44 inline Real operator()(Real x) const; 45 46 inline Real prime(Real x) const; 47 48 std::pair<Real, Real> domain() const; 49 }; 50 51 52 template<class RandomAccessContainer> 53 class cardinal_cubic_hermite_aos { 54 public: 55 using Point = typename RandomAccessContainer::value_type; 56 using Real = typename Point::value_type; 57 58 cardinal_cubic_hermite_aos(RandomAccessContainer && data, Real x0, Real dx); 59 60 inline Real operator()(Real x) const; 61 62 inline Real prime(Real x) const; 63 64 std::pair<Real, Real> domain() const; 65 }; 66 67 } // namespaces 68 69[heading Cubic Hermite Interpolation] 70 71The cubic Hermite interpolant takes non-equispaced data and interpolates between them via cubic Hermite polynomials whose slopes must be provided. 72This is a very nice interpolant for solution skeletons of ODEs steppers, since numerically solving /y/' = /f/(/x/, /y/) produces a list of positions, values, and their derivatives, i.e. (/x/,/y/,/y/')=(/x/,/y/,/f/(/x/,/y/))-which is exactly what is required for the cubic Hermite interpolator. 73The interpolant is /C/[super 1] and evaluation has [bigo](log(/N/)) complexity. 74An example usage is as follows: 75 76 std::vector<double> x{1, 5, 9 , 12}; 77 std::vector<double> y{8,17, 4, -3}; 78 std::vector<double dydx{5, -2, -1}; 79 using boost::math::interpolators::cubic_hermite; 80 auto spline = cubic_hermite(std::move(x), std::move(y), std::move(dydx)); 81 // evaluate at a point: 82 double z = spline(3.4); 83 // evaluate derivative at a point: 84 double zprime = spline.prime(3.4); 85 86Sanity checking of the interpolator can be achieved via: 87 88 std::cout << spline << "\n"; 89 90Note that the interpolator is pimpl'd, so that copying the class is cheap, and hence it can be shared between threads. 91(The call operator and `.prime()` are threadsafe; `push_back` is not.) 92 93This interpolant can be updated in constant time. 94Hence we can use `boost::circular_buffer` to do real-time interpolation: 95 96 #include <boost/circular_buffer.hpp> 97 ... 98 boost::circular_buffer<double> initial_x{1,2,3,4}; 99 boost::circular_buffer<double> initial_y{4,5,6,7}; 100 boost::circular_buffer<double> initial_dydx{1,1,1,1}; 101 auto circular_hermite = cubic_hermite(std::move(initial_x), std::move(initial_y), std::move(initial_dydx)); 102 // interpolate via call operation: 103 double y = circular_hermite(3.5); 104 // add new data: 105 circular_hermite.push_back(5, 8); 106 // interpolate at 4.5: 107 y = circular_hermite(4.5); 108 109For the equispaced case, we can either use `cardinal_cubic_hermite`, which accepts two separate arrays of `y` and `dydx`, or we can use `cardinal_cubic_hermite_aos`, 110which takes a vector of `(y, dydx)`, i.e., and array of structs (`aos`). 111The array of structs should be preferred as it uses cache more effectively. 112 113Usage is as follows: 114 115 using boost::math::interpolators::cardinal_cubic_hermite; 116 double x0 = 0; 117 double dx = 1; 118 std::vector<double> y(128, 1); 119 std::vector<double> dydx(128, 0); 120 auto ch = cardinal_cubic_hermite(std::move(y), std::move(dydx), x0, dx); 121 122For the "array of structs" version: 123 124 using boost::math::interpolators::cardinal_cubic_hermite_aos; 125 std::vector<std::array<double, 2>> data(128); 126 for (auto & datum : data) { 127 datum[0] = 1; // y 128 datum[1] = 0; // y' 129 } 130 131 auto ch = cardinal_cubic_hermite_aos(std::move(data), x0, dx); 132 133For a quick sanity check, we can call `ch.domain()`: 134 135 auto [x_min, x_max] = ch.domain(); 136 137[heading Performance] 138 139Google benchmark was used to evaluate the performance. 140 141``` 142Run on (16 X 4300 MHz CPU s) 143CPU Caches: 144 L1 Data 32 KiB (x8) 145 L1 Instruction 32 KiB (x8) 146 L2 Unified 1024 KiB (x8) 147 L3 Unified 11264 KiB (x1) 148Load Average: 1.16, 1.11, 0.99 149----------------------------------------------------- 150Benchmark Time 151----------------------------------------------------- 152CubicHermite<double>/256 9.67 ns 153CubicHermite<double>/512 10.4 ns 154CubicHermite<double>/1024 10.9 ns 155CubicHermite<double>/2048 12.3 ns 156CubicHermite<double>/4096 13.4 ns 157CubicHermite<double>/8192 14.9 ns 158CubicHermite<double>/16384 15.5 ns 159CubicHermite<double>/32768 18.0 ns 160CubicHermite<double>/65536 19.9 ns 161CubicHermite<double>/131072 23.0 ns 162CubicHermite<double>/262144 25.3 ns 163CubicHermite<double>/524288 28.9 ns 164CubicHermite<double>/1048576 31.0 ns 165CubicHermite<double>_BigO 1.32 lgN 166CubicHermite<double>_RMS 13 % 167CardinalCubicHermite<double>/256 3.14 ns 168CardinalCubicHermite<double>/512 3.13 ns 169CardinalCubicHermite<double>/1024 3.19 ns 170CardinalCubicHermite<double>/2048 3.14 ns 171CardinalCubicHermite<double>/4096 3.38 ns 172CardinalCubicHermite<double>/8192 3.54 ns 173CardinalCubicHermite<double>/16384 3.51 ns 174CardinalCubicHermite<double>/32768 3.76 ns 175CardinalCubicHermite<double>/65536 5.82 ns 176CardinalCubicHermite<double>/131072 5.76 ns 177CardinalCubicHermite<double>/262144 5.85 ns 178CardinalCubicHermite<double>/524288 5.69 ns 179CardinalCubicHermite<double>/1048576 5.77 ns 180CardinalCubicHermite<double>_BigO 4.28 (1) 181CardinalCubicHermite<double>_RMS 28 % 182CardinalCubicHermiteAOS<double>/256 3.22 ns 183CardinalCubicHermiteAOS<double>/512 3.22 ns 184CardinalCubicHermiteAOS<double>/1024 3.26 ns 185CardinalCubicHermiteAOS<double>/2048 3.23 ns 186CardinalCubicHermiteAOS<double>/4096 3.49 ns 187CardinalCubicHermiteAOS<double>/8192 3.57 ns 188CardinalCubicHermiteAOS<double>/16384 3.73 ns 189CardinalCubicHermiteAOS<double>/32768 4.12 ns 190CardinalCubicHermiteAOS<double>/65536 4.13 ns 191CardinalCubicHermiteAOS<double>/131072 4.12 ns 192CardinalCubicHermiteAOS<double>/262144 4.12 ns 193CardinalCubicHermiteAOS<double>/524288 4.12 ns 194CardinalCubicHermiteAOS<double>/1048576 4.19 ns 195CardinalCubicHermiteAOS<double>_BigO 3.73 (1) 196CardinalCubicHermiteAOS<double>_RMS 11 % 197``` 198 199The logarithmic complexity of the non-equispaced version is evident, as is the better cache utilization of the "array of structs" version as the problem size gets larger. 200 201 202[endsect] 203[/section:cubic_hermite] 204