143 lines
3.9 KiB
C++
143 lines
3.9 KiB
C++
/*
|
|
* 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();
|
|
}
|