__all__ = ['matrix', 'bmat', 'mat', 'asmatrix']
import sys
import numeric as N
from numeric import concatenate,isscalar,binary_repr,identity
from numpy.lib.utils import issubdtype
# make translation table
_table = [None]*256
for k in range(256):
_table[k] = chr(k)
_table = ''.join(_table)
_numchars = '0123456789.-+jeEL'
_todelete = []
for k in _table:
if k not in _numchars:
_todelete.append(k)
_todelete = ''.join(_todelete)
del k
def _eval(astr):
return eval(astr.translate(_table,_todelete))
def _convert_from_string(data):
rows = data.split(';')
newdata = []
count = 0
for row in rows:
trow = row.split(',')
newrow = []
for col in trow:
temp = col.split()
newrow.extend(map(_eval,temp))
if count == 0:
Ncols = len(newrow)
elif len(newrow) != Ncols:
raise ValueError, "Rows not the same size."
count += 1
newdata.append(newrow)
return newdata
def asmatrix(data, dtype=None):
""" Returns 'data' as a matrix. Unlike matrix(), no copy is performed
if 'data' is already a matrix or array. Equivalent to:
matrix(data, copy=False)
"""
return matrix(data, dtype=dtype, copy=False)
def matrix_power(M,n):
"""Raise a square matrix to the (integer) power n.
For positive integers n, the power is computed by repeated matrix
squarings and matrix multiplications. If n=0, the identity matrix
of the same type as M is returned. If n<0, the inverse is computed
and raised to the exponent.
Parameters
----------
M : array-like
Must be a square array (that is, of dimension two and with
equal sizes).
n : integer
The exponent can be any integer or long integer, positive
negative or zero.
Returns
-------
M to the power n
The return value is a an array the same shape and size as M;
if the exponent was positive or zero then the type of the
elements is the same as those of M. If the exponent was negative
the elements are floating-point.
Raises
------
LinAlgException
If the matrix is not numerically invertible, an exception is raised.
See Also
--------
The matrix() class provides an equivalent function as the exponentiation
operator.
Examples
--------
>>> matrix_power(array([[0,1],[-1,0]]),10)
array([[-1, 0],
[ 0, -1]])
"""
if len(M.shape) != 2 or M.shape[0] != M.shape[1]:
raise ValueError("input must be a square array")
if not issubdtype(type(n),int):
raise TypeError("exponent must be an integer")
from numpy.linalg import inv
if n==0:
M = M.copy()
M[:] = identity(M.shape[0])
return M
elif n<0:
M = inv(M)
n *= -1
result = M
if n <= 3:
for _ in range(n-1):
result=N.dot(result,M)
return result
# binary decomposition to reduce the number of Matrix
# multiplications for n > 3.
beta = binary_repr(n)
Z,q,t = M,0,len(beta)
while beta[t-q-1] == '0':
Z = N.dot(Z,Z)
q += 1
result = Z
for k in range(q+1,t):
Z = N.dot(Z,Z)
if beta[t-k-1] == '1':
result = N.dot(result,Z)
return result
class matrix(N.ndarray):
"""
mat = matrix(data, dtype=None, copy=True)
Returns a matrix from an array-like object, or a string of
data. A matrix is a specialized 2-d array that retains
its 2-d nature through operations and where '*' means matrix
multiplication and '**' means matrix power.
Parameters
----------
data : array-like or string
If data is a string, then interpret the string as a matrix
with commas or spaces separating columns and semicolons
separating rows.
If data is array-like than convert the array to a matrix.
dtype : data-type
Anything that can be interpreted as a NumPy datatype.
copy : bool
If data is already an ndarray, then this flag determines whether
or not the data will be copied
Examples
--------
>>> import numpy as np
>>> a = np.matrix('1 2; 3 4')
>>> print a
[[1 2]
[3 4]]
"""
__array_priority__ = 10.0
def __new__(subtype, data, dtype=None, copy=True):
if isinstance(data, matrix):
dtype2 = data.dtype
if (dtype is None):
dtype = dtype2
if (dtype2 == dtype) and (not copy):
return data
return data.astype(dtype)
if isinstance(data, N.ndarray):
if dtype is None:
intype = data.dtype
else:
intype = N.dtype(dtype)
new = data.view(subtype)
if intype != data.dtype:
return new.astype(intype)
if copy: return new.copy()
else: return new
if isinstance(data, str):
data = _convert_from_string(data)
# now convert data to an array
arr = N.array(data, dtype=dtype, copy=copy)
ndim = arr.ndim
shape = arr.shape
if (ndim > 2):
raise ValueError, "matrix must be 2-dimensional"
elif ndim == 0:
shape = (1,1)
elif ndim == 1:
shape = (1,shape[0])
order = False
if (ndim == 2) and arr.flags.fortran:
order = True
if not (order or arr.flags.contiguous):
arr = arr.copy()
ret = N.ndarray.__new__(subtype, shape, arr.dtype,
buffer=arr,
order=order)
return ret
def __array_finalize__(self, obj):
self._getitem = False
if (isinstance(obj, matrix) and obj._getitem): return
ndim = self.ndim
if (ndim == 2):
return
if (ndim > 2):
newshape = tuple([x for x in self.shape if x > 1])
ndim = len(newshape)
if ndim == 2:
self.shape = newshape
return
elif (ndim > 2):
raise ValueError, "shape too large to be a matrix."
else:
newshape = self.shape
if ndim == 0:
self.shape = (1,1)
elif ndim == 1:
self.shape = (1,newshape[0])
return
def __getitem__(self, index):
self._getitem = True
try:
out = N.ndarray.__getitem__(self, index)
finally:
self._getitem = False
if not isinstance(out, N.ndarray):
return out
if out.ndim == 0:
return out[()]
if out.ndim == 1:
sh = out.shape[0]
# Determine when we should have a column array
try:
n = len(index)
except:
n = 0
if n > 1 and isscalar(index[1]):
out.shape = (sh,1)
else:
out.shape = (1,sh)
return out
def _get_truendim(self):
shp = self.shape
truend = 0
for val in shp:
if (val > 1): truend += 1
return truend
def __mul__(self, other):
if isinstance(other,(N.ndarray, list, tuple)) :
# This promotes 1-D vectors to row vectors
return N.dot(self, asmatrix(other))
if isscalar(other) or not hasattr(other, '__rmul__') :
return N.dot(self, other)
return NotImplemented
def __rmul__(self, other):
return N.dot(other, self)
def __imul__(self, other):
self[:] = self * other
return self
def __pow__(self, other):
return matrix_power(self, other)
def __rpow__(self, other):
return NotImplemented
def __repr__(self):
s = repr(self.__array__()).replace('array', 'matrix')
# now, 'matrix' has 6 letters, and 'array' 5, so the columns don't
# line up anymore. We need to add a space.
l = s.splitlines()
for i in range(1, len(l)):
if l[i]:
l[i] = ' ' + l[i]
return '\n'.join(l)
def __str__(self):
return str(self.__array__())
def _align(self, axis):
"""A convenience function for operations that need to preserve axis
orientation.
"""
if axis is None:
return self[0,0]
elif axis==0:
return self
elif axis==1:
return self.transpose()
else:
raise ValueError, "unsupported axis"
# Necessary because base-class tolist expects dimension
# reduction by x[0]
def tolist(self):
return self.__array__().tolist()
# To preserve orientation of result...
def sum(self, axis=None, dtype=None, out=None):
"""Sum the matrix over the given axis. If the axis is None, sum
over all dimensions. This preserves the orientation of the
result as a row or column.
"""
return N.ndarray.sum(self, axis, dtype, out)._align(axis)
def mean(self, axis=None, dtype=None, out=None):
"""Compute the mean along the specified axis.
Returns the average of the array elements. The average is taken over
the flattened array by default, otherwise over the specified axis.
Parameters
----------
axis : integer
Axis along which the means are computed. The default is
to compute the standard deviation of the flattened array.
dtype : type
Type to use in computing the means. For arrays of integer type
the default is float32, for arrays of float types it is the
same as the array type.
out : ndarray
Alternative output array in which to place the result. It must
have the same shape as the expected output but the type will be
cast if necessary.
Returns
-------
mean : The return type varies, see above.
A new array holding the result is returned unless out is
specified, in which case a reference to out is returned.
SeeAlso
-------
var : variance
std : standard deviation
Notes
-----
The mean is the sum of the elements along the axis divided by the
number of elements.
"""
return N.ndarray.mean(self, axis, dtype, out)._align(axis)
def std(self, axis=None, dtype=None, out=None, ddof=0):
"""Compute the standard deviation along the specified axis.
Returns the standard deviation of the array elements, a measure of the
spread of a distribution. The standard deviation is computed for the
flattened array by default, otherwise over the specified axis.
Parameters
----------
axis : integer
Axis along which the standard deviation is computed. The
default is to compute the standard deviation of the flattened
array.
dtype : type
Type to use in computing the standard deviation. For arrays of
integer type the default is float32, for arrays of float types
it is the same as the array type.
out : ndarray
Alternative output array in which to place the result. It must
have the same shape as the expected output but the type will be
cast if necessary.
ddof : {0, integer}
Means Delta Degrees of Freedom. The divisor used in calculations
is N-ddof.
Returns
-------
standard deviation : The return type varies, see above.
A new array holding the result is returned unless out is
specified, in which case a reference to out is returned.
SeeAlso
-------
var : variance
mean : average
Notes
-----
The standard deviation is the square root of the
average of the squared deviations from the mean, i.e. var =
sqrt(mean(abs(x - x.mean())**2)). The computed standard
deviation is computed by dividing by the number of elements,
N-ddof. The option ddof defaults to zero, that is, a biased
estimate. Note that for complex numbers std takes the absolute
value before squaring, so that the result is always real
and nonnegative.
"""
return N.ndarray.std(self, axis, dtype, out, ddof)._align(axis)
def var(self, axis=None, dtype=None, out=None, ddof=0):
"""Compute the variance along the specified axis.
Returns the variance of the array elements, a measure of the spread of
a distribution. The variance is computed for the flattened array by
default, otherwise over the specified axis.
Parameters
----------
axis : integer
Axis along which the variance is computed. The default is to
compute the variance of the flattened array.
dtype : data-type
Type to use in computing the variance. For arrays of integer
type the default is float32, for arrays of float types it is
the same as the array type.
out : ndarray
Alternative output array in which to place the result. It must
have the same shape as the expected output but the type will be
cast if necessary.
ddof : {0, integer}
Means Delta Degrees of Freedom. The divisor used in calculations
is N-ddof.
Returns
-------
variance : depends, see above
A new array holding the result is returned unless out is
specified, in which case a reference to out is returned.
SeeAlso
-------
std : standard deviation
mean : average
Notes
-----
The variance is the average of the squared deviations from the
mean, i.e. var = mean(abs(x - x.mean())**2). The mean is
computed by dividing by N-ddof, where N is the number of elements.
The argument ddof defaults to zero; for an unbiased estimate
supply ddof=1. Note that for complex numbers the absolute value
is taken before squaring, so that the result is always real
and nonnegative.
"""
return N.ndarray.var(self, axis, dtype, out)._align(axis)
def prod(self, axis=None, dtype=None, out=None):
return N.ndarray.prod(self, axis, dtype, out)._align(axis)
def any(self, axis=None, out=None):
return N.ndarray.any(self, axis, out)._align(axis)
def all(self, axis=None, out=None):
return N.ndarray.all(self, axis, out)._align(axis)
def max(self, axis=None, out=None):
return N.ndarray.max(self, axis, out)._align(axis)
def argmax(self, axis=None, out=None):
return N.ndarray.argmax(self, axis, out)._align(axis)
def min(self, axis=None, out=None):
return N.ndarray.min(self, axis, out)._align(axis)
def argmin(self, axis=None, out=None):
return N.ndarray.argmin(self, axis, out)._align(axis)
def ptp(self, axis=None, out=None):
return N.ndarray.ptp(self, axis, out)._align(axis)
def getI(self):
M,N = self.shape
if M == N:
from numpy.dual import inv
else:
from numpy.dual import pinv
return asmatrix(func(self))
def getA(self):
return self.__array__()
def getA1(self):
return self.__array__().ravel()
def getT(self):
return self.transpose()
def getH(self):
if issubclass(self.dtype.type, N.complexfloating):
return self.transpose().conjugate()
else:
return self.transpose()
T = property(getT, None, doc="transpose")
A = property(getA, None, doc="base array")
A1 = property(getA1, None, doc="1-d base array")
H = property(getH, None, doc="hermitian (conjugate) transpose")
I = property(getI, None, doc="inverse")
def _from_string(str,gdict,ldict):
rows = str.split(';')
rowtup = []
for row in rows:
trow = row.split(',')
newrow = []
for x in trow:
newrow.extend(x.split())
trow = newrow
coltup = []
for col in trow:
col = col.strip()
try:
thismat = ldict[col]
except KeyError:
try:
thismat = gdict[col]
except KeyError:
raise KeyError, "%s not found" % (col,)
coltup.append(thismat)
rowtup.append(concatenate(coltup,axis=-1))
return concatenate(rowtup,axis=0)
def bmat(obj, ldict=None, gdict=None):
"""
Build a matrix object from string, nested sequence, or array.
Examples
--------
>>> F = bmat('A, B; C, D')
>>> F = bmat([[A,B],[C,D]])
>>> F = bmat(r_[c_[A,B],c_[C,D]])
All of these produce the same matrix::
[ A B ]
[ C D ]
if A, B, C, and D are appropriately shaped 2-d arrays.
"""
if isinstance(obj, str):
if gdict is None:
# get previous frame
frame = sys._getframe().f_back
glob_dict = frame.f_globals
loc_dict = frame.f_locals
else:
glob_dict = gdict
loc_dict = ldict
return matrix(_from_string(obj, glob_dict, loc_dict))
if isinstance(obj, (tuple, list)):
# [[A,B],[C,D]]
arr_rows = []
for row in obj:
if isinstance(row, N.ndarray): # not 2-d
return matrix(concatenate(obj,axis=-1))
else:
arr_rows.append(concatenate(row,axis=-1))
return matrix(concatenate(arr_rows,axis=0))
if isinstance(obj, N.ndarray):
return matrix(obj)
mat = asmatrix
|