def test_known_iscptp(self):
        """
        Superoperator: iscp, istp and iscptp known cases.
        """
        # Check that unitaries are CPTP.
        assert_(identity(2).iscptp)
        assert_(sigmax().iscptp)
        
        # The partial transpose map, whose Choi matrix is SWAP, is TP but not
        # CP.
        W = Qobj(swap(), type='super', superrep='choi')
        assert_(W.istp)
        assert_(not W.iscp)
        assert_(not W.iscptp)
        
        # Subnormalized maps (representing erasure channels, for instance)
        # can be CP but not TP.
        subnorm_map = Qobj(identity(4) * 0.9, type='super', superrep='super')
        assert_(subnorm_map.iscp)
        assert_(not subnorm_map.istp)
        assert_(not subnorm_map.iscptp)

        
        # Check that things which aren't even operators aren't identified as
        # CPTP.
        assert_(not (basis(2).iscptp))

def gate_expand_1toN(U, N, target):
    """
    Create a Qobj representing a one-qubit gate that act on a system with N
    qubits.

    Parameters
    ----------
    U : Qobj
        The one-qubit gate

    N : integer
        The number of qubits in the target space.

    target : integer
        The index of the target qubit.

    Returns
    -------
    gate : qobj
        Quantum object representation of N-qubit gate.
    
    """

    if N < 1:
        raise ValueError("integer N must be larger or equal to 1")

    if target >= N:
        raise ValueError("target must be integer < integer N")

    return tensor([identity(2)] * (target) + [U] +
                  [identity(2)] * (N - target - 1))

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 __init__(self, N, correct_global_phase=True,
                 sx=None, sz=None, sxsy=None):
        """
        Parameters
        ----------
        sx: Integer/List
            The delta for each of the qubits in the system.

        sz: Integer/List
            The epsilon for each of the qubits in the system.

        sxsy: Integer/List
            The interaction strength for each of the qubit pair in the system.
        """

        super(SpinChain, self).__init__(N, correct_global_phase)

        self.sx_ops = [tensor([sigmax() if m == n else identity(2)
                               for n in range(N)])
                       for m in range(N)]
        self.sz_ops = [tensor([sigmaz() if m == n else identity(2)
                               for n in range(N)])
                       for m in range(N)]

        self.sxsy_ops = []
        for n in range(N - 1):
            x = [identity(2)] * N
            x[n] = x[n + 1] = sigmax()
            y = [identity(2)] * N
            y[n] = y[n + 1] = sigmay()
            self.sxsy_ops.append(tensor(x) + tensor(y))

        if sx is None:
            self.sx_coeff = [0.25 * 2 * np.pi] * N
        elif not isinstance(sx, list):
            self.sx_coeff = [sx * 2 * np.pi] * N
        else:
            self.sx_coeff = sx

        if sz is None:
            self.sz_coeff = [1.0 * 2 * np.pi] * N
        elif not isinstance(sz, list):
            self.sz_coeff = [sz * 2 * np.pi] * N
        else:
            self.sz_coeff = sz

        if sxsy is None:
            self.sxsy_coeff = [0.1 * 2 * np.pi] * (N - 1)
        elif not isinstance(sxsy, list):
            self.sxsy_coeff = [sxsy * 2 * np.pi] * (N - 1)
        else:
            self.sxsy_coeff = sxsy

def f(ab, w, ulim, Nt, tmax):
    # Create time series vector
    time = np.linspace(0, tmax, Nt)
    # U = identity(2*2*Nho)
    # create starting unitary
    # could this be wrong ?
    U = tensor(identity(2), identity(2), identity(Nho))

    # Initialize g
    g = 0.0
    #  half the length of ab vector
    # used for decomposing it into a and b
    # possible problems for odd length vectors
    # no because created by concatenating two vectors of same length
    # in whatever program that calls it
    Lab = int(ab.size/2)
    #  at a paticular time t
    # making g
    # making H,Ui
    # and hence U at that time
    # multiplying it to the product of unitatries at previous times
    for t in time:
        for ir in range(Lab):
            # normalizing the A_k and B_k
            norm = np.sqrt(ab[ir]**2 + ab[ir+Lab]**2)
            # numerator built by decomposition of ab just like above
            numer = ab[ir]*np.sin(w[ir]*t) + ab[ir+Lab]*np.cos(w[ir]*t)
            # multiplying by ulim and dividing by norm
            # construction of g_k
            # adding it to rest of g
            g += (ulim*numer)/norm

        # constructing the total hamiltonian
        H = H0 + g*Hi
        # constructing U at the time t
        Ui = (-1j*H*t).expm()
        # multiplying it with the othe unitaries before
        U *= Ui

    # f = hilbert_dist(U.dag(),U_target)
    # fi = np.float64(hilbert_dist(U.dag(),U_target))
    # fif = hilbert_dist(U.dag(),U_target)
    # fis = str(fif)
    # fi = float(fis)
    # stuff inside the trace with normalization
    uu = ((U.dag())*U_target)/(2*2*Nho)
    # taking the trace
    f = uu.tr()
    # absolute value of the trace
    fi = abs(f)
    return fi

 def test_choi_tr(self):
     """
     Superoperator: Checks that the trace of matrices returned by to_choi
     matches that asserted by the docstring for that function.
     """
     for dims in range(2, 5):
         assert_(abs(to_choi(identity(dims)).tr() - dims) <= tol)

def controlled_gate(U, N=2, control=0, target=1, control_value=1):
    """
    Create an N-qubit controlled gate from a single-qubit gate U with the given
    control and target qubits.

    Parameters
    ----------
    U : Qobj
        Arbitrary single-qubit gate.

    N : integer
        The number of qubits in the target space.

    control : integer
        The index of the first control qubit.

    target : integer
        The index of the target qubit.

    control_value : integer (1)
        The state of the control qubit that activates the gate U.

    Returns
    -------
    result : qobj
        Quantum object representing the controlled-U gate.
    
    """

    if [N, control, target] == [2, 0, 1]:
        return (tensor(fock_dm(2, control_value), U) +
                tensor(fock_dm(2, 1 - control_value), identity(2)))
    else:
        U2 = controlled_gate(U, control_value=control_value)
        return gate_expand_2toN(U2, N=N, control=control, target=target)

 def test_average_gate_fidelity(self):
     """
     Metrics: Checks that average gate fidelities are sensible for random
     maps, and are equal to 1 for identity maps.
     """
     for dims in range(2, 5):
         assert_(abs(average_gate_fidelity(identity(dims)) - 1) <= 1e-12)
     assert_(0 <= average_gate_fidelity(self.rand_super()) <= 1)

def test_average_gate_fidelity():
    """
    Metrics: Check avg gate fidelities for random
    maps (equal to 1 for id maps).
    """
    for dims in range(2, 5):
        assert_(abs(average_gate_fidelity(identity(dims)) - 1) <= 1e-12)
    assert_(0 <= average_gate_fidelity(rand_super()) <= 1)

def gate_expand_2toN(U, N, control=None, target=None, targets=None):
    """
    Create a Qobj representing a two-qubit gate that act on a system with N
    qubits.

    Parameters
    ----------
    U : Qobj
        The two-qubit gate

    N : integer
        The number of qubits in the target space.

    control : integer
        The index of the control qubit.

    target : integer
        The index of the target qubit.

    targets : list
        List of target qubits.

    Returns
    -------
    gate : qobj
        Quantum object representation of N-qubit gate.
    
    """

    if targets is not None:
        control, target = targets

    if control is None or target is None:
        raise ValueError("Specify value of control and target")

    if N < 2:
        raise ValueError("integer N must be larger or equal to 2")

    if control >= N or target >= N:
        raise ValueError("control and not target must be integer < integer N")

    if control == target:
        raise ValueError("target and not control cannot be equal")

    p = list(range(N))

    if target == 0 and control == 1:
        p[control], p[target] = p[target], p[control]

    elif target == 0:
        p[1], p[target] = p[target], p[1]
        p[1], p[control] = p[control], p[1]

    else:
        p[1], p[target] = p[target], p[1]
        p[0], p[control] = p[control], p[0]

    return tensor([U] + [identity(2)] * (N - 2)).permute(p)

    def testExpandGate3toN_permutation(self):
        """
        gates: expand 3 to 3 with permuTation (using toffoli)
        """
        for _p in itertools.permutations([0, 1, 2]):
            controls, target = [_p[0], _p[1]], _p[2]
            
            controls = [1, 2]
            target = 0

            p = [1, 2, 3]
            p[controls[0]] = 0
            p[controls[1]] = 1
            p[target] = 2

            U = toffoli(N=3, controls=controls, target=target)
            
            ops = [basis(2, 0).dag(),  basis(2, 0).dag(), identity(2)]
            P = tensor(ops[p[0]], ops[p[1]], ops[p[2]])
            assert_(P * U * P.dag() == identity(2))

            ops = [basis(2, 1).dag(),  basis(2, 0).dag(), identity(2)]
            P = tensor(ops[p[0]], ops[p[1]], ops[p[2]])
            assert_(P * U * P.dag() == identity(2))

            ops = [basis(2, 0).dag(),  basis(2, 1).dag(), identity(2)]
            P = tensor(ops[p[0]], ops[p[1]], ops[p[2]])
            assert_(P * U * P.dag() == identity(2))

            ops = [basis(2, 1).dag(),  basis(2, 1).dag(), identity(2)]
            P = tensor(ops[p[0]], ops[p[1]], ops[p[2]])
            assert_(P * U * P.dag() == sigmax())

    def case_is_clifford(self, U):
        paulis = (identity(2), sigmax(), sigmay(), sigmaz())

        for P in paulis:
            U_P = U * P * U.dag()
            
            out = (np.any(
                np.array([self._prop_identity(U_P * Q) for Q in paulis])
            ))
        return out

    def case_is_clifford(self, U):
        paulis = (identity(2), sigmax(), sigmay(), sigmaz())

        for P in paulis:
            U_P = U * P * U.dag()
            
            assert_(any(
                self._prop_identity(U_P * Q)
                for Q in paulis
            ))

def cphase(theta, N=2, control=0, target=1):
    """
    Returns quantum object representing the phase shift gate.

    Parameters
    ----------
    theta : float
        Phase rotation angle.

    N : integer
        The number of qubits in the target space.

    control : integer
        The index of the control qubit.

    target : integer
        The index of the target qubit.

    Returns
    -------
    U : qobj
        Quantum object representation of controlled phase gate.
    """

    if N < 1 or target < 0 or control < 0:
        raise ValueError("Minimum value: N=1, control=0 and target=0")

    if control >= N or target >= N:
        raise ValueError("control and target need to be smaller than N")

    U_list1 = [identity(2)] * N
    U_list2 = [identity(2)] * N

    U_list1[control] = fock_dm(2, 1)
    U_list1[target] = phasegate(theta)

    U_list2[control] = fock_dm(2, 0)

    U = tensor(U_list1) + tensor(U_list2)
    return U

 def test_choi_tr(self):
     """
     Superoperator: Trace returned by to_choi matches docstring.
     """
     for dims in range(2, 5):
         assert_(abs(to_choi(identity(dims)).tr() - dims) <= tol)

    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())
        yield case, S, True, True, False

        # A single off-diagonal element should not be CP,
        # nor even HP.
        S = sprepost(a, a)
        yield case, S, False, False, False
        
        # Check that unitaries are CPTP and HP.
        yield case, identity(2), True, True, True
        yield case, sigmax(), 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')
        yield 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')
        yield case, subnorm_map, True, True, False

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

    CPTP, expand to arbitrary dimensional systems, etc.
    """
    return Qobj(dims=[[[2], [2]], [[2], [2]]],
                inpt=array([[1. - pe / 2., 0., 0., 1. - pe],
                            [0., pe / 2., 0., 0.],
                            [0., 0., pe / 2., 0.],
                            [1. - pe, 0., 0., 1. - pe / 2.]]),
                superrep='choi')


# CHANGE OF BASIS FUNCTIONS ---------------------------------------------------
# These functions find change of basis matrices, and are useful in converting
# between (for instance) Choi and chi matrices. At some point, these should
# probably be moved out to another module.

_SINGLE_QUBIT_PAULI_BASIS = (identity(2), sigmax(), sigmay(), sigmaz())


def _pauli_basis(nq=1):
    # NOTE: This is slow as can be.
    # TODO: Make this sparse. CSR format was causing problems for the [idx, :]
    #       slicing below.
    B = zeros((4 ** nq, 4 ** nq), dtype=complex)
    dims = [[[2] * nq] * 2] * 2

    for idx, op in enumerate(starmap(tensor,
                                     product(_SINGLE_QUBIT_PAULI_BASIS,
                                             repeat=nq))):
        B[:, idx] = operator_to_vector(op).dag().data.todense()

    return Qobj(B, dims=dims)

#.........这里部分代码省略.........
        p[1], p[controls[1]] = p[controls[1]], p[1]

    elif controls[1] == 1 and target == 2:
        p[0], p[controls[0]] = p[controls[0]], p[0]

    elif controls[0] == 1 and controls[1] == 0:
        p[controls[1]], p[controls[0]] = p[controls[0]], p[controls[1]]
        p2[2], p2[target] = p2[target], p2[2]
        p = [p2[p[k]] for k in range(N)]

    elif controls[0] == 2 and target == 0:
        p[target], p[controls[0]] = p[controls[0]], p[target]
        p1[1], p1[controls[1]] = p1[controls[1]], p1[1]
        p = [p1[p[k]] for k in range(N)]

    elif controls[1] == 2 and target == 1:
        p[target], p[controls[1]] = p[controls[1]], p[target]
        p1[0], p1[controls[0]] = p1[controls[0]], p1[0]
        p = [p1[p[k]] for k in range(N)]

    elif controls[0] == 1 and controls[1] == 2:
        #  controls[0] -> controls[1] -> target -> outside
        p[0], p[1] = p[1], p[0]
        p[0], p[2] = p[2], p[0]
        p[0], p[target] = p[target], p[0]

    elif controls[0] == 2 and target == 1:
        #  controls[0] -> target -> controls[1] -> outside
        p[0], p[2] = p[2], p[0]
        p[0], p[1] = p[1], p[0]
        p[0], p[controls[1]] = p[controls[1]], p[0]

    elif controls[1] == 0 and controls[0] == 2:
        #  controls[1] -> controls[0] -> target -> outside
        p[1], p[0] = p[0], p[1]
        p[1], p[2] = p[2], p[1]
        p[1], p[target] = p[target], p[1]

    elif controls[1] == 2 and target == 0:
        #  controls[1] -> target -> controls[0] -> outside
        p[1], p[2] = p[2], p[1]
        p[1], p[0] = p[0], p[1]
        p[1], p[controls[0]] = p[controls[0]], p[1]

    elif target == 1 and controls[1] == 0:
        #  target -> controls[1] -> controls[0] -> outside
        p[2], p[1] = p[1], p[2]
        p[2], p[0] = p[0], p[2]
        p[2], p[controls[0]] = p[controls[0]], p[2]

    elif target == 0 and controls[0] == 1:
        #  target -> controls[0] -> controls[1] -> outside
        p[2], p[0] = p[0], p[2]
        p[2], p[1] = p[1], p[2]
        p[2], p[controls[1]] = p[controls[1]], p[2]

    elif controls[0] == 0 and controls[1] == 2:
        #  controls[0] -> self, controls[1] -> target -> outside
        p[1], p[2] = p[2], p[1]
        p[1], p[target] = p[target], p[1]

    elif controls[1] == 1 and controls[0] == 2:
        #  controls[1] -> self, controls[0] -> target -> outside
        p[0], p[2] = p[2], p[0]
        p[0], p[target] = p[target], p[0]

    elif target == 2 and controls[0] == 1:
        #  target -> self, controls[0] -> controls[1] -> outside
        p[0], p[1] = p[1], p[0]
        p[0], p[controls[1]] = p[controls[1]], p[0]

    #
    # N > 4 cases
    #

    elif controls[0] == 1 and controls[1] > 2 and target > 2:
        #  controls[0] -> controls[1] -> outside, target -> outside
        p[0], p[1] = p[1], p[0]
        p[0], p[controls[1]] = p[controls[1]], p[0]
        p[2], p[target] = p[target], p[2]

    elif controls[0] == 2 and controls[1] > 2 and target > 2:
        #  controls[0] -> target -> outside, controls[1] -> outside
        p[0], p[2] = p[2], p[0]
        p[0], p[target] = p[target], p[0]
        p[1], p[controls[1]] = p[controls[1]], p[1]

    elif controls[1] == 2 and controls[0] > 2 and target > 2:
        #  controls[1] -> target -> outside, controls[0] -> outside
        p[1], p[2] = p[2], p[1]
        p[1], p[target] = p[target], p[1]
        p[0], p[controls[0]] = p[controls[0]], p[0]

    else:
        p[0], p[controls[0]] = p[controls[0]], p[0]
        p1[1], p1[controls[1]] = p1[controls[1]], p1[1]
        p2[2], p2[target] = p2[target], p2[2]
        p = [p[p1[p2[k]]] for k in range(N)]

    return tensor([U] + [identity(2)] * (N - 3)).permute(p)