def case(S, n_trials=50):        
        S = to_super(S)
        left_dims, right_dims = S.dims
        
        # Assume for the purposes of the test that S maps square operators to square operators.
        in_dim = np.prod(right_dims[0])
        out_dim = np.prod(left_dims[0])
        
        S_dual = to_super(S.dual_chan())
        
        primals = []
        duals = []
    
        for idx_trial in range(n_trials):
            X = rand_dm_ginibre(out_dim)
            X.dims = left_dims
            X = operator_to_vector(X)
            Y = rand_dm_ginibre(in_dim)
            Y.dims = right_dims
            Y = operator_to_vector(Y)

            primals.append((X.dag() * S * Y)[0, 0])
            duals.append((X.dag() * S_dual.dag() * Y)[0, 0])
    
        np.testing.assert_array_almost_equal(primals, duals)

 def test_stinespring_dims(self):
     """
     Stinespring: Check that dims of channels are preserved.
     """
     # FIXME: not the most general test, since this assumes a map
     #        from square matrices to square matrices on the same space.
     chan = super_tensor(to_super(sigmax()), to_super(qeye(3)))
     A, B = to_stinespring(chan)
     assert_equal(A.dims, [[2, 3, 1], [2, 3]])
     assert_equal(B.dims, [[2, 3, 1], [2, 3]])

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_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_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_SuperPreservesSelf(self):
     """
     Superoperator: to_super(q) returns q if q is already a
     supermatrix.
     """
     superop = rand_super()
     assert_(superop is to_super(superop))

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_average_gate_fidelity_target():
    """
    Metrics: Tests that for random unitaries U, AGF(U, U) = 1.
    """
    for _ in range(10):
        U = rand_unitary_haar(13)
        SU = to_super(U)
        assert_almost_equal(average_gate_fidelity(SU, target=U), 1)

def test_tensor_swap_other():
    dims = (2, 3, 4, 5, 7)

    for dim in dims:
        S = to_super(rand_super_bcsz(dim))

        # Swapping the inner indices on a superoperator should give a Choi matrix.
        J = to_choi(S)
        case_tensor_swap(S, [(1, 2)], [[[dim], [dim]], [[dim], [dim]]], J)

def test_tensor_contract_ident():
    qobj = identity([2, 3, 4])
    ans = 3 * identity([2, 4])

    assert_(ans == tensor_contract(qobj, (1, 4)))

    # Now try for superoperators.
    # For now, we just ensure the dims are correct.
    sqobj = to_super(qobj)
    correct_dims = [[[2, 4], [2, 4]], [[2, 4], [2, 4]]]
    assert_equal(correct_dims, tensor_contract(sqobj, (1, 4), (7, 10)).dims)

        def case(map, state):
            S = to_super(map)
            A, B = to_stinespring(map)

            q1 = vector_to_operator(
                S * operator_to_vector(state)
            )
            # FIXME: problem if Kraus index is implicitly
            #        ptraced!
            q2 = (A * state * B.dag()).ptrace((0,))

            assert_((q1 - q2).norm('tr') <= thresh)

    def test_SuperChoiSuper(self):
        """
        Superoperator: Converting superoperator to Choi matrix and back.
        """
        superoperator = rand_super()

        choi_matrix = to_choi(superoperator)
        test_supe = to_super(choi_matrix)

        # Assert both that the result is close to expected, and has the right
        # type.
        assert_((test_supe - superoperator).norm() < tol)
        assert_(choi_matrix.type == "super" and choi_matrix.superrep == "choi")
        assert_(test_supe.type == "super" and test_supe.superrep == "super")

def test_unitarity_known():
    """
    Metrics: Unitarity for known cases.
    """
    def case(q_oper, known_unitarity):
        assert_almost_equal(unitarity(q_oper), known_unitarity)

    yield case, to_super(sigmax()), 1.0
    yield case, sum(map(
        to_super, [qeye(2), sigmax(), sigmay(), sigmaz()]
    )) / 4, 0.0
    yield case, sum(map(
        to_super, [qeye(2), sigmax()]
    )) / 2, 1 / 3.0

def test_dag_preserves_superrep():
    """
    Checks that dag() preserves superrep.
    """

    def case(qobj):
        orig_superrep = qobj.superrep
        assert_equal(qobj.dag().superrep, orig_superrep)

    for dim in (2, 4, 8):
        qobj = rand_super_bcsz(dim)
        yield case, to_super(qobj)
        # These two shouldn't even do anything, since qobj
        # is Hermicity-preserving.
        yield case, to_choi(qobj)
        yield case, to_chi(qobj)

    def test_SuperChoiChiSuper(self):
        """
        Superoperator: Converting two-qubit superoperator through
        Choi and chi representations goes back to right superoperator.
        """
        superoperator = super_tensor(rand_super(2), rand_super(2))

        choi_matrix = to_choi(superoperator)
        chi_matrix = to_chi(choi_matrix)
        test_supe = to_super(chi_matrix)

        # Assert both that the result is close to expected, and has the right
        # type.
        assert_((test_supe - superoperator).norm() < tol)
        assert_(choi_matrix.type == "super" and choi_matrix.superrep == "choi")
        assert_(chi_matrix.type == "super" and chi_matrix.superrep == "chi")
        assert_(test_supe.type == "super" and test_supe.superrep == "super")

    def test_chi_known(self):
        """
        Superoperator: Chi-matrix for known cases is correct.
        """
        def case(S, chi_expected, silent=True):
            chi_actual = to_chi(S)
            chiq = Qobj(chi_expected, dims=[[[2], [2]], [[2], [2]]], superrep='chi')
            if not silent:
                print(chi_actual)
                print(chi_expected)
            assert_almost_equal((chi_actual - chiq).norm('tr'), 0)

        yield case, sigmax(), [
            [0, 0, 0, 0],
            [0, 4, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0]
        ]
        yield case, to_super(sigmax()), [
            [0, 0, 0, 0],
            [0, 4, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0]
        ]
        yield case, qeye(2), [
            [4, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0]
        ]
        yield case, (-1j * sigmax() * pi / 4).expm(), [
            [2, 2j, 0, 0],
            [-2j, 2, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0]
        ]

def test_dnorm_qubit_known_cases():
    """
    Metrics: check agreement for known qubit channels.
    """
    def case(chan1, chan2, expected, significant=4):
        # We again take a generous tolerance so that we don't kill off
        # SCS solvers.
        assert_approx_equal(
            dnorm(chan1, chan2), expected,
            significant=significant
        )

    id_chan = to_choi(qeye(2))
    S_eye = to_super(id_chan)
    X_chan = to_choi(sigmax())
    depol = to_choi(Qobj(
        diag(ones((4,))),
        dims=[[[2], [2]], [[2], [2]]], superrep='chi'
    ))
    S_H = to_super(hadamard_transform())

    W = swap()

    # We need to restrict the number of iterations for things on the boundary,
    # such as perfectly distinguishable channels.
    yield case, id_chan, X_chan, 2
    yield case, id_chan, depol, 1.5

    # Next, we'll generate some test cases based on comparisons to pre-existing
    # dnorm() implementations. In particular, the targets for the following
    # test cases were generated using QuantumUtils for MATLAB (https://goo.gl/oWXhO9).

    def overrotation(x):
        return to_super((1j * np.pi * x * sigmax() / 2).expm())

    for x, target in {
        1.000000e-03: 3.141591e-03,
        3.100000e-03: 9.738899e-03,
        1.000000e-02: 3.141463e-02,
        3.100000e-02: 9.735089e-02,
        1.000000e-01: 3.128689e-01,
        3.100000e-01: 9.358596e-01
    }.items():
        yield case, overrotation(x), id_chan, target

    def had_mixture(x):
        return (1 - x) * S_eye + x * S_H

    for x, target in {
        1.000000e-03: 2.000000e-03,
        3.100000e-03: 6.200000e-03,
        1.000000e-02: 2.000000e-02,
        3.100000e-02: 6.200000e-02,
        1.000000e-01: 2.000000e-01,
        3.100000e-01: 6.200000e-01
    }.items():
        yield case, had_mixture(x), id_chan, target

    def swap_map(x):
        S = (1j * x * W).expm()
        S._type = None
        S.dims = [[[2], [2]], [[2], [2]]]
        S.superrep = 'super'
        return S

    for x, target in {
        1.000000e-03: 2.000000e-03,
        3.100000e-03: 6.199997e-03,
        1.000000e-02: 1.999992e-02,
        3.100000e-02: 6.199752e-02,
        1.000000e-01: 1.999162e-01,
        3.100000e-01: 6.173918e-01
    }.items():
        yield case, swap_map(x), id_chan, target

    # Finally, we add a known case from Johnston's QETLAB documentation,
    # || Phi - I ||,_♢ where Phi(X) = UXU⁺ and U = [[1, 1], [-1, 1]] / sqrt(2).
    yield case, Qobj([[1, 1], [-1, 1]]) / np.sqrt(2), qeye(2), np.sqrt(2)

 def overrotation(x):
     return to_super((1j * np.pi * x * sigmax() / 2).expm())

#.........这里部分代码省略.........
        Diamond norm of q_oper.

    Raises
    ------
    ImportError
        If CVXPY cannot be imported.

    .. _cvxpy: http://www.cvxpy.org/en/latest/
    """
    if cvxpy is None:  # pragma: no cover
        raise ImportError("dnorm() requires CVXPY to be installed.")

    # We follow the strategy of using Watrous' simpler semidefinite
    # program in its primal form. This is the same strategy used,
    # for instance, by both pyGSTi and SchattenNorms.jl. (By contrast,
    # QETLAB uses the dual problem.)

    # Check if A and B are both unitaries. If so, then we can without
    # loss of generality choose B to be the identity by using the
    # unitary invariance of the diamond norm,
    #     || A - B ||_♢ = || A B⁺ - I ||_♢.
    # Then, using the technique mentioned by each of Johnston and
    # da Silva,
    #     || A B⁺ - I ||_♢ = max_{i, j} | \lambda_i(A B⁺) - \lambda_j(A B⁺) |,
    # where \lambda_i(U) is the ith eigenvalue of U.

    if (
        # There's a lot of conditions to check for this path.
        not force_solve and B is not None and
        # Only check if they aren't superoperators.
        A.type == "oper" and B.type == "oper" and
        # The difference of unitaries optimization is currently
        # only implemented for d == 2. Much of the code below is more general,
        # though, in anticipation of generalizing the optimization.
        A.shape[0] == 2
    ):
        # Make an identity the same size as A and B to
        # compare against.
        I = qeye(A.dims[0])
        # Compare to B first, so that an error is raised
        # as soon as possible.
        Bd = B.dag()
        if (
            (B * Bd - I).norm() < 1e-6 and
            (A * A.dag() - I).norm() < 1e-6
        ):
            # Now we are on the fast path, so let's compute the
            # eigenvalues, then find the diameter of the smallest circle
            # containing all of them.
            #
            # For now, this is only implemented for dim = 2, such that
            # generalizing here will allow for generalizing the optimization.
            # A reasonable approach would probably be to use Welzl's algorithm
            # (https://en.wikipedia.org/wiki/Smallest-circle_problem).
            U = A * B.dag()
            eigs = U.eigenenergies()
            eig_distances = np.abs(eigs[:, None] - eigs[None, :])
            return np.max(eig_distances)

    # Force the input superoperator to be a Choi matrix.
    J = to_choi(A)
    
    if B is not None:
        J -= to_choi(B)

    # Watrous 2012 also points out that the diamond norm of Lambda
    # is the same as the completely-bounded operator-norm (∞-norm)
    # of the dual map of Lambda. We can evaluate that norm much more
    # easily if Lambda is completely positive, since then the largest
    # eigenvalue is the same as the largest singular value.
    
    if not force_solve and J.iscp:
        S_dual = to_super(J.dual_chan())
        vec_eye = operator_to_vector(qeye(S_dual.dims[1][1]))
        op = vector_to_operator(S_dual * vec_eye)
        # The 2-norm was not implemented for sparse matrices as of the time
        # of this writing. Thus, we must yet again go dense.
        return la.norm(op.data.todense(), 2)
    
    # If we're still here, we need to actually solve the problem.

    # Assume square...
    dim = np.prod(J.dims[0][0])
    
    # The constraints only depend on the dimension, so
    # we can cache them efficiently.
    problem, Jr, Ji, X, rho0, rho1 = dnorm_problem(dim)
    
    # Load the parameters with the Choi matrix passed in.
    J_dat = J.data
    
    Jr.value = sp.csr_matrix((J_dat.data.real, J_dat.indices, J_dat.indptr), 
                             shape=J_dat.shape)
   
    Ji.value = sp.csr_matrix((J_dat.data.imag, J_dat.indices, J_dat.indptr),
                             shape=J_dat.shape)
    # Finally, set up and run the problem.
    problem.solve(solver=solver, verbose=verbose)
    
    return problem.value

    def test_known_iscptp(self):
        """
        Superoperator: ishp, iscp, istp and iscptp known cases.
        """
        def case(qobj, shouldhp, shouldcp, shouldtp):
            hp = qobj.ishp
            cp = qobj.iscp
            tp = qobj.istp
            cptp = qobj.iscptp

            shouldcptp = shouldcp and shouldtp

            if (
                hp == shouldhp and
                cp == shouldcp and
                tp == shouldtp and
                cptp == shouldcptp
            ):
                return

            fails = []
            if hp != shouldhp:
                fails.append(("ishp", shouldhp, hp))
            if tp != shouldtp:
                fails.append(("istp", shouldtp, tp))
            if cp != shouldcp:
                fails.append(("iscp", shouldcp, cp))
            if cptp != shouldcptp:
                fails.append(("iscptp", shouldcptp, cptp))

            raise AssertionError("Expected {}.".format(" and ".join([
                "{} == {} (got {})".format(fail, expected, got)
                for fail, expected, got in fails
            ])))

        # Conjugation by a creation operator should
        # have be CP (and hence HP), but not TP.
        a = create(2).dag()
        S = sprepost(a, a.dag())
        case(S, True, True, False)

        # A single off-diagonal element should not be CP,
        # nor even HP.
        S = sprepost(a, a)
        case(S, False, False, False)
        
        # Check that unitaries are CPTP and HP.
        case(identity(2), True, True, True)
        case(sigmax(), True, True, True)

        # Check that unitaries on bipartite systems are CPTP and HP.
        case(tensor(sigmax(), identity(2)), True, True, True)

        # Check that a linear combination of bipartitie unitaries is CPTP and HP.
        S = (
            to_super(tensor(sigmax(), identity(2))) + to_super(tensor(identity(2), sigmay()))
        ) / 2
        case(S, True, True, True)

        # The partial transpose map, whose Choi matrix is SWAP, is TP
        # and HP but not CP (one negative eigenvalue).
        W = Qobj(swap(), type='super', superrep='choi')
        case(W, True, False, True)

        # Subnormalized maps (representing erasure channels, for instance)
        # can be CP but not TP.
        subnorm_map = Qobj(identity(4) * 0.9, type='super', superrep='super')
        case(subnorm_map, True, True, False)

        # Check that things which aren't even operators aren't identified as
        # CPTP.
        case(basis(2), False, False, False)