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

```