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