[Numpy-discussion] Rank-0 arrays - reprise
Sat Jan 5 16:10:04 CST 2013
On Sat, Jan 5, 2013 at 3:31 PM, Nathaniel Smith <firstname.lastname@example.org> wrote:
> On 5 Jan 2013 12:16, "Matthew Brett" <email@example.com> wrote:
>> 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:
>> Quoting from an email from Travis in that thread, replying to an email
>> from Tim Hochberg:
>> > 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.
>> 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.
>> There seemed then to be some impetus to remove rank-0 arrays and
>> replace them with Python scalar types with the various numpy
>> precisions :
>> Travis' recent email hints at something that seems similar, but I
>> don't understand what he means:
>> 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.
>> 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 : a = np.array([1, 2, 3])
> In : a
> Out: 1
> In : type(a)
> Out: 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 : int64_dtype = np.dtype(np.int64)
> In : int64_dtype
> Out: dtype('int64')
> In : type(int64_dtype)
> Out: numpy.dtype
> In : int64_dtype.type
> Out: 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 : int64_arr = np.array(1)
> In : int64_arr
> Out: array(1)
> In : int64_arr.dtype
> Out: 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
> 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,
> 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 : scalar = np.int64(1)
> # same as 'scalar = scalar + 2', i.e., creates a new object
> In : scalar += 2
> In : scalar
> Out: 3
> In : zerodim = np.array(1)
> In : zerodim.flags.writeable = False
> In : 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]) 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,
> Q: Well, what about speed? 0-dim arrays are stupidly slow:
> In : x = 1.5
> In : zerodim = np.array(x)
> In : scalar = zerodim[()]
> In : timeit x * x
> 10000000 loops, best of 3: 64.2 ns per loop
> In : timeit scalar * scalar
> 1000000 loops, best of 3: 299 ns per loop
> In : 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
> 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', 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?
>  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?
Thanks for the entertaining explanation.
I don't think 0-dim array being slow is such a big drawback. I would
be really surprised if there was no way to make them faster, and
having unspecified, nearly duplicated type handling code in multiple
places is likely one reason why nobody took time to really make them
Regarding ufunc combination caching, couldn't we do the caching on
demand ? I am not sure how you arrived at a 600 bytes per ufunc, but
in many real world use cases, I would suspect only a few combinations
would be used.
Scalar arrays are ones of the most esoteric feature of numpy, and a
fairly complex one in terms of implementation. Getting rid of it would
be a net plus on that side. Of course, there is the issue of backward
compatibility, whose extend is hard to assess.
More information about the NumPy-Discussion