Matlab vs Python (was RE: Discussion: Introducing new operators for matrix computation)

[David M. Cooke]

At some point, hzhu at localhost.localdomain (Huaiyu Zhu) wrote:

> On 17 Jul 2000 17:39:33 -0400, David M. Cooke <cookedm at physics.mcmaster.ca>
> wrote: 
> >> Consider the proposal of having two classes, Matrixwise and Elementwise,
> >> with methods to cast to each other.  I don't know why
> >> 
> >> (a.E()*b.E()).M()*(c.E()*d.E()).M()
> >> 
> >> would be more preferable to 
> >> 
> >> (a.*b)*(c.*d)
> >
> >Well, it doesn't have to be as complicated as your first example. Assume
> >the default is Matrixwise, and we get rid of some ()'s by overriding
> >__getattr__, so that we can right
> >
> >(a.E*b.E)*(c.E*d.E)
> >
> That's a good point.  Can __getattr__ be made to work only on a selected
> number of attributes?  We could use it on T, H, I, C for transpose,
> Hermitian transpose, inverse and conjugate.  But if it works on all the
> undefind members the error message may be too obscure.  Any examples?

Yep. See code below.

> BTW, you meant
> 
> (a.E*b.E).M*(c.E*d.E).M
> 
> for otherwise it is equivalent to (a.*b).*(c.*d).

Actually, I was thinking that elementwise multiplication would return
an object that does matrixwise multiplication (the idea being that
elementwise is a special case, just as .* is used for it instead of *
in Matlab).

> >But what if we assume that if one operand has .E, the other does also!
> >Then we can write
> >
> >(a.E*b)*(c.E*d)
> 
> This won't work.  It would introduce many hard-to-trace bugs.  Looking at
> 
> (a*b)*(c*d)
> 
> and you have to trace the program all the way down to figure out the
> identities of a,b,c,d.  Note that in concept the identities of the objects
> do not change.  We are changing their identities only as labels to help
> picking different operators.  I can foresee great confusion in practice with
> this approach.

I realize this. Note that you have to trace the program to tell that
a,b,c,d are matrices anyway :-). A lot of this is removed by always
having an operation return a Matrixwise object.

Anyways, here's some code that implements what I was talking about. It
uses the UserArray and Matrix classes included with NumPy:

# matclass.py
from UserArray import UserArray
from Matrix import Matrix

class Matrixwise(Matrix):
    def __init__(self, data, typecode=None):
        if isinstance(data, UserArray):
            UserArray.__init__(self, data.array, typecode=typecode)
        else:
            UserArray.__init__(self, data, typecode)

    def __getattr__(self, name):
        if name == 'E':
            return Elementwise(self)
        else:
            try:
                return self.__dict__[name]
            except KeyError:
                raise AttributeError, name

    def __mul__(self, other):
        if isinstance(other, Elementwise):
            return Elementwise.__rmul__(self, other)
        else:
            return Matrix.__mul__(self, other)

class Elementwise(UserArray):
    def __init__(self, data, typecode=None):
        if isinstance(data, UserArray):
            UserArray.__init__(self, data.array, typecode)
        else:
            UserArray.__init__(self, data, typecode)

    def __mul__(self, other):
        return Matrixwise(UserArray.__mul__(self, other).array)


Here's a session:
>>> import matclass
>>> a = matclass.Matrixwise([[2,3],[4,5]])
>>> b = matclass.Matrixwise([[1,2],[3,4]])
>>> a*b
matclass.Matrixwise([[11, 16],
       [19, 28]])
>>> a.E*b.E
matclass.Matrixwise([[ 2,  6],
       [12, 20]])
>>> a.E*b
matclass.Matrixwise([[ 2,  6],
       [12, 20]])
>>> a*b.E
matclass.Matrixwise([[ 2,  6],
       [12, 20]])

Obviously, this can be extended to include transpose, etc.

-- 
|>|\/|<
----------------------------------------------------------------------------
David M. Cooke
cookedm at mcmaster.ca