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