[PYTHON MATRIX-SIG] Let's get going

Jim Fulton, U.S. Geological Survey jfulton@usgs.gov
Mon, 11 Sep 1995 19:04:39 -0400


Well, now that we're all here, let's get started. :-)

As you know, the purpose of this SIG is to develop a proposal (and
possible an implementation) for a new Python matrix type.  I'd like to
use Jim Hugunin's proposal as a starting point. I give my response to
his proposal below.  Before beginning this, however, I thought I'd
summarize my involvement to date, so you know where I'm coming from.

At the python workshop, I presented some work in progress on a Fortran
(77) Interface Description Language (FIDL) for Python.  This tool
allows one to very quickly create Python modules that call Fortran
libraries given brief high-level descriptions of the library routines
to be called.  The primary data structure in Fortran 77 is the fortran
array, which is a multi-dimensional block of homogenous data elements.
Of course, many interesting C routines, especially in mathematic and
scientific domains, use this same data structure.  To facilitate
interfacing to these sorts of routines, I wanted an an efficient
Python implementation of the same sort of data structure.

In particular:

   o I wanted a python built-in data type containing a data structure
     that I could pass directly to Fortran or C routines.  I did not
     want to have to do alot of data conversion or copying on each
     call.  For example, if a Fortran routine returns a matrix, then
     it should be possible to pass this matrix to another Fortran
     routine without reallocating a data structure or converting the
     entire matrix. 

   o I wanted multi-dimensional access as would be familiar to a
     python programmer. That is, if I have two-dimensional arrays, a
     and b, then:

       a[i][j] = b[k][l]

     should work as expected and should work efficiently.

     In general, given an n-dimensional matrix, m, 

        if n > 1, then m[i] should return an n-1 dimensional matrix,
	and 
	if n = 1, then m[i] should return a scalar.

     Note that for "m[i][j] = spam" to work correctly, assignment to
     the jth element of m[i] must be reflected in m.

   o I need to support *all* data types supported by Fortran 77,
     including strings. Not only do I have to interface to
     mathematical routines, but I need to interface to "legacy"
     systems (Some of which are still being written :).

What I came up with was a fairly simple implementation of a matrix
type that was a pure container type.  My matrix implementation
currently provides no numeric behavior, although a matrix type cries
out for numeric behavior.

My matrix implementation includes the following features:

  o Matrix objects have an internal data pointer that points to a
    homogenous block of memory.  This pointer (or some offset from it,
    of course) can be passed directly to C or Fortran.

  o The internal data structure may be shared among multiple matrix
    objects.  In particular, a __getitem__ from a n-dimensional
    matrix, where n>1, returns a matrix object that shares the same
    data (at an appropriate offset) as the original matrix.  This
    shared data structure is reference counted, so if you have:

      a=Matrix([[1,2,3],[4,5,6],[7,8,9]])
      b=a[1] # b == Matrix([4,5,6])
      del a

    b is not invalidated by the third line.

    Note that if you have:

      a=Matrix([[1,2,3],[4,5,6],[7,8,9]])
      b=Matrix([11,22,33])
      a[1]=b
      a[1][1]=99

    You end up with:

      a == Matrix([[1,2,3],[11,99,33],[7,8,9]])
      b == Matrix([11,22,33])

    That is, assigning to a matrix copies data.  So assignment is by
    copy, even though access is by reference.

    (As you can see, matrices can be created from sequences of
     sequences.  Matrices can also be assigned from arbitrary
     sequences of the right dimension or level of nesting.) 

  o I chose to preserve the copy semantics of slices from lists.  So
    in: 

      a=Matrix([[1,2,3],[4,5,6],[7,8,9]])
      b=a[1][:]  # b == Matrix([4,5,6])
      b[1]=99    # does not affect a

    The third line does not change a.

So that's where I'm coming from. I think this type is generally
useful, which is why I helped start this SIG. I'm pretty open on how
and if matrices should have numeric behavior. With that, here are my
comments on Jim's proposal.


On 18 Aug 1995 17:16:55 GMT 
James Hugunin said:
> I apologize for the length of this posting, but you were warned in the  
> subject header.
> 
> There seems to be a fair amount of interest in the python community  
> concerning the addition of numeric operations to python.  My own desire is  
> to have as large a library of matrix based functions available as possible  
> (linear algebra, eigenfunctions, signal processing, statistics, etc.).  In  
> order to ensure that all of these libraries interoperate, there needs to  
> be agreement on a basic matrix object that can be used to represent arrays  
> of numbers.  The following is a proposal for such an object.  It can be  
> viewed as an extension of the current array object to higher dimensions  
> and numerical operations.
> 
> The wheels have already been set in motion to create a Matrix SIG  
> mailing-list for further discussion of this proposal, as well as other  
> python numerics related issues.  This should come on-line sometime next  
> week and serious discussion of this proposal should probably be deferred  
> until it can be done on that list.  Summaries of that discussion will be  
> posted to the main list as consensus develops.
> 
> Below I specify the behavior of a Matrix object within python.  The C (or  
> FORTRAN) API to such an object is obviously of significant importance as  
> this object is primarily intended as an interface to existing numerical  
> libraries.  I'm currently working on that proposal (at least the C part)  
> as well, but I think everyone will agree that this post is long enough as  
> it is.
> 
> Note: Parts of this proposal are stolen from the python documentation for  
> the array object, and from discussions on the net.  I'd like to  
> particularly thank Jim Fulton for letting me review (and steal ideas from)  
> his current Matrix object and for providing feedback on the first draft of  
> this proposal (though he certainly doesn't agree with everything in it,  
> even now).

We aren't that far apart now, at least not on things that matter to
me. :-)
 
> 
> Matrix Object
> 
> This defines a new object type which can efficiently represent a  
> multidimensional array of basic values (char, unsigned byte, signed byte,  
> int, short, long, float, double, complex float and complex double).   
> Matrices are defined as both a sequence type and a mapping type, and  
> usually (unless it's a matrix of char's) as a number type.  It should be  
> noted that allowing a type to be both a sequence and a number and to  
> behave properly for addition and multiplication, will require some very  
> small patches to the ceval.c code.

While I think the patches are OK, I don't think they're necessary, so
if anyone objects to this proposal on the grounds that changes are
required to ceval.c, you should not be concerned.
 
> All of the examples use the following matrices:
> Assume A = [1,2,3], B = [11,12,13], C = [[1,4,9],[16,25,36]],
> M=[[ 1, 2, 3, 4, 5],[11,12,13,14,15],[21,22,23,24,25],[31,32,33,34,35]]
> 
> 
> **Numeric operations**
> 
> "+", "-", "*", "/", "^"(as power), abs, unary"-", and unary"+" are defined  
> as the corresponding element-wise arithmetic operations.  Because these  

I'd rather see multiplication and division not be elementwise.  But I
don't feel strongly about this.

> assume element-wise correspondence between the two matrices, the operandi  
> must be of the same dimensions (however, note the dimensions-coercion  
> section below).
> 
> A+B = [12,14,16]
> A*B = [11,24,39]
> C*C = [[2,8,18],[32,50,72]]
> -A = [-1,-2,-3]
> 
> Possible suggestions for the other numeric operations:
> "%" represents matrix multiplication (or vector dot product).

If we go with elementwise interpretation for * and /, then I think
that all operators that make sense for floating point numbers should
have an elementwise interpretation, so m%spam should compute a matrix
of mods.  Perhaps one would carry this same argument to bitwise
operators, since you could have matrices of integers.

If we go with elementwise interpretation of numeric operations, then
I'm inclined to use functional, rather than matrix notation for matrix
multiplication, dot products, and so on.

> "~" represents transposition.
> Other suggestions?

Perhaps, as a compromize, there could be a special method or operator
that creates a reference copied version of a matrix that provides a
different semantics.  For example, perhaps the default matrix type
could provide standard matrix interpretation for * and /, but there
could be a method, elementwise that would allow:

   m.elementwise() + a

to do elementwise arithmentic.


IMPORTANT NOTE: Mixed-mode arithmetic between matrices and scalars
		should be defined.
 
> 
> **Sequence operations**
> 
> Concatenation and multiplication sequence operations are not defined, but  
> instead the arithmetic operations are invoked for the "+" and "*"  
> operators. 

Right.  Note if we don't need elementwise bitwise operators, then we
could concatinate with "|", which I find most natural in a matrix
context:

  spam = foo | bar

I suppose one could even repeat with >> or <<, but this is not so natural.

> Indexing by both single indices and by slices is supported.   
> All returned matrices are returned by reference rather than by value.
> len(M) is defined to return the size of the first dimension in in matrix.
> 
> 
> examples:
> len(M) -> 4
>  
> D = M[0]
> D[0] = 66
> print M[0]
> 	[66,2,3,4,5]

Right, for my matrix type.

> The semantics of returning a slice by reference is consistent with the  
> mapping operations given below, however, this is inconsistent with the  
> list objects semantics.  This is open to debate.
> 
> D = M[0:2]
> D[0][0] = 66
> print M[0]
> 	[66,2,3,4,5]

I'm not sure what you were trying to say here, but this example holds
only if M is a list of lists.  If M is one of my matrices, then M is
unchanged. Perhaps that is your point.

Note that I tried to maintain the spirit of list semantics, in which
the slice operator has copy semantics, but with a list, the things you
are copying are, themselves, references.

Anything is open to debate, but how else could this work, and satisfy
other requirements at the same time?
 
> 
> **Mapping operations**
> 
> Recent newsgroup discussions have offered some convincing arguments for  
> treating a matrix as a mapping type where sequences of ints are used as  
> the keys.  The sequence must be of length less than or equal to the numer  
> of dimensions in the matrix.  Each element of the sequence is either an  
> integer or a sequence.  If it is an integer, it returns the corresponding  
> element for that dimension, if it is a sequence then it returns all of the  
> elements given in the sequence.
> 
> ie. A[0] -> 1, A[((1,2))] -> [2,3]
> M[0] -> [1,2,3,4,5]
> M[(0,3)] -> 4
> M[((0,2),0)] -> [1,21]
> M[(range(1,3),range(2,4))] -> [[13,14],[23,24]]

A way of visualizing this is as multi-dimensional slices that let you,
for example, take a rectangular slice out of a 2-d matrix.

> In order to simplify the syntax of these references, it would be valuable  
> to change the syntax of python so that M[0,3] is equivalent to M[(0,3)].   
> This is a small change to the current indexing semantics which will give a  
> reasonable meaning to a currently illegal syntactic item.
> 
> When used as a setvalue, the corresponding elements of the given array  
> will be set to the left hand side.  All of these operations that return a  
> matrix will return it by-reference, meaning that the elements of the new  
> array will be references to the old one and changes to either array will  
> effect the other.  

Note that if one wants a reference, not a copy, one can always apply a
copy operator, [:] as one often does not with lists.

> 
> 
> **Instantiation**
> 
> The following function will be used to instantiate a matrix:
> Matrix(datasource(optional), d1, d2, ..., typecode='d')
> 
> The typecodes are as follows:
> c - char (non-numeric)
> 1 - unsigned char
> b - signed char
> h - short
> i - int
> l - long
> f - float
> d - double
> F - complex float
> D - complex double
> 
> The optional datasource can be a file object or a string object and will  
> fill the matrix with the raw data contained therein.  In this case the  
> dimensions are required.  The first dimension can be input as -1, in which  
> case the first dimension will be set to be as large as necessary to hold  
> the entire contents of the file or string.
> 
> The optional datasource can also be a sequence, in which case the  
> dimensions of the array do not need to be passed in and are instead  
> determined by the hierarchical structure of the sequence.
> 
> examples:
> Matrix([range(3),(4,5,6)], 'i') -> [[0,1,2],[4,5,6]]
> Matrix(2,3) -> [[0,0,0],[0,0,0]]
> Matrix('abcdef', 2,3, 'u') -> [[97,98,99],[100,101,102]]
> Matrix('abcdef', -1,2, 'c') -> ['ab', 'cd', 'ef']

Note that I have found it convenient to implement matrices so that 1-d
matrices of chars behave alot like strings, as I often pass these to
or get these back from Fortran routines in which they are just that.

> 
> **Coercion**
> 
> Dimension coercion:The following is a proposal to make the matrix object  
> behave as closely as possible to matrices in matlab (to make it easy to  
> convert matlab code to python), and to generalize this behavior to  
> arbitrarily high dimensional objects.
> 
> If a matrix of fewer dimensions than that required is provided, then the  
> matrix can be converted to more dimensions by making copies of itself.   
> For purposes of dimension coercion, an int or a float is considered a 0-d  
> matrix.
> 
> ie. A+1 -> [2,3,4], A+C -> [[2,6,12],[17,27,39]]
> 	B^2 -> [121,144,169] vs. B^A -> [11, 144, 2197]

I agree wrt scalars.
 
> Simliar coercion should be performed for most functions acting on  
> matrices.  Helper functions will be provided to make this easy, however  
> there is no way to enforce this.  The following is a proposed standard to  
> allow functions to be applied to higher-dimensional objects in a uniform  
> way.
> 
> If a function is defined to operate on scalars (say sqrt()), then when  
> applied to a matrix it should apply itself to every element in the matrix.
> 
> ie. sqrt(C) -> [[1,2,3],[4,5,6]]
> 
> This is generalized to higher dimensions so that if a function is provided  
> with a matrix of more dimensions that it requires, that function will  
> apply itself to each appropriately sized unit in the matrix.
> 
> ie. sum(A) -> 6 whereas sum(C) -> [14,77]

I am a bit concerned that this will lead to hard to find bugs.  This
also put's extra burden on the function implementor.  I realize (since
you explained it to me in a separate note :) your concern that this is
needed for efficiency, otherwise the map operator, or perhaps a
specialized matrix map operator could be used.

> 
> Type coercion: No type coercion will be performed.  However, it is  
> possible to have a function defined to work on multiple types (operator  
> overloading).  So sqrt might be defined on matrices of both floats and  
> doubles in which case it will work properly when given either matrix as  
> input, but will return a type error if applied to a matrix of ints.
> 
> 
> **Basic Data Items and Methods (non-numeric)**
> 
> typecode - The typecode character used to create the matrix
> 
> itemsize - The length in bytes of one matrix item in the internal  
> representation
> 
> dimensions - A sequence containing the size of the matrix along each of  
  ^^^^^^^^^^ tuple

> its dimensions
> 
> copy - Returns a copy of the matrix.  This is a good way to take an array  
> returned by reference and turn it into one returned by value.

This is not needed.  Use [:].
 
> transpose - The transpose of the matrix (not needed if "~" is used for  
> transpose)

I'm thinking that this (or ~) should return by reference.
 
> byteswap - Useful for reading data from a file written on a machine with a  
> different byte order

Hm. Does this operate in-place or return a value?

You know, there really *should* be a matrix map operator.  This should
function like the standard map operator, but it will generate a matrix
and it will probably want to be told at what level/dimension to apply
the mapping function.  In fact, this would give you arbitrary
elementwise operations given scalar functions.

> 
> typecast(t) - Returns a new matrix of the same dimensions where each item  
> is cast (using C typecasting rules) from the matrix's type to the new type  
> t.

Our constructor already gives you this:

  new_matrix=Matrix(old_matrix,new_type)
 
> tofile(f) - Write all items (as machine values) to the file object f.
> 
> tolist() - Convert the matrix to an ordinary (possibly nested) list with  
> the same items.
> 
> tostring() - Convert the array to an array of machine values and return  
> the
> string representation (the same sequence of bytes that would
> be written to a file by the tofile() method.)
> 
> Note: the append, fromlist, etc. methods are not present because it is  
> assumed that a matrix object is of constant size (there are some strong  
> efficiency reasons for this assumption but this is likely to be another  
> source of arguments).
> 
> **Basic functions taking matrix arguments (non-numeric)**
> 
> Since it isn't possible to append to a matrix (or efficient if it's  
> decided it should be possible) I thought that some sort of function for  
> producing a new matrix from a list of matrices would be useful.  So I'm  
> proposing the following function which is very similar in spirit to the  
> string.join operation.
> 
> concat(l) - l must be a list of matrices, and each matrix in l must have  
> the same number of dimensions, as well as the same size for every  
> dimension except for the first.  This will return a new matrix M whose  
> first dimension is the sum of the first dimensions of all the matrices in  
> the list, and whose data is filled from the data in all of these matrices.
> 
> ie. concat([A,B]) -> [1,2,3,11,12,13]
> 
> This is different from Matrix([A,B]) -> [[1,2,3],[11,12,13]]
> 
> 
> **Basic functions taking matrix arguments (numeric)**
> 
> It is hoped that there will be a large variety of special purpose modules  
> to do fancy numeric operations on matrices.  There should also be a basic  
> set of operations provided in the Matrix module (the module with the  
> instantiation function).  Because type coercion is not going to be  
> performed, it is important to decide both on the list of functions and on  
> the basic types they should support.  The following is a very preliminary  
> list of the functions I consider most fundamental.
> 
> rng(start, stop, step(optional), typecode='d') where all points in the  
> range must be strictly less than the stop point, according to the  
> semantics of range().  This function should possibly be added to the  
> instantiation options.
> 
> ie. rng(0,0.4) = [0,0.1,0.2,0.3]; rng(0,0.1,0.41) = [0,0.1,0.2,0.3,0.4]
> 

Hm.  Why not:

  Matrix.range(start=1,stop,step=1, typecode='d')

and

  Matrix.range([(start=1,stop,step=1), ...], typecode='d')

That is:

  call it range, and give it same semantics, and

  include a multi-dimensional range operator, mrange, that takes a
  sequence of range specs.

> 
> scalar operations (defined on floats, doubles, and all complex):
> everything in mathmodule.c
> round
> 
> scalar operations (defined on all numeric types):
> sign
> 
> vector to scalar operations (defined on all numeric types):
> max, min - Should these return the index of the max or min?
> sum, prod
> 
> vector to vector operations (defined on all numeric types):
> cumsum, cumprod


Comments anyone?




-- Jim Fulton      jfulton@usgs.gov          (703) 648-5622
                   U.S. Geological Survey, Reston VA  22092 
This message is being posted to obtain or provide technical information
relating to my duties at the U.S. Geological Survey.

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

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