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
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 "xaccumulator.hpp"
28#include "xeval.hpp"
29#include "xmanipulation.hpp"
30#include "xoperation.hpp"
31#include "xreducer.hpp"
32#include "xslice.hpp"
33#include "xstrided_view.hpp"
34#include "xtensor_config.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(xtl::negation<is_reducer_options<X>>, xtl::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(xtl::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
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(v < lo, lo, xtl::select(hi < v, hi, v));
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(v < lo, lo, xt_simd::select(hi < v, hi, v));
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
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
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
1091#define XTENSOR_GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
1092
1093// Workaround for MSVC 2015 & GCC 4.9
1094#if (defined(_MSC_VER) && _MSC_VER < 1910) || (defined(__GNUC__) && GCC_VERSION < 49999)
1095#define XTENSOR_DISABLE_LAMBDA_FCT
1096#endif
1097
1098#ifdef XTENSOR_DISABLE_LAMBDA_FCT
1099 struct square_fct
1100 {
1101 template <class T>
1102 auto operator()(T x) const -> decltype(x * x)
1103 {
1104 return x * x;
1105 }
1106 };
1107
1108 struct cube_fct
1109 {
1110 template <class T>
1111 auto operator()(T x) const -> decltype(x * x * x)
1112 {
1113 return x * x * x;
1114 }
1115 };
1116#endif
1117
1127 template <class E1>
1128 inline auto square(E1&& e1) noexcept
1129 {
1130#ifdef XTENSOR_DISABLE_LAMBDA_FCT
1131 return make_lambda_xfunction(square_fct{}, std::forward<E1>(e1));
1132#else
1133 auto fnct = [](auto x) -> decltype(x * x)
1134 {
1135 return x * x;
1136 };
1137 return make_lambda_xfunction(std::move(fnct), std::forward<E1>(e1));
1138#endif
1139 }
1140
1150 template <class E1>
1151 inline auto cube(E1&& e1) noexcept
1152 {
1153#ifdef XTENSOR_DISABLE_LAMBDA_FCT
1154 return make_lambda_xfunction(cube_fct{}, std::forward<E1>(e1));
1155#else
1156 auto fnct = [](auto x) -> decltype(x * x * x)
1157 {
1158 return x * x * x;
1159 };
1160 return make_lambda_xfunction(std::move(fnct), std::forward<E1>(e1));
1161#endif
1162 }
1163
1164#undef XTENSOR_GCC_VERSION
1165#undef XTENSOR_DISABLE_LAMBDA_FCT
1166
1167 namespace detail
1168 {
1169 // Thanks to Matt Pharr in http://pbrt.org/hair.pdf
1170 template <std::size_t N>
1171 struct pow_impl;
1172
1173 template <std::size_t N>
1174 struct pow_impl
1175 {
1176 template <class T>
1177 auto operator()(T v) const -> decltype(v * v)
1178 {
1179 T temp = pow_impl<N / 2>{}(v);
1180 return temp * temp * pow_impl<N & 1>{}(v);
1181 }
1182 };
1183
1184 template <>
1185 struct pow_impl<1>
1186 {
1187 template <class T>
1188 auto operator()(T v) const -> T
1189 {
1190 return v;
1191 }
1192 };
1193
1194 template <>
1195 struct pow_impl<0>
1196 {
1197 template <class T>
1198 auto operator()(T /*v*/) const -> T
1199 {
1200 return T(1);
1201 }
1202 };
1203 }
1204
1222 template <std::size_t N, class E>
1223 inline auto pow(E&& e) noexcept
1224 {
1225 static_assert(N > 0, "integer power cannot be negative");
1226 return make_lambda_xfunction(detail::pow_impl<N>{}, std::forward<E>(e));
1227 }
1228
1238 template <class E>
1239 inline auto sqrt(E&& e) noexcept -> detail::xfunction_type_t<math::sqrt_fun, E>
1240 {
1241 return detail::make_xfunction<math::sqrt_fun>(std::forward<E>(e));
1242 }
1243
1253 template <class E>
1254 inline auto cbrt(E&& e) noexcept -> detail::xfunction_type_t<math::cbrt_fun, E>
1255 {
1256 return detail::make_xfunction<math::cbrt_fun>(std::forward<E>(e));
1257 }
1258
1271 template <class E1, class E2>
1272 inline auto hypot(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::hypot_fun, E1, E2>
1273 {
1274 return detail::make_xfunction<math::hypot_fun>(std::forward<E1>(e1), std::forward<E2>(e2));
1275 }
1276
1277 /***************************
1278 * trigonometric functions *
1279 ***************************/
1280
1294 template <class E>
1295 inline auto sin(E&& e) noexcept -> detail::xfunction_type_t<math::sin_fun, E>
1296 {
1297 return detail::make_xfunction<math::sin_fun>(std::forward<E>(e));
1298 }
1299
1309 template <class E>
1310 inline auto cos(E&& e) noexcept -> detail::xfunction_type_t<math::cos_fun, E>
1311 {
1312 return detail::make_xfunction<math::cos_fun>(std::forward<E>(e));
1313 }
1314
1324 template <class E>
1325 inline auto tan(E&& e) noexcept -> detail::xfunction_type_t<math::tan_fun, E>
1326 {
1327 return detail::make_xfunction<math::tan_fun>(std::forward<E>(e));
1328 }
1329
1339 template <class E>
1340 inline auto asin(E&& e) noexcept -> detail::xfunction_type_t<math::asin_fun, E>
1341 {
1342 return detail::make_xfunction<math::asin_fun>(std::forward<E>(e));
1343 }
1344
1354 template <class E>
1355 inline auto acos(E&& e) noexcept -> detail::xfunction_type_t<math::acos_fun, E>
1356 {
1357 return detail::make_xfunction<math::acos_fun>(std::forward<E>(e));
1358 }
1359
1369 template <class E>
1370 inline auto atan(E&& e) noexcept -> detail::xfunction_type_t<math::atan_fun, E>
1371 {
1372 return detail::make_xfunction<math::atan_fun>(std::forward<E>(e));
1373 }
1374
1387 template <class E1, class E2>
1388 inline auto atan2(E1&& e1, E2&& e2) noexcept -> detail::xfunction_type_t<math::atan2_fun, E1, E2>
1389 {
1390 return detail::make_xfunction<math::atan2_fun>(std::forward<E1>(e1), std::forward<E2>(e2));
1391 }
1392
1393 /************************
1394 * hyperbolic functions *
1395 ************************/
1396
1410 template <class E>
1411 inline auto sinh(E&& e) noexcept -> detail::xfunction_type_t<math::sinh_fun, E>
1412 {
1413 return detail::make_xfunction<math::sinh_fun>(std::forward<E>(e));
1414 }
1415
1425 template <class E>
1426 inline auto cosh(E&& e) noexcept -> detail::xfunction_type_t<math::cosh_fun, E>
1427 {
1428 return detail::make_xfunction<math::cosh_fun>(std::forward<E>(e));
1429 }
1430
1440 template <class E>
1441 inline auto tanh(E&& e) noexcept -> detail::xfunction_type_t<math::tanh_fun, E>
1442 {
1443 return detail::make_xfunction<math::tanh_fun>(std::forward<E>(e));
1444 }
1445
1455 template <class E>
1456 inline auto asinh(E&& e) noexcept -> detail::xfunction_type_t<math::asinh_fun, E>
1457 {
1458 return detail::make_xfunction<math::asinh_fun>(std::forward<E>(e));
1459 }
1460
1470 template <class E>
1471 inline auto acosh(E&& e) noexcept -> detail::xfunction_type_t<math::acosh_fun, E>
1472 {
1473 return detail::make_xfunction<math::acosh_fun>(std::forward<E>(e));
1474 }
1475
1485 template <class E>
1486 inline auto atanh(E&& e) noexcept -> detail::xfunction_type_t<math::atanh_fun, E>
1487 {
1488 return detail::make_xfunction<math::atanh_fun>(std::forward<E>(e));
1489 }
1490
1491 /*****************************
1492 * error and gamma functions *
1493 *****************************/
1494
1508 template <class E>
1509 inline auto erf(E&& e) noexcept -> detail::xfunction_type_t<math::erf_fun, E>
1510 {
1511 return detail::make_xfunction<math::erf_fun>(std::forward<E>(e));
1512 }
1513
1523 template <class E>
1524 inline auto erfc(E&& e) noexcept -> detail::xfunction_type_t<math::erfc_fun, E>
1525 {
1526 return detail::make_xfunction<math::erfc_fun>(std::forward<E>(e));
1527 }
1528
1538 template <class E>
1539 inline auto tgamma(E&& e) noexcept -> detail::xfunction_type_t<math::tgamma_fun, E>
1540 {
1541 return detail::make_xfunction<math::tgamma_fun>(std::forward<E>(e));
1542 }
1543
1553 template <class E>
1554 inline auto lgamma(E&& e) noexcept -> detail::xfunction_type_t<math::lgamma_fun, E>
1555 {
1556 return detail::make_xfunction<math::lgamma_fun>(std::forward<E>(e));
1557 }
1558
1559 /*********************************************
1560 * nearest integer floating point operations *
1561 *********************************************/
1562
1576 template <class E>
1577 inline auto ceil(E&& e) noexcept -> detail::xfunction_type_t<math::ceil_fun, E>
1578 {
1579 return detail::make_xfunction<math::ceil_fun>(std::forward<E>(e));
1580 }
1581
1591 template <class E>
1592 inline auto floor(E&& e) noexcept -> detail::xfunction_type_t<math::floor_fun, E>
1593 {
1594 return detail::make_xfunction<math::floor_fun>(std::forward<E>(e));
1595 }
1596
1606 template <class E>
1607 inline auto trunc(E&& e) noexcept -> detail::xfunction_type_t<math::trunc_fun, E>
1608 {
1609 return detail::make_xfunction<math::trunc_fun>(std::forward<E>(e));
1610 }
1611
1622 template <class E>
1623 inline auto round(E&& e) noexcept -> detail::xfunction_type_t<math::round_fun, E>
1624 {
1625 return detail::make_xfunction<math::round_fun>(std::forward<E>(e));
1626 }
1627
1638 template <class E>
1639 inline auto nearbyint(E&& e) noexcept -> detail::xfunction_type_t<math::nearbyint_fun, E>
1640 {
1641 return detail::make_xfunction<math::nearbyint_fun>(std::forward<E>(e));
1642 }
1643
1654 template <class E>
1655 inline auto rint(E&& e) noexcept -> detail::xfunction_type_t<math::rint_fun, E>
1656 {
1657 return detail::make_xfunction<math::rint_fun>(std::forward<E>(e));
1658 }
1659
1660 /****************************
1661 * classification functions *
1662 ****************************/
1663
1677 template <class E>
1678 inline auto isfinite(E&& e) noexcept -> detail::xfunction_type_t<math::isfinite_fun, E>
1679 {
1680 return detail::make_xfunction<math::isfinite_fun>(std::forward<E>(e));
1681 }
1682
1692 template <class E>
1693 inline auto isinf(E&& e) noexcept -> detail::xfunction_type_t<math::isinf_fun, E>
1694 {
1695 return detail::make_xfunction<math::isinf_fun>(std::forward<E>(e));
1696 }
1697
1707 template <class E>
1708 inline auto isnan(E&& e) noexcept -> detail::xfunction_type_t<math::isnan_fun, E>
1709 {
1710 return detail::make_xfunction<math::isnan_fun>(std::forward<E>(e));
1711 }
1712
1713 namespace detail
1714 {
1715 template <class FUNCTOR, class T, std::size_t... Is>
1716 inline auto get_functor(T&& args, std::index_sequence<Is...>)
1717 {
1718 return FUNCTOR(std::get<Is>(args)...);
1719 }
1720
1721 template <class F, class... A, class... E>
1722 inline auto make_xfunction(std::tuple<A...>&& f_args, E&&... e) noexcept
1723 {
1724 using functor_type = F;
1725 using expression_tag = xexpression_tag_t<E...>;
1726 using type = select_xfunction_expression_t<expression_tag, functor_type, const_xclosure_t<E>...>;
1727 auto functor = get_functor<functor_type>(
1728 std::forward<std::tuple<A...>>(f_args),
1729 std::make_index_sequence<sizeof...(A)>{}
1730 );
1731 return type(std::move(functor), std::forward<E>(e)...);
1732 }
1733
1734 struct isclose
1735 {
1736 using result_type = bool;
1737
1738 isclose(double rtol, double atol, bool equal_nan)
1739 : m_rtol(rtol)
1740 , m_atol(atol)
1741 , m_equal_nan(equal_nan)
1742 {
1743 }
1744
1745 template <class A1, class A2>
1746 bool operator()(const A1& a, const A2& b) const
1747 {
1748 using internal_type = xtl::promote_type_t<A1, A2, double>;
1749 if (math::isnan(a) && math::isnan(b))
1750 {
1751 return m_equal_nan;
1752 }
1753 if (math::isinf(a) && math::isinf(b))
1754 {
1755 // check for both infinity signs equal
1756 return a == b;
1757 }
1758 auto d = math::abs(internal_type(a) - internal_type(b));
1759 return d <= m_atol
1760 || d <= m_rtol
1761 * double((std::max)(math::abs(internal_type(a)), math::abs(internal_type(b)))
1762 );
1763 }
1764
1765 private:
1766
1767 double m_rtol;
1768 double m_atol;
1769 bool m_equal_nan;
1770 };
1771 }
1772
1788 template <class E1, class E2>
1789 inline auto
1790 isclose(E1&& e1, E2&& e2, double rtol = 1e-05, double atol = 1e-08, bool equal_nan = false) noexcept
1791 {
1792 return detail::make_xfunction<detail::isclose>(
1793 std::make_tuple(rtol, atol, equal_nan),
1794 std::forward<E1>(e1),
1795 std::forward<E2>(e2)
1796 );
1797 }
1798
1812 template <class E1, class E2>
1813 inline auto allclose(E1&& e1, E2&& e2, double rtol = 1e-05, double atol = 1e-08) noexcept
1814 {
1815 return xt::all(isclose(std::forward<E1>(e1), std::forward<E2>(e2), rtol, atol));
1816 }
1817
1818 /**********************
1819 * Reducing functions *
1820 **********************/
1821
1841 XTENSOR_REDUCER_FUNCTION(sum, detail::plus, typename std::decay_t<E>::value_type, 0)
1842
1843
1861 XTENSOR_REDUCER_FUNCTION(prod, detail::multiplies, typename std::decay_t<E>::value_type, 1)
1862
1863 namespace detail
1864 {
1865 template <class T, class S, class ST>
1866 inline auto mean_division(S&& s, ST e_size)
1867 {
1868 using value_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
1869 // Avoids floating point exception when s.size is 0
1870 value_type div = s.size() != ST(0) ? static_cast<value_type>(e_size / s.size()) : value_type(0);
1871 return std::move(s) / std::move(div);
1872 }
1873
1874 template <
1875 class T,
1876 class E,
1877 class X,
1878 class D,
1879 class EVS,
1880 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>, xtl::is_integral<D>)>
1881 inline auto mean(E&& e, X&& axes, const D& ddof, EVS es)
1882 {
1883 // sum cannot always be a double. It could be a complex number which cannot operate on
1884 // std::plus<double>.
1885 using size_type = typename std::decay_t<E>::size_type;
1886 const size_type size = e.size();
1887 XTENSOR_ASSERT(static_cast<size_type>(ddof) <= size);
1888 auto s = sum<T>(std::forward<E>(e), std::forward<X>(axes), es);
1889 return mean_division<T>(std::move(s), size - static_cast<size_type>(ddof));
1890 }
1891
1892 template <class T, class E, class I, std::size_t N, class D, class EVS>
1893 inline auto mean(E&& e, const I (&axes)[N], const D& ddof, EVS es)
1894 {
1895 using size_type = typename std::decay_t<E>::size_type;
1896 const size_type size = e.size();
1897 XTENSOR_ASSERT(static_cast<size_type>(ddof) <= size);
1898 auto s = sum<T>(std::forward<E>(e), axes, es);
1899 return mean_division<T>(std::move(s), size - static_cast<size_type>(ddof));
1900 }
1901
1902 template <class T, class E, class D, class EVS, XTL_REQUIRES(is_reducer_options<EVS>, xtl::is_integral<D>)>
1903 inline auto mean_noaxis(E&& e, const D& ddof, EVS es)
1904 {
1905 using value_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
1906 using size_type = typename std::decay_t<E>::size_type;
1907 const size_type size = e.size();
1908 XTENSOR_ASSERT(static_cast<size_type>(ddof) <= size);
1909 auto s = sum<T>(std::forward<E>(e), es);
1910 return std::move(s) / static_cast<value_type>((size - static_cast<size_type>(ddof)));
1911 }
1912 }
1913
1929 template <
1930 class T = void,
1931 class E,
1932 class X,
1933 class EVS = DEFAULT_STRATEGY_REDUCERS,
1934 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>)>
1935 inline auto mean(E&& e, X&& axes, EVS es = EVS())
1936 {
1937 return detail::mean<T>(std::forward<E>(e), std::forward<X>(axes), 0u, es);
1938 }
1939
1940 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
1941 inline auto mean(E&& e, EVS es = EVS())
1942 {
1943 return detail::mean_noaxis<T>(std::forward<E>(e), 0u, es);
1944 }
1945
1946 template <class T = void, class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
1947 inline auto mean(E&& e, const I (&axes)[N], EVS es = EVS())
1948 {
1949 return detail::mean<T>(std::forward<E>(e), axes, 0u, es);
1950 }
1951
1969 template <
1970 class T = void,
1971 class E,
1972 class W,
1973 class X,
1974 class EVS = DEFAULT_STRATEGY_REDUCERS,
1975 XTL_REQUIRES(is_reducer_options<EVS>, xtl::negation<xtl::is_integral<X>>)>
1976 inline auto average(E&& e, W&& weights, X&& axes, EVS ev = EVS())
1977 {
1978 xindex_type_t<typename std::decay_t<E>::shape_type> broadcast_shape;
1979 xt::resize_container(broadcast_shape, e.dimension());
1980 auto ax = normalize_axis(e, axes);
1981 if (weights.dimension() == 1)
1982 {
1983 if (weights.size() != e.shape()[ax[0]])
1984 {
1985 XTENSOR_THROW(std::runtime_error, "Weights need to have the same shape as expression at axes.");
1986 }
1987
1988 std::fill(broadcast_shape.begin(), broadcast_shape.end(), std::size_t(1));
1989 broadcast_shape[ax[0]] = weights.size();
1990 }
1991 else
1992 {
1993 if (!same_shape(e.shape(), weights.shape()))
1994 {
1995 XTENSOR_THROW(
1996 std::runtime_error,
1997 "Weights with dim > 1 need to have the same shape as expression."
1998 );
1999 }
2000
2001 std::copy(e.shape().begin(), e.shape().end(), broadcast_shape.begin());
2002 }
2003
2004 constexpr layout_type L = default_assignable_layout(std::decay_t<W>::static_layout);
2005 auto weights_view = reshape_view<L>(std::forward<W>(weights), std::move(broadcast_shape));
2006 auto scl = sum<T>(weights_view, ax, xt::evaluation_strategy::immediate);
2007 return sum<T>(std::forward<E>(e) * std::move(weights_view), std::move(ax), ev) / std::move(scl);
2008 }
2009
2010 template <
2011 class T = void,
2012 class E,
2013 class W,
2014 class X,
2015 class EVS = DEFAULT_STRATEGY_REDUCERS,
2016 XTL_REQUIRES(is_reducer_options<EVS>, xtl::is_integral<X>)>
2017 inline auto average(E&& e, W&& weights, X axis, EVS ev = EVS())
2018 {
2019 return average(std::forward<E>(e), std::forward<W>(weights), {axis}, std::forward<EVS>(ev));
2020 }
2021
2022 template <class T = void, class E, class W, class X, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2023 inline auto average(E&& e, W&& weights, const X (&axes)[N], EVS ev = EVS())
2024 {
2025 // need to select the X&& overload and forward to different type
2026 using ax_t = std::array<std::size_t, N>;
2027 return average<T>(std::forward<E>(e), std::forward<W>(weights), xt::forward_normalize<ax_t>(e, axes), ev);
2028 }
2029
2030 template <class T = void, class E, class W, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2031 inline auto average(E&& e, W&& weights, EVS ev = EVS())
2032 {
2033 if (weights.dimension() != e.dimension()
2034 || !std::equal(weights.shape().begin(), weights.shape().end(), e.shape().begin()))
2035 {
2036 XTENSOR_THROW(std::runtime_error, "Weights need to have the same shape as expression.");
2037 }
2038
2039 auto div = sum<T>(weights, evaluation_strategy::immediate)();
2040 auto s = sum<T>(std::forward<E>(e) * std::forward<W>(weights), ev) / std::move(div);
2041 return s;
2042 }
2043
2044 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2045 inline auto average(E&& e, EVS ev = EVS())
2046 {
2047 return mean<T>(e, ev);
2048 }
2049
2050 namespace detail
2051 {
2052 template <typename E>
2053 std::enable_if_t<std::is_lvalue_reference<E>::value, E> shared_forward(E e) noexcept
2054 {
2055 return e;
2056 }
2057
2058 template <typename E>
2059 std::enable_if_t<!std::is_lvalue_reference<E>::value, xshared_expression<E>> shared_forward(E e) noexcept
2060 {
2061 return make_xshared(std::move(e));
2062 }
2063 }
2064
2065 template <
2066 class T = void,
2067 class E,
2068 class D,
2069 class EVS = DEFAULT_STRATEGY_REDUCERS,
2070 XTL_REQUIRES(is_reducer_options<EVS>, xtl::is_integral<D>)>
2071 inline auto variance(E&& e, const D& ddof, EVS es = EVS())
2072 {
2073 auto cached_mean = mean<T>(e, es)();
2074 return detail::mean_noaxis<T>(square(std::forward<E>(e) - std::move(cached_mean)), ddof, es);
2075 }
2076
2077 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2078 inline auto variance(E&& e, EVS es = EVS())
2079 {
2080 return variance<T>(std::forward<E>(e), 0u, es);
2081 }
2082
2083 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2084 inline auto stddev(E&& e, EVS es = EVS())
2085 {
2086 return sqrt(variance<T>(std::forward<E>(e), es));
2087 }
2088
2113 template <
2114 class T = void,
2115 class E,
2116 class X,
2117 class D,
2118 class EVS = DEFAULT_STRATEGY_REDUCERS,
2119 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>, xtl::is_integral<D>)>
2120 inline auto variance(E&& e, X&& axes, const D& ddof, EVS es = EVS())
2121 {
2122 decltype(auto) sc = detail::shared_forward<E>(e);
2123 // note: forcing copy of first axes argument -- is there a better solution?
2124 auto axes_copy = axes;
2125 // always eval to prevent repeated evaluations in the next calls
2126 auto inner_mean = eval(mean<T>(sc, std::move(axes_copy), evaluation_strategy::immediate));
2127
2128 // fake keep_dims = 1
2129 // Since the inner_shape might have a reference semantic (e.g. xbuffer_adaptor in bindings)
2130 // We need to map it to another type before modifying it.
2131 // We pragmatically abuse `get_strides_t`
2132 using tmp_shape_t = get_strides_t<typename std::decay_t<E>::shape_type>;
2133 tmp_shape_t keep_dim_shape = xtl::forward_sequence<tmp_shape_t, decltype(e.shape())>(e.shape());
2134 for (const auto& el : axes)
2135 {
2136 keep_dim_shape[el] = 1u;
2137 }
2138
2139 auto mrv = reshape_view<XTENSOR_DEFAULT_LAYOUT>(std::move(inner_mean), std::move(keep_dim_shape));
2140 return detail::mean<T>(square(sc - std::move(mrv)), std::forward<X>(axes), ddof, es);
2141 }
2142
2143 template <
2144 class T = void,
2145 class E,
2146 class X,
2147 class EVS = DEFAULT_STRATEGY_REDUCERS,
2148 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>, xtl::negation<xtl::is_integral<std::decay_t<X>>>, is_reducer_options<EVS>)>
2149 inline auto variance(E&& e, X&& axes, EVS es = EVS())
2150 {
2151 return variance<T>(std::forward<E>(e), std::forward<X>(axes), 0u, es);
2152 }
2153
2175 template <
2176 class T = void,
2177 class E,
2178 class X,
2179 class EVS = DEFAULT_STRATEGY_REDUCERS,
2180 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>)>
2181 inline auto stddev(E&& e, X&& axes, EVS es = EVS())
2182 {
2183 return sqrt(variance<T>(std::forward<E>(e), std::forward<X>(axes), es));
2184 }
2185
2186 template <class T = void, class E, class A, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2187 inline auto stddev(E&& e, const A (&axes)[N], EVS es = EVS())
2188 {
2189 return stddev<T>(
2190 std::forward<E>(e),
2191 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2192 es
2193 );
2194 }
2195
2196 template <
2197 class T = void,
2198 class E,
2199 class A,
2200 std::size_t N,
2201 class EVS = DEFAULT_STRATEGY_REDUCERS,
2202 XTL_REQUIRES(is_reducer_options<EVS>)>
2203 inline auto variance(E&& e, const A (&axes)[N], EVS es = EVS())
2204 {
2205 return variance<T>(
2206 std::forward<E>(e),
2207 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2208 es
2209 );
2210 }
2211
2212 template <class T = void, class E, class A, std::size_t N, class D, class EVS = DEFAULT_STRATEGY_REDUCERS>
2213 inline auto variance(E&& e, const A (&axes)[N], const D& ddof, EVS es = EVS())
2214 {
2215 return variance<T>(
2216 std::forward<E>(e),
2217 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2218 ddof,
2219 es
2220 );
2221 }
2222
2233 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2234 inline auto minmax(E&& e, EVS es = EVS())
2235 {
2236 using std::max;
2237 using std::min;
2238 using value_type = typename std::decay_t<E>::value_type;
2239 using result_type = std::array<value_type, 2>;
2240 using init_value_fct = xt::const_value<result_type>;
2241
2242 auto reduce_func = [](auto r, const auto& v)
2243 {
2244 r[0] = (min) (r[0], v);
2245 r[1] = (max) (r[1], v);
2246 return r;
2247 };
2248
2249 auto init_func = init_value_fct(
2250 result_type{std::numeric_limits<value_type>::max(), std::numeric_limits<value_type>::lowest()}
2251 );
2252
2253 auto merge_func = [](auto r, const auto& s)
2254 {
2255 r[0] = (min) (r[0], s[0]);
2256 r[1] = (max) (r[1], s[1]);
2257 return r;
2258 };
2259 return xt::reduce(
2260 make_xreducer_functor(std::move(reduce_func), std::move(init_func), std::move(merge_func)),
2261 std::forward<E>(e),
2262 arange(e.dimension()),
2263 es
2264 );
2265 }
2266
2285 template <class T = void, class E>
2286 inline auto cumsum(E&& e, std::ptrdiff_t axis)
2287 {
2288 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2289 return accumulate(
2290 make_xaccumulator_functor(detail::plus(), detail::accumulator_identity<init_value_type>()),
2291 std::forward<E>(e),
2292 axis
2293 );
2294 }
2295
2296 template <class T = void, class E>
2297 inline auto cumsum(E&& e)
2298 {
2299 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2300 return accumulate(
2301 make_xaccumulator_functor(detail::plus(), detail::accumulator_identity<init_value_type>()),
2302 std::forward<E>(e)
2303 );
2304 }
2305
2320 template <class T = void, class E>
2321 inline auto cumprod(E&& e, std::ptrdiff_t axis)
2322 {
2323 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2324 return accumulate(
2325 make_xaccumulator_functor(detail::multiplies(), detail::accumulator_identity<init_value_type>()),
2326 std::forward<E>(e),
2327 axis
2328 );
2329 }
2330
2331 template <class T = void, class E>
2332 inline auto cumprod(E&& e)
2333 {
2334 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2335 return accumulate(
2336 make_xaccumulator_functor(detail::multiplies(), detail::accumulator_identity<init_value_type>()),
2337 std::forward<E>(e)
2338 );
2339 }
2340
2341 /*****************
2342 * nan functions *
2343 *****************/
2344
2345 namespace detail
2346 {
2347 struct nan_to_num_functor
2348 {
2349 template <class A>
2350 inline auto operator()(const A& a) const
2351 {
2352 if (math::isnan(a))
2353 {
2354 return A(0);
2355 }
2356 if (math::isinf(a))
2357 {
2358 if (a < 0)
2359 {
2360 return std::numeric_limits<A>::lowest();
2361 }
2362 else
2363 {
2364 return (std::numeric_limits<A>::max)();
2365 }
2366 }
2367 return a;
2368 }
2369 };
2370
2371 struct nan_min
2372 {
2373 template <class T, class U>
2374 constexpr auto operator()(const T lhs, const U rhs) const
2375 {
2376 // Clunky expression for working with GCC 4.9
2377 return math::isnan(lhs)
2378 ? rhs
2379 : (math::isnan(rhs) ? lhs
2380 : std::common_type_t<T, U>(
2381 detail::make_xfunction<math::minimum<void>>(lhs, rhs)
2382 ));
2383 }
2384 };
2385
2386 struct nan_max
2387 {
2388 template <class T, class U>
2389 constexpr auto operator()(const T lhs, const U rhs) const
2390 {
2391 // Clunky expression for working with GCC 4.9
2392 return math::isnan(lhs)
2393 ? rhs
2394 : (math::isnan(rhs) ? lhs
2395 : std::common_type_t<T, U>(
2396 detail::make_xfunction<math::maximum<void>>(lhs, rhs)
2397 ));
2398 }
2399 };
2400
2401 struct nan_plus
2402 {
2403 template <class T, class U>
2404 constexpr auto operator()(const T lhs, const U rhs) const
2405 {
2406 return !math::isnan(rhs) ? lhs + rhs : lhs;
2407 }
2408 };
2409
2410 struct nan_multiplies
2411 {
2412 template <class T, class U>
2413 constexpr auto operator()(const T lhs, const U rhs) const
2414 {
2415 return !math::isnan(rhs) ? lhs * rhs : lhs;
2416 }
2417 };
2418
2419 template <class T, int V>
2420 struct nan_init
2421 {
2422 using value_type = T;
2423 using result_type = T;
2424
2425 constexpr result_type operator()(const value_type lhs) const
2426 {
2427 return math::isnan(lhs) ? result_type(V) : lhs;
2428 }
2429 };
2430 }
2431
2446 template <class E>
2447 inline auto nan_to_num(E&& e)
2448 {
2449 return detail::make_xfunction<detail::nan_to_num_functor>(std::forward<E>(e));
2450 }
2451
2465 XTENSOR_REDUCER_FUNCTION(nanmin, detail::nan_min, typename std::decay_t<E>::value_type, std::nan("0"))
2466
2467
2480 XTENSOR_REDUCER_FUNCTION(nanmax, detail::nan_max, typename std::decay_t<E>::value_type, std::nan("0"))
2481
2497 XTENSOR_REDUCER_FUNCTION(nansum, detail::nan_plus, typename std::decay_t<E>::value_type, 0)
2498
2514 XTENSOR_REDUCER_FUNCTION(nanprod, detail::nan_multiplies, typename std::decay_t<E>::value_type, 1)
2515
2516#define COUNT_NON_ZEROS_CONTENT \
2517 using value_type = typename std::decay_t<E>::value_type; \
2518 using result_type = xt::detail::xreducer_size_type_t<value_type>; \
2519 using init_value_fct = xt::const_value<result_type>; \
2520 \
2521 auto init_fct = init_value_fct(0); \
2522 \
2523 auto reduce_fct = [](const auto& lhs, const auto& rhs) \
2524 { \
2525 using value_t = xt::detail::xreducer_temporary_type_t<std::decay_t<decltype(rhs)>>; \
2526 using result_t = std::decay_t<decltype(lhs)>; \
2527 \
2528 return (rhs != value_t(0)) ? lhs + result_t(1) : lhs; \
2529 }; \
2530 auto merge_func = detail::plus();
2531
2532 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2533 inline auto count_nonzero(E&& e, EVS es = EVS())
2534 {
2535 COUNT_NON_ZEROS_CONTENT;
2536 return xt::reduce(
2537 make_xreducer_functor(std::move(reduce_fct), std::move(init_fct), std::move(merge_func)),
2538 std::forward<E>(e),
2539 es
2540 );
2541 }
2542
2543 template <
2544 class E,
2545 class X,
2546 class EVS = DEFAULT_STRATEGY_REDUCERS,
2547 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>, xtl::negation<xtl::is_integral<X>>)>
2548 inline auto count_nonzero(E&& e, X&& axes, EVS es = EVS())
2549 {
2550 COUNT_NON_ZEROS_CONTENT;
2551 return xt::reduce(
2552 make_xreducer_functor(std::move(reduce_fct), std::move(init_fct), std::move(merge_func)),
2553 std::forward<E>(e),
2554 std::forward<X>(axes),
2555 es
2556 );
2557 }
2558
2559 template <
2560 class E,
2561 class X,
2562 class EVS = DEFAULT_STRATEGY_REDUCERS,
2563 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>, xtl::is_integral<X>)>
2564 inline auto count_nonzero(E&& e, X axis, EVS es = EVS())
2565 {
2566 return count_nonzero(std::forward<E>(e), {axis}, es);
2567 }
2568
2569 template <class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2570 inline auto count_nonzero(E&& e, const I (&axes)[N], EVS es = EVS())
2571 {
2572 COUNT_NON_ZEROS_CONTENT;
2573 return xt::reduce(
2574 make_xreducer_functor(std::move(reduce_fct), std::move(init_fct), std::move(merge_func)),
2575 std::forward<E>(e),
2576 axes,
2577 es
2578 );
2579 }
2580
2581#undef COUNT_NON_ZEROS_CONTENT
2582
2583 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2584 inline auto count_nonnan(E&& e, EVS es = EVS())
2585 {
2586 return xt::count_nonzero(!xt::isnan(std::forward<E>(e)), es);
2587 }
2588
2589 template <
2590 class E,
2591 class X,
2592 class EVS = DEFAULT_STRATEGY_REDUCERS,
2593 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>, xtl::negation<xtl::is_integral<X>>)>
2594 inline auto count_nonnan(E&& e, X&& axes, EVS es = EVS())
2595 {
2596 return xt::count_nonzero(!xt::isnan(std::forward<E>(e)), std::forward<X>(axes), es);
2597 }
2598
2599 template <
2600 class E,
2601 class X,
2602 class EVS = DEFAULT_STRATEGY_REDUCERS,
2603 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>, xtl::is_integral<X>)>
2604 inline auto count_nonnan(E&& e, X&& axes, EVS es = EVS())
2605 {
2606 return xt::count_nonzero(!xt::isnan(std::forward<E>(e)), {axes}, es);
2607 }
2608
2609 template <class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2610 inline auto count_nonnan(E&& e, const I (&axes)[N], EVS es = EVS())
2611 {
2612 return xt::count_nonzero(!xt::isnan(std::forward<E>(e)), axes, es);
2613 }
2614
2629 template <class T = void, class E>
2630 inline auto nancumsum(E&& e, std::ptrdiff_t axis)
2631 {
2632 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2633 return accumulate(
2634 make_xaccumulator_functor(detail::nan_plus(), detail::nan_init<init_value_type, 0>()),
2635 std::forward<E>(e),
2636 axis
2637 );
2638 }
2639
2640 template <class T = void, class E>
2641 inline auto nancumsum(E&& e)
2642 {
2643 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2644 return accumulate(
2645 make_xaccumulator_functor(detail::nan_plus(), detail::nan_init<init_value_type, 0>()),
2646 std::forward<E>(e)
2647 );
2648 }
2649
2664 template <class T = void, class E>
2665 inline auto nancumprod(E&& e, std::ptrdiff_t axis)
2666 {
2667 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2668 return accumulate(
2669 make_xaccumulator_functor(detail::nan_multiplies(), detail::nan_init<init_value_type, 1>()),
2670 std::forward<E>(e),
2671 axis
2672 );
2673 }
2674
2675 template <class T = void, class E>
2676 inline auto nancumprod(E&& e)
2677 {
2678 using init_value_type = std::conditional_t<std::is_same<T, void>::value, typename std::decay_t<E>::value_type, T>;
2679 return accumulate(
2680 make_xaccumulator_functor(detail::nan_multiplies(), detail::nan_init<init_value_type, 1>()),
2681 std::forward<E>(e)
2682 );
2683 }
2684
2685 namespace detail
2686 {
2687 template <class T>
2688 struct diff_impl
2689 {
2690 template <class Arg>
2691 inline void operator()(
2692 Arg& ad,
2693 const std::size_t& n,
2694 xstrided_slice_vector& slice1,
2695 xstrided_slice_vector& slice2,
2696 std::size_t saxis
2697 )
2698 {
2699 for (std::size_t i = 0; i < n; ++i)
2700 {
2701 slice2[saxis] = range(xnone(), ad.shape()[saxis] - 1);
2702 ad = strided_view(ad, slice1) - strided_view(ad, slice2);
2703 }
2704 }
2705 };
2706
2707 template <>
2708 struct diff_impl<bool>
2709 {
2710 template <class Arg>
2711 inline void operator()(
2712 Arg& ad,
2713 const std::size_t& n,
2714 xstrided_slice_vector& slice1,
2715 xstrided_slice_vector& slice2,
2716 std::size_t saxis
2717 )
2718 {
2719 for (std::size_t i = 0; i < n; ++i)
2720 {
2721 slice2[saxis] = range(xnone(), ad.shape()[saxis] - 1);
2722 ad = not_equal(strided_view(ad, slice1), strided_view(ad, slice2));
2723 }
2724 }
2725 };
2726 }
2727
2743 template <
2744 class T = void,
2745 class E,
2746 class X,
2747 class EVS = DEFAULT_STRATEGY_REDUCERS,
2748 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>)>
2749 inline auto nanmean(E&& e, X&& axes, EVS es = EVS())
2750 {
2751 decltype(auto) sc = detail::shared_forward<E>(e);
2752 // note: forcing copy of first axes argument -- is there a better solution?
2753 auto axes_copy = axes;
2754 using value_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
2755 using sum_type = typename std::conditional_t<
2756 std::is_same<T, void>::value,
2757 typename std::common_type_t<typename std::decay_t<E>::value_type, value_type>,
2758 T>;
2759 // sum cannot always be a double. It could be a complex number which cannot operate on
2760 // std::plus<double>.
2761 return nansum<sum_type>(sc, std::forward<X>(axes), es)
2762 / xt::cast<value_type>(count_nonnan(sc, std::move(axes_copy), es));
2763 }
2764
2765 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2766 inline auto nanmean(E&& e, EVS es = EVS())
2767 {
2768 decltype(auto) sc = detail::shared_forward<E>(e);
2769 using value_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
2770 using sum_type = typename std::conditional_t<
2771 std::is_same<T, void>::value,
2772 typename std::common_type_t<typename std::decay_t<E>::value_type, value_type>,
2773 T>;
2774 return nansum<sum_type>(sc, es) / xt::cast<value_type>(count_nonnan(sc, es));
2775 }
2776
2777 template <class T = void, class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2778 inline auto nanmean(E&& e, const I (&axes)[N], EVS es = EVS())
2779 {
2780 return nanmean<T>(
2781 std::forward<E>(e),
2782 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2783 es
2784 );
2785 }
2786
2787 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2788 inline auto nanvar(E&& e, EVS es = EVS())
2789 {
2790 decltype(auto) sc = detail::shared_forward<E>(e);
2791 return nanmean<T>(square(sc - nanmean<T>(sc)), es);
2792 }
2793
2794 template <class T = void, class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_reducer_options<EVS>)>
2795 inline auto nanstd(E&& e, EVS es = EVS())
2796 {
2797 return sqrt(nanvar<T>(std::forward<E>(e), es));
2798 }
2799
2820 template <
2821 class T = void,
2822 class E,
2823 class X,
2824 class EVS = DEFAULT_STRATEGY_REDUCERS,
2825 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>)>
2826 inline auto nanvar(E&& e, X&& axes, EVS es = EVS())
2827 {
2828 decltype(auto) sc = detail::shared_forward<E>(e);
2829 // note: forcing copy of first axes argument -- is there a better solution?
2830 auto axes_copy = axes;
2831 using result_type = typename std::conditional_t<std::is_same<T, void>::value, double, T>;
2832 auto inner_mean = nanmean<result_type>(sc, std::move(axes_copy));
2833
2834 // fake keep_dims = 1
2835 // Since the inner_shape might have a reference semantic (e.g. xbuffer_adaptor in bindings)
2836 // We need to map it to another type before modifying it.
2837 // We pragmatically abuse `get_strides_t`
2838 using tmp_shape_t = get_strides_t<typename std::decay_t<E>::shape_type>;
2839 tmp_shape_t keep_dim_shape = xtl::forward_sequence<tmp_shape_t, decltype(e.shape())>(e.shape());
2840 for (const auto& el : axes)
2841 {
2842 keep_dim_shape[el] = 1;
2843 }
2844 auto mrv = reshape_view<XTENSOR_DEFAULT_LAYOUT>(std::move(inner_mean), std::move(keep_dim_shape));
2845 return nanmean<result_type>(square(cast<result_type>(sc) - std::move(mrv)), std::forward<X>(axes), es);
2846 }
2847
2868 template <
2869 class T = void,
2870 class E,
2871 class X,
2872 class EVS = DEFAULT_STRATEGY_REDUCERS,
2873 XTL_REQUIRES(xtl::negation<is_reducer_options<X>>)>
2874 inline auto nanstd(E&& e, X&& axes, EVS es = EVS())
2875 {
2876 return sqrt(nanvar<T>(std::forward<E>(e), std::forward<X>(axes), es));
2877 }
2878
2879 template <class T = void, class E, class A, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2880 inline auto nanstd(E&& e, const A (&axes)[N], EVS es = EVS())
2881 {
2882 return nanstd<T>(
2883 std::forward<E>(e),
2884 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2885 es
2886 );
2887 }
2888
2889 template <class T = void, class E, class A, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
2890 inline auto nanvar(E&& e, const A (&axes)[N], EVS es = EVS())
2891 {
2892 return nanvar<T>(
2893 std::forward<E>(e),
2894 xtl::forward_sequence<std::array<std::size_t, N>, decltype(axes)>(axes),
2895 es
2896 );
2897 }
2898
2910 template <class T>
2911 auto diff(const xexpression<T>& a, std::size_t n = 1, std::ptrdiff_t axis = -1)
2912 {
2913 typename std::decay_t<T>::temporary_type ad = a.derived_cast();
2914 std::size_t saxis = normalize_axis(ad.dimension(), axis);
2915 if (n <= ad.size())
2916 {
2917 if (n != std::size_t(0))
2918 {
2919 xstrided_slice_vector slice1(ad.dimension(), all());
2920 xstrided_slice_vector slice2(ad.dimension(), all());
2921 slice1[saxis] = range(1, xnone());
2922
2923 detail::diff_impl<typename T::value_type> impl;
2924 impl(ad, n, slice1, slice2, saxis);
2925 }
2926 }
2927 else
2928 {
2929 auto shape = ad.shape();
2930 shape[saxis] = std::size_t(0);
2931 ad.resize(shape);
2932 }
2933 return ad;
2934 }
2935
2947 template <class T>
2948 auto trapz(const xexpression<T>& y, double dx = 1.0, std::ptrdiff_t axis = -1)
2949 {
2950 auto& yd = y.derived_cast();
2951 std::size_t saxis = normalize_axis(yd.dimension(), axis);
2952
2953 xstrided_slice_vector slice1(yd.dimension(), all());
2954 xstrided_slice_vector slice2(yd.dimension(), all());
2955 slice1[saxis] = range(1, xnone());
2956 slice2[saxis] = range(xnone(), yd.shape()[saxis] - 1);
2957
2958 auto trap = dx * (strided_view(yd, slice1) + strided_view(yd, slice2)) * 0.5;
2959
2960 return eval(sum(trap, {saxis}));
2961 }
2962
2974 template <class T, class E>
2975 auto trapz(const xexpression<T>& y, const xexpression<E>& x, std::ptrdiff_t axis = -1)
2976 {
2977 auto& yd = y.derived_cast();
2978 auto& xd = x.derived_cast();
2979 decltype(diff(x)) dx;
2980
2981 std::size_t saxis = normalize_axis(yd.dimension(), axis);
2982
2983 if (xd.dimension() == 1)
2984 {
2985 dx = diff(x);
2986 typename std::decay_t<decltype(yd)>::shape_type shape;
2987 resize_container(shape, yd.dimension());
2988 std::fill(shape.begin(), shape.end(), 1);
2989 shape[saxis] = dx.shape()[0];
2990 dx.reshape(shape);
2991 }
2992 else
2993 {
2994 dx = diff(x, 1, axis);
2995 }
2996
2997 xstrided_slice_vector slice1(yd.dimension(), all());
2998 xstrided_slice_vector slice2(yd.dimension(), all());
2999 slice1[saxis] = range(1, xnone());
3000 slice2[saxis] = range(xnone(), yd.shape()[saxis] - 1);
3001
3002 auto trap = dx * (strided_view(yd, slice1) + strided_view(yd, slice2)) * 0.5;
3003
3004 return eval(sum(trap, {saxis}));
3005 }
3006
3019 template <class E1, class E2, class E3, typename T>
3020 inline auto interp(const E1& x, const E2& xp, const E3& fp, T left, T right)
3021 {
3022 using size_type = common_size_type_t<E1, E2, E3>;
3023 using value_type = typename E3::value_type;
3024
3025 // basic checks
3026 XTENSOR_ASSERT(xp.dimension() == 1);
3027 XTENSOR_ASSERT(std::is_sorted(x.cbegin(), x.cend()));
3028 XTENSOR_ASSERT(std::is_sorted(xp.cbegin(), xp.cend()));
3029
3030 // allocate output
3031 auto f = xtensor<value_type, 1>::from_shape(x.shape());
3032
3033 // counter in "x": from left
3034 size_type i = 0;
3035
3036 // fill f[i] for x[i] <= xp[0]
3037 for (; i < x.size(); ++i)
3038 {
3039 if (x[i] > xp[0])
3040 {
3041 break;
3042 }
3043 f[i] = static_cast<value_type>(left);
3044 }
3045
3046 // counter in "x": from right
3047 // (index counts one right, to terminate the reverse loop, without risking being negative)
3048 size_type imax = x.size();
3049
3050 // fill f[i] for x[-1] >= xp[-1]
3051 for (; imax > 0; --imax)
3052 {
3053 if (x[imax - 1] < xp[xp.size() - 1])
3054 {
3055 break;
3056 }
3057 f[imax - 1] = static_cast<value_type>(right);
3058 }
3059
3060 // catch edge case: all entries are "right"
3061 if (imax == 0)
3062 {
3063 return f;
3064 }
3065
3066 // set "imax" as actual index
3067 // (counted one right, see above)
3068 --imax;
3069
3070 // counter in "xp"
3071 size_type ip = 1;
3072
3073 // fill f[i] for the interior
3074 for (; i <= imax; ++i)
3075 {
3076 // - search next value in "xp"
3077 while (x[i] > xp[ip])
3078 {
3079 ++ip;
3080 }
3081 // - distances as doubles
3082 double dfp = static_cast<double>(fp[ip] - fp[ip - 1]);
3083 double dxp = static_cast<double>(xp[ip] - xp[ip - 1]);
3084 double dx = static_cast<double>(x[i] - xp[ip - 1]);
3085 // - interpolate
3086 f[i] = fp[ip - 1] + static_cast<value_type>(dfp / dxp * dx);
3087 }
3088
3089 return f;
3090 }
3091
3092 namespace detail
3093 {
3094 template <class E1, class E2>
3095 auto calculate_discontinuity(E1&& discontinuity, E2&&)
3096 {
3097 return discontinuity;
3098 }
3099
3100 template <class E2>
3101 auto calculate_discontinuity(xt::placeholders::xtuph, E2&& period)
3102 {
3103 return 0.5 * period;
3104 }
3105
3106 template <class E1, class E2>
3107 auto
3108 calculate_interval(E2&& period, typename std::enable_if<std::is_integral<E1>::value, E1>::type* = 0)
3109 {
3110 auto interval_high = 0.5 * period;
3111 uint64_t remainder = static_cast<uint64_t>(period) % 2;
3112 auto boundary_ambiguous = (remainder == 0);
3113 return std::make_tuple(interval_high, boundary_ambiguous);
3114 }
3115
3116 template <class E1, class E2>
3117 auto
3118 calculate_interval(E2&& period, typename std::enable_if<std::is_floating_point<E1>::value, E1>::type* = 0)
3119 {
3120 auto interval_high = 0.5 * period;
3121 auto boundary_ambiguous = true;
3122 return std::make_tuple(interval_high, boundary_ambiguous);
3123 }
3124 }
3125
3139 template <class E1, class E2 = xt::placeholders::xtuph, class E3 = double>
3140 inline auto unwrap(
3141 E1&& p,
3142 E2 discontinuity = xnone(),
3143 std::ptrdiff_t axis = -1,
3144 E3 period = 2.0 * xt::numeric_constants<double>::PI
3145 )
3146 {
3147 auto discont = detail::calculate_discontinuity(discontinuity, period);
3148 using value_type = typename std::decay_t<E1>::value_type;
3149 std::size_t saxis = normalize_axis(p.dimension(), axis);
3150 auto dd = diff(p, 1, axis);
3151 xstrided_slice_vector slice(p.dimension(), all());
3152 slice[saxis] = range(1, xnone());
3153 auto interval_tuple = detail::calculate_interval<value_type>(period);
3154 auto interval_high = std::get<0>(interval_tuple);
3155 auto boundary_ambiguous = std::get<1>(interval_tuple);
3156 auto interval_low = -interval_high;
3157 auto ddmod = xt::eval(xt::fmod(xt::fmod(dd - interval_low, period) + period, period) + interval_low);
3158 if (boundary_ambiguous)
3159 {
3160 // for `mask = (abs(dd) == period/2)`, the above line made
3161 //`ddmod[mask] == -period/2`. correct these such that
3162 //`ddmod[mask] == sign(dd[mask])*period/2`.
3163 auto boolmap = xt::equal(ddmod, interval_low) && (xt::greater(dd, 0.0));
3164 ddmod = xt::where(boolmap, interval_high, ddmod);
3165 }
3166 auto ph_correct = xt::eval(ddmod - dd);
3167 ph_correct = xt::where(xt::abs(dd) < discont, 0.0, ph_correct);
3168 E1 up(p);
3169 strided_view(up, slice) = strided_view(p, slice)
3170 + xt::cumsum(ph_correct, static_cast<std::ptrdiff_t>(saxis));
3171 return up;
3172 }
3173
3184 template <class E1, class E2, class E3>
3185 inline auto interp(const E1& x, const E2& xp, const E3& fp)
3186 {
3187 return interp(x, xp, fp, fp[0], fp[fp.size() - 1]);
3188 }
3189
3196 template <class E1>
3197 inline auto cov(const E1& x, const E1& y = E1())
3198 {
3199 using value_type = typename E1::value_type;
3200
3201 if (y.dimension() == 0)
3202 {
3203 auto s = x.shape();
3204 using size_type = std::decay_t<decltype(s[0])>;
3205 if (x.dimension() == 1)
3206 {
3207 auto covar = eval(zeros<value_type>({1, 1}));
3208 auto x_norm = x - eval(mean(x));
3209 covar(0, 0) = std::inner_product(x_norm.begin(), x_norm.end(), x_norm.begin(), 0.0)
3210 / value_type(s[0] - 1);
3211 return covar;
3212 }
3213
3214 XTENSOR_ASSERT(x.dimension() == 2);
3215
3216 auto covar = eval(zeros<value_type>({s[0], s[0]}));
3217 auto m = eval(mean(x, {1}));
3218 m.reshape({m.shape()[0], 1});
3219 auto x_norm = x - m;
3220 for (size_type i = 0; i < s[0]; i++)
3221 {
3222 auto xi = strided_view(x_norm, {range(i, i + 1), all()});
3223 for (size_type j = i; j < s[0]; j++)
3224 {
3225 auto xj = strided_view(x_norm, {range(j, j + 1), all()});
3226 covar(j, i) = std::inner_product(xi.begin(), xi.end(), xj.begin(), 0.0)
3227 / value_type(s[1] - 1);
3228 }
3229 }
3230 return eval(covar + transpose(covar) - diag(diagonal(covar)));
3231 }
3232 else
3233 {
3234 return cov(eval(stack(xtuple(x, y))));
3235 }
3236 }
3237
3238 /*
3239 * convolution mode placeholders for selecting the algorithm
3240 * used in computing a 1D convolution.
3241 * Same as NumPy's mode parameter.
3242 */
3243 namespace convolve_mode
3244 {
3245 struct valid
3246 {
3247 };
3248
3249 struct full
3250 {
3251 };
3252 }
3253
3254 namespace detail
3255 {
3256 template <class E1, class E2>
3257 inline auto convolve_impl(E1&& e1, E2&& e2, convolve_mode::valid)
3258 {
3259 using value_type = typename std::decay<E1>::type::value_type;
3260
3261 const std::size_t na = e1.size();
3262 const std::size_t nv = e2.size();
3263 const std::size_t n = na - nv + 1;
3264 xt::xtensor<value_type, 1> out = xt::zeros<value_type>({n});
3265 for (std::size_t i = 0; i < n; i++)
3266 {
3267 for (std::size_t j = 0; j < nv; j++)
3268 {
3269 out(i) += e1(j) * e2(j + i);
3270 }
3271 }
3272 return out;
3273 }
3274
3275 template <class E1, class E2>
3276 inline auto convolve_impl(E1&& e1, E2&& e2, convolve_mode::full)
3277 {
3278 using value_type = typename std::decay<E1>::type::value_type;
3279
3280 const std::size_t na = e1.size();
3281 const std::size_t nv = e2.size();
3282 const std::size_t n = na + nv - 1;
3283 xt::xtensor<value_type, 1> out = xt::zeros<value_type>({n});
3284 for (std::size_t i = 0; i < n; i++)
3285 {
3286 const std::size_t jmn = (i >= nv - 1) ? i - (nv - 1) : 0;
3287 const std::size_t jmx = (i < na - 1) ? i : na - 1;
3288 for (std::size_t j = jmn; j <= jmx; ++j)
3289 {
3290 out(i) += e1(j) * e2(i - j);
3291 }
3292 }
3293 return out;
3294 }
3295 }
3296
3297 /*
3298 * @brief computes the 1D convolution between two 1D expressions
3299 *
3300 * @param a 1D expression
3301 * @param v 1D expression
3302 * @param mode placeholder Select algorithm #convolve_mode
3303 *
3304 * @detail the algorithm convolves a with v and will incur a copy overhead
3305 * should v be longer than a.
3306 */
3307 template <class E1, class E2, class E3>
3308 inline auto convolve(E1&& a, E2&& v, E3 mode)
3309 {
3310 if (a.dimension() != 1 || v.dimension() != 1)
3311 {
3312 XTENSOR_THROW(std::runtime_error, "Invalid dimentions convolution arguments must be 1D expressions");
3313 }
3314
3315 XTENSOR_ASSERT(a.size() > 0 && v.size() > 0);
3316
3317 // swap them so a is always the longest one
3318 if (a.size() < v.size())
3319 {
3320 return detail::convolve_impl(std::forward<E2>(v), std::forward<E1>(a), mode);
3321 }
3322 else
3323 {
3324 return detail::convolve_impl(std::forward<E1>(a), std::forward<E2>(v), mode);
3325 }
3326 }
3327}
3328
3329
3330#endif
xt::cumprod
auto cumprod(E &&e, std::ptrdiff_t axis)
Cumulative product.
Definition xmath.hpp:2321
xt::cumsum
auto cumsum(E &&e, std::ptrdiff_t axis)
Cumulative sum.
Definition xmath.hpp:2286
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:3020
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:3140
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:1813
xt::isfinite
auto isfinite(E &&e) noexcept -> detail::xfunction_type_t< math::isfinite_fun, E >
finite value check
Definition xmath.hpp:1678
xt::isnan
auto isnan(E &&e) noexcept -> detail::xfunction_type_t< math::isnan_fun, E >
NaN check.
Definition xmath.hpp:1708
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:1790
xt::isinf
auto isinf(E &&e) noexcept -> detail::xfunction_type_t< math::isinf_fun, E >
infinity check
Definition xmath.hpp:1693
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:690
xt::equal
auto equal(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< detail::equal_to, E1, E2 >
Element-wise equality.
Definition xoperation.hpp:674
xt::greater
auto greater(E1 &&e1, E2 &&e2) noexcept -> decltype(std::forward< E1 >(e1) > std::forward< E2 >(e2))
Greater than.
Definition xoperation.hpp:741
xt::lgamma
auto lgamma(E &&e) noexcept -> detail::xfunction_type_t< math::lgamma_fun, E >
Natural logarithm of the gamma function.
Definition xmath.hpp:1554
xt::erfc
auto erfc(E &&e) noexcept -> detail::xfunction_type_t< math::erfc_fun, E >
Complementary error function.
Definition xmath.hpp:1524
xt::erf
auto erf(E &&e) noexcept -> detail::xfunction_type_t< math::erf_fun, E >
Error function.
Definition xmath.hpp:1509
xt::tgamma
auto tgamma(E &&e) noexcept -> detail::xfunction_type_t< math::tgamma_fun, E >
Gamma function.
Definition xmath.hpp:1539
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:1456
xt::tanh
auto tanh(E &&e) noexcept -> detail::xfunction_type_t< math::tanh_fun, E >
Hyperbolic tangent function.
Definition xmath.hpp:1441
xt::cosh
auto cosh(E &&e) noexcept -> detail::xfunction_type_t< math::cosh_fun, E >
Hyperbolic cosine function.
Definition xmath.hpp:1426
xt::sinh
auto sinh(E &&e) noexcept -> detail::xfunction_type_t< math::sinh_fun, E >
Hyperbolic sine function.
Definition xmath.hpp:1411
xt::acosh
auto acosh(E &&e) noexcept -> detail::xfunction_type_t< math::acosh_fun, E >
Inverse hyperbolic cosine function.
Definition xmath.hpp:1471
xt::atanh
auto atanh(E &&e) noexcept -> detail::xfunction_type_t< math::atanh_fun, E >
Inverse hyperbolic tangent function.
Definition xmath.hpp:1486
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:777
xt::nanmax
auto nanmax(E &&e, X &&axes, EVS es=EVS())
Maximum element along given axes, ignoring NaNs.
Definition xmath.hpp:2480
xt::nancumsum
auto nancumsum(E &&e, std::ptrdiff_t axis)
Cumulative sum, replacing nan with 0.
Definition xmath.hpp:2630
xt::nancumprod
auto nancumprod(E &&e, std::ptrdiff_t axis)
Cumulative product, replacing nan with 1.
Definition xmath.hpp:2665
xt::nanmean
auto nanmean(E &&e, X &&axes, EVS es=EVS())
Mean of elements over given axes, excluding NaNs.
Definition xmath.hpp:2749
xt::nanmin
auto nanmin(E &&e, X &&axes, EVS es=EVS())
Minimum element over given axes, ignoring NaNs.
Definition xmath.hpp:2465
xt::nanprod
auto nanprod(E &&e, X &&axes, EVS es=EVS())
Product of elements over given axes, replacing NaN with 1.
Definition xmath.hpp:2514
xt::nansum
auto nansum(E &&e, X &&axes, EVS es=EVS())
Sum of elements over given axes, replacing NaN with 0.
Definition xmath.hpp:2497
xt::nan_to_num
auto nan_to_num(E &&e)
Convert nan or +/- inf to numbers.
Definition xmath.hpp:2447
xt::ceil
auto ceil(E &&e) noexcept -> detail::xfunction_type_t< math::ceil_fun, E >
ceil function.
Definition xmath.hpp:1577
xt::trunc
auto trunc(E &&e) noexcept -> detail::xfunction_type_t< math::trunc_fun, E >
trunc function.
Definition xmath.hpp:1607
xt::nearbyint
auto nearbyint(E &&e) noexcept -> detail::xfunction_type_t< math::nearbyint_fun, E >
nearbyint function.
Definition xmath.hpp:1639
xt::floor
auto floor(E &&e) noexcept -> detail::xfunction_type_t< math::floor_fun, E >
floor function.
Definition xmath.hpp:1592
xt::round
auto round(E &&e) noexcept -> detail::xfunction_type_t< math::round_fun, E >
round function.
Definition xmath.hpp:1623
xt::rint
auto rint(E &&e) noexcept -> detail::xfunction_type_t< math::rint_fun, E >
rint function.
Definition xmath.hpp:1655
xt::cube
auto cube(E1 &&e1) noexcept
Cube power function, equivalent to e1 * e1 * e1.
Definition xmath.hpp:1151
xt::sqrt
auto sqrt(E &&e) noexcept -> detail::xfunction_type_t< math::sqrt_fun, E >
Square root function.
Definition xmath.hpp:1239
xt::square
auto square(E1 &&e1) noexcept
Square power function, equivalent to e1 * e1.
Definition xmath.hpp:1128
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:1272
xt::cbrt
auto cbrt(E &&e) noexcept -> detail::xfunction_type_t< math::cbrt_fun, E >
Cubic root function.
Definition xmath.hpp:1254
xt::sum
auto sum(E &&e, X &&axes, EVS es=EVS())
Sum of elements over given axes.
Definition xmath.hpp:1841
xt::prod
auto prod(E &&e, X &&axes, EVS es=EVS())
Product of elements over given axes.
Definition xmath.hpp:1861
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:2948
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:2911
xt::minmax
auto minmax(E &&e, EVS es=EVS())
Minimum and maximum among the elements of an array or expression.
Definition xmath.hpp:2234
xt::average
auto average(E &&e, W &&weights, X &&axes, EVS ev=EVS())
Average of elements over given axes using weights.
Definition xmath.hpp:1976
xt::mean
auto mean(E &&e, X &&axes, EVS es=EVS())
Mean of elements over given axes.
Definition xmath.hpp:1935
xt::atan
auto atan(E &&e) noexcept -> detail::xfunction_type_t< math::atan_fun, E >
Arctangent function.
Definition xmath.hpp:1370
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:1388
xt::asin
auto asin(E &&e) noexcept -> detail::xfunction_type_t< math::asin_fun, E >
Arcsine function.
Definition xmath.hpp:1340
xt::cos
auto cos(E &&e) noexcept -> detail::xfunction_type_t< math::cos_fun, E >
Cosine function.
Definition xmath.hpp:1310
xt::sin
auto sin(E &&e) noexcept -> detail::xfunction_type_t< math::sin_fun, E >
Sine function.
Definition xmath.hpp:1295
xt::tan
auto tan(E &&e) noexcept -> detail::xfunction_type_t< math::tan_fun, E >
Tangent function.
Definition xmath.hpp:1325
xt::acos
auto acos(E &&e) noexcept -> detail::xfunction_type_t< math::acos_fun, E >
Arccosine function.
Definition xmath.hpp:1355
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:237
xt::same_shape
bool same_shape(const S1 &s1, const S2 &s2) noexcept
Check if two objects have the same shape.
Definition xshape.hpp:112
xt
standard mathematical functions for xexpressions
Definition xaccessible.hpp:18
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:818
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:234
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:1002
xt::layout_type
layout_type
Definition xlayout.hpp:24
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:343
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::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:1150
xt::make_xshared
xshared_expression< E > make_xshared(xexpression< E > &&expr)
Helper function to create shared expression from any xexpression.
Definition xexpression.hpp:730
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:728
xt::diag
auto diag(E &&arr, int k=0)
xexpression with values of arr on the diagonal, zeroes otherwise
Definition xbuilder.hpp:1202
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:3197
xt::convolve_mode::full
Definition xmath.hpp:3250
xt::convolve_mode::valid
Definition xmath.hpp:3246
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
xt::xaccumulator_functor
Definition xaccumulator.hpp:43