[Numpy-discussion] Rank-0 arrays - reprise

[Nathaniel Smith]

On 5 Jan 2013 12:16, "Matthew Brett" <matthew.brett at gmail.com> wrote:
>
> Hi,
>
> Following on from Nathaniel's explorations of the scalar - array
> casting rules, some resources on rank-0 arrays.
>
> The discussion that Nathaniel tracked down on "rank-0 arrays"; it also
> makes reference to casting.  The rank-0 arrays seem to have been one
> way of solving the problem of maintaining array dtypes other than bool
> / float / int:
>
> http://mail.scipy.org/pipermail/numpy-discussion/2002-September/001612.html
>
> Quoting from an email from Travis in that thread, replying to an email
> from Tim Hochberg:
>
> http://mail.scipy.org/pipermail/numpy-discussion/2002-September/001647.html
>
> <quote>
> > Frankly, I have no idea what the implimentation details would be, but
> > could we get rid of rank-0 arrays altogether? I have always simply found
> > them strange and confusing... What are they really neccesary for
> > (besides holding scalar values of different precision that standard
> > Pyton scalars)?
>
> With new coercion rules this becomes a possibility.  Arguments against it
> are that  special rank-0 arrays behave as more consistent numbers with the
> rest of Numeric than Python scalars.  In other words they have a length
> and a shape and one can right N-dimensional code that works the same even
> when the result is a scalar.
>
> Another advantage of having a Numeric scalar is that we can control the
> behavior of floating point operations better.
>
> e.g.
>
> if only Python scalars were available and sum(a) returned 0, then
>
>  1 / sum(a)  would behave as Python behaves (always raises error).
>
> while with our own scalars
>
> 1 / sum(a)   could potentially behave however the user wanted.
> </quote>
>
> There seemed then to be some impetus to remove rank-0 arrays and
> replace them with Python scalar types with the various numpy
> precisions :
>
> http://mail.scipy.org/pipermail/numpy-discussion/2002-September/013983.html
>
> Travis' recent email hints at something that seems similar, but I
> don't understand what he means:
>
> http://mail.scipy.org/pipermail/numpy-discussion/2012-December/064795.html
>
> <quote>
> Don't create array-scalars.  Instead, make the data-type object a
> meta-type object whose instances are the items returned from NumPy
> arrays.   There is no need for a separate array-scalar object and in
> fact it's confusing to the type-system.    I understand that now.  I
> did not understand that 5 years ago.
> </quote>
>
> Travis - can you expand?

Numpy has 3 partially overlapping concepts:

A) scalars (what Travis calls "array scalars"): Things like "float64",
"int32". These are ordinary Python classes; usually when you subscript
an array, what you get back is an instance of one of these classes:

In [1]: a = np.array([1, 2, 3])

In [2]: a[0]
Out[2]: 1

In [3]: type(a[0])
Out[3]: numpy.int64

Note that even though they are called "array scalars", they have
nothing to do with the actual ndarray type -- they are totally
separate objects.

B) dtypes: These are instances of class np.dtype. For every scalar
type, there is a corresponding dtype object; plus you can create new
dtype objects for things like record arrays (which correspond to
scalars of type "np.void"; I don't really understand how void scalars
work in detail):

In [8]: int64_dtype = np.dtype(np.int64)

In [9]: int64_dtype
Out[9]: dtype('int64')

In [10]: type(int64_dtype)
Out[10]: numpy.dtype

In [11]: int64_dtype.type
Out[11]: numpy.int64

C) rank-0 arrays: Plain old ndarray objects that happen to have ndim
== 0, shape == (). These are arrays which are scalars, but they are
not array scalars. Arrays HAVE-A dtype.

In [15]: int64_arr = np.array(1)

In [16]: int64_arr
Out[16]: array(1)

In [17]: int64_arr.dtype
Out[17]: dtype('int64')

------------

Okay given that background:

What Travis was saying in that email was that he thought (A) and (B)
should be combined. Instead of having np.float64-the-class and
dtype(np.float64)-the-dtype-object, we should make dtype objects
actually *be* the scalar classes. (They would still be dtype objects,
which means they would be "metaclasses", which is just a fancy way to
say, dtype would be a subclass of the Python class "type", and dtype
objects would be class objects that had extra functionality.)

Those old mailing list threads are debating about (A) versus (C). What
we ended up with is what I described above -- we have "rank-0"
(0-dimensional) arrays, and we have array scalar objects that are a
different set of python types and objects entirely. The actual
implementation is totally different -- to the point that we a 35,000
line auto-generated C file implementing arithmetic for scalars, *and*
a 10,000 line auto-generated C file implementing arithmetic for arrays
(including 0-dim arrays), and these have different functionality and
bugs:
  https://github.com/numpy/numpy/issues/593

However, the actual goal of all this code is to make array scalars and
0-dim arrays entirely indistinguishable. Supposedly they have the same
APIs and generally behave exactly the same, modulo bugs (but surely
there can't be many of those...), and two things:

1) isinstance(scalar, np.int64) is a sorta-legitimate way to do a type
check. But isinstance(zerodim_arr, np.int64) is always false. Instead
you have to use issubdtype(zerodim_arr, np.int64). (I mean, obviously,
right?)

2) Scalars are always read-only, like regular Python scalars. 0-dim
arrays are in general writeable... unless you set them to read-only. I
think the only behavioural difference between an array scalar and a
read-only 0-dim array is that for read-only 0-dim arrays, in-place
operations raise an exception:

In [5]: scalar = np.int64(1)

# same as 'scalar = scalar + 2', i.e., creates a new object
In [6]: scalar += 2

In [7]: scalar
Out[7]: 3

In [10]: zerodim = np.array(1)

In [11]: zerodim.flags.writeable = False

In [12]: zerodim += 2
ValueError: return array is not writeable

Also, scalar indexing of ndarrays returns scalar objects. Except when
it returns a 0-dim array -- I'm pretty sure this can happen when the
moon is right, though I forget the details. ndarray subclasses? custom
dtypes? Maybe someone will remember.

Q: We could make += work on read-only arrays with, like, a 2 line fix.
So wouldn't it be simpler to throw away the tens of thousands of lines
of code used to implement scalars, and just use 0-dim arrays
everywhere instead? So like, np.array([1, 2, 3])[1] would return a
read-only 0-dim array, which acted just like the current scalar
objects in basically every way?

A: Excellent question! So ndarrays would be similar to Python strings
-- indexing an ndarray would return another ndarray, just like
indexing a string returns another string?

Q: Yeah. I mean, I remember that seemed weird when I first learned
Python, but when have you ever felt the Python was really missing a
"character" type like C has?

A: That's true, I don't think I ever have. Plus if you wanted a "real"
float/int/whatever object you could just call float() or int() or use
.item(), just like now. Can you think any problems this would cause,
though?

Q: Well, what about speed? 0-dim arrays are stupidly slow:

In [2]: x = 1.5

In [3]: zerodim = np.array(x)

In [4]: scalar = zerodim[()]

In [5]: timeit x * x
10000000 loops, best of 3: 64.2 ns per loop

In [6]: timeit scalar * scalar
1000000 loops, best of 3: 299 ns per loop

In [7]: timeit zerodim * zerodim
1000000 loops, best of 3: 1.78 us per loop

A: True!

Q: So before we could throw away that code, we'd have to make arrays faster?

A: Is that an objection?

Q: Well, maybe they're already going as fast as they possibly can be?
Part of the motivation for having array scalars in the first place was
that they could be more optimized.

A: It's true, reducing overhead might be hard! For example, with
arrays, you have to look up which ufunc inner loop to use. That
requires considering all kinds of different casts (like it has to
consider, maybe we should cast both arrays to integers and then
multiply those?), and this currently takes up about 700 ns all by
itself!

Q: It takes 700 ns to figure out that to multiply two arrays of
doubles you should use the double-multiplication loop?

A: Well, we support 24 different dtypes out-of-the-box. Caching all
the different combinations so we could skip the ufunc lookup time
would create memory overhead of nearly *600 bytes per ufunc!* So
instead we re-do it from scratch each time.

Q: Uh....

A: C'mon, that's not a question.

Q: Right, okay, how about the isinstance() thing. There are probably
people relying on isinstance(scalar, np.float64) working (even if this
is unwise) -- but if we get rid of scalars, then how could we possibly
make isinstance(zerodim_array, np.float64) work? All 0-dim arrays have
the same type -- ndarray!

A: Well, it turns out that starting in Python 2.6 -- which,
coincidentally, is now our minimum required version! -- you can make
isinstance() and issubclass() do whatever arbitrary checks you want.
Check it out:

class MetaEven(type):
    def __instancecheck__(self, obj):
        return obj % 2 == 0

class Even(object):
     __metaclass__ = MetaEven

assert not isinstance(1, Even)
assert isinstance(2, Even)

So we could just decide that isinstance(foo, some_dtype) returns True
whenever foo is an array with the given dtype, and define np.float64
to be correct dtype. (Thus also fulfilling Travis's idea of getting
rid of the distinction between scalar types and dtypes.)

Q: So basically all the dtypes, including the weird ones like
'np.integer' and 'np.number'[1], would use the standard Python
abstract base class machinery, and we could throw out all the
issubdtype/issubsctype/issctype nonsense, and just use
isinstance/issubclass everywhere instead?
[1] http://docs.scipy.org/doc/numpy/reference/arrays.scalars.html

A: Yeah.

Q: Huh. That does sound nice. I don't know. What other problems can
you think of with this scheme?

-n