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 testOperatorListState(self):
        """
        expect: operator list and state
        """
        res = expect([sigmax(), sigmay(), sigmaz()], fock(2, 0))
        assert_(len(res) == 3)
        assert_(all(abs(res - [0, 0, 1]) < 1e-12))

        res = expect([sigmax(), sigmay(), sigmaz()], fock_dm(2, 1))
        assert_(len(res) == 3)
        assert_(all(abs(res - [0, 0, -1]) < 1e-12))

    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 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 add_annotation(self, state_or_vector, text, **kwargs):
        """Add a text or LaTeX annotation to Bloch sphere,
        parametrized by a qubit state or a vector.

        Parameters
        ----------
        state_or_vector : Qobj/array/list/tuple
            Position for the annotaion.
            Qobj of a qubit or a vector of 3 elements.

        text : str/unicode
            Annotation text.
            You can use LaTeX, but remember to use raw string
            e.g. r"$\\langle x \\rangle$"
            or escape backslashes
            e.g. "$\\\\langle x \\\\rangle$".

        **kwargs :
            Options as for mplot3d.axes3d.text, including:
            fontsize, color, horizontalalignment, verticalalignment.
        """
        if isinstance(state_or_vector, Qobj):
            vec = [expect(sigmax(), state_or_vector),
                   expect(sigmay(), state_or_vector),
                   expect(sigmaz(), state_or_vector)]
        elif isinstance(state_or_vector, (list, ndarray, tuple)) \
                and len(state_or_vector) == 3:
            vec = state_or_vector
        else:
            raise Exception("Position needs to be specified by a qubit " +
                            "state or a 3D vector.")
        self.annotations.append({'position': vec,
                                 'text': text,
                                 'opts': kwargs})

    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 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 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 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 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 qubit_clifford_group(N=None, target=0):
    """
    Generates the Clifford group on a single qubit,
    using the presentation of the group given by Ross and Selinger
    (http://www.mathstat.dal.ca/~selinger/newsynth/).

    Parameters
    -----------

    N : int or None
        Number of qubits on which each operator is to be defined
        (default: 1).
    target : int
        Index of the target qubit on which the single-qubit
        Clifford operators are to act.

    Yields
    ------

    op : Qobj
        Clifford operators, represented as Qobj instances.
    
    """

    # The Ross-Selinger presentation of the single-qubit Clifford
    # group expresses each element in the form C_{ijk} = E^i X^j S^k
    # for gates E, X and S, and for i in range(3), j in range(2) and
    # k in range(4).
    #
    # We start by defining these gates. E is defined in terms of H,
    # \omega and S, so we define \omega and H first.
    w = np.exp(1j * 2 * np.pi / 8)
    H = snot()

    X = sigmax()
    S = phasegate(np.pi / 2)
    E = H * (S ** 3) * w ** 3


    for op in map(
          # partial(reduce, mul) returns a function that takes products of its argument,
          # by analogy to sum. Note that by analogy, sum can be written partial(reduce, add).
          partial(reduce, mul),
          # product(...) yields the Cartesian product of its arguments. Here, each element is
          # a tuple (E**i, X**j, S**k) such that partial(reduce, mul) acting on the tuple
          # yields E**i * X**j * S**k.
          product(_powers(E, 3), _powers(X, 2), _powers(S, 4))
        ):

        # Finally, we optionally expand the gate.
        if N is not None:
            yield gate_expand_1toN(op, N, target)
        else:
            yield op

 def add_states(self, state, kind="vector"):
     """Add a state vector Qobj to Bloch sphere.
     
     Parameters
     ----------
     state : qobj
         Input state vector.
         
     kind : str {'vector','point'}
         Type of object to plot.
     
     """
     if isinstance(state, Qobj):
         state = [state]
     for st in state:
         if kind == "vector":
             vec = [expect(sigmax(), st), expect(sigmay(), st), expect(sigmaz(), st)]
             self.add_vectors(vec)
         elif kind == "point":
             pnt = [expect(sigmax(), st), expect(sigmay(), st), expect(sigmaz(), st)]
             self.add_points(pnt)

def test_QobjHerm():
    "Qobj Hermicity"
    N = 10
    data = np.random.random(
        (N, N)) + 1j * np.random.random((N, N)) - (0.5 + 0.5j)
    q = Qobj(data)
    assert_equal(q.isherm, False)

    data = data + data.conj().T
    q = Qobj(data)
    assert_(q.isherm)

    q_a = destroy(5)
    assert_(not q_a.isherm)

    q_ad = create(5)
    assert_(not q_ad.isherm)

    # test addition of two nonhermitian operators adding up to a hermitian one
    q_x = q_a + q_ad
    assert_hermicity(q_x, True)

    # test addition of one hermitan and one nonhermitian operator
    q = q_x + q_a
    assert_hermicity(q, False)

    # test addition of two hermitan operators
    q = q_x + q_x
    assert_hermicity(q, True)

    # Test multiplication of two Hermitian operators.
    # This results in a skew-Hermitian operator, so
    # we're checking here that __mul__ doesn't set wrong
    # metadata.
    q = sigmax() * sigmay()
    assert_hermicity(q, False, "Expected iZ = X * Y to be skew-Hermitian.")
    # Similarly, we need to check that -Z = X * iY is correctly
    # identified as Hermitian.
    q = sigmax() * (1j * sigmay())
    assert_hermicity(q, True, "Expected -Z = X * iY to be Hermitian.")

    def testExpectSolverCompatibility(self):
        """
        expect: operator list and state list
        """
        c_ops = [0.0001 * sigmaz()]
        e_ops = [sigmax(), sigmay(), sigmaz(), sigmam(), sigmap()]
        times = np.linspace(0, 10, 100)

        res1 = mesolve(sigmax(), fock(2, 0), times, c_ops, e_ops)
        res2 = mesolve(sigmax(), fock(2, 0), times, c_ops, [])

        e1 = res1.expect
        e2 = expect(e_ops, res2.states)

        assert_(len(e1) == len(e2))

        for n in range(len(e1)):
            assert_(len(e1[n]) == len(e2[n]))
            assert_(isinstance(e1[n], np.ndarray))
            assert_(isinstance(e2[n], np.ndarray))
            assert_(e1[n].dtype == e2[n].dtype)
            assert_(all(abs(e1[n] - e2[n]) < 1e-12))

    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 _spin_hamiltonian(N):
    from qutip.tensor import tensor
    from qutip.operators import qeye, sigmax, sigmay, sigmaz
    # array of spin energy splittings and coupling strengths. here we use
    # uniform parameters, but in general we don't have too
    h  = 1.0 * 2 * np.pi * np.ones(N) 
    Jz = 0.1 * 2 * np.pi * np.ones(N)
    Jx = 0.1 * 2 * np.pi * np.ones(N)
    Jy = 0.1 * 2 * np.pi * np.ones(N)
    # dephasing rate
    gamma = 0.01 * np.ones(N)

    si = qeye(2)
    sx = sigmax()
    sy = sigmay()
    sz = sigmaz()

    sx_list = []
    sy_list = []
    sz_list = []

    for n in range(N):
        op_list = []
        for m in range(N):
            op_list.append(si)

        op_list[n] = sx
        sx_list.append(tensor(op_list))

        op_list[n] = sy
        sy_list.append(tensor(op_list))

        op_list[n] = sz
        sz_list.append(tensor(op_list))

    # construct the hamiltonian
    H = 0

    # energy splitting terms
    for n in range(N):
        H += - 0.5 * h[n] * sz_list[n]

    # interaction terms
    for n in range(N-1):
        H += - 0.5 * Jx[n] * sx_list[n] * sx_list[n+1]
        H += - 0.5 * Jy[n] * sy_list[n] * sy_list[n+1]
        H += - 0.5 * Jz[n] * sz_list[n] * sz_list[n+1]
    return H

def test_QobjUnitaryOper():
    "Qobj unitarity"
    # Check some standard operators
    Sx = sigmax()
    Sy = sigmay()
    assert_unitarity(qeye(4), True, "qeye(4) should be unitary.")
    assert_unitarity(Sx, True, "sigmax() should be unitary.")
    assert_unitarity(Sy, True, "sigmax() should be unitary.")
    assert_unitarity(sigmam(), False, "sigmam() should NOT be unitary.")
    assert_unitarity(destroy(10), False, "destroy(10) should NOT be unitary.")
    # Check multiplcation of unitary is unitary
    assert_unitarity(Sx*Sy, True, "sigmax()*sigmay() should be unitary.")
    # Check some other operations clear unitarity
    assert_unitarity(Sx+Sy, False, "sigmax()+sigmay() should NOT be unitary.")
    assert_unitarity(4*Sx, False, "4*sigmax() should NOT be unitary.")
    assert_unitarity(Sx*4, False, "sigmax()*4 should NOT be unitary.")
    assert_unitarity(4+Sx, False, "4+sigmax() should NOT be unitary.")
    assert_unitarity(Sx+4, False, "sigmax()+4 should NOT be unitary.")

    def testOperatorListStateList(self):
        """
        expect: operator list and state list
        """
        operators = [sigmax(), sigmay(), sigmaz(), sigmam(), sigmap()]
        states = [fock(2, 0), fock(2, 1), fock_dm(2, 0), fock_dm(2, 1)]
        res = expect(operators, states)

        assert_(len(res) == len(operators))

        for r_idx, r in enumerate(res):

            assert_(isinstance(r, np.ndarray))

            if operators[r_idx].isherm:
                assert_(r.dtype == np.float64)
            else:
                assert_(r.dtype == np.complex128)

            for s_idx, s in enumerate(states):
                assert_(r[s_idx] == expect(operators[r_idx], states[s_idx]))

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

    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)