def test_diagHamiltonian2():
    """
    Diagonalization of composite systems
    """

    H1 = scipy.rand() * sigmax() + scipy.rand() * sigmay() +\
        scipy.rand() * sigmaz()
    H2 = scipy.rand() * sigmax() + scipy.rand() * sigmay() +\
        scipy.rand() * sigmaz()

    H = tensor(H1, H2)

    evals, ekets = H.eigenstates()

    for n in range(len(evals)):
        # assert that max(H * ket - e * ket) is small
        assert_equal(amax(
            abs((H * ekets[n] - evals[n] * ekets[n]).full())) < 1e-10, True)

    N1 = 10
    N2 = 2

    a1 = tensor(destroy(N1), qeye(N2))
    a2 = tensor(qeye(N1), destroy(N2))
    H = scipy.rand() * a1.dag() * a1 + scipy.rand() * a2.dag() * a2 + \
        scipy.rand() * (a1 + a1.dag()) * (a2 + a2.dag())
    evals, ekets = H.eigenstates()

    for n in range(len(evals)):
        # assert that max(H * ket - e * ket) is small
        assert_equal(amax(
            abs((H * ekets[n] - evals[n] * ekets[n]).full())) < 1e-10, True)

    def test_02_2_qft_bounds(self):
        """
        control.pulseoptim: QFT gate with linear initial pulses (bounds)
        assert that amplitudes remain in bounds
        """
        Sx = sigmax()
        Sy = sigmay()
        Sz = sigmaz()
        Si = 0.5*identity(2)

        H_d = 0.5*(tensor(Sx, Sx) + tensor(Sy, Sy) + tensor(Sz, Sz))
        H_c = [tensor(Sx, Si), tensor(Sy, Si), tensor(Si, Sx), tensor(Si, Sy)]
        U_0 = identity(4)
        # Target for the gate evolution - Quantum Fourier Transform gate
        U_targ = qft.qft(2)

        n_ts = 10
        evo_time = 10

        result = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-9,
                        amp_lbound=-1.0, amp_ubound=1.0,
                        init_pulse_type='LIN',
                        gen_stats=True)
        assert_((result.final_amps >= -1.0).all() and
                    (result.final_amps <= 1.0).all(),
                    msg="Amplitude bounds exceeded for QFT")

    def construct_hamiltonian(self, number_of_spins, alpha):
        """
following example
http://qutip.googlecode.com/svn/doc/2.0.0/html/examples/me/ex-24.html
returns H0 - hamiltonian without the B field
and y_list - list of sigma_y operators
"""
        N = number_of_spins
        si = qeye(2)
        sx = sigmax()
        sy = sigmay()
        # constructing a list of operators sx_list and sy_list where
        # the operator sx_list[i] applies sigma_x on the ith particle and
        # identity to the rest
        sx_list = []
        sy_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))
        # construct the hamiltonian
        H0 = 0
        # ising coupling term, time independent
        for i in range(N):
            for j in range(N):
                if i < j:
                    H0 -= abs(i - j) ** -alpha * sx_list[i] * sx_list[j]
        H1 = 0
        for i in range(N):
            H1 -= sy_list[i]
        return H0, H1

def pauli():
	'''Return the Pauli spin matrices for S=1/2'''
	identity = qutip.qeye(2)
	sx = qutip.sigmax()/2
	sy = qutip.sigmay()/2
	sz = qutip.sigmaz()/2
	return identity, sx, sy, sz

def _qubit_integrate(tlist, psi0, epsilon, delta, g1, g2, solver):

    H = epsilon / 2.0 * sigmaz() + delta / 2.0 * sigmax()

    c_op_list = []

    rate = g1
    if rate > 0.0:
        c_op_list.append(np.sqrt(rate) * sigmam())

    rate = g2
    if rate > 0.0:
        c_op_list.append(np.sqrt(rate) * sigmaz())

    e_ops = [sigmax(), sigmay(), sigmaz()]

    if solver == "me":
        output = mesolve(H, psi0, tlist, c_op_list, e_ops)
    elif solver == "es":
        output = essolve(H, psi0, tlist, c_op_list, e_ops)
    elif solver == "mc":
        output = mcsolve(H, psi0, tlist, c_op_list, e_ops, ntraj=750)
    else:
        raise ValueError("unknown solver")

    return output.expect[0], output.expect[1], output.expect[2]

    def __init__(self, N_field_levels, coupling=None, N_qubits=1):

        # basic parameters
        self.N_field_levels = N_field_levels
        self.N_qubits = N_qubits

        if coupling is None:
            self.g = 0
        else:
            self.g = coupling

        # bare operators
        self.idcavity = qt.qeye(self.N_field_levels)
        self.idqubit = qt.qeye(2)
        self.a_bare = qt.destroy(self.N_field_levels)
        self.sm_bare = qt.sigmam()
        self.sz_bare = qt.sigmaz()
        self.sx_bare = qt.sigmax()
        self.sy_bare = qt.sigmay()

        # 1 atom 1 cavity operators
        self.jc_a = qt.tensor(self.a_bare, self.idqubit)
        self.jc_sm = qt.tensor(self.idcavity, self.sm_bare)
        self.jc_sx = qt.tensor(self.idcavity, self.sx_bare)
        self.jc_sy = qt.tensor(self.idcavity, self.sy_bare)
        self.jc_sz = qt.tensor(self.idcavity, self.sz_bare)

 def setUp(self):
     TestRotatingFrame.setUp(self)
     self.times_to_calc = np.array([0, np.pi/2])
     self.expected_result =np.array(
         np.cos(self.frequency * self.times_to_calc) * q.qeye(2) + \
          1j* np.sin(self.frequency * self.times_to_calc) * q.sigmay()
     )

 def construct_hamiltonian(self, number_of_spins, alpha, B):
     '''
     following example
     http://qutip.googlecode.com/svn/doc/2.0.0/html/examples/me/ex-24.html
     '''
     N = number_of_spins
     si = qeye(2)
     sx = sigmax()
     sy = sigmay()
     #constructing a list of operators sx_list and sy_list where
     #the operator sx_list[i] applies sigma_x on the ith particle and 
     #identity to the rest
     sx_list = []
     sy_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))
     #construct the hamiltonian
     H = 0
     #magnetic field term, hamiltonian is in units of J0
     for i in range(N):
         H-= B * sy_list[i]
     #ising coupling term
     for i in range(N):
         for j in range(N):
             if i < j:
                 H+= abs(i - j)**-alpha * sx_list[i] * sx_list[j]
     return H

def pauli():
    '''Define pauli spin matrices'''
    identity = qutip.qeye(2)
    sx = qutip.sigmax()/2
    sy = qutip.sigmay()/2
    sz = qutip.sigmaz()/2
    return identity, sx, sy, sz

def test_Transformation1():
    "Transform 2-level to eigenbasis and back"
    H1 = scipy.rand() * sigmax() + scipy.rand() * sigmay() + \
        scipy.rand() * sigmaz()
    evals, ekets = H1.eigenstates()
    Heb = H1.transform(ekets)        # eigenbasis (should be diagonal)
    H2 = Heb.transform(ekets, True)  # back to original basis
    assert_equal((H1 - H2).norm() < 1e-6, True)

def test_Transformation6():
    "Check diagonalization via eigenbasis transformation"

    cx, cy, cz = np.random.rand(), np.random.rand(), np.random.rand()
    H = cx * sigmax() + cy * sigmay() + cz * sigmaz()
    evals, evecs = H.eigenstates()
    Heb = H.transform(evecs).tidyup()  # Heb should be diagonal
    assert_(abs(Heb.full() - np.diag(Heb.full().diagonal())).max() < 1e-6)

def Hfred(x,N) : # A possible Hamiltonian for the Fredkin gate
    k = 0
    H = 0
    sx = qt.sigmax()/2
    sy = qt.sigmay()/2
    sz = qt.sigmaz()/2
    Id = qt.qeye(2)

    for q in [sx, sy, sz] :
        temp = 0
        OpChain = [Id]*N
        OpChain[2] = q
        OpChain[1] = q
        temp += x[k]*qt.tensor(OpChain)
        H += temp 
    k+=1        

    for q in [sx,sz]:
        
        temp = 0
        OpChain = [Id]*N
        OpChain[2] = q
        OpChain[0] = q
        temp += x[k]*qt.tensor(OpChain)
    
        OpChain = [Id]*N
        OpChain[1] = q
        OpChain[0] = q
        temp += x[k]*qt.tensor(OpChain)
        k += 1
        H += temp 
    
    for q in [1,2]:
        temp = 0
        OpChain = [Id]*N
        OpChain[q] = sx
        OpChain[3] = sx
        temp += x[k]*qt.tensor(OpChain)
        H += temp 
    k+=1        

    temp = 0
    OpChain = [Id]*N
    OpChain[0] = sz
    temp += x[k]*qt.tensor(OpChain)
    H += temp 
    k += 1
    
    temp = 0
    OpChain = [Id]*N
    OpChain[3] = sx
    temp += x[k]*qt.tensor(OpChain)#last one

    H += temp 

    
    return H

    def local_hamiltonian(self):
        field_vector = self.parent_field.field_vector.values()
        pauli_basis = [qt.sigmax(), qt.sigmay(), qt.sigmaz()]

        b_field = sum(
            [field_vector[i] * pauli_basis[i]
             for i in xrange(0, len(field_vector))
             ]
        )
        return self.gyromagnetic_ratio * b_field * const.HBAR / 2

def test_diagHamiltonian1():
    """
    Diagonalization of random two-level system
    """

    H = scipy.rand() * sigmax() + scipy.rand() * sigmay() +\
        scipy.rand() * sigmaz()

    evals, ekets = H.eigenstates()

    for n in range(len(evals)):
        # assert that max(H * ket - e * ket) is small
        assert_equal(amax(
            abs((H * ekets[n] - evals[n] * ekets[n]).full())) < 1e-10, True)

def test_qpt_snot():
    "quantum process tomography for snot gate"

    U_psi = snot()
    U_rho = spre(U_psi) * spost(U_psi.dag())
    N = 1
    op_basis = [[qeye(2), sigmax(), 1j * sigmay(), sigmaz()] for i in range(N)]
    # op_label = [["i", "x", "y", "z"] for i in range(N)]
    chi1 = qpt(U_rho, op_basis)

    chi2 = np.zeros((2 ** (2 * N), 2 ** (2 * N)), dtype=complex)
    chi2[1, 1] = chi2[1, 3] = chi2[3, 1] = chi2[3, 3] = 0.5

    assert_(norm(chi2 - chi1) < 1e-8)

 def as_qobj(self):
     """
     Returns a representation of the given Pauli operator as a QuTiP
     Qobj instance.
     """
     if qt is None:
         raise RuntimeError("Requires QuTiP.")
     if self == Pauli.I:
         return qt.qeye(2)
     elif self == Pauli.X:
         return qt.sigmax()
     elif self == Pauli.Y:
         return qt.sigmay()
     else:
         return qt.sigmaz()

def ptracetest():
    gamma = 1.
    neq = 2
    psi0 = qt.basis(neq,neq-1)
    psi0 = qt.tensor(psi0,psi0)
    H = qt.tensor(qt.sigmax(),qt.sigmay())
    c1 = np.sqrt(gamma)*qt.sigmax()
    e1 = np.sqrt(gamma)*qt.sigmaz()
    c_ops = [qt.tensor(c1,c1)]
    e_ops = [qt.tensor(e1,e1),qt.tensor(c1,c1)]
    #e_ops = []
    tlist = np.linspace(0,10,100)
    ntraj = 2000
    ptrace_sel = [0]
    sol_f90 = mcf90.mcsolve_f90(H,psi0,tlist,c_ops,e_ops,ntraj=ntraj,
            ptrace_sel=ptrace_sel,calc_entropy=True)

    def qubit_integrate(self, tlist, psi0, epsilon, delta, g1, g2):

        H = epsilon / 2.0 * sigmaz() + delta / 2.0 * sigmax()

        c_op_list = []

        rate = g1
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * sigmam())

        rate = g2
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * sigmaz())

        output = mesolve(H, psi0, tlist, c_op_list, [sigmax(), sigmay(), sigmaz()])
        expt_list = output.expect[0], output.expect[1], output.expect[2]

        return expt_list[0], expt_list[1], expt_list[2]

 def get_reduced_dms(self, states, spin):
     """
     takes a number of states and returns a list of bloch vector of the 0th spin coordinates for each
     """
     sz = sigmaz()
     sy = sigmay()
     sx = sigmax()
     zs = []
     ys = []
     xs = []
     for state in states:
         ptrace = state.ptrace(spin)
         zval = abs(expect(sz, ptrace))
         yval = abs(expect(sy, ptrace))
         xval = abs(expect(sx, ptrace))
         zs.append(zval)
         ys.append(yval)
         xs.append(xval)
     return xs, ys, zs

    def testTLS(self):
        "brmesolve: qubit"

        delta = 0.0 * 2 * np.pi
        epsilon = 0.5 * 2 * np.pi
        gamma = 0.25
        times = np.linspace(0, 10, 100)
        H = delta/2 * sigmax() + epsilon/2 * sigmaz()
        psi0 = (2 * basis(2, 0) + basis(2, 1)).unit()
        c_ops = [np.sqrt(gamma) * sigmam()]
        a_ops = [sigmax()]
        e_ops = [sigmax(), sigmay(), sigmaz()]
        res_me = mesolve(H, psi0, times, c_ops, e_ops)
        res_brme = brmesolve(H, psi0, times, a_ops, e_ops,
                             spectra_cb=[lambda w: gamma * (w >= 0)])

        for idx, e in enumerate(e_ops):
            diff = abs(res_me.expect[idx] - res_brme.expect[idx]).max()
            assert_(diff < 1e-2)