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 testFloquetUnitary(self):
        """
        Floquet: test unitary evolution of time-dependent two-level system
        """

        delta = 1.0 * 2 * np.pi
        eps0 = 1.0 * 2 * np.pi
        A = 0.5 * 2 * np.pi
        omega = np.sqrt(delta ** 2 + eps0 ** 2)
        T = (2 * np.pi) / omega
        tlist = np.linspace(0.0, 2 * T, 101)
        psi0 = rand_ket(2)
        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
        H1 = A / 2.0 * sigmax()
        args = {'w': omega}
        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
        e_ops = [num(2)]

        # solve schrodinger equation with floquet solver
        sol = fsesolve(H, psi0, tlist, e_ops, T, args)

        # compare with results from standard schrodinger equation
        sol_ref = mesolve(H, psi0, tlist, [], e_ops, args)

        assert_(max(abs(sol.expect[0] - sol_ref.expect[0])) < 1e-4)

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 test_crab(self):
        """
        Optimise pulse for Hadamard gate using CRAB algorithm
        Apply guess and ramping pulse
        assert that goal is achieved and fidelity error is below threshold
        assert that starting amplitude is zero
        """
        # Hadamard
        H_d = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 12
        evo_time = 10
        
        # Run the optimisation
        result = cpo.opt_pulse_crab_unitary(H_d, H_c, U_0, U_targ, 
                n_ts, evo_time, 
                fid_err_targ=1e-5, 
                alg_params={'crab_pulse_params':{'randomize_coeffs':False, 
                                                 'randomize_freqs':False}},
                init_coeff_scaling=0.5,
                guess_pulse_type='GAUSSIAN', 
                guess_pulse_params={'variance':0.1*evo_time},
                guess_pulse_scaling=1.0, guess_pulse_offset=1.0,
                amp_lbound=None, amp_ubound=None,
                ramping_pulse_type='GAUSSIAN_EDGE', 
                ramping_pulse_params={'decay_time':evo_time/100.0},
                gen_stats=True)
        assert_(result.goal_achieved, msg="Hadamard goal not achieved")
        assert_almost_equal(result.fid_err, 0.0, decimal=3, 
                            err_msg="Hadamard infidelity too high")
        assert_almost_equal(result.final_amps[0, 0], 0.0, decimal=3, 
                            err_msg="lead in amplitude not zero")

    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 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 test_01_1_unitary_hadamard(self):
        """
        control.pulseoptim: Hadamard gate with linear initial pulses
        assert that goal is achieved and fidelity error is below threshold
        """
        # Hadamard
        H_d = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 10
        evo_time = 10

        # Run the optimisation
        result = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-10,
                        init_pulse_type='LIN',
                        gen_stats=True)
        assert_(result.goal_achieved, msg="Hadamard goal not achieved. "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result.termination_reason, result.fid_err))
        assert_almost_equal(result.fid_err, 0.0, decimal=10,
                            err_msg="Hadamard infidelity too high")

 def count_waves(n_ts, evo_time, ptype, freq=None, num_waves=None):
     
     # Any dyn config will do 
     #Hadamard
     H_d = sigmaz()
     H_c = [sigmax()]
     U_0 = identity(2)
     U_targ = hadamard_transform(1)
     
     pulse_params = {}
     if freq is not None:
         pulse_params['freq'] = freq
     if num_waves is not None:
         pulse_params['num_waves'] = num_waves
     
     optim = cpo.create_pulse_optimizer(H_d, H_c, U_0, U_targ, 
                                 n_ts, evo_time, 
                                 dyn_type='UNIT', 
                                 init_pulse_type=ptype,
                                 init_pulse_params=pulse_params,
                                 gen_stats=False)
     pgen = optim.pulse_generator
     pulse = pgen.gen_pulse()
     
     # count number of waves
     zero_cross = pulse[0:-2]*pulse[1:-1] < 0
     
     return (sum(zero_cross) + 1) / 2

    def test_9_time_dependent_drift(self):
        """
        control.pulseoptim: Hadamard gate with fixed and time varying drift
        assert that goal is achieved for both and that different control
        pulses are produced (only) when they should be
        """
        # Hadamard
        H_0 = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 20
        evo_time = 10
        
        drift_amps_flat = np.ones([n_ts], dtype=float)
        dript_amps_step = [np.round(float(k)/n_ts) for k in range(n_ts)]
        
        # Run the optimisations
        result_fixed = cpo.optimize_pulse_unitary(H_0, H_c, U_0, U_targ, 
                        n_ts, evo_time, 
                        fid_err_targ=1e-10, 
                        init_pulse_type='LIN', 
                        gen_stats=True)
        assert_(result_fixed.goal_achieved, 
                    msg="Fixed drift goal not achieved. "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result_fixed.termination_reason, result_fixed.fid_err))
                    
        H_d = [drift_amps_flat[k]*H_0 for k in range(n_ts)]
        result_flat = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ, 
                        n_ts, evo_time, 
                        fid_err_targ=1e-10, 
                        init_pulse_type='LIN', 
                        gen_stats=True)
        assert_(result_flat.goal_achieved, msg="Flat drift goal not achieved. "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result_flat.termination_reason, result_flat.fid_err))
                    
        # Check fixed and flat produced the same pulse
        assert_almost_equal(result_fixed.final_amps, result_flat.final_amps, 
                            decimal=9, 
                            err_msg="Flat and fixed drift result in "
                                    "different control pules")
                            
        H_d = [dript_amps_step[k]*H_0 for k in range(n_ts)]
        result_step = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ, 
                        n_ts, evo_time, 
                        fid_err_targ=1e-10, 
                        init_pulse_type='LIN', 
                        gen_stats=True)
        assert_(result_step.goal_achieved, msg="Step drift goal not achieved. "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result_step.termination_reason, result_step.fid_err))
                    
        # Check step and flat produced different results
        assert_(np.any(
            np.abs(result_flat.final_amps - result_step.final_amps) > 1e-3), 
                            msg="Flat and step drift result in "
                                    "the same control pules")

    def test_01_6_unitary_hadamard_grad(self):
        """
        control.pulseoptim: Hadamard gate gradient check
        assert that gradient approx and exact gradient match in tolerance
        """
        # Hadamard
        H_d = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 10
        evo_time = 10

        # Create the optim objects
        optim = cpo.create_pulse_optimizer(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-10,
                        dyn_type='UNIT',
                        init_pulse_type='LIN',
                        gen_stats=True)
        dyn = optim.dynamics

        init_amps = optim.pulse_generator.gen_pulse().reshape([-1, 1])
        dyn.initialize_controls(init_amps)

        # Check the exact gradient
        func = optim.fid_err_func_wrapper
        grad = optim.fid_err_grad_wrapper
        x0 = dyn.ctrl_amps.flatten()
        grad_diff = check_grad(func, grad, x0)
        assert_almost_equal(grad_diff, 0.0, decimal=6,
                            err_msg="Unitary gradient outside tolerance")

    def test_01_2_unitary_hadamard_no_stats(self):
        """
        control.pulseoptim: Hadamard gate with linear initial pulses (no stats)
        assert that goal is achieved
        """
        # Hadamard
        H_d = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 10
        evo_time = 10

        # Run the optimisation
        #Try without stats
        result = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-10,
                        init_pulse_type='LIN',
                        gen_stats=False)
        assert_(result.goal_achieved, msg="Hadamard goal not achieved "
                                            "(no stats). "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result.termination_reason, result.fid_err))

    def test_01_4_unitary_hadamard_qobj(self):
        """
        control.pulseoptim: Hadamard gate with linear initial pulses (Qobj)
        assert that goal is achieved
        """
        # Hadamard
        H_d = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 10
        evo_time = 10

        # Run the optimisation
        #Try with Qobj propagation
        result = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-10,
                        init_pulse_type='LIN',
                        dyn_params={'oper_dtype':Qobj},
                        gen_stats=True)
        assert_(result.goal_achieved, msg="Hadamard goal not achieved "
                                            "(Qobj propagation). "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result.termination_reason, result.fid_err))

 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 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 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 __init__(self):
     super(TestBloch, self).__init__()
     self.propagator = propagator(sigmaz(), .1, [])
     self.set_state(qeye(2) + sigmax())
     self.timer = QtCore.QTimer()
     self.timer.setInterval(100)
     self.timer.timeout.connect(self.propagate)
     self.timer.start()

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)