[Numpy-discussion] missing array type

[Travis Oliphant]

Sasha wrote:

>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,
>>right?)
>>
>>    
>>
>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
><http://projects.scipy.org/scipy/numpy/changeset/2054>.  Currently
>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.

[Soap Box]
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? 


-Travis