def qpt(U, op_basis_list):
    """
    Calculate the quantum process tomography chi matrix for a given 
    (possibly nonunitary) transformation matrix U, which transforms a 
    density matrix in vector form according to:

        vec(rho) = U * vec(rho0)

        or

        rho = vec2mat(U * mat2vec(rho0))

    U can be calculated for an open quantum system using the QuTiP propagator
    function.
    """

    E_ops = []
    # loop over all index permutations
    for inds in index_permutations([len(op_list) for op_list in op_basis_list]):
        # loop over all composite systems
        E_op_list = [op_basis_list[k][inds[k]] for k in range(len(op_basis_list))]
        E_ops.append(tensor(E_op_list))

    EE_ops = [spre(E1) * spost(E2.dag()) for E1 in E_ops for E2 in E_ops]

    M = hstack([mat2vec(EE.full()) for EE in EE_ops])

    Uvec = mat2vec(U.full())

    chi_vec = la.solve(M, Uvec)

    return vec2mat(chi_vec)

    def test_ComplexSuperApply(self):
        """
        Superoperator: Efficient numerics and reference return same result,
        acting on non-composite system
        """
        rho_list = list(map(rand_dm, [2, 3, 2, 3, 2]))
        rho_input = tensor(rho_list)
        superop = kraus_to_super(rand_kraus_map(3))

        analytic_result = rho_list
        analytic_result[1] = Qobj(vec2mat(superop.data.todense() *
                                  mat2vec(analytic_result[1].data.todense())))
        analytic_result[3] = Qobj(vec2mat(superop.data.todense() *
                                  mat2vec(analytic_result[3].data.todense())))
        analytic_result = tensor(analytic_result)

        naive_result = subsystem_apply(rho_input, superop,
                                       [False, True, False, True, False],
                                       reference=True)
        naive_diff = (analytic_result - naive_result).data.todense()
        assert_(norm(naive_diff) < 1e-12)

        efficient_result = subsystem_apply(rho_input, superop,
                                           [False, True, False, True, False])
        efficient_diff = (efficient_result - analytic_result).data.todense()
        assert_(norm(efficient_diff) < 1e-12)

def test_vec_to_eigbasis():
    "BR Tools : vector to eigenbasis"
    N = 10
    for kk in range(50):
        H = rand_herm(N,0.5)
        h = H.full('F')
        R = rand_dm(N,0.5)
        r = mat2vec(R.full()).ravel()
        ans = mat2vec(R.transform(H.eigenstates()[1]).full()).ravel()
        out = _test_vec_to_eigbasis(h, r)
        assert_(np.allclose(ans,out))

def test_eigvec_to_fockbasis():
    "BR Tools : eigvector to fockbasis"
    N = 10
    for kk in range(50):
        H = rand_herm(N,0.5)
        h = H.full('F')
        R = rand_dm(N,0.5)
        r = mat2vec(R.full()).ravel()
        eigvals = np.zeros(N,dtype=float)
        Z = _test_zheevr(H.full('F'), eigvals)
        eig_vec = mat2vec(R.transform(H.eigenstates()[1]).full()).ravel()
        out = _test_eigvec_to_fockbasis(eig_vec, Z, N)
        assert_(np.allclose(r,out))

def smesolve_generic(H, rho0, tlist, c_ops, e_ops, rhs, d1, d2, ntraj, nsubsteps):
    """
    internal

    .. note::

        Experimental.

    """
    if debug:
        print(inspect.stack()[0][3])

    N_store = len(tlist)
    N_substeps = nsubsteps
    N = N_store * N_substeps
    dt = (tlist[1] - tlist[0]) / N_substeps

    print("N = %d. dt=%.2e" % (N, dt))

    data = Odedata()

    data.expect = np.zeros((len(e_ops), N_store), dtype=complex)

    # pre-compute collapse operator combinations that are commonly needed
    # when evaluating the RHS of stochastic master equations
    A_ops = []
    for c_idx, c in enumerate(c_ops):

        # xxx: precompute useful operator expressions...
        cdc = c.dag() * c
        Ldt = spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc)
        LdW = spre(c) + spost(c.dag())
        Lm = spre(c) + spost(c.dag())  # currently same as LdW

        A_ops.append([Ldt.data, LdW.data, Lm.data])

    # Liouvillian for the unitary part
    L = -1.0j * (spre(H) - spost(H))  # XXX: should we split the ME in stochastic
    # and deterministic collapse operators here?

    progress_acc = 0.0
    for n in range(ntraj):

        if debug and (100 * float(n) / ntraj) >= progress_acc:
            print("Progress: %.2f" % (100 * float(n) / ntraj))
            progress_acc += 10.0

        rho_t = mat2vec(rho0.full())

        states_list = _smesolve_single_trajectory(
            L, dt, tlist, N_store, N_substeps, rho_t, A_ops, e_ops, data, rhs, d1, d2
        )

        # if average -> average...
        data.states.append(states_list)

    # average
    data.expect = data.expect / ntraj

    return data

def smepdpsolve_generic(ssdata, options, progress_bar):
    """
    For internal use.

    .. note::

        Experimental.

    """
    if debug:
        print(inspect.stack()[0][3])

    N_store = len(ssdata.tlist)
    N_substeps = ssdata.nsubsteps
    N = N_store * N_substeps
    dt = (ssdata.tlist[1] - ssdata.tlist[0]) / N_substeps
    NT = ssdata.ntraj

    data = Odedata()
    data.solver = "smepdpsolve"
    data.times = ssdata.tlist
    data.expect = np.zeros((len(ssdata.e_ops), N_store), dtype=complex)
    data.jump_times = []
    data.jump_op_idx = []

    # Liouvillian for the deterministic part.
    # needs to be modified for TD systems
    L = liouvillian_fast(ssdata.H, ssdata.c_ops)
        
    progress_bar.start(ssdata.ntraj)

    for n in range(ssdata.ntraj):
        progress_bar.update(n)
        rho_t = mat2vec(ssdata.rho0.full()).ravel()

        states_list, jump_times, jump_op_idx = \
            _smepdpsolve_single_trajectory(data, L, dt, ssdata.tlist,
                                           N_store, N_substeps,
                                           rho_t, ssdata.c_ops, ssdata.e_ops)

        data.states.append(states_list)
        data.jump_times.append(jump_times)
        data.jump_op_idx.append(jump_op_idx)

    progress_bar.finished()

    # average density matrices
    if options.average_states and np.any(data.states):
        data.states = [sum(state_list).unit() for state_list in data.states]
    
    # average
    data.expect = data.expect / ssdata.ntraj

    # standard error
    if NT > 1:
        data.se = (data.ss - NT * (data.expect ** 2)) / (NT * (NT - 1))
    else:
        data.se = None

    return data

def smesolve_generic(H, rho0, tlist, c_ops, sc_ops, e_ops,
                     rhs, d1, d2, d2_len, ntraj, nsubsteps,
                     options, progress_bar):
    """
    internal

    .. note::

        Experimental.

    """
    if debug:
        print(inspect.stack()[0][3])

    N_store = len(tlist)
    N_substeps = nsubsteps
    N = N_store * N_substeps
    dt = (tlist[1] - tlist[0]) / N_substeps

    data = Odedata()
    data.solver = "smesolve"
    data.times = tlist
    data.expect = np.zeros((len(e_ops), N_store), dtype=complex)

    # pre-compute collapse operator combinations that are commonly needed
    # when evaluating the RHS of stochastic master equations
    A_ops = []
    for c_idx, c in enumerate(sc_ops):

        # xxx: precompute useful operator expressions...
        cdc = c.dag() * c
        Ldt = spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc)
        LdW = spre(c) + spost(c.dag())
        Lm = spre(c) + spost(c.dag())  # currently same as LdW

        A_ops.append([Ldt.data, LdW.data, Lm.data])

    # Liouvillian for the deterministic part
    L = liouvillian_fast(H, c_ops)  # needs to be modified for TD systems

    progress_bar.start(ntraj)

    for n in range(ntraj):
        progress_bar.update(n)

        rho_t = mat2vec(rho0.full())

        states_list = _smesolve_single_trajectory(
            L, dt, tlist, N_store, N_substeps,
            rho_t, A_ops, e_ops, data, rhs, d1, d2, d2_len)

        # if average -> average...
        data.states.append(states_list)

    progress_bar.finished()

    # average
    data.expect = data.expect / ntraj

    return data

def _mesolve_const(H, rho0, tlist, c_op_list, e_ops, args, opt,
                   progress_bar):
    """
    Evolve the density matrix using an ODE solver, for constant hamiltonian
    and collapse operators.
    """

    if debug:
        print(inspect.stack()[0][3])

    #
    # check initial state
    #
    if isket(rho0):
        # if initial state is a ket and no collapse operator where given,
        # fall back on the unitary schrodinger equation solver
        if len(c_op_list) == 0 and isoper(H):
            return _sesolve_const(H, rho0, tlist, e_ops, args, opt,
                                  progress_bar)

        # Got a wave function as initial state: convert to density matrix.
        rho0 = ket2dm(rho0)

    #
    # construct liouvillian
    #
    if opt.tidy:
        H = H.tidyup(opt.atol)

    L = liouvillian(H, c_op_list)
    

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full()).ravel('F')
    if issuper(rho0):
        r = scipy.integrate.ode(_ode_super_func)
        r.set_f_params(L.data)
    else:
        if opt.use_openmp and L.data.nnz >= qset.openmp_thresh:
            r = scipy.integrate.ode(cy_ode_rhs_openmp)
            r.set_f_params(L.data.data, L.data.indices, L.data.indptr, 
                            opt.openmp_threads)
        else:
            r = scipy.integrate.ode(cy_ode_rhs)
            r.set_f_params(L.data.data, L.data.indices, L.data.indptr)
        # r = scipy.integrate.ode(_ode_rho_test)
        # r.set_f_params(L.data)
    r.set_integrator('zvode', method=opt.method, order=opt.order,
                     atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
                     first_step=opt.first_step, min_step=opt.min_step,
                     max_step=opt.max_step)
    r.set_initial_value(initial_vector, tlist[0])

    #
    # call generic ODE code
    #
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)

def _smepdpsolve_generic(sso, options, progress_bar):
    """
    For internal use. See smepdpsolve.
    """
    if debug:
        logger.debug(inspect.stack()[0][3])

    N_store = len(sso.times)
    N_substeps = sso.nsubsteps
    dt = (sso.times[1] - sso.times[0]) / N_substeps
    nt = sso.ntraj

    data = Result()
    data.solver = "smepdpsolve"
    data.times = sso.times
    data.expect = np.zeros((len(sso.e_ops), N_store), dtype=complex)
    data.jump_times = []
    data.jump_op_idx = []

    # Liouvillian for the deterministic part.
    # needs to be modified for TD systems
    L = liouvillian(sso.H, sso.c_ops)

    progress_bar.start(sso.ntraj)

    for n in range(sso.ntraj):
        progress_bar.update(n)
        rho_t = mat2vec(sso.rho0.full()).ravel()

        states_list, jump_times, jump_op_idx = \
            _smepdpsolve_single_trajectory(data, L, dt, sso.times,
                                           N_store, N_substeps,
                                           rho_t, sso.rho0.dims,
                                           sso.c_ops, sso.e_ops)

        data.states.append(states_list)
        data.jump_times.append(jump_times)
        data.jump_op_idx.append(jump_op_idx)

    progress_bar.finished()

    # average density matrices
    if options.average_states and np.any(data.states):
        data.states = [sum([data.states[m][n] for m in range(nt)]).unit()
                       for n in range(len(data.times))]

    # average
    data.expect = data.expect / sso.ntraj

    # standard error
    if nt > 1:
        data.se = (data.ss - nt * (data.expect ** 2)) / (nt * (nt - 1))
    else:
        data.se = None

    return data

def test_vector_roundtrip():
    "BR Tools : vector roundtrip transform"
    N = 10
    for kk in range(50):
        H = rand_herm(N,0.5)
        h = H.full('F')
        R = rand_dm(N,0.5)
        r = mat2vec(R.full()).ravel()
        out = _test_vector_roundtrip(h,r)
        assert_(np.allclose(r,out))

def countstat_current(L, c_ops=None, rhoss=None, J_ops=None):
    """
    Calculate the current corresponding a system Liouvillian `L` and a list of
    current collapse operators `c_ops` or current superoperators `J_ops`
    (either must be specified). Optionally the steadystate density matrix
    `rhoss` and a list of current superoperators `J_ops` can be specified. If
    either of these are omitted they are computed internally.

    Parameters
    ----------

    L : :class:`qutip.Qobj`
        Qobj representing the system Liouvillian.

    c_ops : array / list (optional)
        List of current collapse operators.

    rhoss : :class:`qutip.Qobj` (optional)
        The steadystate density matrix corresponding the system Liouvillian
        `L`.

    J_ops : array / list (optional)
        List of current superoperators.

    Returns
    --------
    I : array
        The currents `I` corresponding to each current collapse operator
        `c_ops` (or, equivalently, each current superopeator `J_ops`).
    """

    if J_ops is None:
        if c_ops is None:
            raise ValueError("c_ops must be given if J_ops is not")
        J_ops = [sprepost(c, c.dag()) for c in c_ops]

    if rhoss is None:
        if c_ops is None:
            raise ValueError("c_ops must be given if rhoss is not")
        rhoss = steadystate(L, c_ops)

    rhoss_vec = mat2vec(rhoss.full()).ravel()

    N = len(J_ops)
    I = np.zeros(N)

    for i, Ji in enumerate(J_ops):
        I[i] = expect_rho_vec(Ji.data, rhoss_vec, 1)

    return I

def _spectrum_pi(H, wlist, c_ops, a_op, b_op, use_pinv=False):
    """
    Internal function for calculating the spectrum of the correlation function
    :math:`\left<A(\\tau)B(0)\\right>`.
    """

    L = H if issuper(H) else liouvillian(H, c_ops)

    tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
    N = np.prod(L.dims[0][0])

    A = L.full()
    b = spre(b_op).full()
    a = spre(a_op).full()

    tr_vec = np.transpose(mat2vec(tr_mat.full()))

    rho_ss = steadystate(L)
    rho = np.transpose(mat2vec(rho_ss.full()))

    I = np.identity(N * N)
    P = np.kron(np.transpose(rho), tr_vec)
    Q = I - P

    spectrum = np.zeros(len(wlist))

    for idx, w in enumerate(wlist):
        if use_pinv:
            MMR = np.linalg.pinv(-1.0j * w * I + A)
        else:
            MMR = np.dot(Q, np.linalg.solve(-1.0j * w * I + A, Q))

        s = np.dot(tr_vec,
                   np.dot(a, np.dot(MMR, np.dot(b, np.transpose(rho)))))
        spectrum[idx] = -2 * np.real(s[0, 0])

    return spectrum

    def test_SimpleSuperApply(self):
        """
        Non-composite system, operator on Liouville space.
        """
        rho_3 = rand_dm(3)
        superop = kraus_to_super(rand_kraus_map(3))
        analytic_result = vec2mat(superop.data.todense() *
                                  mat2vec(rho_3.data.todense()))

        naive_result = subsystem_apply(rho_3, superop, [True],
                                       reference=True)
        naive_diff = (analytic_result - naive_result).data.todense()
        assert_(norm(naive_diff) < 1e-12)

        efficient_result = subsystem_apply(rho_3, superop, [True])
        efficient_diff = (efficient_result - analytic_result).data.todense()
        assert_(norm(efficient_diff) < 1e-12)

    def test_SimpleSuperApply(self):
        """
        Non-composite system, operator on Liouville space.
        """
        tol = 1e-12
        rho_3 = rand_dm(3)
        superop = kraus_to_super(rand_kraus_map(3))
        analytic_result = vec2mat(superop.data.todense() * mat2vec(rho_3.data.todense()))

        naive_result = subsystem_apply(rho_3, superop, [True], reference=True)
        naive_diff = (analytic_result - naive_result).data.todense()
        naive_diff_norm = norm(naive_diff)
        assert_(
            naive_diff_norm < tol,
            msg="SimpleSuper: naive_diff_norm {} " "is beyond tolerance {}".format(naive_diff_norm, tol),
        )

        efficient_result = subsystem_apply(rho_3, superop, [True])
        efficient_diff = (efficient_result - analytic_result).data.todense()
        efficient_diff_norm = norm(efficient_diff)
        assert_(
            efficient_diff_norm < tol,
            msg="SimpleSuper: efficient_diff_norm {} " "is beyond tolerance {}".format(efficient_diff_norm, tol),
        )

def smesolve_generic(ssdata, options, progress_bar):
    """
    internal

    .. note::

        Experimental.

    """
    if debug:
        print(inspect.stack()[0][3])

    N_store = len(ssdata.tlist)
    N_substeps = ssdata.nsubsteps
    N = N_store * N_substeps
    dt = (ssdata.tlist[1] - ssdata.tlist[0]) / N_substeps
    NT = ssdata.ntraj

    data = Odedata()
    data.solver = "smesolve"
    data.times = ssdata.tlist
    data.expect = np.zeros((len(ssdata.e_ops), N_store), dtype=complex)
    data.ss = np.zeros((len(ssdata.e_ops), N_store), dtype=complex)
    data.noise = []
    data.measurement = []

    # pre-compute suporoperator operator combinations that are commonly needed
    # when evaluating the RHS of stochastic master equations
    A_ops = []
    for c_idx, c in enumerate(ssdata.sc_ops):

        n = c.dag() * c
        A_ops.append([spre(c).data, spost(c).data,
                      spre(c.dag()).data, spost(c.dag()).data,
                      spre(n).data, spost(n).data,
                      (spre(c) * spost(c.dag())).data,
                      lindblad_dissipator(c, data_only=True)])

    s_e_ops = [spre(e) for e in ssdata.e_ops]

    # Liouvillian for the deterministic part.
    # needs to be modified for TD systems
    L = liouvillian_fast(ssdata.H, ssdata.c_ops)

    progress_bar.start(ssdata.ntraj)

    for n in range(ssdata.ntraj):
        progress_bar.update(n)

        rho_t = mat2vec(ssdata.state0.full()).ravel()

        noise = ssdata.noise[n] if ssdata.noise else None

        states_list, dW, m = _smesolve_single_trajectory(
            L, dt, ssdata.tlist, N_store, N_substeps,
            rho_t, A_ops, s_e_ops, data, ssdata.rhs,
            ssdata.d1, ssdata.d2, ssdata.d2_len, ssdata.homogeneous,
            ssdata.distribution, ssdata.args,
            store_measurement=ssdata.store_measurement,
            store_states=ssdata.store_states, noise=noise)

        data.states.append(states_list)
        data.noise.append(dW)
        data.measurement.append(m)

    progress_bar.finished()

    # average density matrices
    if options.average_states and np.any(data.states):
        data.states = [sum(state_list).unit() for state_list in data.states]

    # average
    data.expect = data.expect / NT

    # standard error
    if NT > 1:
        data.se = (data.ss - NT * (data.expect ** 2)) / (NT * (NT - 1))
    else:
        data.se = None

    # convert complex data to real if hermitian
    data.expect = [np.real(data.expect[n,:]) if e.isherm else data.expect[n,:]
                   for n, e in enumerate(ssdata.e_ops)]

    return data

def floquet_markov_mesolve(R, ekets, rho0, tlist, e_ops, f_modes_table=None,
                           options=None, floquet_basis=True):
    """
    Solve the dynamics for the system using the Floquet-Markov master equation.
    """

    if options is None:
        opt = Options()
    else:
        opt = options

    if opt.tidy:
        R.tidyup()

    #
    # check initial state
    #
    if isket(rho0):
        # Got a wave function as initial state: convert to density matrix.
        rho0 = ket2dm(rho0)

    #
    # prepare output array
    #
    n_tsteps = len(tlist)
    dt = tlist[1] - tlist[0]

    output = Result()
    output.solver = "fmmesolve"
    output.times = tlist

    if isinstance(e_ops, FunctionType):
        n_expt_op = 0
        expt_callback = True

    elif isinstance(e_ops, list):

        n_expt_op = len(e_ops)
        expt_callback = False

        if n_expt_op == 0:
            output.states = []
        else:
            if not f_modes_table:
                raise TypeError("The Floquet mode table has to be provided " +
                                "when requesting expectation values.")

            output.expect = []
            output.num_expect = n_expt_op
            for op in e_ops:
                if op.isherm:
                    output.expect.append(np.zeros(n_tsteps))
                else:
                    output.expect.append(np.zeros(n_tsteps, dtype=complex))

    else:
        raise TypeError("Expectation parameter must be a list or a function")

    #
    # transform the initial density matrix to the eigenbasis: from
    # computational basis to the floquet basis
    #
    if ekets is not None:
        rho0 = rho0.transform(ekets)

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full())
    r = scipy.integrate.ode(cy_ode_rhs)
    r.set_f_params(R.data.data, R.data.indices, R.data.indptr)
    r.set_integrator('zvode', method=opt.method, order=opt.order,
                     atol=opt.atol, rtol=opt.rtol, max_step=opt.max_step)
    r.set_initial_value(initial_vector, tlist[0])

    #
    # start evolution
    #
    rho = Qobj(rho0)

    t_idx = 0
    for t in tlist:
        if not r.successful():
            break

        rho.data = vec2mat(r.y)

        if expt_callback:
            # use callback method
            if floquet_basis:
                e_ops(t, Qobj(rho))
            else:
                f_modes_table_t, T = f_modes_table
                f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)
                e_ops(t, Qobj(rho).transform(f_modes_t, True))
        else:
            # calculate all the expectation values, or output rho if
            # no operators
            if n_expt_op == 0:
                if floquet_basis:
#.........这里部分代码省略.........

def _mesolve_list_td(H_func, rho0, tlist, c_op_list, e_ops, args, opt, progress_bar):
    """
    Evolve the density matrix using an ODE solver with time dependent
    Hamiltonian.
    """

    if debug:
        print(inspect.stack()[0][3])

    #
    # check initial state
    #
    if isket(rho0):
        # if initial state is a ket and no collapse operator where given,
        # fall back on the unitary schrodinger equation solver
        if len(c_op_list) == 0:
            return _sesolve_list_td(H_func, rho0, tlist, e_ops, args, opt, progress_bar)

        # Got a wave function as initial state: convert to density matrix.
        rho0 = ket2dm(rho0)

    #
    # construct liouvillian
    #
    if len(H_func) != 2:
        raise TypeError("Time-dependent Hamiltonian list must have two terms.")
    if not isinstance(H_func[0], (list, np.ndarray)) or len(H_func[0]) <= 1:
        raise TypeError("Time-dependent Hamiltonians must be a list " + "with two or more terms")
    if (not isinstance(H_func[1], (list, np.ndarray))) or (len(H_func[1]) != (len(H_func[0]) - 1)):
        raise TypeError(
            "Time-dependent coefficients must be list with "
            + "length N-1 where N is the number of "
            + "Hamiltonian terms."
        )

    if opt.rhs_reuse and config.tdfunc is None:
        rhs_generate(H_func, args)

    lenh = len(H_func[0])
    if opt.tidy:
        H_func[0] = [(H_func[0][k]).tidyup() for k in range(lenh)]
    L_func = [[liouvillian(H_func[0][0], c_op_list)], H_func[1]]
    for m in range(1, lenh):
        L_func[0].append(liouvillian(H_func[0][m], []))

    # create data arrays for time-dependent RHS function
    Ldata = [L_func[0][k].data.data for k in range(lenh)]
    Linds = [L_func[0][k].data.indices for k in range(lenh)]
    Lptrs = [L_func[0][k].data.indptr for k in range(lenh)]
    # setup ode args string
    string = ""
    for k in range(lenh):
        string += "Ldata[%d], Linds[%d], Lptrs[%d]," % (k, k, k)

    if args:
        td_consts = args.items()
        for elem in td_consts:
            string += str(elem[1])
            if elem != td_consts[-1]:
                string += ","

    # run code generator
    if not opt.rhs_reuse or config.tdfunc is None:
        if opt.rhs_filename is None:
            config.tdname = "rhs" + str(os.getpid()) + str(config.cgen_num)
        else:
            config.tdname = opt.rhs_filename
        cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args, config=config)
        cgen.generate(config.tdname + ".pyx")

        code = compile("from " + config.tdname + " import cy_td_ode_rhs", "<string>", "exec")
        exec(code, globals())
        config.tdfunc = cy_td_ode_rhs

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full()).ravel()
    r = scipy.integrate.ode(config.tdfunc)
    r.set_integrator(
        "zvode",
        method=opt.method,
        order=opt.order,
        atol=opt.atol,
        rtol=opt.rtol,
        nsteps=opt.nsteps,
        first_step=opt.first_step,
        min_step=opt.min_step,
        max_step=opt.max_step,
    )
    r.set_initial_value(initial_vector, tlist[0])
    code = compile("r.set_f_params(" + string + ")", "<string>", "exec")
    exec(code)

    #
    # call generic ODE code
    #
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)

def _mesolve_func_td(L_func, rho0, tlist, c_op_list, e_ops, args, opt, progress_bar):
    """
    Evolve the density matrix using an ODE solver with time dependent
    Hamiltonian.
    """

    if debug:
        print(inspect.stack()[0][3])

    #
    # check initial state
    #
    if isket(rho0):
        rho0 = ket2dm(rho0)

    #
    # construct liouvillian
    #
    new_args = None

    if len(c_op_list) > 0:
        L_data = liouvillian(None, c_op_list).data
    else:
        n, m = rho0.shape
        L_data = sp.csr_matrix((n ** 2, m ** 2), dtype=complex)

    if type(args) is dict:
        new_args = {}
        for key in args:
            if isinstance(args[key], Qobj):
                if isoper(args[key]):
                    new_args[key] = (-1j * (spre(args[key]) - spost(args[key]))).data
                else:
                    new_args[key] = args[key].data
            else:
                new_args[key] = args[key]

    elif type(args) is list or type(args) is tuple:
        new_args = []
        for arg in args:
            if isinstance(arg, Qobj):
                if isoper(arg):
                    new_args.append((-1j * (spre(arg) - spost(arg))).data)
                else:
                    new_args.append(arg.data)
            else:
                new_args.append(arg)

        if type(args) is tuple:
            new_args = tuple(new_args)
    else:
        if isinstance(args, Qobj):
            if isoper(args):
                new_args = (-1j * (spre(args) - spost(args))).data
            else:
                new_args = args.data
        else:
            new_args = args

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full()).ravel()
    if not opt.rhs_with_state:
        r = scipy.integrate.ode(cy_ode_rho_func_td)
    else:
        r = scipy.integrate.ode(_ode_rho_func_td_with_state)
    r.set_integrator(
        "zvode",
        method=opt.method,
        order=opt.order,
        atol=opt.atol,
        rtol=opt.rtol,
        nsteps=opt.nsteps,
        first_step=opt.first_step,
        min_step=opt.min_step,
        max_step=opt.max_step,
    )
    r.set_initial_value(initial_vector, tlist[0])
    r.set_f_params(L_data, L_func, new_args)

    #
    # call generic ODE code
    #
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)

#.........这里部分代码省略.........
                Lconst += c
            else:
                raise TypeError(
                    "Incorrect specification of time-dependent "
                    + "Liouvillian (expected operator or "
                    + "superoperator)"
                )

        elif isinstance(c_spec, list):
            c = c_spec[0]
            c_coeff = c_spec[1]

            if isoper(c):
                cdc = c.dag() * c
                L = spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc)
                c_coeff = "(" + c_coeff + ")**2"
            elif issuper(c):
                L = c
            else:
                raise TypeError(
                    "Incorrect specification of time-dependent "
                    + "Liouvillian (expected operator or "
                    + "superoperator)"
                )

            Ldata.append(L.data.data)
            Linds.append(L.data.indices)
            Lptrs.append(L.data.indptr)
            Lcoeff.append(c_coeff)

        else:
            raise TypeError(
                "Incorrect specification of time-dependent " + "collapse operators (expected string format)"
            )

    # add the constant part of the lagrangian
    if Lconst != 0:
        Ldata.append(Lconst.data.data)
        Linds.append(Lconst.data.indices)
        Lptrs.append(Lconst.data.indptr)
        Lcoeff.append("1.0")

    # the total number of liouvillian terms (hamiltonian terms +
    # collapse operators)
    n_L_terms = len(Ldata)

    #
    # setup ode args string: we expand the list Ldata, Linds and Lptrs into
    # and explicit list of parameters
    #
    string_list = []
    for k in range(n_L_terms):
        string_list.append("Ldata[%d], Linds[%d], Lptrs[%d]" % (k, k, k))
    for name, value in args.items():
        if isinstance(value, np.ndarray):
            string_list.append(name)
        else:
            string_list.append(str(value))
    parameter_string = ",".join(string_list)

    #
    # generate and compile new cython code if necessary
    #
    if not opt.rhs_reuse or config.tdfunc is None:
        if opt.rhs_filename is None:
            config.tdname = "rhs" + str(os.getpid()) + str(config.cgen_num)
        else:
            config.tdname = opt.rhs_filename
        cgen = Codegen(h_terms=n_L_terms, h_tdterms=Lcoeff, args=args, config=config)
        cgen.generate(config.tdname + ".pyx")

        code = compile("from " + config.tdname + " import cy_td_ode_rhs", "<string>", "exec")
        exec(code, globals())
        config.tdfunc = cy_td_ode_rhs

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full()).ravel()
    r = scipy.integrate.ode(config.tdfunc)
    r.set_integrator(
        "zvode",
        method=opt.method,
        order=opt.order,
        atol=opt.atol,
        rtol=opt.rtol,
        nsteps=opt.nsteps,
        first_step=opt.first_step,
        min_step=opt.min_step,
        max_step=opt.max_step,
    )
    r.set_initial_value(initial_vector, tlist[0])
    code = compile("r.set_f_params(" + parameter_string + ")", "<string>", "exec")

    exec(code, locals(), args)

    #
    # call generic ODE code
    #
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)

def _mesolve_list_func_td(H_list, rho0, tlist, c_list, e_ops, args, opt, progress_bar):
    """
    Internal function for solving the master equation. See mesolve for usage.
    """

    if debug:
        print(inspect.stack()[0][3])

    #
    # check initial state
    #
    if isket(rho0):
        rho0 = rho0 * rho0.dag()

    #
    # construct liouvillian in list-function format
    #
    L_list = []
    if opt.rhs_with_state:
        constant_func = lambda x, y, z: 1.0
    else:
        constant_func = lambda x, y: 1.0

    # add all hamitonian terms to the lagrangian list
    for h_spec in H_list:

        if isinstance(h_spec, Qobj):
            h = h_spec
            h_coeff = constant_func

        elif isinstance(h_spec, list) and isinstance(h_spec[0], Qobj):
            h = h_spec[0]
            h_coeff = h_spec[1]

        else:
            raise TypeError("Incorrect specification of time-dependent " + "Hamiltonian (expected callback function)")

        if isoper(h):
            L_list.append([(-1j * (spre(h) - spost(h))).data, h_coeff, False])

        elif issuper(h):
            L_list.append([h.data, h_coeff, False])

        else:
            raise TypeError(
                "Incorrect specification of time-dependent " + "Hamiltonian (expected operator or superoperator)"
            )

    # add all collapse operators to the liouvillian list
    for c_spec in c_list:

        if isinstance(c_spec, Qobj):
            c = c_spec
            c_coeff = constant_func
            c_square = False

        elif isinstance(c_spec, list) and isinstance(c_spec[0], Qobj):
            c = c_spec[0]
            c_coeff = c_spec[1]
            c_square = True

        else:
            raise TypeError(
                "Incorrect specification of time-dependent " + "collapse operators (expected callback function)"
            )

        if isoper(c):
            L_list.append([liouvillian(None, [c], data_only=True), c_coeff, c_square])

        elif issuper(c):
            L_list.append([c.data, c_coeff, c_square])

        else:
            raise TypeError(
                "Incorrect specification of time-dependent "
                + "collapse operators (expected operator or "
                + "superoperator)"
            )

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full()).ravel()
    if opt.rhs_with_state:
        r = scipy.integrate.ode(drho_list_td_with_state)
    else:
        r = scipy.integrate.ode(drho_list_td)
    r.set_integrator(
        "zvode",
        method=opt.method,
        order=opt.order,
        atol=opt.atol,
        rtol=opt.rtol,
        nsteps=opt.nsteps,
        first_step=opt.first_step,
        min_step=opt.min_step,
        max_step=opt.max_step,
    )
    r.set_initial_value(initial_vector, tlist[0])
    r.set_f_params(L_list, args)
#.........这里部分代码省略.........