Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <unistd.h>  // IWYU pragma: keep
       8             : #include <vector>
       9             : 
      10             : namespace PUP {
      11             : class er;
      12             : }  // namespace PUP
      13             : 
      14             : namespace intrp {
      15             : /*!
      16             :  * \ingroup NumericalAlgorithmsGroup
      17             :  * \brief A barycentric rational interpolation class
      18             :  *
      19             :  * The class builds a barycentric rational interpolant of a specified order
      20             :  * using the `x_values` and `y_values` passed into the constructor.
      21             :  * Barycentric interpolation requires \f$3N\f$ storage, and costs
      22             :  * \f$\mathcal{O}(Nd)\f$ to construct, where \f$N\f$ is the size of the x- and
      23             :  * y-value vectors and \f$d\f$ is the order of the interpolant. The evaluation
      24             :  * cost is \f$\mathcal{O}(N)\f$ compared to \f$\mathcal{O}(d)\f$ of a spline
      25             :  * method, but constructing the barycentric interpolant does not require any
      26             :  * derivatives of the function to be known.
      27             :  *
      28             :  * The interpolation function is
      29             :  *
      30             :  * \f[
      31             :  *   \mathcal{I}(x)=\frac{\sum_{i=0}^{N-1}w_i y_i /
      32             :  *   (x-x_i)}{\sum_{i=0}^{N-1}w_i/(x-x_i)}
      33             :  * \f]
      34             :  *
      35             :  * where \f$w_i\f$ are the weights. The weights are computed using
      36             :  *
      37             :  * \f[
      38             :  *   w_k=\sum_{i=k-d\\0\le i < N-d}^{k}(-1)^{i}
      39             :  *       \prod_{j=i\\j\ne k}^{i+d}\frac{1}{x_k-x_j} \f]
      40             :  *
      41             :  * \requires `x_values.size() == y_values.size()` and
      42             :  * `x_values_.size() >= order`
      43             :  */
      44           1 : class BarycentricRational {
      45             :  public:
      46           0 :   BarycentricRational() = default;
      47           0 :   BarycentricRational(std::vector<double> x_values,
      48             :                       std::vector<double> y_values, size_t order);
      49             : 
      50           0 :   double operator()(double x_to_interp_to) const;
      51             : 
      52           0 :   const std::vector<double>& x_values() const { return x_values_; }
      53             : 
      54           0 :   const std::vector<double>& y_values() const { return y_values_; }
      55             : 
      56           0 :   size_t order() const;
      57             : 
      58             :   // NOLINTNEXTLINE(google-runtime-references)
      59           0 :   void pup(PUP::er& p);
      60             : 
      61             :  private:
      62           0 :   std::vector<double> x_values_;
      63           0 :   std::vector<double> y_values_;
      64           0 :   std::vector<double> weights_;
      65           0 :   ssize_t order_{0};
      66             : };
      67             : }  // namespace intrp