[Numpy-discussion] A New Coercion Model

[Friedrich Romstedt]

Hi,

there recently were some problems occuring when coercing
numpy.ndarrays with something else.  I would like to try to review the
examples I've seen so far (it are two examples in only ~ one week),
and also would try to summarise the approach taken up to now as far as
I understand.  Furthermore, I would like to start discussion about a
more robust coercion model.

I'm not shure whether I'm in the correct position for starting such a
deep thread, but I simply hope you appreciate and excuse me for trying
to do so.

EXAMPLES

1. numpy.poly1d.

The problem arises when coercing a numpy.poly1d instance as second
operand with a numpy.ndarray (e.g., a scalar array) or a numpy scalar
(like a numpy.int32 instance).  It seems that numpy tries to create an
array out of the poly1d instance, maybe via asarray(), and this
results indeed in an ndarray.

But it isn't clear to me, why.  poly1d supports __len__() and
__getitem__(), but one can easily check that a class supporting this
is _not_ converted into an ndarray by numpy.asarray().  Also int32(2)
* X([1,2]) where X is the resp. class does yield TypeError:
unsupported operand type(s) for *: 'int' and 'instance'.  Same applies
for a class supporting __iter__().  Also same for a class supporting
__iter__, __len__, __getitem__ simultaneously.

instance(numpy.poly1d([1]), numpy.ndarray) returns False.

2. upy.undarray.

As I reported myself, I encountered the same problem with my own code
with also non-numpy.ndarray objects. The class supports also __len__()
and __getitem__(), but I found that it is treated as a scalar.
(Nevertheless, changing this behaviour wouldn't solve the problem.)

numpy.asarray(upy.undarray(...)) returns a scalar array, opposed to
the numpy.asarray(poly1d(...)) case, where it returns a corresponding
ndarray.

It seems to do not matter, that I never tried with numpy.int32
instances on the left hand side of, e.g., multiplication.

3. numpy.polynomial.Polynomial.

I have no idea how this was solved for Polynomial.  Even
numpy.int32(2).__mul__(p) with p being a Polynomial instance works
fine.  Maybe it's a static route?  Polynomial([1,
2]).__array_priority__ does not exist.

CURRENT APPROACH

So, this seems easy: There is __array_priority__, there may be static
routes hardcoded, and there is numpy.set_numeric_ops().  Are there
more ways implemented to treat the problem I'm not aware of?

DISCUSSION

__array_priority__ induces a linear order on all classes in the
namespace.  I think, I guess, this is not appropriate.  At least my
intuition tells me that this easily breaks down.

Consider A and B having precedence over ndarray, for supporting an
expression with A() as right operand.  Now, C is introduced, which
wants to have precedence over A, but B leaving more precedent than the
new C.  Now everything depends on the definition of A and B.  The aim
is unreachable if the __array_priority__s of A and B are by
coincidence not compatible with this new elements of the relation.
(It's maybe an a bit silly example, but I have no better at hand.)

In fact, no one knows the precedence of something over anything else
as long as there is no definition for that.  There may occur even
rings in the relation.  Calling the relation >, although it's no
longer assumed to be linear, there may hold A > B > C > A.

What about implementing a simple class "Relation" in numpy.  Users may
register relations in the numpy.relations instance they are shure to
be existent.  For instance, I would say numpy.relations.set(lower =
numpy.ndarray, higher = upy.undarray).  Same for the numpy.poly1d
thing, and it could also be used for numpy.polynomial.Polynomial.
Though there are subtleties with inheritance.  E.g. when class
X(numpy.ndarray): [...], it should be treated as numpy.ndarray as long
as there is no definition involving X directly.  One could simply say,
the last defined relation rules first, "LIFO".

I could code a Python module for this, but it would slow numpy down.
Maybe a C implementation would be helpful.  I have much experience
with C++ and Python both, but no experience with building numpy ...

But now, I would appreciate any response in discussion.

Friedrich

P.S. To me it's clear that this would apply to 3.0 (or whatever) ...