From 1f9e18e63b7f2bad27b65f9e6b5fad8fee8d15e3 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Wed, 29 May 2019 11:35:54 -0400 Subject: [PATCH 01/11] First pass at vector-valued interpolation in barycentric rational --- .../vector_barycentric_rational_detail.hpp | 212 ++++++++++++++++++ .../vector_barycentric_rational.hpp | 58 +++++ 2 files changed, 270 insertions(+) create mode 100644 include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp create mode 100644 include/boost/math/interpolators/vector_barycentric_rational.hpp diff --git a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp new file mode 100644 index 000000000..b775d009a --- /dev/null +++ b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp @@ -0,0 +1,212 @@ +/* + * Copyright Nick Thompson, 2019 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#ifndef BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_DETAIL_HPP +#define BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_DETAIL_HPP + +#include +#include // for std::move +#include // for std::is_sorted +#include +#include +#include +#include + +namespace boost{ namespace math{ namespace detail{ + +template +class vector_barycentric_rational_imp +{ +public: + vector_barycentric_rational_imp(std::vector&& x, std::vector&& y, size_t approximation_order = 3); + + Real operator()(Real x) const; + + Real prime(Real x) const; + + // The barycentric weights are not really that interesting; except to the unit tests! + Real weight(size_t i) const { return m_w[i]; } + + std::vector&& return_x() + { + return std::move(m_x); + } + + std::vector&& return_y() + { + return std::move(m_y); + } + +private: + + void calculate_weights(size_t approximation_order); + + std::vector m_x; + std::vector m_y; + std::vector m_w; +}; + +template +template +vector_barycentric_rational_imp::vector_barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order) +{ + std::ptrdiff_t n = std::distance(start_x, end_x); + + if (approximation_order >= (std::size_t)n) + { + throw std::domain_error("Approximation order must be < data length."); + } + + // Big sad memcpy. + m_x.resize(n); + m_y.resize(n); + for(unsigned i = 0; start_x != end_x; ++start_x, ++start_y, ++i) + { + // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy: + if(boost::math::isnan(*start_x)) + { + std::string msg = std::string("x[") + boost::lexical_cast(i) + "] is a NAN"; + throw std::domain_error(msg); + } + + if(boost::math::isnan(*start_y)) + { + std::string msg = std::string("y[") + boost::lexical_cast(i) + "] is a NAN"; + throw std::domain_error(msg); + } + + m_x[i] = *start_x; + m_y[i] = *start_y; + } + calculate_weights(approximation_order); +} + +template +barycentric_rational_imp::barycentric_rational_imp(std::vector&& x, std::vector&& y,size_t approximation_order) : m_x(std::move(x)), m_y(std::move(y)) +{ + BOOST_ASSERT_MSG(m_x.size() == m_y.size(), "There must be the same number of abscissas and ordinates."); + BOOST_ASSERT_MSG(approximation_order < m_x.size(), "Approximation order must be < data length."); + BOOST_ASSERT_MSG(std::is_sorted(m_x.begin(), m_x.end()), "The abscissas must be listed in increasing order x[0] < x[1] < ... < x[n-1]."); + calculate_weights(approximation_order); +} + +template +void barycentric_rational_imp::calculate_weights(size_t approximation_order) +{ + using std::abs; + int64_t n = m_x.size(); + m_w.resize(n, 0); + for(int64_t k = 0; k < n; ++k) + { + int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0); + int64_t i_max = k; + if (k >= n - (std::ptrdiff_t)approximation_order) + { + i_max = n - approximation_order - 1; + } + + for(int64_t i = i_min; i <= i_max; ++i) + { + Real inv_product = 1; + int64_t j_max = (std::min)(static_cast(i + approximation_order), static_cast(n - 1)); + for(int64_t j = i; j <= j_max; ++j) + { + if (j == k) + { + continue; + } + + Real diff = m_x[k] - m_x[j]; + using std::numeric_limits; + if (abs(diff) < (numeric_limits::min)()) + { + std::string msg = std::string("Spacing between x[") + + boost::lexical_cast(k) + std::string("] and x[") + + boost::lexical_cast(i) + std::string("] is ") + + boost::lexical_cast(diff) + std::string(", which is smaller than the epsilon of ") + + boost::core::demangle(typeid(Real).name()); + throw std::logic_error(msg); + } + inv_product *= diff; + } + if (i % 2 == 0) + { + m_w[k] += 1/inv_product; + } + else + { + m_w[k] -= 1/inv_product; + } + } + } +} + + +template +Real vector_barycentric_rational_imp::operator()(Real x) const +{ + Real numerator = 0; + Real denominator = 0; + for(size_t i = 0; i < m_x.size(); ++i) + { + // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality. + // However, it has been shown that if x approx x_i, but x != x_i, then inaccuracy in the numerator cancels the inaccuracy in the denominator, + // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715 + if (x == m_x[i]) + { + return m_y[i]; + } + Real t = m_w[i]/(x - m_x[i]); + numerator += t*m_y[i]; + denominator += t; + } + return numerator/denominator; +} + +/* + * A formula for computing the derivative of the barycentric representation is given in + * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner, + * Mathematics of Computation, v47, number 175, 1986. + * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf + * and reviewed in + * Recent developments in barycentric rational interpolation + * Jean–Paul Berrut, Richard Baltensperger and Hans D. Mittelmann + * + * Is it possible to complete this in one pass through the data? + */ + +template +Real vector_barycentric_rational_imp::prime(Real x) const +{ + Real rx = this->operator()(x); + Real numerator = 0; + Real denominator = 0; + for(size_t i = 0; i < m_x.size(); ++i) + { + if (x == m_x[i]) + { + Real sum = 0; + for (size_t j = 0; j < m_x.size(); ++j) + { + if (j == i) + { + continue; + } + sum += m_w[j]*(m_y[i] - m_y[j])/(m_x[i] - m_x[j]); + } + return -sum/m_w[i]; + } + Real t = m_w[i]/(x - m_x[i]); + Real diff = (rx - m_y[i])/(x-m_x[i]); + numerator += t*diff; + denominator += t; + } + + return numerator/denominator; +} +}}} +#endif diff --git a/include/boost/math/interpolators/vector_barycentric_rational.hpp b/include/boost/math/interpolators/vector_barycentric_rational.hpp new file mode 100644 index 000000000..1912b34b1 --- /dev/null +++ b/include/boost/math/interpolators/vector_barycentric_rational.hpp @@ -0,0 +1,58 @@ +/* + * Copyright Nick Thompson, 2019 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + * + * Exactly the same as barycentric_rational.hpp, but delivers values in $\mathbb{R}^n$. + * In some sense this is trivial, since each component of the vector is computed in exactly the same + * as would be computed by barycentric_rational.hpp. But this is a bit more efficient and convenient. + */ + +#ifndef BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_HPP +#define BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_HPP + +#include +#include + +namespace boost{ namespace math{ + +template +class vector_barycentric_rational +{ +public: + using Real = RandomAccessContainer2::value_type; + using Point = RandomAccessContainer1::value_type; + vector_barycentric_rational(RandomAccessContainer1&& positions, RandomAccessContainer2&& times, size_t approximation_order = 3); + + void operator()(Point& x, Real t) const; + + void prime(Point& x, Real t) const; + +private: + std::shared_ptr> m_imp; +}; + + +template +vector_barycentric_rational::vector_barycentric_rational(RandomAccessContainer1&& x, RandomAccessContainer2&& t, size_t approximation_order): + m_imp(std::make_shared>(std::move(x), std::move(y), approximation_order)) +{ + return; +} + +template +void vector_barycentric_rational::operator()(Real x) const +{ + return m_imp->operator()(x); +} + +template +void vector_barycentric_rational::prime(Real x) const +{ + return m_imp->prime(x); +} + + +}} +#endif From 2f725f0299c856a0963912176cb1293b72c83916 Mon Sep 17 00:00:00 2001 From: Nick Date: Wed, 29 May 2019 15:52:23 -0400 Subject: [PATCH 02/11] Tests for vector-valued barycentric rational. [CI SKIP] --- .../vector_barycentric_rational_detail.hpp | 151 +++----- .../vector_barycentric_rational.hpp | 44 ++- test/Jamfile.v2 | 1 + test/test_vector_barycentric_rational.cpp | 321 ++++++++++++++++++ 4 files changed, 405 insertions(+), 112 deletions(-) create mode 100644 test/test_vector_barycentric_rational.cpp diff --git a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp index b775d009a..1e5737008 100644 --- a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp +++ b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp @@ -10,23 +10,22 @@ #include #include // for std::move -#include // for std::is_sorted -#include -#include -#include #include namespace boost{ namespace math{ namespace detail{ -template +template class vector_barycentric_rational_imp { public: - vector_barycentric_rational_imp(std::vector&& x, std::vector&& y, size_t approximation_order = 3); + using Real = typename TimeContainer::value_type; + using Point = typename SpaceContainer::value_type; - Real operator()(Real x) const; + vector_barycentric_rational_imp(TimeContainer&& t, SpaceContainer&& y, size_t approximation_order); - Real prime(Real x) const; + void operator()(Point& p, Real t) const; + + /*Real prime(Real x) const; // The barycentric weights are not really that interesting; except to the unit tests! Real weight(size_t i) const { return m_w[i]; } @@ -39,67 +38,40 @@ public: std::vector&& return_y() { return std::move(m_y); - } + }*/ private: void calculate_weights(size_t approximation_order); - std::vector m_x; - std::vector m_y; - std::vector m_w; + TimeContainer t_; + SpaceContainer y_; + TimeContainer w_; }; -template -template -vector_barycentric_rational_imp::vector_barycentric_rational_imp(InputIterator1 start_x, InputIterator1 end_x, InputIterator2 start_y, size_t approximation_order) +template +vector_barycentric_rational_imp::vector_barycentric_rational_imp(TimeContainer&& t, SpaceContainer&& y, size_t approximation_order) { - std::ptrdiff_t n = std::distance(start_x, end_x); - - if (approximation_order >= (std::size_t)n) - { - throw std::domain_error("Approximation order must be < data length."); - } - - // Big sad memcpy. - m_x.resize(n); - m_y.resize(n); - for(unsigned i = 0; start_x != end_x; ++start_x, ++start_y, ++i) - { - // But if we're going to do a memcpy, we can do some error checking which is inexpensive relative to the copy: - if(boost::math::isnan(*start_x)) - { - std::string msg = std::string("x[") + boost::lexical_cast(i) + "] is a NAN"; - throw std::domain_error(msg); - } - - if(boost::math::isnan(*start_y)) - { - std::string msg = std::string("y[") + boost::lexical_cast(i) + "] is a NAN"; - throw std::domain_error(msg); - } - - m_x[i] = *start_x; - m_y[i] = *start_y; + using Real = typename TimeContainer::value_type; + using std::numeric_limits; + t_ = std::move(t); + y_ = std::move(y); + BOOST_ASSERT_MSG(t_.size() == y_.size(), "There must be the same number of time points as space points."); + BOOST_ASSERT_MSG(approximation_order < y_.size(), "Approximation order must be < data length."); + for (size_t i = 1; i < t_.size(); ++i) { + BOOST_ASSERT_MSG(t_[i] - t_[i-1] > (numeric_limits::min)(), "The abscissas must be listed in strictly increasing order t[0] < t[1] < ... < t[n-1]."); } calculate_weights(approximation_order); } -template -barycentric_rational_imp::barycentric_rational_imp(std::vector&& x, std::vector&& y,size_t approximation_order) : m_x(std::move(x)), m_y(std::move(y)) -{ - BOOST_ASSERT_MSG(m_x.size() == m_y.size(), "There must be the same number of abscissas and ordinates."); - BOOST_ASSERT_MSG(approximation_order < m_x.size(), "Approximation order must be < data length."); - BOOST_ASSERT_MSG(std::is_sorted(m_x.begin(), m_x.end()), "The abscissas must be listed in increasing order x[0] < x[1] < ... < x[n-1]."); - calculate_weights(approximation_order); -} -template -void barycentric_rational_imp::calculate_weights(size_t approximation_order) +template +void vector_barycentric_rational_imp::calculate_weights(size_t approximation_order) { + using Real = typename TimeContainer::value_type; using std::abs; - int64_t n = m_x.size(); - m_w.resize(n, 0); + int64_t n = t_.size(); + w_.resize(n, Real(0)); for(int64_t k = 0; k < n; ++k) { int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0); @@ -119,68 +91,50 @@ void barycentric_rational_imp::calculate_weights(size_t approximation_orde { continue; } - - Real diff = m_x[k] - m_x[j]; - using std::numeric_limits; - if (abs(diff) < (numeric_limits::min)()) - { - std::string msg = std::string("Spacing between x[") - + boost::lexical_cast(k) + std::string("] and x[") - + boost::lexical_cast(i) + std::string("] is ") - + boost::lexical_cast(diff) + std::string(", which is smaller than the epsilon of ") - + boost::core::demangle(typeid(Real).name()); - throw std::logic_error(msg); - } + Real diff = t_[k] - t_[j]; inv_product *= diff; } if (i % 2 == 0) { - m_w[k] += 1/inv_product; + t_[k] += 1/inv_product; } else { - m_w[k] -= 1/inv_product; + t_[k] -= 1/inv_product; } } } } -template -Real vector_barycentric_rational_imp::operator()(Real x) const +template +void vector_barycentric_rational_imp::operator()(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const { - Real numerator = 0; - Real denominator = 0; - for(size_t i = 0; i < m_x.size(); ++i) - { - // Presumably we should see if the accuracy is improved by using ULP distance of say, 5 here, instead of testing for floating point equality. - // However, it has been shown that if x approx x_i, but x != x_i, then inaccuracy in the numerator cancels the inaccuracy in the denominator, - // and the result is fairly accurate. See: http://epubs.siam.org/doi/pdf/10.1137/S0036144502417715 - if (x == m_x[i]) - { - return m_y[i]; - } - Real t = m_w[i]/(x - m_x[i]); - numerator += t*m_y[i]; - denominator += t; + using Real = typename TimeContainer::value_type; + //using Point = typename SpaceContainer::value_type; + for (auto & x : p) { + x = Real(0); } - return numerator/denominator; + Real denominator = 0; + for(size_t i = 0; i < t_.size(); ++i) + { + // See associated commentary in the scalar version of this function. + if (t == t_[i]) + { + p = y_[i]; + return; + } + Real x = w_[i]/(t - t_[i]); + p += x*y_[i]; + denominator += x; + } + p /= denominator; + return; } -/* - * A formula for computing the derivative of the barycentric representation is given in - * "Some New Aspects of Rational Interpolation", by Claus Schneider and Wilhelm Werner, - * Mathematics of Computation, v47, number 175, 1986. - * http://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf - * and reviewed in - * Recent developments in barycentric rational interpolation - * Jean–Paul Berrut, Richard Baltensperger and Hans D. Mittelmann - * - * Is it possible to complete this in one pass through the data? - */ -template -Real vector_barycentric_rational_imp::prime(Real x) const +template +void vector_barycentric_rational_imp::prime(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const { Real rx = this->operator()(x); Real numerator = 0; @@ -208,5 +162,6 @@ Real vector_barycentric_rational_imp::prime(Real x) const return numerator/denominator; } + }}} #endif diff --git a/include/boost/math/interpolators/vector_barycentric_rational.hpp b/include/boost/math/interpolators/vector_barycentric_rational.hpp index 1912b34b1..15b007145 100644 --- a/include/boost/math/interpolators/vector_barycentric_rational.hpp +++ b/include/boost/math/interpolators/vector_barycentric_rational.hpp @@ -17,40 +17,56 @@ namespace boost{ namespace math{ -template +template class vector_barycentric_rational { public: - using Real = RandomAccessContainer2::value_type; - using Point = RandomAccessContainer1::value_type; - vector_barycentric_rational(RandomAccessContainer1&& positions, RandomAccessContainer2&& times, size_t approximation_order = 3); + using Real = typename TimeContainer::value_type; + using Point = typename SpaceContainer::value_type; + vector_barycentric_rational(TimeContainer&& times, SpaceContainer&& points, size_t approximation_order = 3); void operator()(Point& x, Real t) const; + // I have validated using google benchmark that returning a value is no more expensive populating it. + // This is kinda a weird thing to discover since it goes against the advice of basically every high-performance computing book. + Point operator()(Real t) const { + Point p; + this->operator()(p, t); + return p; + } + void prime(Point& x, Real t) const; + Point prime(Real t) const { + Point p; + this->prime(p, t); + return p; + } + private: - std::shared_ptr> m_imp; + std::shared_ptr> m_imp; }; -template -vector_barycentric_rational::vector_barycentric_rational(RandomAccessContainer1&& x, RandomAccessContainer2&& t, size_t approximation_order): - m_imp(std::make_shared>(std::move(x), std::move(y), approximation_order)) +template +vector_barycentric_rational::vector_barycentric_rational(TimeContainer&& times, SpaceContainer&& points, size_t approximation_order): + m_imp(std::make_shared>(std::move(times), std::move(points), approximation_order)) { return; } -template -void vector_barycentric_rational::operator()(Real x) const +template +void vector_barycentric_rational::operator()(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const { - return m_imp->operator()(x); + m_imp->operator()(p, t); + return; } -template -void vector_barycentric_rational::prime(Real x) const +template +void vector_barycentric_rational::prime(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const { - return m_imp->prime(x); + m_imp->prime(p, t); + return; } diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 07795a320..55885bfa2 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -899,6 +899,7 @@ test-suite distribution_tests : test-suite misc : [ run test_print_info_on_type.cpp ] [ run test_barycentric_rational.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_unified_initialization_syntax ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] + [ run test_vector_barycentric_rational.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_smart_ptr cxx11_defaulted_functions cxx11_auto_declarations cxx11_unified_initialization_syntax ] [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] [ run test_constant_generate.cpp : : : release USE_CPP_FLOAT=1 off:no ] [ run test_cubic_b_spline.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_smart_ptr cxx11_defaulted_functions ] off msvc:/bigobj release ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] # does not in fact require C++17 constexpr; requires C++17 std::size. diff --git a/test/test_vector_barycentric_rational.cpp b/test/test_vector_barycentric_rational.cpp new file mode 100644 index 000000000..c369f2dd4 --- /dev/null +++ b/test/test_vector_barycentric_rational.cpp @@ -0,0 +1,321 @@ +// Copyright Nick Thompson, 2017 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_TEST_MODULE vector_barycentric_rational + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef BOOST_HAS_FLOAT128 +#include +#endif + +using std::sqrt; +using std::abs; +using std::numeric_limits; +using boost::multiprecision::cpp_bin_float_50; + +template +void test_interpolation_condition_eigen() +{ + std::cout << "Testing interpolation condition for barycentric interpolation on Eigen vectors of type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::mt19937 gen(4723); + boost::random::uniform_real_distribution dis(0.1f, 1); + std::vector t(100); + std::vector y(100); + t[0] = dis(gen); + y[0][0] = dis(gen); + y[0][1] = dis(gen); + for (size_t i = 1; i < t.size(); ++i) + { + t[i] = t[i-1] + dis(gen); + y[i][0] = dis(gen); + y[i][1] = dis(gen); + } + + boost::math::vector_barycentric_rational interpolator(std::move(t), std::move(y)); + + Eigen::Vector2d z; + for (size_t i = 0; i < t.size(); ++i) + { + interpolator(z, t[i]); + BOOST_CHECK_CLOSE(z[0], y[i][0], 100*numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(z[1], y[i][1], 100*numeric_limits::epsilon()); + } +} + +template +void test_interpolation_condition_ublas() +{ + std::cout << "Testing interpolation condition for vector barycentric interpolation ublas vectors of type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::mt19937 gen(4723); + boost::random::uniform_real_distribution dis(0.1f, 1); + std::vector t(100); + std::vector> y(100); + t[0] = dis(gen); + y[0].resize(2); + y[0][0] = dis(gen); + y[0][1] = dis(gen); + for (size_t i = 1; i < t.size(); ++i) + { + t[i] = t[i-1] + dis(gen); + y[i].resize(2); + y[i][0] = dis(gen); + y[i][1] = dis(gen); + } + + boost::math::vector_barycentric_rational interpolator(std::move(t), std::move(y)); + + boost::numeric::ublas::vector z(2); + for (size_t i = 0; i < t.size(); ++i) + { + interpolator(z, t[i]); + BOOST_CHECK_CLOSE(z[0], y[i][0], 100*numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(z[1], y[i][1], 100*numeric_limits::epsilon()); + } +} + +/*template +void test_interpolation_condition_high_order() +{ + std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::mt19937 gen(5); + boost::random::uniform_real_distribution dis(0.1f, 1); + std::vector x(100); + std::vector y(100); + x[0] = dis(gen); + y[0] = dis(gen); + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = dis(gen); + } + + // Order 5 approximation: + boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size(), 5); + + for (size_t i = 0; i < x.size(); ++i) + { + Real z = interpolator(x[i]); + BOOST_CHECK_CLOSE(z, y[i], 100*numeric_limits::epsilon()); + } +} + + +template +void test_constant() +{ + std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + + std::mt19937 gen(6); + boost::random::uniform_real_distribution dis(0.1f, 1); + std::vector x(100); + std::vector y(100); + Real constant = -8; + x[0] = dis(gen); + y[0] = constant; + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = y[0]; + } + + boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size()); + + for (size_t i = 0; i < x.size(); ++i) + { + // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test. + Real t = x[i] + dis(gen); + Real z = interpolator(t); + BOOST_CHECK_CLOSE(z, constant, 100*sqrt(numeric_limits::epsilon())); + BOOST_CHECK_SMALL(interpolator.prime(t), sqrt(numeric_limits::epsilon())); + } +} + +template +void test_constant_high_order() +{ + std::cout << "Testing that constants are interpolated correctly in high order using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + + std::mt19937 gen(7); + boost::random::uniform_real_distribution dis(0.1f, 1); + std::vector x(100); + std::vector y(100); + Real constant = 5; + x[0] = dis(gen); + y[0] = constant; + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = y[0]; + } + + // Set interpolation order to 7: + boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size(), 7); + + for (size_t i = 0; i < x.size(); ++i) + { + Real t = x[i] + dis(gen); + Real z = interpolator(t); + BOOST_CHECK_CLOSE(z, constant, 1000*sqrt(numeric_limits::epsilon())); + BOOST_CHECK_SMALL(interpolator.prime(t), 100*sqrt(numeric_limits::epsilon())); + } +} + + +template +void test_runge() +{ + std::cout << "Testing interpolation of Runge's 1/(1+25x^2) function using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + + std::mt19937 gen(8); + boost::random::uniform_real_distribution dis(0.005f, 0.01f); + std::vector x(100); + std::vector y(100); + x[0] = -2; + y[0] = 1/(1+25*x[0]*x[0]); + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = 1/(1+25*x[i]*x[i]); + } + + boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size(), 5); + + for (size_t i = 0; i < x.size(); ++i) + { + Real t = x[i]; + Real z = interpolator(t); + BOOST_CHECK_CLOSE(z, y[i], 0.03); + Real z_prime = interpolator.prime(t); + Real num = -50*t; + Real denom = (1+25*t*t)*(1+25*t*t); + if (abs(num/denom) > 0.00001) + { + BOOST_CHECK_CLOSE_FRACTION(z_prime, num/denom, 0.03); + } + } + + + Real tol = 0.0001; + for (size_t i = 0; i < x.size(); ++i) + { + Real t = x[i] + dis(gen); + Real z = interpolator(t); + BOOST_CHECK_CLOSE(z, 1/(1+25*t*t), tol); + Real z_prime = interpolator.prime(t); + Real num = -50*t; + Real denom = (1+25*t*t)*(1+25*t*t); + Real runge_prime = num/denom; + + if (abs(runge_prime) > 0 && abs(z_prime - runge_prime)/abs(runge_prime) > tol) + { + std::cout << "Error too high for t = " << t << " which is a distance " << t - x[i] << " from node " << i << "/" << x.size() << " associated with data (" << x[i] << ", " << y[i] << ")\n"; + BOOST_CHECK_CLOSE_FRACTION(z_prime, runge_prime, tol); + } + } +} + +template +void test_weights() +{ + std::cout << "Testing weights are calculated correctly using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + + std::mt19937 gen(9); + boost::random::uniform_real_distribution dis(0.005, 0.01); + std::vector x(100); + std::vector y(100); + x[0] = -2; + y[0] = 1/(1+25*x[0]*x[0]); + for (size_t i = 1; i < x.size(); ++i) + { + x[i] = x[i-1] + dis(gen); + y[i] = 1/(1+25*x[i]*x[i]); + } + + boost::math::detail::barycentric_rational_imp interpolator(x.data(), x.data() + x.size(), y.data(), 0); + + for (size_t i = 0; i < x.size(); ++i) + { + Real w = interpolator.weight(i); + if (i % 2 == 0) + { + BOOST_CHECK_CLOSE(w, 1, 0.00001); + } + else + { + BOOST_CHECK_CLOSE(w, -1, 0.00001); + } + } + + // d = 1: + interpolator = boost::math::detail::barycentric_rational_imp(x.data(), x.data() + x.size(), y.data(), 1); + + for (size_t i = 1; i < x.size() -1; ++i) + { + Real w = interpolator.weight(i); + Real w_expect = 1/(x[i] - x[i - 1]) + 1/(x[i+1] - x[i]); + if (i % 2 == 0) + { + BOOST_CHECK_CLOSE(w, -w_expect, 0.00001); + } + else + { + BOOST_CHECK_CLOSE(w, w_expect, 0.00001); + } + } + +}*/ + + +BOOST_AUTO_TEST_CASE(vector_barycentric_rational) +{ + // The tests took too long at the higher precisions. + // They still pass, but the CI system is starting to time out, + // so I figured it'd be polite to comment out the most expensive tests. + //test_weights(); + + //test_constant(); + //test_constant(); + //test_constant(); + //test_constant(); + + //test_constant_high_order(); + //test_constant_high_order(); + //test_constant_high_order(); + //test_constant_high_order(); + + //test_interpolation_condition(); + test_interpolation_condition_eigen(); + test_interpolation_condition_ublas(); + //test_interpolation_condition(); + //test_interpolation_condition(); + + //test_interpolation_condition_high_order(); + //test_interpolation_condition_high_order(); + //test_interpolation_condition_high_order(); + //test_interpolation_condition_high_order(); + + //test_runge(); + //test_runge(); + //test_runge(); + +#ifdef BOOST_HAS_FLOAT128 + //test_interpolation_condition(); + //test_constant(); + //test_constant_high_order(); + //test_interpolation_condition_high_order(); + //test_runge(); +#endif +} From 9e21a89675e683c03d1332c4bbfd2c550dd6ca8b Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Thu, 30 May 2019 11:12:38 -0400 Subject: [PATCH 03/11] fix move constructor use [CI SKIP] --- .../vector_barycentric_rational_detail.hpp | 44 +++++++++++++++++-- .../vector_barycentric_rational.hpp | 25 +++++++---- test/test_vector_barycentric_rational.cpp | 27 +++++------- 3 files changed, 68 insertions(+), 28 deletions(-) diff --git a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp index 1e5737008..c5ddfb93b 100644 --- a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp +++ b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp @@ -25,6 +25,8 @@ public: void operator()(Point& p, Real t) const; + void eval_with_prime(Point& x, Point& dxdt, Real t) const; + /*Real prime(Real x) const; // The barycentric weights are not really that interesting; except to the unit tests! @@ -111,7 +113,6 @@ template void vector_barycentric_rational_imp::operator()(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const { using Real = typename TimeContainer::value_type; - //using Point = typename SpaceContainer::value_type; for (auto & x : p) { x = Real(0); } @@ -132,8 +133,45 @@ void vector_barycentric_rational_imp::operator()( return; } - template +void vector_barycentric_rational_imp::eval_with_prime(typename SpaceContainer::value_type& x, typename SpaceContainer::value_type& dxdt, typename TimeContainer::value_type t) const +{ + using Point = typename SpaceContainer::value_type; + using Real = typename TimeContainer::value_type; + this->operator()(x, t); + Point numerator; + for (size_t i = 0; i < x.size(); ++i) { + numerator[i] = 0; + } + Real denominator = 0; + for(size_t i = 0; i < t_.size(); ++i) + { + if (t == t_[i]) + { + Real sum = 0; + for (size_t j = 0; j < t_.size(); ++j) + { + if (j == i) + { + continue; + } + sum += w_[j]*(y_[i] - y_[j])/(t_[i] - t_[j]); + } + dxdt = -sum/w_[i]; + return; + } + Real tw = w_[i]/(t - t_[i]); + Real diff = (x - y_[i])/(t-t_[i]); + numerator += tw*diff; + denominator += tw; + } + + dxdt = numerator/denominator; + return; +} + + +/*template void vector_barycentric_rational_imp::prime(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const { Real rx = this->operator()(x); @@ -161,7 +199,7 @@ void vector_barycentric_rational_imp::prime(typen } return numerator/denominator; -} +}*/ }}} #endif diff --git a/include/boost/math/interpolators/vector_barycentric_rational.hpp b/include/boost/math/interpolators/vector_barycentric_rational.hpp index 15b007145..e9ce2ef17 100644 --- a/include/boost/math/interpolators/vector_barycentric_rational.hpp +++ b/include/boost/math/interpolators/vector_barycentric_rational.hpp @@ -35,7 +35,10 @@ public: return p; } - void prime(Point& x, Real t) const; + void prime(Point& dxdt, Real t) const { + Point x; + m_imp->eval_with_prime(x, dxdt, t); + } Point prime(Real t) const { Point p; @@ -43,6 +46,18 @@ public: return p; } + void eval_with_prime(Point& x, Point& dxdt, Real t) const { + m_imp->eval_with_prime(x, dxdt, t); + return; + } + + std::pair eval_with_prime(Real t) const { + Point x; + Point dxdt; + m_imp->eval_with_prime(x, dxdt, t); + return {x, dxdt}; + } + private: std::shared_ptr> m_imp; }; @@ -62,13 +77,5 @@ void vector_barycentric_rational::operator()(type return; } -template -void vector_barycentric_rational::prime(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const -{ - m_imp->prime(p, t); - return; -} - - }} #endif diff --git a/test/test_vector_barycentric_rational.cpp b/test/test_vector_barycentric_rational.cpp index c369f2dd4..e6330e2f7 100644 --- a/test/test_vector_barycentric_rational.cpp +++ b/test/test_vector_barycentric_rational.cpp @@ -15,16 +15,11 @@ #include #include #include -#include -#ifdef BOOST_HAS_FLOAT128 -#include -#endif using std::sqrt; using std::abs; using std::numeric_limits; -using boost::multiprecision::cpp_bin_float_50; template void test_interpolation_condition_eigen() @@ -44,14 +39,16 @@ void test_interpolation_condition_eigen() y[i][1] = dis(gen); } + std::vector y_copy = y; + std::vector t_copy = t; boost::math::vector_barycentric_rational interpolator(std::move(t), std::move(y)); Eigen::Vector2d z; - for (size_t i = 0; i < t.size(); ++i) + for (size_t i = 0; i < t_copy.size(); ++i) { - interpolator(z, t[i]); - BOOST_CHECK_CLOSE(z[0], y[i][0], 100*numeric_limits::epsilon()); - BOOST_CHECK_CLOSE(z[1], y[i][1], 100*numeric_limits::epsilon()); + interpolator(z, t_copy[i]); + BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits::epsilon()); } } @@ -225,7 +222,7 @@ void test_runge() BOOST_CHECK_CLOSE_FRACTION(z_prime, runge_prime, tol); } } -} +}*/ template void test_weights() @@ -275,8 +272,7 @@ void test_weights() BOOST_CHECK_CLOSE(w, w_expect, 0.00001); } } - -}*/ +} BOOST_AUTO_TEST_CASE(vector_barycentric_rational) @@ -284,7 +280,7 @@ BOOST_AUTO_TEST_CASE(vector_barycentric_rational) // The tests took too long at the higher precisions. // They still pass, but the CI system is starting to time out, // so I figured it'd be polite to comment out the most expensive tests. - //test_weights(); + test_weights(); //test_constant(); //test_constant(); @@ -296,11 +292,10 @@ BOOST_AUTO_TEST_CASE(vector_barycentric_rational) //test_constant_high_order(); //test_constant_high_order(); - //test_interpolation_condition(); + test_interpolation_condition_eigen(); test_interpolation_condition_eigen(); test_interpolation_condition_ublas(); - //test_interpolation_condition(); - //test_interpolation_condition(); + test_interpolation_condition_ublas(); //test_interpolation_condition_high_order(); //test_interpolation_condition_high_order(); From 18feb0fc2a0c793c7312c42e181e72b93294a893 Mon Sep 17 00:00:00 2001 From: Nick Date: Thu, 30 May 2019 13:25:53 -0400 Subject: [PATCH 04/11] Documentation, more unit tests [CI SKIP] --- .../vector_barycentric_rational.qbk | 75 +++++ doc/math.qbk | 1 + .../vector_barycentric_rational_detail.hpp | 66 +---- .../vector_barycentric_rational.hpp | 5 +- test/test_vector_barycentric_rational.cpp | 259 +++++++----------- 5 files changed, 189 insertions(+), 217 deletions(-) create mode 100644 doc/interpolators/vector_barycentric_rational.qbk diff --git a/doc/interpolators/vector_barycentric_rational.qbk b/doc/interpolators/vector_barycentric_rational.qbk new file mode 100644 index 000000000..2647e9118 --- /dev/null +++ b/doc/interpolators/vector_barycentric_rational.qbk @@ -0,0 +1,75 @@ +[/ + Copyright 2019 Nick Thompson + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or copy at + http://www.boost.org/LICENSE_1_0.txt). +] + +[section:vector_barycentric Vector-valued Barycentric Rational Interpolation] + +[heading Synopsis] + +`` +#include + +namespace boost{ namespace math{ + +template +class vector_barycentric_rational +{ +public: + using Real = typename TimeContainer::value_type; + using Point = typename SpaceContainer::value_type; + vector_barycentric_rational(TimeContainer&& times, SpaceContainer&& points, size_t approximation_order = 3); + + void operator()(Point& x, Real t) const; + + Point operator()(Real t) const; + + void prime(Point& dxdt, Real t) const; + + Point prime(Real t); + + void eval_with_prime(Point& x, Point& dxdt, Real t) const; + + std::pair eval_with_prime(Real t) const; +}; + +}} +`` + +[heading Description] + +The /n/ dimensional vector-valued barycentric rational interpolator is exactly the same as /n/ scalar-valued barycentric rational interpolators. +This is provided primarily for convenience and a slight improvement in efficiency over using /n/ different rational interpolators and combining their results. + +Use of the class requires a `Point`-type with the classical vector space operations defined on it, such as multiplication by scalars and vector additions. +These requirements are satisfied by (for example) `Eigen::Vector`s and `boost::numeric::ublas::vector` classes. +The call to the constructor computes the weights: + + std::vector t(100); + std::vector y(100); + + boost::math::vector_barycentric_rational interpolant(std::move(t), std::move(y)); + +To evaluate the interpolant, use + + double t = 2.3; + Eigen::Vector2d y = interpolant(t); + +If you want to populate a vector passed into the interpolant, rather than get it returned, that syntax is supported: + + Eigen::Vector2d y; + interpolant(y, t); + +We tested this with `Eigen::Vector`s and found no performance benefit, but other `Point`-types might not be the same. + +To evaluate the derivative of the interpolant use + + auto [y, y_prime] = interpolant.eval_with_prime(x); + +Computation of the derivative requires evaluation, so if you can try to use both values at once. + + +[endsect] [/section:vector_barycentric Vector Barycentric Rational Interpolation] diff --git a/doc/math.qbk b/doc/math.qbk index 3be47a250..6a485b63a 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -659,6 +659,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [mathpart interpolation Interpolation] [include interpolators/cubic_b_spline.qbk] [include interpolators/barycentric_rational_interpolation.qbk] +[include interpolators/vector_barycentric_rational.qbk] [include interpolators/catmull_rom.qbk] [endmathpart] diff --git a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp index c5ddfb93b..58de6bb5e 100644 --- a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp +++ b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp @@ -27,20 +27,8 @@ public: void eval_with_prime(Point& x, Point& dxdt, Real t) const; - /*Real prime(Real x) const; - - // The barycentric weights are not really that interesting; except to the unit tests! - Real weight(size_t i) const { return m_w[i]; } - - std::vector&& return_x() - { - return std::move(m_x); - } - - std::vector&& return_y() - { - return std::move(m_y); - }*/ + // The barycentric weights are only interesting to the unit tests: + Real weight(size_t i) const { return w_[i]; } private: @@ -58,6 +46,7 @@ vector_barycentric_rational_imp::vector_barycentr using std::numeric_limits; t_ = std::move(t); y_ = std::move(y); + BOOST_ASSERT_MSG(t_.size() == y_.size(), "There must be the same number of time points as space points."); BOOST_ASSERT_MSG(approximation_order < y_.size(), "Approximation order must be < data length."); for (size_t i = 1; i < t_.size(); ++i) { @@ -98,11 +87,11 @@ void vector_barycentric_rational_imp::calculate_w } if (i % 2 == 0) { - t_[k] += 1/inv_product; + w_[k] += 1/inv_product; } else { - t_[k] -= 1/inv_product; + w_[k] -= 1/inv_product; } } } @@ -140,16 +129,19 @@ void vector_barycentric_rational_imp::eval_with_p using Real = typename TimeContainer::value_type; this->operator()(x, t); Point numerator; - for (size_t i = 0; i < x.size(); ++i) { + for (decltype(x.size()) i = 0; i < x.size(); ++i) { numerator[i] = 0; } Real denominator = 0; - for(size_t i = 0; i < t_.size(); ++i) + for(decltype(t_.size()) i = 0; i < t_.size(); ++i) { if (t == t_[i]) { - Real sum = 0; - for (size_t j = 0; j < t_.size(); ++j) + Point sum; + for (decltype(x.size()) i = 0; i < x.size(); ++i) { + numerator[i] = 0; + } + for (decltype(t_.size()) j = 0; j < t_.size(); ++j) { if (j == i) { @@ -161,45 +153,13 @@ void vector_barycentric_rational_imp::eval_with_p return; } Real tw = w_[i]/(t - t_[i]); - Real diff = (x - y_[i])/(t-t_[i]); + Point diff = (x - y_[i])/(t-t_[i]); numerator += tw*diff; denominator += tw; } - dxdt = numerator/denominator; return; } - -/*template -void vector_barycentric_rational_imp::prime(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const -{ - Real rx = this->operator()(x); - Real numerator = 0; - Real denominator = 0; - for(size_t i = 0; i < m_x.size(); ++i) - { - if (x == m_x[i]) - { - Real sum = 0; - for (size_t j = 0; j < m_x.size(); ++j) - { - if (j == i) - { - continue; - } - sum += m_w[j]*(m_y[i] - m_y[j])/(m_x[i] - m_x[j]); - } - return -sum/m_w[i]; - } - Real t = m_w[i]/(x - m_x[i]); - Real diff = (rx - m_y[i])/(x-m_x[i]); - numerator += t*diff; - denominator += t; - } - - return numerator/denominator; -}*/ - }}} #endif diff --git a/include/boost/math/interpolators/vector_barycentric_rational.hpp b/include/boost/math/interpolators/vector_barycentric_rational.hpp index e9ce2ef17..8289799bc 100644 --- a/include/boost/math/interpolators/vector_barycentric_rational.hpp +++ b/include/boost/math/interpolators/vector_barycentric_rational.hpp @@ -27,7 +27,8 @@ public: void operator()(Point& x, Real t) const; - // I have validated using google benchmark that returning a value is no more expensive populating it. + // I have validated using google benchmark that returning a value is no more expensive populating it, + // at least for Eigen vectors with known size at compile-time. // This is kinda a weird thing to discover since it goes against the advice of basically every high-performance computing book. Point operator()(Real t) const { Point p; @@ -50,7 +51,7 @@ public: m_imp->eval_with_prime(x, dxdt, t); return; } - + std::pair eval_with_prime(Real t) const { Point x; Point dxdt; diff --git a/test/test_vector_barycentric_rational.cpp b/test/test_vector_barycentric_rational.cpp index e6330e2f7..326c11434 100644 --- a/test/test_vector_barycentric_rational.cpp +++ b/test/test_vector_barycentric_rational.cpp @@ -1,4 +1,4 @@ -// Copyright Nick Thompson, 2017 +// Copyright Nick Thompson, 2019 // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt @@ -16,7 +16,6 @@ #include #include - using std::sqrt; using std::abs; using std::numeric_limits; @@ -72,40 +71,48 @@ void test_interpolation_condition_ublas() y[i][1] = dis(gen); } + std::vector t_copy = t; + std::vector> y_copy = y; + boost::math::vector_barycentric_rational interpolator(std::move(t), std::move(y)); boost::numeric::ublas::vector z(2); - for (size_t i = 0; i < t.size(); ++i) + for (size_t i = 0; i < t_copy.size(); ++i) { - interpolator(z, t[i]); - BOOST_CHECK_CLOSE(z[0], y[i][0], 100*numeric_limits::epsilon()); - BOOST_CHECK_CLOSE(z[1], y[i][1], 100*numeric_limits::epsilon()); + interpolator(z, t_copy[i]); + BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits::epsilon()); } } -/*template +template void test_interpolation_condition_high_order() { std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; std::mt19937 gen(5); boost::random::uniform_real_distribution dis(0.1f, 1); - std::vector x(100); - std::vector y(100); - x[0] = dis(gen); - y[0] = dis(gen); - for (size_t i = 1; i < x.size(); ++i) + std::vector t(100); + std::vector y(100); + t[0] = dis(gen); + y[0][0] = dis(gen); + y[0][1] = dis(gen); + for (size_t i = 1; i < t.size(); ++i) { - x[i] = x[i-1] + dis(gen); - y[i] = dis(gen); + t[i] = t[i-1] + dis(gen); + y[i][0] = dis(gen); + y[i][1] = dis(gen); } - // Order 5 approximation: - boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size(), 5); + std::vector y_copy = y; + std::vector t_copy = t; + boost::math::vector_barycentric_rational interpolator(std::move(t), std::move(y), 5); - for (size_t i = 0; i < x.size(); ++i) + Eigen::Vector2d z; + for (size_t i = 0; i < t_copy.size(); ++i) { - Real z = interpolator(x[i]); - BOOST_CHECK_CLOSE(z, y[i], 100*numeric_limits::epsilon()); + interpolator(z, t_copy[i]); + BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits::epsilon()); } } @@ -117,113 +124,83 @@ void test_constant() std::mt19937 gen(6); boost::random::uniform_real_distribution dis(0.1f, 1); - std::vector x(100); - std::vector y(100); - Real constant = -8; - x[0] = dis(gen); - y[0] = constant; - for (size_t i = 1; i < x.size(); ++i) + std::vector t(100); + std::vector y(100); + t[0] = dis(gen); + Real constant0 = dis(gen); + Real constant1 = dis(gen); + y[0][0] = constant0; + y[0][1] = constant1; + for (size_t i = 1; i < t.size(); ++i) { - x[i] = x[i-1] + dis(gen); - y[i] = y[0]; + t[i] = t[i-1] + dis(gen); + y[i][0] = constant0; + y[i][1] = constant1; } - boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size()); + std::vector y_copy = y; + std::vector t_copy = t; + boost::math::vector_barycentric_rational interpolator(std::move(t), std::move(y)); - for (size_t i = 0; i < x.size(); ++i) + Eigen::Vector2d z; + for (size_t i = 0; i < t_copy.size(); ++i) { // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test. - Real t = x[i] + dis(gen); - Real z = interpolator(t); - BOOST_CHECK_CLOSE(z, constant, 100*sqrt(numeric_limits::epsilon())); - BOOST_CHECK_SMALL(interpolator.prime(t), sqrt(numeric_limits::epsilon())); + Real t = t_copy[i] + dis(gen); + z = interpolator(t); + BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits::epsilon())); + BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits::epsilon())); + Eigen::Vector2d zprime = interpolator.prime(t); + Real zero_0 = zprime[0]; + Real zero_1 = zprime[1]; + BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits::epsilon())); + BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits::epsilon())); } } + template void test_constant_high_order() { - std::cout << "Testing that constants are interpolated correctly in high order using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; - std::mt19937 gen(7); + std::mt19937 gen(6); boost::random::uniform_real_distribution dis(0.1f, 1); - std::vector x(100); - std::vector y(100); - Real constant = 5; - x[0] = dis(gen); - y[0] = constant; - for (size_t i = 1; i < x.size(); ++i) + std::vector t(100); + std::vector y(100); + t[0] = dis(gen); + Real constant0 = dis(gen); + Real constant1 = dis(gen); + y[0][0] = constant0; + y[0][1] = constant1; + for (size_t i = 1; i < t.size(); ++i) { - x[i] = x[i-1] + dis(gen); - y[i] = y[0]; + t[i] = t[i-1] + dis(gen); + y[i][0] = constant0; + y[i][1] = constant1; } - // Set interpolation order to 7: - boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size(), 7); + std::vector y_copy = y; + std::vector t_copy = t; + boost::math::vector_barycentric_rational interpolator(std::move(t), std::move(y), 5); - for (size_t i = 0; i < x.size(); ++i) + Eigen::Vector2d z; + for (size_t i = 0; i < t_copy.size(); ++i) { - Real t = x[i] + dis(gen); - Real z = interpolator(t); - BOOST_CHECK_CLOSE(z, constant, 1000*sqrt(numeric_limits::epsilon())); - BOOST_CHECK_SMALL(interpolator.prime(t), 100*sqrt(numeric_limits::epsilon())); + // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test. + Real t = t_copy[i] + dis(gen); + z = interpolator(t); + BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits::epsilon())); + BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits::epsilon())); + Eigen::Vector2d zprime = interpolator.prime(t); + Real zero_0 = zprime[0]; + Real zero_1 = zprime[1]; + BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits::epsilon())); + BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits::epsilon())); } } -template -void test_runge() -{ - std::cout << "Testing interpolation of Runge's 1/(1+25x^2) function using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; - - std::mt19937 gen(8); - boost::random::uniform_real_distribution dis(0.005f, 0.01f); - std::vector x(100); - std::vector y(100); - x[0] = -2; - y[0] = 1/(1+25*x[0]*x[0]); - for (size_t i = 1; i < x.size(); ++i) - { - x[i] = x[i-1] + dis(gen); - y[i] = 1/(1+25*x[i]*x[i]); - } - - boost::math::barycentric_rational interpolator(x.data(), y.data(), y.size(), 5); - - for (size_t i = 0; i < x.size(); ++i) - { - Real t = x[i]; - Real z = interpolator(t); - BOOST_CHECK_CLOSE(z, y[i], 0.03); - Real z_prime = interpolator.prime(t); - Real num = -50*t; - Real denom = (1+25*t*t)*(1+25*t*t); - if (abs(num/denom) > 0.00001) - { - BOOST_CHECK_CLOSE_FRACTION(z_prime, num/denom, 0.03); - } - } - - - Real tol = 0.0001; - for (size_t i = 0; i < x.size(); ++i) - { - Real t = x[i] + dis(gen); - Real z = interpolator(t); - BOOST_CHECK_CLOSE(z, 1/(1+25*t*t), tol); - Real z_prime = interpolator.prime(t); - Real num = -50*t; - Real denom = (1+25*t*t)*(1+25*t*t); - Real runge_prime = num/denom; - - if (abs(runge_prime) > 0 && abs(z_prime - runge_prime)/abs(runge_prime) > tol) - { - std::cout << "Error too high for t = " << t << " which is a distance " << t - x[i] << " from node " << i << "/" << x.size() << " associated with data (" << x[i] << ", " << y[i] << ")\n"; - BOOST_CHECK_CLOSE_FRACTION(z_prime, runge_prime, tol); - } - } -}*/ - template void test_weights() { @@ -231,38 +208,26 @@ void test_weights() std::mt19937 gen(9); boost::random::uniform_real_distribution dis(0.005, 0.01); - std::vector x(100); - std::vector y(100); - x[0] = -2; - y[0] = 1/(1+25*x[0]*x[0]); - for (size_t i = 1; i < x.size(); ++i) + std::vector t(100); + std::vector y(100); + t[0] = dis(gen); + y[0][0] = dis(gen); + y[0][1] = dis(gen); + for (size_t i = 1; i < t.size(); ++i) { - x[i] = x[i-1] + dis(gen); - y[i] = 1/(1+25*x[i]*x[i]); + t[i] = t[i-1] + dis(gen); + y[i][0] = dis(gen); + y[i][1] = dis(gen); } - boost::math::detail::barycentric_rational_imp interpolator(x.data(), x.data() + x.size(), y.data(), 0); + std::vector y_copy = y; + std::vector t_copy = t; + boost::math::detail::vector_barycentric_rational_imp interpolator(std::move(t), std::move(y), 1); - for (size_t i = 0; i < x.size(); ++i) + for (size_t i = 1; i < t_copy.size() - 1; ++i) { Real w = interpolator.weight(i); - if (i % 2 == 0) - { - BOOST_CHECK_CLOSE(w, 1, 0.00001); - } - else - { - BOOST_CHECK_CLOSE(w, -1, 0.00001); - } - } - - // d = 1: - interpolator = boost::math::detail::barycentric_rational_imp(x.data(), x.data() + x.size(), y.data(), 1); - - for (size_t i = 1; i < x.size() -1; ++i) - { - Real w = interpolator.weight(i); - Real w_expect = 1/(x[i] - x[i - 1]) + 1/(x[i+1] - x[i]); + Real w_expect = 1/(t_copy[i] - t_copy[i - 1]) + 1/(t_copy[i+1] - t_copy[i]); if (i % 2 == 0) { BOOST_CHECK_CLOSE(w, -w_expect, 0.00001); @@ -277,40 +242,10 @@ void test_weights() BOOST_AUTO_TEST_CASE(vector_barycentric_rational) { - // The tests took too long at the higher precisions. - // They still pass, but the CI system is starting to time out, - // so I figured it'd be polite to comment out the most expensive tests. test_weights(); - - //test_constant(); - //test_constant(); - //test_constant(); - //test_constant(); - - //test_constant_high_order(); - //test_constant_high_order(); - //test_constant_high_order(); - //test_constant_high_order(); - - test_interpolation_condition_eigen(); + test_constant(); + test_constant_high_order(); test_interpolation_condition_eigen(); test_interpolation_condition_ublas(); - test_interpolation_condition_ublas(); - - //test_interpolation_condition_high_order(); - //test_interpolation_condition_high_order(); - //test_interpolation_condition_high_order(); - //test_interpolation_condition_high_order(); - - //test_runge(); - //test_runge(); - //test_runge(); - -#ifdef BOOST_HAS_FLOAT128 - //test_interpolation_condition(); - //test_constant(); - //test_constant_high_order(); - //test_interpolation_condition_high_order(); - //test_runge(); -#endif + test_interpolation_condition_high_order(); } From aa5a63135cea5cfb356c65ad08d6466fef800b60 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Thu, 30 May 2019 13:42:58 -0400 Subject: [PATCH 05/11] Remove typo from docs and kick off build. --- doc/interpolators/vector_barycentric_rational.qbk | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/interpolators/vector_barycentric_rational.qbk b/doc/interpolators/vector_barycentric_rational.qbk index 2647e9118..08fa9f945 100644 --- a/doc/interpolators/vector_barycentric_rational.qbk +++ b/doc/interpolators/vector_barycentric_rational.qbk @@ -50,7 +50,7 @@ The call to the constructor computes the weights: std::vector t(100); std::vector y(100); - + // initialize t and y . . . boost::math::vector_barycentric_rational interpolant(std::move(t), std::move(y)); To evaluate the interpolant, use From a97edb11395244ff7c9feeb6a8566bb661807bb6 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Tue, 4 Jun 2019 11:47:17 -0400 Subject: [PATCH 06/11] Allow use of std::array [CI SKIP] --- .../vector_barycentric_rational_detail.hpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp index 58de6bb5e..ef9845277 100644 --- a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp +++ b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp @@ -115,10 +115,14 @@ void vector_barycentric_rational_imp::operator()( return; } Real x = w_[i]/(t - t_[i]); - p += x*y_[i]; + for (size_t j = 0; j < p.size(); ++j){ + p[j] += x*y_[i][j]; + } denominator += x; } - p /= denominator; + for (size_t j = 0; j < p.size(); ++j){ + p[j] /= denominator; + } return; } @@ -147,9 +151,13 @@ void vector_barycentric_rational_imp::eval_with_p { continue; } - sum += w_[j]*(y_[i] - y_[j])/(t_[i] - t_[j]); + for (size_t k = 0; k < sum.size(); ++k) { + sum[k] += w_[j]*(y_[i][k] - y_[j][k])/(t_[i] - t_[j]); + } + } + for (size_t k = 0; k < sum.size(); ++k) { + dxdt[k] = -sum[k]/w_[i]; } - dxdt = -sum/w_[i]; return; } Real tw = w_[i]/(t - t_[i]); From ce816b3d0111e8389f5ecd2da6175d492941bb38 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 15 Jun 2019 12:59:25 -0400 Subject: [PATCH 07/11] Attempt to fix travis failure with clang++-6: The following packages have unmet dependencies: clang-6.0 : Depends: libjsoncpp0 (>= 0.6.0~rc2) but it is not installable. --- .travis.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.travis.yml b/.travis.yml index e61d94154..2a6e77425 100644 --- a/.travis.yml +++ b/.travis.yml @@ -568,6 +568,7 @@ matrix: addons: apt: packages: + - libjsoncpp0 - clang-6.0 sources: - ubuntu-toolchain-r-test @@ -579,6 +580,7 @@ matrix: addons: apt: packages: + - libjsoncpp0 - clang-6.0 sources: - ubuntu-toolchain-r-test @@ -590,6 +592,7 @@ matrix: addons: apt: packages: + - libjsoncpp0 - clang-6.0 sources: - ubuntu-toolchain-r-test @@ -601,6 +604,7 @@ matrix: addons: apt: packages: + - libjsoncpp0 - clang-6.0 sources: - ubuntu-toolchain-r-test @@ -612,6 +616,7 @@ matrix: addons: apt: packages: + - libjsoncpp0 - clang-6.0 sources: - ubuntu-toolchain-r-test From 69d5ebc0bea0c8d1ab09c5fb8b52e5dfba523e5f Mon Sep 17 00:00:00 2001 From: Nick Date: Sat, 15 Jun 2019 16:19:58 -0400 Subject: [PATCH 08/11] Second attempt to fix travis build with clang 6.0. --- .travis.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.travis.yml b/.travis.yml index 2a6e77425..298f835ba 100644 --- a/.travis.yml +++ b/.travis.yml @@ -568,7 +568,7 @@ matrix: addons: apt: packages: - - libjsoncpp0 + - libjsoncpp-dev - clang-6.0 sources: - ubuntu-toolchain-r-test @@ -580,7 +580,7 @@ matrix: addons: apt: packages: - - libjsoncpp0 + - libjsoncpp-dev - clang-6.0 sources: - ubuntu-toolchain-r-test @@ -592,7 +592,7 @@ matrix: addons: apt: packages: - - libjsoncpp0 + - libjsoncpp-dev - clang-6.0 sources: - ubuntu-toolchain-r-test @@ -604,7 +604,7 @@ matrix: addons: apt: packages: - - libjsoncpp0 + - libjsoncpp-dev - clang-6.0 sources: - ubuntu-toolchain-r-test @@ -616,7 +616,7 @@ matrix: addons: apt: packages: - - libjsoncpp0 + - libjsoncpp-dev - clang-6.0 sources: - ubuntu-toolchain-r-test From 2d770b9d32a01b02ec4b3d893ca5b6fa765024ec Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Mon, 17 Jun 2019 08:28:45 -0400 Subject: [PATCH 09/11] Vector barycentric rational: Enable compilation with containers which do not have multiplication/division/addition defined on them. --- .../vector_barycentric_rational.qbk | 9 +- .../vector_barycentric_rational_detail.hpp | 44 ++++++--- test/test_vector_barycentric_rational.cpp | 92 ++++++++++++++++++- 3 files changed, 124 insertions(+), 21 deletions(-) diff --git a/doc/interpolators/vector_barycentric_rational.qbk b/doc/interpolators/vector_barycentric_rational.qbk index 08fa9f945..335ecf372 100644 --- a/doc/interpolators/vector_barycentric_rational.qbk +++ b/doc/interpolators/vector_barycentric_rational.qbk @@ -44,14 +44,15 @@ public: The /n/ dimensional vector-valued barycentric rational interpolator is exactly the same as /n/ scalar-valued barycentric rational interpolators. This is provided primarily for convenience and a slight improvement in efficiency over using /n/ different rational interpolators and combining their results. -Use of the class requires a `Point`-type with the classical vector space operations defined on it, such as multiplication by scalars and vector additions. -These requirements are satisfied by (for example) `Eigen::Vector`s and `boost::numeric::ublas::vector` classes. +Use of the class requires a `Point`-type which has size known at compile time. +These requirements are satisfied by (for example) `Eigen::Vector2d`s and `std::array` classes. The call to the constructor computes the weights: + using boost::math::vector_barycentric_rational; std::vector t(100); std::vector y(100); - // initialize t and y . . . - boost::math::vector_barycentric_rational interpolant(std::move(t), std::move(y)); + // initialize t and y . . . + vector_barycentric_rational interpolant(std::move(t), std::move(y)); To evaluate the interpolant, use diff --git a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp index ef9845277..2f8d177a3 100644 --- a/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp +++ b/include/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp @@ -49,7 +49,8 @@ vector_barycentric_rational_imp::vector_barycentr BOOST_ASSERT_MSG(t_.size() == y_.size(), "There must be the same number of time points as space points."); BOOST_ASSERT_MSG(approximation_order < y_.size(), "Approximation order must be < data length."); - for (size_t i = 1; i < t_.size(); ++i) { + for (size_t i = 1; i < t_.size(); ++i) + { BOOST_ASSERT_MSG(t_[i] - t_[i-1] > (numeric_limits::min)(), "The abscissas must be listed in strictly increasing order t[0] < t[1] < ... < t[n-1]."); } calculate_weights(approximation_order); @@ -102,7 +103,8 @@ template void vector_barycentric_rational_imp::operator()(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const { using Real = typename TimeContainer::value_type; - for (auto & x : p) { + for (auto & x : p) + { x = Real(0); } Real denominator = 0; @@ -115,12 +117,14 @@ void vector_barycentric_rational_imp::operator()( return; } Real x = w_[i]/(t - t_[i]); - for (size_t j = 0; j < p.size(); ++j){ + for (decltype(p.size()) j = 0; j < p.size(); ++j) + { p[j] += x*y_[i][j]; } denominator += x; } - for (size_t j = 0; j < p.size(); ++j){ + for (decltype(p.size()) j = 0; j < p.size(); ++j) + { p[j] /= denominator; } return; @@ -133,7 +137,8 @@ void vector_barycentric_rational_imp::eval_with_p using Real = typename TimeContainer::value_type; this->operator()(x, t); Point numerator; - for (decltype(x.size()) i = 0; i < x.size(); ++i) { + for (decltype(x.size()) i = 0; i < x.size(); ++i) + { numerator[i] = 0; } Real denominator = 0; @@ -142,30 +147,45 @@ void vector_barycentric_rational_imp::eval_with_p if (t == t_[i]) { Point sum; - for (decltype(x.size()) i = 0; i < x.size(); ++i) { - numerator[i] = 0; + for (decltype(x.size()) i = 0; i < x.size(); ++i) + { + sum[i] = 0; } + for (decltype(t_.size()) j = 0; j < t_.size(); ++j) { if (j == i) { continue; } - for (size_t k = 0; k < sum.size(); ++k) { + for (decltype(sum.size()) k = 0; k < sum.size(); ++k) + { sum[k] += w_[j]*(y_[i][k] - y_[j][k])/(t_[i] - t_[j]); } } - for (size_t k = 0; k < sum.size(); ++k) { + for (decltype(sum.size()) k = 0; k < sum.size(); ++k) + { dxdt[k] = -sum[k]/w_[i]; } return; } Real tw = w_[i]/(t - t_[i]); - Point diff = (x - y_[i])/(t-t_[i]); - numerator += tw*diff; + Point diff; + for (decltype(diff.size()) j = 0; j < diff.size(); ++j) + { + diff[j] = (x[j] - y_[i][j])/(t-t_[i]); + } + for (decltype(diff.size()) j = 0; j < diff.size(); ++j) + { + numerator[j] += tw*diff[j]; + } denominator += tw; } - dxdt = numerator/denominator; + + for (decltype(dxdt.size()) j = 0; j < dxdt.size(); ++j) + { + dxdt[j] = numerator[j]/denominator; + } return; } diff --git a/test/test_vector_barycentric_rational.cpp b/test/test_vector_barycentric_rational.cpp index 326c11434..2ff5a4793 100644 --- a/test/test_vector_barycentric_rational.cpp +++ b/test/test_vector_barycentric_rational.cpp @@ -8,6 +8,7 @@ #include #include +#include #include #include #include @@ -23,7 +24,8 @@ using std::numeric_limits; template void test_interpolation_condition_eigen() { - std::cout << "Testing interpolation condition for barycentric interpolation on Eigen vectors of type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::cout << "Testing interpolation condition for barycentric interpolation on Eigen vectors of type " + << boost::typeindex::type_id().pretty_name() << "\n"; std::mt19937 gen(4723); boost::random::uniform_real_distribution dis(0.1f, 1); std::vector t(100); @@ -51,10 +53,44 @@ void test_interpolation_condition_eigen() } } +template +void test_interpolation_condition_std_array() +{ + std::cout << "Testing interpolation condition for barycentric interpolation on std::array vectors of type " + << boost::typeindex::type_id().pretty_name() << "\n"; + std::mt19937 gen(4723); + boost::random::uniform_real_distribution dis(0.1f, 1); + std::vector t(100); + std::vector> y(100); + t[0] = dis(gen); + y[0][0] = dis(gen); + y[0][1] = dis(gen); + for (size_t i = 1; i < t.size(); ++i) + { + t[i] = t[i-1] + dis(gen); + y[i][0] = dis(gen); + y[i][1] = dis(gen); + } + + std::vector> y_copy = y; + std::vector t_copy = t; + boost::math::vector_barycentric_rational interpolator(std::move(t), std::move(y)); + + std::array z; + for (size_t i = 0; i < t_copy.size(); ++i) + { + interpolator(z, t_copy[i]); + BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits::epsilon()); + } +} + + template void test_interpolation_condition_ublas() { - std::cout << "Testing interpolation condition for vector barycentric interpolation ublas vectors of type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::cout << "Testing interpolation condition for barycentric interpolation ublas vectors of type " + << boost::typeindex::type_id().pretty_name() << "\n"; std::mt19937 gen(4723); boost::random::uniform_real_distribution dis(0.1f, 1); std::vector t(100); @@ -118,9 +154,10 @@ void test_interpolation_condition_high_order() template -void test_constant() +void test_constant_eigen() { - std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on Eigen vectors of type " + << boost::typeindex::type_id().pretty_name() << "\n"; std::mt19937 gen(6); boost::random::uniform_real_distribution dis(0.1f, 1); @@ -159,6 +196,49 @@ void test_constant() } +template +void test_constant_std_array() +{ + std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on std::array vectors of type " + << boost::typeindex::type_id().pretty_name() << "\n"; + + std::mt19937 gen(6); + boost::random::uniform_real_distribution dis(0.1f, 1); + std::vector t(100); + std::vector> y(100); + t[0] = dis(gen); + Real constant0 = dis(gen); + Real constant1 = dis(gen); + y[0][0] = constant0; + y[0][1] = constant1; + for (size_t i = 1; i < t.size(); ++i) + { + t[i] = t[i-1] + dis(gen); + y[i][0] = constant0; + y[i][1] = constant1; + } + + std::vector> y_copy = y; + std::vector t_copy = t; + boost::math::vector_barycentric_rational interpolator(std::move(t), std::move(y)); + + std::array z; + for (size_t i = 0; i < t_copy.size(); ++i) + { + // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test. + Real t = t_copy[i] + dis(gen); + z = interpolator(t); + BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits::epsilon())); + BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits::epsilon())); + std::array zprime = interpolator.prime(t); + Real zero_0 = zprime[0]; + Real zero_1 = zprime[1]; + BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits::epsilon())); + BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits::epsilon())); + } +} + + template void test_constant_high_order() { @@ -243,9 +323,11 @@ void test_weights() BOOST_AUTO_TEST_CASE(vector_barycentric_rational) { test_weights(); - test_constant(); + test_constant_eigen(); + test_constant_std_array(); test_constant_high_order(); test_interpolation_condition_eigen(); test_interpolation_condition_ublas(); + test_interpolation_condition_std_array(); test_interpolation_condition_high_order(); } From 676a584c10bdac4b99b865695756ed4f632490cd Mon Sep 17 00:00:00 2001 From: Nick Date: Mon, 17 Jun 2019 16:40:40 -0400 Subject: [PATCH 10/11] Update .travis.yml Try Peter Dimov's suggestion to fix build from https://lists.boost.org/Archives/boost/2019/06/246359.php --- .travis.yml | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/.travis.yml b/.travis.yml index 298f835ba..5d46c76a2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -568,11 +568,10 @@ matrix: addons: apt: packages: - - libjsoncpp-dev - clang-6.0 sources: - ubuntu-toolchain-r-test - - llvm-toolchain-trusty-6.0 + - llvm-toolchain-xenial-6.0 - os: linux compiler: clang++-6.0 @@ -580,11 +579,10 @@ matrix: addons: apt: packages: - - libjsoncpp-dev - clang-6.0 sources: - ubuntu-toolchain-r-test - - llvm-toolchain-trusty-6.0 + - llvm-toolchain-xenial-6.0 - os: linux compiler: clang++-6.0 @@ -592,11 +590,10 @@ matrix: addons: apt: packages: - - libjsoncpp-dev - clang-6.0 sources: - ubuntu-toolchain-r-test - - llvm-toolchain-trusty-6.0 + - llvm-toolchain-xenial-6.0 - os: linux compiler: clang++-6.0 @@ -604,11 +601,10 @@ matrix: addons: apt: packages: - - libjsoncpp-dev - clang-6.0 sources: - ubuntu-toolchain-r-test - - llvm-toolchain-trusty-6.0 + - llvm-toolchain-xenial-6.0 - os: linux compiler: clang++-6.0 @@ -616,11 +612,10 @@ matrix: addons: apt: packages: - - libjsoncpp-dev - clang-6.0 sources: - ubuntu-toolchain-r-test - - llvm-toolchain-trusty-6.0 + - llvm-toolchain-xenial-6.0 - os: osx env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=special_fun From 442fa8c2ad1dc35991861feec6d2fa1c10579653 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Wed, 19 Jun 2019 07:21:00 -0400 Subject: [PATCH 11/11] Make sure that vector interpolant agrees with scalar interpolant. [CI SKIP] --- test/test_vector_barycentric_rational.cpp | 52 +++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/test/test_vector_barycentric_rational.cpp b/test/test_vector_barycentric_rational.cpp index 2ff5a4793..6250fa897 100644 --- a/test/test_vector_barycentric_rational.cpp +++ b/test/test_vector_barycentric_rational.cpp @@ -15,12 +15,63 @@ #include #include #include +#include #include using std::sqrt; using std::abs; using std::numeric_limits; +template +void test_agreement_with_1d() +{ + std::cout << "Testing with 1D interpolation on type " + << boost::typeindex::type_id().pretty_name() << "\n"; + std::mt19937 gen(4723); + boost::random::uniform_real_distribution dis(0.1f, 1); + std::vector t(100); + std::vector y(100); + t[0] = dis(gen); + y[0][0] = dis(gen); + y[0][1] = dis(gen); + for (size_t i = 1; i < t.size(); ++i) + { + t[i] = t[i-1] + dis(gen); + y[i][0] = dis(gen); + y[i][1] = dis(gen); + } + + std::vector y_copy = y; + std::vector t_copy = t; + std::vector t_copy0 = t; + std::vector t_copy1 = t; + + std::vector y_copy0(y.size()); + std::vector y_copy1(y.size()); + for (size_t i = 0; i < y.size(); ++i) { + y_copy0[i] = y[i][0]; + y_copy1[i] = y[i][1]; + } + + boost::random::uniform_real_distribution dis2(t[0], t[t.size()-1]); + boost::math::vector_barycentric_rational interpolator(std::move(t), std::move(y)); + boost::math::barycentric_rational scalar_interpolator0(std::move(t_copy0), std::move(y_copy0)); + boost::math::barycentric_rational scalar_interpolator1(std::move(t_copy1), std::move(y_copy1)); + + + Eigen::Vector2d z; + + size_t samples = 0; + while (samples++ < 1000) + { + Real t = dis2(gen); + interpolator(z, t); + BOOST_CHECK_CLOSE(z[0], scalar_interpolator0(t), 2*numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(z[1], scalar_interpolator1(t), 2*numeric_limits::epsilon()); + } +} + + template void test_interpolation_condition_eigen() { @@ -330,4 +381,5 @@ BOOST_AUTO_TEST_CASE(vector_barycentric_rational) test_interpolation_condition_ublas(); test_interpolation_condition_std_array(); test_interpolation_condition_high_order(); + test_agreement_with_1d(); }