interpolators/detail/barycentric_rational_detail.hpp

0001 /*
0002  *  Copyright Nick Thompson, 2017
0003  *  Use, modification and distribution are subject to the
0004  *  Boost Software License, Version 1.0. (See accompanying file
0005  *  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006  */
0007
0008 #ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
0009 #define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
0010
0011 #include <vector>
0012 #include <utility> // for std::move
0013 #include <algorithm> // for std::is_sorted
0014 #include <string>
0015 #include <boost/math/special_functions/fpclassify.hpp>
0016 #include <boost/math/tools/assert.hpp>
0017
0018 namespace boost{ namespace math{ namespace interpolators { namespace detail{
0019
0020 template<class Real>
0021 class barycentric_rational_imp
0022 {
0023 public:
0024     template <class InputIterator1, class InputIterator2>
0025     barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3);
0026
0027     barycentric_rational_imp(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order = 3);
0028
0029     Real operator()(Real x) const;
0030
0031     Real prime(Real x) const;
0032
0033     // The barycentric weights are not really that interesting; except to the unit tests!
0034     Real weight(size_t i) const { return m_w[i]; }
0035
0036     std::vector<Real>&& return_x()
0037     {
0038         return std::move(m_x);
0039     }
0040
0041     std::vector<Real>&& return_y()
0042     {
0043         return std::move(m_y);
0044     }
0045
0046 private:
0047
0048     void calculate_weights(size_t approximation_order);
0049
0050     std::vector<Real> m_x;
0051     std::vector<Real> m_y;
0052     std::vector<Real> m_w;
0053 };
0054
0055 template <class Real>
0056 template <class InputIterator1, class InputIterator2>
0057 barycentric_rational_imp<Real>::barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order)
0058 {
0059     std::ptrdiff_t n = std::distance(start_x, end_x);
0060
0061     if (approximation_order >= (std::size_t)n)
0062     {
0063         throw std::domain_error("Approximation order must be < data length.");
0064     }
0065
0066     // Big sad memcpy.
0067     m_x.resize(n);
0068     m_y.resize(n);
0069     for(unsigned i = 0; start_x != end_x; ++start_x, ++start_y, ++i)
0070     {
0071         // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy:
0072         if(boost::math::isnan(*start_x))
0073         {
0074             std::string msg = std::string("x[") + std::to_string(i) + "] is a NAN";
0075             throw std::domain_error(msg);
0076         }
0077
0078         if(boost::math::isnan(*start_y))
0079         {
0080            std::string msg = std::string("y[") + std::to_string(i) + "] is a NAN";
0081            throw std::domain_error(msg);
0082         }
0083
0084         m_x[i] = *start_x;
0085         m_y[i] = *start_y;
0086     }
0087     calculate_weights(approximation_order);
0088 }
0089
0090 template <class Real>
0091 barycentric_rational_imp<Real>::barycentric_rational_imp(std::vector<Real>&& x, std::vector<Real>&& y,size_t approximation_order) : m_x(std::move(x)), m_y(std::move(y))
0092 {
0093     BOOST_MATH_ASSERT_MSG(m_x.size() == m_y.size(), "There must be the same number of abscissas and ordinates.");
0094     BOOST_MATH_ASSERT_MSG(approximation_order < m_x.size(), "Approximation order must be < data length.");
0095     BOOST_MATH_ASSERT_MSG(std::is_sorted(m_x.begin(), m_x.end()), "The abscissas must be listed in increasing order x[0] < x[1] < ... < x[n-1].");
0096     calculate_weights(approximation_order);
0097 }
0098
0099 template<class Real>
0100 void barycentric_rational_imp<Real>::calculate_weights(size_t approximation_order)
0101 {
0102     using std::abs;
0103     int64_t n = m_x.size();
0104     m_w.resize(n, 0);
0105     for(int64_t k = 0; k < n; ++k)
0106     {
0107         int64_t i_min = (std::max)(k - static_cast<int64_t>(approximation_order), static_cast<int64_t>(0));
0108         int64_t i_max = k;
0109         if (k >= n - (std::ptrdiff_t)approximation_order)
0110         {
0111             i_max = n - approximation_order - 1;
0112         }
0113
0114         for(int64_t i = i_min; i <= i_max; ++i)
0115         {
0116             Real inv_product = 1;
0117             int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1));
0118             for(int64_t j = i; j <= j_max; ++j)
0119             {
0120                 if (j == k)
0121                 {
0122                     continue;
0123                 }
0124
0125                 Real diff = m_x[k] - m_x[j];
0126                 using std::numeric_limits;
0127                 if (abs(diff) < (numeric_limits<Real>::min)())
0128                 {
0129                    std::string msg = std::string("Spacing between  x[")
0130                       + std::to_string(k) + std::string("] and x[")
0131                       + std::to_string(i) + std::string("] is ")
0132                       + std::string("smaller than the epsilon of ")
0133                       + std::string(typeid(Real).name());
0134                     throw std::logic_error(msg);
0135                 }
0136                 inv_product *= diff;
0137             }
0138             if (i % 2 == 0)
0139             {
0140                 m_w[k] += 1/inv_product;
0141             }
0142             else
0143             {
0144                 m_w[k] -= 1/inv_product;
0145             }
0146         }
0147     }
0148 }
0149
0150
0151 template<class Real>
0152 Real barycentric_rational_imp<Real>::operator()(Real x) const
0153 {
0154     Real numerator = 0;
0155     Real denominator = 0;
0156     for(size_t i = 0; i < m_x.size(); ++i)
0157     {
0158         // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality.
0159         // However, it has been shown that if x approx x_i, but x != x_i, then inaccuracy in the numerator cancels the inaccuracy in the denominator,
0160         // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715
0161         if (x == m_x[i])
0162         {
0163             return m_y[i];
0164         }
0165         Real t = m_w[i]/(x - m_x[i]);
0166         numerator += t*m_y[i];
0167         denominator += t;
0168     }
0169     return numerator/denominator;
0170 }
0171
0172 /*
0173  * A formula for computing the derivative of the barycentric representation is given in
0174  * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner,
0175  * Mathematics of Computation, v47, number 175, 1986.
0176  * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf
0177  * and reviewed in
0178  * Recent developments in barycentric rational interpolation
0179  * Jean-Paul Berrut, Richard Baltensperger and Hans D. Mittelmann
0180  *
0181  * Is it possible to complete this in one pass through the data?
0182  */
0183
0184 template<class Real>
0185 Real barycentric_rational_imp<Real>::prime(Real x) const
0186 {
0187     Real rx = this->operator()(x);
0188     Real numerator = 0;
0189     Real denominator = 0;
0190     for(size_t i = 0; i < m_x.size(); ++i)
0191     {
0192         if (x == m_x[i])
0193         {
0194             Real sum = 0;
0195             for (size_t j = 0; j < m_x.size(); ++j)
0196             {
0197                 if (j == i)
0198                 {
0199                     continue;
0200                 }
0201                 sum += m_w[j]*(m_y[i] - m_y[j])/(m_x[i] - m_x[j]);
0202             }
0203             return -sum/m_w[i];
0204         }
0205         Real t = m_w[i]/(x - m_x[i]);
0206         Real diff = (rx - m_y[i])/(x-m_x[i]);
0207         numerator += t*diff;
0208         denominator += t;
0209     }
0210
0211     return numerator/denominator;
0212 }
0213 }}}}
0214 #endif