interpolators/detail/cubic_hermite_detail.hpp

0001 // Copyright Nick Thompson, 2020
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_INTERPOLATORS_DETAIL_CUBIC_HERMITE_DETAIL_HPP
0008 #define BOOST_MATH_INTERPOLATORS_DETAIL_CUBIC_HERMITE_DETAIL_HPP
0009 #include <stdexcept>
0010 #include <algorithm>
0011 #include <cmath>
0012 #include <iostream>
0013 #include <sstream>
0014 #include <limits>
0015
0016 namespace boost {
0017 namespace math {
0018 namespace interpolators {
0019 namespace detail {
0020
0021 template<class RandomAccessContainer>
0022 class cubic_hermite_detail {
0023 public:
0024     using Real = typename RandomAccessContainer::value_type;
0025     using Size = typename RandomAccessContainer::size_type;
0026
0027     cubic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer dydx)
0028      : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}
0029     {
0030         using std::abs;
0031         using std::isnan;
0032         if (x_.size() != y_.size())
0033         {
0034             throw std::domain_error("There must be the same number of ordinates as abscissas.");
0035         }
0036         if (x_.size() != dydx_.size())
0037         {
0038             throw std::domain_error("There must be the same number of ordinates as derivative values.");
0039         }
0040         if (x_.size() < 2)
0041         {
0042             throw std::domain_error("Must be at least two data points.");
0043         }
0044         Real x0 = x_[0];
0045         for (size_t i = 1; i < x_.size(); ++i)
0046         {
0047             Real x1 = x_[i];
0048             if (x1 <= x0)
0049             {
0050                 std::ostringstream oss;
0051                 oss.precision(std::numeric_limits<Real>::digits10+3);
0052                 oss << "Abscissas must be listed in strictly increasing order x0 < x1 < ... < x_{n-1}, ";
0053                 oss << "but at x[" << i - 1 << "] = " << x0 << ", and x[" << i << "] = " << x1 << ".\n";
0054                 throw std::domain_error(oss.str());
0055             }
0056             x0 = x1;
0057         }
0058     }
0059
0060     void push_back(Real x, Real y, Real dydx)
0061     {
0062         using std::abs;
0063         using std::isnan;
0064         if (x <= x_.back())
0065         {
0066              throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
0067         }
0068         x_.push_back(x);
0069         y_.push_back(y);
0070         dydx_.push_back(dydx);
0071     }
0072
0073     Real operator()(Real x) const
0074     {
0075         if  (x < x_[0] || x > x_.back())
0076         {
0077             std::ostringstream oss;
0078             oss.precision(std::numeric_limits<Real>::digits10+3);
0079             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0080                 << x_[0] << ", " << x_.back() << "]";
0081             throw std::domain_error(oss.str());
0082         }
0083         // We need t := (x-x_k)/(x_{k+1}-x_k) \in [0,1) for this to work.
0084         // Sadly this neccessitates this loathesome check, otherwise we get t = 1 at x = xf.
0085         if (x == x_.back())
0086         {
0087             return y_.back();
0088         }
0089
0090         auto it = std::upper_bound(x_.begin(), x_.end(), x);
0091         auto i = std::distance(x_.begin(), it) -1;
0092         Real x0 = *(it-1);
0093         Real x1 = *it;
0094         Real y0 = y_[i];
0095         Real y1 = y_[i+1];
0096         Real s0 = dydx_[i];
0097         Real s1 = dydx_[i+1];
0098         Real dx = (x1-x0);
0099         Real t = (x-x0)/dx;
0100
0101         // See the section 'Representations' in the page
0102         // https://en.wikipedia.org/wiki/Cubic_Hermite_spline
0103         Real y = (1-t)*(1-t)*(y0*(1+2*t) + s0*(x-x0))
0104               + t*t*(y1*(3-2*t) + dx*s1*(t-1));
0105         return y;
0106     }
0107
0108     Real prime(Real x) const
0109     {
0110         if  (x < x_[0] || x > x_.back())
0111         {
0112             std::ostringstream oss;
0113             oss.precision(std::numeric_limits<Real>::digits10+3);
0114             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0115                 << x_[0] << ", " << x_.back() << "]";
0116             throw std::domain_error(oss.str());
0117         }
0118         if (x == x_.back())
0119         {
0120             return dydx_.back();
0121         }
0122         auto it = std::upper_bound(x_.begin(), x_.end(), x);
0123         auto i = std::distance(x_.begin(), it) -1;
0124         Real x0 = *(it-1);
0125         Real x1 = *it;
0126         Real y0 = y_[i];
0127         Real y1 = y_[i+1];
0128         Real s0 = dydx_[i];
0129         Real s1 = dydx_[i+1];
0130         Real dx = (x1-x0);
0131
0132         Real d1 = (y1 - y0 - s0*dx)/(dx*dx);
0133         Real d2 = (s1 - s0)/(2*dx);
0134         Real c2 = 3*d1 - 2*d2;
0135         Real c3 = 2*(d2 - d1)/dx;
0136         return s0 + 2*c2*(x-x0) + 3*c3*(x-x0)*(x-x0);
0137     }
0138
0139
0140     friend std::ostream& operator<<(std::ostream & os, const cubic_hermite_detail & m)
0141     {
0142         os << "(x,y,y') = {";
0143         for (size_t i = 0; i < m.x_.size() - 1; ++i)
0144         {
0145             os << "(" << m.x_[i] << ", " << m.y_[i] << ", " << m.dydx_[i] << "),  ";
0146         }
0147         auto n = m.x_.size()-1;
0148         os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ")}";
0149         return os;
0150     }
0151
0152     Size size() const
0153     {
0154         return x_.size();
0155     }
0156
0157     int64_t bytes() const
0158     {
0159         return 3*x_.size()*sizeof(Real) + 3*sizeof(x_);
0160     }
0161
0162     std::pair<Real, Real> domain() const
0163     {
0164         return {x_.front(), x_.back()};
0165     }
0166
0167     RandomAccessContainer x_;
0168     RandomAccessContainer y_;
0169     RandomAccessContainer dydx_;
0170 };
0171
0172 template<class RandomAccessContainer>
0173 class cardinal_cubic_hermite_detail {
0174 public:
0175     using Real = typename RandomAccessContainer::value_type;
0176     using Size = typename RandomAccessContainer::size_type;
0177
0178     cardinal_cubic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer dydx, Real x0, Real dx)
0179     : y_{std::move(y)}, dy_{std::move(dydx)}, x0_{x0}, inv_dx_{1/dx}
0180     {
0181         using std::abs;
0182         using std::isnan;
0183         if (y_.size() != dy_.size())
0184         {
0185             throw std::domain_error("There must be the same number of derivatives as ordinates.");
0186         }
0187         if (y_.size() < 2)
0188         {
0189             throw std::domain_error("Must be at least two data points.");
0190         }
0191         if (dx <= 0)
0192         {
0193             throw std::domain_error("dx > 0 is required.");
0194         }
0195
0196         for (auto & dy : dy_)
0197         {
0198             dy *= dx;
0199         }
0200     }
0201
0202     // Why not implement push_back? It's awkward: If the buffer is circular, x0_ += dx_.
0203     // If the buffer is not circular, x0_ is unchanged.
0204     // We need a concept for circular_buffer!
0205
0206     inline Real operator()(Real x) const
0207     {
0208         const Real xf = x0_ + (y_.size()-1)/inv_dx_;
0209         if  (x < x0_ || x > xf)
0210         {
0211             std::ostringstream oss;
0212             oss.precision(std::numeric_limits<Real>::digits10+3);
0213             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0214                 << x0_ << ", " << xf << "]";
0215             throw std::domain_error(oss.str());
0216         }
0217         if (x == xf)
0218         {
0219             return y_.back();
0220         }
0221         return this->unchecked_evaluation(x);
0222     }
0223
0224     inline Real unchecked_evaluation(Real x) const
0225     {
0226         using std::floor;
0227         Real s = (x-x0_)*inv_dx_;
0228         Real ii = floor(s);
0229         auto i = static_cast<decltype(y_.size())>(ii);
0230         Real t = s - ii;
0231         Real y0 = y_[i];
0232         Real y1 = y_[i+1];
0233         Real dy0 = dy_[i];
0234         Real dy1 = dy_[i+1];
0235
0236         Real r = 1-t;
0237         return r*r*(y0*(1+2*t) + dy0*t)
0238               + t*t*(y1*(3-2*t) - dy1*r);
0239     }
0240
0241     inline Real prime(Real x) const
0242     {
0243         const Real xf = x0_ + (y_.size()-1)/inv_dx_;
0244         if  (x < x0_ || x > xf)
0245         {
0246             std::ostringstream oss;
0247             oss.precision(std::numeric_limits<Real>::digits10+3);
0248             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0249                 << x0_ << ", " << xf << "]";
0250             throw std::domain_error(oss.str());
0251         }
0252         if (x == xf)
0253         {
0254             return dy_.back()*inv_dx_;
0255         }
0256         return this->unchecked_prime(x);
0257     }
0258
0259     inline Real unchecked_prime(Real x) const
0260     {
0261         using std::floor;
0262         Real s = (x-x0_)*inv_dx_;
0263         Real ii = floor(s);
0264         auto i = static_cast<decltype(y_.size())>(ii);
0265         Real t = s - ii;
0266         Real y0 = y_[i];
0267         Real y1 = y_[i+1];
0268         Real dy0 = dy_[i];
0269         Real dy1 = dy_[i+1];
0270
0271         Real dy = 6*t*(1-t)*(y1 - y0)  + (3*t*t - 4*t+1)*dy0 + t*(3*t-2)*dy1;
0272         return dy*inv_dx_;
0273     }
0274
0275
0276     Size size() const
0277     {
0278         return y_.size();
0279     }
0280
0281     int64_t bytes() const
0282     {
0283         return 2*y_.size()*sizeof(Real) + 2*sizeof(y_) + 2*sizeof(Real);
0284     }
0285
0286     std::pair<Real, Real> domain() const
0287     {
0288         Real xf = x0_ + (y_.size()-1)/inv_dx_;
0289         return {x0_, xf};
0290     }
0291
0292 private:
0293
0294     RandomAccessContainer y_;
0295     RandomAccessContainer dy_;
0296     Real x0_;
0297     Real inv_dx_;
0298 };
0299
0300
0301 template<class RandomAccessContainer>
0302 class cardinal_cubic_hermite_detail_aos {
0303 public:
0304     using Point = typename RandomAccessContainer::value_type;
0305     using Real = typename Point::value_type;
0306     using Size = typename RandomAccessContainer::size_type;
0307
0308     cardinal_cubic_hermite_detail_aos(RandomAccessContainer && dat, Real x0, Real dx)
0309     : dat_{std::move(dat)}, x0_{x0}, inv_dx_{1/dx}
0310     {
0311         if (dat_.size() < 2)
0312         {
0313             throw std::domain_error("Must be at least two data points.");
0314         }
0315         if (dat_[0].size() != 2)
0316         {
0317             throw std::domain_error("Each datum must contain (y, y'), and nothing else.");
0318         }
0319         if (dx <= 0)
0320         {
0321             throw std::domain_error("dx > 0 is required.");
0322         }
0323
0324         for (auto & d : dat_)
0325         {
0326             d[1] *= dx;
0327         }
0328     }
0329
0330     inline Real operator()(Real x) const
0331     {
0332         const Real xf = x0_ + (dat_.size()-1)/inv_dx_;
0333         if  (x < x0_ || x > xf)
0334         {
0335             std::ostringstream oss;
0336             oss.precision(std::numeric_limits<Real>::digits10+3);
0337             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0338                 << x0_ << ", " << xf << "]";
0339             throw std::domain_error(oss.str());
0340         }
0341         if (x == xf)
0342         {
0343             return dat_.back()[0];
0344         }
0345         return this->unchecked_evaluation(x);
0346     }
0347
0348     inline Real unchecked_evaluation(Real x) const
0349     {
0350         using std::floor;
0351         Real s = (x-x0_)*inv_dx_;
0352         Real ii = floor(s);
0353         auto i = static_cast<decltype(dat_.size())>(ii);
0354
0355         Real t = s - ii;
0356         // If we had infinite precision, this would never happen.
0357         // But we don't have infinite precision.
0358         if (t == 0)
0359         {
0360             return dat_[i][0];
0361         }
0362         Real y0 = dat_[i][0];
0363         Real y1 = dat_[i+1][0];
0364         Real dy0 = dat_[i][1];
0365         Real dy1 = dat_[i+1][1];
0366
0367         Real r = 1-t;
0368         return r*r*(y0*(1+2*t) + dy0*t)
0369               + t*t*(y1*(3-2*t) - dy1*r);
0370     }
0371
0372     inline Real prime(Real x) const
0373     {
0374         const Real xf = x0_ + (dat_.size()-1)/inv_dx_;
0375         if  (x < x0_ || x > xf)
0376         {
0377             std::ostringstream oss;
0378             oss.precision(std::numeric_limits<Real>::digits10+3);
0379             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0380                 << x0_ << ", " << xf << "]";
0381             throw std::domain_error(oss.str());
0382         }
0383         if (x == xf)
0384         {
0385             return dat_.back()[1]*inv_dx_;
0386         }
0387         return this->unchecked_prime(x);
0388     }
0389
0390     inline Real unchecked_prime(Real x) const
0391     {
0392         using std::floor;
0393         Real s = (x-x0_)*inv_dx_;
0394         Real ii = floor(s);
0395         auto i = static_cast<decltype(dat_.size())>(ii);
0396         Real t = s - ii;
0397         if (t == 0)
0398         {
0399             return dat_[i][1]*inv_dx_;
0400         }
0401         Real y0 = dat_[i][0];
0402         Real dy0 = dat_[i][1];
0403         Real y1 = dat_[i+1][0];
0404         Real dy1 = dat_[i+1][1];
0405
0406         Real dy = 6*t*(1-t)*(y1 - y0)  + (3*t*t - 4*t+1)*dy0 + t*(3*t-2)*dy1;
0407         return dy*inv_dx_;
0408     }
0409
0410
0411     Size size() const
0412     {
0413         return dat_.size();
0414     }
0415
0416     int64_t bytes() const
0417     {
0418         return dat_.size()*dat_[0].size()*sizeof(Real) + sizeof(dat_) + 2*sizeof(Real);
0419     }
0420
0421     std::pair<Real, Real> domain() const
0422     {
0423         Real xf = x0_ + (dat_.size()-1)/inv_dx_;
0424         return {x0_, xf};
0425     }
0426
0427
0428 private:
0429     RandomAccessContainer dat_;
0430     Real x0_;
0431     Real inv_dx_;
0432 };
0433
0434 }
0435 }
0436 }
0437 }
0438 #endif