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