diff --git a/include/boost/math/interpolators/cubic_hermite.hpp b/include/boost/math/interpolators/cubic_hermite.hpp index 48346ab1ec..09081c6ba5 100644 --- a/include/boost/math/interpolators/cubic_hermite.hpp +++ b/include/boost/math/interpolators/cubic_hermite.hpp @@ -8,18 +8,19 @@ #define BOOST_MATH_INTERPOLATORS_CUBIC_HERMITE_HPP #include #include +#include namespace boost { namespace math { namespace interpolators { -template +template> class cubic_hermite { public: using Real = typename RandomAccessContainer::value_type; cubic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx) - : impl_(std::make_shared>(std::move(x), std::move(y), std::move(dydx))) + : impl_(std::make_shared>(std::move(x), std::move(y), std::move(dydx))) {} inline Real operator()(Real x) const { @@ -43,6 +44,7 @@ class cubic_hermite { int64_t bytes() const { + if ( ! valid()) return 0; return impl_->bytes() + sizeof(impl_); } @@ -51,17 +53,25 @@ class cubic_hermite { return impl_->domain(); } + bool valid() const { + return impl_->valid(); + } + + std::string const& error_msg() const { + return impl_->error_msg(); + } + private: std::shared_ptr> impl_; }; -template +template> class cardinal_cubic_hermite { public: using Real = typename RandomAccessContainer::value_type; cardinal_cubic_hermite(RandomAccessContainer && y, RandomAccessContainer && dydx, Real x0, Real dx) - : impl_(std::make_shared>(std::move(y), std::move(dydx), x0, dx)) + : impl_(std::make_shared>(std::move(y), std::move(dydx), x0, dx)) {} inline Real operator()(Real x) const @@ -82,6 +92,7 @@ class cardinal_cubic_hermite { int64_t bytes() const { + if ( ! valid()) return 0; return impl_->bytes() + sizeof(impl_); } @@ -90,19 +101,26 @@ class cardinal_cubic_hermite { return impl_->domain(); } + bool valid() const { + return impl_->valid(); + } + + std::string const& error_msg() const { + return impl_->error_msg(); + } + private: std::shared_ptr> impl_; }; - -template +template> class cardinal_cubic_hermite_aos { public: using Point = typename RandomAccessContainer::value_type; using Real = typename Point::value_type; cardinal_cubic_hermite_aos(RandomAccessContainer && data, Real x0, Real dx) - : impl_(std::make_shared>(std::move(data), x0, dx)) + : impl_(std::make_shared>(std::move(data), x0, dx)) {} inline Real operator()(Real x) const @@ -123,6 +141,7 @@ class cardinal_cubic_hermite_aos { int64_t bytes() const { + if ( ! valid()) return 0; return impl_->bytes() + sizeof(impl_); } @@ -131,6 +150,14 @@ class cardinal_cubic_hermite_aos { return impl_->domain(); } + bool valid() const { + return impl_->valid(); + } + + std::string const& error_msg() const { + return impl_->error_msg(); + } + private: std::shared_ptr> impl_; }; diff --git a/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp b/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp index 921902e1f2..048dbbf056 100644 --- a/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp +++ b/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp @@ -13,12 +13,14 @@ #include #include +#include + namespace boost { namespace math { namespace interpolators { namespace detail { -template +template> class cubic_hermite_detail { public: using Real = typename RandomAccessContainer::value_type; @@ -27,19 +29,29 @@ class cubic_hermite_detail { cubic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer dydx) : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)} { + static const char* function = "boost::math::interpolators::detail::cubic_hermite_detail::cubic_hermite_detail"; using std::abs; using std::isnan; if (x_.size() != y_.size()) { - throw std::domain_error("There must be the same number of ordinates as abscissas."); + error_msg_ = "There must be the same number of ordinates as abscissas."; + boost::math::policies::raise_domain_error(function, error_msg_.c_str(), false, Policy()); + valid_ = false; + return; } if (x_.size() != dydx_.size()) { - throw std::domain_error("There must be the same number of ordinates as derivative values."); + error_msg_ = "There must be the same number of ordinates as derivative values."; + boost::math::policies::raise_domain_error(function, error_msg_.c_str(), false, Policy()); + valid_ = false; + return; } if (x_.size() < 2) { - throw std::domain_error("Must be at least two data points."); + error_msg_ = "Must be at least two data points."; + boost::math::policies::raise_domain_error(function, error_msg_.c_str(), false, Policy()); + valid_ = false; + return; } Real x0 = x_[0]; for (size_t i = 1; i < x_.size(); ++i) @@ -51,19 +63,25 @@ class cubic_hermite_detail { oss.precision(std::numeric_limits::digits10+3); oss << "Abscissas must be listed in strictly increasing order x0 < x1 < ... < x_{n-1}, "; oss << "but at x[" << i - 1 << "] = " << x0 << ", and x[" << i << "] = " << x1 << ".\n"; - throw std::domain_error(oss.str()); + error_msg_ = oss.str(); + valid_ = false; + return; } x0 = x1; } + valid_ = true; } void push_back(Real x, Real y, Real dydx) { + if ( ! valid_) return; // We do not realize if the object is constructed correctly here. + + static const char* function = "boost::math::interpolators::detail::cubic_hermite_detail::push_back"; using std::abs; using std::isnan; if (x <= x_.back()) { - throw std::domain_error("Calling push_back must preserve the monotonicity of the x's"); + return boost::math::policies::raise_domain_error(function, "Calling push_back must preserve the monotonicity of the x's", false, Policy()); } x_.push_back(x); y_.push_back(y); @@ -72,13 +90,16 @@ class cubic_hermite_detail { Real operator()(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); + + static const char* function = "boost::math::interpolators::detail::cubic_hermite_detail::operator()"; if (x < x_[0] || x > x_.back()) { std::ostringstream oss; oss.precision(std::numeric_limits::digits10+3); oss << "Requested abscissa x = " << x << ", which is outside of allowed range [" << x_[0] << ", " << x_.back() << "]"; - throw std::domain_error(oss.str()); + return boost::math::policies::raise_domain_error(function, oss.str().c_str(), x, Policy()); } // We need t := (x-x_k)/(x_{k+1}-x_k) \in [0,1) for this to work. // Sadly this neccessitates this loathesome check, otherwise we get t = 1 at x = xf. @@ -107,13 +128,16 @@ class cubic_hermite_detail { Real prime(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); + + static const char* function = "boost::math::interpolators::detail::cubic_hermite_detail::prime"; if (x < x_[0] || x > x_.back()) { std::ostringstream oss; oss.precision(std::numeric_limits::digits10+3); oss << "Requested abscissa x = " << x << ", which is outside of allowed range [" << x_[0] << ", " << x_.back() << "]"; - throw std::domain_error(oss.str()); + return boost::math::policies::raise_domain_error(function, oss.str().c_str(), x, Policy()); } if (x == x_.back()) { @@ -136,9 +160,10 @@ class cubic_hermite_detail { return s0 + 2*c2*(x-x0) + 3*c3*(x-x0)*(x-x0); } - friend std::ostream& operator<<(std::ostream & os, const cubic_hermite_detail & m) { + if ( ! m.valid_) return os; + os << "(x,y,y') = {"; for (size_t i = 0; i < m.x_.size() - 1; ++i) { @@ -151,25 +176,39 @@ class cubic_hermite_detail { Size size() const { + if ( ! valid_) return 0; return x_.size(); } int64_t bytes() const { + if ( ! valid_) return 0; return 3*x_.size()*sizeof(Real) + 3*sizeof(x_); } std::pair domain() const { + if ( ! valid_) return {0, 0}; return {x_.front(), x_.back()}; } + bool valid() const { + return valid_; + } + + std::string const& error_msg() const { + return error_msg_; + } + RandomAccessContainer x_; RandomAccessContainer y_; RandomAccessContainer dydx_; +private: + bool valid_ = false; + std::string error_msg_; }; -template +template> class cardinal_cubic_hermite_detail { public: using Real = typename RandomAccessContainer::value_type; @@ -178,25 +217,36 @@ class cardinal_cubic_hermite_detail { cardinal_cubic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer dydx, Real x0, Real dx) : y_{std::move(y)}, dy_{std::move(dydx)}, x0_{x0}, inv_dx_{1/dx} { + static const char* function = "boost::math::interpolators::detail::cardinal_cubic_hermite_detail::cardinal_cubic_hermite_detail"; using std::abs; using std::isnan; if (y_.size() != dy_.size()) { - throw std::domain_error("There must be the same number of derivatives as ordinates."); + error_msg_ = "There must be the same number of derivatives as ordinates."; + boost::math::policies::raise_domain_error(function, error_msg_.c_str(), false, Policy()); + valid_ = false; + return; } if (y_.size() < 2) { - throw std::domain_error("Must be at least two data points."); + error_msg_ = "Must be at least two data points."; + boost::math::policies::raise_domain_error(function, error_msg_.c_str(), false, Policy()); + valid_ = false; + return; } if (dx <= 0) { - throw std::domain_error("dx > 0 is required."); + error_msg_ = "dx > 0 is required."; + boost::math::policies::raise_domain_error(function, error_msg_.c_str(), false, Policy()); + valid_ = false; + return; } for (auto & dy : dy_) { dy *= dx; } + valid_ = true; } // Why not implement push_back? It's awkward: If the buffer is circular, x0_ += dx_. @@ -205,6 +255,9 @@ class cardinal_cubic_hermite_detail { inline Real operator()(Real x) const { + static const char* function = "boost::math::interpolators::detail::cardinal_cubic_hermite_detail::operator()"; + if ( ! valid_) return std::numeric_limits::quiet_NaN(); + const Real xf = x0_ + (y_.size()-1)/inv_dx_; if (x < x0_ || x > xf) { @@ -212,7 +265,7 @@ class cardinal_cubic_hermite_detail { oss.precision(std::numeric_limits::digits10+3); oss << "Requested abscissa x = " << x << ", which is outside of allowed range [" << x0_ << ", " << xf << "]"; - throw std::domain_error(oss.str()); + return boost::math::policies::raise_domain_error(function, oss.str().c_str(), x, Policy()); } if (x == xf) { @@ -223,6 +276,7 @@ class cardinal_cubic_hermite_detail { inline Real unchecked_evaluation(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); using std::floor; Real s = (x-x0_)*inv_dx_; Real ii = floor(s); @@ -240,6 +294,8 @@ class cardinal_cubic_hermite_detail { inline Real prime(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); + static const char* function = "boost::math::interpolators::detail::cardinal_cubic_hermite_detail::prime"; const Real xf = x0_ + (y_.size()-1)/inv_dx_; if (x < x0_ || x > xf) { @@ -247,7 +303,7 @@ class cardinal_cubic_hermite_detail { oss.precision(std::numeric_limits::digits10+3); oss << "Requested abscissa x = " << x << ", which is outside of allowed range [" << x0_ << ", " << xf << "]"; - throw std::domain_error(oss.str()); + return boost::math::policies::raise_domain_error(function, oss.str().c_str(), x, Policy()); } if (x == xf) { @@ -258,6 +314,8 @@ class cardinal_cubic_hermite_detail { inline Real unchecked_prime(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); + using std::floor; Real s = (x-x0_)*inv_dx_; Real ii = floor(s); @@ -272,33 +330,44 @@ class cardinal_cubic_hermite_detail { return dy*inv_dx_; } - Size size() const { + if ( ! valid_) return 0; return y_.size(); } int64_t bytes() const { + if ( ! valid_) return 0; return 2*y_.size()*sizeof(Real) + 2*sizeof(y_) + 2*sizeof(Real); } std::pair domain() const { + if ( ! valid_) return {0, 0}; Real xf = x0_ + (y_.size()-1)/inv_dx_; return {x0_, xf}; } -private: + bool valid() const { + return valid_; + } + std::string const& error_msg() const { + return error_msg_; + } +private: RandomAccessContainer y_; RandomAccessContainer dy_; Real x0_; Real inv_dx_; + + bool valid_ = false; + std::string error_msg_; }; -template +template> class cardinal_cubic_hermite_detail_aos { public: using Point = typename RandomAccessContainer::value_type; @@ -306,29 +375,43 @@ class cardinal_cubic_hermite_detail_aos { using Size = typename RandomAccessContainer::size_type; cardinal_cubic_hermite_detail_aos(RandomAccessContainer && dat, Real x0, Real dx) - : dat_{std::move(dat)}, x0_{x0}, inv_dx_{1/dx} + : dat_{std::move(dat)}, x0_{x0}, inv_dx_{1/dx} { + static const char* function = "boost::math::interpolators::detail::cardinal_cubic_hermite_detail_aos::cardinal_cubic_hermite_detail_aos"; if (dat_.size() < 2) { - throw std::domain_error("Must be at least two data points."); + error_msg_ = "Must be at least two data points."; + boost::math::policies::raise_domain_error(function, error_msg_.c_str(), false, Policy()); + valid_ = false; + return; } if (dat_[0].size() != 2) { - throw std::domain_error("Each datum must contain (y, y'), and nothing else."); + error_msg_ = "Each datum must contain (y, y'), and nothing else."; + boost::math::policies::raise_domain_error(function, error_msg_.c_str(), false, Policy()); + valid_ = false; + return; } if (dx <= 0) { - throw std::domain_error("dx > 0 is required."); + error_msg_ = "dx > 0 is required."; + boost::math::policies::raise_domain_error(function, error_msg_.c_str(), false, Policy()); + valid_ = false; + return; } for (auto & d : dat_) { d[1] *= dx; } + valid_ = true; } inline Real operator()(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); + + static const char* function = "boost::math::interpolators::detail::cardinal_cubic_hermite_detail_aos::operator()"; const Real xf = x0_ + (dat_.size()-1)/inv_dx_; if (x < x0_ || x > xf) { @@ -336,7 +419,7 @@ class cardinal_cubic_hermite_detail_aos { oss.precision(std::numeric_limits::digits10+3); oss << "Requested abscissa x = " << x << ", which is outside of allowed range [" << x0_ << ", " << xf << "]"; - throw std::domain_error(oss.str()); + return boost::math::policies::raise_domain_error(function, oss.str().c_str(), x, Policy()); } if (x == xf) { @@ -347,6 +430,8 @@ class cardinal_cubic_hermite_detail_aos { inline Real unchecked_evaluation(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); + using std::floor; Real s = (x-x0_)*inv_dx_; Real ii = floor(s); @@ -371,6 +456,9 @@ class cardinal_cubic_hermite_detail_aos { inline Real prime(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); + + static const char* function = "boost::math::interpolators::detail::cardinal_cubic_hermite_detail_aos::prime"; const Real xf = x0_ + (dat_.size()-1)/inv_dx_; if (x < x0_ || x > xf) { @@ -378,7 +466,7 @@ class cardinal_cubic_hermite_detail_aos { oss.precision(std::numeric_limits::digits10+3); oss << "Requested abscissa x = " << x << ", which is outside of allowed range [" << x0_ << ", " << xf << "]"; - throw std::domain_error(oss.str()); + return boost::math::policies::raise_domain_error(function, oss.str().c_str(), x, Policy()); } if (x == xf) { @@ -389,6 +477,8 @@ class cardinal_cubic_hermite_detail_aos { inline Real unchecked_prime(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); + using std::floor; Real s = (x-x0_)*inv_dx_; Real ii = floor(s); @@ -407,28 +497,40 @@ class cardinal_cubic_hermite_detail_aos { return dy*inv_dx_; } - Size size() const { + if ( ! valid_) return 0; return dat_.size(); } int64_t bytes() const { + if ( ! valid_) return 0; return dat_.size()*dat_[0].size()*sizeof(Real) + sizeof(dat_) + 2*sizeof(Real); } std::pair domain() const { + if ( ! valid_) return {0, 0}; Real xf = x0_ + (dat_.size()-1)/inv_dx_; return {x0_, xf}; } + bool valid() const { + return valid_; + } + + std::string const& error_msg() const { + return error_msg_; + } private: RandomAccessContainer dat_; Real x0_; Real inv_dx_; + + bool valid_ = false; + std::string error_msg_; }; } @@ -436,3 +538,4 @@ class cardinal_cubic_hermite_detail_aos { } } #endif + diff --git a/include/boost/math/interpolators/makima.hpp b/include/boost/math/interpolators/makima.hpp index f69fe40733..43b854fb01 100644 --- a/include/boost/math/interpolators/makima.hpp +++ b/include/boost/math/interpolators/makima.hpp @@ -17,7 +17,7 @@ namespace boost { namespace math { namespace interpolators { -template +template> class makima { public: using Real = typename RandomAccessContainer::value_type; @@ -26,11 +26,15 @@ class makima { Real left_endpoint_derivative = std::numeric_limits::quiet_NaN(), Real right_endpoint_derivative = std::numeric_limits::quiet_NaN()) { + static const char* function = "boost::math::interpolators::makima::makima"; using std::isnan; using std::abs; if (x.size() < 4) { - throw std::domain_error("Must be at least four data points."); + error_msg_ = "Must be at least four data points."; + boost::math::policies::raise_domain_error(function, error_msg_.c_str(), false, Policy()); + valid_ = false; + return; } RandomAccessContainer s(x.size(), std::numeric_limits::quiet_NaN()); Real m2 = (y[3]-y[2])/(x[3]-x[2]); @@ -54,14 +58,12 @@ class makima { { s[0] = left_endpoint_derivative; } - w1 = abs(m2-m1) + abs(m2+m1)/2; w2 = abs(m0-mm1) + abs(m0+mm1)/2; s[1] = (w1*m0 + w2*m1)/(w1+w2); if (isnan(s[1])) { s[1] = 0; } - for (decltype(s.size()) i = 2; i < s.size()-2; ++i) { Real mim2 = (y[i-1]-y[i-2])/(x[i-1]-x[i-2]); Real mim1 = (y[i ]-y[i-1])/(x[i ]-x[i-1]); @@ -84,16 +86,12 @@ class makima { Real mn = 2*mnm1 - mnm2; w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2; w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2; - s[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2); if (isnan(s[n-2])) { s[n-2] = 0; } - w1 = abs(mn - mnm1) + abs(mn+mnm1)/2; w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2; - - if (isnan(right_endpoint_derivative)) { s[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2); @@ -107,27 +105,36 @@ class makima { } impl_ = std::make_shared>(std::move(x), std::move(y), std::move(s)); + valid_ = impl_->valid(); + if ( ! valid_) { + error_msg_ = impl_->error_msg(); + } } Real operator()(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); return impl_->operator()(x); } Real prime(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); return impl_->prime(x); } friend std::ostream& operator<<(std::ostream & os, const makima & m) { + if ( ! m.valid_) return os; os << *m.impl_; return os; } void push_back(Real x, Real y) { + if ( ! valid_) return; // We do not realize if the object is constructed properly here. + static const char* function = "boost::math::interpolators::makima::push_back"; using std::abs; using std::isnan; if (x <= impl_->x_.back()) { - throw std::domain_error("Calling push_back must preserve the monotonicity of the x's"); + boost::math::policies::raise_domain_error(function, "Calling push_back must preserve the monotonicity of the x's", x, Policy()); } impl_->x_.push_back(x); impl_->y_.push_back(y); @@ -170,6 +177,8 @@ class makima { private: std::shared_ptr> impl_; + bool valid_ = false; + std::string error_msg_; }; } diff --git a/include/boost/math/interpolators/pchip.hpp b/include/boost/math/interpolators/pchip.hpp index 6697166df0..efe0f7328c 100644 --- a/include/boost/math/interpolators/pchip.hpp +++ b/include/boost/math/interpolators/pchip.hpp @@ -14,7 +14,7 @@ namespace boost { namespace math { namespace interpolators { -template +template> class pchip { public: using Real = typename RandomAccessContainer::value_type; @@ -29,7 +29,9 @@ class pchip { std::ostringstream oss; oss << __FILE__ << ":" << __LINE__ << ":" << __func__; oss << " This interpolator requires at least four data points."; - throw std::domain_error(oss.str()); + error_msg_ = oss.str(); + valid_ = false; + return; } RandomAccessContainer s(x.size(), std::numeric_limits::quiet_NaN()); if (isnan(left_endpoint_derivative)) @@ -43,11 +45,9 @@ class pchip { { s[0] = left_endpoint_derivative; } - for (decltype(s.size()) k = 1; k < s.size()-1; ++k) { Real hkm1 = x[k] - x[k-1]; Real dkm1 = (y[k] - y[k-1])/hkm1; - Real hk = x[k+1] - x[k]; Real dk = (y[k+1] - y[k])/hk; Real w1 = 2*hk + hkm1; @@ -63,7 +63,6 @@ class pchip { // Un-numbered equation just before Section 3.5: s[k] = (w1+w2)/(w1/dkm1 + w2/dk); } - } auto n = s.size(); if (isnan(right_endpoint_derivative)) @@ -75,27 +74,36 @@ class pchip { s[n-1] = right_endpoint_derivative; } impl_ = std::make_shared>(std::move(x), std::move(y), std::move(s)); + valid_ = impl_->valid(); + if ( ! valid_) { + error_msg_ = impl_->error_msg(); + } } Real operator()(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); return impl_->operator()(x); } Real prime(Real x) const { + if ( ! valid_) return std::numeric_limits::quiet_NaN(); return impl_->prime(x); } friend std::ostream& operator<<(std::ostream & os, const pchip & m) { + if ( ! m.valid_) return os; os << *m.impl_; return os; } void push_back(Real x, Real y) { + if ( ! valid_) return; // We do not realize if the object is constructed properly here. + static const char* function = "boost::math::interpolators::pchip::push_back"; using std::abs; using std::isnan; if (x <= impl_->x_.back()) { - throw std::domain_error("Calling push_back must preserve the monotonicity of the x's"); + boost::math::policies::raise_domain_error(function, "Calling push_back must preserve the monotonicity of the x's", x, Policy()); } impl_->x_.push_back(x); impl_->y_.push_back(y); @@ -123,6 +131,8 @@ class pchip { private: std::shared_ptr> impl_; + bool valid_ = false; + std::string error_msg_; }; }