#!/usr/bin/env python
# (C) 2000 Huaiyu Zhu <hzhu@users.sourceforge.net>. Licence: GPL
# $Id: Matrix.py,v 1.17 2001/08/26 12:40:43 hzhu Exp $
"""
Matrix class with matlab-like interface. Related functions.
Row vec = 1xn mat. Col vec = nx1 mat.
1x1 mat is different from a number
Matrix multiplication follows linear algebra rules
"""
from MatPy import MatrixType
from MatPy import NaN,nan,Inf,inf
import Numeric, MLab, RandomArray, LinearAlgebra, string, __builtin__
from Numeric import Int
_eps = 2.2204e-16
_eps_approx = 1e-13
_mul = MLab.matrixmultiply
_array = MLab.array
_Axis = MLab.NewAxis
_abs = MLab.absolute
_sum = MLab.sum
_sqrt = MLab.sqrt
_pow = MLab.power
_min = MLab.min
_max = MLab.max
_minimum = Numeric.minimum
_maximum = Numeric.maximum
_trsp = MLab.transpose
_conj = MLab.conjugate
_diag = MLab.diag
_sol = LinearAlgebra.solve_linear_equations
_inv = LinearAlgebra.inverse
_eig = LinearAlgebra.eigenvalues
_ones = MLab.ones
_zeros = MLab.zeros
_eye = MLab.eye
_rand = RandomArray.random
_gt = Numeric.greater
_ge = Numeric.greater_equal
_lt = Numeric.less
_le = Numeric.less_equal
_eq = Numeric.equal
#------------------------------------------------------------------
class Matrix(MatrixType):
"class Matrix : main matrix class"
def __init__(self, data):
"__init__: Initialize to a new Matrix of given data"
if isinstance(data, Matrix):
self.data = data.data
else:
self.data = _array(data)
m = self.data.shape
if len(m) == 2: pass
elif len(m) == 1: self.data.shape = (1, m[0])
elif len(m) == 0: self.data.shape = (1, 1)
else : raise ValueError, "Matrix of shape %s"%`m`
self.shape = self.data.shape
self.typecode = self.data.typecode()
def copy(self):
"A new matrix with a copy of self.data"
import copy
data = copy.copy(self.data)
return Matrix(data)
#------------------------------------------------------------------
# Methods that look like attributes
def __getattr__(self, name):
"Short cut to methods T H C I eig real imag"
if name == 'T': return Matrix(_trsp(self.data))
elif name == 'H': return Matrix(_conj(_trsp(self.data)))
elif name == 'C': return Matrix(_conj(self.data))
elif name == 'I': return Matrix(_inv(self.data))
elif name == 'eig': return Matrix(_eig(self.data))
elif name == 'real':
if self.typecode == "D": return Matrix(self.data.real)
else: return self
elif name == 'imag':
if self.typecode == "D": return Matrix(self.data.imag)
else: return zeros(self.shape)
else:
try:
return self.__dict__[name]
except KeyError:
raise AttributeError, name
def __coerce__(self, other): return self, other
def checkreal(self):
if self.typecode=='D':
print "norms of this complex value:",
print "%s+j%s"%(norm(self.real),norm(self.imag))
else:
return 1
#------------------------------------------------------------------
# Indices and sclices
def __getitem__(self, i):
#print "getitem", i
if type(i) is type(1):
if self.shape[0] == 1: return self.data[0,i]
elif self.shape[1] == 1: return self.data[i,0]
else: return self.data[i]
elif type(i) is type((1,)):
result = self.data[i]
if type(i[0]) == type(1) and type(i[1]) == type(slice(1)):
result = result[_Axis,:]
elif type(i[1]) == type(1) and type(i[0]) == type(slice(1)):
result = result[:,_Axis]
else:
raise TypeError, i, type(i)
return to_Matrix(result)
def __setitem__(self, i, item):
#print "setitem"
item = to_array(item)
if type(i) is type(1):
if self.shape[0] == 1: self.data[0,i] = item
elif self.shape[1] == 1: self.data[i,0] = item
else: self.data[i] = item
def __getslice__(self, i, j):
#print "getslice", i, j, self.shape
if self.shape[0] == 1: result = self.data[0,i:j]
elif self.shape[1] == 1: result = self.data[i:j,0][:, _Axis]
else: result = self.data[i:j]
return to_Matrix(result)
def __setslice__(self, i, j, other):
#print "setslice"
if isinstance(other, Matrix): tmp = other.data
else: tmp = other
if self.shape[0] == 1: self.data[0,i:j] = tmp
elif self.shape[1] == 1: self.data[i:j,0] = tmp
else: self.data[i:j] = tmp
def __len__(self): return len(self.data)
#------------------------------------------------------------------
# Elementwise binary operations
def e_add (self, other): return Matrix(self.data + to_array(other))
def re_add(self, other): return Matrix(to_array(other) + self.data)
def e_sub (self, other): return Matrix(self.data - to_array(other))
def re_sub(self, other): return Matrix(to_array(other) - self.data)
def e_mul (self, other): return Matrix(self.data * to_array(other))
def re_mul(self, other): return Matrix(to_array(other) * self.data)
def e_div (self, other): return Matrix(self.data / to_array(other))
def re_div(self, other): return Matrix(to_array(other) / self.data)
def e_pow (self, n): return Matrix(self.data ** n)
def re_pow(self, n): return Matrix(n ** self.data)
# Operator names
__tadd__ = e_add
__rtadd__ = re_add
__tsub__ = e_sub
__rtsub__ = re_sub
__tmul__ = e_mul
__rtmul__ = re_mul
__tdiv__ = e_div
__rtdiv__ = re_div
__tpow__ = e_pow
__rtpow__ = re_pow
def __neg__(self): return Matrix(-self.data)
def __int__(self): return Matrix(self.data.astype(Int))
#------------------------------------------------------------------
# Elementwise comparison
def e_cmp(self, other):
"e_cmp(self, other): elementwise comparison (other may be scalar)"
if not isinstance(other, Matrix):
other = ones(size(self))*other
return mmap2(cmp, self, other)
def gt(self, other): return Matrix(_gt(self.data, to_array(other)))
def ge(self, other): return Matrix(_ge(self.data, to_array(other)))
def lt(self, other): return Matrix(_lt(self.data, to_array(other)))
def le(self, other): return Matrix(_le(self.data, to_array(other)))
def eq(self, other): return Matrix(_eq(self.data, to_array(other)))
#------------------------------------------------------------------
# Matrixwise binary operations
def checkshape(self, other, op):
"Check dimension is both are Matrix"
if isinstance(other, Matrix) and self.shape != other.shape:
raise ValueError, "%s shapes %s and %s" %(
op, self.shape, other.shape)
def m_add(self, other):
self.checkshape(other, "add")
return Matrix(self.data + to_array(other))
def rm_add(self, other):
#self.checkshape(other, "radd")
return Matrix(to_array(other) + self.data)
def m_sub(self, other):
self.checkshape(other, "sub")
return Matrix(self.data - to_array(other))
def rm_sub(self, other):
#self.checkshape(other, "rsub")
return Matrix(to_array(other) - self.data)
def m_mul(self, other):
#if isinstance(other, MatrixType): other = other.data
return to_Matrix(_mul(self.data, to_array(other)))
def rm_mul(self, other):
"Somehow this gets transposed by MLab"
#if isinstance(other, MatrixType): other = other.data
return to_Matrix(_mul(to_array(other), self.data)).T
def m_div(self, other):
if isinstance(other, MatrixType):
return solve(other.T, self.T).T
else:
return to_Matrix(self.data / other)
def rm_div(self, other):
if isinstance(other, MatrixType):
return solve(self.T, other.T).T
else:
return other * self.I
def m_sol(self, other):
return solve(self, other)
def rm_sol(self, other):
return solve(other, self)
def m_pow(self, n):
return mfuncC(lambda x,n=n:_pow(x,n), self)
def rm_pow(self, n):
return mfuncC(lambda x,n=n:_pow(n,x), self)
#------------------------------------------------------------------
# Matrixwise augmented assignments
def m_iadd(self, other):
self.checkshape(other, "iadd")
self.data = to_array(self.data + other)
self.shape = self.data.shape
self.typecode = self.data.typecode()
return self
def m_isub(self, other):
self.checkshape(other, "isub")
self.data = to_array(self.data - other)
self.shape = self.data.shape
self.typecode = self.data.typecode()
return self
def m_imul(self, other):
if isinstance(other, MatrixType): other = other.data
self.data = _mul(self.data, other)
self.shape = self.data.shape
self.typecode = self.data.typecode()
return self
def m_idiv(self, other):
"Temporary hack"
if isinstance(other, MatrixType):
self.data = _trsp(_sol(_trsp(other.data), _trsp(self.data)))
else:
self.data = self.data / other
self.shape = self.data.shape
self.typecode = self.data.typecode()
return self
def m_ipow(self, n):
self.data = _pow(self.data, n)
self.shape = self.data.shape
self.typecode = self.data.typecode()
return self
# Operator names
__add__ = m_add
__radd__ = rm_add
__sub__ = m_sub
__rsub__ = rm_sub
__mul__ = m_mul
__rmul__ = rm_mul
__div__ = m_div
__rdiv__ = rm_div
__sol__ = m_sol
__rsol__ = rm_sol
__pow__ = m_pow
__rpow__ = rm_pow
__iadd__ = m_iadd
__isub__ = m_isub
__imul__ = m_imul
__idiv__ = m_idiv
__ipow__ = m_ipow
def __cmp__(self, other): return norm(self-other) > 0
#------------------------------------------------------------------
# Casting
def to_list(data):
"to_list(data) : convert to double list"
data = to_array(data)
return map(list, data)
def to_Matrix(data):
"to_Matrix(data) : covert to Matrix if possible"
if isinstance(data, Matrix): return Matrix(data.data)
elif type(data) == type(_array([])): return Matrix(data)
elif type(data) == type([]): return Matrix(data)
else: return data
def to_array(data):
"to_array(data) : convert to array"
if isinstance(data, Matrix): return data.data
else: return _array(data)
def to_arrayC(data):
"to_arrayC(data) : convert to complex array"
if isinstance(data, Matrix): return data.data+0j
else: return _array(data)+0j
def to_number(data):
"to_number(data) : conbert to number if possible"
if hasattr(data, 'shape'):
if data.shape == (1,1): return data[0,0]
elif data.shape == (1): return data[0]
else: raise ValueError, "Can't convert to number"
else: return data
def Matrix_c(x):
"Matrix_c(x) : return a matrix from a list as a column"
if type(x) not in [type([]), type(())]:
raise TypeError, "argument must be a sequence"
lists = map(list, map(to_array,x))
import operator
return Matrix(reduce (operator.add, lists))
def Matrix_r(x):
"Matrix_r(x) : return a matrix from a list as a row"
if type(x) not in [type([]), type(())]:
raise TypeError, "argument must be a sequence"
return Matrix_c(map(lambda a:a.T, x)).T
def Matrix_cr(x):
"Matrix_cr(x) : return a matrix from a double list as a column of rows"
return Matrix_c(map(Matrix_r, x))
def rows(x):
"rows(x) : split x into a list of rows"
#print x
return map(to_Matrix, tuple(x.data))
def cols(x):
"cols(x) : split x into a list of columns"
return map(lambda x:x.T, rows(x.T))
def size(data):
"size(data) : return data.shape"
return to_Matrix(data).shape
#------------------------------------------------------------------
# General elementwise operations
# These are used when there is no NumPy functions.
def map1(f=None, *args):
"map1(f, *args): a costly reimplementation of map"
return apply(map, [f] + list(args))
def map2(f=None, *args):
"map2(f, *args): elementwise f on matrices args"
matlist = map(to_list, args)
return apply(map, [map] + [[f]*len(args[0])] + matlist)
def mmap2(*args):
"mmap2(f, *args): return Matrix as elementwise f on matrices args"
return Matrix(apply(map2, args))
def reduce1(f, vec):
"reduce1(f, list): a costly reimplementation of reduce"
return reduce(f, vec)
def reduce2(f, mat):
"reduce2(f, *args): elementwise f on matrices args"
return reduce(f, map(lambda x,f=f:reduce(f,x), to_list(mat)))
def mreduce2(f, mat):
"mreduce2(f, *args): return Matrix as elementwise f on matrices args"
return reduce2(f, map(list, args))
def all(mat):
"all(x) : return row 0-1 vector, indicates each column being all true"
#return reduce2(operator.__and__, mat)
return efunc(MLab.alltrue, mat)
def any(mat):
"any(x) : return row 0-1 vector, indicates each column being partly true"
#return reduce2(operator.__or__, mat)
return efunc(MLab.sometrue, mat)
def find(mat):
"find(x) : return list of indices (i,j) for which x[i,j] is true"
m,n = mat.shape
result = []
for i in xrange(m):
for j in xrange(n):
if mat[i,j]: result.append((i,j))
return result
#------------------------------------------------------------------
# Linear algebra
# These used in efuncs.py and mfuncs.py to general functions.
def plainsolve(a, b):
"solve(a,b) : return solution of a*x==b"
aT = a.T
return to_Matrix(_sol(a.data, b.data))
def solve(a, b):
"solve(a,b) : return LMS solution of a*x==b"
aT = a.T
return to_Matrix(_sol((aT*a).data, (aT*b).data))
sol = solve
def eig(x):
"""eig(x) : return eigenvalues as an array, eigenvectors as a Matrix
Definition: x * u.T = u.T * v
"""
s = x.shape
assert s[0] == s[1]
v, u = MLab.eig(x.data)
return v, to_Matrix(u)
def svd(a):
"svd(x) : return Matrix u, array s, Matrix v as svd decomposition"
u, s, v = MLab.svd(a.data)
return Matrix(u), s, Matrix(v)
def efunc(f, *args):
"efunc(f, *args) : elementwise function f of args"
return to_Matrix(apply(f, map(to_array, args)))
def efuncC(f, *args):
"""efuncC(f, *args) : elementwise function f of args with possible complex elements"""
y = to_Matrix(apply(f, map(to_arrayC, args)))
return approx_real(y)
def efuncRC(*args):
"""efuncRC(*args) : elementwise function of matrices args,
allow complex if one arg is complex"""
for x in args[1:]:
if x.typecode == "D":
return apply(efuncC, args)
else:
return apply(efunc, args)
def mfunc(f, x):
"""mfunc(f, x) : matrix function f of matrix x.
Note: Numeric defines (v,u) = eig(x) => x*u.T = u.T * Diag(v)
"""
(v, u) = eig(x)
uT = u.T
V = Matrix(_diag(f(v)))
y = uT * V * uT.I
return approx_real(y)
def mfuncC(f, x):
"""mfuncC(f, x) : matrix function with possibly complex eigenvalues.
Note: Numeric defines (v,u) = eig(x) => x*u.T = u.T * Diag(v)
"""
(v, u) = eig(x)
uT = u.T
V = Matrix(_diag(f(v+0j)))
y = uT * V * uT.I
return approx_real(y)
def mfunc_p(f, x):
"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 _min(v) + _eps < 0:
raise TypeError, "%s only applicable to positive numbers" % `f`
else:
uT = u.T
V = diag(Matrix(f(__builtin__.max(v,0))))
y = u.T * V * uT.I
return approx_real(y)
def approx_real(x):
"approx_real(x) : return x.real if |x.imag| < |x.real| * _eps_approx"
try:
xtype = x.typecode
except AttributeError:
if abs(x.imag) <= abs(x.real) * _eps_approx:
return x.real
else:
return x
if xtype == 'D' and norm(x.imag) <= norm(x.real) * _eps_approx:
return x.real
else:
return x
#------------------------------------------------------------------
# New matrices
def Mrange(*i):
"Mrange(*i) : row vector of range(i) - deprecated"
return to_Matrix(apply(range,i))
def r_range(*i):
"r_range(*i) : row vector of range(i)"
return to_Matrix(apply(range,i))
def c_range(*i):
"c_range(*i) : col vector of range(i)"
return to_Matrix(apply(range,i)).T
def zeros(x, typecode='d'):
"zeros(size, *typecode) : zero matrix of given size and type"
return to_Matrix(_zeros(x, typecode))
def ones(x):
"ones(size) : matrix of ones of given size"
return to_Matrix(_ones(x))
def eye(x):
"eye(size) : identity matrix of given size"
return to_Matrix(_eye(x))
def rand(*s):
"rand(size) : standard uniformly distributed random matrix of given size"
return to_Matrix(apply(_rand, s))
def unit(n, m):
"unit(n, m) : unit vector of dimension n, pointing at mth direction"
x = zeros(n).T
x[m] = 1.
return x
# Matrices of the same shape
def sort(x):
"sort(x) : return columnwise sorted"
return efunc(MLab.msort, x)
def diff(x):
"diff(x) : return columnwise diff"
return efunc(MLab.diff, x.T).T
def diff0(x):
"diff(x) : return columnwise diff, keeping first row. Inverse of cumsum"
return Matrix_c((a[0,:], diff(a)))
def trapz(x):
"trapz(x) : return columnwise integration using trapzoidal rule"
return efunc(MLab.trapz, x)
def cumsum(x):
"cumsum(x) : return columnwise cumulative sum"
return efunc(MLab.cumsum, x)
def cumprod(x):
"cumprod(x) : return columnwise cumulative prod"
return efunc(MLab.cumprod, x)
def flipud(x):
"flipud(x) : return up-down flipped"
return efunc(MLab.flipud, x)
def fliplr(x):
"flipud(x) : return left-right flipped"
return efunc(MLab.fliplr, x)
# Reduce to row vector
def norm1(x):
"norm1(x) : return columnwise norms"
return _sqrt(_sum(to_array(x)**2))
def sum(x):
"sum(x) : return columnwise sum"
return efunc(_sum, x)
def prod(x):
"prod(x) : return columnwise product"
return efunc(MLab.prod, x)
def mean(x):
"mean(x) : return columnwise mean"
return efunc(MLab.mean, x)
def median(x):
"median(x) : return columnwise median"
return efunc(MLab.median, x)
def std(x):
"std(x) : return columnwise standard deviation"
return efunc(MLab.std, x)
def ptp(x):
"ptp(x) : return columnwise ptp?"
return efunc(MLab.ptp, x)
def min(*args):
"""min(x) : return columnwise minimum if one argument,
return elementwise min across multiple arguments"""
if len(args) == 1:
return efunc(_min, args[0])
else:
return to_Matrix(reduce(_minimum, map(to_array,args)))
def max(*args):
"""max(x) : return columnwise maximum if one argument,
return elementwise max across multiple arguements"""
if len(args) == 1:
return efunc(_max, args[0])
else:
return to_Matrix(reduce(_maximum, map(to_array,args)))
# Different objects
def sum2(x):
"sum2(x) : return sum of all elements"
return _sum(_sum(to_array(x)))
def min2(x):
"min2(x) : return min of all elements"
return _min(_min(to_array(x)))
def max2(x):
"max2(x) : return max of all elements"
return _max(_max(to_array(x)))
def cov(x):
"std(x) : return covariance matrix for columns of x"
return efunc(MLab.cov, x)
def norm(x):
"norm(x) : return Frobenius norm of x"
return _sqrt(sum2(_abs(to_array(x))**2))
def mnorm(x):
"norm(x) : return Hilbert-Schmidt (?) norm of x"
return _sqrt(sum2(_abs(to_array(x.H * x))))
def Diag(x):
"Diag(x) : return diag matrix if x is row or col vector"
if x.shape[1] == 1: x = x.T
a = to_array(x)
if a.shape[0] == 1:
return Matrix(_diag(a[0]))
else:
raise ValueError, "x must be a vector"
def diag(x):
"diag(x) : return diagonal of x as row vector"
if (x.shape[1] == 1 or x.shape[0] == 1) and x.shape[1]*x.shape[0] != 1:
raise ValueError, "x cannot be a vector" + `x.shape`
a = to_array(x)
return Matrix(_diag(a))
|