MORE INFO - Array assigment with Numeric.py

Huaiyu Zhu hzhu at rocket.knowledgetrack.com
Wed May 31 19:31:56 CEST 2000


On Wed, 31 May 2000 08:42:21 -0400, Morris, Brad [WDLN2:2W40:EXCH]
<confused at americasm01.nt.com> wrote: 
[snip]
>i.e. there appear to be three possible shapes for a vector where x
>represents the number of elements:
>1) (x,)
>2) (x,1)
>3) (1,x)
>The difference between (2) and (3) are obvious (row vs column) but I
>don't
>understand how (1) is different from (2) or why sometimes (1) is
>returned
>and other times (2) or (3). I think this is part of the problem.

Try the Matrix module enclosed below that has a Matlab-like interface.
The usage is ilusstrated by the self-tests.

I had this problem with NumPy and struggled for weeks.  It's not worth it.
The new module only took me less than one day to write and it worked
reasonably well.  I plan to put it on sourceforge as many people seem to
have similar problems.  Right now this is just a wrapper around NumPy but
someday someone might implement some parts directly in C.

Huaiyu


#!/usr/bin/env python
# (C) 2000 Huaiyu Zhu <huaiyu_zhu at yahoo.com>.  Licence: GPL
# $Id: Matrix.py,v 1.3 2000-05-30 19:02:10-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, __builtin__
eps = 2.2204e-16
mul = MLab.matrixmultiply
sol = LinearAlgebra.solve_linear_equations
_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 MatrixType: pass

class Matrix(MatrixType):
	"Matrix class similar to matlab"
	
	def __init__(self, data):
		if type(data) == type(MLab.array([])) or \
		   type(data) == type([]) or \
		   type(data) == type((1,)):
			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, "Matrix of shape %s"%`m`
		elif isinstance(data, Matrix):
			import copy
			self.data = copy.copy(data.data)
		elif isinstance(data, Scalar):
			self.data = MLab.array([[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):		return Matrix(self.data + other)
	def __radd__(self, other):		return Matrix(other + self.data)

	def __sub__(self, other):		return Matrix(self.data - other)
	def __rsub__(self, other):		return Matrix(other - self.data)

	def __mul__(self, other):
		if isinstance(other, MatrixType):		other = other.data
		result = mul(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 isinstance(other, MatrixType):		other = other.data
		if type(other) == type(MLab.array([])):
			result = mul(other, self.data)
		else:			result = other * self.data
		if MLab.sum(result.shape) == 2:			return Scalar(result[0,0])
		else:			return Matrix(result)

	def __div__(self, other):
		if isinstance(other, MatrixType):
			result = solve(other.T(), self.T()).T()
		else:
			result = Matrix(self.data / other)
		if MLab.sum(result.shape) == 2:		return Scalar(result[0,0])
		else:			return Matrix(result)

	def __rdiv__(self, other):
		if isinstance(other, MatrixType):
			result = solve(self.T(), other.T()).T()
		else:
			result = other * self.I()
		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 __add__(self, other): return self.data + other
	def __radd__(self, other): return other + self.data
	def __sub__(self, other): return self.data - other
	def __rsub__(self, other): return other - self.data
	def __mul__(self, other): return self.data * other
	def __rmul__(self, other): return other * self.data
	def __div__(self, other): return self.data / other
	def __rdiv__(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
	def __int__(self):		return int(self.data)
	def __float__(self):		return float(self.data)

#------------------------------------------------------------------

"New matrices"
def Mrange(*i):	return Matrix(apply(range,i))
def zeros(a, typecode='d'):	return Matrix(MLab.zeros(a, typecode))
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(Matrix(a).data))
def exp(a):		return Matrix(MLab.exp(Matrix(a).data))
def sin(a):		return Matrix(MLab.sin(Matrix(a).data))
def sinc(a):	return Matrix(MLab.sinc(Matrix(a).data))
def sqrt(a):	return Matrix(MLab.sqrt(Matrix(a).data))
def tanh(a):	return Matrix(MLab.tanh(Matrix(a).data))

"Matrices of the same shape"
def sort(a):	return Matrix(MLab.msort(Matrix(a).data))
def diff(a):	return Matrix(MLab.diff(Matrix(a).data))
def cumsum(a):	return Matrix(MLab.cumsum(Matrix(a).data))
def cumprod(a):	return Matrix(MLab.cumprod(Matrix(a).data))

"Reduce to row vector"
def sum(a):		return Matrix(MLab.sum(Matrix(a).data))
def prod(a):	return Matrix(MLab.prod(Matrix(a).data))
def mean(a):	return Matrix(MLab.mean(Matrix(a).data))
def median(a):	return Matrix(MLab.median(Matrix(a).data))
def std(a):		return Matrix(MLab.std(Matrix(a).data))
def min(a):		return Matrix(MLab.min(Matrix(a).data))
def max(a):		return Matrix(MLab.max(Matrix(a).data))
def ptp(a):		return Matrix(MLab.ptp(Matrix(a).data))

"Matrix functions"
def expm(a):	return mfunc(MLab.exp, a)
def logm(a):	return mfuncC(MLab.log, a)
def sqrtm(a):	return mfuncC(MLab.sqrt, a)
def powm(a, n):	return mfuncC(lambda x,n=n:MLab.power(x,n), a)

"Different objects"
def diag(a):	return Matrix(MLab.diag(a))
def sums(a):	return Scalar(MLab.sum(MLab.sum(Matrix(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(Matrix(a).data))

#------------------------------------------------------------------
"Linear Algebra"

def solve(A, b):
	AT = A.T()
	return Matrix(sol((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)
	if v.typecode() == 'D':
		if norm(Matrix(v.imag)) > norm(Matrix(v.real))*1e-14:
			raise TypeError, "%s only applicable to real numbers" % `f`
		else:
			v = v.real
	print v, v.typecode()
	
	if MLab.min(v) + eps < 0:
		raise TypeError,  "%s only applicable to positive numbers" % `f`
	else:
		uT = u.T()
		y = u.T() * diag(f(__builtin__.max(v,0))) * 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 mfuncC(f, x):
	"Matrix function with possibly complex eigenvalues"
	(v, u) = eig(x)
	uT = u.T()
	y = uT * diag(f(v+0j)) * 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 vector"
	a = Mrange(-1, 3)
	print a
	a [2] = 2
	print a.shape,  a[0,2], a[2]

	print "-"*40, "basic matrix"
	b = zeros((2,2))
	b [1,0] = 1; b[0,1] = 2
	print b
	print b.shape,  b[0,1], b[1,0]
	print 1 + b + 1
	print 2 - b - 2
	print 3 * b * 3
	print b / 4
	print 4 / b * b

	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