diff options
Diffstat (limited to 'numpy/matrixlib/defmatrix.py')
| -rw-r--r-- | numpy/matrixlib/defmatrix.py | 677 |
1 files changed, 677 insertions, 0 deletions
diff --git a/numpy/matrixlib/defmatrix.py b/numpy/matrixlib/defmatrix.py new file mode 100644 index 000000000..dfa7e4a8d --- /dev/null +++ b/numpy/matrixlib/defmatrix.py @@ -0,0 +1,677 @@ +__all__ = ['matrix', 'bmat', 'mat', 'asmatrix'] + +import sys +import numpy.core.numeric as N +from numpy.core.numeric import concatenate, isscalar, binary_repr, identity, asanyarray +from numpy.core.numerictypes 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): + """ + Interpret the input as a matrix. + + Unlike `matrix`, `asmatrix` does not make a copy if the input is already + a matrix or an ndarray. Equivalent to ``matrix(data, copy=False)``. + + Parameters + ---------- + data : array_like + Input data. + + Returns + ------- + mat : matrix + `data` interpreted as a matrix. + + Examples + -------- + >>> x = np.array([[1, 2], [3, 4]]) + + >>> m = np.asmatrix(x) + + >>> x[0,0] = 5 + + >>> m + matrix([[5, 2], + [3, 4]]) + + """ + 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 + -------- + >>> np.linalg.matrix_power(np.array([[0,1],[-1,0]]),10) + array([[-1, 0], + [ 0, -1]]) + + """ + M = asanyarray(M) + 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): + """ + matrix(data, dtype=None, copy=True) + + Returns a matrix from an array-like object, or from a string + of data. A matrix is a specialized 2-d array that retains + its 2-d nature through operations. It has certain special + operators, such as ``*`` (matrix multiplication) and + ``**`` (matrix power). + + Parameters + ---------- + data : array_like or string + If data is a string, the string is interpreted as a matrix + with commas or spaces separating columns, and semicolons + separating rows. + dtype : data-type + Data-type of the output matrix. + copy : bool + If data is already an ndarray, then this flag determines whether + the data is copied, or whether a view is constructed. + + See Also + -------- + array + + Examples + -------- + >>> a = np.matrix('1 2; 3 4') + >>> print a + [[1 2] + [3 4]] + + >>> np.matrix([[1, 2], [3, 4]]) + matrix([[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 __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 __ipow__(self, other): + self[:] = self ** other + return self + + 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, ddof)._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): + """ + Test whether any array element along a given axis evaluates to True. + + Refer to `numpy.any` for full documentation. + + Parameters + ---------- + axis: int, optional + Axis along which logical OR is performed + out: ndarray, optional + Output to existing array instead of creating new one, must have + same shape as expected output + + Returns + ------- + any : bool, ndarray + Returns a single bool if `axis` is ``None``; otherwise, + returns `ndarray` + + """ + 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 as func + else: + from numpy.dual import pinv as func + return asmatrix(func(self)) + + def getA(self): + return self.__array__() + + def getA1(self): + return self.__array__().ravel() + + def getT(self): + """ + m.T + + Returns the transpose of m. + + Examples + -------- + >>> m = np.matrix('[1, 2; 3, 4]') + >>> m + matrix([[1, 2], + [3, 4]]) + >>> m.T + matrix([[1, 3], + [2, 4]]) + + See Also + -------- + transpose + + """ + 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 a string, nested sequence, or array. + + Parameters + ---------- + obj : string, sequence or array + Input data. Variables names in the current scope may + be referenced, even if `obj` is a string. + + Returns + ------- + out : matrix + Returns a matrix object, which is a specialized 2-D array. + + See Also + -------- + matrix + + Examples + -------- + >>> A = np.mat('1 1; 1 1') + >>> B = np.mat('2 2; 2 2') + >>> C = np.mat('3 4; 5 6') + >>> D = np.mat('7 8; 9 0') + + All the following expressions construct the same block matrix: + + >>> np.bmat([[A, B], [C, D]]) + matrix([[1, 1, 2, 2], + [1, 1, 2, 2], + [3, 4, 7, 8], + [5, 6, 9, 0]]) + >>> np.bmat(np.r_[np.c_[A, B], np.c_[C, D]]) + matrix([[1, 1, 2, 2], + [1, 1, 2, 2], + [3, 4, 7, 8], + [5, 6, 9, 0]]) + >>> np.bmat('A,B; C,D') + matrix([[1, 1, 2, 2], + [1, 1, 2, 2], + [3, 4, 7, 8], + [5, 6, 9, 0]]) + + """ + 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 |