Re: [eigen] Matrix2i mean

[Petr Kubánek <pkubanek@xxxxxxxxx>]

Patch attached. Hopefully clarifies the documentation. I will see if I
will find time to work on better solution.

Petr

On Tue, 2019-11-05 at 17:11 +0100, Christoph Hertzberg wrote:
> Adding an additional template method would not be a problem (minimal 
> example: https://godbolt.org/z/ZB9WDI), though I don't see any real 
> advantage vs writing
> 
>      m.cast<int>().sum();
> 
> Sure, `m.sum<int>();` would be shorter, but also more confusing.
> 
> Your "another template parameter for the matrices" formulation
> sounded 
> to me like you wanted to add another template parameter to 
> `Eigen::Matrix` (which is not an option).
> 
> Any patches clarifying the documentation are welcome!
> 
> Christoph
> 
> 
> On 05/11/2019 16.44, Petr Kubánek wrote:
> > Hi,
> > 
> > I know what's the problem. I am just looking either to document it
> > or
> > for a solution.
> > 
> > Why would:
> > 
> > template <typename dt> dt sum()
> > 
> > break API/ABI? You will be able to call either sum() or
> > sum<int32_t>(),
> > shouldn't you? I will try that and submit a patch.
> > 
> > Having sum<dt>() working, one can hopefully create mean<dt_sum,
> > dt_mean>(), so one can code:
> > 
> > mean<int32_t,double>()
> > 
> > Petr
> > 
> > On Tue, 2019-11-05 at 16:36 +0100, Christoph Hertzberg wrote:
> > > On 05/11/2019 16.11, Peter wrote:
> > > > [...]
> > > > 
> > > > actually, this would also be interesting for the scalar
> > > > products
> > > > in
> > > > general, namely a different
> > > > type for accumulating the sums within a scalar product, e.g. as
> > > > yet
> > > > another template parameter for the matrices.
> > > 
> > > Adding another template parameter to basic types is not an
> > > option.
> > > This
> > > would break ABI and API compatibility (even if the parameter has
> > > a
> > > default).
> > > 
> > > You could create your own custom type `my_int16` for which
> > > `my_int16
> > > *
> > > my_int16` results in a `my_int32` (this needs to be told to
> > > Eigen,
> > > similar to how real*complex products are handled).
> > > 
> > > > > [...]
> > > > 
> > > > I think it's more subtle than that.
> > > > Even
> > > > 
> > > >     int16_t  Two =  2;
> > > >     int16_t  Max =  INT16_MAX;
> > > >     int16_t mean = ( Max/2 + Max + Two + Two ) / int16_t(4);
> > > > 
> > > > doesn't produce an overflow.
> > > 
> > > Yes, because `Max/2` gets implicitly converted to `int`.
> > > Actually,
> > > even
> > > adding two `int16_t` get implicitly converted to `int` (search
> > > for
> > > integer promotion rules -- I don't know them entirely either).
> > > And
> > > dividing an `int` by an `int16_t` results in an `int` which is
> > > only
> > > afterwards converted to `int16_t`.
> > > 
> > > In contrast, Eigen::DenseBase::sum() does something (more or
> > > less)
> > > equivalent to:
> > > 
> > >       T sum = T(0);
> > >       for(Index i=0; i< size(); ++i)
> > >           sum+=coeff(i);
> > >       return sum;
> > > 
> > > i.e., after each addition the result gets reduced to the scalar
> > > type
> > > of
> > > the matrix, thus it will immediately overflow.
> > > 
> > > And mean() essentially just takes the return value of `sum()` and
> > > divides it by the size.
> > > 
> > > 
> > > Christoph
> > > 
> > > 
> > 
> > 
> > 
> 
> 

# HG changeset patch
# User Petr Kubánek <pkubanek@xxxxxxxxx>
# Date 1573065305 25200
#      Wed Nov 06 11:35:05 2019 -0700
# Node ID 42bca2d654395d5356bc053f03f5e7e558273cf7
# Parent  afc120bc03bdc4265858d6f86218eb1fed51b1b9
improves reductions documentation

diff --git a/Eigen/src/Core/Redux.h b/Eigen/src/Core/Redux.h
--- a/Eigen/src/Core/Redux.h
+++ b/Eigen/src/Core/Redux.h
@@ -436,65 +436,69 @@ DenseBase<Derived>::minCoeff() const
 template<typename Derived>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
 DenseBase<Derived>::maxCoeff() const
 {
   return derived().redux(Eigen::internal::scalar_max_op<Scalar,Scalar>());
 }
 
 /** \returns the sum of all coefficients of \c *this
+  * \warning internal accumulator type and returned type matches base type. Overflow might occur for e.g. int values. Use \link MatrixBase::cast() .cast() \endlink to cast your variable to proper type.
   *
   * If \c *this is empty, then the value 0 is returned.
   *
   * \sa trace(), prod(), mean()
   */
 template<typename Derived>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
 DenseBase<Derived>::sum() const
 {
   if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0))
     return Scalar(0);
   return derived().redux(Eigen::internal::scalar_sum_op<Scalar,Scalar>());
 }
 
 /** \returns the mean of all coefficients of *this
-*
-* \sa trace(), prod(), sum()
-*/
+  * \warning internal accumulator type and returned type matches base type. Overflow might occur for e.g. int values. Use \link MatrixBase::cast() .cast() \endlink to cast your variable to proper type.
+  *
+  * \sa trace(), prod(), sum()
+  */
 template<typename Derived>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
 DenseBase<Derived>::mean() const
 {
 #ifdef __INTEL_COMPILER
   #pragma warning push
   #pragma warning ( disable : 2259 )
 #endif
   return Scalar(derived().redux(Eigen::internal::scalar_sum_op<Scalar,Scalar>())) / Scalar(this->size());
 #ifdef __INTEL_COMPILER
   #pragma warning pop
 #endif
 }
 
 /** \returns the product of all coefficients of *this
+  * \warning internal accumulator type and returned type matches base type. Overflow might occur for e.g. int values. Use \link MatrixBase::cast() .cast() \endlink to cast your variable to proper type.
   *
   * Example: \include MatrixBase_prod.cpp
   * Output: \verbinclude MatrixBase_prod.out
   *
   * \sa sum(), mean(), trace()
   */
 template<typename Derived>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
 DenseBase<Derived>::prod() const
 {
   if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0))
     return Scalar(1);
   return derived().redux(Eigen::internal::scalar_product_op<Scalar>());
 }
 
 /** \returns the trace of \c *this, i.e. the sum of the coefficients on the main diagonal.
+  * \warning internal accumulator type and returned type matches base type. Overflow might occur for e.g. int values. Use \link MatrixBase::cast() .cast() \endlink to cast your variable to proper type.
   *
   * \c *this can be any matrix, not necessarily square.
   *
   * \sa diagonal(), sum()
   */
 template<typename Derived>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
 MatrixBase<Derived>::trace() const
diff --git a/doc/TutorialReductionsVisitorsBroadcasting.dox b/doc/TutorialReductionsVisitorsBroadcasting.dox
--- a/doc/TutorialReductionsVisitorsBroadcasting.dox
+++ b/doc/TutorialReductionsVisitorsBroadcasting.dox
@@ -16,16 +16,18 @@ returning the sum of all the coefficient
 <tr><th>Example:</th><th>Output:</th></tr>
 <tr><td>
 \include tut_arithmetic_redux_basic.cpp
 </td>
 <td>
 \verbinclude tut_arithmetic_redux_basic.out
 </td></tr></table>
 
+Please note \link DenseBase::sum() .sum() \endlink calculates sum using same storage type as the matrix or array.
+
 The \em trace of a matrix, as returned by the function \c trace(), is the sum of the diagonal coefficients and can equivalently be computed <tt>a.diagonal().sum()</tt>.
 
 
 \subsection TutorialReductionsVisitorsBroadcastingReductionsNorm Norm computations
 
 The (Euclidean a.k.a. \f$\ell^2\f$) squared norm of a vector can be obtained \link MatrixBase::squaredNorm() squaredNorm() \endlink. It is equal to the dot product of the vector by itself, and equivalently to the sum of squared absolute values of its coefficients.
 
 Eigen also provides the \link MatrixBase::norm() norm() \endlink method, which returns the square root of \link MatrixBase::squaredNorm() squaredNorm() \endlink.