]>
Commit | Line | Data |
---|---|---|
1e59de90 TL |
1 | /* |
2 | * Copyright Nick Thompson, 2021 | |
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 <random> | |
11 | #include <array> | |
12 | #include <boost/core/demangle.hpp> | |
13 | #include <boost/math/interpolators/bilinear_uniform.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::bilinear_uniform; | |
20 | ||
21 | template<class Real> | |
22 | void test_four_values() | |
23 | { | |
24 | Real x0 = 0; | |
25 | Real y0 = 0; | |
26 | Real dx = 1; | |
27 | Real dy = 1; | |
28 | Real value = 1.5; | |
29 | std::vector<Real> v(2*2, value); | |
30 | auto v_copy = v; | |
31 | auto ub = bilinear_uniform<decltype(v)>(std::move(v_copy), 2, 2, dx, dy, x0, y0); | |
32 | for (Real x = x0; x <= x0 + dx; x += dx/8) { | |
33 | for (Real y = y0; y <= y0 + dx; y += dy/8) { | |
34 | CHECK_ULP_CLOSE(value, ub(x, y), 1); | |
35 | } | |
36 | } | |
37 | ||
38 | // Now we test the unit square: | |
39 | std::random_device rd; | |
40 | std::uniform_real_distribution<Real> dis(1,2); | |
41 | ||
42 | int i = 0; | |
43 | while (i++ < 300) { | |
44 | v[0] = dis(rd); | |
45 | v[1] = dis(rd); | |
46 | v[2] = dis(rd); | |
47 | v[3] = dis(rd); | |
48 | ||
49 | // See https://en.wikipedia.org/wiki/Bilinear_interpolation, section: Unit square | |
50 | auto f = [&v](Real x, Real y) { | |
51 | return v[0]*(1-x)*(1-y) + v[1]*x*(1-y) + v[2]*(1-x)*y + v[3]*x*y; | |
52 | }; | |
53 | ||
54 | v_copy = v; | |
55 | ub = bilinear_uniform<decltype(v_copy)>(std::move(v_copy), 2, 2, dx, dy, x0, y0); | |
56 | for (Real x = x0; x <= x0 + dx; x += dx/16) { | |
57 | for (Real y = y0; y <= y0 + dx; y += dy/16) { | |
58 | CHECK_ULP_CLOSE(f(x,y), ub(x, y), 3); | |
59 | } | |
60 | } | |
61 | } | |
62 | } | |
63 | ||
64 | template<typename Real> | |
65 | void test_linear() | |
66 | { | |
67 | std::random_device rd; | |
68 | std::uniform_real_distribution<Real> dis(1,2); | |
69 | std::array<Real, 4> a{dis(rd), dis(rd), dis(rd), dis(rd)}; | |
70 | auto f = [&a](Real x, Real y) { | |
71 | return a[0] + a[1]*x + a[2]*y + a[3]*x*y; | |
72 | }; | |
73 | ||
74 | for (int rows = 2; rows < 20; ++rows) { | |
75 | for (int cols = 2; cols < 20; ++cols) { | |
76 | Real dx = dis(rd); | |
77 | Real dy = dis(rd); | |
78 | Real x0 = dis(rd); | |
79 | Real y0 = dis(rd); | |
80 | std::vector<Real> v(rows*cols, std::numeric_limits<Real>::quiet_NaN()); | |
81 | for (int i = 0; i < cols; ++i) { | |
82 | for (int j = 0; j < rows; ++j) { | |
83 | v[j*cols + i] = f(x0 + i*dx, y0 + j*dy); | |
84 | } | |
85 | } | |
86 | auto ub = bilinear_uniform<decltype(v)>(std::move(v), rows, cols, dx, dy, x0, y0); | |
87 | ||
88 | for (Real x = x0; x < x0 + (cols-1)*dx; x += dx/8) { | |
89 | for (Real y = y0; y < y0 + (rows-1)*dy; y += dy/8) { | |
90 | if (!CHECK_ULP_CLOSE(f(x,y), ub(x, y), 13)) { | |
91 | std::cerr << " f(" << x << ", " << y << ") = " << f(x,y) << "\n"; | |
92 | std::cerr << "ub(" << x << ", " << y << ") = " << ub(x,y) << "\n"; | |
93 | } | |
94 | } | |
95 | } | |
96 | } | |
97 | } | |
98 | } | |
99 | ||
100 | ||
101 | ||
102 | int main() | |
103 | { | |
104 | test_four_values<float>(); | |
105 | test_four_values<double>(); | |
106 | test_four_values<long double>(); | |
107 | test_linear<double>(); | |
108 | return boost::math::test::report_errors(); | |
109 | } |