xtensor
 
Loading...
Searching...
No Matches
xmath.hpp
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
13
14#ifndef XTENSOR_MATH_HPP
15#define XTENSOR_MATH_HPP
16
17#include <algorithm>
18#include <array>
19#include <cmath>
20#include <complex>
21#include <type_traits>
22
23#include <xtl/xcomplex.hpp>
24#include <xtl/xsequence.hpp>
25#include <xtl/xtype_traits.hpp>
26
27#include "../core/xeval.hpp"
28#include "../core/xoperation.hpp"
29#include "../core/xtensor_config.hpp"
30#include "../misc/xmanipulation.hpp"
31#include "../reducers/xaccumulator.hpp"
32#include "../reducers/xreducer.hpp"
33#include "../views/xslice.hpp"
34#include "../views/xstrided_view.hpp"
35
36namespace xt
37{
38 template <class T = double>
39 struct numeric_constants
40 {
41 static constexpr T PI = 3.141592653589793238463;
42 static constexpr T PI_2 = 1.57079632679489661923;
43 static constexpr T PI_4 = 0.785398163397448309616;
44 static constexpr T D_1_PI = 0.318309886183790671538;
45 static constexpr T D_2_PI = 0.636619772367581343076;
46 static constexpr T D_2_SQRTPI = 1.12837916709551257390;
47 static constexpr T SQRT2 = 1.41421356237309504880;
48 static constexpr T SQRT1_2 = 0.707106781186547524401;
49 static constexpr T E = 2.71828182845904523536;
50 static constexpr T LOG2E = 1.44269504088896340736;
51 static constexpr T LOG10E = 0.434294481903251827651;
52 static constexpr T LN2 = 0.693147180559945309417;
53 };
54
55 /***********
56 * Helpers *
57 ***********/
58
59#define XTENSOR_UNSIGNED_ABS_FUNC(T) \
60 constexpr inline T abs(const T& x) \
61 { \
62 return x; \
63 }
64
65#define XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, T) \
66 constexpr inline bool FUNC_NAME(const T& /*x*/) noexcept \
67 { \
68 return RETURN_VAL; \
69 }
70
71#define XTENSOR_INT_SPECIALIZATION(FUNC_NAME, RETURN_VAL) \
72 XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, char); \
73 XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, short); \
74 XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, int); \
75 XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, long); \
76 XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, long long); \
77 XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, unsigned char); \
78 XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, unsigned short); \
79 XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, unsigned int); \
80 XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, unsigned long); \
81 XTENSOR_INT_SPECIALIZATION_IMPL(FUNC_NAME, RETURN_VAL, unsigned long long);
82
83
84#define XTENSOR_UNARY_MATH_FUNCTOR(NAME) \
85 struct NAME##_fun \
86 { \
87 template <class T> \
88 constexpr auto operator()(const T& arg) const \
89 { \
90 using math::NAME; \
91 return NAME(arg); \
92 } \
93 template <class B> \
94 constexpr auto simd_apply(const B& arg) const \
95 { \
96 using math::NAME; \
97 return NAME(arg); \
98 } \
99 }
100
101#define XTENSOR_UNARY_MATH_FUNCTOR_COMPLEX_REDUCING(NAME) \
102 struct NAME##_fun \
103 { \
104 template <class T> \
105 constexpr auto operator()(const T& arg) const \
106 { \
107 using math::NAME; \
108 return NAME(arg); \
109 } \
110 template <class B> \
111 constexpr auto simd_apply(const B& arg) const \
112 { \
113 using math::NAME; \
114 return NAME(arg); \
115 } \
116 }
117
118#define XTENSOR_BINARY_MATH_FUNCTOR(NAME) \
119 struct NAME##_fun \
120 { \
121 template <class T1, class T2> \
122 constexpr auto operator()(const T1& arg1, const T2& arg2) const \
123 { \
124 using math::NAME; \
125 return NAME(arg1, arg2); \
126 } \
127 template <class B> \
128 constexpr auto simd_apply(const B& arg1, const B& arg2) const \
129 { \
130 using math::NAME; \
131 return NAME(arg1, arg2); \
132 } \
133 }
134
135#define XTENSOR_TERNARY_MATH_FUNCTOR(NAME) \
136 struct NAME##_fun \
137 { \
138 template <class T1, class T2, class T3> \
139 constexpr auto operator()(const T1& arg1, const T2& arg2, const T3& arg3) const \
140 { \
141 using math::NAME; \
142 return NAME(arg1, arg2, arg3); \
143 } \
144 template <class B> \
145 auto simd_apply(const B& arg1, const B& arg2, const B& arg3) const \
146 { \
147 using math::NAME; \
148 return NAME(arg1, arg2, arg3); \
149 } \
150 }
151
152 namespace math
153 {
154 using std::abs;
155 using std::fabs;
156
157 using std::acos;
158 using std::asin;
159 using std::atan;
160 using std::cos;
161 using std::sin;
162 using std::tan;
163
164 using std::acosh;
165 using std::asinh;
166 using std::atanh;
167 using std::cosh;
168 using std::sinh;
169 using std::tanh;
170
171 using std::cbrt;
172 using std::sqrt;
173
174 using std::exp;
175 using std::exp2;
176 using std::expm1;
177 using std::ilogb;
178 using std::log;
179 using std::log10;
180 using std::log1p;
181 using std::log2;
182 using std::logb;
183
184 using std::ceil;
185 using std::floor;
186 using std::llround;
187 using std::lround;
188 using std::nearbyint;
189 using std::remainder;
190 using std::rint;
191 using std::round;
192 using std::trunc;
193
194 using std::erf;
195 using std::erfc;
196 using std::lgamma;
197 using std::tgamma;
198
199 using std::arg;
200 using std::conj;
201 using std::imag;
202 using std::real;
203
204 using std::atan2;
205
206// copysign is not in the std namespace for MSVC
207#if !defined(_MSC_VER)
208 using std::copysign;
209#endif
210 using std::fdim;
211 using std::fmax;
212 using std::fmin;
213 using std::fmod;
214 using std::hypot;
215 using std::pow;
216
217 using std::fma;
218 using std::fpclassify;
219
220 // Overload isinf, isnan and isfinite because glibc implementation
221 // might return int instead of bool and the SIMD detection requires
222 // bool return type.
223 template <class T>
224 inline std::enable_if_t<xtl::is_arithmetic<T>::value, bool> isinf(const T& t)
225 {
226 return bool(std::isinf(t));
227 }
228
229 template <class T>
230 inline std::enable_if_t<xtl::is_arithmetic<T>::value, bool> isnan(const T& t)
231 {
232 return bool(std::isnan(t));
233 }
234
235 template <class T>
236 inline std::enable_if_t<xtl::is_arithmetic<T>::value, bool> isfinite(const T& t)
237 {
238 return bool(std::isfinite(t));
239 }
240
241 // Overload isinf, isnan and isfinite for complex datatypes,
242 // following the Python specification:
243 template <class T>
244 inline bool isinf(const std::complex<T>& c)
245 {
246 return std::isinf(std::real(c)) || std::isinf(std::imag(c));
247 }
248
249 template <class T>
250 inline bool isnan(const std::complex<T>& c)
251 {
252 return std::isnan(std::real(c)) || std::isnan(std::imag(c));
253 }
254
255 template <class T>
256 inline bool isfinite(const std::complex<T>& c)
257 {
258 return !isinf(c) && !isnan(c);
259 }
260
261 // VS2015 STL defines isnan, isinf and isfinite as template
262 // functions, breaking ADL.
263#if defined(_WIN32) && defined(XTENSOR_USE_XSIMD)
264 /*template <class T, class A>
265 inline xsimd::batch_bool<T, A> isinf(const xsimd::batch<T, A>& b)
266 {
267 return xsimd::isinf(b);
268 }
269 template <class T, class A>
270 inline xsimd::batch_bool<T, A> isnan(const xsimd::batch<T, A>& b)
271 {
272 return xsimd::isnan(b);
273 }
274 template <class T, class A>
275 inline xsimd::batch_bool<T, A> isfinite(const xsimd::batch<T, A>& b)
276 {
277 return xsimd::isfinite(b);
278 }*/
279#endif
280 // The following specializations are needed to avoid 'ambiguous overload' errors,
281 // whereas 'unsigned char' and 'unsigned short' are automatically converted to 'int'.
282 // we're still adding those functions to silence warnings
283 XTENSOR_UNSIGNED_ABS_FUNC(unsigned char)
284 XTENSOR_UNSIGNED_ABS_FUNC(unsigned short)
285 XTENSOR_UNSIGNED_ABS_FUNC(unsigned int)
286 XTENSOR_UNSIGNED_ABS_FUNC(unsigned long)
287 XTENSOR_UNSIGNED_ABS_FUNC(unsigned long long)
288
289#ifdef _WIN32
290 XTENSOR_INT_SPECIALIZATION(isinf, false);
291 XTENSOR_INT_SPECIALIZATION(isnan, false);
292 XTENSOR_INT_SPECIALIZATION(isfinite, true);
293#endif
294
295 XTENSOR_UNARY_MATH_FUNCTOR_COMPLEX_REDUCING(abs);
296
297 XTENSOR_UNARY_MATH_FUNCTOR(fabs);
298 XTENSOR_BINARY_MATH_FUNCTOR(fmod);
299 XTENSOR_BINARY_MATH_FUNCTOR(remainder);
300 XTENSOR_TERNARY_MATH_FUNCTOR(fma);
301 XTENSOR_BINARY_MATH_FUNCTOR(fmax);
302 XTENSOR_BINARY_MATH_FUNCTOR(fmin);
303 XTENSOR_BINARY_MATH_FUNCTOR(fdim);
304 XTENSOR_UNARY_MATH_FUNCTOR(exp);
305 XTENSOR_UNARY_MATH_FUNCTOR(exp2);
306 XTENSOR_UNARY_MATH_FUNCTOR(expm1);
307 XTENSOR_UNARY_MATH_FUNCTOR(log);
308 XTENSOR_UNARY_MATH_FUNCTOR(log10);
309 XTENSOR_UNARY_MATH_FUNCTOR(log2);
310 XTENSOR_UNARY_MATH_FUNCTOR(log1p);
311 XTENSOR_BINARY_MATH_FUNCTOR(pow);
312 XTENSOR_UNARY_MATH_FUNCTOR(sqrt);
313 XTENSOR_UNARY_MATH_FUNCTOR(cbrt);
314 XTENSOR_BINARY_MATH_FUNCTOR(hypot);
315 XTENSOR_UNARY_MATH_FUNCTOR(sin);
316 XTENSOR_UNARY_MATH_FUNCTOR(cos);
317 XTENSOR_UNARY_MATH_FUNCTOR(tan);
318 XTENSOR_UNARY_MATH_FUNCTOR(asin);
319 XTENSOR_UNARY_MATH_FUNCTOR(acos);
320 XTENSOR_UNARY_MATH_FUNCTOR(atan);
321 XTENSOR_BINARY_MATH_FUNCTOR(atan2);
322 XTENSOR_UNARY_MATH_FUNCTOR(sinh);
323 XTENSOR_UNARY_MATH_FUNCTOR(cosh);
324 XTENSOR_UNARY_MATH_FUNCTOR(tanh);
325 XTENSOR_UNARY_MATH_FUNCTOR(asinh);
326 XTENSOR_UNARY_MATH_FUNCTOR(acosh);
327 XTENSOR_UNARY_MATH_FUNCTOR(atanh);
328 XTENSOR_UNARY_MATH_FUNCTOR(erf);
329 XTENSOR_UNARY_MATH_FUNCTOR(erfc);
330 XTENSOR_UNARY_MATH_FUNCTOR(tgamma);
331 XTENSOR_UNARY_MATH_FUNCTOR(lgamma);
332 XTENSOR_UNARY_MATH_FUNCTOR(ceil);
333 XTENSOR_UNARY_MATH_FUNCTOR(floor);
334 XTENSOR_UNARY_MATH_FUNCTOR(trunc);
335 XTENSOR_UNARY_MATH_FUNCTOR(round);
336 XTENSOR_UNARY_MATH_FUNCTOR(nearbyint);
337 XTENSOR_UNARY_MATH_FUNCTOR(rint);
338 XTENSOR_UNARY_MATH_FUNCTOR(isfinite);
339 XTENSOR_UNARY_MATH_FUNCTOR(isinf);
340 XTENSOR_UNARY_MATH_FUNCTOR(isnan);
341 XTENSOR_UNARY_MATH_FUNCTOR(conj);
342 }
343
344#undef XTENSOR_UNARY_MATH_FUNCTOR
345#undef XTENSOR_BINARY_MATH_FUNCTOR
346#undef XTENSOR_TERNARY_MATH_FUNCTOR
347#undef XTENSOR_UNARY_MATH_FUNCTOR_COMPLEX_REDUCING
348#undef XTENSOR_UNSIGNED_ABS_FUNC
349
350 namespace detail
351 {
352 template <class R, class T>
353 std::enable_if_t<!has_iterator_interface<R>::value, R> fill_init(T init)
354 {
355 return R(init);
356 }
357
358 template <class R, class T>
359 std::enable_if_t<has_iterator_interface<R>::value, R> fill_init(T init)
360 {
361 R result;
362 std::fill(std::begin(result), std::end(result), init);
363 return result;
364 }
365 }
366
367#define XTENSOR_REDUCER_FUNCTION(NAME, FUNCTOR, INIT_VALUE_TYPE, INIT) \
368 template < \
369 class T = void, \
370 class E, \
371 class X, \
372 class EVS = DEFAULT_STRATEGY_REDUCERS, \
373 XTL_REQUIRES(std::negation<is_reducer_options<X>>, std::negation<xtl::is_integral<std::decay_t<X>>>)> \
374 inline auto NAME(E&& e, X&& axes, EVS es = EVS()) \
375 { \
376 using init_value_type = std::conditional_t<std::is_same<T, void>::value, INIT_VALUE_TYPE, T>; \
377 using functor_type = FUNCTOR; \
378 using init_value_fct = xt::const_value<init_value_type>; \
379 return xt::reduce( \
380 make_xreducer_functor(functor_type(), init_value_fct(detail::fill_init<init_value_type>(INIT))), \
381 std::forward<E>(e), \
382 std::forward<X>(axes), \
383 es \
384 ); \
385 } \
386 \
387 template < \
388 class T = void, \
389 class E, \
390 class X, \
391 class EVS = DEFAULT_STRATEGY_REDUCERS, \
392 XTL_REQUIRES(std::negation<is_reducer_options<X>>, xtl::is_integral<std::decay_t<X>>)> \
393 inline auto NAME(E&& e, X axis, EVS es = EVS()) \
394 { \
395 return NAME(std::forward<E>(e), {axis}, es); \
396 } \
397 \
398 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)> \
399 inline auto NAME(E&& e, EVS es = EVS()) \
400 { \
401 using init_value_type = std::conditional_t<std::is_same<T, void>::value, INIT_VALUE_TYPE, T>; \
402 using functor_type = FUNCTOR; \
403 using init_value_fct = xt::const_value<init_value_type>; \
404 return xt::reduce( \
405 make_xreducer_functor(functor_type(), init_value_fct(detail::fill_init<init_value_type>(INIT))), \
406 std::forward<E>(e), \
407 es \
408 ); \
409 } \
410 \
411 template <class T = void, class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS> \
412 inline auto NAME(E&& e, const I(&axes)[N], EVS es = EVS()) \
413 { \
414 using init_value_type = std::conditional_t<std::is_same<T, void>::value, INIT_VALUE_TYPE, T>; \
415 using functor_type = FUNCTOR; \
416 using init_value_fct = xt::const_value<init_value_type>; \
417 return xt::reduce( \
418 make_xreducer_functor(functor_type(), init_value_fct(detail::fill_init<init_value_type>(INIT))), \
419 std::forward<E>(e), \
420 axes, \
421 es \
422 ); \
423 }
424
425 /*******************
426 * basic functions *
427 *******************/
428
432
442 template <class E>
443 inline auto abs(E&& e) noexcept -> detail::xfunction_type_t<math::abs_fun, E>
444 {
445 return detail::make_xfunction<math::abs_fun>(std::forward<E>(e));
446 }
447
457 template <class E>
458 inline auto fabs(E&& e) noexcept -> detail::xfunction_type_t<math::fabs_fun, E>
459 {
460 return detail::make_xfunction<math::fabs_fun>(std::forward<E>(e));
461 }
462
474 template <class E1, class E2>
475 inline auto fmod(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::fmod_fun, E1, E2>
476 {
477 return detail::make_xfunction<math::fmod_fun>(std::forward<E1>(e1), std::forward<E2>(e2));
478 }
479
491 template <class E1, class E2>
492 inline auto remainder(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::remainder_fun, E1, E2>
493 {
494 return detail::make_xfunction<math::remainder_fun>(std::forward<E1>(e1), std::forward<E2>(e2));
495 }
496
509 template <class E1, class E2, class E3>
510 inline auto fma(E1&& e1, E2&& e2, E3&& e3) noexcept -> detail::xfunction_type_t<math::fma_fun, E1, E2, E3>
511 {
512 return detail::make_xfunction<math::fma_fun>(
513 std::forward<E1>(e1),
514 std::forward<E2>(e2),
515 std::forward<E3>(e3)
516 );
517 }
518
530 template <class E1, class E2>
531 inline auto fmax(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::fmax_fun, E1, E2>
532 {
533 return detail::make_xfunction<math::fmax_fun>(std::forward<E1>(e1), std::forward<E2>(e2));
534 }
535
547 template <class E1, class E2>
548 inline auto fmin(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::fmin_fun, E1, E2>
549 {
550 return detail::make_xfunction<math::fmin_fun>(std::forward<E1>(e1), std::forward<E2>(e2));
551 }
552
564 template <class E1, class E2>
565 inline auto fdim(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::fdim_fun, E1, E2>
566 {
567 return detail::make_xfunction<math::fdim_fun>(std::forward<E1>(e1), std::forward<E2>(e2));
568 }
569
570 namespace math
571 {
572 template <class T = void>
573 struct minimum
574 {
575 template <class A1, class A2>
576 constexpr auto operator()(const A1& t1, const A2& t2) const noexcept
577 {
578 return xtl::select(t1 < t2, t1, t2);
579 }
580
581 template <class A1, class A2>
582 constexpr auto simd_apply(const A1& t1, const A2& t2) const noexcept
583 {
584 return xt_simd::select(t1 < t2, t1, t2);
585 }
586 };
587
588 template <class T = void>
589 struct maximum
590 {
591 template <class A1, class A2>
592 constexpr auto operator()(const A1& t1, const A2& t2) const noexcept
593 {
594 return xtl::select(t1 > t2, t1, t2);
595 }
596
597 template <class A1, class A2>
598 constexpr auto simd_apply(const A1& t1, const A2& t2) const noexcept
599 {
600 return xt_simd::select(t1 > t2, t1, t2);
601 }
602 };
603
604 struct clamp_fun
605 {
606 template <class A1, class A2, class A3>
607 constexpr auto operator()(const A1& v, const A2& lo, const A3& hi) const
608 {
609 return xtl::select(lo < hi, xtl::select(v < lo, lo, xtl::select(hi < v, hi, v)), hi);
610 }
611
612 template <class A1, class A2, class A3>
613 constexpr auto simd_apply(const A1& v, const A2& lo, const A3& hi) const
614 {
615 return xt_simd::select(lo < hi, xt_simd::select(v < lo, lo, xt_simd::select(hi < v, hi, v)), hi);
616 }
617 };
618
619 struct deg2rad
620 {
621 template <class A, std::enable_if_t<xtl::is_integral<A>::value, int> = 0>
622 constexpr double operator()(const A& a) const noexcept
623 {
624 return a * xt::numeric_constants<double>::PI / 180.0;
625 }
626
627 template <class A, std::enable_if_t<std::is_floating_point<A>::value, int> = 0>
628 constexpr auto operator()(const A& a) const noexcept
629 {
630 return a * xt::numeric_constants<A>::PI / A(180.0);
631 }
632
633 template <class A, std::enable_if_t<xtl::is_integral<A>::value, int> = 0>
634 constexpr double simd_apply(const A& a) const noexcept
635 {
636 return a * xt::numeric_constants<double>::PI / 180.0;
637 }
638
639 template <class A, std::enable_if_t<std::is_floating_point<A>::value, int> = 0>
640 constexpr auto simd_apply(const A& a) const noexcept
641 {
642 return a * xt::numeric_constants<A>::PI / A(180.0);
643 }
644 };
645
646 struct rad2deg
647 {
648 template <class A, std::enable_if_t<xtl::is_integral<A>::value, int> = 0>
649 constexpr double operator()(const A& a) const noexcept
650 {
651 return a * 180.0 / xt::numeric_constants<double>::PI;
652 }
653
654 template <class A, std::enable_if_t<std::is_floating_point<A>::value, int> = 0>
655 constexpr auto operator()(const A& a) const noexcept
656 {
657 return a * A(180.0) / xt::numeric_constants<A>::PI;
658 }
659
660 template <class A, std::enable_if_t<xtl::is_integral<A>::value, int> = 0>
661 constexpr double simd_apply(const A& a) const noexcept
662 {
663 return a * 180.0 / xt::numeric_constants<double>::PI;
664 }
665
666 template <class A, std::enable_if_t<std::is_floating_point<A>::value, int> = 0>
667 constexpr auto simd_apply(const A& a) const noexcept
668 {
669 return a * A(180.0) / xt::numeric_constants<A>::PI;
670 }
671 };
672 }
673
683 template <class E>
684 inline auto deg2rad(E&& e) noexcept -> detail::xfunction_type_t<math::deg2rad, E>
685 {
686 return detail::make_xfunction<math::deg2rad>(std::forward<E>(e));
687 }
688
698 template <class E>
699 inline auto radians(E&& e) noexcept -> detail::xfunction_type_t<math::deg2rad, E>
700 {
701 return detail::make_xfunction<math::deg2rad>(std::forward<E>(e));
702 }
703
713 template <class E>
714 inline auto rad2deg(E&& e) noexcept -> detail::xfunction_type_t<math::rad2deg, E>
715 {
716 return detail::make_xfunction<math::rad2deg>(std::forward<E>(e));
717 }
718
728 template <class E>
729 inline auto degrees(E&& e) noexcept -> detail::xfunction_type_t<math::rad2deg, E>
730 {
731 return detail::make_xfunction<math::rad2deg>(std::forward<E>(e));
732 }
733
744 template <class E1, class E2>
745 inline auto maximum(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::maximum<void>, E1, E2>
746 {
747 return detail::make_xfunction<math::maximum<void>>(std::forward<E1>(e1), std::forward<E2>(e2));
748 }
749
760 template <class E1, class E2>
761 inline auto minimum(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::minimum<void>, E1, E2>
762 {
763 return detail::make_xfunction<math::minimum<void>>(std::forward<E1>(e1), std::forward<E2>(e2));
764 }
765
777 XTENSOR_REDUCER_FUNCTION(
778 amax,
779 math::maximum<void>,
780 typename std::decay_t<E>::value_type,
781 std::numeric_limits<xvalue_type_t<std::decay_t<E>>>::lowest()
782 )
783
784
795 XTENSOR_REDUCER_FUNCTION(
796 amin,
797 math::minimum<void>,
798 typename std::decay_t<E>::value_type,
799 std::numeric_limits<xvalue_type_t<std::decay_t<E>>>::max()
800 )
801
814 template <class E1, class E2, class E3>
815 inline auto clip(E1&& e1, E2&& lo, E3&& hi) noexcept
816 -> detail::xfunction_type_t<math::clamp_fun, E1, E2, E3>
817 {
818 return detail::make_xfunction<math::clamp_fun>(
819 std::forward<E1>(e1),
820 std::forward<E2>(lo),
821 std::forward<E3>(hi)
822 );
823 }
824
825 namespace math
826 {
827 template <class T>
828 struct sign_impl
829 {
830 template <class XT = T>
831 static constexpr std::enable_if_t<xtl::is_signed<XT>::value, T> run(T x)
832 {
833 return std::isnan(x) ? std::numeric_limits<T>::quiet_NaN()
834 : x == 0 ? T(copysign(T(0), x))
835 : T(copysign(T(1), x));
836 }
837
838 template <class XT = T>
839 static constexpr std::enable_if_t<xtl::is_complex<XT>::value, T> run(T x)
840 {
841 return T(
842 sign_impl<typename T::value_type>::run(
843 (x.real() != typename T::value_type(0)) ? x.real() : x.imag()
844 ),
845 0
846 );
847 }
848
849 template <class XT = T>
850 static constexpr std::enable_if_t<std::is_unsigned<XT>::value, T> run(T x)
851 {
852 return T(x > T(0));
853 }
854 };
855
856 struct sign_fun
857 {
858 template <class T>
859 constexpr auto operator()(const T& x) const
860 {
861 return sign_impl<T>::run(x);
862 }
863 };
864 }
865
876 template <class E>
877 inline auto sign(E&& e) noexcept -> detail::xfunction_type_t<math::sign_fun, E>
878 {
879 return detail::make_xfunction<math::sign_fun>(std::forward<E>(e));
880 }
881
882 /*************************
883 * exponential functions *
884 *************************/
885
889
899 template <class E>
900 inline auto exp(E&& e) noexcept -> detail::xfunction_type_t<math::exp_fun, E>
901 {
902 return detail::make_xfunction<math::exp_fun>(std::forward<E>(e));
903 }
904
914 template <class E>
915 inline auto exp2(E&& e) noexcept -> detail::xfunction_type_t<math::exp2_fun, E>
916 {
917 return detail::make_xfunction<math::exp2_fun>(std::forward<E>(e));
918 }
919
929 template <class E>
930 inline auto expm1(E&& e) noexcept -> detail::xfunction_type_t<math::expm1_fun, E>
931 {
932 return detail::make_xfunction<math::expm1_fun>(std::forward<E>(e));
933 }
934
944 template <class E>
945 inline auto log(E&& e) noexcept -> detail::xfunction_type_t<math::log_fun, E>
946 {
947 return detail::make_xfunction<math::log_fun>(std::forward<E>(e));
948 }
949
959 template <class E>
960 inline auto log10(E&& e) noexcept -> detail::xfunction_type_t<math::log10_fun, E>
961 {
962 return detail::make_xfunction<math::log10_fun>(std::forward<E>(e));
963 }
964
974 template <class E>
975 inline auto log2(E&& e) noexcept -> detail::xfunction_type_t<math::log2_fun, E>
976 {
977 return detail::make_xfunction<math::log2_fun>(std::forward<E>(e));
978 }
979
989 template <class E>
990 inline auto log1p(E&& e) noexcept -> detail::xfunction_type_t<math::log1p_fun, E>
991 {
992 return detail::make_xfunction<math::log1p_fun>(std::forward<E>(e));
993 }
994
995 /*******************
996 * power functions *
997 *******************/
998
1002
1014 template <class E1, class E2>
1015 inline auto pow(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::pow_fun, E1, E2>
1016 {
1017 return detail::make_xfunction<math::pow_fun>(std::forward<E1>(e1), std::forward<E2>(e2));
1018 }
1019
1020 namespace detail
1021 {
1022 template <class F, class... T, typename = decltype(std::declval<F>()(std::declval<T>()...))>
1023 std::true_type supports_test(const F&, const T&...);
1024 std::false_type supports_test(...);
1025
1026 template <class... T>
1027 struct supports;
1028
1029 template <class F, class... T>
1030 struct supports<F(T...)> : decltype(supports_test(std::declval<F>(), std::declval<T>()...))
1031 {
1032 };
1033
1034 template <class F>
1035 struct lambda_adapt
1036 {
1037 explicit lambda_adapt(F&& lmbd)
1038 : m_lambda(std::move(lmbd))
1039 {
1040 }
1041
1042 template <class... T>
1043 auto operator()(T... args) const
1044 {
1045 return m_lambda(args...);
1046 }
1047
1048 template <class... T, XTL_REQUIRES(detail::supports<F(T...)>)>
1049 auto simd_apply(T... args) const
1050 {
1051 return m_lambda(args...);
1052 }
1053
1054 F m_lambda;
1055 };
1056 }
1057
1084 template <class F, class... E>
1085 inline auto make_lambda_xfunction(F&& lambda, E&&... args)
1086 {
1087 using xfunction_type = typename detail::xfunction_type<detail::lambda_adapt<F>, E...>::type;
1088 return xfunction_type(detail::lambda_adapt<F>(std::forward<F>(lambda)), std::forward<E>(args)...);
1089 }
1090
1100 template <class E1>
1101 inline auto square(E1&& e1) noexcept
1102 {
1103 auto fnct = [](auto x) -> decltype(x * x)
1104 {
1105 return x * x;
1106 };
1107 return make_lambda_xfunction(std::move(fnct), std::forward<E1>(e1));
1108 }
1109
1119 template <class E1>
1120 inline auto cube(E1&& e1) noexcept
1121 {
1122 auto fnct = [](auto x) -> decltype(x * x * x)
1123 {
1124 return x * x * x;
1125 };
1126 return make_lambda_xfunction(std::move(fnct), std::forward<E1>(e1));
1127 }
1128
1129 namespace detail
1130 {
1131 // Thanks to Matt Pharr in http://pbrt.org/hair.pdf
1132 template <std::size_t N>
1133 struct pow_impl;
1134
1135 template <std::size_t N>
1136 struct pow_impl
1137 {
1138 template <class T>
1139 auto operator()(T v) const -> decltype(v * v)
1140 {
1141 T temp = pow_impl<N / 2>{}(v);
1142 return temp * temp * pow_impl<N & 1>{}(v);
1143 }
1144 };
1145
1146 template <>
1147 struct pow_impl<1>
1148 {
1149 template <class T>
1150 auto operator()(T v) const -> T
1151 {
1152 return v;
1153 }
1154 };
1155
1156 template <>
1157 struct pow_impl<0>
1158 {
1159 template <class T>
1160 auto operator()(T /*v*/) const -> T
1161 {
1162 return T(1);
1163 }
1164 };
1165 }
1166
1184 template <std::size_t N, class E>
1185 inline auto pow(E&& e) noexcept
1186 {
1187 static_assert(N > 0, "integer power cannot be negative");
1188 return make_lambda_xfunction(detail::pow_impl<N>{}, std::forward<E>(e));
1189 }
1190
1200 template <class E>
1201 inline auto sqrt(E&& e) noexcept -> detail::xfunction_type_t<math::sqrt_fun, E>
1202 {
1203 return detail::make_xfunction<math::sqrt_fun>(std::forward<E>(e));
1204 }
1205
1215 template <class E>
1216 inline auto cbrt(E&& e) noexcept -> detail::xfunction_type_t<math::cbrt_fun, E>
1217 {
1218 return detail::make_xfunction<math::cbrt_fun>(std::forward<E>(e));
1219 }
1220
1233 template <class E1, class E2>
1234 inline auto hypot(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::hypot_fun, E1, E2>
1235 {
1236 return detail::make_xfunction<math::hypot_fun>(std::forward<E1>(e1), std::forward<E2>(e2));
1237 }
1238
1239 /***************************
1240 * trigonometric functions *
1241 ***************************/
1242
1246
1256 template <class E>
1257 inline auto sin(E&& e) noexcept -> detail::xfunction_type_t<math::sin_fun, E>
1258 {
1259 return detail::make_xfunction<math::sin_fun>(std::forward<E>(e));
1260 }
1261
1271 template <class E>
1272 inline auto cos(E&& e) noexcept -> detail::xfunction_type_t<math::cos_fun, E>
1273 {
1274 return detail::make_xfunction<math::cos_fun>(std::forward<E>(e));
1275 }
1276
1286 template <class E>
1287 inline auto tan(E&& e) noexcept -> detail::xfunction_type_t<math::tan_fun, E>
1288 {
1289 return detail::make_xfunction<math::tan_fun>(std::forward<E>(e));
1290 }
1291
1301 template <class E>
1302 inline auto asin(E&& e) noexcept -> detail::xfunction_type_t<math::asin_fun, E>
1303 {
1304 return detail::make_xfunction<math::asin_fun>(std::forward<E>(e));
1305 }
1306
1316 template <class E>
1317 inline auto acos(E&& e) noexcept -> detail::xfunction_type_t<math::acos_fun, E>
1318 {
1319 return detail::make_xfunction<math::acos_fun>(std::forward<E>(e));
1320 }
1321
1331 template <class E>
1332 inline auto atan(E&& e) noexcept -> detail::xfunction_type_t<math::atan_fun, E>
1333 {
1334 return detail::make_xfunction<math::atan_fun>(std::forward<E>(e));
1335 }
1336
1349 template <class E1, class E2>
1350 inline auto atan2(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::atan2_fun, E1, E2>
1351 {
1352 return detail::make_xfunction<math::atan2_fun>(std::forward<E1>(e1), std::forward<E2>(e2));
1353 }
1354
1355 /************************
1356 * hyperbolic functions *
1357 ************************/
1358
1362
1372 template <class E>
1373 inline auto sinh(E&& e) noexcept -> detail::xfunction_type_t<math::sinh_fun, E>
1374 {
1375 return detail::make_xfunction<math::sinh_fun>(std::forward<E>(e));
1376 }
1377
1387 template <class E>
1388 inline auto cosh(E&& e) noexcept -> detail::xfunction_type_t<math::cosh_fun, E>
1389 {
1390 return detail::make_xfunction<math::cosh_fun>(std::forward<E>(e));
1391 }
1392
1402 template <class E>
1403 inline auto tanh(E&& e) noexcept -> detail::xfunction_type_t<math::tanh_fun, E>
1404 {
1405 return detail::make_xfunction<math::tanh_fun>(std::forward<E>(e));
1406 }
1407
1417 template <class E>
1418 inline auto asinh(E&& e) noexcept -> detail::xfunction_type_t<math::asinh_fun, E>
1419 {
1420 return detail::make_xfunction<math::asinh_fun>(std::forward<E>(e));
1421 }
1422
1432 template <class E>
1433 inline auto acosh(E&& e) noexcept -> detail::xfunction_type_t<math::acosh_fun, E>
1434 {
1435 return detail::make_xfunction<math::acosh_fun>(std::forward<E>(e));
1436 }
1437
1447 template <class E>
1448 inline auto atanh(E&& e) noexcept -> detail::xfunction_type_t<math::atanh_fun, E>
1449 {
1450 return detail::make_xfunction<math::atanh_fun>(std::forward<E>(e));
1451 }
1452
1453 /*****************************
1454 * error and gamma functions *
1455 *****************************/
1456
1460
1470 template <class E>
1471 inline auto erf(E&& e) noexcept -> detail::xfunction_type_t<math::erf_fun, E>
1472 {
1473 return detail::make_xfunction<math::erf_fun>(std::forward<E>(e));
1474 }
1475
1485 template <class E>
1486 inline auto erfc(E&& e) noexcept -> detail::xfunction_type_t<math::erfc_fun, E>
1487 {
1488 return detail::make_xfunction<math::erfc_fun>(std::forward<E>(e));
1489 }
1490
1500 template <class E>
1501 inline auto tgamma(E&& e) noexcept -> detail::xfunction_type_t<math::tgamma_fun, E>
1502 {
1503 return detail::make_xfunction<math::tgamma_fun>(std::forward<E>(e));
1504 }
1505
1515 template <class E>
1516 inline auto lgamma(E&& e) noexcept -> detail::xfunction_type_t<math::lgamma_fun, E>
1517 {
1518 return detail::make_xfunction<math::lgamma_fun>(std::forward<E>(e));
1519 }
1520
1521 /*********************************************
1522 * nearest integer floating point operations *
1523 *********************************************/
1524
1528
1538 template <class E>
1539 inline auto ceil(E&& e) noexcept -> detail::xfunction_type_t<math::ceil_fun, E>
1540 {
1541 return detail::make_xfunction<math::ceil_fun>(std::forward<E>(e));
1542 }
1543
1553 template <class E>
1554 inline auto floor(E&& e) noexcept -> detail::xfunction_type_t<math::floor_fun, E>
1555 {
1556 return detail::make_xfunction<math::floor_fun>(std::forward<E>(e));
1557 }
1558
1568 template <class E>
1569 inline auto trunc(E&& e) noexcept -> detail::xfunction_type_t<math::trunc_fun, E>
1570 {
1571 return detail::make_xfunction<math::trunc_fun>(std::forward<E>(e));
1572 }
1573
1584 template <class E>
1585 inline auto round(E&& e) noexcept -> detail::xfunction_type_t<math::round_fun, E>
1586 {
1587 return detail::make_xfunction<math::round_fun>(std::forward<E>(e));
1588 }
1589
1600 template <class E>
1601 inline auto nearbyint(E&& e) noexcept -> detail::xfunction_type_t<math::nearbyint_fun, E>
1602 {
1603 return detail::make_xfunction<math::nearbyint_fun>(std::forward<E>(e));
1604 }
1605
1616 template <class E>
1617 inline auto rint(E&& e) noexcept -> detail::xfunction_type_t<math::rint_fun, E>
1618 {
1619 return detail::make_xfunction<math::rint_fun>(std::forward<E>(e));
1620 }
1621
1622 /****************************
1623 * classification functions *
1624 ****************************/
1625
1629
1639 template <class E>
1640 inline auto isfinite(E&& e) noexcept -> detail::xfunction_type_t<math::isfinite_fun, E>
1641 {
1642 return detail::make_xfunction<math::isfinite_fun>(std::forward<E>(e));
1643 }
1644
1654 template <class E>
1655 inline auto isinf(E&& e) noexcept -> detail::xfunction_type_t<math::isinf_fun, E>
1656 {
1657 return detail::make_xfunction<math::isinf_fun>(std::forward<E>(e));
1658 }
1659
1669 template <class E>
1670 inline auto isnan(E&& e) noexcept -> detail::xfunction_type_t<math::isnan_fun, E>
1671 {
1672 return detail::make_xfunction<math::isnan_fun>(std::forward<E>(e));
1673 }
1674
1675 namespace detail
1676 {
1677 template <class FUNCTOR, class T, std::size_t... Is>
1678 inline auto get_functor(T&& args, std::index_sequence<Is...>)
1679 {
1680 return FUNCTOR(std::get<Is>(args)...);
1681 }
1682
1683 template <class F, class... A, class... E>
1684 inline auto make_xfunction(std::tuple<A...>&& f_args, E&&... e) noexcept
1685 {
1686 using functor_type = F;
1687 using expression_tag = xexpression_tag_t<E...>;
1688 using type = select_xfunction_expression_t<expression_tag, functor_type, const_xclosure_t<E>...>;
1689 auto functor = get_functor<functor_type>(
1690 std::forward<std::tuple<A...>>(f_args),
1691 std::make_index_sequence<sizeof...(A)>{}
1692 );
1693 return type(std::move(functor), std::forward<E>(e)...);
1694 }
1695
1696 struct isclose
1697 {
1698 using result_type = bool;
1699
1700 isclose(double rtol, double atol, bool equal_nan)
1701 : m_rtol(rtol)
1702 , m_atol(atol)
1703 , m_equal_nan(equal_nan)
1704 {
1705 }
1706
1707 template <class A1, class A2>
1708 bool operator()(const A1& a, const A2& b) const
1709 {
1710 using internal_type = xtl::promote_type_t<A1, A2, double>;
1711 if (math::isnan(a) && math::isnan(b))
1712 {
1713 return m_equal_nan;
1714 }
1715 if (math::isinf(a) && math::isinf(b))
1716 {
1717 // check for both infinity signs equal
1718 return a == b;
1719 }
1720 auto d = math::abs(internal_type(a) - internal_type(b));
1721 return d <= m_atol
1722 || d <= m_rtol
1723 * double((std::max)(math::abs(internal_type(a)), math::abs(internal_type(b)))
1724 );
1725 }
1726
1727 private:
1728
1729 double m_rtol;
1730 double m_atol;
1731 bool m_equal_nan;
1732 };
1733 }
1734
1750 template <class E1, class E2>
1751 inline auto
1752 isclose(E1&& e1, E2&& e2, double rtol = 1e-05, double atol = 1e-08, bool equal_nan = false) noexcept
1753 {
1754 return detail::make_xfunction<detail::isclose>(
1755 std::make_tuple(rtol, atol, equal_nan),
1756 std::forward<E1>(e1),
1757 std::forward<E2>(e2)
1758 );
1759 }
1760
1774 template <class E1, class E2>
1775 inline auto allclose(E1&& e1, E2&& e2, double rtol = 1e-05, double atol = 1e-08) noexcept
1776 {
1777 return xt::all(isclose(std::forward<E1>(e1), std::forward<E2>(e2), rtol, atol));
1778 }
1779
1780 /**********************
1781 * Reducing functions *
1782 **********************/
1783
1787
1803 XTENSOR_REDUCER_FUNCTION(sum, detail::plus, typename std::decay_t<E>::value_type, 0)
1804
1805
1823 XTENSOR_REDUCER_FUNCTION(prod, detail::multiplies, typename std::decay_t<E>::value_type, 1)
1824
1825 namespace detail
1826 {
1827 template <class T, class S, class ST>
1828 inline auto mean_division(S&& s, ST e_size)
1829 {
1830 using value_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
1831 // Avoids floating point exception when s.size is 0
1832 value_type div = s.size() != ST(0) ? static_cast<value_type>(e_size / s.size()) : value_type(0);
1833 return std::move(s) / std::move(div);
1834 }
1835
1836 template <
1837 class T,
1838 class E,
1839 class X,
1840 class D,
1841 class EVS,
1842 XTL_REQUIRES(std::negation<is_reducer_options<X>>, xtl::is_integral<D>)>
1843 inline auto mean(E&& e, X&& axes, const D& ddof, EVS es)
1844 {
1845 // sum cannot always be a double. It could be a complex number which cannot operate on
1846 // std::plus<double>.
1847 using size_type = typename std::decay_t<E>::size_type;
1848 const size_type size = e.size();
1849 XTENSOR_ASSERT(static_cast<size_type>(ddof) <= size);
1850 auto s = sum<T>(std::forward<E>(e), std::forward<X>(axes), es);
1851 return mean_division<T>(std::move(s), size - static_cast<size_type>(ddof));
1852 }
1853
1854 template <class T, class E, class I, std::size_t N, class D, class EVS>
1855 inline auto mean(E&& e, const I (&axes)[N], const D& ddof, EVS es)
1856 {
1857 using size_type = typename std::decay_t<E>::size_type;
1858 const size_type size = e.size();
1859 XTENSOR_ASSERT(static_cast<size_type>(ddof) <= size);
1860 auto s = sum<T>(std::forward<E>(e), axes, es);
1861 return mean_division<T>(std::move(s), size - static_cast<size_type>(ddof));
1862 }
1863
1864 template <class T, class E, class D, class EVS, XTL_REQUIRES(is_reducer_options<EVS>, xtl::is_integral<D>)>
1865 inline auto mean_noaxis(E&& e, const D& ddof, EVS es)
1866 {
1867 using value_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
1868 using size_type = typename std::decay_t<E>::size_type;
1869 const size_type size = e.size();
1870 XTENSOR_ASSERT(static_cast<size_type>(ddof) <= size);
1871 auto s = sum<T>(std::forward<E>(e), es);
1872 return std::move(s) / static_cast<value_type>((size - static_cast<size_type>(ddof)));
1873 }
1874 }
1875
1891 template <
1892 class T = void,
1893 class E,
1894 class X,
1895 class EVS = DEFAULT_STRATEGY_REDUCERS,
1896 XTL_REQUIRES(std::negation<is_reducer_options<X>>)>
1897 inline auto mean(E&& e, X&& axes, EVS es = EVS())
1898 {
1899 return detail::mean<T>(std::forward<E>(e), std::forward<X>(axes), 0u, es);
1900 }
1901
1902 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
1903 inline auto mean(E&& e, EVS es = EVS())
1904 {
1905 return detail::mean_noaxis<T>(std::forward<E>(e), 0u, es);
1906 }
1907
1908 template <class T = void, class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
1909 inline auto mean(E&& e, const I (&axes)[N], EVS es = EVS())
1910 {
1911 return detail::mean<T>(std::forward<E>(e), axes, 0u, es);
1912 }
1913
1931 template <
1932 class T = void,
1933 class E,
1934 class W,
1935 class X,
1936 class EVS = DEFAULT_STRATEGY_REDUCERS,
1937 XTL_REQUIRES(is_reducer_options<EVS>, std::negation<xtl::is_integral<X>>)>
1938 inline auto average(E&& e, W&& weights, X&& axes, EVS ev = EVS())
1939 {
1940 xindex_type_t<typename std::decay_t<E>::shape_type> broadcast_shape;
1941 xt::resize_container(broadcast_shape, e.dimension());
1942 auto ax = normalize_axis(e, axes);
1943 if (weights.dimension() == 1)
1944 {
1945 if (weights.size() != e.shape()[ax[0]])
1946 {
1947 XTENSOR_THROW(std::runtime_error, "Weights need to have the same shape as expression at axes.");
1948 }
1949
1950 std::fill(broadcast_shape.begin(), broadcast_shape.end(), std::size_t(1));
1951 broadcast_shape[ax[0]] = weights.size();
1952 }
1953 else
1954 {
1955 if (!same_shape(e.shape(), weights.shape()))
1956 {
1957 XTENSOR_THROW(
1958 std::runtime_error,
1959 "Weights with dim > 1 need to have the same shape as expression."
1960 );
1961 }
1962
1963 std::copy(e.shape().begin(), e.shape().end(), broadcast_shape.begin());
1964 }
1965
1966 constexpr layout_type L = default_assignable_layout(std::decay_t<W>::static_layout);
1967 auto weights_view = reshape_view<L>(std::forward<W>(weights), std::move(broadcast_shape));
1968 auto scl = sum<T>(weights_view, ax, xt::evaluation_strategy::immediate);
1969 return sum<T>(std::forward<E>(e) * std::move(weights_view), std::move(ax), ev) / std::move(scl);
1970 }
1971
1972 template <
1973 class T = void,
1974 class E,
1975 class W,
1976 class X,
1977 class EVS = DEFAULT_STRATEGY_REDUCERS,
1978 XTL_REQUIRES(is_reducer_options<EVS>, xtl::is_integral<X>)>
1979 inline auto average(E&& e, W&& weights, X axis, EVS ev = EVS())
1980 {
1981 return average(std::forward<E>(e), std::forward<W>(weights), {axis}, std::forward<EVS>(ev));
1982 }
1983
1984 template <class T = void, class E, class W, class X, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
1985 inline auto average(E&& e, W&& weights, const X (&axes)[N], EVS ev = EVS())
1986 {
1987 // need to select the X&& overload and forward to different type
1988 using ax_t = std::array<std::size_t, N>;
1989 return average<T>(std::forward<E>(e), std::forward<W>(weights), xt::forward_normalize<ax_t>(e, axes), ev);
1990 }
1991
1992 template <class T = void, class E, class W, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
1993 inline auto average(E&& e, W&& weights, EVS ev = EVS())
1994 {
1995 if (weights.dimension() != e.dimension()
1996 || !std::equal(weights.shape().begin(), weights.shape().end(), e.shape().begin()))
1997 {
1998 XTENSOR_THROW(std::runtime_error, "Weights need to have the same shape as expression.");
1999 }
2000
2001 auto div = sum<T>(weights, evaluation_strategy::immediate)();
2002 auto s = sum<T>(std::forward<E>(e) * std::forward<W>(weights), ev) / std::move(div);
2003 return s;
2004 }
2005
2006 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2007 inline auto average(E&& e, EVS ev = EVS())
2008 {
2009 return mean<T>(e, ev);
2010 }
2011
2012 namespace detail
2013 {
2014 template <typename E>
2015 std::enable_if_t<std::is_lvalue_reference<E>::value, E> shared_forward(E e) noexcept
2016 {
2017 return e;
2018 }
2019
2020 template <typename E>
2021 std::enable_if_t<!std::is_lvalue_reference<E>::value, xshared_expression<E>> shared_forward(E e) noexcept
2022 {
2023 return make_xshared(std::move(e));
2024 }
2025 }
2026
2027 template <
2028 class T = void,
2029 class E,
2030 class D,
2031 class EVS = DEFAULT_STRATEGY_REDUCERS,
2032 XTL_REQUIRES(is_reducer_options<EVS>, xtl::is_integral<D>)>
2033 inline auto variance(E&& e, const D& ddof, EVS es = EVS())
2034 {
2035 auto cached_mean = mean<T>(e, es)();
2036 return detail::mean_noaxis<T>(square(std::forward<E>(e) - std::move(cached_mean)), ddof, es);
2037 }
2038
2039 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2040 inline auto variance(E&& e, EVS es = EVS())
2041 {
2042 return variance<T>(std::forward<E>(e), 0u, es);
2043 }
2044
2045 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2046 inline auto stddev(E&& e, EVS es = EVS())
2047 {
2048 return sqrt(variance<T>(std::forward<E>(e), es));
2049 }
2050
2075 template <
2076 class T = void,
2077 class E,
2078 class X,
2079 class D,
2080 class EVS = DEFAULT_STRATEGY_REDUCERS,
2081 XTL_REQUIRES(std::negation<is_reducer_options<X>>, xtl::is_integral<D>)>
2082 inline auto variance(E&& e, X&& axes, const D& ddof, EVS es = EVS())
2083 {
2084 decltype(auto) sc = detail::shared_forward<E>(e);
2085 // note: forcing copy of first axes argument -- is there a better solution?
2086 auto axes_copy = axes;
2087 // always eval to prevent repeated evaluations in the next calls
2088 auto inner_mean = eval(mean<T>(sc, std::move(axes_copy), evaluation_strategy::immediate));
2089
2090 // fake keep_dims = 1
2091 // Since the inner_shape might have a reference semantic (e.g. xbuffer_adaptor in bindings)
2092 // We need to map it to another type before modifying it.
2093 // We pragmatically abuse `get_strides_t`
2094 using tmp_shape_t = get_strides_t<typename std::decay_t<E>::shape_type>;
2095 tmp_shape_t keep_dim_shape = xtl::forward_sequence<tmp_shape_t, decltype(e.shape())>(e.shape());
2096 for (const auto& el : axes)
2097 {
2098 keep_dim_shape[el] = 1u;
2099 }
2100
2101 auto mrv = reshape_view<XTENSOR_DEFAULT_LAYOUT>(std::move(inner_mean), std::move(keep_dim_shape));
2102 return detail::mean<T>(square(sc - std::move(mrv)), std::forward<X>(axes), ddof, es);
2103 }
2104
2105 template <
2106 class T = void,
2107 class E,
2108 class X,
2109 class EVS = DEFAULT_STRATEGY_REDUCERS,
2110 XTL_REQUIRES(std::negation<is_reducer_options<X>>, std::negation<xtl::is_integral<std::decay_t<X>>>, is_reducer_options<EVS>)>
2111 inline auto variance(E&& e, X&& axes, EVS es = EVS())
2112 {
2113 return variance<T>(std::forward<E>(e), std::forward<X>(axes), 0u, es);
2114 }
2115
2137 template <
2138 class T = void,
2139 class E,
2140 class X,
2141 class EVS = DEFAULT_STRATEGY_REDUCERS,
2142 XTL_REQUIRES(std::negation<is_reducer_options<X>>)>
2143 inline auto stddev(E&& e, X&& axes, EVS es = EVS())
2144 {
2145 return sqrt(variance<T>(std::forward<E>(e), std::forward<X>(axes), es));
2146 }
2147
2148 template <class T = void, class E, class A, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2149 inline auto stddev(E&& e, const A (&axes)[N], EVS es = EVS())
2150 {
2151 return stddev<T>(
2152 std::forward<E>(e),
2153 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2154 es
2155 );
2156 }
2157
2158 template <
2159 class T = void,
2160 class E,
2161 class A,
2162 std::size_t N,
2163 class EVS = DEFAULT_STRATEGY_REDUCERS,
2164 XTL_REQUIRES(is_reducer_options<EVS>)>
2165 inline auto variance(E&& e, const A (&axes)[N], EVS es = EVS())
2166 {
2167 return variance<T>(
2168 std::forward<E>(e),
2169 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2170 es
2171 );
2172 }
2173
2174 template <class T = void, class E, class A, std::size_t N, class D, class EVS = DEFAULT_STRATEGY_REDUCERS>
2175 inline auto variance(E&& e, const A (&axes)[N], const D& ddof, EVS es = EVS())
2176 {
2177 return variance<T>(
2178 std::forward<E>(e),
2179 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2180 ddof,
2181 es
2182 );
2183 }
2184
2195 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2196 inline auto minmax(E&& e, EVS es = EVS())
2197 {
2198 using std::max;
2199 using std::min;
2200 using value_type = typename std::decay_t<E>::value_type;
2201 using result_type = std::array<value_type, 2>;
2202 using init_value_fct = xt::const_value<result_type>;
2203
2204 auto reduce_func = [](auto r, const auto& v)
2205 {
2206 r[0] = (min) (r[0], v);
2207 r[1] = (max) (r[1], v);
2208 return r;
2209 };
2210
2211 auto init_func = init_value_fct(
2212 result_type{std::numeric_limits<value_type>::max(), std::numeric_limits<value_type>::lowest()}
2213 );
2214
2215 auto merge_func = [](auto r, const auto& s)
2216 {
2217 r[0] = (min) (r[0], s[0]);
2218 r[1] = (max) (r[1], s[1]);
2219 return r;
2220 };
2221 return xt::reduce(
2222 make_xreducer_functor(std::move(reduce_func), std::move(init_func), std::move(merge_func)),
2223 std::forward<E>(e),
2224 arange(e.dimension()),
2225 es
2226 );
2227 }
2228
2232
2247 template <class T = void, class E>
2248 inline auto cumsum(E&& e, std::ptrdiff_t axis)
2249 {
2250 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2251 return accumulate(
2252 make_xaccumulator_functor(detail::plus(), detail::accumulator_identity<init_value_type>()),
2253 std::forward<E>(e),
2254 axis
2255 );
2256 }
2257
2258 template <class T = void, class E>
2259 inline auto cumsum(E&& e)
2260 {
2261 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2262 return accumulate(
2263 make_xaccumulator_functor(detail::plus(), detail::accumulator_identity<init_value_type>()),
2264 std::forward<E>(e)
2265 );
2266 }
2267
2282 template <class T = void, class E>
2283 inline auto cumprod(E&& e, std::ptrdiff_t axis)
2284 {
2285 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2286 return accumulate(
2287 make_xaccumulator_functor(detail::multiplies(), detail::accumulator_identity<init_value_type>()),
2288 std::forward<E>(e),
2289 axis
2290 );
2291 }
2292
2293 template <class T = void, class E>
2294 inline auto cumprod(E&& e)
2295 {
2296 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2297 return accumulate(
2298 make_xaccumulator_functor(detail::multiplies(), detail::accumulator_identity<init_value_type>()),
2299 std::forward<E>(e)
2300 );
2301 }
2302
2303 /*****************
2304 * nan functions *
2305 *****************/
2306
2307 namespace detail
2308 {
2309 struct nan_to_num_functor
2310 {
2311 template <class A>
2312 inline auto operator()(const A& a) const
2313 {
2314 if (math::isnan(a))
2315 {
2316 return A(0);
2317 }
2318 if (math::isinf(a))
2319 {
2320 if (a < 0)
2321 {
2322 return std::numeric_limits<A>::lowest();
2323 }
2324 else
2325 {
2326 return (std::numeric_limits<A>::max)();
2327 }
2328 }
2329 return a;
2330 }
2331 };
2332
2333 struct nan_min
2334 {
2335 template <class T, class U>
2336 constexpr auto operator()(const T lhs, const U rhs) const
2337 {
2338 // Clunky expression for working with GCC 4.9
2339 return math::isnan(lhs)
2340 ? rhs
2341 : (math::isnan(rhs) ? lhs
2342 : std::common_type_t<T, U>(
2343 detail::make_xfunction<math::minimum<void>>(lhs, rhs)
2344 ));
2345 }
2346 };
2347
2348 struct nan_max
2349 {
2350 template <class T, class U>
2351 constexpr auto operator()(const T lhs, const U rhs) const
2352 {
2353 // Clunky expression for working with GCC 4.9
2354 return math::isnan(lhs)
2355 ? rhs
2356 : (math::isnan(rhs) ? lhs
2357 : std::common_type_t<T, U>(
2358 detail::make_xfunction<math::maximum<void>>(lhs, rhs)
2359 ));
2360 }
2361 };
2362
2363 struct nan_plus
2364 {
2365 template <class T, class U>
2366 constexpr auto operator()(const T lhs, const U rhs) const
2367 {
2368 return !math::isnan(rhs) ? lhs + rhs : lhs;
2369 }
2370 };
2371
2372 struct nan_multiplies
2373 {
2374 template <class T, class U>
2375 constexpr auto operator()(const T lhs, const U rhs) const
2376 {
2377 return !math::isnan(rhs) ? lhs * rhs : lhs;
2378 }
2379 };
2380
2381 template <class T, int V>
2382 struct nan_init
2383 {
2384 using value_type = T;
2385 using result_type = T;
2386
2387 constexpr result_type operator()(const value_type lhs) const
2388 {
2389 return math::isnan(lhs) ? result_type(V) : lhs;
2390 }
2391 };
2392 }
2393
2397
2408 template <class E>
2409 inline auto nan_to_num(E&& e)
2410 {
2411 return detail::make_xfunction<detail::nan_to_num_functor>(std::forward<E>(e));
2412 }
2413
2427 XTENSOR_REDUCER_FUNCTION(nanmin, detail::nan_min, typename std::decay_t<E>::value_type, std::nan("0"))
2428
2429
2442 XTENSOR_REDUCER_FUNCTION(nanmax, detail::nan_max, typename std::decay_t<E>::value_type, std::nan("0"))
2443
2459 XTENSOR_REDUCER_FUNCTION(nansum, detail::nan_plus, typename std::decay_t<E>::value_type, 0)
2460
2476 XTENSOR_REDUCER_FUNCTION(nanprod, detail::nan_multiplies, typename std::decay_t<E>::value_type, 1)
2477
2478#define COUNT_NON_ZEROS_CONTENT \
2479 using value_type = typename std::decay_t<E>::value_type; \
2480 using result_type = xt::detail::xreducer_size_type_t<value_type>; \
2481 using init_value_fct = xt::const_value<result_type>; \
2482 \
2483 auto init_fct = init_value_fct(0); \
2484 \
2485 auto reduce_fct = [](const auto& lhs, const auto& rhs) \
2486 { \
2487 using value_t = xt::detail::xreducer_temporary_type_t<std::decay_t<decltype(rhs)>>; \
2488 using result_t = std::decay_t<decltype(lhs)>; \
2489 \
2490 return (rhs != value_t(0)) ? lhs + result_t(1) : lhs; \
2491 }; \
2492 auto merge_func = detail::plus();
2493
2494 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2495 inline auto count_nonzero(E&& e, EVS es = EVS())
2496 {
2497 COUNT_NON_ZEROS_CONTENT;
2498 return xt::reduce(
2499 make_xreducer_functor(std::move(reduce_fct), std::move(init_fct), std::move(merge_func)),
2500 std::forward<E>(e),
2501 es
2502 );
2503 }
2504
2505 template <
2506 class E,
2507 class X,
2508 class EVS = DEFAULT_STRATEGY_REDUCERS,
2509 XTL_REQUIRES(std::negation<is_reducer_options<X>>, std::negation<xtl::is_integral<X>>)>
2510 inline auto count_nonzero(E&& e, X&& axes, EVS es = EVS())
2511 {
2512 COUNT_NON_ZEROS_CONTENT;
2513 return xt::reduce(
2514 make_xreducer_functor(std::move(reduce_fct), std::move(init_fct), std::move(merge_func)),
2515 std::forward<E>(e),
2516 std::forward<X>(axes),
2517 es
2518 );
2519 }
2520
2521 template <
2522 class E,
2523 class X,
2524 class EVS = DEFAULT_STRATEGY_REDUCERS,
2525 XTL_REQUIRES(std::negation<is_reducer_options<X>>, xtl::is_integral<X>)>
2526 inline auto count_nonzero(E&& e, X axis, EVS es = EVS())
2527 {
2528 return count_nonzero(std::forward<E>(e), {axis}, es);
2529 }
2530
2531 template <class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2532 inline auto count_nonzero(E&& e, const I (&axes)[N], EVS es = EVS())
2533 {
2534 COUNT_NON_ZEROS_CONTENT;
2535 return xt::reduce(
2536 make_xreducer_functor(std::move(reduce_fct), std::move(init_fct), std::move(merge_func)),
2537 std::forward<E>(e),
2538 axes,
2539 es
2540 );
2541 }
2542
2543#undef COUNT_NON_ZEROS_CONTENT
2544
2545 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2546 inline auto count_nonnan(E&& e, EVS es = EVS())
2547 {
2548 return xt::count_nonzero(!xt::isnan(std::forward<E>(e)), es);
2549 }
2550
2551 template <
2552 class E,
2553 class X,
2554 class EVS = DEFAULT_STRATEGY_REDUCERS,
2555 XTL_REQUIRES(std::negation<is_reducer_options<X>>, std::negation<xtl::is_integral<X>>)>
2556 inline auto count_nonnan(E&& e, X&& axes, EVS es = EVS())
2557 {
2558 return xt::count_nonzero(!xt::isnan(std::forward<E>(e)), std::forward<X>(axes), es);
2559 }
2560
2561 template <
2562 class E,
2563 class X,
2564 class EVS = DEFAULT_STRATEGY_REDUCERS,
2565 XTL_REQUIRES(std::negation<is_reducer_options<X>>, xtl::is_integral<X>)>
2566 inline auto count_nonnan(E&& e, X&& axes, EVS es = EVS())
2567 {
2568 return xt::count_nonzero(!xt::isnan(std::forward<E>(e)), {axes}, es);
2569 }
2570
2571 template <class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2572 inline auto count_nonnan(E&& e, const I (&axes)[N], EVS es = EVS())
2573 {
2574 return xt::count_nonzero(!xt::isnan(std::forward<E>(e)), axes, es);
2575 }
2576
2591 template <class T = void, class E>
2592 inline auto nancumsum(E&& e, std::ptrdiff_t axis)
2593 {
2594 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2595 return accumulate(
2596 make_xaccumulator_functor(detail::nan_plus(), detail::nan_init<init_value_type, 0>()),
2597 std::forward<E>(e),
2598 axis
2599 );
2600 }
2601
2602 template <class T = void, class E>
2603 inline auto nancumsum(E&& e)
2604 {
2605 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2606 return accumulate(
2607 make_xaccumulator_functor(detail::nan_plus(), detail::nan_init<init_value_type, 0>()),
2608 std::forward<E>(e)
2609 );
2610 }
2611
2626 template <class T = void, class E>
2627 inline auto nancumprod(E&& e, std::ptrdiff_t axis)
2628 {
2629 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2630 return accumulate(
2631 make_xaccumulator_functor(detail::nan_multiplies(), detail::nan_init<init_value_type, 1>()),
2632 std::forward<E>(e),
2633 axis
2634 );
2635 }
2636
2637 template <class T = void, class E>
2638 inline auto nancumprod(E&& e)
2639 {
2640 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2641 return accumulate(
2642 make_xaccumulator_functor(detail::nan_multiplies(), detail::nan_init<init_value_type, 1>()),
2643 std::forward<E>(e)
2644 );
2645 }
2646
2647 namespace detail
2648 {
2649 template <class T>
2650 struct diff_impl
2651 {
2652 template <class Arg>
2653 inline void operator()(
2654 Arg& ad,
2655 const std::size_t& n,
2656 xstrided_slice_vector& slice1,
2657 xstrided_slice_vector& slice2,
2658 std::size_t saxis
2659 )
2660 {
2661 for (std::size_t i = 0; i < n; ++i)
2662 {
2663 slice2[saxis] = range(xnone(), ad.shape()[saxis] - 1);
2664 ad = strided_view(ad, slice1) - strided_view(ad, slice2);
2665 }
2666 }
2667 };
2668
2669 template <>
2670 struct diff_impl<bool>
2671 {
2672 template <class Arg>
2673 inline void operator()(
2674 Arg& ad,
2675 const std::size_t& n,
2676 xstrided_slice_vector& slice1,
2677 xstrided_slice_vector& slice2,
2678 std::size_t saxis
2679 )
2680 {
2681 for (std::size_t i = 0; i < n; ++i)
2682 {
2683 slice2[saxis] = range(xnone(), ad.shape()[saxis] - 1);
2684 ad = not_equal(strided_view(ad, slice1), strided_view(ad, slice2));
2685 }
2686 }
2687 };
2688 }
2689
2705 template <
2706 class T = void,
2707 class E,
2708 class X,
2709 class EVS = DEFAULT_STRATEGY_REDUCERS,
2710 XTL_REQUIRES(std::negation<is_reducer_options<X>>)>
2711 inline auto nanmean(E&& e, X&& axes, EVS es = EVS())
2712 {
2713 decltype(auto) sc = detail::shared_forward<E>(e);
2714 // note: forcing copy of first axes argument -- is there a better solution?
2715 auto axes_copy = axes;
2716 using value_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
2717 using sum_type = typename std::conditional_t<
2718 std::is_same<T, void>::value,
2719 typename std::common_type_t<typename std::decay_t<E>::value_type, value_type>,
2720 T>;
2721 // sum cannot always be a double. It could be a complex number which cannot operate on
2722 // std::plus<double>.
2723 return nansum<sum_type>(sc, std::forward<X>(axes), es)
2724 / xt::cast<value_type>(count_nonnan(sc, std::move(axes_copy), es));
2725 }
2726
2727 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2728 inline auto nanmean(E&& e, EVS es = EVS())
2729 {
2730 decltype(auto) sc = detail::shared_forward<E>(e);
2731 using value_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
2732 using sum_type = typename std::conditional_t<
2733 std::is_same<T, void>::value,
2734 typename std::common_type_t<typename std::decay_t<E>::value_type, value_type>,
2735 T>;
2736 return nansum<sum_type>(sc, es) / xt::cast<value_type>(count_nonnan(sc, es));
2737 }
2738
2739 template <class T = void, class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2740 inline auto nanmean(E&& e, const I (&axes)[N], EVS es = EVS())
2741 {
2742 return nanmean<T>(
2743 std::forward<E>(e),
2744 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2745 es
2746 );
2747 }
2748
2749 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2750 inline auto nanvar(E&& e, EVS es = EVS())
2751 {
2752 decltype(auto) sc = detail::shared_forward<E>(e);
2753 return nanmean<T>(square(sc - nanmean<T>(sc)), es);
2754 }
2755
2756 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2757 inline auto nanstd(E&& e, EVS es = EVS())
2758 {
2759 return sqrt(nanvar<T>(std::forward<E>(e), es));
2760 }
2761
2782 template <
2783 class T = void,
2784 class E,
2785 class X,
2786 class EVS = DEFAULT_STRATEGY_REDUCERS,
2787 XTL_REQUIRES(std::negation<is_reducer_options<X>>)>
2788 inline auto nanvar(E&& e, X&& axes, EVS es = EVS())
2789 {
2790 decltype(auto) sc = detail::shared_forward<E>(e);
2791 // note: forcing copy of first axes argument -- is there a better solution?
2792 auto axes_copy = axes;
2793 using result_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
2794 auto inner_mean = nanmean<result_type>(sc, std::move(axes_copy));
2795
2796 // fake keep_dims = 1
2797 // Since the inner_shape might have a reference semantic (e.g. xbuffer_adaptor in bindings)
2798 // We need to map it to another type before modifying it.
2799 // We pragmatically abuse `get_strides_t`
2800 using tmp_shape_t = get_strides_t<typename std::decay_t<E>::shape_type>;
2801 tmp_shape_t keep_dim_shape = xtl::forward_sequence<tmp_shape_t, decltype(e.shape())>(e.shape());
2802 for (const auto& el : axes)
2803 {
2804 keep_dim_shape[el] = 1;
2805 }
2806 auto mrv = reshape_view<XTENSOR_DEFAULT_LAYOUT>(std::move(inner_mean), std::move(keep_dim_shape));
2807 return nanmean<result_type>(square(cast<result_type>(sc) - std::move(mrv)), std::forward<X>(axes), es);
2808 }
2809
2830 template <
2831 class T = void,
2832 class E,
2833 class X,
2834 class EVS = DEFAULT_STRATEGY_REDUCERS,
2835 XTL_REQUIRES(std::negation<is_reducer_options<X>>)>
2836 inline auto nanstd(E&& e, X&& axes, EVS es = EVS())
2837 {
2838 return sqrt(nanvar<T>(std::forward<E>(e), std::forward<X>(axes), es));
2839 }
2840
2841 template <class T = void, class E, class A, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2842 inline auto nanstd(E&& e, const A (&axes)[N], EVS es = EVS())
2843 {
2844 return nanstd<T>(
2845 std::forward<E>(e),
2846 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2847 es
2848 );
2849 }
2850
2851 template <class T = void, class E, class A, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2852 inline auto nanvar(E&& e, const A (&axes)[N], EVS es = EVS())
2853 {
2854 return nanvar<T>(
2855 std::forward<E>(e),
2856 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2857 es
2858 );
2859 }
2860
2872 template <class T>
2873 auto diff(const xexpression<T>& a, std::size_t n = 1, std::ptrdiff_t axis = -1)
2874 {
2875 typename std::decay_t<T>::temporary_type ad = a.derived_cast();
2876 std::size_t saxis = normalize_axis(ad.dimension(), axis);
2877 if (n <= ad.size())
2878 {
2879 if (n != std::size_t(0))
2880 {
2881 xstrided_slice_vector slice1(ad.dimension(), all());
2882 xstrided_slice_vector slice2(ad.dimension(), all());
2883 slice1[saxis] = range(1, xnone());
2884
2885 detail::diff_impl<typename T::value_type> impl;
2886 impl(ad, n, slice1, slice2, saxis);
2887 }
2888 }
2889 else
2890 {
2891 auto shape = ad.shape();
2892 shape[saxis] = std::size_t(0);
2893 ad.resize(shape);
2894 }
2895 return ad;
2896 }
2897
2909 template <class T>
2910 auto trapz(const xexpression<T>& y, double dx = 1.0, std::ptrdiff_t axis = -1)
2911 {
2912 auto& yd = y.derived_cast();
2913 std::size_t saxis = normalize_axis(yd.dimension(), axis);
2914
2915 xstrided_slice_vector slice1(yd.dimension(), all());
2916 xstrided_slice_vector slice2(yd.dimension(), all());
2917 slice1[saxis] = range(1, xnone());
2918 slice2[saxis] = range(xnone(), yd.shape()[saxis] - 1);
2919
2920 auto trap = dx * (strided_view(yd, slice1) + strided_view(yd, slice2)) * 0.5;
2921
2922 return eval(sum(trap, {saxis}));
2923 }
2924
2936 template <class T, class E>
2937 auto trapz(const xexpression<T>& y, const xexpression<E>& x, std::ptrdiff_t axis = -1)
2938 {
2939 auto& yd = y.derived_cast();
2940 auto& xd = x.derived_cast();
2941 decltype(diff(x)) dx;
2942
2943 std::size_t saxis = normalize_axis(yd.dimension(), axis);
2944
2945 if (xd.dimension() == 1)
2946 {
2947 dx = diff(x);
2948 typename std::decay_t<decltype(yd)>::shape_type shape;
2949 resize_container(shape, yd.dimension());
2950 std::fill(shape.begin(), shape.end(), 1);
2951 shape[saxis] = dx.shape()[0];
2952 dx.reshape(shape);
2953 }
2954 else
2955 {
2956 dx = diff(x, 1, axis);
2957 }
2958
2959 xstrided_slice_vector slice1(yd.dimension(), all());
2960 xstrided_slice_vector slice2(yd.dimension(), all());
2961 slice1[saxis] = range(1, xnone());
2962 slice2[saxis] = range(xnone(), yd.shape()[saxis] - 1);
2963
2964 auto trap = dx * (strided_view(yd, slice1) + strided_view(yd, slice2)) * 0.5;
2965
2966 return eval(sum(trap, {saxis}));
2967 }
2968
2981 template <class E1, class E2, class E3, typename T>
2982 inline auto interp(const E1& x, const E2& xp, const E3& fp, T left, T right)
2983 {
2984 using size_type = common_size_type_t<E1, E2, E3>;
2985 using value_type = typename E3::value_type;
2986
2987 // basic checks
2988 XTENSOR_ASSERT(xp.dimension() == 1);
2989 XTENSOR_ASSERT(std::is_sorted(x.cbegin(), x.cend()));
2990 XTENSOR_ASSERT(std::is_sorted(xp.cbegin(), xp.cend()));
2991
2992 // allocate output
2993 auto f = xtensor<value_type, 1>::from_shape(x.shape());
2994
2995 // counter in "x": from left
2996 size_type i = 0;
2997
2998 // fill f[i] for x[i] <= xp[0]
2999 for (; i < x.size(); ++i)
3000 {
3001 if (x[i] > xp[0])
3002 {
3003 break;
3004 }
3005 f[i] = static_cast<value_type>(left);
3006 }
3007
3008 // counter in "x": from right
3009 // (index counts one right, to terminate the reverse loop, without risking being negative)
3010 size_type imax = x.size();
3011
3012 // fill f[i] for x[-1] >= xp[-1]
3013 for (; imax > 0; --imax)
3014 {
3015 if (x[imax - 1] < xp[xp.size() - 1])
3016 {
3017 break;
3018 }
3019 f[imax - 1] = static_cast<value_type>(right);
3020 }
3021
3022 // catch edge case: all entries are "right"
3023 if (imax == 0)
3024 {
3025 return f;
3026 }
3027
3028 // set "imax" as actual index
3029 // (counted one right, see above)
3030 --imax;
3031
3032 // counter in "xp"
3033 size_type ip = 1;
3034
3035 // fill f[i] for the interior
3036 for (; i <= imax; ++i)
3037 {
3038 // - search next value in "xp"
3039 while (x[i] > xp[ip])
3040 {
3041 ++ip;
3042 }
3043 // - distances as doubles
3044 double dfp = static_cast<double>(fp[ip] - fp[ip - 1]);
3045 double dxp = static_cast<double>(xp[ip] - xp[ip - 1]);
3046 double dx = static_cast<double>(x[i] - xp[ip - 1]);
3047 // - interpolate
3048 f[i] = fp[ip - 1] + static_cast<value_type>(dfp / dxp * dx);
3049 }
3050
3051 return f;
3052 }
3053
3054 namespace detail
3055 {
3056 template <class E1, class E2>
3057 auto calculate_discontinuity(E1&& discontinuity, E2&&)
3058 {
3059 return discontinuity;
3060 }
3061
3062 template <class E2>
3063 auto calculate_discontinuity(xt::placeholders::xtuph, E2&& period)
3064 {
3065 return 0.5 * period;
3066 }
3067
3068 template <class E1, class E2>
3069 auto
3070 calculate_interval(E2&& period, typename std::enable_if<std::is_integral<E1>::value, E1>::type* = 0)
3071 {
3072 auto interval_high = 0.5 * period;
3073 uint64_t remainder = static_cast<uint64_t>(period) % 2;
3074 auto boundary_ambiguous = (remainder == 0);
3075 return std::make_tuple(interval_high, boundary_ambiguous);
3076 }
3077
3078 template <class E1, class E2>
3079 auto
3080 calculate_interval(E2&& period, typename std::enable_if<std::is_floating_point<E1>::value, E1>::type* = 0)
3081 {
3082 auto interval_high = 0.5 * period;
3083 auto boundary_ambiguous = true;
3084 return std::make_tuple(interval_high, boundary_ambiguous);
3085 }
3086 }
3087
3100
3101 template <class E1, class E2 = xt::placeholders::xtuph, class E3 = double>
3102 inline auto unwrap(
3103 E1&& p,
3104 E2 discontinuity = xnone(),
3105 std::ptrdiff_t axis = -1,
3106 E3 period = 2.0 * xt::numeric_constants<double>::PI
3107 )
3108 {
3109 auto discont = detail::calculate_discontinuity(discontinuity, period);
3110 using value_type = typename std::decay_t<E1>::value_type;
3111 std::size_t saxis = normalize_axis(p.dimension(), axis);
3112 auto dd = diff(p, 1, axis);
3113 xstrided_slice_vector slice(p.dimension(), all());
3114 slice[saxis] = range(1, xnone());
3115 auto interval_tuple = detail::calculate_interval<value_type>(period);
3116 auto interval_high = std::get<0>(interval_tuple);
3117 auto boundary_ambiguous = std::get<1>(interval_tuple);
3118 auto interval_low = -interval_high;
3119 auto ddmod = xt::eval(xt::fmod(xt::fmod(dd - interval_low, period) + period, period) + interval_low);
3120 if (boundary_ambiguous)
3121 {
3122 // for `mask = (abs(dd) == period/2)`, the above line made
3123 //`ddmod[mask] == -period/2`. correct these such that
3124 //`ddmod[mask] == sign(dd[mask])*period/2`.
3125 auto boolmap = xt::equal(ddmod, interval_low) && (xt::greater(dd, 0.0));
3126 ddmod = xt::where(boolmap, interval_high, ddmod);
3127 }
3128 auto ph_correct = xt::eval(ddmod - dd);
3129 ph_correct = xt::where(xt::abs(dd) < discont, 0.0, ph_correct);
3130 E1 up(p);
3131 strided_view(up, slice) = strided_view(p, slice)
3132 + xt::cumsum(ph_correct, static_cast<std::ptrdiff_t>(saxis));
3133 return up;
3134 }
3135
3146 template <class E1, class E2, class E3>
3147 inline auto interp(const E1& x, const E2& xp, const E3& fp)
3148 {
3149 return interp(x, xp, fp, fp[0], fp[fp.size() - 1]);
3150 }
3151
3158 template <class E1>
3159 inline auto cov(const E1& x, const E1& y = E1())
3160 {
3161 using value_type = typename E1::value_type;
3162
3163 if (y.dimension() == 0)
3164 {
3165 auto s = x.shape();
3166 using size_type = std::decay_t<decltype(s[0])>;
3167 if (x.dimension() == 1)
3168 {
3169 auto covar = eval(zeros<value_type>({1, 1}));
3170 auto x_norm = x - eval(mean(x));
3171 covar(0, 0) = std::inner_product(x_norm.begin(), x_norm.end(), x_norm.begin(), 0.0)
3172 / value_type(s[0] - 1);
3173 return covar;
3174 }
3175
3176 XTENSOR_ASSERT(x.dimension() == 2);
3177
3178 auto covar = eval(zeros<value_type>({s[0], s[0]}));
3179 auto m = eval(mean(x, {1}));
3180 m.reshape({m.shape()[0], 1});
3181 auto x_norm = x - m;
3182 for (size_type i = 0; i < s[0]; i++)
3183 {
3184 auto xi = strided_view(x_norm, {range(i, i + 1), all()});
3185 for (size_type j = i; j < s[0]; j++)
3186 {
3187 auto xj = strided_view(x_norm, {range(j, j + 1), all()});
3188 covar(j, i) = std::inner_product(xi.begin(), xi.end(), xj.begin(), 0.0)
3189 / value_type(s[1] - 1);
3190 }
3191 }
3192 return eval(covar + transpose(covar) - diag(diagonal(covar)));
3193 }
3194 else
3195 {
3196 return cov(eval(stack(xtuple(x, y))));
3197 }
3198 }
3199
3200 /*
3201 * convolution mode placeholders for selecting the algorithm
3202 * used in computing a 1D convolution.
3203 * Same as NumPy's mode parameter.
3204 */
3205 namespace convolve_mode
3206 {
3207 struct valid
3208 {
3209 };
3210
3211 struct full
3212 {
3213 };
3214 }
3215
3216 namespace detail
3217 {
3218 template <class E1, class E2>
3219 inline auto convolve_impl(E1&& e1, E2&& e2, convolve_mode::valid)
3220 {
3221 using value_type = typename std::decay<E1>::type::value_type;
3222
3223 const std::size_t na = e1.size();
3224 const std::size_t nv = e2.size();
3225 const std::size_t n = na - nv + 1;
3226 xt::xtensor<value_type, 1> out = xt::zeros<value_type>({n});
3227 for (std::size_t i = 0; i < n; i++)
3228 {
3229 for (std::size_t j = 0; j < nv; j++)
3230 {
3231 out(i) += e1(j) * e2(j + i);
3232 }
3233 }
3234 return out;
3235 }
3236
3237 template <class E1, class E2>
3238 inline auto convolve_impl(E1&& e1, E2&& e2, convolve_mode::full)
3239 {
3240 using value_type = typename std::decay<E1>::type::value_type;
3241
3242 const std::size_t na = e1.size();
3243 const std::size_t nv = e2.size();
3244 const std::size_t n = na + nv - 1;
3245 xt::xtensor<value_type, 1> out = xt::zeros<value_type>({n});
3246 for (std::size_t i = 0; i < n; i++)
3247 {
3248 const std::size_t jmn = (i >= nv - 1) ? i - (nv - 1) : 0;
3249 const std::size_t jmx = (i < na - 1) ? i : na - 1;
3250 for (std::size_t j = jmn; j <= jmx; ++j)
3251 {
3252 out(i) += e1(j) * e2(i - j);
3253 }
3254 }
3255 return out;
3256 }
3257 }
3258
3259 /*
3260 * @brief computes the 1D convolution between two 1D expressions
3261 *
3262 * @param a 1D expression
3263 * @param v 1D expression
3264 * @param mode placeholder Select algorithm #convolve_mode
3265 *
3266 * @detail the algorithm convolves a with v and will incur a copy overhead
3267 * should v be longer than a.
3268 */
3269 template <class E1, class E2, class E3>
3270 inline auto convolve(E1&& a, E2&& v, E3 mode)
3271 {
3272 if (a.dimension() != 1 || v.dimension() != 1)
3273 {
3274 XTENSOR_THROW(std::runtime_error, "Invalid dimentions convolution arguments must be 1D expressions");
3275 }
3276
3277 XTENSOR_ASSERT(a.size() > 0 && v.size() > 0);
3278
3279 // swap them so a is always the longest one
3280 if (a.size() < v.size())
3281 {
3282 return detail::convolve_impl(std::forward<E2>(v), std::forward<E1>(a), mode);
3283 }
3284 else
3285 {
3286 return detail::convolve_impl(std::forward<E1>(a), std::forward<E2>(v), mode);
3287 }
3288 }
3289}
3290
3291
3292#endif
xt::xexpression
Base class for xexpressions.
Definition xexpression.hpp:45
xt::xexpression::derived_cast
derived_type & derived_cast() &noexcept
Returns a reference to the actual derived type of the xexpression.
Definition xexpression.hpp:115
xt::cumprod
auto cumprod(E &&e, std::ptrdiff_t axis)
Cumulative product.
Definition xmath.hpp:2283
xt::cumsum
auto cumsum(E &&e, std::ptrdiff_t axis)
Cumulative sum.
Definition xmath.hpp:2248
xt::fma
auto fma(E1 &&e1, E2 &&e2, E3 &&e3) noexcept -> detail::xfunction_type_t< math::fma_fun, E1, E2, E3 >
Fused multiply-add operation.
Definition xmath.hpp:510
xt::deg2rad
auto deg2rad(E &&e) noexcept -> detail::xfunction_type_t< math::deg2rad, E >
Convert angles from degrees to radians.
Definition xmath.hpp:684
xt::amax
auto amax(E &&e, X &&axes, EVS es=EVS())
Maximum element along given axis.
Definition xmath.hpp:782
xt::remainder
auto remainder(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::remainder_fun, E1, E2 >
Signed remainder of the division operation.
Definition xmath.hpp:492
xt::degrees
auto degrees(E &&e) noexcept -> detail::xfunction_type_t< math::rad2deg, E >
Convert angles from radians to degrees.
Definition xmath.hpp:729
xt::interp
auto interp(const E1 &x, const E2 &xp, const E3 &fp, T left, T right)
Returns the one-dimensional piecewise linear interpolant to a function with given discrete data point...
Definition xmath.hpp:2982
xt::fmod
auto fmod(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::fmod_fun, E1, E2 >
Remainder of the floating point division operation.
Definition xmath.hpp:475
xt::abs
auto abs(E &&e) noexcept -> detail::xfunction_type_t< math::abs_fun, E >
Absolute value function.
Definition xmath.hpp:443
xt::fabs
auto fabs(E &&e) noexcept -> detail::xfunction_type_t< math::fabs_fun, E >
Absolute value function.
Definition xmath.hpp:458
xt::minimum
auto minimum(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::minimum< void >, E1, E2 >
Elementwise minimum.
Definition xmath.hpp:761
xt::maximum
auto maximum(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::maximum< void >, E1, E2 >
Elementwise maximum.
Definition xmath.hpp:745
xt::fmax
auto fmax(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::fmax_fun, E1, E2 >
Maximum function.
Definition xmath.hpp:531
xt::clip
auto clip(E1 &&e1, E2 &&lo, E3 &&hi) noexcept -> detail::xfunction_type_t< math::clamp_fun, E1, E2, E3 >
Clip values between hi and lo.
Definition xmath.hpp:815
xt::radians
auto radians(E &&e) noexcept -> detail::xfunction_type_t< math::deg2rad, E >
Convert angles from degrees to radians.
Definition xmath.hpp:699
xt::fdim
auto fdim(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::fdim_fun, E1, E2 >
Positive difference function.
Definition xmath.hpp:565
xt::amin
auto amin(E &&e, X &&axes, EVS es=EVS())
Minimum element along given axis.
Definition xmath.hpp:800
xt::rad2deg
auto rad2deg(E &&e) noexcept -> detail::xfunction_type_t< math::rad2deg, E >
Convert angles from radians to degrees.
Definition xmath.hpp:714
xt::fmin
auto fmin(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::fmin_fun, E1, E2 >
Minimum function.
Definition xmath.hpp:548
xt::sign
auto sign(E &&e) noexcept -> detail::xfunction_type_t< math::sign_fun, E >
Returns an element-wise indication of the sign of a number.
Definition xmath.hpp:877
xt::unwrap
auto unwrap(E1 &&p, E2 discontinuity=xnone(), std::ptrdiff_t axis=-1, E3 period=2.0 *xt::numeric_constants< double >::PI)
Unwrap by taking the complement of large deltas with respect to the period.
Definition xmath.hpp:3102
xt::cast
auto cast(E &&e) noexcept -> detail::xfunction_type_t< typename detail::cast< R >::functor, E >
Element-wise static_cast.
Definition xoperation.hpp:989
xt::allclose
auto allclose(E1 &&e1, E2 &&e2, double rtol=1e-05, double atol=1e-08) noexcept
Check if all elements in e1 are close to the corresponding elements in e2.
Definition xmath.hpp:1775
xt::isfinite
auto isfinite(E &&e) noexcept -> detail::xfunction_type_t< math::isfinite_fun, E >
finite value check
Definition xmath.hpp:1640
xt::isnan
auto isnan(E &&e) noexcept -> detail::xfunction_type_t< math::isnan_fun, E >
NaN check.
Definition xmath.hpp:1670
xt::isclose
auto isclose(E1 &&e1, E2 &&e2, double rtol=1e-05, double atol=1e-08, bool equal_nan=false) noexcept
Element-wise closeness detection.
Definition xmath.hpp:1752
xt::isinf
auto isinf(E &&e) noexcept -> detail::xfunction_type_t< math::isinf_fun, E >
infinity check
Definition xmath.hpp:1655
xt::not_equal
auto not_equal(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< detail::not_equal_to, E1, E2 >
Element-wise inequality.
Definition xoperation.hpp:689
xt::equal
auto equal(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< detail::equal_to, E1, E2 >
Element-wise equality.
Definition xoperation.hpp:673
xt::greater
auto greater(E1 &&e1, E2 &&e2) noexcept -> decltype(std::forward< E1 >(e1) > std::forward< E2 >(e2))
Greater than.
Definition xoperation.hpp:740
xt::lgamma
auto lgamma(E &&e) noexcept -> detail::xfunction_type_t< math::lgamma_fun, E >
Natural logarithm of the gamma function.
Definition xmath.hpp:1516
xt::erfc
auto erfc(E &&e) noexcept -> detail::xfunction_type_t< math::erfc_fun, E >
Complementary error function.
Definition xmath.hpp:1486
xt::erf
auto erf(E &&e) noexcept -> detail::xfunction_type_t< math::erf_fun, E >
Error function.
Definition xmath.hpp:1471
xt::tgamma
auto tgamma(E &&e) noexcept -> detail::xfunction_type_t< math::tgamma_fun, E >
Gamma function.
Definition xmath.hpp:1501
xt::log1p
auto log1p(E &&e) noexcept -> detail::xfunction_type_t< math::log1p_fun, E >
Natural logarithm of one plus function.
Definition xmath.hpp:990
xt::expm1
auto expm1(E &&e) noexcept -> detail::xfunction_type_t< math::expm1_fun, E >
Natural exponential minus one function.
Definition xmath.hpp:930
xt::exp2
auto exp2(E &&e) noexcept -> detail::xfunction_type_t< math::exp2_fun, E >
Base 2 exponential function.
Definition xmath.hpp:915
xt::log
auto log(E &&e) noexcept -> detail::xfunction_type_t< math::log_fun, E >
Natural logarithm function.
Definition xmath.hpp:945
xt::log2
auto log2(E &&e) noexcept -> detail::xfunction_type_t< math::log2_fun, E >
Base 2 logarithm function.
Definition xmath.hpp:975
xt::exp
auto exp(E &&e) noexcept -> detail::xfunction_type_t< math::exp_fun, E >
Natural exponential function.
Definition xmath.hpp:900
xt::log10
auto log10(E &&e) noexcept -> detail::xfunction_type_t< math::log10_fun, E >
Base 10 logarithm function.
Definition xmath.hpp:960
xt::asinh
auto asinh(E &&e) noexcept -> detail::xfunction_type_t< math::asinh_fun, E >
Inverse hyperbolic sine function.
Definition xmath.hpp:1418
xt::tanh
auto tanh(E &&e) noexcept -> detail::xfunction_type_t< math::tanh_fun, E >
Hyperbolic tangent function.
Definition xmath.hpp:1403
xt::cosh
auto cosh(E &&e) noexcept -> detail::xfunction_type_t< math::cosh_fun, E >
Hyperbolic cosine function.
Definition xmath.hpp:1388
xt::sinh
auto sinh(E &&e) noexcept -> detail::xfunction_type_t< math::sinh_fun, E >
Hyperbolic sine function.
Definition xmath.hpp:1373
xt::acosh
auto acosh(E &&e) noexcept -> detail::xfunction_type_t< math::acosh_fun, E >
Inverse hyperbolic cosine function.
Definition xmath.hpp:1433
xt::atanh
auto atanh(E &&e) noexcept -> detail::xfunction_type_t< math::atanh_fun, E >
Inverse hyperbolic tangent function.
Definition xmath.hpp:1448
xt::all
bool all(E &&e)
Any.
Definition xoperation.hpp:959
xt::where
auto where(E1 &&e1, E2 &&e2, E3 &&e3) noexcept -> detail::xfunction_type_t< detail::conditional_ternary, E1, E2, E3 >
Ternary selection.
Definition xoperation.hpp:776
xt::nanmax
auto nanmax(E &&e, X &&axes, EVS es=EVS())
Maximum element along given axes, ignoring NaNs.
Definition xmath.hpp:2442
xt::nancumsum
auto nancumsum(E &&e, std::ptrdiff_t axis)
Cumulative sum, replacing nan with 0.
Definition xmath.hpp:2592
xt::nancumprod
auto nancumprod(E &&e, std::ptrdiff_t axis)
Cumulative product, replacing nan with 1.
Definition xmath.hpp:2627
xt::nanmean
auto nanmean(E &&e, X &&axes, EVS es=EVS())
Mean of elements over given axes, excluding NaNs.
Definition xmath.hpp:2711
xt::nanmin
auto nanmin(E &&e, X &&axes, EVS es=EVS())
Minimum element over given axes, ignoring NaNs.
Definition xmath.hpp:2427
xt::nanprod
auto nanprod(E &&e, X &&axes, EVS es=EVS())
Product of elements over given axes, replacing NaN with 1.
Definition xmath.hpp:2476
xt::nansum
auto nansum(E &&e, X &&axes, EVS es=EVS())
Sum of elements over given axes, replacing NaN with 0.
Definition xmath.hpp:2459
xt::nan_to_num
auto nan_to_num(E &&e)
Convert nan or +/- inf to numbers.
Definition xmath.hpp:2409
xt::ceil
auto ceil(E &&e) noexcept -> detail::xfunction_type_t< math::ceil_fun, E >
ceil function.
Definition xmath.hpp:1539
xt::trunc
auto trunc(E &&e) noexcept -> detail::xfunction_type_t< math::trunc_fun, E >
trunc function.
Definition xmath.hpp:1569
xt::nearbyint
auto nearbyint(E &&e) noexcept -> detail::xfunction_type_t< math::nearbyint_fun, E >
nearbyint function.
Definition xmath.hpp:1601
xt::floor
auto floor(E &&e) noexcept -> detail::xfunction_type_t< math::floor_fun, E >
floor function.
Definition xmath.hpp:1554
xt::round
auto round(E &&e) noexcept -> detail::xfunction_type_t< math::round_fun, E >
round function.
Definition xmath.hpp:1585
xt::rint
auto rint(E &&e) noexcept -> detail::xfunction_type_t< math::rint_fun, E >
rint function.
Definition xmath.hpp:1617
xt::cube
auto cube(E1 &&e1) noexcept
Cube power function, equivalent to e1 * e1 * e1.
Definition xmath.hpp:1120
xt::sqrt
auto sqrt(E &&e) noexcept -> detail::xfunction_type_t< math::sqrt_fun, E >
Square root function.
Definition xmath.hpp:1201
xt::square
auto square(E1 &&e1) noexcept
Square power function, equivalent to e1 * e1.
Definition xmath.hpp:1101
xt::pow
auto pow(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::pow_fun, E1, E2 >
Power function.
Definition xmath.hpp:1015
xt::hypot
auto hypot(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::hypot_fun, E1, E2 >
Hypotenuse function.
Definition xmath.hpp:1234
xt::cbrt
auto cbrt(E &&e) noexcept -> detail::xfunction_type_t< math::cbrt_fun, E >
Cubic root function.
Definition xmath.hpp:1216
xt::sum
auto sum(E &&e, X &&axes, EVS es=EVS())
Sum of elements over given axes.
Definition xmath.hpp:1803
xt::prod
auto prod(E &&e, X &&axes, EVS es=EVS())
Product of elements over given axes.
Definition xmath.hpp:1823
xt::trapz
auto trapz(const xexpression< T > &y, double dx=1.0, std::ptrdiff_t axis=-1)
Integrate along the given axis using the composite trapezoidal rule.
Definition xmath.hpp:2910
xt::diff
auto diff(const xexpression< T > &a, std::size_t n=1, std::ptrdiff_t axis=-1)
Calculate the n-th discrete difference along the given axis.
Definition xmath.hpp:2873
xt::minmax
auto minmax(E &&e, EVS es=EVS())
Minimum and maximum among the elements of an array or expression.
Definition xmath.hpp:2196
xt::average
auto average(E &&e, W &&weights, X &&axes, EVS ev=EVS())
Average of elements over given axes using weights.
Definition xmath.hpp:1938
xt::mean
auto mean(E &&e, X &&axes, EVS es=EVS())
Mean of elements over given axes.
Definition xmath.hpp:1897
xt::atan
auto atan(E &&e) noexcept -> detail::xfunction_type_t< math::atan_fun, E >
Arctangent function.
Definition xmath.hpp:1332
xt::atan2
auto atan2(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::atan2_fun, E1, E2 >
Artangent function, using signs to determine quadrants.
Definition xmath.hpp:1350
xt::asin
auto asin(E &&e) noexcept -> detail::xfunction_type_t< math::asin_fun, E >
Arcsine function.
Definition xmath.hpp:1302
xt::cos
auto cos(E &&e) noexcept -> detail::xfunction_type_t< math::cos_fun, E >
Cosine function.
Definition xmath.hpp:1272
xt::sin
auto sin(E &&e) noexcept -> detail::xfunction_type_t< math::sin_fun, E >
Sine function.
Definition xmath.hpp:1257
xt::tan
auto tan(E &&e) noexcept -> detail::xfunction_type_t< math::tan_fun, E >
Tangent function.
Definition xmath.hpp:1287
xt::acos
auto acos(E &&e) noexcept -> detail::xfunction_type_t< math::acos_fun, E >
Arccosine function.
Definition xmath.hpp:1317
xt::conj
auto conj(E &&e) noexcept
Return an xt::xfunction evaluating to the complex conjugate of the given expression.
Definition xcomplex.hpp:207
xt::eval
auto eval(T &&t) -> std::enable_if_t< detail::is_container< std::decay_t< T > >::value, T && >
Force evaluation of xexpression.
Definition xeval.hpp:46
xt::transpose
auto transpose(E &&e) noexcept
Returns a transpose view by reversing the dimensions of xexpression e.
Definition xmanipulation.hpp:239
xt::same_shape
bool same_shape(const S1 &s1, const S2 &s2) noexcept
Check if two objects have the same shape.
Definition xshape.hpp:109
xt
standard mathematical functions for xexpressions
Definition xchunked_array.hpp:20
xt::stack
auto stack(std::tuple< CT... > &&t, std::size_t axis=0)
Stack xexpressions along axis.
Definition xbuilder.hpp:883
xt::range
auto range(A start_val, B stop_val)
Select a range from start_val to stop_val (excluded).
Definition xslice.hpp:734
xt::arange
auto arange(T start, T stop, S step=1) noexcept
Generates numbers evenly spaced within given half-open interval [start, stop).
Definition xbuilder.hpp:432
xt::all
auto all() noexcept
Returns a slice representing a full dimension, to be used as an argument of view function.
Definition xslice.hpp:221
xt::xstrided_slice_vector
std::vector< xstrided_slice< std::ptrdiff_t > > xstrided_slice_vector
vector of slices used to build a xstrided_view
Definition xstrided_view.hpp:341
xt::make_lambda_xfunction
auto make_lambda_xfunction(F &&lambda, E &&... args)
Create a xfunction from a lambda.
Definition xmath.hpp:1085
xt::reduce
auto reduce(F &&f, E &&e, X &&axes, EVS &&options=EVS())
Returns an xexpression applying the specified reducing function to an expression over the given axes.
Definition xreducer.hpp:993
xt::layout_type
layout_type
Definition xlayout.hpp:24
xt::zeros
auto zeros(S shape) noexcept
Returns an xexpression containing zeros of the specified shape.
Definition xbuilder.hpp:66
xt::accumulate
auto accumulate(F &&f, E &&e, EVS evaluation_strategy=EVS())
Accumulate and flatten array NOTE This function is not lazy!
Definition xaccumulator.hpp:344
xt::xtensor
xtensor_container< uvector< T, A >, N, L > xtensor
Alias template on xtensor_container with default parameters for data container type.
Definition xtensor_forward.hpp:137
xt::diagonal
auto diagonal(E &&arr, int offset=0, std::size_t axis_1=0, std::size_t axis_2=1)
Returns the elements on the diagonal of arr If arr has more than two dimensions, then the axes specif...
Definition xbuilder.hpp:1143
xt::make_xshared
xshared_expression< E > make_xshared(xexpression< E > &&expr)
Helper function to create shared expression from any xexpression.
Definition xexpression.hpp:731
xt::strided_view
auto strided_view(E &&e, S &&shape, X &&stride, std::size_t offset=0, layout_type layout=L) noexcept
Construct a strided view from an xexpression, shape, strides and offset.
Definition xstrided_view.hpp:726
xt::diag
auto diag(E &&arr, int k=0)
xexpression with values of arr on the diagonal, zeroes otherwise
Definition xbuilder.hpp:1195
xt::xtuple
auto xtuple(Types &&... args)
Creates tuples from arguments for concatenate and stack.
Definition xbuilder.hpp:707
xt::cov
auto cov(const E1 &x, const E1 &y=E1())
Returns the covariance matrix.
Definition xmath.hpp:3159
xt::const_value
Definition xreducer.hpp:573
xt::convolve_mode::full
Definition xmath.hpp:3212
xt::convolve_mode::valid
Definition xmath.hpp:3208
xt::is_reducer_options
Definition xreducer.hpp:175
xt::math::clamp_fun
Definition xmath.hpp:605
xt::math::sign_fun
Definition xmath.hpp:857
xt::math::sign_impl
Definition xmath.hpp:829
xt::numeric_constants
Definition xmath.hpp:40