Matlab-like Matrix module (Was Re: NumPy and Octave)
Huaiyu Zhu
hzhu at rocket.knowledgetrack.com
Tue May 30 16:04:01 EDT 2000
Matrix module that behaves like matlab
Here it is, enclosed below. Comments, patches and ideas are welcome!
>Anyway, implementing basic matrix and vector classes should be an
>afternoon job. Perhaps many people have written some and not considered
>them worth publishing!
Well, it took an afternoon plus an evening, mainly to devise enough tests).
Caution: This is a first try. The only tests are those after the line
if __name__ == "__main__". More tests are needed.
Note: I'm looking for a home for this stuff. I can't figure out how to get
onto starship. I'm also trying to find a name for this, and my current
choice is MatPy.
Plans:
1. Implement as much matlab stuff as possible.
2. Reduce the cost of wrapper class when necessary.
3. (Future) cooperation with octave.
::::::::::::::
MatPy.py
::::::::::::::
#!/usr/bin/env python
# (C) 2000 Huaiyu Zhu <huaiyu_zhu at yahoo.com>. Licence: GPL
# $Id: Matrix.py,v 1.1 2000-05-28 01:11:32-07 hzhu Exp hzhu $
"""
Matrix class with matlab-like interface.
Row vec = 1xn mat. Col vec = nx1 mat. Scalar behaves as number.
Matrix multiplication follows linear algebra rules
Todo: sub matrices and block matrices
"""
import MLab, RandomArray, LinearAlgebra, string
eps = 2.2204e-16
_time = 0
def tic(): global _time; _time = time.time()
def toc(): print "time lapse =", time.time() - _time
def short(x): return "%-6.4g" % x
def shortC(x): return "%-6.4g%+.4gj" % (x.real, x.imag)
def long(x): return "%-11.8g" % x
def longC(x): return "%-11.8g%+.8gj" % (x.real, x.imag)
def array2str(a, format=short): return " " + string.join(map(format, list(a)))
def array2strC(a): return array2str(a, format=shortC)
#------------------------------------------------------------------
class Matrix:
"Matrix class similar to matlab"
def __init__(self, data):
if type(data) == type(MLab.array([])) or type(data) == type([]):
self.data = MLab.array(data)
m = self.data.shape
if len(m) == 1: self.data.shape = (1, m[0])
elif len(m) != 2: raise ValueError, "dimension wrong"
elif isinstance(data, Matrix):
import copy
self.data = copy.copy(data.data)
else:
try:
assert data == data + 0
except:
print "data is %s" % data
raise
self.data = MLab.array([[data]])
self.shape = self.data.shape
self.typecode = self.data.typecode()
def __repr__(self):
if MLab.sum(self.shape) == 2:
return `self.data[0][0]`
else:
if self.typecode == 'd' or self.typecode == 'l':
s = map(array2str, self.data)
elif self.typecode == 'D':
s = map(array2strC, self.data)
else:
raise TypeError, "typecode =" + self.typecode
return "[" + string.join(s,'\n ') + " ]"
def __getitem__(self, i):
if type(i) is type(1):
if self.data.shape[0] == 1: return self.data[0,i]
elif self.data.shape[1] == 1: return self.data[i,0]
else: return self.data[i]
def __setitem__(self, i, item):
if type(i) is type(1):
if self.data.shape[0] == 1: self.data[0,i] = item
elif self.data.shape[1] == 1: self.data[i,0] = item
else: self.data[i] = item
def __add__(self, other):
if type(other) == type(self): return Matrix(self.data + other.data)
elif type(other) == type(1): return Matrix(self.data + other)
def __sub__(self, other):
if type(other) == type(self): return Matrix(self.data - other.data)
elif type(other) == type(1): return Matrix(self.data - other)
def __mul__(self, other):
if type(other) == type(self):
other = other.data
result = MLab.matrixmultiply(self.data, other)
if MLab.sum(result.shape) == 2:
return Scalar(result[0,0])
else:
return Matrix(result)
def __rmul__(self, other):
"There is a bug in MLab.matrixmultiply(scalar, matrix)"
if type(other) == type(self):
other = other.data
if type(other) == type(MLab.array([])):
result = MLab.matrixmultiply(other, self.data)
else:
result = other * self.data
if MLab.sum(result.shape) == 2:
return Scalar(result[0,0])
else:
return Matrix(result)
def __neg__(self): return Matrix(-self.data)
def __pow__(self, n): return Matrix(self.data ** n)
def T(self): return Matrix(MLab.transpose(self.data))
def I(self): return Matrix(LinearAlgebra.inverse(self.data))
def eig(self): return Matrix(LinearAlgebra.eigenvalues(self.data))
def real(self):
if self.typecode == "D": return Matrix(self.data.real)
else: return Matrix(self)
def imag(self):
if self.typecode == "D": return Matrix(self.data.imag)
else: return zeros(self.shape)
#------------------------------------------------------------------
class Scalar:
"Scalar class"
def __init__(self, data):
if isinstance(data, Matrix): self.data = data[0]
else: self.data = data
def __repr__(self):
if type(self.data) == type(1j): return shortC(self.data)
else: return short(self.data)
def __mul__(self, other): return self.data * other
def __rmul__(self, other): return other * self.data
def __add__(self, other): return self.data + other
def __radd__(self, other): return other + self.data
def __cmp__(self, other):
if self.data < other: return -1
elif self.data == other: return 0
else: return 1
#------------------------------------------------------------------
"New matrices"
def Mrange(*i): return Matrix(apply(range,i))
def zeros(a): return Matrix(MLab.zeros(a))
def ones(a): return Matrix(MLab.ones(a))
def eye(a): return Matrix(MLab.eye(a))
def rand(a): return Matrix(RandomArray.random(a))
"Element-wise functions"
def abs(a): return Matrix(MLab.abs(a.data))
def exp(a): return Matrix(MLab.exp(a.data))
def sin(a): return Matrix(MLab.sin(a.data))
def sinc(a): return Matrix(MLab.sinc(a.data))
def sqrt(a): return Matrix(MLab.sqrt(a.data))
def tanh(a): return Matrix(MLab.tanh(a.data))
"Matrices of the same shape"
def sort(a): return Matrix(MLab.msort(a.data))
def diff(a): return Matrix(MLab.diff(a.data))
def cumsum(a): return Matrix(MLab.cumsum(a.data))
def cumprod(a): return Matrix(MLab.cumprod(a.data))
"Reduce to row vector"
def sum(a): return Matrix(MLab.sum(a.data))
def prod(a): return Matrix(MLab.prod(a.data))
def mean(a): return Matrix(MLab.mean(a.data))
def median(a): return Matrix(MLab.median(a.data))
def std(a): return Matrix(MLab.std(a.data))
def min(a): return Matrix(MLab.min(a.data))
def max(a): return Matrix(MLab.max(a.data))
def ptp(a): return Matrix(MLab.ptp(a.data))
"Matrix functions"
def expm(a): return mfunc(MLab.exp, a)
def logm(a): return mfunc_p(MLab.log, a)
def sqrtm(a): return mfunc_p(MLab.sqrt, a)
def powm(a, n): return mfunc_p(lambda x,n=n:MLab.power(x,n), a)
"Different objects"
def diag(a): return Matrix(MLab.diag(a))
def sums(a): return Matrix(MLab.sum(MLab.sum(a.data)))
def norm(a): return Scalar(sqrt(sums(a ** 2)))
def mnorm(a): return Scalar(sqrt(a.T() * a))
def cov(a): return Matrix(MLab.cov(a.data))
#------------------------------------------------------------------
"Linear Algebra"
def solve(A, b):
AT = A.T()
return Matrix(LinearAlgebra.solve_linear_equations((AT*A).data, (AT*b).data))
def eig(a):
"""Return eigenvalues as an array, eigenvectors as a Matrix
Definition: a * u.T = u.T * v
"""
s = a.shape
assert s[0] == s[1]
v, u = MLab.eig(a.data)
return v, Matrix(u)
def svd(a):
"Return Matrix u, array x, Matrix v as svd decomposition"
u, x, v = MLab.eig(a.data)
return Matrix(u), x, Matrix(v)
def mfunc_p(f, x):
"Matrix function with positive eigenvalues"
(v, u) = eig(x)
uT = u.T()
v = MLab.absolute(v)
y = u.T() * diag(f(v-__builtins__.min(MLab.min(v),0))+eps) * uT.I()
return approx_real(y)
def mfunc(f, x):
"Matrix function"
(v, u) = eig(x)
uT = u.T()
y = uT * diag(f(v)) * uT.I()
return approx_real(y)
def approx_real(x):
if x.typecode == 'D' and norm(x.imag()) < norm(x.real())*1e-14:
return x.real()
else:
return x
#------------------------------------------------------------------
if __name__ == "__main__":
import time, sys
print "-"*40, "basic matrix"
a = Mrange(-1, 3)
print a
a [2] = 2
print a.shape, a[0,2], a[2]
print "-"*40, "arithmatic, norm"
print 1.1 * Matrix(a)*2 + 2 + a
print Scalar(-2) * a * norm(a)
print Scalar(3)*Scalar(0.99) < Scalar(1) + Scalar(2) == \
Scalar(3.) < Scalar(3)+eps*3
print "-"*40, "special matrices"
c = ones(a.shape) * 2 + zeros(a.shape) + rand(a.shape)
print c
print "-"*40, "transpose, multiplication"
d = -a.T()
print d, d.shape, a
print d*c, c*d
print "-"*40, "max, min, sum, cumsum"
X = rand((9,4))
print X
print max(X)
print min(X)
print sum(X)
print cumsum(X)
print "-"*40, "linear transform"
y = X * d
print y.T()
print "-"*40, "eigenvalues"
C = rand((5,5))
print C
print C.eig()
print "-"*40, "matrix function"
print norm(mfunc(lambda x:x, C) - C)
print "-"*40, "sqrtm, powm"
A = X.T() * X
B = sqrtm(A)
print B, norm(powm(B,2)-A)
print "-"*40, "expm, logm"
print expm(C)
print logm(expm(C)), norm(logm(expm(C)) - C)
print "-"*40, "inverse, eye"
print A.I(), norm(A.I() * A - eye(4))
print "-"*40, "solve by inverse"
tic()
A = X.T() * X
b = X.T() * y
d1 = (A.I() * b)
toc()
e = d1 - d
print norm(e), mnorm(e)
print "-"*40, "linear equation"
tic()
d2 = solve(X , y)
toc()
print norm(d2 - d), norm(d2 - d1), norm(X*d2 - y)
sys.exit()
--
Huaiyu Zhu hzhu at knowledgetrack.com
KnowledgeTrack Corporation Tel: 925 738 1907
7020 Koll Center Parkway, Suite 110 Fax: 925 738 1039
Pleasanton, CA 94566
More information about the Python-list
mailing list