Source code for sympy.matrices.expressions.matexpr

from __future__ import print_function, division

from functools import wraps

from sympy.core import S, Symbol, Tuple, Integer, Basic, Expr, Eq
from sympy.core.decorators import call_highest_priority
from sympy.core.compatibility import range
from sympy.core.sympify import SympifyError, sympify
from sympy.functions import conjugate, adjoint
from sympy.functions.special.tensor_functions import KroneckerDelta
from sympy.matrices import ShapeError
from sympy.simplify import simplify
from sympy.utilities.misc import filldedent


def _sympifyit(arg, retval=None):
    # This version of _sympifyit sympifies MutableMatrix objects
    def deco(func):
        @wraps(func)
        def __sympifyit_wrapper(a, b):
            try:
                b = sympify(b, strict=True)
                return func(a, b)
            except SympifyError:
                return retval

        return __sympifyit_wrapper

    return deco


[docs]class MatrixExpr(Basic): """ Superclass for Matrix Expressions MatrixExprs represent abstract matrices, linear transformations represented within a particular basis. Examples ======== >>> from sympy import MatrixSymbol >>> A = MatrixSymbol('A', 3, 3) >>> y = MatrixSymbol('y', 3, 1) >>> x = (A.T*A).I * A * y See Also ======== MatrixSymbol MatAdd MatMul Transpose Inverse """ # Should not be considered iterable by the # sympy.core.compatibility.iterable function. Subclass that actually are # iterable (i.e., explicit matrices) should set this to True. _iterable = False _op_priority = 11.0 is_Matrix = True is_MatrixExpr = True is_Identity = None is_Inverse = False is_Transpose = False is_ZeroMatrix = False is_MatAdd = False is_MatMul = False is_commutative = False def __new__(cls, *args, **kwargs): args = map(sympify, args) return Basic.__new__(cls, *args, **kwargs) # The following is adapted from the core Expr object def __neg__(self): return MatMul(S.NegativeOne, self).doit() def __abs__(self): raise NotImplementedError @_sympifyit('other', NotImplemented) @call_highest_priority('__radd__') def __add__(self, other): return MatAdd(self, other).doit() @_sympifyit('other', NotImplemented) @call_highest_priority('__add__') def __radd__(self, other): return MatAdd(other, self).doit() @_sympifyit('other', NotImplemented) @call_highest_priority('__rsub__') def __sub__(self, other): return MatAdd(self, -other).doit() @_sympifyit('other', NotImplemented) @call_highest_priority('__sub__') def __rsub__(self, other): return MatAdd(other, -self).doit() @_sympifyit('other', NotImplemented) @call_highest_priority('__rmul__') def __mul__(self, other): return MatMul(self, other).doit() @_sympifyit('other', NotImplemented) @call_highest_priority('__rmul__') def __matmul__(self, other): return MatMul(self, other).doit() @_sympifyit('other', NotImplemented) @call_highest_priority('__mul__') def __rmul__(self, other): return MatMul(other, self).doit() @_sympifyit('other', NotImplemented) @call_highest_priority('__mul__') def __rmatmul__(self, other): return MatMul(other, self).doit() @_sympifyit('other', NotImplemented) @call_highest_priority('__rpow__') def __pow__(self, other): if not self.is_square: raise ShapeError("Power of non-square matrix %s" % self) elif self.is_Identity: return self elif other is S.NegativeOne: return Inverse(self) elif other is S.Zero: return Identity(self.rows) elif other is S.One: return self return MatPow(self, other) @_sympifyit('other', NotImplemented) @call_highest_priority('__pow__') def __rpow__(self, other): raise NotImplementedError("Matrix Power not defined") @_sympifyit('other', NotImplemented) @call_highest_priority('__rdiv__') def __div__(self, other): return self * other**S.NegativeOne @_sympifyit('other', NotImplemented) @call_highest_priority('__div__') def __rdiv__(self, other): raise NotImplementedError() #return MatMul(other, Pow(self, S.NegativeOne)) __truediv__ = __div__ __rtruediv__ = __rdiv__ @property def rows(self): return self.shape[0] @property def cols(self): return self.shape[1] @property def is_square(self): return self.rows == self.cols def _eval_conjugate(self): from sympy.matrices.expressions.adjoint import Adjoint from sympy.matrices.expressions.transpose import Transpose return Adjoint(Transpose(self)) def as_real_imag(self): from sympy import I real = (S(1)/2) * (self + self._eval_conjugate()) im = (self - self._eval_conjugate())/(2*I) return (real, im) def _eval_inverse(self): from sympy.matrices.expressions.inverse import Inverse return Inverse(self) def _eval_transpose(self): return Transpose(self) def _eval_power(self, exp): return MatPow(self, exp) def _eval_simplify(self, **kwargs): if self.is_Atom: return self else: return self.__class__(*[simplify(x, **kwargs) for x in self.args]) def _eval_adjoint(self): from sympy.matrices.expressions.adjoint import Adjoint return Adjoint(self) def _entry(self, i, j): raise NotImplementedError( "Indexing not implemented for %s" % self.__class__.__name__) def adjoint(self): return adjoint(self)
[docs] def as_coeff_Mul(self, rational=False): """Efficiently extract the coefficient of a product. """ return S.One, self
def conjugate(self): return conjugate(self) def transpose(self): from sympy.matrices.expressions.transpose import transpose return transpose(self) T = property(transpose, None, None, 'Matrix transposition.') def inverse(self): return self._eval_inverse() @property def I(self): return self.inverse() def valid_index(self, i, j): def is_valid(idx): return isinstance(idx, (int, Integer, Symbol, Expr)) return (is_valid(i) and is_valid(j) and (self.rows is None or (0 <= i) != False and (i < self.rows) != False) and (0 <= j) != False and (j < self.cols) != False) def __getitem__(self, key): if not isinstance(key, tuple) and isinstance(key, slice): from sympy.matrices.expressions.slice import MatrixSlice return MatrixSlice(self, key, (0, None, 1)) if isinstance(key, tuple) and len(key) == 2: i, j = key if isinstance(i, slice) or isinstance(j, slice): from sympy.matrices.expressions.slice import MatrixSlice return MatrixSlice(self, i, j) i, j = sympify(i), sympify(j) if self.valid_index(i, j) != False: return self._entry(i, j) else: raise IndexError("Invalid indices (%s, %s)" % (i, j)) elif isinstance(key, (int, Integer)): # row-wise decomposition of matrix rows, cols = self.shape # allow single indexing if number of columns is known if not isinstance(cols, Integer): raise IndexError(filldedent(''' Single indexing is only supported when the number of columns is known.''')) key = sympify(key) i = key // cols j = key % cols if self.valid_index(i, j) != False: return self._entry(i, j) else: raise IndexError("Invalid index %s" % key) elif isinstance(key, (Symbol, Expr)): raise IndexError(filldedent(''' Only integers may be used when addressing the matrix with a single index.''')) raise IndexError("Invalid index, wanted %s[i,j]" % self)
[docs] def as_explicit(self): """ Returns a dense Matrix with elements represented explicitly Returns an object of type ImmutableDenseMatrix. Examples ======== >>> from sympy import Identity >>> I = Identity(3) >>> I I >>> I.as_explicit() Matrix([ [1, 0, 0], [0, 1, 0], [0, 0, 1]]) See Also ======== as_mutable: returns mutable Matrix type """ from sympy.matrices.immutable import ImmutableDenseMatrix return ImmutableDenseMatrix([[ self[i, j] for j in range(self.cols)] for i in range(self.rows)])
[docs] def as_mutable(self): """ Returns a dense, mutable matrix with elements represented explicitly Examples ======== >>> from sympy import Identity >>> I = Identity(3) >>> I I >>> I.shape (3, 3) >>> I.as_mutable() Matrix([ [1, 0, 0], [0, 1, 0], [0, 0, 1]]) See Also ======== as_explicit: returns ImmutableDenseMatrix """ return self.as_explicit().as_mutable()
def __array__(self): from numpy import empty a = empty(self.shape, dtype=object) for i in range(self.rows): for j in range(self.cols): a[i, j] = self[i, j] return a
[docs] def equals(self, other): """ Test elementwise equality between matrices, potentially of different types >>> from sympy import Identity, eye >>> Identity(3).equals(eye(3)) True """ return self.as_explicit().equals(other)
def canonicalize(self): return self def as_coeff_mmul(self): return 1, MatMul(self)
class MatrixElement(Expr): parent = property(lambda self: self.args[0]) i = property(lambda self: self.args[1]) j = property(lambda self: self.args[2]) _diff_wrt = True is_symbol = True is_commutative = True def __new__(cls, name, n, m): n, m = map(sympify, (n, m)) from sympy import MatrixBase if isinstance(name, (MatrixBase,)): if n.is_Integer and m.is_Integer: return name[n, m] name = sympify(name) obj = Expr.__new__(cls, name, n, m) return obj def doit(self, **kwargs): deep = kwargs.get('deep', True) if deep: args = [arg.doit(**kwargs) for arg in self.args] else: args = self.args return args[0][args[1], args[2]] def _eval_derivative(self, v): if not isinstance(v, MatrixElement): from sympy import MatrixBase if isinstance(self.parent, MatrixBase): return self.parent.diff(v)[self.i, self.j] return S.Zero if self.args[0] != v.args[0]: return S.Zero return KroneckerDelta(self.args[1], v.args[1])*KroneckerDelta(self.args[2], v.args[2])
[docs]class MatrixSymbol(MatrixExpr): """Symbolic representation of a Matrix object Creates a SymPy Symbol to represent a Matrix. This matrix has a shape and can be included in Matrix Expressions >>> from sympy import MatrixSymbol, Identity >>> A = MatrixSymbol('A', 3, 4) # A 3 by 4 Matrix >>> B = MatrixSymbol('B', 4, 3) # A 4 by 3 Matrix >>> A.shape (3, 4) >>> 2*A*B + Identity(3) I + 2*A*B """ is_commutative = False def __new__(cls, name, n, m): n, m = sympify(n), sympify(m) obj = Basic.__new__(cls, name, n, m) return obj def _hashable_content(self): return(self.name, self.shape) @property def shape(self): return self.args[1:3] @property def name(self): return self.args[0] def _eval_subs(self, old, new): # only do substitutions in shape shape = Tuple(*self.shape)._subs(old, new) return MatrixSymbol(self.name, *shape) def __call__(self, *args): raise TypeError( "%s object is not callable" % self.__class__ ) def _entry(self, i, j): return MatrixElement(self, i, j) @property def free_symbols(self): return set((self,)) def doit(self, **hints): if hints.get('deep', True): return type(self)(self.name, self.args[1].doit(**hints), self.args[2].doit(**hints)) else: return self def _eval_simplify(self, **kwargs): return self
[docs]class Identity(MatrixExpr): """The Matrix Identity I - multiplicative identity >>> from sympy.matrices import Identity, MatrixSymbol >>> A = MatrixSymbol('A', 3, 5) >>> I = Identity(3) >>> I*A A """ is_Identity = True def __new__(cls, n): return super(Identity, cls).__new__(cls, sympify(n)) @property def rows(self): return self.args[0] @property def cols(self): return self.args[0] @property def shape(self): return (self.args[0], self.args[0]) def _eval_transpose(self): return self def _eval_trace(self): return self.rows def _eval_inverse(self): return self def conjugate(self): return self def _entry(self, i, j): eq = Eq(i, j) if eq is S.true: return S.One elif eq is S.false: return S.Zero return KroneckerDelta(i, j) def _eval_determinant(self): return S.One
[docs]class ZeroMatrix(MatrixExpr): """The Matrix Zero 0 - additive identity >>> from sympy import MatrixSymbol, ZeroMatrix >>> A = MatrixSymbol('A', 3, 5) >>> Z = ZeroMatrix(3, 5) >>> A+Z A >>> Z*A.T 0 """ is_ZeroMatrix = True def __new__(cls, m, n): return super(ZeroMatrix, cls).__new__(cls, m, n) @property def shape(self): return (self.args[0], self.args[1]) @_sympifyit('other', NotImplemented) @call_highest_priority('__rpow__') def __pow__(self, other): if other != 1 and not self.is_square: raise ShapeError("Power of non-square matrix %s" % self) if other == 0: return Identity(self.rows) if other < 1: raise ValueError("Matrix det == 0; not invertible.") return self def _eval_transpose(self): return ZeroMatrix(self.cols, self.rows) def _eval_trace(self): return S.Zero def _eval_determinant(self): return S.Zero def conjugate(self): return self def _entry(self, i, j): return S.Zero def __nonzero__(self): return False __bool__ = __nonzero__
def matrix_symbols(expr): return [sym for sym in expr.free_symbols if sym.is_Matrix] from .matmul import MatMul from .matadd import MatAdd from .matpow import MatPow from .transpose import Transpose from .inverse import Inverse