boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp

   1 /*
   2  *  Copyright Nick Thompson, 2019
   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_VECTOR_BARYCENTRIC_RATIONAL_DETAIL_HPP
   9 #define BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_DETAIL_HPP
  10
  11 #include <vector>
  12 #include <utility> // for std::move
  13 #include <boost/assert.hpp>
  14
  15 namespace boost{ namespace math{ namespace detail{
  16
  17 template <class TimeContainer, class SpaceContainer>
  18 class vector_barycentric_rational_imp
  19 {
  20 public:
  21     using Real = typename TimeContainer::value_type;
  22     using Point = typename SpaceContainer::value_type;
  23
  24     vector_barycentric_rational_imp(TimeContainer&& t, SpaceContainer&& y, size_t approximation_order);
  25
  26     void operator()(Point& p, Real t) const;
  27
  28     void eval_with_prime(Point& x, Point& dxdt, Real t) const;
  29
  30     // The barycentric weights are only interesting to the unit tests:
  31     Real weight(size_t i) const { return w_[i]; }
  32
  33 private:
  34
  35     void calculate_weights(size_t approximation_order);
  36
  37     TimeContainer t_;
  38     SpaceContainer y_;
  39     TimeContainer w_;
  40 };
  41
  42 template <class TimeContainer, class SpaceContainer>
  43 vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::vector_barycentric_rational_imp(TimeContainer&& t, SpaceContainer&& y, size_t approximation_order)
  44 {
  45     using std::numeric_limits;
  46     t_ = std::move(t);
  47     y_ = std::move(y);
  48
  49     BOOST_ASSERT_MSG(t_.size() == y_.size(), "There must be the same number of time points as space points.");
  50     BOOST_ASSERT_MSG(approximation_order < y_.size(), "Approximation order must be < data length.");
  51     for (size_t i = 1; i < t_.size(); ++i)
  52     {
  53         BOOST_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].");
  54     }
  55     calculate_weights(approximation_order);
  56 }
  57
  58
  59 template<class TimeContainer, class SpaceContainer>
  60 void vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::calculate_weights(size_t approximation_order)
  61 {
  62     using Real = typename TimeContainer::value_type;
  63     using std::abs;
  64     int64_t n = t_.size();
  65     w_.resize(n, Real(0));
  66     for(int64_t k = 0; k < n; ++k)
  67     {
  68         int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0);
  69         int64_t i_max = k;
  70         if (k >= n - (std::ptrdiff_t)approximation_order)
  71         {
  72             i_max = n - approximation_order - 1;
  73         }
  74
  75         for(int64_t i = i_min; i <= i_max; ++i)
  76         {
  77             Real inv_product = 1;
  78             int64_t j_max = (std::min)(static_cast<int64_t>(i + approximation_order), static_cast<int64_t>(n - 1));
  79             for(int64_t j = i; j <= j_max; ++j)
  80             {
  81                 if (j == k)
  82                 {
  83                     continue;
  84                 }
  85                 Real diff = t_[k] - t_[j];
  86                 inv_product *= diff;
  87             }
  88             if (i % 2 == 0)
  89             {
  90                 w_[k] += 1/inv_product;
  91             }
  92             else
  93             {
  94                 w_[k] -= 1/inv_product;
  95             }
  96         }
  97     }
  98 }
  99
 100
 101 template<class TimeContainer, class SpaceContainer>
 102 void vector_barycentric_rational_imp<TimeContainer, SpaceContainer>::operator()(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const
 103 {
 104     using Real = typename TimeContainer::value_type;
 105     for (auto & x : p)
 106     {
 107         x = Real(0);
 108     }
 109     Real denominator = 0;
 110     for(size_t i = 0; i < t_.size(); ++i)
 111     {
 112         // See associated commentary in the scalar version of this function.
 113         if (t == t_[i])
 114         {
 115             p = y_[i];
 116             return;
 117         }
 118         Real x = w_[i]/(t - t_[i]);
 119         for (decltype(p.size()) j = 0; j < p.size(); ++j)
 120         {
 121             p[j] += x*y_[i][j];
 122         }
 123         denominator += x;
 124     }
 125     for (decltype(p.size()) j = 0; j < p.size(); ++j)
 126     {
 127         p[j] /= denominator;
 128     }
 129     return;
 130 }
 131
 132 template<class TimeContainer, class SpaceContainer>
 133 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
 134 {
 135     using Point = typename SpaceContainer::value_type;
 136     using Real = typename TimeContainer::value_type;
 137     this->operator()(x, t);
 138     Point numerator;
 139     for (decltype(x.size()) i = 0; i < x.size(); ++i)
 140     {
 141         numerator[i] = 0;
 142     }
 143     Real denominator = 0;
 144     for(decltype(t_.size()) i = 0; i < t_.size(); ++i)
 145     {
 146         if (t == t_[i])
 147         {
 148             Point sum;
 149             for (decltype(x.size()) i = 0; i < x.size(); ++i)
 150             {
 151                 sum[i] = 0;
 152             }
 153
 154             for (decltype(t_.size()) j = 0; j < t_.size(); ++j)
 155             {
 156                 if (j == i)
 157                 {
 158                     continue;
 159                 }
 160                 for (decltype(sum.size()) k = 0; k < sum.size(); ++k)
 161                 {
 162                     sum[k] += w_[j]*(y_[i][k] - y_[j][k])/(t_[i] - t_[j]);
 163                 }
 164             }
 165             for (decltype(sum.size()) k = 0; k < sum.size(); ++k)
 166             {
 167                 dxdt[k] = -sum[k]/w_[i];
 168             }
 169             return;
 170         }
 171         Real tw = w_[i]/(t - t_[i]);
 172         Point diff;
 173         for (decltype(diff.size()) j = 0; j < diff.size(); ++j)
 174         {
 175             diff[j] = (x[j] - y_[i][j])/(t-t_[i]);
 176         }
 177         for (decltype(diff.size()) j = 0; j < diff.size(); ++j)
 178         {
 179             numerator[j] += tw*diff[j];
 180         }
 181         denominator += tw;
 182     }
 183
 184     for (decltype(dxdt.size()) j = 0; j < dxdt.size(); ++j)
 185     {
 186         dxdt[j] = numerator[j]/denominator;
 187     }
 188     return;
 189 }
 190
 191 }}}
 192 #endif