1 /***************************************************************************
2 * Copyright (c) Johan Mabille, Sylvain Corlay and Wolf Vollprecht          *
3 * Copyright (c) QuantStack                                                 *
4 *                                                                          *
5 * Distributed under the terms of the BSD 3-Clause License.                 *
6 *                                                                          *
7 * The full license is in the file LICENSE, distributed with this software. *
8 ****************************************************************************/
9 
10 #include "test_common_macros.hpp"
11 
12 #if (defined(__GNUC__) && !defined(__clang__))
13 #pragma GCC diagnostic push
14 #pragma GCC diagnostic ignored "-Wconversion"
15 #include "xtensor/xmath.hpp"
16 #pragma GCC diagnostic pop
17 #endif
18 
19 #include "xtensor/xarray.hpp"
20 #include "xtensor/xtensor.hpp"
21 #include "xtensor/xmath.hpp"
22 #include "xtensor/xstrided_view.hpp"
23 #include <xtensor/xindex_view.hpp>
24 
25 #include "xtl/xtype_traits.hpp"
26 
27 namespace xt
28 {
29     using namespace std::complex_literals;
30     static const double nanv = std::nan("0");
31     static const double d_min = std::numeric_limits<double>::min();
32     static const double d_max = std::numeric_limits<double>::max();
33 
34 #define NAN_SENSITIVE_EQ(E1, E2, PLACE_HOLDER)                                                          \
35         EXPECT_EQ(xt::isnan(E1), xt::equal(E2, PLACE_HOLDER));                                          \
36         EXPECT_EQ(xt::filter(E1, !xt::isnan(E1)), xt::filter(E2, xt::not_equal(E2, PLACE_HOLDER)));
37 
38     namespace nantest
39     {
40         xarray<double> aN = {{ nanv, nanv, 123, 3 }, { 1, 2, nanv, 3 }, { 1, 1, nanv, 3 }};
41         xarray<double> aR = {{ 0, 0, 123 , 3 }, { 1, 2, 0 , 3}, { 1, 1, 0, 3 }};
42         xarray<double> aP = {{ 1, 1, 123 , 3}, { 1, 2, 1, 3 }, { 1, 1, 1, 3 }};
43         xarray<double> aI = {{ d_max, d_max, 123 , 3}, { 1, 2, d_max, 3 }, { 1, 1, d_max, 3 }};
44         xarray<double> aA = {{ d_min, d_min, 123 , 3}, { 1, 2, d_min, 3 }, { 1, 1, d_min, 3 }};
45 
46         xarray<double> xN = {{{nanv, nanv}, {1, 2}}, {{3, nanv}, {nanv, 5}}};
47         xarray<double> xR = {{{0, 0}, {1,2}}, {{3, 0}, {0, 5}}};
48         xarray<double> xP = {{{1, 1}, {1,2}}, {{3, 1}, {1, 5}}};
49         xarray<double> xI = {{{d_max, d_max}, {1, 2}}, {{3, d_max}, {d_max, 5}}};
50         xarray<double> xA = {{{d_min, d_min}, {1, 2}}, {{3, d_min}, {d_min, 5}}};
51 
52         xarray<std::complex<double>> cN = {
53             {1.0 + 1.0i, 1.0 + 1.0i, nanv      },
54             {1.0 - 1.0i, 1.0       , 3.0 + 2.0i}
55         };
56 
57     }
58 
TEST(xnanfunctions,count_nonnan)59     TEST(xnanfunctions, count_nonnan)
60     {
61         xarray<double> a = {{0, 1, 2, 3}, {nanv, nanv, nanv, nanv}, {3, nanv, 1, nanv}};
62         std::size_t as = count_nonnan(a)();
63         std::size_t ase = count_nonnan(a, evaluation_strategy::immediate)();
64         EXPECT_EQ(as, 6u);
65         EXPECT_EQ(ase, 6u);
66 
67         xarray<std::size_t> ea0 = {2, 1, 2, 1};
68         xarray<std::size_t> ea1 = {4, 0, 2};
69 
70         EXPECT_EQ(count_nonnan(a, {0}), ea0);
71         EXPECT_EQ(count_nonnan(a, {1}), ea1);
72 
73         EXPECT_EQ(count_nonnan(a, {0}, evaluation_strategy::immediate), ea0);
74         EXPECT_EQ(count_nonnan(a, {1}, evaluation_strategy::immediate), ea1);
75     }
76 
TEST(xnanfunctions,nan_to_num)77     TEST(xnanfunctions, nan_to_num)
78     {
79         double neg_inf = -std::numeric_limits<double>::infinity();
80         double inf = std::numeric_limits<double>::infinity();
81         xarray<double> a = {{ nanv, nanv, 123}, {0.5123, neg_inf, inf}};
82 
83         auto expr = nan_to_num(a);
84         EXPECT_EQ(expr(0, 0), 0);
85         EXPECT_EQ(expr(0, 1), 0);
86         EXPECT_EQ(expr(0, 2), 123);
87         EXPECT_EQ(expr(1, 0), 0.5123);
88         EXPECT_EQ(expr(1, 1), std::numeric_limits<double>::lowest());
89         EXPECT_TRUE(expr(1, 1) < 0);
90         EXPECT_EQ(expr(1, 2), std::numeric_limits<double>::max());
91 
92         xarray<double> exp = {{ 0, 0, 123 }, {0.5123, std::numeric_limits<double>::lowest(), std::numeric_limits<double>::max() }};
93         xarray<double> assigned = exp;
94         EXPECT_EQ(assigned, exp);
95     }
96 
TEST(xnanfunctions,nanmin)97     TEST(xnanfunctions, nanmin)
98     {
99         EXPECT_EQ(nanmin(nantest::aN), amin(nantest::aI));
100         EXPECT_EQ(nanmin(nantest::aN, {0}), amin(nantest::aI, {0}));
101         EXPECT_EQ(nanmin(nantest::aN, {1}), amin(nantest::aI, {1}));
102 
103         for (size_t i = 0; i < 3; ++i)
104         {
105             auto ret = nanmin(nantest::xN, {i});
106             auto reference = amin(nantest::xI, {i});
107             NAN_SENSITIVE_EQ(ret, reference, d_max)
108         }
109     }
110 
TEST(xnanfunctions,nanmax)111     TEST(xnanfunctions, nanmax)
112     {
113         EXPECT_EQ(nanmax(nantest::aN), amax(nantest::aA));
114         EXPECT_EQ(nanmax(nantest::aN, {0}), amax(nantest::aA, {0}));
115         EXPECT_EQ(nanmax(nantest::aN, {1}), amax(nantest::aA, {1}));
116 
117         for (size_t i = 0; i < 3; ++i)
118         {
119             auto ret = nanmax(nantest::xN, {i});
120             auto reference = amax(nantest::xA, {i});
121             NAN_SENSITIVE_EQ(ret, reference, d_min)
122         }
123     }
124 
TEST(xnanfunctions,nansum)125     TEST(xnanfunctions, nansum)
126     {
127         xarray<double> res = nansum(nantest::aN);
128         xarray<double> res0 = sum(nantest::aR);
129         EXPECT_EQ(res(0), 137);
130 
131         EXPECT_EQ(nansum(nantest::aN, {0}), sum(nantest::aR, {0}));
132         EXPECT_EQ(nansum(nantest::aN, {1}), sum(nantest::aR, {1}));
133         EXPECT_EQ(nansum(nantest::xN, {0}), sum(nantest::xR, {0}));
134         EXPECT_EQ(nansum(nantest::xN, {1}), sum(nantest::xR, {1}));
135         EXPECT_EQ(nansum(nantest::xN, {2}), sum(nantest::xR, {2}));
136     }
137 
TEST(xnanfunctions,nanprod)138     TEST(xnanfunctions, nanprod)
139     {
140         xarray<double> res = nanprod(nantest::aN);
141         xarray<double> res0 = prod(nantest::aP);
142         EXPECT_EQ(res(0), 6642);
143 
144         EXPECT_EQ(nanprod(nantest::aN, {0}), prod(nantest::aP, {0}));
145         EXPECT_EQ(nanprod(nantest::aN, {1}), prod(nantest::aP, {1}));
146         EXPECT_EQ(nanprod(nantest::xN, {0}), prod(nantest::xP, {0}));
147         EXPECT_EQ(nanprod(nantest::xN, {1}), prod(nantest::xP, {1}));
148         EXPECT_EQ(nanprod(nantest::xN, {2}), prod(nantest::xP, {2}));
149     }
150 
TEST(xnanfunctions,nancumsum)151     TEST(xnanfunctions, nancumsum)
152     {
153         EXPECT_EQ(nancumsum(nantest::aN), cumsum(nantest::aR));
154         EXPECT_EQ(nancumsum(nantest::aN, 0), cumsum(nantest::aR, 0));
155         EXPECT_EQ(nancumsum(nantest::aN, 1), cumsum(nantest::aR, 1));
156         EXPECT_EQ(nancumsum(nantest::xN, 0), cumsum(nantest::xR, 0));
157         EXPECT_EQ(nancumsum(nantest::xN, 1), cumsum(nantest::xR, 1));
158         EXPECT_EQ(nancumsum(nantest::xN, 2), cumsum(nantest::xR, 2));
159     }
160 
TEST(xnanfunctions,multid)161     TEST(xnanfunctions, multid)
162     {
163         xarray<double> arr = xt::arange(3 * 4 * 2 * 8 * 7);
164         arr.reshape({3, 4, 2, 8, 7});
165         xarray<double> carr = arr;
166         strided_view(arr, {0, xt::ellipsis()}) = nanv;
167         strided_view(carr, {0, xt::ellipsis()}) = 0;
168 
169         EXPECT_EQ(nancumsum(arr, 0), cumsum(carr, 0));
170         EXPECT_EQ(nancumsum(arr, 1), cumsum(carr, 1));
171         EXPECT_EQ(nancumsum(arr, 2), cumsum(carr, 2));
172         EXPECT_EQ(nancumsum(arr, 3), cumsum(carr, 3));
173         EXPECT_EQ(nancumsum(arr, 4), cumsum(carr, 4));
174     }
175 
TEST(xnanfunctions,nancumprod)176     TEST(xnanfunctions, nancumprod)
177     {
178         EXPECT_EQ(nancumprod(nantest::aN), cumprod(nantest::aP));
179         EXPECT_EQ(nancumprod(nantest::aN, 0), cumprod(nantest::aP, 0));
180         EXPECT_EQ(nancumprod(nantest::aN, 1), cumprod(nantest::aP, 1));
181         EXPECT_EQ(nancumprod(nantest::xN, 0), cumprod(nantest::xP, 0));
182         EXPECT_EQ(nancumprod(nantest::xN, 1), cumprod(nantest::xP, 1));
183         EXPECT_EQ(nancumprod(nantest::xN, 2), cumprod(nantest::xP, 2));
184     }
185 
TEST(xnanfunctions,nanmean)186     TEST(xnanfunctions, nanmean)
187     {
188         auto as = nanmean(nantest::aN)();
189         auto ase = nanmean(nantest::aN, evaluation_strategy::immediate)();
190         EXPECT_DOUBLE_EQ(as, 17.125);
191         EXPECT_DOUBLE_EQ(ase, 17.125);
192 
193         xarray<double> eaN0 = {1.0, 1.5, 123, 3};
194         xarray<double> eaN1 = {63.0, 2.0, 5.0/3.0};
195 
196         EXPECT_TENSOR_EQ(nanmean(nantest::aN, {0}), eaN0);
197         EXPECT_TENSOR_EQ(nanmean(nantest::aN, {1}), eaN1);
198 
199         std::array<std::size_t, 1> axis{0};
200         EXPECT_EQ(nanmean(nantest::aN, axis), eaN0);
201 
202         EXPECT_TENSOR_EQ(nanmean(nantest::aN, {0}, evaluation_strategy::immediate), eaN0);
203         EXPECT_TENSOR_EQ(nanmean(nantest::aN, {1}, evaluation_strategy::immediate), eaN1);
204 
205         auto cs = nanmean(nantest::cN)();
206         auto cse = nanmean(nantest::cN, evaluation_strategy::immediate)();
207         EXPECT_DOUBLE_EQ(cs, std::complex<double>(1.4, 0.6));
208         EXPECT_DOUBLE_EQ(cse, std::complex<double>(1.4, 0.6));
209 
210         xarray<std::complex<double>> ecN0 = {1.0 + 0.0i, 1.0+0.5i, 3.0+2.0i};
211         xarray<std::complex<double>> ecN1 = {1.0 + 1.0i, (5.0 + 1.0i) / 3.0};
212 
213         EXPECT_TENSOR_EQ(nanmean(nantest::cN, {0}), ecN0);
214         EXPECT_TENSOR_EQ(nanmean(nantest::cN, {1}), ecN1);
215 
216         EXPECT_TENSOR_EQ(nanmean(nantest::cN, {0}, evaluation_strategy::immediate), ecN0);
217         EXPECT_TENSOR_EQ(nanmean(nantest::cN, {1}, evaluation_strategy::immediate), ecN1);
218     }
219 
220 
TEST(xnanfunctions,nanvar)221     TEST(xnanfunctions, nanvar)
222     {
223         auto as = nanvar(nantest::aN)();
224         auto ase = nanvar(nantest::aN, evaluation_strategy::immediate)();
225         EXPECT_EQ(as, 1602.109375);
226         EXPECT_EQ(ase, 1602.109375);
227 
228         xarray<double> eaN0 = {0.0, 0.25, 0.0, 0.0};
229         xarray<double> eaN1 = {3600.0, 2.0/3.0, 8.0/9.0};
230 
231         EXPECT_EQ(nanvar(nantest::aN, {0}), eaN0);
232         EXPECT_TRUE(allclose(nanvar(nantest::aN, {1}), eaN1));
233 
234         std::array<std::size_t, 1> axis{0};
235         EXPECT_EQ(nanvar(nantest::aN, axis), eaN0);
236 
237         EXPECT_EQ(nanvar(nantest::aN, {0}, evaluation_strategy::immediate), eaN0);
238         EXPECT_TRUE(allclose(nanvar(nantest::aN, {1}, evaluation_strategy::immediate), eaN1));
239     }
240 
241     using shape_type = dynamic_shape<size_t>;
242 
243     /*******************
244      * type conversion *
245      *******************/
246 
247 #define CHECK_RESULT_TYPE(EXPRESSION, EXPECTED_TYPE)                                 \
248     {                                                                                \
249         using result_type = typename std::decay_t<decltype(EXPRESSION)>::value_type; \
250         EXPECT_TRUE((std::is_same<result_type, EXPECTED_TYPE>::value));              \
251     }
252 
TEST(xnanfunctions,result_type)253     TEST(xnanfunctions, result_type) {
254         shape_type shape = {4, 3, 2};
255         xarray<short> ashort(shape);
256         xarray<unsigned short> aushort(shape);
257         xarray<int> aint(shape);
258         xarray<unsigned int> auint(shape);
259         xarray<long long> along(shape);
260         xarray<unsigned long long> aulong(shape);
261         xarray<float> afloat(shape);
262         xarray<double> adouble(shape);
263 
264 #define CHECK_RESULT_TYPE_FOR_ALL(INPUT, RESULT_TYPE, MINMAX_TYPE)  \
265         CHECK_RESULT_TYPE(nansum(INPUT, {1, 2}), RESULT_TYPE);      \
266         CHECK_RESULT_TYPE(nanmean(INPUT, {1, 2}), double);          \
267         CHECK_RESULT_TYPE(nanvar(INPUT, {1, 2}), double);           \
268         CHECK_RESULT_TYPE(nanstd(INPUT, {1, 2}), double);           \
269         CHECK_RESULT_TYPE(nanmin(INPUT, {1, 2}), MINMAX_TYPE);      \
270         CHECK_RESULT_TYPE(nanmax(INPUT, {1, 2}), MINMAX_TYPE);
271 
272 #define CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(INPUT, TEMPLATE_TYPE, RESULT_TYPE, STD_TYPE, MINMAX_TYPE)  \
273         CHECK_RESULT_TYPE(nansum<TEMPLATE_TYPE>(INPUT, {1, 2}), RESULT_TYPE)                           \
274         CHECK_RESULT_TYPE(nanmean<TEMPLATE_TYPE>(INPUT, {1, 2}), RESULT_TYPE)                          \
275         CHECK_RESULT_TYPE(nanvar<TEMPLATE_TYPE>(INPUT, {1, 2}), RESULT_TYPE)                           \
276         CHECK_RESULT_TYPE(nanstd<TEMPLATE_TYPE>(INPUT, {1, 2}), STD_TYPE)                              \
277         CHECK_RESULT_TYPE(nanmin<TEMPLATE_TYPE>(INPUT, {1, 2}), MINMAX_TYPE)                           \
278         CHECK_RESULT_TYPE(nanmax<TEMPLATE_TYPE>(INPUT, {1, 2}), MINMAX_TYPE)
279 
280         /*********
281          * short *
282          *********/
283         CHECK_RESULT_TYPE_FOR_ALL(ashort, int, short);
284         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(ashort, short, int, double, short);
285         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(ashort, int, int, double, int);
286         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(ashort, double, double, double, double);
287 
288         /******************
289          * unsigned short *
290          ******************/
291         CHECK_RESULT_TYPE_FOR_ALL(aushort, int, unsigned short);
292         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(aushort, unsigned short, int, double, unsigned short);
293         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(aushort, int, int, double, int);
294         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(aushort, double, double, double, double);
295 
296         /*********
297          * int *
298          *********/
299         CHECK_RESULT_TYPE_FOR_ALL(aint, int, int);
300         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(aint, unsigned short, int, double, int);
301         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(aint, int, int, double, int);
302         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(aint, double, double, double, double);
303 
304         /****************
305          * unsigned int *
306          ****************/
307         CHECK_RESULT_TYPE_FOR_ALL(auint, unsigned int, unsigned int);
308         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(auint, unsigned int, unsigned int, double, unsigned int);
309         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(auint, unsigned int, unsigned int, double, unsigned int);
310         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(auint, double, double, double, double);
311 
312         /**********************
313          * long long *
314          **********************/
315 #ifndef SKIP_ON_WERROR
316         // intermediate computation done in double may imply precision loss
317         CHECK_RESULT_TYPE_FOR_ALL(along, signed long long, signed long long);
318 #endif
319         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(along, int, signed long long, double, signed long long);
320         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(along, long, signed long long, double, signed long long);
321         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(along, signed long long, signed long long, double, signed long long);
322 
323         /**********************
324          * unsigned long long *
325          **********************/
326 #ifndef SKIP_ON_WERROR
327         // intermediate computation done in double may imply precision loss
328         CHECK_RESULT_TYPE_FOR_ALL(aulong, unsigned long long, unsigned long long);
329 #endif
330         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(aulong, unsigned int, unsigned long long, double, unsigned long long);
331         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(aulong, unsigned long, unsigned long long, double, unsigned long long);
332         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(aulong, unsigned long long, unsigned long long, double, unsigned long long);
333 
334         /*********
335          * float *
336          *********/
337         CHECK_RESULT_TYPE_FOR_ALL(afloat, float, float);
338 #ifndef SKIP_ON_WERROR
339         // final conversion to int may imply conversion loss
340         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(afloat, int, float, float, float);
341 #endif
342         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(afloat, float, float, float, float);
343         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(afloat, double, double, double, double);
344 
345         /**********
346          * double *
347          **********/
348         CHECK_RESULT_TYPE_FOR_ALL(adouble, double, double);
349         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(adouble, float, double, double, double);
350         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(adouble, double, double, double, double);
351         CHECK_TEMPLATED_RESULT_TYPE_FOR_ALL(adouble, long double, long double, long double, long double);
352 
353     }
354 }
355