def test_scalar_numpy():
    if not np:
        skip("numpy not installed or Python too old.")

    assert represent(Integer(1), format='numpy') == 1
    assert represent(Float(1.0), format='numpy') == 1.0
    assert represent(1.0+I, format='numpy') == 1.0+1.0j

def test_QubitBra():
    assert Qubit(0).dual_class() == QubitBra
    assert QubitBra(0).dual_class() == Qubit
    assert represent(Qubit(1,1,0), nqubits=3).H ==\
           represent(QubitBra(1,1,0), nqubits=3)
    assert Qubit(0,1)._eval_innerproduct_QubitBra(QubitBra(1,0)) == Integer(0)
    assert Qubit(0,1)._eval_innerproduct_QubitBra(QubitBra(0,1)) == Integer(1)

 def _rewrite_basis(self, basis, evect, **options):
     from sympy.physics.quantum.represent import represent
     j = sympify(self.j)
     jvals = self.jvals
     if j.is_number:
         if j == int(j):
             start = j**2
         else:
             start = (2*j-1)*(2*j+1)/4
         vect = represent(self, basis=basis, **options)
         result = Add(*[vect[start+i] * evect(j,j-i,*jvals) for i in range(2*j+1)])
         if options.get('coupled') is False:
             return uncouple(result)
         return result
     else:
         # TODO: better way to get angles of rotation
         mi = symbols('mi')
         angles = represent(self.__class__(0,mi),basis=basis)[0].args[3:6]
         if angles == (0,0,0):
             return self
         else:
             state = evect(j, mi, *jvals)
             lt = Rotation.D(j, mi, self.m, *angles)
             result = lt * state
             return Sum(lt * state, (mi,-j,j))

def test_RaisingOp():
    assert Dagger(ad) == a
    assert Commutator(ad, a).doit() == Integer(-1)
    assert Commutator(ad, N).doit() == Integer(-1)*ad
    assert qapply(ad*k) == (sqrt(k.n + 1)*SHOKet(k.n + 1)).expand()
    assert qapply(ad*kz) == (sqrt(kz.n + 1)*SHOKet(kz.n + 1)).expand()
    assert qapply(ad*kf) == (sqrt(kf.n + 1)*SHOKet(kf.n + 1)).expand()
    assert ad.rewrite('xp').doit() == \
        (Integer(1)/sqrt(Integer(2)*hbar*m*omega))*(Integer(-1)*I*Px + m*omega*X)
    assert ad.hilbert_space == ComplexSpace(S.Infinity)
    for i in range(ndim - 1):
        assert ad_rep_sympy[i + 1,i] == sqrt(i + 1)

    if not np:
        skip("numpy not installed or Python too old.")

    ad_rep_numpy = represent(ad, basis=N, ndim=4, format='numpy')
    for i in range(ndim - 1):
        assert ad_rep_numpy[i + 1,i] == float(sqrt(i + 1))

    if not np:
        skip("numpy not installed or Python too old.")
    if not scipy:
        skip("scipy not installed.")
    else:
        sparse = scipy.sparse

    ad_rep_scipy = represent(ad, basis=N, ndim=4, format='scipy.sparse', spmatrix='lil')
    for i in range(ndim - 1):
        assert ad_rep_scipy[i + 1,i] == float(sqrt(i + 1))

    assert ad_rep_numpy.dtype == 'float64'
    assert ad_rep_scipy.dtype == 'float64'

def test_entropy():
    up = JzKet(S(1)/2, S(1)/2)
    down = JzKet(S(1)/2, -S(1)/2)
    d = Density((up, 0.5), (down, 0.5))

    # test for density object
    ent = entropy(d)
    assert entropy(d) == 0.5*log(2)
    assert d.entropy() == 0.5*log(2)

    np = import_module('numpy', min_module_version='1.4.0')
    if np:
        #do this test only if 'numpy' is available on test machine
        np_mat = represent(d, format='numpy')
        ent = entropy(np_mat)
        assert isinstance(np_mat, np.matrixlib.defmatrix.matrix)
        assert ent.real == 0.69314718055994529
        assert ent.imag == 0

    scipy = import_module('scipy', __import__kwargs={'fromlist': ['sparse']})
    if scipy and np:
        #do this test only if numpy and scipy are available
        mat = represent(d, format="scipy.sparse")
        assert isinstance(mat, scipy_sparse_matrix)
        assert ent.real == 0.69314718055994529
        assert ent.imag == 0

def fidelity(state1, state2):
    """ Computes the fidelity [1]_ between two quantum states

    The arguments provided to this function should be a square matrix or a
    Density object. If it is a square matrix, it is assumed to be diagonalizable.

    Parameters
    ==========

    state1, state2 : a density matrix or Matrix


    Examples
    ========

    >>> from sympy import S, sqrt
    >>> from sympy.physics.quantum.dagger import Dagger
    >>> from sympy.physics.quantum.spin import JzKet
    >>> from sympy.physics.quantum.density import Density, fidelity
    >>> from sympy.physics.quantum.represent import represent
    >>>
    >>> up = JzKet(S(1)/2,S(1)/2)
    >>> down = JzKet(S(1)/2,-S(1)/2)
    >>> amp = 1/sqrt(2)
    >>> updown = (amp * up) + (amp * down)
    >>>
    >>> # represent turns Kets into matrices
    >>> up_dm = represent(up * Dagger(up))
    >>> down_dm = represent(down * Dagger(down))
    >>> updown_dm = represent(updown * Dagger(updown))
    >>>
    >>> fidelity(up_dm, up_dm)
    1
    >>> fidelity(up_dm, down_dm) #orthogonal states
    0
    >>> fidelity(up_dm, updown_dm).evalf().round(3)
    0.707

    References
    ==========

    .. [1] https://en.wikipedia.org/wiki/Fidelity_of_quantum_states

    """
    state1 = represent(state1) if isinstance(state1, Density) else state1
    state2 = represent(state2) if isinstance(state2, Density) else state2

    if (not isinstance(state1, Matrix) or
            not isinstance(state2, Matrix)):
        raise ValueError("state1 and state2 must be of type Density or Matrix "
                         "received type=%s for state1 and type=%s for state2" %
                         (type(state1), type(state2)))

    if ( state1.shape != state2.shape and state1.is_square):
        raise ValueError("The dimensions of both args should be equal and the "
                         "matrix obtained should be a square matrix")

    sqrt_state1 = state1**Rational(1, 2)
    return Tr((sqrt_state1 * state2 * sqrt_state1)**Rational(1, 2)).doit()

def test_cnot_gate():
    """Test the CNOT gate."""
    circuit = CNotGate(1,0)
    assert represent(circuit, nqubits=2) ==\
        Matrix([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]])
    circuit = circuit*Qubit('111')
    assert matrix_to_qubit(represent(circuit, nqubits=3)) ==\
        apply_operators(circuit)

def test_scalar_scipy_sparse():
    if not np:
        skip("numpy not installed or Python too old.")
    if not scipy:
        skip("scipy not installed.")

    assert represent(Integer(1), format='scipy.sparse') == 1
    assert represent(Float(1.0), format='scipy.sparse') == 1.0
    assert represent(1.0+I, format='scipy.sparse') == 1.0+1.0j

def test_swap_gate():
    """Test the SWAP gate."""
    swap_gate_matrix = Matrix(((1, 0, 0, 0), (0, 0, 1, 0), (0, 1, 0, 0), (0, 0, 0, 1)))
    assert represent(SwapGate(1, 0).decompose(), nqubits=2) == swap_gate_matrix
    assert qapply(SwapGate(1, 3) * Qubit("0010")) == Qubit("1000")
    nqubits = 4
    for i in range(nqubits):
        for j in range(i):
            assert represent(SwapGate(i, j), nqubits=nqubits) == represent(SwapGate(i, j).decompose(), nqubits=nqubits)

def test_swap_gate():
    """Test the SWAP gate."""
    swap_gate_matrix = Matrix(((1,0,0,0),(0,0,1,0),(0,1,0,0),(0,0,0,1)))
    assert represent(SwapGate(1,0).decompose(), nqubits=2) == swap_gate_matrix
    assert apply_operators(SwapGate(1,3)*Qubit('0010')) == Qubit('1000')
    nqubits = 4
    for i in range(nqubits):
        for j in range(i):
            assert represent(SwapGate(i,j), nqubits=nqubits) ==\
                represent(SwapGate(i,j).decompose(), nqubits=nqubits)

def test_one_qubit_commutators():
    """Test single qubit gate commutation relations."""
    for g1 in (IdentityGate, X, Y, Z, H, T, S):
        for g2 in (IdentityGate, X, Y, Z, H, T, S):
            e = Commutator(g1(0),g2(0))
            a = matrix_to_zero(represent(e, nqubits=1, format='sympy'))
            b = matrix_to_zero(represent(e.doit(), nqubits=1, format='sympy'))
            assert a == b
            e = Commutator(g1(0),g2(1))
            assert e.doit() == 0

def test_cnot_gate():
    """Test the CNOT gate."""
    circuit = CNotGate(1, 0)
    assert represent(circuit, nqubits=2) == Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
    circuit = circuit * Qubit("111")
    assert matrix_to_qubit(represent(circuit, nqubits=3)) == qapply(circuit)

    circuit = CNotGate(1, 0)
    assert Dagger(circuit) == circuit
    assert Dagger(Dagger(circuit)) == circuit
    assert circuit * circuit == 1

def test_one_qubit_anticommutators():
    """Test single qubit gate anticommutation relations."""
    for g1 in (IdentityGate, X, Y, Z, H):
        for g2 in (IdentityGate, X, Y, Z, H):
            e = AntiCommutator(g1(0), g2(0))
            a = matrix_to_zero(represent(e, nqubits=1, format="sympy"))
            b = matrix_to_zero(represent(e.doit(), nqubits=1, format="sympy"))
            assert a == b
            e = AntiCommutator(g1(0), g2(1))
            a = matrix_to_zero(represent(e, nqubits=2, format="sympy"))
            b = matrix_to_zero(represent(e.doit(), nqubits=2, format="sympy"))
            assert a == b

def bures_angle(state1, state2):
    """ Computes the Bures angle [1], [2] between two quantum states
    The arguments provided to this function should be a square matrix or a
    Density object. If it is a square matrix, it is assumed to be diagonalizable.
    Parameters:
    ==========
    state1, state2 : a density matrix or Matrix
    Examples:
    =========
    >>> from sympy.physics.quantum import TensorProduct, Ket, Dagger
    >>> from sympy.physics.quantum.density import bures_angle
    >>> from sympy import Matrix
    >>> from math import sqrt
    >>> # define qubits |0>, |1>, |00>, and |11>
    >>> q0 = Matrix([1,0])
    >>> q1 = Matrix([0,1])
    >>> q00 = TensorProduct(q0,q0)
    >>> q11 = TensorProduct(q1,q1)
    >>> # create set of maximally entangled Bell states 
    >>> phip = 1/sqrt(2) * ( q00 + q11 )
    >>> phim = 1/sqrt(2) * ( q00 - q11 )
    >>> # create the corresponding density matrices for the Bell states
    >>> phip_dm = phip * Dagger(phip)
    >>> phim_dm = phim * Dagger(phim)
    >>> # calculates the Bures angle between two orthogonal states (yields: 0)
    >>> print bures_angle(phip_dm, phim_dm)
    1.0
    >>> # calculates the Bures angle between two identitcal states (yields: 1)
    >>> print bures_angle(phip_dm, phip_dm)
    0.0

    References
    ==========
    .. [1] http://en.wikipedia.org/wiki/Bures_metric
    .. [2] Quantum Computation and Quantum Information. M. Nielsen, I. Chuang, 
           Cambridge University Press, (2001) (Eq. 9.82, pg 413). 
    """
    state1 = represent(state1) if isinstance(state1, Density) else state1
    state2 = represent(state2) if isinstance(state2, Density) else state2

    if (not isinstance(state1, Matrix) or
            not isinstance(state2, Matrix)):
        raise ValueError("state1 and state2 must be of type Density or Matrix "
                         "received type=%s for state1 and type=%s for state2" %
                         (type(state1), type(state2)))

    if ( state1.shape != state2.shape and state1.is_square):
        raise ValueError("The dimensions of both args should be equal and the "
                         "matrix obtained should be a square matrix")

    # the pi/2 is for normalization
    return acos( fidelity(state1, state2) ) / (pi/2)

def test_UGate_OneQubitGate_combo():
    v, w, f, g = symbols('v w f g')
    uMat1 = ImmutableMatrix([[v, w], [f, g]])
    cMat1 = Matrix([[v, w + 1, 0, 0], [f + 1, g, 0, 0], [0, 0, v, w + 1], [0, 0, f + 1, g]])
    u1 = X(0) + UGate(0, uMat1)
    assert represent(u1, nqubits=2) == cMat1

    uMat2 = ImmutableMatrix([[1/sqrt(2), 1/sqrt(2)], [I/sqrt(2), -I/sqrt(2)]])
    cMat2_1 = Matrix([[1/2 + I/2, 1/2 - I/2], [1/2 - I/2, 1/2 + I/2]])
    cMat2_2 = Matrix([[1, 0], [0, I]])
    u2 = UGate(0, uMat2)
    assert represent(H(0)*u2, nqubits=1) == cMat2_1
    assert represent(u2*H(0), nqubits=1) == cMat2_2

def test_UGate():
    a, b, c, d = symbols("a,b,c,d")
    uMat = Matrix([[a, b], [c, d]])

    # Test basic case where gate exists in 1-qubit space
    u1 = UGate((0,), uMat)
    assert represent(u1, nqubits=1) == uMat
    assert qapply(u1 * Qubit("0")) == a * Qubit("0") + c * Qubit("1")
    assert qapply(u1 * Qubit("1")) == b * Qubit("0") + d * Qubit("1")

    # Test case where gate exists in a larger space
    u2 = UGate((1,), uMat)
    u2Rep = represent(u2, nqubits=2)
    for i in range(4):
        assert u2Rep * qubit_to_matrix(IntQubit(i, 2)) == qubit_to_matrix(qapply(u2 * IntQubit(i, 2)))

def test_OracleGate():
    v = OracleGate(1, lambda qubits: qubits == IntQubit(0))
    assert qapply(v*IntQubit(0)) == -IntQubit(0)
    assert qapply(v*IntQubit(1)) == IntQubit(1)

    nbits = 2
    v = OracleGate(2, return_one_on_two)
    assert qapply(v*IntQubit(0, nbits)) == IntQubit(0, nqubits=nbits)
    assert qapply(v*IntQubit(1, nbits)) == IntQubit(1, nqubits=nbits)
    assert qapply(v*IntQubit(2, nbits)) == -IntQubit(2, nbits)
    assert qapply(v*IntQubit(3, nbits)) == IntQubit(3, nbits)

    assert represent(OracleGate(1, lambda qubits: qubits == IntQubit(0)), nqubits=1) == \
           Matrix([[-1, 0], [0, 1]])
    assert represent(v, nqubits=2) == Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]])

def test_superposition_of_states():
    assert qapply(CNOT(0, 1)*HadamardGate(0)*(1/sqrt(2)*Qubit('01') + 1/sqrt(2)*Qubit('10'))).expand() == (Qubit('01')/2 + Qubit('00')/2 - Qubit('11')/2 +
     Qubit('10')/2)

    assert matrix_to_qubit(represent(CNOT(0, 1)*HadamardGate(0)
        *(1/sqrt(2)*Qubit('01') + 1/sqrt(2)*Qubit('10')), nqubits=2)) == \
        (Qubit('01')/2 + Qubit('00')/2 - Qubit('11')/2 + Qubit('10')/2)

def qubit_to_matrix(qubit, format='sympy'):
    """Converts an Add/Mul of Qubit objects into it's matrix representation

    This function is the inverse of ``matrix_to_qubit`` and is a shorthand
    for ``represent(qubit)``.
    """
    return represent(qubit, format=format)

def test_apply_represent_equality():
    gates = [
        HadamardGate(int(3 * random.random())),
        XGate(int(3 * random.random())),
        ZGate(int(3 * random.random())),
        YGate(int(3 * random.random())),
        ZGate(int(3 * random.random())),
        PhaseGate(int(3 * random.random())),
    ]

    circuit = Qubit(
        int(random.random() * 2),
        int(random.random() * 2),
        int(random.random() * 2),
        int(random.random() * 2),
        int(random.random() * 2),
        int(random.random() * 2),
    )
    for i in range(int(random.random() * 6)):
        circuit = gates[int(random.random() * 6)] * circuit

    mat = represent(circuit, nqubits=6)
    states = qapply(circuit)
    state_rep = matrix_to_qubit(mat)
    states = states.expand()
    state_rep = state_rep.expand()
    assert state_rep == states