interpolators/detail/vector_barycentric_rational_detail.hpp

0001 /*
0002  *  Copyright Nick Thompson, 2019
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_VECTOR_BARYCENTRIC_RATIONAL_DETAIL_HPP
0009 #define BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_DETAIL_HPP
0010
0011 #include <cmath>
0012 #include <vector>
0013 #include <utility> // for std::move
0014 #include <limits>
0015 #include <algorithm>
0016 #include <boost/math/tools/assert.hpp>
0017
0018 namespace boost{ namespace math{ namespace interpolators{ namespace detail{
0019
0020 template <class TimeContainer, class SpaceContainer>
0021 class vector_barycentric_rational_imp
0022 {
0023 public:
0024     using Real = typename TimeContainer::value_type;
0025     using Point = typename SpaceContainer::value_type;
0026
0027     vector_barycentric_rational_imp(TimeContainer&& t, SpaceContainer&& y, size_t approximation_order);
0028
0029     void operator()(Point& p, Real t) const;
0030
0031     void eval_with_prime(Point& x, Point& dxdt, Real t) const;
0032
0033     // The barycentric weights are only interesting to the unit tests:
0034     Real weight(size_t i) const { return w_[i]; }
0035
0036 private:
0037
0038     void calculate_weights(size_t approximation_order);
0039
0040     TimeContainer t_;
0041     SpaceContainer y_;
0042     TimeContainer w_;
0043 };
0044
0045 template <class TimeContainer, class SpaceContainer>
0046 vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::vector_barycentric_rational_imp(TimeContainer&& t, SpaceContainer&& y, size_t approximation_order)
0047 {
0048     using std::numeric_limits;
0049     t_ = std::move(t);
0050     y_ = std::move(y);
0051
0052     BOOST_MATH_ASSERT_MSG(t_.size() == y_.size(), "There must be the same number of time points as space points.");
0053     BOOST_MATH_ASSERT_MSG(approximation_order < y_.size(), "Approximation order must be < data length.");
0054     for (size_t i = 1; i < t_.size(); ++i)
0055     {
0056         BOOST_MATH_ASSERT_MSG(t_[i] - t_[i-1] >  (numeric_limits<typename TimeContainer::value_type>::min)(), "The abscissas must be listed in strictly increasing order t[0] < t[1] < ... < t[n-1].");
0057     }
0058     calculate_weights(approximation_order);
0059 }
0060
0061
0062 template<class TimeContainer, class SpaceContainer>
0063 void vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::calculate_weights(size_t approximation_order)
0064 {
0065     using Real = typename TimeContainer::value_type;
0066     using std::abs;
0067     int64_t n = t_.size();
0068     w_.resize(n, Real(0));
0069     for(int64_t k = 0; k < n; ++k)
0070     {
0071         int64_t i_min = (std::max)(k - static_cast<int64_t>(approximation_order), static_cast<int64_t>(0));
0072         int64_t i_max = k;
0073         if (k >= n - (std::ptrdiff_t)approximation_order)
0074         {
0075             i_max = n - approximation_order - 1;
0076         }
0077
0078         for(int64_t i = i_min; i <= i_max; ++i)
0079         {
0080             Real inv_product = 1;
0081             int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1));
0082             for(int64_t j = i; j <= j_max; ++j)
0083             {
0084                 if (j == k)
0085                 {
0086                     continue;
0087                 }
0088                 Real diff = t_[k] - t_[j];
0089                 inv_product *= diff;
0090             }
0091             if (i % 2 == 0)
0092             {
0093                 w_[k] += 1/inv_product;
0094             }
0095             else
0096             {
0097                 w_[k] -= 1/inv_product;
0098             }
0099         }
0100     }
0101 }
0102
0103
0104 template<class TimeContainer, class SpaceContainer>
0105 void vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::operator()(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const
0106 {
0107     using Real = typename TimeContainer::value_type;
0108     for (auto & x : p)
0109     {
0110         x = Real(0);
0111     }
0112     Real denominator = 0;
0113     for(size_t i = 0; i < t_.size(); ++i)
0114     {
0115         // See associated commentary in the scalar version of this function.
0116         if (t == t_[i])
0117         {
0118             p = y_[i];
0119             return;
0120         }
0121         Real x = w_[i]/(t - t_[i]);
0122         for (decltype(p.size()) j = 0; j < p.size(); ++j)
0123         {
0124             p[j] += x*y_[i][j];
0125         }
0126         denominator += x;
0127     }
0128     for (decltype(p.size()) j = 0; j < p.size(); ++j)
0129     {
0130         p[j] /= denominator;
0131     }
0132     return;
0133 }
0134
0135 template<class TimeContainer, class SpaceContainer>
0136 void vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::eval_with_prime(typename SpaceContainer::value_type& x, typename SpaceContainer::value_type& dxdt, typename TimeContainer::value_type t) const
0137 {
0138     using Point = typename SpaceContainer::value_type;
0139     using Real = typename TimeContainer::value_type;
0140     this->operator()(x, t);
0141     Point numerator;
0142     for (decltype(x.size()) i = 0; i < x.size(); ++i)
0143     {
0144         numerator[i] = 0;
0145     }
0146     Real denominator = 0;
0147     for(decltype(t_.size()) i = 0; i < t_.size(); ++i)
0148     {
0149         if (t == t_[i])
0150         {
0151             Point sum;
0152             for (decltype(x.size()) i = 0; i < x.size(); ++i)
0153             {
0154                 sum[i] = 0;
0155             }
0156
0157             for (decltype(t_.size()) j = 0; j < t_.size(); ++j)
0158             {
0159                 if (j == i)
0160                 {
0161                     continue;
0162                 }
0163                 for (decltype(sum.size()) k = 0; k < sum.size(); ++k)
0164                 {
0165                     sum[k] += w_[j]*(y_[i][k] - y_[j][k])/(t_[i] - t_[j]);
0166                 }
0167             }
0168             for (decltype(sum.size()) k = 0; k < sum.size(); ++k)
0169             {
0170                 dxdt[k] = -sum[k]/w_[i];
0171             }
0172             return;
0173         }
0174         Real tw = w_[i]/(t - t_[i]);
0175         Point diff;
0176         for (decltype(diff.size()) j = 0; j < diff.size(); ++j)
0177         {
0178             diff[j] = (x[j] - y_[i][j])/(t-t_[i]);
0179         }
0180         for (decltype(diff.size()) j = 0; j < diff.size(); ++j)
0181         {
0182             numerator[j] += tw*diff[j];
0183         }
0184         denominator += tw;
0185     }
0186
0187     for (decltype(dxdt.size()) j = 0; j < dxdt.size(); ++j)
0188     {
0189         dxdt[j] = numerator[j]/denominator;
0190     }
0191     return;
0192 }
0193
0194 }}}}
0195 #endif