Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <boost/math/interpolators/cardinal_cubic_b_spline.hpp>
8 : #include <cstddef>
9 : #include <optional>
10 : #include <utility>
11 : #include <vector>
12 :
13 : /// \cond
14 : namespace PUP {
15 : class er;
16 : } // namespace PUP
17 : /// \endcond
18 :
19 : namespace intrp {
20 :
21 : /*!
22 : * \ingroup NumericalAlgorithmsGroup
23 : * \brief A cubic B-spline interpolant of uniformly spaced samples
24 : *
25 : * Interpolates samples \f$f_i = f(t_0 + i \Delta t)\f$, \f$i = 0, \ldots,
26 : * N-1\f$, with a cubic B-spline. The interpolation error decreases as
27 : * \f$\mathcal{O}(\Delta t^4)\f$ for smooth data. This class wraps
28 : * `boost::math::interpolators::cardinal_cubic_b_spline` and adds:
29 : *
30 : * - Serialization: the samples, start time, and time step fully determine the
31 : * interpolant, so it can be sent in a `pup` and rebuilt on the receiving
32 : * side.
33 : * - Clamped evaluation: evaluation is clamped to the bounds of the sampled
34 : * interval, so times that fall outside the interval by roundoff evaluate
35 : * to the boundary values instead of extrapolating the spline. Times
36 : * outside the bounds beyond roundoff trigger an `ASSERT` in debug builds.
37 : *
38 : * Here is an example how to use this class:
39 : *
40 : * \snippet Test_UniformCardinalBSpline.cpp uniform_cardinal_b_spline_example
41 : *
42 : * \note Requires Boost 1.81 or newer at runtime because of an accuracy fix
43 : * for the estimate of the derivative at the right endpoint in the Boost
44 : * implementation, see
45 : * https://github.com/boostorg/math/commit/4809e714d4806c07da3a3def0c4550daa0529b8d.
46 : * The constructor raises an error for older Boost versions.
47 : */
48 1 : class UniformCardinalBSpline {
49 : public:
50 : /*!
51 : * \brief Construct from uniformly spaced samples.
52 : *
53 : * \param values The sampled function values \f$f(t_0 + i \Delta t)\f$. At
54 : * least 5 samples are required.
55 : * \param start_time The time \f$t_0\f$ of the first sample.
56 : * \param time_step The (positive) spacing \f$\Delta t\f$ between samples.
57 : */
58 1 : UniformCardinalBSpline(std::vector<double> values, double start_time,
59 : double time_step);
60 :
61 0 : UniformCardinalBSpline() = default;
62 :
63 : /// Evaluate the interpolant at `time`, clamped to the `bounds()`
64 1 : double operator()(double time) const;
65 :
66 : /// The sampled function values
67 1 : const std::vector<double>& values() const { return values_; }
68 :
69 : /// The time of the first sample
70 1 : double start_time() const { return start_time_; }
71 :
72 : /// The spacing between samples
73 1 : double time_step() const { return time_step_; }
74 :
75 : /// The first and last sample time
76 1 : std::array<double, 2> bounds() const;
77 :
78 : // NOLINTNEXTLINE(google-runtime-references)
79 0 : void pup(PUP::er& p);
80 :
81 : private:
82 0 : void initialize_interpolant();
83 :
84 0 : std::vector<double> values_{};
85 0 : double start_time_{};
86 0 : double time_step_{};
87 : std::optional<boost::math::interpolators::cardinal_cubic_b_spline<double>>
88 0 : interpolant_{};
89 : };
90 :
91 0 : bool operator==(const UniformCardinalBSpline& lhs,
92 : const UniformCardinalBSpline& rhs);
93 :
94 0 : bool operator!=(const UniformCardinalBSpline& lhs,
95 : const UniformCardinalBSpline& rhs);
96 :
97 : /*!
98 : * \ingroup NumericalAlgorithmsGroup
99 : * \brief Estimate the maximum interpolation error of a
100 : * `intrp::UniformCardinalBSpline` through the given samples.
101 : *
102 : * Measures the deviation of an interpolant through every other sample from
103 : * the sampled values at the sample times, where the error of that coarser
104 : * interpolant peaks. For smooth data the interpolation error of a cubic
105 : * spline decreases by a factor of 16 when the step size is halved, so the
106 : * error of the interpolant through all samples is about 1/16 of the measured
107 : * deviation for smooth data in the convergent regime. The deviation is divided
108 : * by only 8 to obtain a conservative estimate, accounting for nonsmooth data
109 : * that converges slower.
110 : *
111 : * \param values The sampled function values. At least 9 samples are required
112 : * so the coarser interpolant has at least 5.
113 : * \param start_time The time of the first sample.
114 : * \param time_step The (positive) spacing between samples.
115 : */
116 1 : double estimate_interpolation_error(const std::vector<double>& values,
117 : double start_time, double time_step);
118 :
119 : /*!
120 : * \ingroup NumericalAlgorithmsGroup
121 : * \brief Compress uniformly spaced samples into a
122 : * `intrp::UniformCardinalBSpline` with fewer points that reproduces the
123 : * samples to the given tolerance.
124 : *
125 : * Builds a reference interpolant through all samples and resamples it on
126 : * coarser uniform grids over the same time interval, starting with 6 points
127 : * and doubling until the resampled interpolant reproduces the input samples
128 : * to within `absolute_tolerance` at the sample times. This reduces memory
129 : * for smooth, densely sampled data, e.g. time series of spectral modes in a
130 : * simulation.
131 : *
132 : * \returns The compressed interpolant and its maximum absolute deviation
133 : * from the input samples at the sample times.
134 : *
135 : * If the tolerance cannot be met with fewer points than the input, the
136 : * returned interpolant holds the original samples. Its deviation from the
137 : * samples at the sample times vanishes by construction, so the returned
138 : * error is the error of the last (coarser) candidate in the doubling search
139 : * as a conservative estimate, and may exceed the tolerance. Inputs with 6 or
140 : * fewer samples are returned unchanged with zero error.
141 : *
142 : * \param values The sampled function values. At least 5 samples are required.
143 : * \param start_time The time of the first sample.
144 : * \param time_step The (positive) spacing between samples.
145 : * \param absolute_tolerance Maximum allowed absolute deviation of the
146 : * compressed interpolant from the input samples.
147 : */
148 1 : std::pair<UniformCardinalBSpline, double> compress_to_tolerance(
149 : const std::vector<double>& values, double start_time, double time_step,
150 : double absolute_tolerance);
151 :
152 : } // namespace intrp
|
