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