def test_super_tensor_property():
    """
    Tensor: Super_tensor correctly tensors on underlying spaces.
    """
    U1 = rand_unitary(3)
    U2 = rand_unitary(5)

    U = tensor(U1, U2)
    S_tens = to_super(U)

    S_supertens = super_tensor(to_super(U1), to_super(U2))

    assert_(S_tens == S_supertens)
    assert_equal(S_supertens.superrep, 'super')

def test_call():
    """
    Test Qobj: Call
    """
    # Make test objects.
    psi = rand_ket(3)
    rho = rand_dm_ginibre(3)
    U = rand_unitary(3)
    S = rand_super_bcsz(3)

    # Case 0: oper(ket).
    assert U(psi) == U * psi

    # Case 1: oper(oper). Should raise TypeError.
    with expect_exception(TypeError):
        U(rho)

    # Case 2: super(ket).
    assert S(psi) == vector_to_operator(S * operator_to_vector(ket2dm(psi)))

    # Case 3: super(oper).
    assert S(rho) == vector_to_operator(S * operator_to_vector(rho))

    # Case 4: super(super). Should raise TypeError.
    with expect_exception(TypeError):
        S(S)

def test_trunc_neg():
    """
    Test Qobj: Checks trunc_neg in several different cases.
    """

    @has_description
    def case(qobj, method, expected=None):
        pos_qobj = qobj.trunc_neg(method=method)
        assert(all([energy > -1e-8 for energy in pos_qobj.eigenenergies()]))
        assert_almost_equal(pos_qobj.tr(), 1)

        if expected is not None:
            assert_almost_equal(pos_qobj.data.todense(), expected.data.todense())

    for method in ('clip', 'sgs'):
        # Make sure that it works for operators that are already positive.
        yield case("Test Qobj: trunc_neg works for positive opers."), \
            rand_dm(5), method
        # Make sure that it works for a diagonal matrix.
        yield case("Test Qobj: trunc_neg works for diagonal opers."), \
            Qobj(np.diag([1.1, -0.1])), method, Qobj(np.diag([1.0, 0.0]))
        # Make sure that it works for a non-diagonal matrix.
        U = rand_unitary(3)
        yield case("Test Qobj: trunc_neg works for non-diagonal opers."), \
            U * Qobj(np.diag([1.1, 0, -0.1])) * U.dag(), \
            method, \
            U * Qobj(np.diag([1.0, 0.0, 0.0])) * U.dag()

    # Check the test case in SGS.
    yield (
        case("Test Qobj: trunc_neg works for SGS known-good test case."),
        Qobj(np.diag([3. / 5, 1. / 2, 7. / 20, 1. / 10, -11. / 20])), 'sgs',
        Qobj(np.diag([9. / 20, 7. / 20, 1. / 5, 0, 0]))
    )

def test_composite_oper():
    """
    Composite: Tests compositing unitaries and superoperators.
    """
    U1 = rand_unitary(3)
    U2 = rand_unitary(5)
    S1 = to_super(U1)
    S2 = to_super(U2)

    S3 = rand_super(4)
    S4 = rand_super(7)

    assert_(composite(U1, U2) == tensor(U1, U2))
    assert_(composite(S3, S4) == super_tensor(S3, S4))
    assert_(composite(U1, S4) == super_tensor(S1, S4))
    assert_(composite(S3, U2) == super_tensor(S3, S2))

    def test_ComplexSingleApply(self):
        """
        Composite system, operator on Hilbert space.
        """
        tol = 1e-12
        rho_list = list(map(rand_dm, [2, 3, 2, 3, 2]))
        rho_input = tensor(rho_list)
        single_op = rand_unitary(3)

        analytic_result = rho_list
        analytic_result[1] = single_op * analytic_result[1] * single_op.dag()
        analytic_result[3] = single_op * analytic_result[3] * single_op.dag()
        analytic_result = tensor(analytic_result)

        naive_result = subsystem_apply(rho_input, single_op, [False, True, False, True, False], reference=True)
        naive_diff = (analytic_result - naive_result).data.todense()
        naive_diff_norm = norm(naive_diff)
        assert_(
            naive_diff_norm < tol,
            msg="ComplexSingle: naive_diff_norm {} " "is beyond tolerance {}".format(naive_diff_norm, tol),
        )

        efficient_result = subsystem_apply(rho_input, single_op, [False, True, False, True, False])
        efficient_diff = (efficient_result - analytic_result).data.todense()
        efficient_diff_norm = norm(efficient_diff)
        assert_(
            efficient_diff_norm < tol,
            msg="ComplexSingle: efficient_diff_norm {} " "is beyond tolerance {}".format(efficient_diff_norm, tol),
        )

def test_CheckMulType():
    "Qobj multiplication type"

    # ket-bra and bra-ket multiplication
    psi = basis(5)
    dm = psi * psi.dag()
    assert_(dm.isoper)
    assert_(dm.isherm)

    nrm = psi.dag() * psi
    assert_equal(np.prod(nrm.shape), 1)
    assert_((abs(nrm) == 1)[0, 0])

    # operator-operator multiplication
    H1 = rand_herm(3)
    H2 = rand_herm(3)
    out = H1 * H2
    assert_(out.isoper)
    out = H1 * H1
    assert_(out.isoper)
    assert_(out.isherm)
    out = H2 * H2
    assert_(out.isoper)
    assert_(out.isherm)

    U = rand_unitary(5)
    out = U.dag() * U
    assert_(out.isoper)
    assert_(out.isherm)

    N = num(5)

    out = N * N
    assert_(out.isoper)
    assert_(out.isherm)

    # operator-ket and bra-operator multiplication
    op = sigmax()
    ket1 = basis(2)
    ket2 = op * ket1
    assert_(ket2.isket)

    bra1 = basis(2).dag()
    bra2 = bra1 * op
    assert_(bra2.isbra)

    assert_(bra2.dag() == ket2)

    # superoperator-operket and operbra-superoperator multiplication
    sop = to_super(sigmax())
    opket1 = operator_to_vector(fock_dm(2))
    opket2 = sop * opket1
    assert(opket2.isoperket)

    opbra1 = operator_to_vector(fock_dm(2)).dag()
    opbra2 = opbra1 * sop
    assert(opbra2.isoperbra)

    assert_(opbra2.dag() == opket2)

def test_QobjPermute():
    "Qobj permute"
    A = basis(3, 0)
    B = basis(5, 4)
    C = basis(4, 2)
    psi = tensor(A, B, C)
    psi2 = psi.permute([2, 0, 1])
    assert_(psi2 == tensor(C, A, B))
    
    psi_bra = psi.dag()
    psi2_bra = psi_bra.permute([2, 0, 1])
    assert_(psi2_bra == tensor(C, A, B).dag())

    A = fock_dm(3, 0)
    B = fock_dm(5, 4)
    C = fock_dm(4, 2)
    rho = tensor(A, B, C)
    rho2 = rho.permute([2, 0, 1])
    assert_(rho2 == tensor(C, A, B))

    for ii in range(3):
        A = rand_ket(3)
        B = rand_ket(4)
        C = rand_ket(5)
        psi = tensor(A, B, C)
        psi2 = psi.permute([1, 0, 2])
        assert_(psi2 == tensor(B, A, C))
        
        psi_bra = psi.dag()
        psi2_bra = psi_bra.permute([1, 0, 2])
        assert_(psi2_bra == tensor(B, A, C).dag())

    for ii in range(3):
        A = rand_dm(3)
        B = rand_dm(4)
        C = rand_dm(5)
        rho = tensor(A, B, C)
        rho2 = rho.permute([1, 0, 2])
        assert_(rho2 == tensor(B, A, C))
        
        rho_vec = operator_to_vector(rho)
        rho2_vec = rho_vec.permute([[1, 0, 2],[4,3,5]])
        assert_(rho2_vec == operator_to_vector(tensor(B, A, C)))
        
        rho_vec_bra = operator_to_vector(rho).dag()
        rho2_vec_bra = rho_vec_bra.permute([[1, 0, 2],[4,3,5]])
        assert_(rho2_vec_bra == operator_to_vector(tensor(B, A, C)).dag())
        
    for ii in range(3):
        super_dims = [3, 5, 4]
        U = rand_unitary(np.prod(super_dims), density=0.02, dims=[super_dims, super_dims])
        Unew = U.permute([2,1,0])
        S_tens = to_super(U)
        S_tens_new = to_super(Unew)
        assert_(S_tens_new == S_tens.permute([[2,1,0],[5,4,3]]))

def test_fidelity2():
    """
    Metrics: Fidelity, invariance under unitary trans.
    """
    for k in range(10):
        rho1 = rand_dm(25, 0.25)
        rho2 = rand_dm(25, 0.25)
        U = rand_unitary(25, 0.25)
        F = fidelity(rho1, rho2)
        FU = fidelity(U*rho1*U.dag(), U*rho2*U.dag())
        assert_(abs((F-FU)/F) < 1e-5)

def test_tracedist1():
    """
    Metrics: Trace dist., invariance under unitary trans.
    """
    for k in range(10):
        rho1 = rand_dm(25, 0.25)
        rho2 = rand_dm(25, 0.25)
        U = rand_unitary(25, 0.25)
        D = tracedist(rho1, rho2)
        DU = tracedist(U*rho1*U.dag(), U*rho2*U.dag())
        assert_(abs((D-DU)/D) < 1e-5)

 def test_SimpleSingleApply(self):
     """
     Non-composite system, operator on Hilbert space.
     """
     rho_3 = rand_dm(3)
     single_op = rand_unitary(3)
     analytic_result = single_op * rho_3 * single_op.dag()
     naive_result = subsystem_apply(rho_3, single_op, [True],
                                    reference=True)
     efficient_result = subsystem_apply(rho_3, single_op, [True])
     naive_diff = (analytic_result - naive_result).data.todense()
     efficient_diff = (efficient_result - analytic_result).data.todense()
     assert_(norm(naive_diff) < 1e-12 and norm(efficient_diff) < 1e-12)

def test_zcsr_transpose():
    "spmath: zcsr_transpose"
    for k in range(50):
        ra = np.random.randint(2,100)
        A = rand_ket(ra,0.5).data
        B = A.T.tocsr()
        C = A.trans()
        x = np.all(B.data == C.data)
        y = np.all(B.indices == C.indices)
        z = np.all(B.indptr == C.indptr)
        assert_(x*y*z)
    
    for k in range(50):
        ra = np.random.randint(2,100)
        A = rand_herm(5,1.0/ra).data
        B = A.T.tocsr()
        C = A.trans()
        x = np.all(B.data == C.data)
        y = np.all(B.indices == C.indices)
        z = np.all(B.indptr == C.indptr)
        assert_(x*y*z)
        
    for k in range(50):
        ra = np.random.randint(2,100)
        A = rand_dm(ra,1.0/ra).data
        B = A.T.tocsr()
        C = A.trans()
        x = np.all(B.data == C.data)
        y = np.all(B.indices == C.indices)
        z = np.all(B.indptr == C.indptr)
        assert_(x*y*z)
        
    for k in range(50):
        ra = np.random.randint(2,100)
        A = rand_unitary(ra,1.0/ra).data
        B = A.T.tocsr()
        C = A.trans()
        x = np.all(B.data == C.data)
        y = np.all(B.indices == C.indices)
        z = np.all(B.indptr == C.indptr)
        assert_(x*y*z)

    def test_ComplexSingleApply(self):
        """
        Composite system, operator on Hilbert space.
        """
        rho_list = list(map(rand_dm, [2, 3, 2, 3, 2]))
        rho_input = tensor(rho_list)
        single_op = rand_unitary(3)

        analytic_result = rho_list
        analytic_result[1] = single_op * analytic_result[1] * single_op.dag()
        analytic_result[3] = single_op * analytic_result[3] * single_op.dag()
        analytic_result = tensor(analytic_result)

        naive_result = subsystem_apply(rho_input, single_op,
                                       [False, True, False, True, False],
                                       reference=True)
        naive_diff = (analytic_result - naive_result).data.todense()
        assert_(norm(naive_diff) < 1e-12)

        efficient_result = subsystem_apply(rho_input, single_op,
                                           [False, True, False, True, False])
        efficient_diff = (efficient_result - analytic_result).data.todense()
        assert_(norm(efficient_diff) < 1e-12)

    def test_SimpleSingleApply(self):
        """
        Non-composite system, operator on Hilbert space.
        """
        tol = 1e-12
        rho_3 = rand_dm(3)
        single_op = rand_unitary(3)
        analytic_result = single_op * rho_3 * single_op.dag()
        naive_result = subsystem_apply(rho_3, single_op, [True], reference=True)
        naive_diff = (analytic_result - naive_result).data.todense()
        naive_diff_norm = norm(naive_diff)
        assert_(
            naive_diff_norm < tol,
            msg="SimpleSingle: naive_diff_norm {} " "is beyond tolerance {}".format(naive_diff_norm, tol),
        )

        efficient_result = subsystem_apply(rho_3, single_op, [True])
        efficient_diff = (efficient_result - analytic_result).data.todense()
        efficient_diff_norm = norm(efficient_diff)
        assert_(
            efficient_diff_norm < tol,
            msg="SimpleSingle: efficient_diff_norm {} " "is beyond tolerance {}".format(efficient_diff_norm, tol),
        )