Re: problem with sum (...)

[David Bateman]

John W. Eaton wrote:
> On  5-Nov-2007, David Bateman wrote:
> 
> | In any case consider the attached patch that adds this
> | functionality to "sum"..
> 
> I applied the patch for OPERATORS/op-bm-bm.cc but not the other part
> because it doesn't seem to implement the correctly functionality.
> If I'm understanding correctly, your patch computes the sum using
> doubles internally, then casts the result back to the original type.
> That's not quite the same as doing the sum in the native type, which
> is what the Matlab docs say they are doing.  For example, with the
> saturation rules for integer arithmetic, Matlab computes
> 
>   sum (int8 ([127, 10, -20]), 'native')
>   == (127 + 10) + -20
>   == 127 + -20
>   == 107

Not that I don't believe you, but I checked and with MatlabR2007a this
returns 107 as you stated.. Damn.

> 
> but with your patch, we get 117 because the saturation limit is only
> applied at the end instead of after each operation.

Ok, then what about the attached..

> 
> Note that with 'native' it is critical to sum the elements in a
> specific order.  Ugh.

The attached patch also adds a test for this

> 
> | Perhaps the other data reduction operators
> | should also include this functionality?
> 
> Hmm.  Maybe only if it is needed for compatibility.

Ok, then lets not add it to the other functions.

D.

> 
> jwe
> 


*** ./liboctave/intNDArray.cc.orig4     2007-11-06 22:22:29.686756124 +0100
--- ./liboctave/intNDArray.cc   2007-11-06 22:37:53.161623350 +0100
***************
*** 222,227 ****
--- 222,234 ----
    return is;
  }
  
+ template <class T>
+ intNDArray<T>
+ intNDArray<T>::sum (int dim) const
+ {
+   MX_ND_REDUCTION (retval(result_idx) += intNDArray<T>::elem (iter_idx), 0, 
intNDArray<T>);
+ }
+ 
  /*
  ;;; Local Variables: ***
  ;;; mode: C++ ***
*** ./liboctave/intNDArray.h.orig4      2007-11-06 22:22:43.340044496 +0100
--- ./liboctave/intNDArray.h    2007-11-06 22:21:46.070029488 +0100
***************
*** 73,78 ****
--- 73,80 ----
    boolNDArray all (int dim = -1) const;
    boolNDArray any (int dim = -1) const;
  
+   intNDArray sum (int dim) const;
+ 
    intNDArray squeeze (void) const
      { return intNDArray<T> (MArrayN<T>::squeeze ()); }
  
*** ./liboctave/boolNDArray.h.orig4     2007-11-06 22:25:56.903955680 +0100
--- ./liboctave/boolNDArray.h   2007-11-06 22:24:40.850919669 +0100
***************
*** 66,71 ****
--- 66,73 ----
    boolNDArray all (int dim = -1) const;
    boolNDArray any (int dim = -1) const;
  
+   boolNDArray sum (int dim = -1) const;
+ 
    boolNDArray concat (const boolNDArray& rb, const Array<octave_idx_type>& 
ra_idx);
  
    boolNDArray& insert (const boolNDArray& a, octave_idx_type r, 
octave_idx_type c);
*** ./liboctave/boolNDArray.cc.orig4    2007-11-06 22:26:06.176472384 +0100
--- ./liboctave/boolNDArray.cc  2007-11-06 22:34:42.813544553 +0100
***************
*** 57,62 ****
--- 57,68 ----
    MX_ND_ANY_ALL_REDUCTION (MX_ND_ANY_EVAL (MX_ND_ANY_EXPR), false);
  }
  
+ boolNDArray 
+ boolNDArray::sum (int dim) const
+ {
+   MX_ND_REDUCTION (retval(result_idx) |= elem (iter_idx), true, boolNDArray);
+ }
+ 
  boolNDArray
  boolNDArray::concat (const boolNDArray& rb, const Array<octave_idx_type>& 
ra_idx)
  {
*** ./src/data.cc.orig4 2007-11-05 14:51:04.238672315 +0100
--- ./src/data.cc       2007-11-06 23:25:09.056812611 +0100
***************
*** 370,375 ****
--- 370,481 ----
    return retval;
  }
  
+ #define NATIVE_REDUCTION_1(FCN, TYPE, DIM) \
+   (arg.is_ ## TYPE ## _type ()) \
+     { \
+       TYPE ## NDArray tmp = arg. TYPE ##_array_value (); \
+       \
+       if (! error_state) \
+         retval = tmp.FCN (DIM); \
+     }
+ 
+ #define NATIVE_REDUCTION(FCN) \
+  \
+   octave_value retval; \
+  \
+   int nargin = args.length (); \
+  \
+   bool isnative = false; \
+   \
+   if (nargin > 1 && args(nargin - 1).is_string ()) \
+     { \
+       std::string str = args(nargin - 1).string_value (); \
+       \
+       if (! error_state) \
+       { \
+         if (str == "native") \
+           isnative = true; \
+         else if (str != "double") /* Ignore double as no single type */ \
+           error ("sum: unrecognized string argument"); \
+           nargin --; \
+       } \
+     } \
+   \
+   if (nargin == 1 || nargin == 2) \
+     { \
+       octave_value arg = args(0); \
+  \
+       int dim = (nargin == 1 ? -1 : args(1).int_value (true) - 1); \
+  \
+       if (! error_state) \
+       { \
+         if (dim >= -1) \
+           { \
+               if (isnative) \
+                 { \
+                   if NATIVE_REDUCTION_1 (FCN, uint8, dim) \
+                   else if NATIVE_REDUCTION_1 (FCN, uint16, dim) \
+                   else if NATIVE_REDUCTION_1 (FCN, uint32, dim) \
+                   else if NATIVE_REDUCTION_1 (FCN, uint64, dim) \
+                   else if NATIVE_REDUCTION_1 (FCN, int8, dim) \
+                   else if NATIVE_REDUCTION_1 (FCN, int16, dim) \
+                   else if NATIVE_REDUCTION_1 (FCN, int32, dim) \
+                   else if NATIVE_REDUCTION_1 (FCN, int64, dim) \
+                   else if NATIVE_REDUCTION_1 (FCN, bool, dim) \
+                   else if (arg.is_char_matrix ()) \
+                     { \
+                        error (#FCN, ": invalid char type"); \
+                     } \
+                 else if (arg.is_complex_type ()) \
+                   { \
+                     ComplexNDArray tmp = arg.complex_array_value (); \
+                       \
+                     if (! error_state) \
+                       retval = tmp.FCN (dim); \
+                   } \
+                 else if (arg.is_real_type ()) \
+                   { \
+                     NDArray tmp = arg.array_value (); \
+                       \
+                     if (! error_state) \
+                       retval = tmp.FCN (dim); \
+                   } \
+                   else \
+                   { \
+                     gripe_wrong_type_arg (#FCN, arg); \
+                     return retval; \
+                   } \
+                 } \
+             else if (arg.is_real_type ()) \
+               { \
+                 NDArray tmp = arg.array_value (); \
+                   \
+                 if (! error_state) \
+                   retval = tmp.FCN (dim); \
+               } \
+             else if (arg.is_complex_type ()) \
+               { \
+                 ComplexNDArray tmp = arg.complex_array_value (); \
+                   \
+                 if (! error_state) \
+                   retval = tmp.FCN (dim); \
+               } \
+             else \
+               { \
+                 gripe_wrong_type_arg (#FCN, arg); \
+                 return retval; \
+               } \
+           } \
+         else \
+           error (#FCN ": invalid dimension argument = %d", dim + 1); \
+       } \
+       \
+     } \
+   else \
+     print_usage (); \
+  \
+   return retval
+ 
  #define DATA_REDUCTION(FCN) \
   \
    octave_value retval; \
***************
*** 1213,1228 ****
  DEFUN (sum, args, ,
    "-*- texinfo -*-\n\
  @deftypefn {Built-in Function} {} sum (@var{x}, @var{dim})\n\
  Sum of elements along dimension @var{dim}.  If @var{dim} is\n\
  omitted, it defaults to 1 (column-wise sum).\n\
  \n\
  As a special case, if @var{x} is a vector and @var{dim} is omitted,\n\
  return the sum of the elements.\n\
  @end deftypefn")
  {
!   DATA_REDUCTION (sum);
  }
  
  DEFUN (sumsq, args, ,
    "-*- texinfo -*-\n\
  @deftypefn {Built-in Function} {} sumsq (@var{x}, @var{dim})\n\
--- 1319,1355 ----
  DEFUN (sum, args, ,
    "-*- texinfo -*-\n\
  @deftypefn {Built-in Function} {} sum (@var{x}, @var{dim})\n\
+ @deftypefnx {Built-in Function} {} sum (@dots{}, 'native')\n\
  Sum of elements along dimension @var{dim}.  If @var{dim} is\n\
  omitted, it defaults to 1 (column-wise sum).\n\
  \n\
  As a special case, if @var{x} is a vector and @var{dim} is omitted,\n\
  return the sum of the elements.\n\
+ \n\
+ If the optional argument 'native' is given, then the sum is performed\n\
+ in the same type as the original argument, rather than in the default\n\
+ double type. For example\n\
+ \n\
+ @example\n\
+ sum ([true, true])\n\
+   @result{} 2\n\
+ sum ([true, true], 'native')\n\
+   @result{} true\n\
+ @end example\n\
  @end deftypefn")
  {
!   NATIVE_REDUCTION (sum);
  }
  
+ /*
+ 
+ %!assert (sum([true,true]), 2)
+ %!assert (sum([true,true],'native'), true)
+ %!assert (sum(int8([127,10,-20])), 117);
+ %!assert (sum(int8([127,10,-20]),'native'), int8(107));
+ 
+ */
+ 
  DEFUN (sumsq, args, ,
    "-*- texinfo -*-\n\
  @deftypefn {Built-in Function} {} sumsq (@var{x}, @var{dim})\n\
***************
*** 2613,2619 ****
  %!assert(norm(x,Inf), 7);
  %!assert(norm(x,-Inf), 1);
  %!assert(norm(x,"inf"), 7);
! %!assert(norm(x,"fro"), 10);
  %!assert(norm(x), 10);
  %!assert(norm([1e200, 1]), 1e200);
  %!assert(norm([3+4i, 3-4i, sqrt(31)]), 9, -4*eps);
--- 2740,2746 ----
  %!assert(norm(x,Inf), 7);
  %!assert(norm(x,-Inf), 1);
  %!assert(norm(x,"inf"), 7);
! %!assert(norm(x,"fro"), 10, -eps);
  %!assert(norm(x), 10);
  %!assert(norm([1e200, 1]), 1e200);
  %!assert(norm([3+4i, 3-4i, sqrt(31)]), 9, -4*eps);
***************
*** 2623,2628 ****
--- 2750,2761 ----
  %!assert(norm(m,2), 34, -eps);
  %!assert(norm(m,Inf), 34);
  %!assert(norm(m,"inf"), 34);
+ %!shared m2, flo, fhi
+ %! m2 = [1,2;3,4];
+ %! flo = 1e-300;
+ %! fhi = 1e+300;
+ %!assert (norm(flo*m2,"fro"), sqrt(30)*flo, -eps)
+ %!assert (norm(fhi*m2,"fro"), sqrt(30)*fhi, -eps)
  */
  
  // Compute various norms of the vector X.
***************
*** 2684,2715 ****
          if ((x_arg.rows () == 1 || x_arg.columns () == 1)
              && ! (x_arg.is_sparse_type () || x_arg.is_integer_type ()))
            {
!             double p_val;
  
!             octave_value p_arg;
! 
!             if (nargin == 1)
!               p_arg = 2;
!             else
!               p_arg = args(1);
! 
!             if (p_arg.is_string ())
                {
!                 std::string p = args(1).string_value ();
  
!                 if (p == "inf")
!                   p_val = octave_Inf;
!                 else if (p == "fro")
!                   p_val = -1;
                  else
!                   error ("norm: unrecognized norm `%s'", p.c_str ());
!               }
!             else
!               {
!                 p_val = p_arg.double_value ();
  
!                 if (error_state)
!                   error ("norm: unrecognized norm value");
                }
  
              if (! error_state)
--- 2817,2846 ----
          if ((x_arg.rows () == 1 || x_arg.columns () == 1)
              && ! (x_arg.is_sparse_type () || x_arg.is_integer_type ()))
            {
!             double p_val = 2;
  
!             if (nargin == 2)
                {
!                 octave_value p_arg = args(1);
! 
!                 if (p_arg.is_string ())
!                   {
!                     std::string p = args(1).string_value ();
  
!                     if (p == "inf")
!                       p_val = octave_Inf;
!                     else if (p == "fro")
!                       p_val = -1;
!                     else
!                       error ("norm: unrecognized norm `%s'", p.c_str ());
!                   }
                  else
!                   {
!                     p_val = p_arg.double_value ();
  
!                     if (error_state)
!                       error ("norm: unrecognized norm value");
!                   }
                }
  
              if (! error_state)

2007-11-06  David Bateman  <address@hidden>

        * data.cc (DATA_REDUCTION): Handle the 'native' and 'double'
        arguments of the Fsum function.