[PYTHON MATRIX-SIG] Mutli-dimensional indexing and other comments

[Chris Chase S1A]

This is a long note.

I have tried to go back and read through the various matrix/array
proposals from the start of this list.  I have some comments on
various proposals and opinions.  I also would like to make some
suggestions about multi-dimensional indexing. 

----

Hinsen Konrad's proposal:

     - Like arrays in J, my arrays are immutable, i.e. there is no
       provision for changing individual elements.

I find that changinng individual elements a necessity in almost all
the computations that I do.  For example, this kind of operation is
pervasive in imaging applications (e.g. segmentation and
region-of-interest manipulations) and block-matrix linear algebra
computations.  For me, setting elements is a must.

-----

Array functions:

Hinsen Konrad supports calculator sytle array usage, i.e. interactive
computation:

     I am typing things like

       sqrt(Array("file1"))*Array("file2)

     which I prefer a lot to

       Array.sqrt(Array.Array("file1"))*Array.Array("file2)

I definitely prefer the first sytle for interactive computations.

I also like Konrad's array function prototype that allows rank
specification, e.g. product.over[1](b).  Guido did not like the rank
specification via an index argument.  Could the rank be overriden
using a keyword argument instead of an index argument,
e.g. product.over(b,rank=1)?  Since keywords are commonly used to
override default arguments this would seem a natural possibility.

------

Heirarchical (list-like) indexing:

Jim Fulton's proposal thinks of multi-dimensional arrays as being
heirarchical so that a[i] returns a N-1 dimensioned array.  This
allows him to use list style indexing, e.g. a[i][j].  This approach
would seem to rule out flatten indexing (one-dimensional indexing) and
multi-dimensional slicing, e.g. would a[2:9][1:5] just end up being
equivalent to a[3:7]?  Perhaps that is what is wanted.  If the
heirarchical approach is used then a different mechanism would be
necessary for more general multi-dimensional indexing.  I could see
that heirarchical indexing like list indexing, would be useful.  My
opinion is that arbitrary multi-dimensional slicing is more useful for
the homogeneous arrays we are discussing (they are not lists).
Perhaps both could be supported.

----

I am confused about proposals for array slices to be references.  In
Jim's proposal, assignment is by copy and access is by reference.
Then am I correct in understanding that for b=a[i] (a is
two-dimensional) changing elements of b will not change a?

-----

General Multi-dimensional indexing:

In general, I have not been able to keep straight the various
proposals.  It seems that most issues have been concerned with
implementation issues or with limiting extensions so they will not
break the current language.

Instead of diving into the implementation questions I would like to
present my views on multi-dimensional indexing for arrays.  Then I
would like to know how the current proposals fit into this view.

I have worked with or investigated a number of interactive
array-oriented computational languages: IDL, Matlab, Tela, scilab,
Octave, rlab, Mathematica, APL.  I have seen the following indexing
concepts all of which are very useful and completely generalize array
access.

Like Konrad, I will use the term "rank" to refer to the number of
dimensions of an array and "shape" to refer to the ordered list of
dimension sizes of an array.

The arrays are homogeneous and can be viewed as one-dimensional (as in
C or Fortran) or multi-dimensional.

index vector: a scalar, slice (with optional stride), or array.  An
array index could be our new array (matrix) object, tuple or list.

Let A be a multi-dimensional array (rank is 3 in the examples) to be
indexed and let a1, a2, a3 be arbitrary index vectors.

one-dimensional indexing: 

a) flattened indexing.  Takes a single index vector and the array is
   viewed like contiguously-stored arrays in Fortran or C (one or the
   other, but consistent in that always the first or last dimension
   changes most rapidly.)  The shape of the result is the same as the
   shape of the index vector.

b (not really 1-D?) heirachical indexing (Jim Fulton's list-like indexing): 
    Mentioned above.  The index vector is used only for the first
    dimension.  I have never used this type of indexing for
    homogeneous numeric arrays (because it is a special case of the
    product indexing), but I see that it could be useful when viewing
    the array as a list.  For a scalar index the result has rank one
    less than A.  If a non-scalar is used as an index vector then the
    index vector would be treated in flattened form and the result
    would have have the same rank as A.

Multi-dimensional indexing:

   The number of index vectors is equal to the rank A with one index
   vector for each dimension.  There are two types of indexing that I
   have seen:

a) product indexing: You might call this arbitrary slice indexing.  I
   sometimes call this Cartesian product indexing because all
   combinations of the index vector elements are used for indexing.
   All index vectors are used in flattened form.  The result has the
   same rank but the size of each dimension is equal to the length of
   the corresponding index vector.  The elements of the result are
   taken from all possible ordered combinations of the index vector
   elements, e.g. the result element at index i,j,k is taken from A at
   index a1(i),a2(j),a3(k).

b) mapped indexing (as in Tela): The index vectors all have the same
   shape and are used in flattened form.  Indexes into array A are
   generated from ordered elementwise grouping, i.e. the result at 1-D
   index i is taken from A at index a1(i),a2(i),a3(i).  The shape of
   the result is the same as the shape of the index vectors.
   
Selection indexing (as in Octave or APL): 
  A type of 1-D indexing where the index vector is a {0,1} vector that
  has the same number of elements as A.  The result contains only
  those elements of A where the corresponding index vector is
  non-zero.

  Rather than support this directly, many languages, e.g. IDL and
  Matlab, support this indirectly with a "where" or "find" function.
  where(A) returns a rank 1 array containing in increasing order the
  one-dimensional indexes where A is nonzero.  The result of where(A)
  can then be used for 1-D indexing.  This is more powerful then
  supporting selection directly.

Insertion indexing (as in IDL):
  Used in setting blocks of items.  Takes a 1-D or multi-dimensional
  scalar index that specifies the starting position for the object
  insertion, overwriting existing items.  For a 1-D scalar index the
  object inserted is used in flattened form.  For a multi-dim index
  the object inserted must have the same rank as A.  When the inserted
  object would extend beyond the dimension bounds of A there are
  several possible behaviors: signal error, index wrap-around, or
  truncation of A to fit.


In most of the languages that I have used, the syntax A[] is used for
both 1-D and multi-dimensional product/slicing indexing.  Sometimes,
1-D indexing is used when only a single index vector is given.  Konrad
suggests that this is just syntactic sugar in place of function calls
like "take" and "ravel".  But it is extremely useful for writing
clear, readable code (even APL acknowledges this by supporting "[]"
syntax for indexing).

A natural solution syntactically would be:
1. "[]" subscripting to support 1-D indexing and product indexing
   depending on the number of index vectors given.
2. allowing ":" slice notation for index vectors.

Mapped indexing, heirarchical indexing, insertion and where/selection
could be access methods.

It does not seem possible without major internal Python changes to
implement both 1) and 2) because of the ingrained and non-extensible
1-D indexing built-in to Python.

Someone (?) suggested that arbitrary (product) indexing did not seem
useful and slicing was sufficient.  For my personal interests, product
indexing is necessary for a huge number of applications.  I use mapped
indexing less often, but it is similarly necessary in many
applications.  Of course all indexing can be converted manually into
1-D indexing, but this would make an array extension almost unuseable
for interactive use.

------

Question: Do any of the proposals completely support these types of
indexing in some form?

The closest I have seen to supporting product indexing is
James Hugunin proposal of using a mapping type:

     M[(range(1,3),range(2,4))] -> [[13,14],[23,24]]

The idea suggests wrapping the index arguments of [] into a tuple.
However, this would not allow both 1-D and multi-dimensional indexing
for the same array without possible ambiguity.

To avoid the extra parentheses Guido suggests:
 
    - Allowing M[i, j] for (multidimensional) sequence types would
      also mean that D[i, j] would be equivalent to D[(i, j)] for
      dictionaries.

But Jim Fulton:

     I see no reason to support M[i,j] for arbitrary sequence types.  I'd
     say that if a type wants to support multiple arguments to [], then it
     should provide mapping behavior and have the mapping implementation
     sniff for either an integer or a tuple argument and do the right
     thing.  

     I am *very much* against a language change to support this.

One complex solution that preserves the current Python language:

Treat M[i,j] as multiple dimension indexing and M[(i,j)] as 1-D
indexing.  Think of M[i,j] as a function call with 2 arguments and
M[(i,j)] as a function call with a single argument.  As specified by
the current language, M[(i,j)] is a mapping type, but (i,j) is treated
as a one-dimensional index vector selecting items i and j.  On the
other hand, when multiple index arguments are used, e.g. M[i,j], a
different multi-dimension index method taking a variable length
argument list is called, e.g. __msetitem__ or __mgetitem__.  This
allows both 1-D and multi-dim indexing for the same object that
_preserves_ existing language behavior.  

This also overcomes the discussion of about bundling the index
arguments into a tuple, e.g. a[1,], which I think could be a source of
errors and confusion.

Of course, this solution requires large changes to Python internals.

------------- 

Multi-dimensional slicing:

If possible I would prefer slice notation over range().

Using range() is certainly not as clean as ":".  ":" is more than
syntaic sugar since it will use dimension length information for the
object that is not available to range().  With range() how do I
specify the entire dimension or dimension length for the upper bound?
For large array expressions the repetitive range() can generate a lot
of clutter.

For example:

a[1:3,:,5:] == a[range(1:3),range(0,shape(a)[1])),range(5,shape(a)[2])]

The left side is much more readable at a glance.  Then imagine if this
was part of a larger expression.

Of course this could have been made a little simpler by storing the
shape in an auxillary variable:

s = shape(a)
then:
a[1:3,:,5:] == a[range(1:3),range(0,s[1])),range(5,len(s[2]))]

To support slice ":" for multi-dimensions Guido suggests an expression
like (2:3).  But Jim Fulton replies:

     Can't be, (2:3) is not a valid expression, so it can't yield a valid
     element of the tuple.

Additionally, allowing (2:3) expressions in place of range() expressions
would break old slice usage.

A possible complex solution to allow ":" multi-dimensional expressions
with backwards compatibility:

To make slice expressions valid there would have to be some kind of
built-in slice type.  For example, it could just be a tuple subclass
with the only difference being a type name of "slice".  When used for
multi-dim indexing the indexing function would have to support
scalars, arrays or slices for the indexes by checking the type of the
index.  For backwards compatibility, when using one-dimensional
indexing M[2:3] would call __getslice__(2,3) by unpacking the slice
tuple.

When using slices for multi-dimensional indexes, supporting negative
upper bounds would have to be a built-in dimension length function so
that the dimension upper bound would be available, e.g. when computing
M[2,3:-1].  If negative upper bounds could be forfeited for multi-dim
slices then the dimension length function would not be necessary ("None"
could be use for the upper bound in "2:").

If such additions to Python are just too much work, then I think Jim's
range() workaround is the next best thing.

------

The solutions above to allow ":" and 1-D plus multi-dimensional
indexing are not simple to implement.  They would require major
internal changes to the compiler in addition to other areas.  But they
would _preserve_ backwards compatibility with the current language.
They are complex because the original language was not designed with
extensibility to completely general multi-dimensional indexing in
mind.

General opinion about possible Python array/matrix extensions: 

Trying to implement multi-dimensional arrays without any changes to
Python internals could end up being cumbersome, inelegant workarounds
that may be confusing and discouraging to end users (especially those
coming from numeric array languages like APL, Matlab, IDL, Tela, etc).
The extensions would most likely not be completely general.  Why
should we limit the hamstring the capabilities of an array extension
just to avoid internal Python changes?  Additionally, while they might
be useable in scripts, they will be inefficient for extensive
interactive use.

Chris Chase

=================
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
=================