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