[Numpy-discussion] untenable matrix behavior in SVN

[Travis E. Oliphant]

Charles R Harris wrote:
>
>
> On Tue, Apr 29, 2008 at 7:24 AM, Stéfan van der Walt <stefan at sun.ac.za 
> <mailto:stefan at sun.ac.za>> wrote:
>
>     Hi Charles
>
>     2008/4/29 Charles R Harris <charlesr.harris at gmail.com
>     <mailto:charlesr.harris at gmail.com>>:
>     > May I add that if I edit defmatrix.py to act like an array for
>     scalar
>     > indexing, then the following works.
>     >
>     > In [1]: a = matrix(eye(2))
>     >
>     > In [2]: array([a,a])
>     > Out[2]:
>     > array([[[ 1.,  0.],
>     >          [ 0.,  1.]],
>     >
>     >        [[ 1.,  0.],
>     >         [ 0.,  1.]]])
>     >
>     > This generates an error with the current version of matrix and,
>     frankly, I
>     > am not going to be bothered going all through the numpy c
>     sources to special
>     > case matrices to fix that. Someone else can do it if they wish.
>     There are
>     > recursive routines that expect the dimensions to decrease on
>     each call.
>
>     I'd also like to see matrices become proper hierarchical containers --
>     the question is just how to do that.  Thus far, I'm most convinced by
>     the arguments for RowVectors/Columns, which leaves us with a sane
>     model for doing linear algebra, while providing the enhancements you
>     mentioned here and in comments to another ticket.
>
>     We were thinking of raising a warning on scalar indexing for 1.1, but
>     given the above, would that be sensical?
>
>
> Hi Stefan,
>
> The numpy c routines call PySequence_GetItem(s,i) as well as ask for 
> the length (first index), so I think these should be left as they are 
> for arrays in order to guarantee that matrices are compatible with all 
> the normal array operations. This means either returning special row 
> objects that index as expected
> or returning 1D arrays. I don't think the '*' operator has these 
> problems, but in any case that is a well documented feature of matrices.

Thanks for looking in to this Chuck,

I'm quite persuaded now that a[i] should return a 1-d object for 
arrays.    In addition to the places Chuck identified,  there are at 
least 2 other places where special code was written to work-around the 
expectation that item selection returns an object with reduced 
dimensionality (a special-cased .tolist for matrices and a special-cased 
getitem in the fancy indexing code).  

As the number of special-case work-arounds grows the more I'm convinced 
the conceptualization is wrong.   So, I now believe we should change the 
a[i] for matrices to return a 1-d array.  

The only down-side I see is that a[i] != a[i,:] for matrices.   
However,  matrix(a[i]) == a[i,:], and so I'm not sure there is really a 
problem, there.  I also don't think that whatever problem may actually 
be buried in the fact that type(a[i]) != type(a[i,:]) is worse than the 
problem that several pieces of NumPy code actually expect hierarchical 
container behavior of multi-dimensional sequences.  

I don't think making the small change to have a[i] return 1-d arrays 
precludes us from that 1-d array being a bit more formalized in 1.2 as a 
RowVector should we choose that direction.   It does, however, fix 
several real bugs right now.

So, I'm

+1 on making the following change:

def __getitem__(self, index):
+         if isscalar(index):
+             return self.__array__()[index]

-Travis