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>
       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