boost/math/interpolators/detail/barycentric_rational_detail.hpp

   1 /*
   2  *  Copyright Nick Thompson, 2017
   3  *  Use, modification and distribution are subject to the
   4  *  Boost Software License, Version 1.0. (See accompanying file
   5  *  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7
   8 #ifndef BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
   9 #define BOOST_MATH_INTERPOLATORS_BARYCENTRIC_RATIONAL_DETAIL_HPP
  10
  11 #include <vector>
  12 #include <utility> // for std::move
  13 #include <algorithm> // for std::is_sorted
  14 #include <boost/lexical_cast.hpp>
  15 #include <boost/math/special_functions/fpclassify.hpp>
  16 #include <boost/core/demangle.hpp>
  17 #include <boost/assert.hpp>
  18
  19 namespace boost{ namespace math{ namespace detail{
  20
  21 template<class Real>
  22 class barycentric_rational_imp
  23 {
  24 public:
  25     template <class InputIterator1, class InputIterator2>
  26     barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order = 3);
  27
  28     barycentric_rational_imp(std::vector<Real>&& x, std::vector<Real>&& y, size_t approximation_order = 3);
  29
  30     Real operator()(Real x) const;
  31
  32     Real prime(Real x) const;
  33
  34     // The barycentric weights are not really that interesting; except to the unit tests!
  35     Real weight(size_t i) const { return m_w[i]; }
  36
  37     std::vector<Real>&& return_x()
  38     {
  39         return std::move(m_x);
  40     }
  41
  42     std::vector<Real>&& return_y()
  43     {
  44         return std::move(m_y);
  45     }
  46
  47 private:
  48
  49     void calculate_weights(size_t approximation_order);
  50
  51     std::vector<Real> m_x;
  52     std::vector<Real> m_y;
  53     std::vector<Real> m_w;
  54 };
  55
  56 template <class Real>
  57 template <class InputIterator1, class InputIterator2>
  58 barycentric_rational_imp<Real>::barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order)
  59 {
  60     std::ptrdiff_t n = std::distance(start_x, end_x);
  61
  62     if (approximation_order >= (std::size_t)n)
  63     {
  64         throw std::domain_error("Approximation order must be < data length.");
  65     }
  66
  67     // Big sad memcpy.
  68     m_x.resize(n);
  69     m_y.resize(n);
  70     for(unsigned i = 0; start_x != end_x; ++start_x, ++start_y, ++i)
  71     {
  72         // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy:
  73         if(boost::math::isnan(*start_x))
  74         {
  75             std::string msg = std::string("x[") + boost::lexical_cast<std::string>(i) + "] is a NAN";
  76             throw std::domain_error(msg);
  77         }
  78
  79         if(boost::math::isnan(*start_y))
  80         {
  81            std::string msg = std::string("y[") + boost::lexical_cast<std::string>(i) + "] is a NAN";
  82            throw std::domain_error(msg);
  83         }
  84
  85         m_x[i] = *start_x;
  86         m_y[i] = *start_y;
  87     }
  88     calculate_weights(approximation_order);
  89 }
  90
  91 template <class Real>
  92 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))
  93 {
  94     BOOST_ASSERT_MSG(m_x.size() == m_y.size(), "There must be the same number of abscissas and ordinates.");
  95     BOOST_ASSERT_MSG(approximation_order < m_x.size(), "Approximation order must be < data length.");
  96     BOOST_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].");
  97     calculate_weights(approximation_order);
  98 }
  99
 100 template<class Real>
 101 void barycentric_rational_imp<Real>::calculate_weights(size_t approximation_order)
 102 {
 103     using std::abs;
 104     int64_t n = m_x.size();
 105     m_w.resize(n, 0);
 106     for(int64_t k = 0; k < n; ++k)
 107     {
 108         int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0);
 109         int64_t i_max = k;
 110         if (k >= n - (std::ptrdiff_t)approximation_order)
 111         {
 112             i_max = n - approximation_order - 1;
 113         }
 114
 115         for(int64_t i = i_min; i <= i_max; ++i)
 116         {
 117             Real inv_product = 1;
 118             int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1));
 119             for(int64_t j = i; j <= j_max; ++j)
 120             {
 121                 if (j == k)
 122                 {
 123                     continue;
 124                 }
 125
 126                 Real diff = m_x[k] - m_x[j];
 127                 using std::numeric_limits;
 128                 if (abs(diff) < (numeric_limits<Real>::min)())
 129                 {
 130                    std::string msg = std::string("Spacing between  x[")
 131                       + boost::lexical_cast<std::string>(k) + std::string("] and x[")
 132                       + boost::lexical_cast<std::string>(i) + std::string("] is ")
 133                       + boost::lexical_cast<std::string>(diff) + std::string(", which is smaller than the epsilon of ")
 134                       + boost::core::demangle(typeid(Real).name());
 135                     throw std::logic_error(msg);
 136                 }
 137                 inv_product *= diff;
 138             }
 139             if (i % 2 == 0)
 140             {
 141                 m_w[k] += 1/inv_product;
 142             }
 143             else
 144             {
 145                 m_w[k] -= 1/inv_product;
 146             }
 147         }
 148     }
 149 }
 150
 151
 152 template<class Real>
 153 Real barycentric_rational_imp<Real>::operator()(Real x) const
 154 {
 155     Real numerator = 0;
 156     Real denominator = 0;
 157     for(size_t i = 0; i < m_x.size(); ++i)
 158     {
 159         // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality.
 160         // 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,
 161         // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715
 162         if (x == m_x[i])
 163         {
 164             return m_y[i];
 165         }
 166         Real t = m_w[i]/(x - m_x[i]);
 167         numerator += t*m_y[i];
 168         denominator += t;
 169     }
 170     return numerator/denominator;
 171 }
 172
 173 /*
 174  * A formula for computing the derivative of the barycentric representation is given in
 175  * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner,
 176  * Mathematics of Computation, v47, number 175, 1986.
 177  * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf
 178  * and reviewed in
 179  * Recent developments in barycentric rational interpolation
 180  * Jean-Paul Berrut, Richard Baltensperger and Hans D. Mittelmann
 181  *
 182  * Is it possible to complete this in one pass through the data?
 183  */
 184
 185 template<class Real>
 186 Real barycentric_rational_imp<Real>::prime(Real x) const
 187 {
 188     Real rx = this->operator()(x);
 189     Real numerator = 0;
 190     Real denominator = 0;
 191     for(size_t i = 0; i < m_x.size(); ++i)
 192     {
 193         if (x == m_x[i])
 194         {
 195             Real sum = 0;
 196             for (size_t j = 0; j < m_x.size(); ++j)
 197             {
 198                 if (j == i)
 199                 {
 200                     continue;
 201                 }
 202                 sum += m_w[j]*(m_y[i] - m_y[j])/(m_x[i] - m_x[j]);
 203             }
 204             return -sum/m_w[i];
 205         }
 206         Real t = m_w[i]/(x - m_x[i]);
 207         Real diff = (rx - m_y[i])/(x-m_x[i]);
 208         numerator += t*diff;
 209         denominator += t;
 210     }
 211
 212     return numerator/denominator;
 213 }
 214 }}}
 215 #endif