123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385 |
- // 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)
- #define BOOST_TEST_MODULE vector_barycentric_rational
- #include <cmath>
- #include <random>
- #include <array>
- #include <Eigen/Dense>
- #include <boost/numeric/ublas/vector.hpp>
- #include <boost/random/uniform_real_distribution.hpp>
- #include <boost/type_index.hpp>
- #include <boost/test/included/unit_test.hpp>
- #include <boost/test/tools/floating_point_comparison.hpp>
- #include <boost/math/interpolators/barycentric_rational.hpp>
- #include <boost/math/interpolators/vector_barycentric_rational.hpp>
- using std::sqrt;
- using std::abs;
- using std::numeric_limits;
- template<class Real>
- void test_agreement_with_1d()
- {
- std::cout << "Testing with 1D interpolation on type "
- << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- std::mt19937 gen(4723);
- boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
- std::vector<Real> t(100);
- std::vector<Eigen::Vector2d> 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<Eigen::Vector2d> y_copy = y;
- std::vector<Real> t_copy = t;
- std::vector<Real> t_copy0 = t;
- std::vector<Real> t_copy1 = t;
- std::vector<Real> y_copy0(y.size());
- std::vector<Real> 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<Real> dis2(t[0], t[t.size()-1]);
- boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
- boost::math::barycentric_rational<Real> scalar_interpolator0(std::move(t_copy0), std::move(y_copy0));
- boost::math::barycentric_rational<Real> 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), 10000*numeric_limits<Real>::epsilon());
- BOOST_CHECK_CLOSE(z[1], scalar_interpolator1(t), 10000*numeric_limits<Real>::epsilon());
- }
- }
- template<class Real>
- void test_interpolation_condition_eigen()
- {
- std::cout << "Testing interpolation condition for barycentric interpolation on Eigen vectors of type "
- << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- std::mt19937 gen(4723);
- boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
- std::vector<Real> t(100);
- std::vector<Eigen::Vector2d> 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<Eigen::Vector2d> y_copy = y;
- std::vector<Real> t_copy = t;
- boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
- Eigen::Vector2d 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<Real>::epsilon());
- BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
- }
- }
- template<class Real>
- void test_interpolation_condition_std_array()
- {
- std::cout << "Testing interpolation condition for barycentric interpolation on std::array vectors of type "
- << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- std::mt19937 gen(4723);
- boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
- std::vector<Real> t(100);
- std::vector<std::array<Real, 2>> 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<std::array<Real, 2>> y_copy = y;
- std::vector<Real> t_copy = t;
- boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
- std::array<Real, 2> 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<Real>::epsilon());
- BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
- }
- }
- template<class Real>
- void test_interpolation_condition_ublas()
- {
- std::cout << "Testing interpolation condition for barycentric interpolation ublas vectors of type "
- << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- std::mt19937 gen(4723);
- boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
- std::vector<Real> t(100);
- std::vector<boost::numeric::ublas::vector<Real>> 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);
- }
- std::vector<Real> t_copy = t;
- std::vector<boost::numeric::ublas::vector<Real>> y_copy = y;
- boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
- boost::numeric::ublas::vector<Real> z(2);
- 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<Real>::epsilon());
- BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
- }
- }
- template<class Real>
- void test_interpolation_condition_high_order()
- {
- std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- std::mt19937 gen(5);
- boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
- std::vector<Real> t(100);
- std::vector<Eigen::Vector2d> 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<Eigen::Vector2d> y_copy = y;
- std::vector<Real> t_copy = t;
- boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 5);
- Eigen::Vector2d 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<Real>::epsilon());
- BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
- }
- }
- template<class Real>
- void test_constant_eigen()
- {
- std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on Eigen vectors of type "
- << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- std::mt19937 gen(6);
- boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
- std::vector<Real> t(100);
- std::vector<Eigen::Vector2d> 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<Eigen::Vector2d> y_copy = y;
- std::vector<Real> t_copy = t;
- boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
- 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 = t_copy[i] + dis(gen);
- z = interpolator(t);
- BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon()));
- BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::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<Real>::epsilon()));
- BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon()));
- }
- }
- template<class Real>
- 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<Real>().pretty_name() << "\n";
- std::mt19937 gen(6);
- boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
- std::vector<Real> t(100);
- std::vector<std::array<Real, 2>> 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<std::array<Real,2>> y_copy = y;
- std::vector<Real> t_copy = t;
- boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
- std::array<Real, 2> 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<Real>::epsilon()));
- BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::epsilon()));
- std::array<Real, 2> zprime = interpolator.prime(t);
- Real zero_0 = zprime[0];
- Real zero_1 = zprime[1];
- BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits<Real>::epsilon()));
- BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon()));
- }
- }
- template<class Real>
- void test_constant_high_order()
- {
- std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- std::mt19937 gen(6);
- boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
- std::vector<Real> t(100);
- std::vector<Eigen::Vector2d> 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<Eigen::Vector2d> y_copy = y;
- std::vector<Real> t_copy = t;
- boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 5);
- 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 = t_copy[i] + dis(gen);
- z = interpolator(t);
- BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon()));
- BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::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<Real>::epsilon()));
- BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon()));
- }
- }
- template<class Real>
- void test_weights()
- {
- std::cout << "Testing weights are calculated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- std::mt19937 gen(9);
- boost::random::uniform_real_distribution<Real> dis(0.005, 0.01);
- std::vector<Real> t(100);
- std::vector<Eigen::Vector2d> 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<Eigen::Vector2d> y_copy = y;
- std::vector<Real> t_copy = t;
- boost::math::detail::vector_barycentric_rational_imp<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 1);
- for (size_t i = 1; i < t_copy.size() - 1; ++i)
- {
- Real w = interpolator.weight(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);
- }
- else
- {
- BOOST_CHECK_CLOSE(w, w_expect, 0.00001);
- }
- }
- }
- BOOST_AUTO_TEST_CASE(vector_barycentric_rational)
- {
- test_weights<double>();
- test_constant_eigen<double>();
- test_constant_std_array<double>();
- test_constant_high_order<double>();
- test_interpolation_condition_eigen<double>();
- test_interpolation_condition_ublas<double>();
- test_interpolation_condition_std_array<double>();
- test_interpolation_condition_high_order<double>();
- test_agreement_with_1d<double>();
- }
|