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