diff options
Diffstat (limited to 'Eigen/src/Core/VectorwiseOp.h')
| -rw-r--r-- | Eigen/src/Core/VectorwiseOp.h | 229 |
1 files changed, 159 insertions, 70 deletions
diff --git a/Eigen/src/Core/VectorwiseOp.h b/Eigen/src/Core/VectorwiseOp.h index 4fe267e9f..870f4f1e4 100644 --- a/Eigen/src/Core/VectorwiseOp.h +++ b/Eigen/src/Core/VectorwiseOp.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2008-2019 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> // // This Source Code Form is subject to the terms of the Mozilla @@ -65,10 +65,10 @@ class PartialReduxExpr : public internal::dense_xpr_base< PartialReduxExpr<Matri explicit PartialReduxExpr(const MatrixType& mat, const MemberOp& func = MemberOp()) : m_matrix(mat), m_functor(func) {} - EIGEN_DEVICE_FUNC - Index rows() const { return (Direction==Vertical ? 1 : m_matrix.rows()); } - EIGEN_DEVICE_FUNC - Index cols() const { return (Direction==Horizontal ? 1 : m_matrix.cols()); } + EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR + Index rows() const EIGEN_NOEXCEPT { return (Direction==Vertical ? 1 : m_matrix.rows()); } + EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR + Index cols() const EIGEN_NOEXCEPT { return (Direction==Horizontal ? 1 : m_matrix.cols()); } EIGEN_DEVICE_FUNC typename MatrixType::Nested nestedExpression() const { return m_matrix; } @@ -81,39 +81,46 @@ class PartialReduxExpr : public internal::dense_xpr_base< PartialReduxExpr<Matri const MemberOp m_functor; }; -#define EIGEN_MEMBER_FUNCTOR(MEMBER,COST) \ - template <typename ResultType> \ - struct member_##MEMBER { \ - EIGEN_EMPTY_STRUCT_CTOR(member_##MEMBER) \ - typedef ResultType result_type; \ - template<typename Scalar, int Size> struct Cost \ - { enum { value = COST }; }; \ - template<typename XprType> \ - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \ - ResultType operator()(const XprType& mat) const \ - { return mat.MEMBER(); } \ +template<typename A,typename B> struct partial_redux_dummy_func; + +#define EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(MEMBER,COST,VECTORIZABLE,BINARYOP) \ + template <typename ResultType,typename Scalar> \ + struct member_##MEMBER { \ + EIGEN_EMPTY_STRUCT_CTOR(member_##MEMBER) \ + typedef ResultType result_type; \ + typedef BINARYOP<Scalar,Scalar> BinaryOp; \ + template<int Size> struct Cost { enum { value = COST }; }; \ + enum { Vectorizable = VECTORIZABLE }; \ + template<typename XprType> \ + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \ + ResultType operator()(const XprType& mat) const \ + { return mat.MEMBER(); } \ + BinaryOp binaryFunc() const { return BinaryOp(); } \ } +#define EIGEN_MEMBER_FUNCTOR(MEMBER,COST) \ + EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(MEMBER,COST,0,partial_redux_dummy_func) + namespace internal { -EIGEN_MEMBER_FUNCTOR(squaredNorm, Size * NumTraits<Scalar>::MulCost + (Size-1)*NumTraits<Scalar>::AddCost); EIGEN_MEMBER_FUNCTOR(norm, (Size+5) * NumTraits<Scalar>::MulCost + (Size-1)*NumTraits<Scalar>::AddCost); EIGEN_MEMBER_FUNCTOR(stableNorm, (Size+5) * NumTraits<Scalar>::MulCost + (Size-1)*NumTraits<Scalar>::AddCost); EIGEN_MEMBER_FUNCTOR(blueNorm, (Size+5) * NumTraits<Scalar>::MulCost + (Size-1)*NumTraits<Scalar>::AddCost); EIGEN_MEMBER_FUNCTOR(hypotNorm, (Size-1) * functor_traits<scalar_hypot_op<Scalar> >::Cost ); -EIGEN_MEMBER_FUNCTOR(sum, (Size-1)*NumTraits<Scalar>::AddCost); -EIGEN_MEMBER_FUNCTOR(mean, (Size-1)*NumTraits<Scalar>::AddCost + NumTraits<Scalar>::MulCost); -EIGEN_MEMBER_FUNCTOR(minCoeff, (Size-1)*NumTraits<Scalar>::AddCost); -EIGEN_MEMBER_FUNCTOR(maxCoeff, (Size-1)*NumTraits<Scalar>::AddCost); EIGEN_MEMBER_FUNCTOR(all, (Size-1)*NumTraits<Scalar>::AddCost); EIGEN_MEMBER_FUNCTOR(any, (Size-1)*NumTraits<Scalar>::AddCost); EIGEN_MEMBER_FUNCTOR(count, (Size-1)*NumTraits<Scalar>::AddCost); -EIGEN_MEMBER_FUNCTOR(prod, (Size-1)*NumTraits<Scalar>::MulCost); -template <int p, typename ResultType> +EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(sum, (Size-1)*NumTraits<Scalar>::AddCost, 1, internal::scalar_sum_op); +EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(minCoeff, (Size-1)*NumTraits<Scalar>::AddCost, 1, internal::scalar_min_op); +EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(maxCoeff, (Size-1)*NumTraits<Scalar>::AddCost, 1, internal::scalar_max_op); +EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(prod, (Size-1)*NumTraits<Scalar>::MulCost, 1, internal::scalar_product_op); + +template <int p, typename ResultType,typename Scalar> struct member_lpnorm { typedef ResultType result_type; - template<typename Scalar, int Size> struct Cost + enum { Vectorizable = 0 }; + template<int Size> struct Cost { enum { value = (Size+5) * NumTraits<Scalar>::MulCost + (Size-1)*NumTraits<Scalar>::AddCost }; }; EIGEN_DEVICE_FUNC member_lpnorm() {} template<typename XprType> @@ -121,17 +128,20 @@ struct member_lpnorm { { return mat.template lpNorm<p>(); } }; -template <typename BinaryOp, typename Scalar> +template <typename BinaryOpT, typename Scalar> struct member_redux { + typedef BinaryOpT BinaryOp; typedef typename result_of< BinaryOp(const Scalar&,const Scalar&) >::type result_type; - template<typename _Scalar, int Size> struct Cost - { enum { value = (Size-1) * functor_traits<BinaryOp>::Cost }; }; + + enum { Vectorizable = functor_traits<BinaryOp>::PacketAccess }; + template<int Size> struct Cost { enum { value = (Size-1) * functor_traits<BinaryOp>::Cost }; }; EIGEN_DEVICE_FUNC explicit member_redux(const BinaryOp func) : m_functor(func) {} template<typename Derived> EIGEN_DEVICE_FUNC inline result_type operator()(const DenseBase<Derived>& mat) const { return mat.redux(m_functor); } + const BinaryOp& binaryFunc() const { return m_functor; } const BinaryOp m_functor; }; } @@ -139,18 +149,38 @@ struct member_redux { /** \class VectorwiseOp * \ingroup Core_Module * - * \brief Pseudo expression providing partial reduction operations + * \brief Pseudo expression providing broadcasting and partial reduction operations * * \tparam ExpressionType the type of the object on which to do partial reductions - * \tparam Direction indicates the direction of the redux (#Vertical or #Horizontal) + * \tparam Direction indicates whether to operate on columns (#Vertical) or rows (#Horizontal) * - * This class represents a pseudo expression with partial reduction features. + * This class represents a pseudo expression with broadcasting and partial reduction features. * It is the return type of DenseBase::colwise() and DenseBase::rowwise() - * and most of the time this is the only way it is used. + * and most of the time this is the only way it is explicitly used. + * + * To understand the logic of rowwise/colwise expression, let's consider a generic case `A.colwise().foo()` + * where `foo` is any method of `VectorwiseOp`. This expression is equivalent to applying `foo()` to each + * column of `A` and then re-assemble the outputs in a matrix expression: + * \code [A.col(0).foo(), A.col(1).foo(), ..., A.col(A.cols()-1).foo()] \endcode * * Example: \include MatrixBase_colwise.cpp * Output: \verbinclude MatrixBase_colwise.out * + * The begin() and end() methods are obviously exceptions to the previous rule as they + * return STL-compatible begin/end iterators to the rows or columns of the nested expression. + * Typical use cases include for-range-loop and calls to STL algorithms: + * + * Example: \include MatrixBase_colwise_iterator_cxx11.cpp + * Output: \verbinclude MatrixBase_colwise_iterator_cxx11.out + * + * For a partial reduction on an empty input, some rules apply. + * For the sake of clarity, let's consider a vertical reduction: + * - If the number of columns is zero, then a 1x0 row-major vector expression is returned. + * - Otherwise, if the number of rows is zero, then + * - a row vector of zeros is returned for sum-like reductions (sum, squaredNorm, norm, etc.) + * - a row vector of ones is returned for a product reduction (e.g., <code>MatrixXd(n,0).colwise().prod()</code>) + * - an assert is triggered for all other reductions (minCoeff,maxCoeff,redux(bin_op)) + * * \sa DenseBase::colwise(), DenseBase::rowwise(), class PartialReduxExpr */ template<typename ExpressionType, int Direction> class VectorwiseOp @@ -163,11 +193,11 @@ template<typename ExpressionType, int Direction> class VectorwiseOp typedef typename internal::ref_selector<ExpressionType>::non_const_type ExpressionTypeNested; typedef typename internal::remove_all<ExpressionTypeNested>::type ExpressionTypeNestedCleaned; - template<template<typename _Scalar> class Functor, - typename Scalar_=Scalar> struct ReturnType + template<template<typename OutScalar,typename InputScalar> class Functor, + typename ReturnScalar=Scalar> struct ReturnType { typedef PartialReduxExpr<ExpressionType, - Functor<Scalar_>, + Functor<ReturnScalar,Scalar>, Direction > Type; }; @@ -187,23 +217,6 @@ template<typename ExpressionType, int Direction> class VectorwiseOp protected: - typedef typename internal::conditional<isVertical, - typename ExpressionType::ColXpr, - typename ExpressionType::RowXpr>::type SubVector; - /** \internal - * \returns the i-th subvector according to the \c Direction */ - EIGEN_DEVICE_FUNC - SubVector subVector(Index i) - { - return SubVector(m_matrix.derived(),i); - } - - /** \internal - * \returns the number of subvectors in the direction \c Direction */ - EIGEN_DEVICE_FUNC - Index subVectors() const - { return isVertical?m_matrix.cols():m_matrix.rows(); } - template<typename OtherDerived> struct ExtendedType { typedef Replicate<OtherDerived, isVertical ? 1 : ExpressionType::RowsAtCompileTime, @@ -258,42 +271,101 @@ template<typename ExpressionType, int Direction> class VectorwiseOp EIGEN_DEVICE_FUNC inline const ExpressionType& _expression() const { return m_matrix; } + #ifdef EIGEN_PARSED_BY_DOXYGEN + /** STL-like <a href="https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator">RandomAccessIterator</a> + * iterator type over the columns or rows as returned by the begin() and end() methods. + */ + random_access_iterator_type iterator; + /** This is the const version of iterator (aka read-only) */ + random_access_iterator_type const_iterator; + #else + typedef internal::subvector_stl_iterator<ExpressionType, DirectionType(Direction)> iterator; + typedef internal::subvector_stl_iterator<const ExpressionType, DirectionType(Direction)> const_iterator; + typedef internal::subvector_stl_reverse_iterator<ExpressionType, DirectionType(Direction)> reverse_iterator; + typedef internal::subvector_stl_reverse_iterator<const ExpressionType, DirectionType(Direction)> const_reverse_iterator; + #endif + + /** returns an iterator to the first row (rowwise) or column (colwise) of the nested expression. + * \sa end(), cbegin() + */ + iterator begin() { return iterator (m_matrix, 0); } + /** const version of begin() */ + const_iterator begin() const { return const_iterator(m_matrix, 0); } + /** const version of begin() */ + const_iterator cbegin() const { return const_iterator(m_matrix, 0); } + + /** returns a reverse iterator to the last row (rowwise) or column (colwise) of the nested expression. + * \sa rend(), crbegin() + */ + reverse_iterator rbegin() { return reverse_iterator (m_matrix, m_matrix.template subVectors<DirectionType(Direction)>()-1); } + /** const version of rbegin() */ + const_reverse_iterator rbegin() const { return const_reverse_iterator (m_matrix, m_matrix.template subVectors<DirectionType(Direction)>()-1); } + /** const version of rbegin() */ + const_reverse_iterator crbegin() const { return const_reverse_iterator (m_matrix, m_matrix.template subVectors<DirectionType(Direction)>()-1); } + + /** returns an iterator to the row (resp. column) following the last row (resp. column) of the nested expression + * \sa begin(), cend() + */ + iterator end() { return iterator (m_matrix, m_matrix.template subVectors<DirectionType(Direction)>()); } + /** const version of end() */ + const_iterator end() const { return const_iterator(m_matrix, m_matrix.template subVectors<DirectionType(Direction)>()); } + /** const version of end() */ + const_iterator cend() const { return const_iterator(m_matrix, m_matrix.template subVectors<DirectionType(Direction)>()); } + + /** returns a reverse iterator to the row (resp. column) before the first row (resp. column) of the nested expression + * \sa begin(), cend() + */ + reverse_iterator rend() { return reverse_iterator (m_matrix, -1); } + /** const version of rend() */ + const_reverse_iterator rend() const { return const_reverse_iterator (m_matrix, -1); } + /** const version of rend() */ + const_reverse_iterator crend() const { return const_reverse_iterator (m_matrix, -1); } + /** \returns a row or column vector expression of \c *this reduxed by \a func * * The template parameter \a BinaryOp is the type of the functor * of the custom redux operator. Note that func must be an associative operator. * + * \warning the size along the reduction direction must be strictly positive, + * otherwise an assertion is triggered. + * * \sa class VectorwiseOp, DenseBase::colwise(), DenseBase::rowwise() */ template<typename BinaryOp> EIGEN_DEVICE_FUNC const typename ReduxReturnType<BinaryOp>::Type redux(const BinaryOp& func = BinaryOp()) const - { return typename ReduxReturnType<BinaryOp>::Type(_expression(), internal::member_redux<BinaryOp,Scalar>(func)); } + { + eigen_assert(redux_length()>0 && "you are using an empty matrix"); + return typename ReduxReturnType<BinaryOp>::Type(_expression(), internal::member_redux<BinaryOp,Scalar>(func)); + } typedef typename ReturnType<internal::member_minCoeff>::Type MinCoeffReturnType; typedef typename ReturnType<internal::member_maxCoeff>::Type MaxCoeffReturnType; - typedef typename ReturnType<internal::member_squaredNorm,RealScalar>::Type SquaredNormReturnType; - typedef typename ReturnType<internal::member_norm,RealScalar>::Type NormReturnType; + typedef PartialReduxExpr<const CwiseUnaryOp<internal::scalar_abs2_op<Scalar>, const ExpressionTypeNestedCleaned>,internal::member_sum<RealScalar,RealScalar>,Direction> SquaredNormReturnType; + typedef CwiseUnaryOp<internal::scalar_sqrt_op<RealScalar>, const SquaredNormReturnType> NormReturnType; typedef typename ReturnType<internal::member_blueNorm,RealScalar>::Type BlueNormReturnType; typedef typename ReturnType<internal::member_stableNorm,RealScalar>::Type StableNormReturnType; typedef typename ReturnType<internal::member_hypotNorm,RealScalar>::Type HypotNormReturnType; typedef typename ReturnType<internal::member_sum>::Type SumReturnType; - typedef typename ReturnType<internal::member_mean>::Type MeanReturnType; + typedef EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(SumReturnType,Scalar,quotient) MeanReturnType; typedef typename ReturnType<internal::member_all>::Type AllReturnType; typedef typename ReturnType<internal::member_any>::Type AnyReturnType; - typedef PartialReduxExpr<ExpressionType, internal::member_count<Index>, Direction> CountReturnType; + typedef PartialReduxExpr<ExpressionType, internal::member_count<Index,Scalar>, Direction> CountReturnType; typedef typename ReturnType<internal::member_prod>::Type ProdReturnType; typedef Reverse<const ExpressionType, Direction> ConstReverseReturnType; typedef Reverse<ExpressionType, Direction> ReverseReturnType; template<int p> struct LpNormReturnType { - typedef PartialReduxExpr<ExpressionType, internal::member_lpnorm<p,RealScalar>,Direction> Type; + typedef PartialReduxExpr<ExpressionType, internal::member_lpnorm<p,RealScalar,Scalar>,Direction> Type; }; /** \returns a row (or column) vector expression of the smallest coefficient * of each column (or row) of the referenced expression. * + * \warning the size along the reduction direction must be strictly positive, + * otherwise an assertion is triggered. + * * \warning the result is undefined if \c *this contains NaN. * * Example: \include PartialRedux_minCoeff.cpp @@ -302,11 +374,17 @@ template<typename ExpressionType, int Direction> class VectorwiseOp * \sa DenseBase::minCoeff() */ EIGEN_DEVICE_FUNC const MinCoeffReturnType minCoeff() const - { return MinCoeffReturnType(_expression()); } + { + eigen_assert(redux_length()>0 && "you are using an empty matrix"); + return MinCoeffReturnType(_expression()); + } /** \returns a row (or column) vector expression of the largest coefficient * of each column (or row) of the referenced expression. * + * \warning the size along the reduction direction must be strictly positive, + * otherwise an assertion is triggered. + * * \warning the result is undefined if \c *this contains NaN. * * Example: \include PartialRedux_maxCoeff.cpp @@ -315,7 +393,10 @@ template<typename ExpressionType, int Direction> class VectorwiseOp * \sa DenseBase::maxCoeff() */ EIGEN_DEVICE_FUNC const MaxCoeffReturnType maxCoeff() const - { return MaxCoeffReturnType(_expression()); } + { + eigen_assert(redux_length()>0 && "you are using an empty matrix"); + return MaxCoeffReturnType(_expression()); + } /** \returns a row (or column) vector expression of the squared norm * of each column (or row) of the referenced expression. @@ -327,7 +408,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp * \sa DenseBase::squaredNorm() */ EIGEN_DEVICE_FUNC const SquaredNormReturnType squaredNorm() const - { return SquaredNormReturnType(_expression()); } + { return SquaredNormReturnType(m_matrix.cwiseAbs2()); } /** \returns a row (or column) vector expression of the norm * of each column (or row) of the referenced expression. @@ -339,7 +420,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp * \sa DenseBase::norm() */ EIGEN_DEVICE_FUNC const NormReturnType norm() const - { return NormReturnType(_expression()); } + { return NormReturnType(squaredNorm()); } /** \returns a row (or column) vector expression of the norm * of each column (or row) of the referenced expression. @@ -404,7 +485,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp * \sa DenseBase::mean() */ EIGEN_DEVICE_FUNC const MeanReturnType mean() const - { return MeanReturnType(_expression()); } + { return sum() / Scalar(Direction==Vertical?m_matrix.rows():m_matrix.cols()); } /** \returns a row (or column) vector expression representing * whether \b all coefficients of each respective column (or row) are \c true. @@ -500,7 +581,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) //eigen_assert((m_matrix.isNull()) == (other.isNull())); FIXME - return const_cast<ExpressionType&>(m_matrix = extendedTo(other.derived())); + return m_matrix = extendedTo(other.derived()); } /** Adds the vector \a other to each subvector of \c *this */ @@ -510,7 +591,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) - return const_cast<ExpressionType&>(m_matrix += extendedTo(other.derived())); + return m_matrix += extendedTo(other.derived()); } /** Substracts the vector \a other to each subvector of \c *this */ @@ -520,7 +601,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) - return const_cast<ExpressionType&>(m_matrix -= extendedTo(other.derived())); + return m_matrix -= extendedTo(other.derived()); } /** Multiples each subvector of \c *this by the vector \a other */ @@ -532,7 +613,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType) EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) m_matrix *= extendedTo(other.derived()); - return const_cast<ExpressionType&>(m_matrix); + return m_matrix; } /** Divides each subvector of \c *this by the vector \a other */ @@ -544,7 +625,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType) EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived) m_matrix /= extendedTo(other.derived()); - return const_cast<ExpressionType&>(m_matrix); + return m_matrix; } /** Returns the expression of the sum of the vector \a other to each subvector of \c *this */ @@ -609,7 +690,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp EIGEN_DEVICE_FUNC CwiseBinaryOp<internal::scalar_quotient_op<Scalar>, const ExpressionTypeNestedCleaned, - const typename OppositeExtendedType<typename ReturnType<internal::member_norm,RealScalar>::Type>::Type> + const typename OppositeExtendedType<NormReturnType>::Type> normalized() const { return m_matrix.cwiseQuotient(extendedToOpposite(this->norm())); } @@ -658,7 +739,15 @@ template<typename ExpressionType, int Direction> class VectorwiseOp EIGEN_DEVICE_FUNC const HNormalizedReturnType hnormalized() const; +# ifdef EIGEN_VECTORWISEOP_PLUGIN +# include EIGEN_VECTORWISEOP_PLUGIN +# endif + protected: + Index redux_length() const + { + return Direction==Vertical ? m_matrix.rows() : m_matrix.cols(); + } ExpressionTypeNested m_matrix; }; @@ -670,7 +759,7 @@ template<typename ExpressionType, int Direction> class VectorwiseOp * \sa rowwise(), class VectorwiseOp, \ref TutorialReductionsVisitorsBroadcasting */ template<typename Derived> -inline typename DenseBase<Derived>::ColwiseReturnType +EIGEN_DEVICE_FUNC inline typename DenseBase<Derived>::ColwiseReturnType DenseBase<Derived>::colwise() { return ColwiseReturnType(derived()); @@ -684,7 +773,7 @@ DenseBase<Derived>::colwise() * \sa colwise(), class VectorwiseOp, \ref TutorialReductionsVisitorsBroadcasting */ template<typename Derived> -inline typename DenseBase<Derived>::RowwiseReturnType +EIGEN_DEVICE_FUNC inline typename DenseBase<Derived>::RowwiseReturnType DenseBase<Derived>::rowwise() { return RowwiseReturnType(derived()); |