element-wise exponent operator

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

element-wise exponent operator

Rik-4
jwe, maintainers,

I started trying to understand bug #57671 (Power operator operating on
complex numbers results in NaN when broadcasting) a week ago.  The problem
behavior is

([1;0]*j) .^ [1, 0]
ans =

   0.0000 + 1.0000i   1.0000 +      0i
        0 +      0i      NaN -    NaNi

But, if the matrices are full so that there is no broadcasting then Octave
works correctly

([1,1; 0,0]*j) .^ [1, 0; 1, 0]
ans =

   0 + 1i   1 + 0i
   0 + 0i   1 + 0i

After obtuse debugging sessions with gdb, I found that it is not
broadcasting per se which is the problem, but how Octave is calling
std::pow from the C++ library.  The bug manifests in many different ways. 
For example, besides broadcasting, real scalars behave differently than
complex scalars.

0.^0
ans = 1

But, for complex scalars,

complex (0,0) .^ complex (0,0)
ans =  NaN + NaNi

Bugs #57908, #57671, #57476, #54302, #52706 are all manifestations of the
same issue.

My first request is for an explanation of the convoluted mapping strategy
that Octave uses for operators.  There are tons of C++ files in
libinterp/operators/ which seem to map all of the various combinations of a
scalar, vector, or matrix object with another scalar, vector, or matrix
object.  And this structure is also echoed in liboctave/operators.  And
then there are exceptions to this structure such as
libinterp/corefcn/xpow.cc which slots in to the architecture as it gets
called from libinterp/operators, but isn't stored there.  Why aren't we
just mapping the '+' operator to a BLAS call?  Or maybe we are and I don't
understand how we get there.

For bug #57671, the code which ends up getting called is in xpow.cc case #10.

// -*- 10 -*-
octave_value
elem_xpow (const ComplexNDArray& a, const NDArray& b)
{
  dim_vector a_dims = a.dims ();
  dim_vector b_dims = b.dims ();

  if (a_dims != b_dims)
    {
      if (! is_valid_bsxfun ("operator .^", a_dims, b_dims))
        octave::err_nonconformant ("operator .^", a_dims, b_dims);

      return bsxfun_pow (a, b);
    }

  ComplexNDArray result (a_dims);

  for (octave_idx_type i = 0; i < a.numel (); i++)
    {
      octave_quit ();
      double btmp = b(i);
      if (xisint (btmp))
        result(i) = std::pow (a(i), static_cast<int> (btmp));
      else
        result(i) = std::pow (a(i), btmp);
    }

  return result;
}

For the non-broadcasting case, which works, the key is to note the check
for an integer power (xisint (btmp)).  If this is found it calls std::pow
from the C++ library with a different function signature.  When
broadcasting is used it calls bsxfun_pow and this ends up calling code from
liboctave.  In particular, the individual operation is in
liboctave/operators/mx-inlines.cc

template <typename R, typename X, typename Y>
inline void
mx_inline_pow (size_t n, R *r, const X *x, const Y *y)
{
  using std::pow;

  for (size_t i = 0; i < n; i++)
    r[i] = pow (x[i], y[i]);
}

This calls std::pow with a function signature of two reals which results in
NaN + NaNi.

At a minimum, all the uses of std::pow in xpow.cc and mx-inlines.cc should
be checked to see if they need to install a test for an integer and then
perform the same cast as is done for case 10.

Handling Complex values will be a little more difficult.  The function
signatures available for pow are

template<class T> complex<T> pow (const complex<T>& x, int y);
template<class T> complex<T> pow (const complex<T>& x, const complex<T>& y);
template<class T> complex<T> pow (const complex<T>& x, const T& y);
template<class T> complex<T> pow (const T& x, const complex<T>& y);

There is no signature for pow (complex<real>, complex<int>) which means we
likely need to check first for whether the exponent can be narrowed to
real, and then checked subsequently for whether it is an integer.  As an
example,

[0, j] .^ complex (0,0)
ans =

   NaN + NaNi     1 +   0i

But, if the code in xpow.cc is modified as shown below

// -*- 11 -*-
octave_value
elem_xpow (const ComplexNDArray& a, const Complex& b)
{
  ComplexNDArray result (a.dims ());

  bool isint = (b.imag() == 0) && xisint (b.real ());
  int bint = static_cast<int> (b.real ());

  for (octave_idx_type i = 0; i < a.numel (); i++)
    {
      octave_quit ();
      if (isint)
        result(i) = std::pow (a(i), bint);
      else
        result(i) = std::pow (a(i), b);
    }

  return result;
}


[0, j] .^ complex (0,0)
ans =

   1   1

All of this will be tedious, because there are so many cases, which makes
one wonder whether the code in xpow.cc can't make better use of templates.

Regardless, once fixed we also need to add BIST tests for all the weird
corner cases.

--Rik


Reply | Threaded
Open this post in threaded view
|

Re: element-wise exponent operator

John W. Eaton
Administrator
On 2/28/20 1:09 PM, Rik wrote:

> I started trying to understand bug #57671 (Power operator operating on
> complex numbers results in NaN when broadcasting) a week ago.  The problem
> behavior is
>
> ([1;0]*j) .^ [1, 0]
> ans =
>
>     0.0000 + 1.0000i   1.0000 +      0i
>          0 +      0i      NaN -    NaNi
>
> [...]
 >
> My first request is for an explanation of the convoluted mapping strategy
> that Octave uses for operators.  There are tons of C++ files in
> libinterp/operators/ which seem to map all of the various combinations of a
> scalar, vector, or matrix object with another scalar, vector, or matrix
> object.

This is a very old part of Octave from before C++ had templates and
before Octave hadt classes (for dispatching based on type, which could
also be used for built-in types) or broadcasting.

Since

   scalar .* matrix

behaves differently from

   matrix .* matrix

I decided to have actual scalar types in Octave so that we could write
operators that did special things for mixed-type operations like the
above.  Then each operator function could be simpler (no special check
is needed for scalar-by-matrix operations if matrices are always
narrowed to scalars by the interpreter before dispatching to specific
operators.  But you have more operator functions, and when you add in
complex versions, a LOT more.  But I thought it was worth it to avoid
wasting memory and time expanding (possibly large) matrices from real to
complex if we had special operations for mixed types where possible and
only resort to type conversions when mixed type operations won't work.

To determine which function to call, we have lookup tables for operators
based on the type of both operands and we may also apply type conversion
operators.  I agree this is all a bit complicated.  Maybe we can
simplify some of it.

> And this structure is also echoed in liboctave/operators.

I've been thinking about that a little bit recently while looking at the
MEX interface.  We have things like Matrix, ComplexMatrix, etc. in
liboctave and a wrappers for them so they can be octave_value objects
with a common interface for the interpreter.  Matlab appears to just
have mxArray that contains various types.  We could do the same by
moving the array storage directly into an object in the octave_value
hierarchy, but I thought it would be more convenient for interfacing
with and writing numerical code if we had separate array types that were
independent of the Octave interpreter.

> And
> then there are exceptions to this structure such as
> libinterp/corefcn/xpow.cc which slots in to the architecture as it gets
> called from libinterp/operators, but isn't stored there.

The operators directory was originally intended to be a place to store
the wrapper functions that all had the same signatures (for the table of
operators in the interpreter) and xpow.cc was for the type-specific
operations with Matrix, ComplexMatrix, etc. arguments.  It's in
libinterp instead of liboctave because the return type is octave_value
because pow operations on two real values can sometimes produce complex.

> Why aren't we
> just mapping the '+' operator to a BLAS call?  Or maybe we are and I don't
> understand how we get there.

Are there BLAS functions for vector addition, subtraction,
multiplication, or division, or for mixed-type operations?

> For bug #57671, the code which ends up getting called is in xpow.cc case #10.
>
> // -*- 10 -*-
> octave_value
> elem_xpow (const ComplexNDArray& a, const NDArray& b)
> {
>    dim_vector a_dims = a.dims ();
>    dim_vector b_dims = b.dims ();
>
>    if (a_dims != b_dims)
>      {
>        if (! is_valid_bsxfun ("operator .^", a_dims, b_dims))
>          octave::err_nonconformant ("operator .^", a_dims, b_dims);
>
>        return bsxfun_pow (a, b);
>      }
>
>    ComplexNDArray result (a_dims);
>
>    for (octave_idx_type i = 0; i < a.numel (); i++)
>      {
>        octave_quit ();
>        double btmp = b(i);
>        if (xisint (btmp))
>          result(i) = std::pow (a(i), static_cast<int> (btmp));
>        else
>          result(i) = std::pow (a(i), btmp);
>      }
>
>    return result;
> }
>
> For the non-broadcasting case, which works, the key is to note the check
> for an integer power (xisint (btmp)).  If this is found it calls std::pow
> from the C++ library with a different function signature.

Should we be using a wrapper for pow and burying this logic inside a
template function?

> When
> broadcasting is used it calls bsxfun_pow and this ends up calling code from
> liboctave.  In particular, the individual operation is in
> liboctave/operators/mx-inlines.cc
>
> template <typename R, typename X, typename Y>
> inline void
> mx_inline_pow (size_t n, R *r, const X *x, const Y *y)
> {
>    using std::pow;
>
>    for (size_t i = 0; i < n; i++)
>      r[i] = pow (x[i], y[i]);
> }
>  > This calls std::pow with a function signature of two reals which
results in
> NaN + NaNi.
>
> At a minimum, all the uses of std::pow in xpow.cc and mx-inlines.cc should
> be checked to see if they need to install a test for an integer and then
> perform the same cast as is done for case 10.

It seems like we should be doing this with a wrapper function for pow to
get the behavior we want.

There may also be better ways to write the templates for brodcasting.

> Handling Complex values will be a little more difficult.  The function
> signatures available for pow are
>
> template<class T> complex<T> pow (const complex<T>& x, int y);
> template<class T> complex<T> pow (const complex<T>& x, const complex<T>& y);
> template<class T> complex<T> pow (const complex<T>& x, const T& y);
> template<class T> complex<T> pow (const T& x, const complex<T>& y);
>
> There is no signature for pow (complex<real>, complex<int>) which means we
> likely need to check first for whether the exponent can be narrowed to
> real, and then checked subsequently for whether it is an integer.

To do this correctly everywhere would probably be easiest if we provided
our own Complex class and used it consistently or at least had a wrapper
function for pow that we used instead of using the ones from
std::complex directly.

> All of this will be tedious, because there are so many cases, which makes
> one wonder whether the code in xpow.cc can't make better use of templates.

I'm sure it can.  The functions there were originally written long
before C++ had useful templates.

> Regardless, once fixed we also need to add BIST tests for all the weird
> corner cases.

Yes.

jwe



Reply | Threaded
Open this post in threaded view
|

Re: element-wise exponent operator

Rik-4
On 03/10/2020 10:10 AM, John W. Eaton wrote:
On 2/28/20 1:09 PM, Rik wrote:

I started trying to understand bug #57671 (Power operator operating on
complex numbers results in NaN when broadcasting) a week ago.  The problem
behavior is

([1;0]*j) .^ [1, 0]
ans =

    0.0000 + 1.0000i   1.0000 +      0i
         0 +      0i      NaN -    NaNi

[...]
>
My first request is for an explanation of the convoluted mapping strategy
that Octave uses for operators.  There are tons of C++ files in
libinterp/operators/ which seem to map all of the various combinations of a
scalar, vector, or matrix object with another scalar, vector, or matrix
object.

This is a very old part of Octave from before C++ had templates and before Octave hadt classes (for dispatching based on type, which could also be used for built-in types) or broadcasting.

Since

  scalar .* matrix

behaves differently from

  matrix .* matrix

I decided to have actual scalar types in Octave so that we could write operators that did special things for mixed-type operations like the above.  Then each operator function could be simpler (no special check is needed for scalar-by-matrix operations if matrices are always narrowed to scalars by the interpreter before dispatching to specific operators.  But you have more operator functions, and when you add in complex versions, a LOT more.  But I thought it was worth it to avoid wasting memory and time expanding (possibly large) matrices from real to complex if we had special operations for mixed types where possible and only resort to type conversions when mixed type operations won't work.

To determine which function to call, we have lookup tables for operators based on the type of both operands and we may also apply type conversion operators.  I agree this is all a bit complicated.  Maybe we can simplify some of it.

I added a project to the wiki about this:

Re-implement operators using templates and modern C++. Current system evolved before templates and makes extensive use of macros to define interactions between scalar<->scalar, scalar<->matrix, scalar<->float, etc., etc.
  • In liboctave, the directory to work on is liboctave/operators
  • In libinterp, the directory to work on is libinterp/operators
  • In libinterp, there is also xpow.cc, xdiv.cc in libinterp/corefcn



And this structure is also echoed in liboctave/operators.

I've been thinking about that a little bit recently while looking at the MEX interface.  We have things like Matrix, ComplexMatrix, etc. in liboctave and a wrappers for them so they can be octave_value objects with a common interface for the interpreter.  Matlab appears to just have mxArray that contains various types.  We could do the same by moving the array storage directly into an object in the octave_value hierarchy, but I thought it would be more convenient for interfacing with and writing numerical code if we had separate array types that were independent of the Octave interpreter.

And
then there are exceptions to this structure such as
libinterp/corefcn/xpow.cc which slots in to the architecture as it gets
called from libinterp/operators, but isn't stored there.

The operators directory was originally intended to be a place to store the wrapper functions that all had the same signatures (for the table of operators in the interpreter) and xpow.cc was for the type-specific operations with Matrix, ComplexMatrix, etc. arguments.  It's in libinterp instead of liboctave because the return type is octave_value because pow operations on two real values can sometimes produce complex.

Why aren't we
just mapping the '+' operator to a BLAS call?  Or maybe we are and I don't
understand how we get there.

Are there BLAS functions for vector addition, subtraction, multiplication, or division, or for mixed-type operations?

No, there aren't.  I just assumed there must be, but checking the BLAS library (http://www.netlib.org/blas/) shows nothing of the sort.


For bug #57671, the code which ends up getting called is in xpow.cc case #10.

// -*- 10 -*-
octave_value
elem_xpow (const ComplexNDArray& a, const NDArray& b)
{
   dim_vector a_dims = a.dims ();
   dim_vector b_dims = b.dims ();

   if (a_dims != b_dims)
     {
       if (! is_valid_bsxfun ("operator .^", a_dims, b_dims))
         octave::err_nonconformant ("operator .^", a_dims, b_dims);

       return bsxfun_pow (a, b);
     }

   ComplexNDArray result (a_dims);

   for (octave_idx_type i = 0; i < a.numel (); i++)
     {
       octave_quit ();
       double btmp = b(i);
       if (xisint (btmp))
         result(i) = std::pow (a(i), static_cast<int> (btmp));
       else
         result(i) = std::pow (a(i), btmp);
     }

   return result;
}

For the non-broadcasting case, which works, the key is to note the check
for an integer power (xisint (btmp)).  If this is found it calls std::pow
from the C++ library with a different function signature.

Should we be using a wrapper for pow and burying this logic inside a template function?

Probably, but I'm going to lump that change in with converting most of the system to use templates.  For the time being, I will patch the code we have in xpow.cc and mx-inlines.cc to handle these exceptional cases.


When
broadcasting is used it calls bsxfun_pow and this ends up calling code from
liboctave.  In particular, the individual operation is in
liboctave/operators/mx-inlines.cc

template <typename R, typename X, typename Y>
inline void
mx_inline_pow (size_t n, R *r, const X *x, const Y *y)
{
   using std::pow;

   for (size_t i = 0; i < n; i++)
     r[i] = pow (x[i], y[i]);
}
 > This calls std::pow with a function signature of two reals which
results in
NaN + NaNi.

At a minimum, all the uses of std::pow in xpow.cc and mx-inlines.cc should
be checked to see if they need to install a test for an integer and then
perform the same cast as is done for case 10.

It seems like we should be doing this with a wrapper function for pow to get the behavior we want.

There may also be better ways to write the templates for brodcasting.

Handling Complex values will be a little more difficult.  The function
signatures available for pow are

template<class T> complex<T> pow (const complex<T>& x, int y);
template<class T> complex<T> pow (const complex<T>& x, const complex<T>& y);
template<class T> complex<T> pow (const complex<T>& x, const T& y);
template<class T> complex<T> pow (const T& x, const complex<T>& y);

There is no signature for pow (complex<real>, complex<int>) which means we
likely need to check first for whether the exponent can be narrowed to
real, and then checked subsequently for whether it is an integer.

To do this correctly everywhere would probably be easiest if we provided our own Complex class and used it consistently or at least had a wrapper function for pow that we used instead of using the ones from std::complex directly.

All of this will be tedious, because there are so many cases, which makes
one wonder whether the code in xpow.cc can't make better use of templates.

I'm sure it can.  The functions there were originally written long before C++ had useful templates.

Regardless, once fixed we also need to add BIST tests for all the weird
corner cases.


I added this to the wiki Projects page under the Testing category.

*** thorough tests for power operator including corner cases and strange combinations such as complex .^ range.

--Rik