xtensor
Loading...
Searching...
No Matches
xnorm.hpp
1/***************************************************************************
2 * Copyright (c) Johan Mabille, Sylvain Corlay and Wolf Vollprecht *
3 * Copyright (c) Ullrich Koethe
4 * Copyright (c) QuantStack *
5 * *
6 * Distributed under the terms of the BSD 3-Clause License. *
7 * *
8 * The full license is in the file LICENSE, distributed with this software. *
9 ****************************************************************************/
10
11#ifndef XTENSOR_NORM_HPP
12#define XTENSOR_NORM_HPP
13
14#include <cmath>
15// std::abs(int) prior to C++ 17
16#include <complex>
17#include <cstdlib>
18
19#include <xtl/xtype_traits.hpp>
20
21#include "xmath.hpp"
22#include "xoperation.hpp"
23#include "xutils.hpp"
24
25namespace xt
26{
27 /********************************************
28 * type inference for norm and squared norm *
29 ********************************************/
30
31 template <class T>
32 struct norm_type;
33
34 template <class T>
35 struct squared_norm_type;
36
37 namespace traits_detail
38 {
39
40 template <class T, bool scalar = xtl::is_arithmetic<T>::value>
42
43 template <class T>
45 {
46 static const bool value = false;
47 using norm_type = void*;
48 using squared_norm_type = void*;
49 };
50
51 template <class T>
53 {
54 static const bool value = true;
55 using norm_type = xtl::promote_type_t<T>;
56 using squared_norm_type = xtl::promote_type_t<T>;
57 };
58
59 template <class T, bool integral = xtl::is_integral<T>::value, bool floating = std::is_floating_point<T>::value>
61
62 template <>
63 struct norm_of_array_elements_impl<void*, false, false>
64 {
65 using norm_type = void*;
66 using squared_norm_type = void*;
67 };
68
69 template <class T>
71 {
72 using norm_type = typename norm_type<T>::type;
73 using squared_norm_type = typename squared_norm_type<T>::type;
74 };
75
76 template <class T>
78 {
79 static_assert(
80 !std::is_same<T, char>::value,
81 "'char' is not a numeric type, use 'signed char' or 'unsigned char'."
82 );
83
84 using norm_type = double;
86 };
87
88 template <class T>
90 {
91 using norm_type = double;
93 };
94
95 template <>
97 {
98 using norm_type = long double;
99 using squared_norm_type = long double;
100 };
101
102 template <class ARRAY>
104 {
105 static void* test(...);
106
107 template <class U>
108 static typename U::value_type test(U*, typename U::value_type* = 0);
109
110 using T = decltype(test(std::declval<ARRAY*>()));
111
112 static const bool value = !std::is_same<T, void*>::value;
113
114 using norm_type = typename norm_of_array_elements_impl<T>::norm_type;
115 using squared_norm_type = typename norm_of_array_elements_impl<T>::squared_norm_type;
116 };
117
118 template <class U>
120 {
121 using T = std::decay_t<U>;
122
123 static_assert(
124 !std::is_same<T, char>::value,
125 "'char' is not a numeric type, use 'signed char' or 'unsigned char'."
126 );
127
130
131 static const bool value = norm_of_scalar::value || norm_of_vector::value;
132
133 static_assert(value, "norm_type<T> are undefined for type U.");
134 };
135 } // namespace traits_detail
136
152 template <class T>
154 {
156
157 using type = typename std::conditional<
158 base_type::norm_of_vector::value,
159 typename base_type::norm_of_vector::norm_type,
160 typename base_type::norm_of_scalar::norm_type>::type;
161 };
162
166 template <class T>
167 using norm_type_t = typename norm_type<T>::type;
168
185 template <class T>
187 {
189
190 using type = typename std::conditional<
191 base_type::norm_of_vector::value,
192 typename base_type::norm_of_vector::squared_norm_type,
193 typename base_type::norm_of_scalar::squared_norm_type>::type;
194 };
195
199 template <class T>
200 using squared_norm_type_t = typename squared_norm_type<T>::type;
201
202 /*************************************
203 * norm functions for built-in types *
204 *************************************/
205
207#define XTENSOR_DEFINE_SIGNED_NORMS(T) \
208 inline auto norm_lp(T t, double p) noexcept \
209 { \
210 using rt = decltype(std::abs(t)); \
211 return p == 0.0 ? static_cast<rt>(t != 0) : std::abs(t); \
212 } \
213 inline auto norm_lp_to_p(T t, double p) noexcept \
214 { \
215 using rt = xtl::real_promote_type_t<T>; \
216 return p == 0.0 ? static_cast<rt>(t != 0) \
217 : std::pow(static_cast<rt>(std::abs(t)), static_cast<rt>(p)); \
218 } \
219 inline std::size_t norm_l0(T t) noexcept \
220 { \
221 return (t != 0); \
222 } \
223 inline auto norm_l1(T t) noexcept \
224 { \
225 return std::abs(t); \
226 } \
227 inline auto norm_l2(T t) noexcept \
228 { \
229 return std::abs(t); \
230 } \
231 inline auto norm_linf(T t) noexcept \
232 { \
233 return std::abs(t); \
234 } \
235 inline auto norm_sq(T t) noexcept \
236 { \
237 return t * t; \
238 }
239
240 XTENSOR_DEFINE_SIGNED_NORMS(signed char)
247 XTENSOR_DEFINE_SIGNED_NORMS(long double)
248
249#undef XTENSOR_DEFINE_SIGNED_NORMS
250
251#define XTENSOR_DEFINE_UNSIGNED_NORMS(T) \
252 inline T norm_lp(T t, double p) noexcept \
253 { \
254 return p == 0.0 ? (t != 0) : t; \
255 } \
256 inline auto norm_lp_to_p(T t, double p) noexcept \
257 { \
258 using rt = xtl::real_promote_type_t<T>; \
259 return p == 0.0 ? static_cast<rt>(t != 0) : std::pow(static_cast<rt>(t), static_cast<rt>(p)); \
260 } \
261 inline T norm_l0(T t) noexcept \
262 { \
263 return t != 0 ? 1 : 0; \
264 } \
265 inline T norm_l1(T t) noexcept \
266 { \
267 return t; \
268 } \
269 inline T norm_l2(T t) noexcept \
270 { \
271 return t; \
272 } \
273 inline T norm_linf(T t) noexcept \
274 { \
275 return t; \
276 } \
277 inline auto norm_sq(T t) noexcept \
278 { \
279 return t * t; \
280 }
281
282 XTENSOR_DEFINE_UNSIGNED_NORMS(unsigned char)
283 XTENSOR_DEFINE_UNSIGNED_NORMS(unsigned short)
285 XTENSOR_DEFINE_UNSIGNED_NORMS(unsigned long)
286 XTENSOR_DEFINE_UNSIGNED_NORMS(unsigned long long)
287
288#undef XTENSOR_DEFINE_UNSIGNED_NORMS
289
290 /***********************************
291 * norm functions for std::complex *
292 ***********************************/
293
298 template <class T>
299 inline uint64_t norm_l0(const std::complex<T>& t) noexcept
300 {
301 return t.real() != 0 || t.imag() != 0;
302 }
303
307 template <class T>
308 inline auto norm_l1(const std::complex<T>& t) noexcept
309 {
310 return std::abs(t.real()) + std::abs(t.imag());
311 }
312
317 template <class T>
318 inline auto norm_l2(const std::complex<T>& t) noexcept
319 {
320 return std::abs(t);
321 }
322
328 template <class T>
329 inline auto norm_sq(const std::complex<T>& t) noexcept
330 {
331 // Does not use std::norm since it returns a std::complex on OSX
332 return t.real() * t.real() + t.imag() * t.imag();
333 }
334
338 template <class T>
339 inline auto norm_linf(const std::complex<T>& t) noexcept
340 {
341 return (std::max)(std::abs(t.real()), std::abs(t.imag()));
342 }
343
347 template <class T>
348 inline auto norm_lp_to_p(const std::complex<T>& t, double p) noexcept
349 {
350 using rt = decltype(std::pow(std::abs(t.real()), static_cast<T>(p)));
351 return p == 0 ? static_cast<rt>(t.real() != 0 || t.imag() != 0)
352 : std::pow(std::abs(t.real()), static_cast<T>(p))
353 + std::pow(std::abs(t.imag()), static_cast<T>(p));
354 }
355
359 template <class T>
360 inline auto norm_lp(const std::complex<T>& t, double p) noexcept
361 {
362 return p == 0 ? norm_lp_to_p(t, p) : std::pow(norm_lp_to_p(t, p), 1.0 / p);
363 }
364
365 /***********************************
366 * norm functions for xexpressions *
367 ***********************************/
368
369#define XTENSOR_NORM_FUNCTION_AXES(NAME) \
370 template <class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS> \
371 inline auto NAME(E&& e, const I(&axes)[N], EVS es = EVS()) noexcept \
372 { \
373 using axes_type = std::array<typename std::decay_t<E>::size_type, N>; \
374 return NAME(std::forward<E>(e), xtl::forward_sequence<axes_type, decltype(axes)>(axes), es); \
375 }
376
377 namespace detail
378 {
379 template <class T>
380 struct norm_value_type
381 {
382 using type = T;
383 };
384
385 template <class T>
386 struct norm_value_type<std::complex<T>>
387 {
388 using type = T;
389 };
390
391 template <class T>
392 using norm_value_type_t = typename norm_value_type<T>::type;
393 }
394
395#define XTENSOR_EMPTY
396#define XTENSOR_COMMA ,
397#define XTENSOR_NORM_FUNCTION(NAME, RESULT_TYPE, REDUCE_EXPR, REDUCE_OP, MERGE_FUNC) \
398 template <class E, class X, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(xtl::negation<is_reducer_options<X>>)> \
399 inline auto NAME(E&& e, X&& axes, EVS es = EVS()) noexcept \
400 { \
401 using value_type = typename std::decay_t<E>::value_type; \
402 using result_type = detail::norm_value_type_t<RESULT_TYPE>; \
403 \
404 auto reduce_func = [](result_type const& r, value_type const& v) \
405 { \
406 return REDUCE_EXPR(r REDUCE_OP NAME(v)); \
407 }; \
408 \
409 return xt::reduce( \
410 make_xreducer_functor(std::move(reduce_func), const_value<result_type>(0), MERGE_FUNC<result_type>()), \
411 std::forward<E>(e), \
412 std::forward<X>(axes), \
413 es \
414 ); \
415 } \
416 \
417 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_xexpression<E>)> \
418 inline auto NAME(E&& e, EVS es = EVS()) noexcept \
419 { \
420 return NAME(std::forward<E>(e), arange(e.dimension()), es); \
421 } \
422 XTENSOR_NORM_FUNCTION_AXES(NAME)
423
424 XTENSOR_NORM_FUNCTION(norm_l0, unsigned long long, XTENSOR_EMPTY, +, std::plus)
425 XTENSOR_NORM_FUNCTION(norm_l1, xtl::big_promote_type_t<value_type>, XTENSOR_EMPTY, +, std::plus)
426 XTENSOR_NORM_FUNCTION(norm_sq, xtl::big_promote_type_t<value_type>, XTENSOR_EMPTY, +, std::plus)
427 XTENSOR_NORM_FUNCTION(
428 norm_linf,
429 decltype(norm_linf(std::declval<value_type>())),
430 (std::max<result_type>),
431 XTENSOR_COMMA,
432 math::maximum
433 )
434
435#undef XTENSOR_EMPTY
436#undef XTENSOR_COMMA
437#undef XTENSOR_NORM_FUNCTION
438#undef XTENSOR_NORM_FUNCTION_AXES
440
452 template <class E, class X, class EVS, class>
453 auto norm_l0(E&& e, X&& axes, EVS es) noexcept;
454
467 template <class E, class X, class EVS, class>
468 auto norm_l1(E&& e, X&& axes, EVS es) noexcept;
469
482 template <class E, class X, class EVS, class>
483 auto norm_sq(E&& e, X&& axes, EVS es) noexcept;
484
493 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_xexpression<E>)>
494 inline auto norm_l2(E&& e, EVS es = EVS()) noexcept
495 {
496 using std::sqrt;
497 return sqrt(norm_sq(std::forward<E>(e), es));
498 }
499
511 template <
512 class E,
513 class X,
514 class EVS = DEFAULT_STRATEGY_REDUCERS,
515 XTL_REQUIRES(is_xexpression<E>, xtl::negation<is_reducer_options<X>>)>
516 inline auto norm_l2(E&& e, X&& axes, EVS es = EVS()) noexcept
517 {
518 return sqrt(norm_sq(std::forward<E>(e), std::forward<X>(axes), es));
519 }
520
521 template <class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
522 inline auto norm_l2(E&& e, const I (&axes)[N], EVS es = EVS()) noexcept
523 {
524 using axes_type = std::array<typename std::decay_t<E>::size_type, N>;
525 return sqrt(norm_sq(std::forward<E>(e), xtl::forward_sequence<axes_type, decltype(axes)>(axes), es));
526 }
527
540 template <class E, class X, class EVS, class>
541 auto norm_linf(E&& e, X&& axes, EVS es) noexcept;
542
556 template <class E, class X, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(xtl::negation<is_reducer_options<X>>)>
557 inline auto norm_lp_to_p(E&& e, double p, X&& axes, EVS es = EVS()) noexcept
558 {
559 using value_type = typename std::decay_t<E>::value_type;
560 using result_type = norm_type_t<std::decay_t<E>>;
561
562 auto reduce_func = [p](const result_type& r, const value_type& v)
563 {
564 return r + norm_lp_to_p(v, p);
565 };
566 return xt::reduce(
567 make_xreducer_functor(std::move(reduce_func), xt::const_value<result_type>(0), std::plus<result_type>()),
568 std::forward<E>(e),
569 std::forward<X>(axes),
570 es
571 );
572 }
573
574 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_xexpression<E>)>
575 inline auto norm_lp_to_p(E&& e, double p, EVS es = EVS()) noexcept
576 {
577 return norm_lp_to_p(std::forward<E>(e), p, arange(e.dimension()), es);
578 }
579
580 template <class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
581 inline auto norm_lp_to_p(E&& e, double p, const I (&axes)[N], EVS es = EVS()) noexcept
582 {
583 using axes_type = std::array<typename std::decay_t<E>::size_type, N>;
584 return norm_lp_to_p(std::forward<E>(e), p, xtl::forward_sequence<axes_type, decltype(axes)>(axes), es);
585 }
586
600 template <class E, class X, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(xtl::negation<is_reducer_options<X>>)>
601 inline auto norm_lp(E&& e, double p, X&& axes, EVS es = EVS())
602 {
603 XTENSOR_PRECONDITION(p != 0, "norm_lp(): p must be nonzero, use norm_l0() instead.");
604 return pow(norm_lp_to_p(std::forward<E>(e), p, std::forward<X>(axes), es), 1.0 / p);
605 }
606
607 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_xexpression<E>)>
608 inline auto norm_lp(E&& e, double p, EVS es = EVS())
609 {
610 return norm_lp(std::forward<E>(e), p, arange(e.dimension()), es);
611 }
612
613 template <class E, class I, std::size_t N, class EVS = DEFAULT_STRATEGY_REDUCERS>
614 inline auto norm_lp(E&& e, double p, const I (&axes)[N], EVS es = EVS())
615 {
616 using axes_type = std::array<typename std::decay_t<E>::size_type, N>;
617 return norm_lp(std::forward<E>(e), p, xtl::forward_sequence<axes_type, decltype(axes)>(axes), es);
618 }
619
629 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_xexpression<E>)>
630 inline auto norm_induced_l1(E&& e, EVS es = EVS())
631 {
632 XTENSOR_PRECONDITION(
633 e.dimension() == 2,
634 "norm_induced_l1(): only applicable to matrices (e.dimension() must be 2)."
635 );
636 return norm_linf(norm_l1(std::forward<E>(e), {0}, es), es);
637 }
638
649 template <class E, class EVS = DEFAULT_STRATEGY_REDUCERS, XTL_REQUIRES(is_xexpression<E>)>
650 inline auto norm_induced_linf(E&& e, EVS es = EVS())
651 {
652 XTENSOR_PRECONDITION(
653 e.dimension() == 2,
654 "norm_induced_linf(): only applicable to matrices (e.dimension() must be 2)."
655 );
656 return norm_linf(norm_l1(std::forward<E>(e), {1}, es), es);
657 }
658
659} // namespace xt
660
661#endif
auto sqrt(E &&e) noexcept -> detail::xfunction_type_t< math::sqrt_fun, E >
Square root function.
Definition xmath.hpp:1239
auto pow(E1 &&e1, E2 &&e2) noexcept -> detail::xfunction_type_t< math::pow_fun, E1, E2 >
Power function.
Definition xmath.hpp:1015
auto norm_sq(E &&e, X &&axes, EVS es) noexcept
Squared L2 norm of an array-like argument over given axes.
auto norm_induced_linf(E &&e, EVS es=EVS())
Induced L-infinity norm of a matrix.
Definition xnorm.hpp:650
auto norm_lp(E &&e, double p, X &&axes, EVS es=EVS())
Lp norm of an array-like argument over given axes.
Definition xnorm.hpp:601
auto norm_l2(E &&e, EVS es=EVS()) noexcept
L2 norm of a scalar or array-like argument.
Definition xnorm.hpp:494
auto norm_induced_l1(E &&e, EVS es=EVS())
Induced L1 norm of a matrix.
Definition xnorm.hpp:630
auto norm_l1(E &&e, X &&axes, EVS es) noexcept
L1 norm of an array-like argument over given axes.
auto norm_lp_to_p(E &&e, double p, X &&axes, EVS es=EVS()) noexcept
p-th power of the Lp norm of an array-like argument over given axes.
Definition xnorm.hpp:557
auto norm_l0(E &&e, X &&axes, EVS es) noexcept
L0 (count) pseudo-norm of an array-like argument over given axes.
auto norm_linf(E &&e, X &&axes, EVS es) noexcept
Infinity (maximum) norm of an array-like argument over given axes.
standard mathematical functions for xexpressions
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
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.
typename squared_norm_type< T >::type squared_norm_type_t
Abbreviation of 'typename squared_norm_type<T>::type'.
Definition xnorm.hpp:200
typename norm_type< T >::type norm_type_t
Abbreviation of 'typename norm_type<T>::type'.
Definition xnorm.hpp:167
Traits class for the result type of the norm_l2() function.
Definition xnorm.hpp:154
Traits class for the result type of the norm_sq() function.
Definition xnorm.hpp:187