1 /*
2 * Copyright Nick Thompson, 2019
3 * Use, modification and distribution are subject to the
4 * Boost Software License, Version 1.0. (See accompanying file
5 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 */
7
8 #include "math_unit_test.hpp"
9 #include <numeric>
10 #include <utility>
11 #include <random>
12 #include <boost/core/demangle.hpp>
13 #include <boost/math/interpolators/whittaker_shannon.hpp>
14 #ifdef BOOST_HAS_FLOAT128
15 #include <boost/multiprecision/float128.hpp>
16 using boost::multiprecision::float128;
17 #endif
18
19 using boost::math::interpolators::whittaker_shannon;
20
21 template<class Real>
test_trivial()22 void test_trivial()
23 {
24 Real t0 = 0;
25 Real h = Real(1)/Real(16);
26 std::vector<Real> v{1.5};
27 std::vector<Real> v_copy = v;
28 auto ws = whittaker_shannon<decltype(v)>(std::move(v), t0, h);
29
30
31 Real expected = 0;
32 if(!CHECK_MOLLIFIED_CLOSE(expected, ws.prime(0), 10*std::numeric_limits<Real>::epsilon())) {
33 std::cerr << " Problem occurred at abscissa " << 0 << "\n";
34 }
35
36 expected = -v_copy[0]/h;
37 if(!CHECK_MOLLIFIED_CLOSE(expected, ws.prime(h), 10*std::numeric_limits<Real>::epsilon())) {
38 std::cerr << " Problem occurred at abscissa " << 0 << "\n";
39 }
40 }
41
42 template<class Real>
test_knots()43 void test_knots()
44 {
45 Real t0 = 0;
46 Real h = Real(1)/Real(16);
47 size_t n = 512;
48 std::vector<Real> v(n);
49 std::mt19937 gen(323723);
50 std::uniform_real_distribution<Real> dis(1.0, 2.0);
51
52 for(size_t i = 0; i < n; ++i) {
53 v[i] = static_cast<Real>(dis(gen));
54 }
55 auto ws = whittaker_shannon<decltype(v)>(std::move(v), t0, h);
56
57 size_t i = 0;
58 while (i < n) {
59 Real t = t0 + i*h;
60 Real expected = ws[i];
61 Real computed = ws(t);
62 CHECK_ULP_CLOSE(expected, computed, 16);
63 ++i;
64 }
65 }
66
67 template<class Real>
test_bump()68 void test_bump()
69 {
70 using std::exp;
71 using std::abs;
72 using std::sqrt;
73 auto bump = [](Real x) { if (abs(x) >= 1) { return Real(0); } return exp(-Real(1)/(Real(1)-x*x)); };
74
75 auto bump_prime = [&bump](Real x) { Real z = 1-x*x; return -2*x*bump(x)/(z*z); };
76
77 Real t0 = -1;
78 size_t n = 2049;
79 Real h = Real(2)/Real(n-1);
80
81 std::vector<Real> v(n);
82 for(size_t i = 0; i < n; ++i) {
83 Real t = t0 + i*h;
84 v[i] = bump(t);
85 }
86
87
88 std::vector<Real> v_copy = v;
89 auto ws = whittaker_shannon<decltype(v)>(std::move(v), t0, h);
90
91 // Test the knots:
92 for(size_t i = v_copy.size()/4; i < 3*v_copy.size()/4; ++i) {
93 Real t = t0 + i*h;
94 Real expected = v_copy[i];
95 Real computed = ws(t);
96 if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 10*std::numeric_limits<Real>::epsilon())) {
97 std::cerr << " Problem occurred at abscissa " << t << "\n";
98 }
99
100 Real expected_prime = bump_prime(t);
101 Real computed_prime = ws.prime(t);
102 if(!CHECK_MOLLIFIED_CLOSE(expected_prime, computed_prime, 1000*std::numeric_limits<Real>::epsilon())) {
103 std::cerr << " Problem occurred at abscissa " << t << "\n";
104 }
105
106 }
107
108 std::mt19937 gen(323723);
109 std::uniform_real_distribution<long double> dis(-0.85, 0.85);
110
111 size_t i = 0;
112 while (i++ < 1000)
113 {
114 Real t = static_cast<Real>(dis(gen));
115 Real expected = bump(t);
116 Real computed = ws(t);
117 if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 10*std::numeric_limits<Real>::epsilon())) {
118 std::cerr << " Problem occurred at abscissa " << t << "\n";
119 }
120
121 Real expected_prime = bump_prime(t);
122 Real computed_prime = ws.prime(t);
123 if(!CHECK_MOLLIFIED_CLOSE(expected_prime, computed_prime, sqrt(std::numeric_limits<Real>::epsilon()))) {
124 std::cerr << " Problem occurred at abscissa " << t << "\n";
125 }
126 }
127 }
128
129
main()130 int main()
131 {
132 test_knots<float>();
133 test_knots<double>();
134 test_knots<long double>();
135
136 test_bump<double>();
137 test_bump<long double>();
138
139 test_trivial<float>();
140 test_trivial<double>();
141 return boost::math::test::report_errors();
142 }
143