[Numpy-discussion] A case for rank-0 arrays

[Sasha <nd...@ma...>]

I have reviewed mailing list discussions of rank-0 arrays vs. scalars
and I concluded that the current implementation that contains both is
(almost) correct.  I will address the "almost" part with a concrete
proposal at the end of this post (search for PROPOSALS if you are only
interested in the practical part).

The main criticism of supporting both scalars and rank-0 arrays is
that it is "unpythonic" in the sense that it provides two almost
equivalent ways to achieve the same result.  However, I am now
convinced that this is the case where practicality beats purity.

If you take the "one way" rule to it's logical conclusion, you will
find that once your language has functions, it does not need numbers
or any other data type because they all can be represented by
functions (see http://en.wikipedia.org/wiki/Church_numeral).   Another
example of core python violating the "one way rule" is the presence of
scalars and length-1 tuples.  In S+, for example, scalars are
represented by single element lists.

The situation with ndarrays is somewhat similar. A rank-N array is
very similar to a function with N arguments, where each argument has a
finite domain (i-th domain of a is range(a.shape[i])).  A rank-0 array
is just a function with no arguments and as such it is quite different
from a scalar.  Just as a function with no arguments cannot be
replaced by a constant in the case when a value returned may change
during the run of the program, rank-0 array cannot be replaced by an
array scalar because it is mutable.  (See
http://projects.scipy.org/scipy/numpy/wiki/ZeroRankArray for use
cases).

Rather than trying to hide rank-0 arrays from the end-user and treat
it as an implementation artifact, I believe numpy should emphasize the
difference between rank-0 arrays and scalars and have clear rules on
when to use what.

PROPOSALS
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D

Here are three suggestions:

1. Probably the most controversial question is what getitem should
return. I believe that most of the confusion comes from the fact that
the same syntax implements two different operations: indexing and
projection (for the lack of better name).  Using the analogy between
ndarrays and functions, indexing is just the application of the
function to its arguments and projection is the function projection
((f, x) -> lambda (*args): f(x, *args)).

The problem is that the same syntax results in different operations
depending on the rank of the array.

Let

>>> x =3D ones((2,2))
>>> y =3D ones(2)

then x[1] is projection and type(x[1]) is ndarray, but y[1] is
indexing and type(y[1]) is int32.  Similarly, y[1,...] is indexing,
while x[1,...] is projection.

I propose to change numpy rules so that if ellipsis is present inside
[], the operation is always projection and both y[1,...] and
x[1,1,...] return zero-rank arrays.  Note that I have previously
rejected Francesc's idea that x[...] and x[()] should have different
meaning for zero-rank arrays.  I was wrong.

2. Another source of ambiguity is the various "reduce" operations such
as sum or max.  Using the previous example, type(x.sum(axis=3D0)) is
ndarray, but type(y.sum(axis=3D0)) is int32.  I propose two changes:

   a. Make x.sum(axis)  return ndarray unless axis is None, making
type(y.sum(axis=3D0)) is ndarray true in the example.

   b. Allow axis to be a sequence of ints and make
x.sum(axis=3Drange(rank(x))) return rank-0 array according to the rule
2.a above.

   c. Make x.sum() raise an error for rank-0 arrays and scalars, but
allow x.sum(axis=3D()) to return x.  This will make numpy sum consistent
with the built-in sum that does not work on scalars.

3. This is a really small change currently

>>> empty(())
array(0)

but

>>> ndarray(())
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
ValueError: need to give a valid shape as the first argument


I propose to make shape=3D() valid in ndarray constructor.