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