MORE INFO - Array assigment with Numeric.py
Huaiyu Zhu
hzhu at rocket.knowledgetrack.com
Wed May 31 13:31:56 EDT 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