123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142 |
- /*
- * 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)
- */
- #include "math_unit_test.hpp"
- #include <numeric>
- #include <utility>
- #include <random>
- #include <boost/core/demangle.hpp>
- #include <boost/math/interpolators/whittaker_shannon.hpp>
- #ifdef BOOST_HAS_FLOAT128
- #include <boost/multiprecision/float128.hpp>
- using boost::multiprecision::float128;
- #endif
- using boost::math::interpolators::whittaker_shannon;
- template<class Real>
- void test_trivial()
- {
- Real t0 = 0;
- Real h = Real(1)/Real(16);
- std::vector<Real> v{1.5};
- std::vector<Real> v_copy = v;
- auto ws = whittaker_shannon<decltype(v)>(std::move(v), t0, h);
- Real expected = 0;
- if(!CHECK_MOLLIFIED_CLOSE(expected, ws.prime(0), 10*std::numeric_limits<Real>::epsilon())) {
- std::cerr << " Problem occured at abscissa " << 0 << "\n";
- }
- expected = -v_copy[0]/h;
- if(!CHECK_MOLLIFIED_CLOSE(expected, ws.prime(h), 10*std::numeric_limits<Real>::epsilon())) {
- std::cerr << " Problem occured at abscissa " << 0 << "\n";
- }
- }
- template<class Real>
- void test_knots()
- {
- Real t0 = 0;
- Real h = Real(1)/Real(16);
- size_t n = 512;
- std::vector<Real> v(n);
- std::mt19937 gen(323723);
- std::uniform_real_distribution<Real> dis(1.0, 2.0);
- for(size_t i = 0; i < n; ++i) {
- v[i] = static_cast<Real>(dis(gen));
- }
- auto ws = whittaker_shannon<decltype(v)>(std::move(v), t0, h);
- size_t i = 0;
- while (i < n) {
- Real t = t0 + i*h;
- Real expected = ws[i];
- Real computed = ws(t);
- CHECK_ULP_CLOSE(expected, computed, 16);
- ++i;
- }
- }
- template<class Real>
- void test_bump()
- {
- using std::exp;
- using std::abs;
- using std::sqrt;
- auto bump = [](Real x) { if (abs(x) >= 1) { return Real(0); } return exp(-Real(1)/(Real(1)-x*x)); };
- auto bump_prime = [&bump](Real x) { Real z = 1-x*x; return -2*x*bump(x)/(z*z); };
- Real t0 = -1;
- size_t n = 2049;
- Real h = Real(2)/Real(n-1);
- std::vector<Real> v(n);
- for(size_t i = 0; i < n; ++i) {
- Real t = t0 + i*h;
- v[i] = bump(t);
- }
- std::vector<Real> v_copy = v;
- auto ws = whittaker_shannon<decltype(v)>(std::move(v), t0, h);
- // Test the knots:
- for(size_t i = v_copy.size()/4; i < 3*v_copy.size()/4; ++i) {
- Real t = t0 + i*h;
- Real expected = v_copy[i];
- Real computed = ws(t);
- if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 10*std::numeric_limits<Real>::epsilon())) {
- std::cerr << " Problem occured at abscissa " << t << "\n";
- }
- Real expected_prime = bump_prime(t);
- Real computed_prime = ws.prime(t);
- if(!CHECK_MOLLIFIED_CLOSE(expected_prime, computed_prime, 1000*std::numeric_limits<Real>::epsilon())) {
- std::cerr << " Problem occured at abscissa " << t << "\n";
- }
- }
- std::mt19937 gen(323723);
- std::uniform_real_distribution<long double> dis(-0.85, 0.85);
- size_t i = 0;
- while (i++ < 1000)
- {
- Real t = static_cast<Real>(dis(gen));
- Real expected = bump(t);
- Real computed = ws(t);
- if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 10*std::numeric_limits<Real>::epsilon())) {
- std::cerr << " Problem occured at abscissa " << t << "\n";
- }
- Real expected_prime = bump_prime(t);
- Real computed_prime = ws.prime(t);
- if(!CHECK_MOLLIFIED_CLOSE(expected_prime, computed_prime, sqrt(std::numeric_limits<Real>::epsilon()))) {
- std::cerr << " Problem occured at abscissa " << t << "\n";
- }
- }
- }
- int main()
- {
- test_knots<float>();
- test_knots<double>();
- test_knots<long double>();
- test_bump<double>();
- test_bump<long double>();
- test_trivial<float>();
- test_trivial<double>();
- return boost::math::test::report_errors();
- }
|