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 cnot():
    """
    Quantum object representing the CNOT gate.

    Returns
    -------
    cnot_gate : qobj
        Quantum object representation of CNOT gate

    Examples
    --------
    >>> cnot()
    Quantum object: dims = [[2, 2], [2, 2]], \
shape = [4, 4], type = oper, isHerm = True
    Qobj data =
        [[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
         [ 0.+0.j  1.+0.j  0.+0.j  0.+0.j]
         [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]
         [ 0.+0.j  0.+0.j  1.+0.j  0.+0.j]]

    """
    uu = tensor(basis(2),basis(2))
    ud = tensor(basis(2),basis(2,1))
    du = tensor(basis(2,1),basis(2))
    dd = tensor(basis(2,1),basis(2,1))
    Q = uu * uu.dag() + ud * ud.dag() + dd * du.dag() + du * dd.dag()
    return Qobj(Q)

def _subsystem_apply_reference(state, channel, mask):
    if isket(state):
        state = ket2dm(state)

    if isoper(channel):
        full_oper = tensor([channel if mask[j]
                            else qeye(state.dims[0][j])
                            for j in range(len(state.dims[0]))])
        return full_oper * state * full_oper.dag()
    else:
        # Go to Choi, then Kraus
        # chan_mat = array(channel.data.todense())
        choi_matrix = super_to_choi(channel)
        vals, vecs = eig(choi_matrix.full())
        vecs = list(map(array, zip(*vecs)))
        kraus_list = [sqrt(vals[j]) * vec2mat(vecs[j])
                      for j in range(len(vals))]
        # Kraus operators to be padded with identities:
        k_qubit_kraus_list = product(kraus_list, repeat=sum(mask))
        rho_out = Qobj(inpt=zeros(state.shape), dims=state.dims)
        for operator_iter in k_qubit_kraus_list:
            operator_iter = iter(operator_iter)
            op_iter_list = [next(operator_iter).conj().T if mask[j]
                            else qeye(state.dims[0][j])
                            for j in range(len(state.dims[0]))]
            full_oper = tensor(list(map(Qobj, op_iter_list)))
            rho_out = rho_out + full_oper * state * full_oper.dag()
        return Qobj(rho_out)

    def testExpandGate2toNSwap(self):
        """
        gates: expand 2 to N (using swap)
        """

        a, b = np.random.rand(), np.random.rand()
        k1 = (a * basis(2, 0) + b * basis(2, 1)).unit()

        c, d = np.random.rand(), np.random.rand()
        k2 = (c * basis(2, 0) + d * basis(2, 1)).unit()

        N = 6
        kets = [rand_ket(2) for k in range(N)]

        for m in range(N):
            for n in set(range(N)) - {m}:

                psi_in = tensor([k1 if k == m else k2 if k == n else kets[k]
                                 for k in range(N)])
                psi_out = tensor([k2 if k == m else k1 if k == n else kets[k]
                                  for k in range(N)])

                targets = [m, n]
                G = swap(N, targets)

                psi_out = G * psi_in

                assert_((psi_out - G * psi_in).norm() < 1e-12)

def entangling_power(U):
    """
    Calculate the entangling power of a two-qubit gate U, which
    is zero of nonentangling gates and 1 and 2/9 for maximally
    entangling gates.

    Parameters
    ----------
    U : qobj
        Qobj instance representing a two-qubit gate.

    Returns
    -------
    ep : float
        The entanglement power of U (real number between 0 and 1)

    References:

        Explorations in Quantum Computing, Colin P. Williams (Springer, 2011)
    """

    if not U.isoper:
        raise Exception("U must be an operator.")

    if U.dims != [[2, 2], [2, 2]]:
        raise Exception("U must be a two-qubit gate.")

    a = (tensor(U, U).dag() * swap(N=4, targets=[1, 3]) *
         tensor(U, U) * swap(N=4, targets=[1, 3]))
    b = (tensor(swap() * U, swap() * U).dag() * swap(N=4, targets=[1, 3]) *
         tensor(swap() * U, swap() * U) * swap(N=4, targets=[1, 3]))

    return 5.0/9 - 1.0/36 * (a.tr() + b.tr()).real

def _jc_liouvillian(N):
    from qutip.tensor import tensor
    from qutip.operators import destroy, qeye
    from qutip.superoperator import liouvillian
    wc = 1.0  * 2 * np.pi  # cavity frequency
    wa = 1.0  * 2 * np.pi  # atom frequency
    g  = 0.05 * 2 * np.pi  # coupling strength
    kappa = 0.005          # cavity dissipation rate
    gamma = 0.05           # atom dissipation rate
    n_th_a = 1           # temperature in frequency units
    use_rwa = 0
    # operators
    a  = tensor(destroy(N), qeye(2))
    sm = tensor(qeye(N), destroy(2))
    # Hamiltonian
    if use_rwa:
        H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
    else:
        H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() + a) * (sm + sm.dag())
    c_op_list = []

    rate = kappa * (1 + n_th_a)
    if rate > 0.0:
        c_op_list.append(np.sqrt(rate) * a)

    rate = kappa * n_th_a
    if rate > 0.0:
        c_op_list.append(np.sqrt(rate) * a.dag())

    rate = gamma
    if rate > 0.0:
        c_op_list.append(np.sqrt(rate) * sm)

    return liouvillian(H, c_op_list)

def choi_to_stinespring(q_oper, thresh=1e-10):
    # TODO: document!
    kU, kV = _generalized_kraus(q_oper, thresh=thresh)

    assert(len(kU) == len(kV))
    dK = len(kU)
    dL = kU[0].shape[0]
    dR = kV[0].shape[1]
    # Also remember the dims breakout.
    out_dims, in_dims = q_oper.dims
    out_left, out_right = out_dims
    in_left, in_right = in_dims

    A = Qobj(zeros((dK * dL, dL)), dims=[out_left + [dK], out_right + [1]])
    B = Qobj(zeros((dK * dR, dR)), dims=[in_left + [dK], in_right + [1]])

    for idx_kraus, (KL, KR) in enumerate(zip(kU, kV)):
        A += tensor(KL, basis(dK, idx_kraus))
        B += tensor(KR, basis(dK, idx_kraus))
        
    # There is no input (right) Kraus index, so strip that off.
    del A.dims[1][-1]
    del B.dims[1][-1]

    return A, B

    def test_ComplexSuperApply(self):
        """
        Superoperator: Efficient numerics and reference return same result,
        acting on non-composite system
        """
        rho_list = list(map(rand_dm, [2, 3, 2, 3, 2]))
        rho_input = tensor(rho_list)
        superop = kraus_to_super(rand_kraus_map(3))

        analytic_result = rho_list
        analytic_result[1] = Qobj(vec2mat(superop.data.todense() *
                                  mat2vec(analytic_result[1].data.todense())))
        analytic_result[3] = Qobj(vec2mat(superop.data.todense() *
                                  mat2vec(analytic_result[3].data.todense())))
        analytic_result = tensor(analytic_result)

        naive_result = subsystem_apply(rho_input, superop,
                                       [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, superop,
                                           [False, True, False, True, False])
        efficient_diff = (efficient_result - analytic_result).data.todense()
        assert_(norm(efficient_diff) < 1e-12)

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_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 testExpandGate1toN(self):
        """
        gates: expand 1 to N
        """
        N = 7

        for g in [rx, ry, rz, phasegate]:

            theta = np.random.rand() * 2 * 3.1415
            a, b = np.random.rand(), np.random.rand()
            psi1 = (a * basis(2, 0) + b * basis(2, 1)).unit()
            psi2 = g(theta) * psi1

            psi_rand_list = [rand_ket(2) for k in range(N)]

            for m in range(N):

                psi_in = tensor([psi1 if k == m else psi_rand_list[k]
                                 for k in range(N)])
                psi_out = tensor([psi2 if k == m else psi_rand_list[k]
                                  for k in range(N)])

                G = g(theta, N, m)
                psi_res = G * psi_in

                assert_((psi_out - psi_res).norm() < 1e-12)

def cnot(N=None, control=None, target=None):
    """
    Quantum object representing the CNOT gate.

    Returns
    -------
    cnot_gate : qobj
        Quantum object representation of CNOT gate

    Examples
    --------
    >>> cnot()
    Quantum object: dims = [[2, 2], [2, 2]], \
shape = [4, 4], type = oper, isHerm = True
    Qobj data =
        [[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
         [ 0.+0.j  1.+0.j  0.+0.j  0.+0.j]
         [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]
         [ 0.+0.j  0.+0.j  1.+0.j  0.+0.j]]

    """
    if not N is None and not control is None and not target is None:
        return gate_expand_2toN(cnot(), N, control, target)
    else:
        uu = tensor(basis(2), basis(2))
        ud = tensor(basis(2), basis(2, 1))
        du = tensor(basis(2, 1), basis(2))
        dd = tensor(basis(2, 1), basis(2, 1))
        Q = uu * uu.dag() + ud * ud.dag() + dd * du.dag() + du * dd.dag()
        return Qobj(Q)

 def __init__(self,N,state=None):
     #construct register intial state
     if state==None:
         reg=tensor([basis(2) for k in range(N)])
     if isinstance(state,str):
         state=_reg_str2array(state,N)
         reg=tensor([basis(2,state[k]) for k in state])
     Qobj.__init__(self, reg.data, reg.dims, reg.shape,
              reg.type, reg.isherm, fast=False)

    def __init__(self, N, state=None):
        # construct register intial state
        if state is None:
            reg = tensor([basis(2) for k in range(N)])
        elif isinstance(state, str):
            state = _reg_str2array(state, N)
            reg = tensor([basis(2, k) for k in state])
        else:
            raise ValueError("state must be a string")

        Qobj.__init__(self, reg.data, reg.dims, reg.isherm)

    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 testExpandGate3toN(self):
        """
        gates: expand 3 to N (using toffoli, fredkin, and random 3 qubit gate)
        """

        a, b = np.random.rand(), np.random.rand()
        psi1 = (a * basis(2, 0) + b * basis(2, 1)).unit()

        c, d = np.random.rand(), np.random.rand()
        psi2 = (c * basis(2, 0) + d * basis(2, 1)).unit()

        e, f = np.random.rand(), np.random.rand()
        psi3 = (e * basis(2, 0) + f * basis(2, 1)).unit()

        N = 4
        psi_rand_list = [rand_ket(2) for k in range(N)]

        _rand_gate_U = tensor([rand_herm(2, density=1) for k in range(3)])

        def _rand_3qubit_gate(N=None, controls=None, k=None):
            if N is None:
                return _rand_gate_U
            else:
                return gate_expand_3toN(_rand_gate_U, N, controls, k)

        for g in [fredkin, toffoli, _rand_3qubit_gate]:

            psi_ref_in = tensor(psi1, psi2, psi3)
            psi_ref_out = g() * psi_ref_in

            for m in range(N):
                for n in set(range(N)) - {m}:
                    for k in set(range(N)) - {m, n}:

                        psi_list = [psi_rand_list[p] for p in range(N)]
                        psi_list[m] = psi1
                        psi_list[n] = psi2
                        psi_list[k] = psi3
                        psi_in = tensor(psi_list)

                        if g == fredkin:
                            targets = [n, k]
                            G = g(N, control=m, targets=targets)
                        else:
                            controls = [m, n]
                            G = g(N, controls, k)

                        psi_out = G * psi_in

                        o1 = psi_out.overlap(psi_in)
                        o2 = psi_ref_out.overlap(psi_ref_in)
                        assert_(abs(o1 - o2) < 1e-12)

def bell_state(state='00'):
    """
    Returns the Bell state:

        |B00> = 1 / sqrt(2)*[|0>|0>+|1>|1>]
        |B01> = 1 / sqrt(2)*[|0>|0>-|1>|1>]
        |B10> = 1 / sqrt(2)*[|0>|1>+|1>|0>]
        |B11> = 1 / sqrt(2)*[|0>|1>-|1>|0>]

    Returns
    -------
    Bell_state : qobj
        Bell state

    """
    if state == '00':
        Bell_state = tensor(
            basis(2), basis(2))+tensor(basis(2, 1), basis(2, 1))
    elif state == '01':
        Bell_state = tensor(
            basis(2), basis(2))-tensor(basis(2, 1), basis(2, 1))
    elif state == '10':
        Bell_state = tensor(
            basis(2), basis(2, 1))+tensor(basis(2, 1), basis(2))
    elif state == '11':
        Bell_state = tensor(
            basis(2), basis(2, 1))-tensor(basis(2, 1), basis(2))

    return Bell_state.unit()

def bell_state11():
    """
    Returns the B11 Bell state:
    
        |B11>=1/sqrt(2)*[|0>|1>-|1>|0>]
    
    Returns
    -------
    B11 : qobj
        B11 Bell state
    
    """
    B=tensor(basis(2),basis(2,1))-tensor(basis(2,1),basis(2))
    return B.unit()

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

    A = fock_dm(5, 0)
    B = fock_dm(5, 4)
    C = fock_dm(5, 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(5)
        B = rand_ket(5)
        C = rand_ket(5)
        psi = tensor(A, B, C)
        psi2 = psi.permute([1, 0, 2])
        assert_(psi2 == tensor(B, A, C))

    for ii in range(3):
        A = rand_dm(5)
        B = rand_dm(5)
        C = rand_dm(5)
        rho = tensor(A, B, C)
        rho2 = rho.permute([1, 0, 2])
        assert_(rho2 == tensor(B, A, C))

def bell_state10():
    """
    Returns the B10 Bell state:
    
        |B10>=1/sqrt(2)*[|0>|1>+|1>|0>]
    
    Returns
    -------
    B10 : qobj
        B10 Bell state
    
    """
    B=tensor(basis(2),basis(2,1))+tensor(basis(2,1),basis(2))
    return B.unit()