def test_inverse():
    n, m, l = symbols("n m l", integer=True)
    A = MatrixSymbol("A", n, m)
    C = MatrixSymbol("C", n, n)
    D = MatrixSymbol("D", n, n)
    E = MatrixSymbol("E", m, n)

    raises(ShapeError, lambda: Inverse(A))
    assert Inverse(Inverse(C)) == C

    assert Inverse(C) * C == Identity(C.rows)

    assert Inverse(eye(3)) == eye(3)

    assert Inverse(S(3)) == S(1) / 3

    assert Inverse(Identity(n)) == Identity(n)

    # Simplifies Muls if possible (i.e. submatrices are square)
    assert Inverse(C * D) == D.I * C.I
    # But still works when not possible
    assert Inverse(A * E).is_Inverse

    # We play nice with traditional explicit matrices
    assert Inverse(Matrix([[1, 2], [3, 4]])) == Matrix([[1, 2], [3, 4]]).inv()

    def _represent_ZGate(self, basis, **options):
        """
            Represents the (I)QFT In the Z Basis
        """
        nqubits = options.get('nqubits', 0)
        if nqubits == 0:
            raise QuantumError(
                'The number of qubits must be given as nqubits.')
        if nqubits < self.min_qubits:
            raise QuantumError(
                'The number of qubits %r is too small for the gate.' % nqubits
            )
        size = self.size
        omega = self.omega

        #Make a matrix that has the basic Fourier Transform Matrix
        arrayFT = [[omega**(
            i*j % size)/sqrt(size) for i in range(size)] for j in range(size)]
        matrixFT = Matrix(arrayFT)

        #Embed the FT Matrix in a higher space, if necessary
        if self.label[0] != 0:
            matrixFT = matrix_tensor_product(eye(2**self.label[0]), matrixFT)
        if self.min_qubits < nqubits:
            matrixFT = matrix_tensor_product(
                matrixFT, eye(2**(nqubits - self.min_qubits)))

        return matrixFT

def test_inverse():
    raises(ShapeError, lambda: Inverse(A))
    raises(ShapeError, lambda: Inverse(A*B))

    assert Inverse(C).args == (C, S(-1))
    assert Inverse(C).shape == (n, n)
    assert Inverse(A*E).shape == (n, n)
    assert Inverse(E*A).shape == (m, m)
    assert Inverse(C).inverse() == C
    assert isinstance(Inverse(Inverse(C)), Inverse)

    assert Inverse(*Inverse(E*A).args) == Inverse(E*A)

    assert C.inverse().inverse() == C

    assert C.inverse()*C == Identity(C.rows)

    assert Identity(n).inverse() == Identity(n)
    assert (3*Identity(n)).inverse() == Identity(n)/3

    # Simplifies Muls if possible (i.e. submatrices are square)
    assert (C*D).inverse() == D.I*C.I
    # But still works when not possible
    assert isinstance((A*E).inverse(), Inverse)
    assert Inverse(C*D).doit(inv_expand=False) == Inverse(C*D)

    assert Inverse(eye(3)).doit() == eye(3)
    assert Inverse(eye(3)).doit(deep=False) == eye(3)

def test(As,Bs,Cs,Ds,Us,sI_A,Js,Ys,B0,C0,C,Psi,Lamda):

    D1=BlockMatrix([[Us.inv()*Psi,zeros(Us.rows,1)],[-Lamda,eye(Lamda.rows)]]).as_mutable()
    D2=BlockMatrix([[sI_A,B0],[-C,Js]]).as_mutable()
    D3=BlockMatrix([[As,Bs],[-Cs,Ds]]).as_mutable()
    D4=BlockMatrix([[C0,-Ys],[zeros(Ys.cols,C0.cols),eye(Ys.cols)]]).as_mutable()
    return expand(simplify(D1*D2))==expand(D3*D4)

def test_inverse():
    n, m, l = symbols('n m l', integer=True)
    A = MatrixSymbol('A', n, m)
    B = MatrixSymbol('B', m, l)
    C = MatrixSymbol('C', n, n)
    D = MatrixSymbol('D', n, n)
    E = MatrixSymbol('E', m, n)

    raises(ShapeError, "Inverse(A)")
    assert Inverse(Inverse(C)) == C

    assert Inverse(C)*C == Identity(C.n)

    assert Inverse(eye(3)) == eye(3)

    assert Inverse(S(3)) == S(1)/3

    assert Inverse(Identity(n)) == Identity(n)

    # Simplifies Muls if possible (i.e. submatrices are square)
    assert Inverse(C*D) == D.I*C.I
    # But still works when not possible
    assert Inverse(A*E).is_Inverse

    # We play nice with traditional explicit matrices
    assert Inverse(Matrix([[1,2],[3,4]])) == Matrix([[1,2],[3,4]]).inv()

def test(B0,Ts,E,A,B,C,U,V,Rc,Rs,Ss,p,m,rt):
       
    D1=BlockDiagMatrix(B0,eye(p)).as_mutable()
    D2=BlockMatrix([[Ts,U],[-V,zeros(p,m)]]).as_mutable()
    D3=BlockMatrix([[s*E-A,B],[-C,zeros(C.rows,B.cols)]]).as_mutable()
    D4Block11=simplify((Ss*Rs.inv())).as_mutable()
    D4=BlockDiagMatrix(D4Block11,eye(m)).as_mutable()
    return simplify(D1*D2)==simplify(D3*D4) 

def test(rE, r, p, m, PE, PA, PC, As, Bs, Cs, Ds, E, A, B, C, D):
    D1Block2 = simplify(PC * s * (PE - s * PA).inv())
    D1Block3 = BlockMatrix([[Bs], [Ds]]).as_mutable()
    D1Block4 = BlockMatrix([[zeros(r, p)], [eye(p)]]).as_mutable()
    D1 = BlockMatrix([[eye(r + p), D1Block2, D1Block3, D1Block4]]).as_mutable()
    D2 = BlockMatrix(2, 2, [s * E - A, B, -C, D]).as_mutable()
    D3 = BlockMatrix(2, 2, [As, Bs, -Cs, Ds]).as_mutable()
    D4 = BlockMatrix(
        [[eye(r), zeros(r, rE + 2 * p), zeros(r, m)], [zeros(m, r), zeros(m, rE + 2 * p), eye(m)]]
    ).as_mutable()
    return simplify(D1 * D2) == simplify(D3 * D4)

def test(B0,Ts,E,E0,A0,A,B,C,U,V,Rc,Rs,Qs,Ss,p,m,rt,n):

    D1=BlockMatrix([[zeros(n,rt+p)],[eye(rt +p)]]).as_mutable()
    D2=BlockMatrix([[Ts,U],[-V,zeros(p,m)]]).as_mutable()
    D3=BlockMatrix([[s*E-A,B],[-C,zeros(C.rows,B.cols)]]).as_mutable()
    D4Block1=simplify((-Ss*(Qs.inv()))).as_mutable()
    D4Block1=D4Block1.col_join(eye(rt))
    D4=BlockDiagMatrix(D4Block1,eye(m)).as_mutable()
    

    return simplify(D1*D2) == simplify(D3*D4)

def test_det():
    assert isinstance(Determinant(A), Determinant)
    assert not isinstance(Determinant(A), MatrixExpr)
    raises(ShapeError, lambda: Determinant(C))
    print det(eye(3))
    assert det(eye(3)) == 1
    assert det(Matrix(3, 3, [1, 3, 2, 4, 1, 3, 2, 5, 2])) == 17
    A / det(A)  # Make sure this is possible

    raises(TypeError, lambda: Determinant(S.One))

    assert Determinant(A).arg is A

def make_collision(choice):
    model_switcher = {
        # Relax 2nd moments for hydro, SOI
        'hydro': (eye(q) - S_Relax) * temp_populations
                 + S_Relax * hardcoded_cm_eq
                 + (eye(q) - S_Relax / 2) * hardcoded_F_cm,
        # Relax 1st moments for ADE, SOI without force
        'ade': None  # TODO: write the collision for advection-diffusion equation
    }
    # Get the function from switcher dictionary
    cm_after_collision = model_switcher.get(choice, lambda: "Invalid argument")
    print_as_vector(cm_after_collision, print_symbol=pop_in_str)

def test_transpose():
    Sq = MatrixSymbol('Sq', n, n)

    assert Transpose(A).shape == (m, n)
    assert Transpose(A*B).shape == (l, n)
    assert Transpose(Transpose(A)) == A

    assert Transpose(eye(3)) == eye(3)

    assert Transpose(S(5)) == S(5)

    assert Transpose(Matrix([[1, 2], [3, 4]])) == Matrix([[1, 3], [2, 4]])

    assert Transpose(Trace(Sq)) == Trace(Sq)

def test_transpose():
    n, m, l = symbols('n m l', integer=True)
    A = MatrixSymbol('A', n, m)
    B = MatrixSymbol('B', m, l)

    assert Transpose(A).shape == (m,n)
    assert Transpose(A*B).shape == (l,n)
    assert Transpose(Transpose(A)) == A

    assert Transpose(eye(3)) == eye(3)

    assert Transpose(S(5)) == S(5)

    assert Transpose(Matrix([[1,2],[3,4]])) == Matrix([[1,3],[2,4]])

def test_sparse_solve():
    from sympy.matrices import SparseMatrix
    A = SparseMatrix(((25, 15, -5), (15, 18, 0), (-5, 0, 11)))
    assert A.cholesky() == Matrix([
        [ 5, 0, 0],
        [ 3, 3, 0],
        [-1, 1, 3]])
    assert A.cholesky() * A.cholesky().T == Matrix([
        [25, 15, -5],
        [15, 18, 0],
        [-5, 0, 11]])

    A = SparseMatrix(((25, 15, -5), (15, 18, 0), (-5, 0, 11)))
    L, D = A.LDLdecomposition()
    assert 15*L == Matrix([
        [15, 0, 0],
        [ 9, 15, 0],
        [-3, 5, 15]])
    assert D == Matrix([
        [25, 0, 0],
        [ 0, 9, 0],
        [ 0, 0, 9]])
    assert L * D * L.T == A

    A = SparseMatrix(((3, 0, 2), (0, 0, 1), (1, 2, 0)))
    assert A.inv() * A == SparseMatrix(eye(3))

    A = SparseMatrix([
        [ 2, -1, 0],
        [-1, 2, -1],
        [ 0, 0, 2]])
    ans = SparseMatrix([
        [S(2)/3, S(1)/3, S(1)/6],
        [S(1)/3, S(2)/3, S(1)/3],
        [     0,      0, S(1)/2]])
    assert A.inv(method='CH') == ans
    assert A.inv(method='LDL') == ans
    assert A * ans == SparseMatrix(eye(3))

    s = A.solve(A[:, 0], 'LDL')
    assert A*s == A[:, 0]
    s = A.solve(A[:, 0], 'CH')
    assert A*s == A[:, 0]
    A = A.col_join(A)
    s = A.solve_least_squares(A[:, 0], 'CH')
    assert A*s == A[:, 0]
    s = A.solve_least_squares(A[:, 0], 'LDL')
    assert A*s == A[:, 0]

    def _syzygies(self):
        """Compute syzygies. See [SCA, algorithm 2.5.4]."""
        # NOTE if self.gens is a standard basis, this can be done more
        #      efficiently using Schreyer's theorem
        from sympy.matrices import eye

        # First bullet point
        k = len(self.gens)
        r = self.rank
        im = eye(k)
        Rkr = self.ring.free_module(r + k)
        newgens = []
        for j, f in enumerate(self.gens):
            m = [0]*(r + k)
            for i, v in enumerate(f):
                m[i] = f[i]
            for i in range(k):
                m[r + i] = im[j, i]
            newgens.append(Rkr.convert(m))
        # Note: we need *descending* order on module index, and TOP=False to
        #       get an eliminetaion order
        F = Rkr.submodule(*newgens, order='ilex', TOP=False)

        # Second bullet point: standard basis of F
        G = F._groebner_vec()

        # Third bullet point: G0 = G intersect the new k components
        G0 = [x[r:] for x in G if all(y == self.ring.convert(0)
                                      for y in x[:r])]

        # Fourth and fifth bullet points: we are done
        return G0

def test_singular_values():
    x = Symbol('x', real=True)

    A = EigenOnlyMatrix([[0, 1*I], [2, 0]])
    # if singular values can be sorted, they should be in decreasing order
    assert A.singular_values() == [2, 1]

    A = eye(3)
    A[1, 1] = x
    A[2, 2] = 5
    vals = A.singular_values()
    # since Abs(x) cannot be sorted, test set equality
    assert set(vals) == set([5, 1, Abs(x)])

    A = EigenOnlyMatrix([[sin(x), cos(x)], [-cos(x), sin(x)]])
    vals = [sv.trigsimp() for sv in A.singular_values()]
    assert vals == [S(1), S(1)]

    A = EigenOnlyMatrix([
        [2, 4],
        [1, 3],
        [0, 0],
        [0, 0]
        ])
    assert A.singular_values() == \
        [sqrt(sqrt(221) + 15), sqrt(15 - sqrt(221))]
    assert A.T.singular_values() == \
        [sqrt(sqrt(221) + 15), sqrt(15 - sqrt(221)), 0, 0]

def test_trace():
    n, m, l = symbols("n m l", integer=True)
    A = MatrixSymbol("A", n, n)
    B = MatrixSymbol("B", n, n)
    assert isinstance(Trace(A), Trace)
    assert not isinstance(Trace(A), MatrixExpr)
    raises(ShapeError, lambda: Trace(MatrixSymbol("B", 3, 4)))
    assert Trace(eye(3)) == 3
    assert Trace(Matrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 9])) == 15

    A / Trace(A)  # Make sure this is possible

    # Some easy simplifications
    assert Trace(Identity(5)) == 5
    assert Trace(ZeroMatrix(5, 5)) == 0
    assert Trace(2 * A * B) == 2 * Trace(A * B)
    assert Trace(A.T) == Trace(A)

    i, j = symbols("i,j")
    F = FunctionMatrix(3, 3, Lambda((i, j), i + j))
    assert Trace(F).doit() == (0 + 0) + (1 + 1) + (2 + 2)

    raises(TypeError, lambda: Trace(S.One))

    assert Trace(A).arg is A

 def sigma(k):
     s = zeros(3)
     e = eye(3)
     for i in range(3):
         for j in range(3):
             s[i, j] = (im if i == k else e)[i,j]
     return s

def demo_corr_bound(n):
    def corr(i, j):
        if i == j:
            return 1
        else:
            return rho
    Rho = Matrix(n, n, corr)
    print "what's the bound of rho to make below a correlation matrix?\n"
    pprint(Rho)
    print "\ncan be viewed as the sum of 2 matrices Rho = A + B:\n"
    A = (1-rho)*eye(n)
    B = rho*ones(n,n)
    pprint(A)
    pprint(B)
    print "\nthe eigen value and its dimention of first matrix A is:"
    pprint(A.eigenvects())
    print "\nas for the seconde matrix B\n"\
            "it's product of any vector v:"
    v = IndexedBase('v')
    vec = Matrix(n, 1, lambda i, j: v[i+1])
    pprint(vec)
    pprint(B*vec)
    print "\nin order for it's to equal to a linear transform of v\n"\
            "we can see its eigen values, vectors and dimentions are:"
    pprint(B.eigenvects())
    print "\nfor any eigen vector of B v, we have Rho*v :\n"
    pprint(Rho*vec)
    print "\nwhich means v is also an eigen vector of Rho,\n"\
            "the eigen values, vectors and dimentions of Rho are:\n"
    pprint(Rho.eigenvects())
    print "\nsince have no negative eigen values <=> positive semidefinite:\n"\
            "the boundaries for rho are: [1/(%d-1),1]" %n

def test_Trace():
    assert isinstance(Trace(A), Trace)
    assert not isinstance(Trace(A), MatrixExpr)
    raises(ShapeError, lambda: Trace(C))
    assert trace(eye(3)) == 3
    assert trace(Matrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 9])) == 15

    assert adjoint(Trace(A)) == trace(Adjoint(A))
    assert conjugate(Trace(A)) == trace(Adjoint(A))
    assert transpose(Trace(A)) == Trace(A)

    A / Trace(A)  # Make sure this is possible

    # Some easy simplifications
    assert trace(Identity(5)) == 5
    assert trace(ZeroMatrix(5, 5)) == 0
    assert trace(OneMatrix(1, 1)) == 1
    assert trace(OneMatrix(2, 2)) == 2
    assert trace(OneMatrix(n, n)) == n
    assert trace(2*A*B) == 2*Trace(A*B)
    assert trace(A.T) == trace(A)

    i, j = symbols('i j')
    F = FunctionMatrix(3, 3, Lambda((i, j), i + j))
    assert trace(F) == (0 + 0) + (1 + 1) + (2 + 2)

    raises(TypeError, lambda: Trace(S.One))

    assert Trace(A).arg is A

    assert str(trace(A)) == str(Trace(A).doit())

    assert Trace(A).is_commutative is True

    def cartan_matrix(self):
        """
        Returns the Cartan matrix for A_n.
        The Cartan matrix matrix for a Lie algebra is
        generated by assigning an ordering to the simple
        roots, (alpha[1], ...., alpha[l]).  Then the ijth
        entry of the Cartan matrix is (<alpha[i],alpha[j]>).

        Examples
        ========

        >>> from sympy.liealgebras.cartan_type import CartanType
        >>> c = CartanType('A4')
        >>> c.cartan_matrix()
        Matrix([
        [ 2, -1,  0,  0],
        [-1,  2, -1,  0],
        [ 0, -1,  2, -1],
        [ 0,  0, -1,  2]])

        """

        n = self.n
        m = 2 * eye(n)
        i = 1
        while i < n-1:
            m[i, i+1] = -1
            m[i, i-1] = -1
            i += 1
        m[0,1] = -1
        m[n-1, n-2] = -1
        return m