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

James Hugunin jjh@mama-bear.lcs.mit.edu
Wed, 1 Nov 95 15:28:58 -0500

```> 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

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