Python sympy.Matrix() Examples

def get_equations(self):
        """
        :return: Functions to calculate A, B and f given state x and input u
        """
        f = sp.zeros(3, 1)

        x = sp.Matrix(sp.symbols('x y theta', real=True))
        u = sp.Matrix(sp.symbols('v w', real=True))

        f[0, 0] = u[0, 0] * sp.cos(x[2, 0])
        f[1, 0] = u[0, 0] * sp.sin(x[2, 0])
        f[2, 0] = u[1, 0]

        f = sp.simplify(f)
        A = sp.simplify(f.jacobian(x))
        B = sp.simplify(f.jacobian(u))

        f_func = sp.lambdify((x, u), f, 'numpy')
        A_func = sp.lambdify((x, u), A, 'numpy')
        B_func = sp.lambdify((x, u), B, 'numpy')

        return f_func, A_func, B_func 

def _generateTransitionMatrix(self, A=None):#, transitionExpressionList=None):
        '''
        Finds the transition matrix from the set of ode.  It is
        important to note that although some of the functions used
        in this method appear to be the same as _getReactantMatrix
        and _getStateChangeMatrix, they are different in the sense
        that the functions called here is focused on the terms of
        the equation rather than the states.
        '''
        if A is None:
            if not ode_utils.none_or_empty_list(self._odeList):
                eqn_list = [t.equation for t in self._odeList]
                A = sympy.Matrix(checkEquation(eqn_list,
                                               *self._getListOfVariablesDict(),
                                               subs_derived=False))
            else:
                raise Exception("Object was not initialized using a set of ode")

        bdList, _term = _ode_composition.getUnmatchedExpressionVector(A, True)
        fx = _ode_composition.stripBDFromOde(A, bdList)
        states = [s for s in self._iterStateList()]
        M, _remain = _ode_composition.odeToPureTransition(fx, states, True)
        return M 

def prepare_channel_operator_list(ops_list):
        """
        Prepares a list of channel operators.

        Args:
            ops_list (List): The list of operators to prepare

        Returns:
            List: The channel operator list
        """
        # convert to sympy matrices and verify that each singleton is
        # in a tuple; also add identity matrix
        from sympy import Matrix, eye
        result = []
        for ops in ops_list:
            if not isinstance(ops, tuple) and not isinstance(ops, list):
                ops = [ops]
            result.append([Matrix(op) for op in ops])
        n = result[0][0].shape[0]  # grab the dimensions from the first element
        result = [[eye(n)]] + result
        return result

    # pylint: disable=invalid-name 

def as_matrix(self):
        """Returns the data of the table in Matrix form.

        Examples
        ========
        >>> from sympy import TableForm
        >>> t = TableForm([[5, 7], [4, 2], [10, 3]], headings='automatic')
        >>> t
          | 1  2
        --------
        1 | 5  7
        2 | 4  2
        3 | 10 3
        >>> t.as_matrix()
        Matrix([
        [ 5, 7],
        [ 4, 2],
        [10, 3]])
        """
        from sympy import Matrix
        return Matrix(self._lines) 

def diff_inference():
    
    r,s=symbols('r,s')
    
    la1,la2,lb1,lb2=symbols('l^a_1,l^a_2,l^b_1,l^b_2')
    
    la1=(1-s)/2
    la2=(1+s)/2
    lb1=(1-r)/2
    lb2=(1+r)/2
    N1=la1*lb1
    N2=la1*lb2
    N3=la2*lb1
    N4=la2*lb2
    
    N=Matrix([[N1,0,N2,0,N3,0,N4,0],
              [0,N1,0,N2,0,N3,0,N4]]) 

def test_conv10b():
    A = sympy.Matrix([[sympy.Symbol("x"), sympy.Symbol("y")],
                     [sympy.Symbol("z"), sympy.Symbol("t")]])
    assert sympify(A) == DenseMatrix(2, 2, [Symbol("x"), Symbol("y"),
                                            Symbol("z"), Symbol("t")])

    B = sympy.Matrix([[1, 2], [3, 4]])
    assert sympify(B) == DenseMatrix(2, 2, [Integer(1), Integer(2), Integer(3),
                                            Integer(4)])

    C = sympy.Matrix([[7, sympy.Symbol("y")],
                     [sympy.Function("g")(sympy.Symbol("z")), 3 + 2*sympy.I]])
    assert sympify(C) == DenseMatrix(2, 2, [Integer(7), Symbol("y"),
                                            function_symbol("g",
                                                            Symbol("z")),
                                            3 + 2*I]) 

def block_matrix(blocks):
    """
    Construct a matrix where the elements are specified by the block structure
    by joining the blocks appropriately.

    Parameters
    ----------
    blocks : two level deep iterable of sympy Matrix objects
        The block specification of the matrices used to construct the block
        matrix.

    Returns
    -------
    matrix : sympy Matrix
        A matrix whose elements are the elements of the blocks with the
        specified block structure.
    """
    return sp.Matrix.col_join(
        *tuple(
            sp.Matrix.row_join(
                *tuple(mat for mat in row)) for row in blocks
        )
    ) 

def generate_channel_quadratic_programming_matrices(
            self, channel, symbols):
        """
        Generate matrices for quadratic program.

        Args:
             channel (Matrix): a 4x4 symbolic matrix
             symbols (list): the symbols x1, ..., xn which may occur in the matrix

        Returns:
            list: A list of 4x4 complex matrices ([D1, D2, ..., Dn], E) such that:
            channel == x1*D1 + ... + xn*Dn + E
        """
        return (
            [self.get_matrix_from_channel(channel, symbol) for symbol in symbols],
            self.get_const_matrix_from_channel(channel, symbols)
        ) 

def main():
    Format()
    a = Matrix ( 2, 2, ( 1, 2, 3, 4 ) )
    b = Matrix ( 2, 1, ( 5, 6 ) )
    c = a * b
    print(a,b,'=',c)

    x, y = symbols ( 'x, y' )

    d = Matrix ( 1, 2, ( x ** 3, y ** 3 ))
    e = Matrix ( 2, 2, ( x ** 2, 2 * x * y, 2 * x * y, y ** 2 ) )
    f = d * e

    print('%',d,e,'=',f)

    # xpdf()
    xpdf(pdfprog=None)
    return 

def main():
    Format()
    a = Matrix ( 2, 2, ( 1, 2, 3, 4 ) )
    b = Matrix ( 2, 1, ( 5, 6 ) )
    c = a * b
    print(a,b,'=',c)

    x, y = symbols ('x, y')

    d = Matrix ( 1, 2, ( x ** 3, y ** 3 ))
    e = Matrix ( 2, 2, ( x ** 2, 2 * x * y, 2 * x * y, y ** 2 ) )
    f = d * e

    print('%',d,e,'=',f)

    # xpdf()
    xpdf(pdfprog=None)
    return 

def test_conv10():
    A = DenseMatrix(1, 4, [Integer(1), Integer(2), Integer(3), Integer(4)])
    assert (A._sympy_() == sympy.Matrix(1, 4,
                                        [sympy.Integer(1), sympy.Integer(2),
                                         sympy.Integer(3), sympy.Integer(4)]))

    B = DenseMatrix(4, 1, [Symbol("x"), Symbol("y"), Symbol("z"), Symbol("t")])
    assert (B._sympy_() == sympy.Matrix(4, 1,
                                        [sympy.Symbol("x"), sympy.Symbol("y"),
                                         sympy.Symbol("z"), sympy.Symbol("t")])
            )

    C = DenseMatrix(2, 2,
                    [Integer(5), Symbol("x"),
                     function_symbol("f", Symbol("x")), 1 + I])

    assert (C._sympy_() ==
            sympy.Matrix([[5, sympy.Symbol("x")],
                          [sympy.Function("f")(sympy.Symbol("x")),
                           1 + sympy.I]])) 

def fcts(self):

        j = sympy.Matrix([
            r**2 * sympy.sin(t) * z,
            r * sympy.sin(t) * z,
            r * sympy.sin(t) * z**2,
        ])

        # Create an isotropic sigma vector
        Sig = sympy.Matrix([
            [1120/(69*sympy.pi)*(r*z)**2 * sympy.sin(t)**2, 0, 0],
            [0, 1120/(69*sympy.pi)*(r*z)**2 * sympy.sin(t)**2, 0],
            [0, 0, 1120/(69*sympy.pi)*(r*z)**2 * sympy.sin(t)**2]
        ])

        return j, Sig 

def _qsympify_sequence(seq):
    """Convert elements of a sequence to standard form.

    This is like sympify, but it performs special logic for arguments passed
    to QExpr. The following conversions are done:

    * (list, tuple, Tuple) => _qsympify_sequence each element and convert
      sequence to a Tuple.
    * basestring => Symbol
    * Matrix => Matrix
    * other => sympify

    Strings are passed to Symbol, not sympify to make sure that variables like
    'pi' are kept as Symbols, not the SymPy built-in number subclasses.

    Examples
    ========

    >>> from sympsi.qexpr import _qsympify_sequence
    >>> _qsympify_sequence((1,2,[3,4,[1,]]))
    (1, 2, (3, 4, (1,)))

    """

    return tuple(__qsympify_sequence_helper(seq)) 

def linearize_symbolic(self,
                            zeros=False) -> List[sympy.MutableDenseMatrix]:
        nx = len(self.x)
        nu = len(self.u)
        ny = len(self.y)
        nf = len(self.f)
        ng = len(self.g)
        A = sympy.Matrix([]) if zeros == False else sympy.Matrix.zeros(nx, nx)
        B = sympy.Matrix([]) if zeros == False else sympy.Matrix.zeros(nx, nu)
        C = sympy.Matrix([]) if zeros == False else sympy.Matrix.zeros(ny, nx)
        D = sympy.Matrix([]) if zeros == False else sympy.Matrix.zeros(ny, nu)
        if nx > 0:
            if nf > 0:
                A = self.f.jacobian(self.x)
            if ng > 0:
                C = self.g.jacobian(self.x)
        if nu > 0:
            if nf > 0:
                B = self.f.jacobian(self.u)
            if ng > 0:
                D = self.g.jacobian(self.u)
        return [A, B, C, D] 

def _findIndex(eq_vec, expr):
    """
    Given a vector of expressions, find where you will locate the
    input term.

    Parameters
    ----------
    eq_vec: :class:`sympy.Matrix`
        vector of sympy equation
    expr: sympy type
        An expression that we would like to find

    Returns
    -------
    list:
        of index that contains the expression.  Can be an empty list
        or with multiple integer
    """
    out = list()
    for i, a in enumerate(eq_vec):
        j = _hasExpression(a, expr)
        if j is True:
            out.append(i)
    return out 

def pureTransitionToOde(A):
    """
    Get the ode from a pure transition matrix

    Parameters
    ----------
    A: `sympy.Matrix`
        a transition matrix of size [n \times n]

    Returns
    -------
    b: `sympy.Matrix`
        a matrix of size [n \times 1] which is the ode
    """
    nrow, ncol = A.shape
    assert nrow == ncol, "Need a square matrix"
    B = [sum(A[:, i]) - sum(A[i, :]) for i in range(nrow)]
    return sympy.simplify(sympy.Matrix(B)) 

def symbolic_rotation_matrix(phi, theta, symbolic_psi):
    """Retourne une matrice de rotation où psi est symbolique"""
    return sympy.Matrix(rz_matrix(phi)) * sympy.Matrix(rx_matrix(theta)) * symbolic_rz_matrix(symbolic_psi) 

def _interpolation_coeffs(self):
        """
        Symbolic expression for the coefficients for sparse point interpolation
        according to:

            https://en.wikipedia.org/wiki/Bilinear_interpolation.

        Returns
        -------
        Matrix of coefficient expressions.
        """
        # Grid indices corresponding to the corners of the cell ie x1, y1, z1
        indices1 = tuple(sympy.symbols('%s1' % d) for d in self.grid.dimensions)
        indices2 = tuple(sympy.symbols('%s2' % d) for d in self.grid.dimensions)
        # 1, x1, y1, z1, x1*y1, ...
        indices = list(powerset(indices1))
        indices[0] = (1,)
        point_sym = list(powerset(self.sfunction._point_symbols))
        point_sym[0] = (1,)
        # 1, px. py, pz, px*py, ...
        A = []
        ref_A = [np.prod(ind) for ind in indices]
        # Create the matrix with the same increment order as the point increment
        for i in self.sfunction._point_increments:
            # substitute x1 by x2 if increment in that dimension
            subs = dict((indices1[d], indices2[d] if i[d] == 1 else indices1[d])
                        for d in range(len(i)))
            A += [[1] + [a.subs(subs) for a in ref_A[1:]]]

        A = sympy.Matrix(A)
        # Coordinate values of the sparse point
        p = sympy.Matrix([[np.prod(ind)] for ind in point_sym])

        # reference cell x1:0, x2:h_x
        left = dict((a, 0) for a in indices1)
        right = dict((b, dim.spacing) for b, dim in zip(indices2, self.grid.dimensions))
        reference_cell = {**left, **right}
        # Substitute in interpolation matrix
        A = A.subs(reference_cell)
        return A.inv().T * p 

def quat_matrix_r(p):
  return sp.Matrix([[p[0], -p[1], -p[2], -p[3]],
                    [p[1],  p[0],  p[3], -p[2]],
                    [p[2], -p[3],  p[0],  p[1]],
                    [p[3],  p[2], -p[1],  p[0]]]) 

def _sympy_tridiag(a, b):
    """Creates the tridiagonal sympy matrix tridiag(b, a, b).
    """
    n = len(a)
    assert n == len(b)
    A = [[0 for _ in range(n)] for _ in range(n)]
    for i in range(n):
        A[i][i] = a[i]
    for i in range(n - 1):
        A[i][i + 1] = b[i + 1]
        A[i + 1][i] = b[i + 1]
    return sympy.Matrix(A) 

def central_diff_weights(points, divisions=1):
    # Check the cache
    if (divisions, points) in central_diff_weights_precomputed:
        return central_diff_weights_precomputed[(divisions, points)]

    if points < divisions + 1:
        raise ValueError("Points < divisions + 1, cannot compute")
    if points % 2 == 0:
        raise ValueError("Odd number of points required")
    ho = points >> 1 
    
    x = [[xi] for xi in range(-ho, ho+1)]
    X = []
    for xi in x:
        line = [1.0] + [xi[0]**k for k in range(1, points)]
        X.append(line)
    factor = product(range(1, divisions + 1))
#    from scipy.linalg import inv
    from sympy import Matrix
    # The above coefficients were generated from a routine which used Fractions
    # Additional coefficients cannot reliable be computed with numpy or floating
    # point numbers - the error is too great
    # sympy must be used for reliability
    inverted = [[float(j) for j in i] for i in Matrix(X).inv().tolist()]
    w = [i*factor for i in inverted[divisions]]
#    w = [i*factor for i in (inv(X)[divisions]).tolist()]
    central_diff_weights_precomputed[(divisions, points)] = w
    return w 

def _interpret_measurement(self, measurements):
        # reverse measurement bitstrings and remove all zero entry
        linear = [(k[::-1], v) for k, v in measurements.items()
                  if k != "0" * len(self._oracle.variable_register)]
        # sort bitstrings by their probabilities
        linear.sort(key=lambda x: x[1], reverse=True)

        # construct matrix
        equations = []
        for k, _ in linear:
            equations.append([int(c) for c in k])
        y = Matrix(equations)

        # perform gaussian elimination
        y_transformed = y.rref(iszerofunc=lambda x: x % 2 == 0)

        def mod(x, modulus):
            numer, denom = x.as_numer_denom()
            return numer * mod_inverse(denom, modulus) % modulus
        y_new = y_transformed[0].applyfunc(lambda x: mod(x, 2))

        # determine hidden string from final matrix
        rows, _ = y_new.shape
        hidden = [0] * len(self._oracle.variable_register)
        for r in range(rows):
            yi = [i for i, v in enumerate(list(y_new[r, :])) if v == 1]
            if len(yi) == 2:
                hidden[yi[0]] = '1'
                hidden[yi[1]] = '1'
        return "".join(str(x) for x in hidden)[::-1] 

def symbolic_rz_matrix(symbolic_theta):
    """Matrice symbolique de rotation autour de l'axe Z"""
    return sympy.Matrix([
        [sympy.cos(symbolic_theta), -sympy.sin(symbolic_theta), 0],
        [sympy.sin(symbolic_theta), sympy.cos(symbolic_theta), 0],
        [0, 0, 1]
    ]) 

def symbolic_axis_rotation_matrix(axis, symbolic_theta):
    """Returns a rotation matrix around the given axis"""
    [x, y, z] = axis
    c = sympy.cos(symbolic_theta)
    s = sympy.sin(symbolic_theta)
    return sympy.Matrix(_axis_rotation_matrix_formula(x, y, z, c, s)) 

def _apply_geometric_transformations(self, theta, symbolic):

        if symbolic:
            # Angle symbolique qui paramètre la rotation du joint en cours
            symbolic_frame_matrix = np.eye(4)

            # Apply translation matrix
            symbolic_frame_matrix = symbolic_frame_matrix * sympy.Matrix(geometry.homogeneous_translation_matrix(*self.translation_vector))

            # Apply orientation matrix
            symbolic_frame_matrix = symbolic_frame_matrix * geometry.cartesian_to_homogeneous(geometry.rpy_matrix(*self.orientation))

            # Apply rotation matrix
            if self.rotation is not None:
                symbolic_frame_matrix = symbolic_frame_matrix * geometry.cartesian_to_homogeneous(geometry.symbolic_axis_rotation_matrix(self.rotation, theta), matrix_type="sympy")

            symbolic_frame_matrix = sympy.lambdify(theta, symbolic_frame_matrix, "numpy")

            return symbolic_frame_matrix

        else:
            # Init the transformation matrix
            frame_matrix = np.eye(4)

            # First, apply translation matrix
            frame_matrix = np.dot(frame_matrix, geometry.homogeneous_translation_matrix(*self.translation_vector))

            # Apply orientation
            frame_matrix = np.dot(frame_matrix, geometry.cartesian_to_homogeneous(geometry.rpy_matrix(*self.orientation)))

            # Apply rotation matrix
            if self.rotation is not None:
                frame_matrix = np.dot(frame_matrix, geometry.cartesian_to_homogeneous(geometry.axis_rotation_matrix(self.rotation, theta)))

            return frame_matrix 

def get_generating_function(f, T=None):
    """
    Get the generating function

    :param f: function addr
    :param T: edges matrix of the function f
    :return: the generating function
    """

    if not T:
        edges = [(a[0].addr, a[1].addr) for a in f.graph.edges()]
        N = sorted([x for x in f.block_addrs])
        T = []
        for b1 in N:
            T.append([])
            for b2 in N:
                if b1 == b2 or (b1, b2) in edges:
                    T[-1].append(1)
                else:
                    T[-1].append(0)
    else:
        N = T[0]

    T = sympy.Matrix(T)
    z = sympy.var('z')
    I = sympy.eye(len(N))

    tmp = I - z * T
    tmp.row_del(len(N) - 1)
    tmp.col_del(0)
    det_num = tmp.det()
    det_den = (I - z * T).det()
    quot = det_num / det_den
    g_z = ((-1) ** (len(N) + 1)) * quot
    return g_z 

def __init__(self):
        self.t = sympy.symbols('t')
        self.x = sympy.Matrix([])
        self.u = sympy.Matrix([])
        self.y = sympy.Matrix([])
        self.p = sympy.Matrix([])
        self.c = sympy.Matrix([])
        self.v = sympy.Matrix([])
        self.x0 = {}
        self.u0 = {}
        self.p0 = {}
        self.c0 = {}
        self.eqs = []
        self.f = None
        self.g = None 

def quat_matrix_l(p):
  return sp.Matrix([[p[0], -p[1], -p[2], -p[3]],
                    [p[1],  p[0], -p[3],  p[2]],
                    [p[2],  p[3],  p[0], -p[1]],
                    [p[3], -p[2],  p[1],  p[0]]]) 

def quat_rotate(q0, q1, q2, q3):
  # make symbolic rotation matrix from quat
  return sp.Matrix([[q0**2 + q1**2 - q2**2 - q3**2, 2*(q1*q2 + q0*q3), 2*(q1*q3 - q0*q2)],
                    [2*(q1*q2 - q0*q3), q0**2 - q1**2 + q2**2 - q3**2, 2*(q2*q3 + q0*q1)],
                    [2*(q1*q3 + q0*q2), 2*(q2*q3 - q0*q1), q0**2 - q1**2 - q2**2 + q3**2]]).T 

def __init__(self):

        super(Aircraft, self).__init__()

        # states
        self.x = sympy.Matrix([])
        self.x0 = {
            }

        # variables
        self.v = sympy.Matrix([])

        # constants
        self.c = sympy.Matrix([])
        self.c0 = {
            }

        # parameters
        self.p = sympy.Matrix([])
        self.p0 = {
            }

        # inputs
        self.u = sympy.Matrix([])
        self.u0 = {
            }

        # outputs
        self.y = sympy.Matrix([])

        # equations
        self.eqs = [
            ,
            ]

        self.compute_fg()