[Numpy-discussion] missing array type
oliphant.travis at ieee.org
Mon Feb 27 22:02:04 CST 2006
>On 2/27/06, Travis Oliphant <oliphant at ee.byu.edu> wrote:
>>.... I think 0-stride arrays are acceptable (I
>>think you can make them now but you have to provide your own memory,
>Not really. Ndarray constructor has never allowed zeros in strides. It
>was possible to set strides to a tuple containing zeros after
>construction in some cases. I've changed that in r2054
>zero strides are not allowed.
Ah, right. It was only possible to do it in C-code. But, it is
possible to do it in C-code.
Since Colin has expressed some reservations, it's probably a good idea
to continue the discussion before doing anything.
One issue I have with zero-stride arrays is that essentially is what
broadcasting is all about. Recently there has been a discussion about
bringing repmat functionality over. The repmat function is used in some
array languages largely because there is no such thing as broadcasting
and arrays are not ND.
Perhaps what is desired instead is rather than play games with indexing
on a two dimensional array you simply define the appropriate
4-dimensional array. Currently you can define the size of the new
dimensions to be 1 and they will act like 0-strided arrays when you
operate with other arrays of any desired shape. Zero-strided arrays are
actually quite fundamental to the notion of broadcasting.
I've been annoyed for several years that the idea of linear operators is
constrained in most libraries to 2 dimensions. There are many times I
want to find an inverse of an operator that is most naturally expressed
with 6 dimensions. I have to myself play games with indexing to give
the computer a matrix it can understand. Why is that? I think the
computer should be doing the work of raveling and unraveling those
indices for me. I think we have the opportunity in NumPy/SciPy to be
much more general. A tensor class that handles the "index-raveling"
that so many people have become conditioned to think is necessary could
and should be handled by the class. If you've ever written
finite-element code you should know exactly what I mean.
[End Soap Box]
On the one hand, we could just tell people to try and use broadcasting
so that zero-strided arrays show up in Python in definitive ways. On
the other we can just expose the power of zero-strided arrays to Python
and let people come up with their own rules. I lean toward giving
people the capability and letting them show me what it can do.
The only thing controversial, I think is the behavior of outputs on
ufuncs for strided arrays. Currently ufunc outputs always have full
strides unless an output array is given. Changing this default
behavior would require some justification (not to mention some code
tweaking). I'm not immediately inclined to change it even if
zero-strided arrays are allowed to be created from Python.
>In order to make zero stride arrays really useful, they should survive
>transformation by ufunc. With my patch if x is a zero-stride array of
>length N, then exp(x) is a regular array and exp is called N times to
>compute the result. That would be a much bigger project. As a first
>step, I would just disallow using zero-stride arrays as output to
>avoid problems with inplace operations.
Hmm.. Could you show us again what you mean by these problems and the
better behavior that could happen if ufuncs were changed?
More information about the Numpy-discussion