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