Python qutip.expect函数代码示例

def testHOFiniteTemperatureStates():
    """
    brmesolve: harmonic oscillator, finite temperature, states
    """

    N = 10
    w0 = 1.0 * 2 * np.pi
    g = 0.05 * w0
    kappa = 0.25
    times = np.linspace(0, 25, 1000)
    a = destroy(N)
    H = w0 * a.dag() * a + g * (a + a.dag())
    psi0 = ket2dm((basis(N, 4) + basis(N, 2) + basis(N, 0)).unit())

    n_th = 1.5
    w_th = w0/np.log(1 + 1/n_th)

    def S_w(w):
        if w >= 0:
            return (n_th + 1) * kappa
        else:
            return (n_th + 1) * kappa * np.exp(w / w_th)

    c_ops = [np.sqrt(kappa * (n_th + 1)) * a, np.sqrt(kappa * n_th) * a.dag()]
    a_ops = [a + a.dag()]
    e_ops = []

    res_me = mesolve(H, psi0, times, c_ops, e_ops)
    res_brme = brmesolve(H, psi0, times, a_ops, e_ops, [S_w])

    n_me = expect(a.dag() * a, res_me.states)
    n_brme = expect(a.dag() * a, res_brme.states)

    diff = abs(n_me - n_brme).max()
    assert_(diff < 1e-2)

    def compute(N, wc, wa, glist, use_rwa):

        # Pre-compute operators for the hamiltonian
        a  = tensor(destroy(N), qeye(2))
        sm = tensor(qeye(N), destroy(2))
        nc = a.dag() * a
        na = sm.dag() * sm

        idx = 0
        na_expt = zeros(shape(glist))
        nc_expt = zeros(shape(glist))
        for g in glist:

            # recalculate the hamiltonian for each value of g
            if use_rwa:
                H = wc * nc + wa * na + g * (a.dag() * sm + a * sm.dag())
            else:
                H = wc * nc + wa * na + g * (a.dag() + a) * (sm + sm.dag())

            # find the groundstate of the composite system
            evals, ekets = H.eigenstates()
            psi_gnd = ekets[0]
            na_expt[idx] = expect(na, psi_gnd)
            nc_expt[idx] = expect(nc, psi_gnd)

            idx += 1

        return nc_expt, na_expt, ket2dm(psi_gnd)

 def correlator(self):
     """correlator
     Measure of quantum vs semiclassical"""
     if not self.precalc:
         self._calculate()
     return np.abs(
         np.asarray([qt.expect(self.a * self.sm, rho) for rho in self.rhos_ss])
         - np.asarray([qt.expect(self.a, rho) for rho in self.rhos_ss])
         * np.asarray([qt.expect(self.sm, rho) for rho in self.rhos_ss])
     ).T

 def g2(self):
     if not self.precalc:
         self._calculate()
     return np.abs(
         np.asarray(
             [
                 qt.expect(self.a.dag() * self.a.dag() * self.a * self.a, rho)
                 / qt.expect(self.a.dag() * self.a, rho) ** 2
                 for rho in self.rhos_ss
             ]
         )
     ).T

def test_ho_lgmres():
    "Steady state: Thermal HO - iterative-lgmres solver"
    # thermal steadystate of an oscillator: compare numerics with analytical
    # formula
    a = destroy(40)
    H = 0.5 * 2 * np.pi * a.dag() * a
    gamma1 = 0.05

    wth_vec = np.linspace(0.1, 3, 20)
    p_ss = np.zeros(np.shape(wth_vec))

    for idx, wth in enumerate(wth_vec):

        n_th = 1.0 / (np.exp(1.0 / wth) - 1)  # bath temperature
        c_op_list = []
        rate = gamma1 * (1 + n_th)
        c_op_list.append(np.sqrt(rate) * a)
        rate = gamma1 * n_th
        c_op_list.append(np.sqrt(rate) * a.dag())
        rho_ss = steadystate(H, c_op_list, method='iterative-lgmres')
        p_ss[idx] = np.real(expect(a.dag() * a, rho_ss))

    p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1)
    delta = sum(abs(p_ss_analytic - p_ss))
    assert_equal(delta < 1e-3, True)

def test_qubit_power():
    "Steady state: Thermal qubit - power solver"
    # thermal steadystate of a qubit: compare numerics with analytical formula
    sz = sigmaz()
    sm = destroy(2)

    H = 0.5 * 2 * np.pi * sz
    gamma1 = 0.05

    wth_vec = np.linspace(0.1, 3, 20)
    p_ss = np.zeros(np.shape(wth_vec))

    for idx, wth in enumerate(wth_vec):

        n_th = 1.0 / (np.exp(1.0 / wth) - 1)  # bath temperature
        c_op_list = []
        rate = gamma1 * (1 + n_th)
        c_op_list.append(np.sqrt(rate) * sm)
        rate = gamma1 * n_th
        c_op_list.append(np.sqrt(rate) * sm.dag())
        rho_ss = steadystate(H, c_op_list, method='power')
        p_ss[idx] = expect(sm.dag() * sm, rho_ss)

    p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec))
    delta = sum(abs(p_ss_analytic - p_ss))
    assert_equal(delta < 1e-5, True)

    def steady(N = 20):  # number of basis states to consider
        n=num(N)
        a = destroy(N)
        H = a.dag() * a
        print H.eigenstates()
        #psi0 = basis(N, 10)  # initial state
        kappa = 0.1  # coupling to oscillator
        c_op_list = []
        n_th_a = 2  # temperature with average of 2 excitations
        rate = kappa * (1 + n_th_a)
        c_op_list.append(sqrt(rate) * a)  # decay operators
        rate = kappa * n_th_a
        c_op_list.append(sqrt(rate) * a.dag())  # excitation operators
        final_state = steadystate(H, c_op_list)
        fexpt = expect(a.dag() * a, final_state)

        #tlist = linspace(0, 100, 100)

        #mcdata = mcsolve(H, psi0, tlist, c_op_list, [a.dag() * a], ntraj=100)

        #medata = mesolve(H, psi0, tlist, c_op_list, [a.dag() * a])
        #plot(tlist, mcdata.expect[0],
        #plt.plot(tlist, medata.expect[0], lw=2)
        plt.axhline(y=fexpt, color='r', lw=1.5) # ss expt. value as horiz line (= 2)
        plt.ylim([0, 10])
        plt.show()

    def find_expect(self, vg, pwr, fd):
        if self.power_plot:
            phi, pwr=vg
        else:
            phi, fd=vg
        pwr_fridge=pwr-self.atten
        lin_pwr=0.001*10**(pwr_fridge/10.0)
        Omega=sqrt(lin_pwr/h*2.0)

        gamma, Delta=self._get_GammaDelta(fd=fd, f0=self.f0, Np=self.Np, gamma=self.gamma)
        g_el=self._get_Gamma_C(fd=fd)
        wTvec=self._get_fTvec(phi=phi, gamma=gamma, Delta=Delta, fd=fd, Psaw=lin_pwr)

        if self.acoustic_plot:
            Om=Omega*sqrt(gamma/fd)
        else:
            Om=Omega*sqrt(g_el/fd)
        wT = wTvec-fd*self.nvec #rotating frame of gate drive \omega_m-m*\omega_\gate
        transmon_levels = Qobj(diag(wT[range(self.N_dim)]))
        rate1 = (gamma+g_el)*(1.0 + self.N_gamma)
        rate2 = (gamma+g_el)*self.N_gamma
        c_ops=[sqrt(rate1)*self.a_op, sqrt(rate2)*self.a_dag]#, sqrt(rate3)*self.a_op, sqrt(rate4)*self.a_dag]
        Omega_vec=-0.5j*(Om*self.a_dag - conj(Om)*self.a_op)
        H=transmon_levels +Omega_vec
        final_state = steadystate(H, c_ops) #solve master equation
        fexpt=expect(self.a_op, final_state) #expectation value of relaxation operator
        #return fexpt
        if self.acoustic_plot:
            return 1.0*gamma/Om*fexpt
        else:
            return 1.0*sqrt(g_el*gamma)/Om*fexpt

 def e_ops_func(t, rho, transformation=Utrans, e_ops=results.e_ops):
     """Transform the density matrix into the lab frame, and calculate the expectation values. 
     TODO: this could probably be streamlined... need to look into qutip's code
     """
     rho_lab_frame=Utrans(t).dag()*q.Qobj(rho)*Utrans(t)
     for i, e_operator in enumerate(e_ops):
         expectation_values[i][e_ops_func.idx]=q.expect(e_operator, rho_lab_frame).real
     e_ops_func.idx+=1

def essolve(H, rho0, tlist, c_op_list, expt_op_list):
    """
    Evolution of a state vector or density matrix (`rho0`) for a given
    Hamiltonian (`H`) and set of collapse operators (`c_op_list`), by expressing
    the ODE as an exponential series. The output is either the state vector at
    arbitrary points in time (`tlist`), or the expectation values of the supplied
    operators (`expt_op_list`). 
   
    Parameters
    ----------
    H : qobj/function_type 
        System Hamiltonian.
    
    rho0 : qobj 
        Initial state density matrix.
    
    tlist : list/array
        ``list`` of times for :math:`t`.
    
    c_op_list : list 
        ``list`` of ``qobj`` collapse operators.
    
    expt_op_list : list
        ``list`` of ``qobj`` operators for which to evaluate expectation values.


    Returns
    -------
     expt_array : array
        Expectation values of wavefunctions/density matrices for the times specified in ``tlist``.        

    
    .. note:: This solver does not support time-dependent Hamiltonians.

    """
    n_expt_op = len(expt_op_list)
    n_tsteps  = len(tlist)

    # Calculate the Liouvillian
    if c_op_list == None or len(c_op_list) == 0:
        L = H
    else:
        L = liouvillian(H, c_op_list)

    es = ode2es(L, rho0)

    # evaluate the expectation values      
    if n_expt_op == 0:
        result_list = [Qobj() for k in range(n_tsteps)]
    else:
        result_list = zeros([n_expt_op, n_tsteps], dtype=complex)

    for n in range(0, n_expt_op):
        result_list[n,:] = esval(expect(expt_op_list[n],es),tlist)

    return result_list

 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 testCoherentState(self):
     """
     states: coherent state
     """
     N = 10
     alpha = 0.5
     c1 = coherent(N, alpha) # displacement method
     c2 = coherent(7, alpha, offset=3) # analytic method
     assert_(abs(expect(destroy(N), c1) - alpha) < 1e-10)
     assert_((c1[3:]-c2).norm() < 1e-7)

def get_dms(states):
    '''
    takes a number of states and returns a list of bloch vector coordinates for each
    '''
    si = qeye(2)
    sz = sigmaz()
    sy = sigmay()
    sx = sigmax()
    zs = []
    ys = []
    xs = []
    for state in states:
        ptrace = state.ptrace(0)
        zval = expect(sz, ptrace )
        yval = expect(sy, ptrace )
        xval = expect(sx, ptrace )
        zs.append(zval)
        ys.append(yval)
        xs.append(xval)
    return xs,ys,zs

def test_SparseHermValsVecs():
    """
    Sparse eigs Hermitian
    """

    # check using number operator
    N = num(10)
    spvals, spvecs = N.eigenstates(sparse=True)
    for k in range(10):
        # check that eigvals are in proper order
        assert_equal(abs(spvals[k] - k) <= 1e-13, True)
        # check that eigenvectors are right and in right order
        assert_equal(abs(expect(N, spvecs[k]) - spvals[k]) < 5e-14, True)

    # check ouput of only a few eigenvals/vecs
    spvals, spvecs = N.eigenstates(sparse=True, eigvals=7)
    assert_equal(len(spvals), 7)
    assert_equal(spvals[0] <= spvals[-1], True)
    for k in range(7):
        assert_equal(abs(spvals[k] - k) < 1e-12, True)

    spvals, spvecs = N.eigenstates(sparse=True, sort='high', eigvals=5)
    assert_equal(len(spvals), 5)
    assert_equal(spvals[0] >= spvals[-1], True)
    vals = np.arange(9, 4, -1)
    for k in range(5):
        # check that eigvals are ordered from high to low
        assert_equal(abs(spvals[k] - vals[k]) < 5e-14, True)
        assert_equal(abs(expect(N, spvecs[k]) - vals[k]) < 1e-14, True)
    # check using random Hermitian
    H = rand_herm(10)
    spvals, spvecs = H.eigenstates(sparse=True)
    # check that sorting is lowest eigval first
    assert_equal(spvals[0] <= spvals[-1], True)
    # check that spvals equal expect vals
    for k in range(10):
        assert_equal(abs(expect(H, spvecs[k]) - spvals[k]) < 5e-14, True)
        # check that ouput is real for Hermitian operator
        assert_equal(np.isreal(spvals[k]), True)

    def jc_steadystate(self, N, wc, wa, g, kappa, gamma,
                       pump, psi0, use_rwa, tlist):

        # Hamiltonian
        a = tensor(destroy(N), identity(2))
        sm = tensor(identity(N), destroy(2))

        if use_rwa:
            # use the rotating wave approxiation
            H = wc * a.dag(
            ) * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
        else:
            H = wc * a.dag() * a + wa * sm.dag() * sm + g * (
                a.dag() + a) * (sm + sm.dag())

        # collapse operators
        c_op_list = []

        n_th_a = 0.0  # zero temperature

        rate = kappa * (1 + n_th_a)
        c_op_list.append(np.sqrt(rate) * a)

        rate = kappa * n_th_a
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * a.dag())

        rate = gamma
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * sm)

        rate = pump
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * sm.dag())

        # find the steady state
        rho_ss = steadystate(H, c_op_list)

        return expect(a.dag() * a, rho_ss), expect(sm.dag() * sm, rho_ss)