[PYTHON MATRIX-SIG] Some notes on the matrix class

[James Hugunin]

> Date: Wed, 1 Nov 1995 10:56:19 +0800
> From: "P. Dubois" <dubois@kristen.llnl.gov>
>
> This is a very nice facility already and has an amazing amount of
> functionality.

Thanks for the compliment and the feedback.

>    1. A function equal to max(m) - min(m), often used as a reduction operator
>       in higher dimensions. We call it ptp (peak-to-peak).

I'll add the following to the standard functions
def ptp(m): return max(m)-min(m)

where within Matrix.py, max(m) <--> maximum.reduce(m)

>    2. ranf(m) fills m with uniform random numbers on (0.0, 1.0)

I just added a function for this, (I needed to initialize a weight matrix  
for a NN) it has the form, m = rand(d1,...,dn).  where m will be a new  
matrix with the given dimensions.  It uses whrandom.py for getting its  
random numbers so is slow.  If somebody with a good random number generator  
in C would like to add this to matrixmodule.c as a new instantiation method  
in C, I'd use it.

>    3. A function reshape(x,n,m,...) equal to a copy of x reshaped as (n,m,...)

Already in my working copy of Matrix.py (requested by Konrad Hinsen).  I  
use the syntax reshape(x, (d1,...dn)) so that reshape(a, b.shape) works.  I  
don't like the idea of automatically making copies of the matrix in calls to  
reshape, but I have added the method m.copy(n=1) to make an arbitray number  
of copies of a matrix (this is much like the multiply operator for lists).   
(I also added a concat method noticing that I was missing this other  
functionality of lists).

>    4. The existing ones(n,m) will return a matrix of shape (n,m)
>       filled with 1.0; equally useful would be one where the i,j element is 
>       1.0 iff i==j, 0.0 otherwise.

This is now given by the eye(n,m,k=0) function (stolen from matlab), where  
the i,j element is 1.0 iff i-j==k, 0.0 otherwise.  This is easy enough to  
contruct yourself out of outer products (which is what eye does).

ie.
identity = equal.outer(mrange(n), mrange(m))

Aren't outer products cool?

> b. compress(condition, x) = those elements of x corresponding to those
>    elements of condition that are "true". condition and x must have
>    the same length (or x scalar), result is a one-dimensional vector with
>    length equal to the number of true elements in condition. Naturally this 
>    is most useful if x < y returns a vector instead of a scalar.
>    Most used form however is with a condition of the form x > 0. or some
>    similar comparison to a scalar, which are expressible already in Python 
>    if we take "true" in this sense to be the normal Python meaning.

This can be expressed as reshape(x, None)[condition.nonzero()].

My personal take on this is that if condition is not a vector, then this  
function is meaningless.  If it is a vector, it should be able to be applied  
to any dimension of x.

I'll add the following to Matrix.py to support this:

def compress(condition, m, dimension=-1):
	if dimension < 0: dimension = len(m.shape)+dimension
	if len(condition.shape)!=1 or len(condition)!=m.shape[dimension]:
 		raise ValueError
	index = [All]*len(m.shape)
	index[dimension] = condition.nonzero()
	print index
	return m[index]

BTW - while x < y returns a (meaningless) scalar, Matrix.less(x, y) returns  
a matrix of booleans.

There has been a small amount of talk about the possibilty of having  
"boolean types" in addition to number, sequence, and mapping types that  
would support >, <, ==, etc.  After writing a lot of code with matrices, I  
feel that many things could be expressed a lot more elegantly if this was  
allowed.  Unfortunately, Guido appears to be skeptical of this proposal, so  
don't hold your breath.

I've run into particular collisions with and, or, and not.  These functions  
when applied to matrices are renamed to the ugly "Matrix.andBoolean", etc.   
The problem is that the former are reserved words in python, so you can't  
create functions with those names (even in the "safe" context of a module).

> c. My Basis users love:
>    where(condition,x,y) is shaped like condition and has elements of x and
>    y where condition is respectively true or false; we allow broadcast of
>    scalars x or y, so some thought about ranks should reveal the equivalent
>    notion here. Often used in the form where(z, x/z, 1.e99); in this case
>    we rely on a convention that something/0.0 = Infinity, where Infinity is
>    some variable the users can set to suit them, but anyway the idea is that
>    the computation x/z goes at compiled speed and does not cause an exception.

I had something like this in my original matrix proposal, but was advised  
that it could be expressed as without any extra C-code:

def where(condition, x, y):
	c = notEqual(condition, 0)
	return c*x+(1-c)*y

where the I am relying on the fact that all comparision operators return 1  
for true.

This is in fact about 1/5th the speed of writing a where function in C, so  
it should be dropped to compiled code at some time.  I'll add this to the  
standard set of functions in Matrix.py.

>    We haven't found this perfect and of course would prefer that the
>    computation is simply not done at components where z is false. Given
>    the genius you guys have already displayed ... (Yes, you can do this
>    with a for loop, this is just a speed-freak need).

I'm a speed freak too, if you ever have to use a for loop to iterate over  
the contents of a matrix, feel free to report it as a bug in the matrix  
package.
--
You raise another interesting issue here that I've completely ignored in  
the matrix object so far.  Numeric exceptions are not handled at all; in  
return, matrix math executes at the speed of compiled code.  However, it  
wouldn't be too hard to add functions like safeDivide to the ofunc  
collection that would signal exceptions if someone considers this important.

I don't even have the hack that x/0.0 = some user settable Infinity.  x/0.0  
is set to the IEEE value for Infinity.

> d. One curiosity is the notation something[(a,b,c)]. This is required
>    because something[a,b,c] is a syntax error. But perhaps Python could
>    allow the latter interpreted as the former.

Guido beat me to this one.  I'd like to add that if anybody decides to go  
rummaging around in Grammer and compile.c, Guido also mentioned that he'd be  
willing to let a**b <--> pow(a, b) into the syntax.  This also seems like  
an easy hack for someone with the right know-how, and would be a nice bonus  
for numerical code.

There have been so many things added to the matrix object this week, that  
it seems I should get a new version out to those interested testers (so you  
won't need to ask me for features that are already there).  I'll put out a  
new object on Friday (if all goes well).  It will have all of the above  
additions (and many more), slightly different syntax for specifying the rank  
of an operator, working complex numbers, possibly working generic Python  
Objects as matrix elements and almost all of the memory leaks that I  
identified as clearly belonging to my code squashed (thanks to Purify).   
There are still a few memory leaks that appear to be in the python  
internals, and not due to my code (I think).  I'll check those out when I  
have more time.

-Jim

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

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