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 snot(N=None, target=None):
    """Quantum object representing the SNOT (2-qubit Hadamard) gate.

    Returns
    -------
    snot_gate : qobj
        Quantum object representation of SNOT gate.

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

    """
    if not N is None and not target is None:
        return gate_expand_1toN(snot(), N, target)
    else:
        u = basis(2, 0)
        d = basis(2, 1)
        Q = 1.0 / sqrt(2.0) * (sigmax() + sigmaz())
        return Q

def phasegate(theta):
    """
    Returns quantum object representing the phase shift gate.

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

    Returns
    -------
    phase_gate : qobj
        Quantum object representation of phase shift gate.

    Examples
    --------
    >>> phasegate(pi/4)
    Quantum object: dims = [[2], [2]], \
shape = [2, 2], type = oper, isHerm = False
    Qobj data =
    [[ 1.00000000+0.j          0.00000000+0.j        ]
     [ 0.00000000+0.j          0.70710678+0.70710678j]]

    """
    u = basis(2)
    d = basis(2,1)
    Q = u * u.dag() + (exp(1.0j * theta) * d * d.dag())
    return Qobj(Q)

def qubit_states(N=1, states=[0]):
    """
    Function to define initial state of the qubits.

    Parameters
    ----------
    N : Integer
        Number of qubits in the register.
    states : List
        Initial state of each qubit.

    Returns
    ----------
    qstates : Qobj
        List of qubits.

    """
    state_list = []
    for i in range(N):
        if N > len(states) and i >= len(states):
            state_list.append(0)
        else:
            state_list.append(states[i])

    return tensor([alpha * basis(2, 1) + sqrt(1 - alpha**2) * basis(2, 0)
                  for alpha in state_list])

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 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 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 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_unit():
    """
    Test Qobj: unit
    """
    psi = 10*np.random.randn()*basis(2,0)-10*np.random.randn()*1j*basis(2,1)
    psi2 = psi.unit()
    psi.unit(inplace=True)
    assert_(psi == psi2)
    assert_almost_equal(np.linalg.norm(psi.full()), 1.0)

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 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 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 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()

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 testQubitStates(self):
        """
        Tests the qubit_states function.
        """
        psi0_a = basis(2, 0)
        psi0_b = qubit_states()
        assert_(psi0_a == psi0_b)

        psi1_a = basis(2, 1)
        psi1_b = qubit_states(states=[1])
        assert_(psi1_a == psi1_b)

        psi01_a = tensor(psi0_a, psi1_a)
        psi01_b = qubit_states(N=2, states=[0, 1])
        assert_(psi01_a == psi01_b)

def test_KetType():
    "Qobj ket type"

    psi = basis(2, 1)

    assert_(psi.isket)
    assert_(not psi.isbra)
    assert_(not psi.isoper)
    assert_(not psi.issuper)

    psi = tensor(basis(2, 1), basis(2, 0))

    assert_(psi.isket)
    assert_(not psi.isbra)
    assert_(not psi.isoper)
    assert_(not psi.issuper)

def test_BraType():
    "Qobj bra type"

    psi = basis(2, 1).dag()

    assert_equal(psi.isket, False)
    assert_equal(psi.isbra, True)
    assert_equal(psi.isoper, False)
    assert_equal(psi.issuper, False)

    psi = tensor(basis(2, 1).dag(), basis(2, 0).dag())

    assert_equal(psi.isket, False)
    assert_equal(psi.isbra, True)
    assert_equal(psi.isoper, False)
    assert_equal(psi.issuper, False)

def rand_ket_haar(N=2, dims=None):
    """
    Returns a Haar random pure state of dimension ``dim`` by
    applying a Haar random unitary to a fixed pure state.


    Parameters
    ----------
    N : int
        Dimension of the state vector to be returned.

    dims : list of ints, or None
        Left-dimensions of the resultant quantum object.
        If None, [N] is used.

    Returns
    -------
    psi : Qobj
        A random state vector drawn from the Haar measure.
    """
    if dims:
        _check_ket_dims(dims, N)
    else:
        dims = [N]
    psi = rand_unitary_haar(N) * basis(N, 0)
    psi.dims = [dims, [1]]
    return psi

    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))