]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
1 | /* |
2 | * (C) Copyright Nick Thompson 2018. | |
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 <vector> | |
9 | #include <array> | |
10 | #include <forward_list> | |
11 | #include <algorithm> | |
12 | #include <random> | |
13 | #include <boost/core/lightweight_test.hpp> | |
14 | #include <boost/numeric/ublas/vector.hpp> | |
15 | #include <boost/math/constants/constants.hpp> | |
16 | #include <boost/math/statistics/univariate_statistics.hpp> | |
17 | #include <boost/math/statistics/signal_statistics.hpp> | |
18 | #include <boost/multiprecision/cpp_bin_float.hpp> | |
19 | #include <boost/multiprecision/cpp_complex.hpp> | |
20 | ||
21 | using std::abs; | |
22 | using boost::multiprecision::cpp_bin_float_50; | |
23 | using boost::multiprecision::cpp_complex_50; | |
24 | using boost::math::constants::two_pi; | |
25 | ||
26 | /* | |
27 | * Test checklist: | |
28 | * 1) Does it work with multiprecision? | |
29 | * 2) Does it work with .cbegin()/.cend() if the data is not altered? | |
30 | * 3) Does it work with ublas and std::array? (Checking Eigen and Armadillo will make the CI system really unhappy.) | |
31 | * 4) Does it work with std::forward_list if a forward iterator is all that is required? | |
32 | * 5) Does it work with complex data if complex data is sensible? | |
33 | * 6) Does it work with integer data if sensible? | |
34 | */ | |
35 | ||
36 | template<class Real> | |
37 | void test_hoyer_sparsity() | |
38 | { | |
39 | using std::sqrt; | |
40 | Real tol = 5*std::numeric_limits<Real>::epsilon(); | |
41 | std::vector<Real> v{1,0,0}; | |
42 | Real hs = boost::math::statistics::hoyer_sparsity(v.begin(), v.end()); | |
43 | BOOST_TEST(abs(hs - 1) < tol); | |
44 | ||
45 | hs = boost::math::statistics::hoyer_sparsity(v); | |
46 | BOOST_TEST(abs(hs - 1) < tol); | |
47 | ||
48 | // Does it work with constant iterators? | |
49 | hs = boost::math::statistics::hoyer_sparsity(v.cbegin(), v.cend()); | |
50 | BOOST_TEST(abs(hs - 1) < tol); | |
51 | ||
52 | v[0] = 1; | |
53 | v[1] = 1; | |
54 | v[2] = 1; | |
55 | hs = boost::math::statistics::hoyer_sparsity(v.cbegin(), v.cend()); | |
56 | BOOST_TEST(abs(hs) < tol); | |
57 | ||
58 | std::array<Real, 3> w{1,1,1}; | |
59 | hs = boost::math::statistics::hoyer_sparsity(w); | |
60 | BOOST_TEST(abs(hs) < tol); | |
61 | ||
62 | // Now some statistics: | |
63 | // If x_i ~ Unif(0,1), E[x_i] = 1/2, E[x_i^2] = 1/3. | |
64 | // Therefore, E[||x||_1] = N/2, E[||x||_2] = sqrt(N/3), | |
65 | // and hoyer_sparsity(x) is close to (1-sqrt(3)/2)/(1-1/sqrt(N)) | |
66 | std::mt19937 gen(82); | |
67 | std::uniform_real_distribution<long double> dis(0, 1); | |
68 | v.resize(5000); | |
69 | for (size_t i = 0; i < v.size(); ++i) { | |
70 | v[i] = dis(gen); | |
71 | } | |
72 | hs = boost::math::statistics::hoyer_sparsity(v); | |
73 | Real expected = (1.0 - boost::math::constants::root_three<Real>()/2)/(1.0 - 1.0/sqrt(v.size())); | |
74 | BOOST_TEST(abs(expected - hs) < 0.01); | |
75 | ||
76 | // Does it work with a forward list? | |
77 | std::forward_list<Real> u1{1, 1, 1}; | |
78 | hs = boost::math::statistics::hoyer_sparsity(u1); | |
79 | BOOST_TEST(abs(hs) < tol); | |
80 | ||
81 | // Does it work with a boost ublas vector? | |
82 | boost::numeric::ublas::vector<Real> u2(3); | |
83 | u2[0] = 1; | |
84 | u2[1] = 1; | |
85 | u2[2] = 1; | |
86 | hs = boost::math::statistics::hoyer_sparsity(u2); | |
87 | BOOST_TEST(abs(hs) < tol); | |
88 | ||
89 | } | |
90 | ||
91 | template<class Z> | |
92 | void test_integer_hoyer_sparsity() | |
93 | { | |
94 | using std::sqrt; | |
95 | double tol = 5*std::numeric_limits<double>::epsilon(); | |
96 | std::vector<Z> v{1,0,0}; | |
97 | double hs = boost::math::statistics::hoyer_sparsity(v); | |
98 | BOOST_TEST(abs(hs - 1) < tol); | |
99 | ||
100 | v[0] = 1; | |
101 | v[1] = 1; | |
102 | v[2] = 1; | |
103 | hs = boost::math::statistics::hoyer_sparsity(v); | |
104 | BOOST_TEST(abs(hs) < tol); | |
105 | } | |
106 | ||
107 | ||
108 | template<class Complex> | |
109 | void test_complex_hoyer_sparsity() | |
110 | { | |
111 | typedef typename Complex::value_type Real; | |
112 | using std::sqrt; | |
113 | Real tol = 5*std::numeric_limits<Real>::epsilon(); | |
114 | std::vector<Complex> v{{0,1}, {0, 0}, {0,0}}; | |
115 | Real hs = boost::math::statistics::hoyer_sparsity(v.begin(), v.end()); | |
116 | BOOST_TEST(abs(hs - 1) < tol); | |
117 | ||
118 | hs = boost::math::statistics::hoyer_sparsity(v); | |
119 | BOOST_TEST(abs(hs - 1) < tol); | |
120 | ||
121 | // Does it work with constant iterators? | |
122 | hs = boost::math::statistics::hoyer_sparsity(v.cbegin(), v.cend()); | |
123 | BOOST_TEST(abs(hs - 1) < tol); | |
124 | ||
125 | // All are the same magnitude: | |
126 | v[0] = {0, 1}; | |
127 | v[1] = {1, 0}; | |
128 | v[2] = {0,-1}; | |
129 | hs = boost::math::statistics::hoyer_sparsity(v.cbegin(), v.cend()); | |
130 | BOOST_TEST(abs(hs) < tol); | |
131 | } | |
132 | ||
133 | ||
134 | template<class Real> | |
135 | void test_absolute_gini_coefficient() | |
136 | { | |
137 | using boost::math::statistics::absolute_gini_coefficient; | |
138 | using boost::math::statistics::sample_absolute_gini_coefficient; | |
139 | Real tol = std::numeric_limits<Real>::epsilon(); | |
140 | std::vector<Real> v{-1,0,0}; | |
141 | Real gini = sample_absolute_gini_coefficient(v.begin(), v.end()); | |
142 | BOOST_TEST(abs(gini - 1) < tol); | |
143 | ||
144 | gini = absolute_gini_coefficient(v); | |
145 | BOOST_TEST(abs(gini - Real(2)/Real(3)) < tol); | |
146 | ||
147 | v[0] = 1; | |
148 | v[1] = -1; | |
149 | v[2] = 1; | |
150 | gini = absolute_gini_coefficient(v.begin(), v.end()); | |
151 | BOOST_TEST(abs(gini) < tol); | |
152 | gini = sample_absolute_gini_coefficient(v.begin(), v.end()); | |
153 | BOOST_TEST(abs(gini) < tol); | |
154 | ||
155 | std::vector<std::complex<Real>> w(128); | |
156 | std::complex<Real> i{0,1}; | |
157 | for(size_t k = 0; k < w.size(); ++k) | |
158 | { | |
159 | w[k] = exp(i*static_cast<Real>(k)/static_cast<Real>(w.size())); | |
160 | } | |
161 | gini = absolute_gini_coefficient(w.begin(), w.end()); | |
162 | BOOST_TEST(abs(gini) < tol); | |
163 | gini = sample_absolute_gini_coefficient(w.begin(), w.end()); | |
164 | BOOST_TEST(abs(gini) < tol); | |
165 | ||
166 | // The population Gini index is invariant under "cloning": If w = v \oplus v, then G(w) = G(v). | |
167 | // We use the sample Gini index, so we need to rescale | |
168 | std::vector<Real> u(1000); | |
169 | std::mt19937 gen(35); | |
170 | std::uniform_real_distribution<long double> dis(0, 50); | |
171 | for (size_t i = 0; i < u.size()/2; ++i) | |
172 | { | |
173 | u[i] = dis(gen); | |
174 | } | |
175 | for (size_t i = 0; i < u.size()/2; ++i) | |
176 | { | |
177 | u[i + u.size()/2] = u[i]; | |
178 | } | |
179 | Real population_gini1 = absolute_gini_coefficient(u.begin(), u.begin() + u.size()/2); | |
180 | Real population_gini2 = absolute_gini_coefficient(u.begin(), u.end()); | |
181 | ||
182 | BOOST_TEST(abs(population_gini1 - population_gini2) < 10*tol); | |
183 | ||
184 | // The Gini coefficient of a uniform distribution is (b-a)/(3*(b+a)), see https://en.wikipedia.org/wiki/Gini_coefficient | |
185 | Real expected = (dis.b() - dis.a() )/(3*(dis.a() + dis.b())); | |
186 | ||
187 | BOOST_TEST(abs(expected - population_gini1) < 0.01); | |
188 | ||
189 | std::exponential_distribution<long double> exp_dis(1); | |
190 | for (size_t i = 0; i < u.size(); ++i) | |
191 | { | |
192 | u[i] = exp_dis(gen); | |
193 | } | |
194 | population_gini2 = absolute_gini_coefficient(u); | |
195 | ||
196 | BOOST_TEST(abs(population_gini2 - 0.5) < 0.01); | |
197 | } | |
198 | ||
199 | ||
200 | template<class Real> | |
201 | void test_oracle_snr() | |
202 | { | |
203 | using std::abs; | |
204 | Real tol = 100*std::numeric_limits<Real>::epsilon(); | |
205 | size_t length = 100; | |
206 | std::vector<Real> signal(length, 1); | |
207 | std::vector<Real> noisy_signal = signal; | |
208 | ||
209 | noisy_signal[0] += 1; | |
210 | Real snr = boost::math::statistics::oracle_snr(signal, noisy_signal); | |
211 | Real snr_db = boost::math::statistics::oracle_snr_db(signal, noisy_signal); | |
212 | BOOST_TEST(abs(snr - length) < tol); | |
213 | BOOST_TEST(abs(snr_db - 10*log10(length)) < tol); | |
214 | } | |
215 | ||
216 | template<class Z> | |
217 | void test_integer_oracle_snr() | |
218 | { | |
219 | double tol = std::numeric_limits<double>::epsilon(); | |
220 | size_t length = 100; | |
221 | std::vector<Z> signal(length, 1); | |
222 | std::vector<Z> noisy_signal = signal; | |
223 | ||
224 | noisy_signal[0] += 1; | |
225 | double snr = boost::math::statistics::oracle_snr(signal, noisy_signal); | |
226 | double snr_db = boost::math::statistics::oracle_snr_db(signal, noisy_signal); | |
227 | BOOST_TEST(abs(snr - length) < tol); | |
228 | BOOST_TEST(abs(snr_db - 10*log10(length)) < tol); | |
229 | } | |
230 | ||
231 | template<class Complex> | |
232 | void test_complex_oracle_snr() | |
233 | { | |
234 | using Real = typename Complex::value_type; | |
235 | using std::abs; | |
236 | using std::log10; | |
237 | Real tol = 100*std::numeric_limits<Real>::epsilon(); | |
238 | size_t length = 100; | |
239 | std::vector<Complex> signal(length, {1,0}); | |
240 | std::vector<Complex> noisy_signal = signal; | |
241 | ||
242 | noisy_signal[0] += Complex(1,0); | |
243 | Real snr = boost::math::statistics::oracle_snr(signal, noisy_signal); | |
244 | Real snr_db = boost::math::statistics::oracle_snr_db(signal, noisy_signal); | |
245 | BOOST_TEST(abs(snr - length) < tol); | |
246 | BOOST_TEST(abs(snr_db - 10*log10(length)) < tol); | |
247 | } | |
248 | ||
249 | template<class Real> | |
250 | void test_m2m4_snr_estimator() | |
251 | { | |
252 | Real tol = std::numeric_limits<Real>::epsilon(); | |
253 | std::vector<Real> signal(5000, 1); | |
254 | std::vector<Real> x(signal.size()); | |
255 | std::mt19937 gen(18); | |
256 | std::normal_distribution<Real> dis{0, 1.0}; | |
257 | ||
258 | for (size_t i = 0; i < x.size(); ++i) | |
259 | { | |
260 | signal[i] = 5*sin(100*6.28*i/x.size()); | |
261 | x[i] = signal[i] + dis(gen); | |
262 | } | |
263 | ||
264 | // Kurtosis of a sine wave is 1.5: | |
265 | auto m2m4_db = boost::math::statistics::m2m4_snr_estimator_db(x, 1.5); | |
266 | auto oracle_snr_db = boost::math::statistics::mean_invariant_oracle_snr_db(signal, x); | |
267 | BOOST_TEST(abs(m2m4_db - oracle_snr_db) < 0.2); | |
268 | ||
269 | std::uniform_real_distribution<Real> uni_dis{-1,1}; | |
270 | for (size_t i = 0; i < x.size(); ++i) | |
271 | { | |
272 | x[i] = signal[i] + uni_dis(gen); | |
273 | } | |
274 | ||
275 | // Kurtosis of continuous uniform distribution over [-1,1] is 1.8: | |
276 | m2m4_db = boost::math::statistics::m2m4_snr_estimator_db(x, 1.5, 1.8); | |
277 | oracle_snr_db = boost::math::statistics::mean_invariant_oracle_snr_db(signal, x); | |
278 | // The performance depends on the exact numbers generated by the distribution, but this isn't bad: | |
279 | BOOST_TEST(abs(m2m4_db - oracle_snr_db) < 0.2); | |
280 | ||
281 | // The SNR estimator should be scale invariant. | |
282 | // If x has snr y, then kx should have snr y. | |
283 | Real ka = 1.5; | |
284 | Real kw = 1.8; | |
285 | auto m2m4 = boost::math::statistics::m2m4_snr_estimator(x.begin(), x.end(), ka, kw); | |
286 | for(size_t i = 0; i < x.size(); ++i) | |
287 | { | |
288 | x[i] *= 4096; | |
289 | } | |
290 | auto m2m4_2 = boost::math::statistics::m2m4_snr_estimator(x.begin(), x.end(), ka, kw); | |
291 | BOOST_TEST(abs(m2m4 - m2m4_2) < tol); | |
292 | } | |
293 | ||
294 | int main() | |
295 | { | |
296 | test_absolute_gini_coefficient<float>(); | |
297 | test_absolute_gini_coefficient<double>(); | |
298 | test_absolute_gini_coefficient<long double>(); | |
299 | ||
300 | test_hoyer_sparsity<float>(); | |
301 | test_hoyer_sparsity<double>(); | |
302 | test_hoyer_sparsity<long double>(); | |
303 | test_hoyer_sparsity<cpp_bin_float_50>(); | |
304 | ||
305 | test_integer_hoyer_sparsity<int>(); | |
306 | test_integer_hoyer_sparsity<unsigned>(); | |
307 | ||
308 | test_complex_hoyer_sparsity<std::complex<float>>(); | |
309 | test_complex_hoyer_sparsity<std::complex<double>>(); | |
310 | test_complex_hoyer_sparsity<std::complex<long double>>(); | |
311 | test_complex_hoyer_sparsity<cpp_complex_50>(); | |
312 | ||
313 | test_oracle_snr<float>(); | |
314 | test_oracle_snr<double>(); | |
315 | test_oracle_snr<long double>(); | |
316 | test_oracle_snr<cpp_bin_float_50>(); | |
317 | ||
318 | test_integer_oracle_snr<int>(); | |
319 | test_integer_oracle_snr<unsigned>(); | |
320 | ||
321 | test_complex_oracle_snr<std::complex<float>>(); | |
322 | test_complex_oracle_snr<std::complex<double>>(); | |
323 | test_complex_oracle_snr<std::complex<long double>>(); | |
324 | test_complex_oracle_snr<cpp_complex_50>(); | |
325 | ||
326 | test_m2m4_snr_estimator<float>(); | |
327 | test_m2m4_snr_estimator<double>(); | |
328 | test_m2m4_snr_estimator<long double>(); | |
329 | ||
330 | return boost::report_errors(); | |
331 | } |