[Matrix-SIG] Zeros

[Charles G Waldman]

Paul F. Dubois writes:

 > I question seriously whether Numerical should have any support for
 > non-numeric types, actually. It causes serious complications for very little
 > real gain. There is an array object in the core that might suit your needs.

I would like to see even more support for non-numeric types, for the
following reason: I could use the same NumPy code to approach problems
from either a symbolic and a numerical approach.  E.g. I would like to
create a symbolic "Polynomial" class, populate a matrix with these
objects, then call some Numeric functions to perform symbolic algebra.
Then when I'm happy with that, I can plug in numeric values and have a
fast number-crunching engine.

For instance:

>>> from Numeric import *
>>> from LinearAlgebra import determinant
>>> det=determinant
>>> class Poly:
...     def __init__(self,val):
...             self.val=val
...     def __add__(self,other):
...             return Poly("%s+%s"%(self.val,other.val))
...     def __mul__(self,other):
...             return Poly("%s*%s"%(self.val,other.val))
...     def __repr__(self):
...             return self.val
...     def __str__(self):
...             return self.val
... 
>>> a,b,c,d=map(Poly,'abcd')
>>> a
a
>>> a+b
a+b
>>> m=array([a,b,c,d],'O')
>>> m.shape=2,2
>>> m
array([[a , b ],
       [c , d ]],'O')
>>> m*m
array([[a*a , b*b ],
       [c*c , d*d ]],'O')
>>> m2=array([w,x,y,z],'O')
>>> m2.shape=2,2
>>> dot(m,m2)
array([[a*w+b*y , a*x+b*z ],
       [c*w+d*y , c*x+d*z ]],'O')


This, IMO, is a great way to test code and see that your matrix
calculations are actually doing what you want them to.  Seems to me
like the combination of unlimited Python dynamism and speed when you
need it, is something really unique to NumPy.... the possibility that
you could create arrays of rational numbers, or elements of a finite
field, etc, and compute with those is something I'd not like to see go
away.

Unfortunately, as has been pointed out many times,
 the support for type 'O' arrays is not uniform:

>>> det(m)
Traceback (innermost last):
  File "<stdin>", line 1, in ?
  File "/usr/lib/python1.5/site-packages/Numeric/LinearAlgebra.py", line 228, in determinant
    t =_commonType(a)
  File "/usr/lib/python1.5/site-packages/Numeric/LinearAlgebra.py", line 27, in _commonType
    kind = max(kind, _array_kind[t])
KeyError: O



I think the solution to this is not to eliminate type 'O' but to do
some work to make it work better.

Of course I realize whom the burden of such work rests upon...