Python Matrix.rref方法代码示例

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import rref [as 别名]
def test_nullspace():
    # first test reduced row-ech form
    R = Rational

    M = Matrix([[5,7,2,1],
               [1,6,2,-1]])
    out, tmp = M.rref()
    assert out == Matrix([[1,0,-R(2)/23,R(13)/23],
                              [0,1,R(8)/23, R(-6)/23]])

    M = Matrix([[-5,-1, 4,-3,-1],
                [ 1,-1,-1, 1, 0],
                [-1, 0, 0, 0, 0],
                [ 4, 1,-4, 3, 1],
                [-2, 0, 2,-2,-1]])
    assert M*M.nullspace()[0] == Matrix(5,1,[0]*5)

    M = Matrix([[1,3,0,2,6,3,1],
                [-2,-6,0,-2,-8,3,1],
                [3,9,0,0,6,6,2],
                [-1,-3,0,1,0,9,3]])
    out, tmp = M.rref()
    assert out == Matrix([[1,3,0,0,2,0,0],
                               [0,0,0,1,2,0,0],
                               [0,0,0,0,0,1,R(1)/3],
                               [0,0,0,0,0,0,0]])

    # now check the vectors
    basis = M.nullspace()
    assert basis[0] == Matrix([-3,1,0,0,0,0,0])
    assert basis[1] == Matrix([0,0,1,0,0,0,0])
    assert basis[2] == Matrix([-2,0,0,-2,1,0,0])
    assert basis[3] == Matrix([0,0,0,0,0,R(-1)/3, 1])

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import rref [as 别名]
def linear_rref(A, b, Matrix=None, S=None):
    """ Transform a linear system to reduced row-echelon form

    Transforms both the matrix and right-hand side of a linear
    system of equations to reduced row echelon form

    Parameters
    ----------
    A: Matrix-like
        iterable of rows
    b: iterable

    Returns
    -------
    A', b' - transformed versions

    """
    if Matrix is None:
        from sympy import Matrix
    if S is None:
        from sympy import S
    mat_rows = [_map2l(S, list(row) + [v]) for row, v in zip(A, b)]
    aug = Matrix(mat_rows)
    raug, pivot = aug.rref()
    nindep = len(pivot)
    return raug[:nindep, :-1], raug[:nindep, -1]

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import rref [as 别名]
 def balance(self):
     # Get list of unique elements
     elements = set()
     [elements.add(element)
      for reactant in self.left for element in reactant.count().keys()]
     elements = tuple(elements)
     # Build the matrix
     rows = []
     for element in elements:
         row = []
         for reactant in self.left:
             row.append(reactant.count(element))
         for reactant in self.right:
             row.append(-reactant.count(element))
         rows.append(row)
     # Balance equation with linear algebra
     # http://www.ctroms.com/blog/math/2011/05/01/balancing-chemical-equations-with-linear-algebra/
     mat = Matrix(rows)
     solution, pivots = mat.rref()
     values = [solution.row(r)[-1] for r in range(solution.rows)]
     factor = lcm([value.as_numer_denom()[1] for value in values])
     coeffs = [-solution.row(i)[i] * solution.row(i)[-1]
               * factor for i in pivots] + [factor]
     for reactant, coeff in zip(self.left + self.right, coeffs):
         reactant.coefficient = coeff
     return self

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import rref [as 别名]
def find_dependent_rational_basis(x):
    """
    Find linearly dependent columns similar to MATLAB null(x, 'r').
    :param x: The input numpy array
    :return: The linear combination matrix (Null Space of Matrix).
    """
    m, n = x.shape
    x = Matrix(x)
    R = x.rref()
    r = len(R[1])
    nopiv = np.arange(0, n)
    p = np.array(R[1])
    a = np.array(R[0])
    nopiv = np.delete(nopiv, p, axis = 0)
    Z = np.zeros([n, n-r])
    if n > r:
        Z[nopiv, :] = np.eye(n-r, n-r)
        if r > 0:
            Z[p, :] = -a[0:r, nopiv]
    return(Z)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import rref [as 别名]
from sympy import symbols, Matrix
from numpy import exp

A = Matrix([[2, -1, 0, 2], [-1, 2, -1, 4], [0, -1, 2, 6]])
print A.rref()

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import rref [as 别名]

#.........这里部分代码省略.........

    # test_inverse
    A = eye(4)
    assert A.inv() == eye(4)
    assert A.inv("LU") == eye(4)
    assert A.inv("ADJ") == eye(4)
    A = SMatrix([[2,3,5],
                [3,6,2],
                [8,3,6]])
    Ainv = A.inv()
    assert A*Ainv == eye(3)
    assert A.inv("LU") == Ainv
    assert A.inv("ADJ") == Ainv

    # test_cross
    v1 = Matrix(1,3,[1,2,3])
    v2 = Matrix(1,3,[3,4,5])
    assert v1.cross(v2) == Matrix(1,3,[-2,4,-2])
    assert v1.norm(v1) == 14

    # test_cofactor
    assert eye(3) == eye(3).cofactorMatrix()
    test = SMatrix([[1,3,2],[2,6,3],[2,3,6]])
    assert test.cofactorMatrix() == SMatrix([[27,-6,-6],[-12,2,3],[-3,1,0]])
    test = SMatrix([[1,2,3],[4,5,6],[7,8,9]])
    assert test.cofactorMatrix() == SMatrix([[-3,6,-3],[6,-12,6],[-3,6,-3]])

    # test_jacobian
    x = Symbol('x')
    y = Symbol('y')
    L = SMatrix(1,2,[x**2*y, 2*y**2 + x*y])
    syms = [x,y]
    assert L.jacobian(syms) == Matrix([[2*x*y, x**2],[y, 4*y+x]])

    L = SMatrix(1,2,[x, x**2*y**3])
    assert L.jacobian(syms) == SMatrix([[1, 0], [2*x*y**3, x**2*3*y**2]])

    # test_QR
    A = Matrix([[1,2],[2,3]])
    Q, S = A.QRdecomposition()
    R = Rational
    assert Q == Matrix([[5**R(-1,2), (R(2)/5)*(R(1)/5)**R(-1,2)], [2*5**R(-1,2), (-R(1)/5)*(R(1)/5)**R(-1,2)]])
    assert S == Matrix([[5**R(1,2), 8*5**R(-1,2)], [0, (R(1)/5)**R(1,2)]])
    assert Q*S == A
    assert Q.T * Q == eye(2)

    # test nullspace
    # first test reduced row-ech form
    R = Rational

    M = Matrix([[5,7,2,1],
               [1,6,2,-1]])
    out, tmp = M.rref()
    assert out == Matrix([[1,0,-R(2)/23,R(13)/23],
                              [0,1,R(8)/23, R(-6)/23]])

    M = Matrix([[1,3,0,2,6,3,1],
                [-2,-6,0,-2,-8,3,1],
                [3,9,0,0,6,6,2],
                [-1,-3,0,1,0,9,3]])
    out, tmp = M.rref()
    assert out == Matrix([[1,3,0,0,2,0,0],
                               [0,0,0,1,2,0,0],
                               [0,0,0,0,0,1,R(1)/3],
                               [0,0,0,0,0,0,0]])
    # now check the vectors
    basis = M.nullspace()
    assert basis[0] == Matrix([[-3,1,0,0,0,0,0]])
    assert basis[1] == Matrix([[0,0,1,0,0,0,0]])
    assert basis[2] == Matrix([[-2,0,0,-2,1,0,0]])
    assert basis[3] == Matrix([[0,0,0,0,0,R(-1)/3, 1]])


    # test eigen
    x = Symbol('x')
    y = Symbol('y')
    eye3 = eye(3)
    assert eye3.charpoly(x) == (1-x)**3
    assert eye3.charpoly(y) == (1-y)**3
    # test values
    M = Matrix([(0,1,-1),
                (1,1,0),
                (-1,0,1) ])
    vals = M.eigenvals()
    vals.sort()
    assert vals == [-1, 1, 2]

    R = Rational
    M = Matrix([ [1,0,0],
                 [0,1,0],
                 [0,0,1]])
    assert M.eigenvects() == [[1, 3, [Matrix(1,3,[1,0,0]), Matrix(1,3,[0,1,0]), Matrix(1,3,[0,0,1])]]]
    M = Matrix([ [5,0,2],
                 [3,2,0],
                 [0,0,1]])
    assert M.eigenvects() == [[1, 1, [Matrix(1,3,[R(-1)/2,R(3)/2,1])]],
                              [2, 1, [Matrix(1,3,[0,1,0])]],
                              [5, 1, [Matrix(1,3,[1,1,0])]]]

    assert M.zeros((3, 5)) == SMatrix(3, 5, {})

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import rref [as 别名]
def find_basis_set(vecs, **kwargs):
    A = Matrix(vecs)
    basiscols = A.rref()[1]

    return PY.double(vecs[:, PY.array(basiscols)])

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import rref [as 别名]

#.........这里部分代码省略.........
                                raise_keyerror_if_weight_absent = (
                                    states[new_level][tuple(new_weight)])
                                states[new_level][tuple(new_weight)][
                                    "states"].append(new_state)
                            except KeyError:
                                states[new_level][tuple(new_weight)] = {
                                    "states" : [new_state],
                                    "rotation_to_ob" : Matrix([[1]]),
                                }
                        except KeyError:
                            states[new_level] = {
                                tuple(new_weight) : {
                                    "states" : [new_state],
                                    "rotation_to_ob" : Matrix([[1]]),
                                }
                            }
            #resolve degeneracy
            new_level = level + 1
            for weight in weights[new_level]:
                states_for_this_weight = states[new_level][weight]["states"]
                degeneracy = len(states_for_this_weight)
                if degeneracy == 1:
                    continue
                scalar_product_matrix = Matrix(degeneracy, degeneracy,
                    lambda i,j : states_for_this_weight[i]["norm"]*
                    states_for_this_weight[j]["norm"]*self.scalar_product(
                        states_for_this_weight[i]["lowering_chain"],
                        states_for_this_weight[j]["lowering_chain"],
                        simple_root_length_squared_list, cartan_matrix,
                        weights)
                    if i > j else 0)
                scalar_product_matrix += scalar_product_matrix.T
                scalar_product_matrix += eye(degeneracy)
                rref = scalar_product_matrix.rref()
                dependents = [index for index in range(degeneracy)
                              if index not in rref[1]]
                # calculate additional matrix elements (if any)
                for independent in rref[1]:
                    for state_num in range(degeneracy):
                        if state_num != independent:
                            # state_num's parent might be linked to independent
                            state_num_norm = states_for_this_weight[state_num][
                                "norm"]
                            state_num_matrix_element_information = (
                                states_for_this_weight[state_num][
                                    "matrix_element_information"])
                            parent_level = (
                                state_num_matrix_element_information[0][
                                    "level"])
                            parent_weight = (
                                state_num_matrix_element_information[0][
                                    "weight"])
                            parent_state_num = (
                                state_num_matrix_element_information[0][
                                    "state_num"])
                            direction = (
                                state_num_matrix_element_information[0][
                                    "direction"])
                            parent_norm = states[parent_level][parent_weight][
                                "states"][parent_state_num]["norm"]
                            matrix_element = (
                                Rational(1,1)*scalar_product_matrix[
                                    state_num, independent]*parent_norm/
                                state_num_norm)
                            if matrix_element == 0:
                                continue