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 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 _generate_rho_A_ops(sc, L, dt):
    """
    pre-compute superoperator operator combinations that are commonly needed
    when evaluating the RHS of stochastic master equations
    """
    out = []
    for c_idx, c in enumerate(sc):
        n = c.dag() * c
        out.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)])

    return out

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 _generate_A_ops_Euler(sc, L, dt):
	"""
	combine precomputed operators in one long operator for the Euler method
	"""
	A_len = len(sc)
	out = []
	out += [spre(c).data + spost(c.dag()).data for c in sc]
	out += [(L + np.sum([lindblad_dissipator(c, data_only=True) for c in sc], axis=0))*dt]
	out1 = [[sp.vstack(out).tocsr(), sc[0].shape[0]]]
	#the following hack is required for compatibility with old A_ops
	out1 += [[] for n in xrange(A_len-1)]
	return out1

def _generate_A_ops_Milstein(sc, L, dt):
	"""
	combine precomputed operators in one long operator for the Milstein method
	with commuting stochastic jump operators.
	"""
	A_len = len(sc)
	temp = [spre(c).data + spost(c.dag()).data for c in sc]
	out = []
	out += temp
	out += [temp[n]*temp[n] for n in xrange(A_len)]
	out += [temp[n]*temp[m] for (n,m) in np.ndindex(A_len,A_len) if n > m]
	out += [(L + np.sum([lindblad_dissipator(c, data_only=True) for c in sc], axis=0))*dt]
	out1 = [[sp.vstack(out).tocsr(), sc[0].shape[0]]]
	#the following hack is required for compatibility with old A_ops
	out1 += [[] for n in xrange(A_len-1)]
	return out1

def to_super(q_oper):
    """
    Converts a Qobj representing a quantum map to the supermatrix (Liouville)
    representation.

    Parameters
    ----------
    q_oper : Qobj
        Superoperator to be converted to supermatrix representation. If
        ``q_oper`` is ``type="oper"``, then it is taken to act by conjugation,
        such that ``to_super(A) == sprepost(A, A.dag())``.

    Returns
    -------
    superop : Qobj
        A quantum object representing the same map as ``q_oper``, such that
        ``superop.superrep == "super"``.

    Raises
    ------
    TypeError
        If the given quantum object is not a map, or cannot be converted
        to supermatrix representation.
    """
    if q_oper.type == 'super':
        # Case 1: Already done.
        if q_oper.superrep == "super":
            return q_oper
        # Case 2: Can directly convert.
        elif q_oper.superrep == 'choi':
            return choi_to_super(q_oper)
        # Case 3: Need to go through Choi.
        elif q_oper.superrep == 'chi':
            return to_super(to_choi(q_oper))
        # Case 4: Something went wrong.
        else:
            raise ValueError(
                "Unrecognized superrep '{}'.".format(q_oper.superrep))
    elif q_oper.type == 'oper':  # Assume unitary
        return spre(q_oper) * spost(q_oper.dag())
    else:
        raise TypeError(
            "Conversion of Qobj with type = {0.type} "
            "and superrep = {0.superrep} to supermatrix not "
            "supported.".format(q_oper)
        )

def test_SuperType():
    "Qobj superoperator type"

    psi = basis(2, 1)
    rho = psi * psi.dag()

    sop = spre(rho)

    assert_equal(sop.isket, False)
    assert_equal(sop.isbra, False)
    assert_equal(sop.isoper, False)
    assert_equal(sop.issuper, True)

    sop = spost(rho)

    assert_equal(sop.isket, False)
    assert_equal(sop.isbra, False)
    assert_equal(sop.isoper, False)
    assert_equal(sop.issuper, True)

def to_chi(q_oper):
    """
    Converts a Qobj representing a quantum map to a representation as a chi
    (process) matrix in the Pauli basis, such that the trace of the returned
    operator is equal to the dimension of the system.

    Parameters
    ----------
    q_oper : Qobj
        Superoperator to be converted to Chi representation. If
        ``q_oper`` is ``type="oper"``, then it is taken to act by conjugation,
        such that ``to_chi(A) == to_chi(sprepost(A, A.dag()))``.

    Returns
    -------
    chi : Qobj
        A quantum object representing the same map as ``q_oper``, such that
        ``chi.superrep == "chi"``.

    Raises
    ------
    TypeError: if the given quantum object is not a map, or cannot be converted
        to Chi representation.
    """
    if q_oper.type == 'super':
        # Case 1: Already done.
        if q_oper.superrep == 'chi':
            return q_oper
        # Case 2: Can directly convert.
        elif q_oper.superrep == 'choi':
            return choi_to_chi(q_oper)
        # Case 3: Need to go through Choi.
        elif q_oper.superrep == 'super':
            return to_chi(to_choi(q_oper))
        else:
            raise TypeError(q_oper.superrep)
    elif q_oper.type == 'oper':
        return to_chi(spre(q_oper) * spost(q_oper.dag()))
    else:
        raise TypeError(
            "Conversion of Qobj with type = {0.type} "
            "and superrep = {0.choi} to Choi not supported.".format(q_oper)
        )

def to_choi(q_oper):
    """
    Converts a Qobj representing a quantum map to the Choi representation,
    such that the trace of the returned operator is equal to the dimension
    of the system.

    Parameters
    ----------
    q_oper : Qobj
        Superoperator to be converted to Choi representation. If
        ``q_oper`` is ``type="oper"``, then it is taken to act by conjugation,
        such that ``to_choi(A) == to_choi(sprepost(A, A.dag()))``.

    Returns
    -------
    choi : Qobj
        A quantum object representing the same map as ``q_oper``, such that
        ``choi.superrep == "choi"``.

    Raises
    ------
    TypeError: if the given quantum object is not a map, or cannot be converted
        to Choi representation.
    """
    if q_oper.type == 'super':
        if q_oper.superrep == 'choi':
            return q_oper
        if q_oper.superrep == 'super':
            return super_to_choi(q_oper)
        if q_oper.superrep == 'chi':
            return chi_to_choi(q_oper)
        else:
            raise TypeError(q_oper.superrep)
    elif q_oper.type == 'oper':
        return super_to_choi(spre(q_oper) * spost(q_oper.dag()))
    else:
        raise TypeError(
            "Conversion of Qobj with type = {0.type} "
            "and superrep = {0.choi} to Choi not supported.".format(q_oper)
        )

def rhs_generate(H,c_ops,args={},options=Odeoptions(),name=None):
    """
    Generates the Cython functions needed for solving the dynamics of a
    given system using the mesolve function inside a parfor loop.  
    
    Parameters
    ----------
    H : qobj
        System Hamiltonian.
    c_ops : list
        ``list`` of collapse operators.
    args : dict
        Arguments for time-dependent Hamiltonian and collapse operator terms.
    options : Odeoptions
        Instance of ODE solver options.
    name: str
        Name of generated RHS
    
    Notes
    -----
    Using this function with any solver other than the mesolve function
    will result in an error.
    
    """
    _reset_odeconfig() #clear odeconfig
    if name:
        odeconfig.tdname=name
    else:
        odeconfig.tdname="rhs"+str(odeconfig.cgen_num)
    
    n_op = len(c_ops)

    Lconst = 0        

    Ldata = []
    Linds = []
    Lptrs = []
    Lcoeff = []
    
    # loop over all hamiltonian terms, convert to superoperator form and 
    # add the data of sparse matrix represenation to 
    for h_spec in H:
        if isinstance(h_spec, Qobj):
            h = h_spec
            Lconst += -1j*(spre(h) - spost(h)) 
        
        elif isinstance(h_spec, list): 
            h = h_spec[0]
            h_coeff = h_spec[1]

            L = -1j*(spre(h) - spost(h))

            Ldata.append(L.data.data)
            Linds.append(L.data.indices)
            Lptrs.append(L.data.indptr)
            Lcoeff.append(h_coeff)
            
        else:
            raise TypeError("Incorrect specification of time-dependent " + 
                             "Hamiltonian (expected string format)")
    
    # loop over all collapse operators        
    for c_spec in c_ops:
        if isinstance(c_spec, Qobj):
            c = c_spec
            cdc = c.dag() * c
            Lconst += spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc) 

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

            cdc = c.dag() * c
            L = spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc) 

            Ldata.append(L.data.data)
            Linds.append(L.data.indices)
            Lptrs.append(L.data.indptr)
            Lcoeff.append("("+c_coeff+")**2")

        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)
    
    cgen=Codegen(h_terms=n_L_terms,h_tdterms=Lcoeff, args=args)
    cgen.generate(odeconfig.tdname+".pyx")
    os.environ['CFLAGS'] = '-O3 -w'
    import pyximport
    pyximport.install(setup_args={'include_dirs':[numpy.get_include()]})
#.........这里部分代码省略.........

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)

def _mesolve_list_str_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: must be a density matrix
    #
    if isket(rho0):
        rho0 = rho0 * rho0.dag()

    #
    # construct liouvillian
    #
    Lconst = 0

    Ldata = []
    Linds = []
    Lptrs = []
    Lcoeff = []

    # loop over all hamiltonian terms, convert to superoperator form and
    # add the data of sparse matrix representation to
    for h_spec in H_list:

        if isinstance(h_spec, Qobj):
            h = h_spec

            if isoper(h):
                Lconst += -1j * (spre(h) - spost(h))
            elif issuper(h):
                Lconst += h
            else:
                raise TypeError(
                    "Incorrect specification of time-dependent "
                    + "Hamiltonian (expected operator or "
                    + "superoperator)"
                )

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

            if isoper(h):
                L = -1j * (spre(h) - spost(h))
            elif issuper(h):
                L = h
            else:
                raise TypeError(
                    "Incorrect specification of time-dependent "
                    + "Hamiltonian (expected operator or "
                    + "superoperator)"
                )

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

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

    # loop over all collapse operators
    for c_spec in c_list:

        if isinstance(c_spec, Qobj):
            c = c_spec

            if isoper(c):
                cdc = c.dag() * c
                Lconst += spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc)
            elif issuper(c):
                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)"
                )

#.........这里部分代码省略.........

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)
#.........这里部分代码省略.........

def rhs_generate(H, c_ops, args={}, options=Options(), name=None,
                 cleanup=True):
    """
    Generates the Cython functions needed for solving the dynamics of a
    given system using the mesolve function inside a parfor loop.

    Parameters
    ----------
    H : qobj
        System Hamiltonian.

    c_ops : list
        ``list`` of collapse operators.

    args : dict
        Arguments for time-dependent Hamiltonian and collapse operator terms.

    options : Options
        Instance of ODE solver options.

    name: str
        Name of generated RHS

    cleanup: bool
        Whether the generated cython file should be automatically removed or
        not.

    Notes
    -----
    Using this function with any solver other than the mesolve function
    will result in an error.

    """
    config.reset()
    config.options = options

    if name:
        config.tdname = name
    else:
        config.tdname = "rhs" + str(os.getpid()) + str(config.cgen_num)

    Lconst = 0

    Ldata = []
    Linds = []
    Lptrs = []
    Lcoeff = []

    # loop over all hamiltonian terms, convert to superoperator form and
    # add the data of sparse matrix represenation to

    msg = "Incorrect specification of time-dependence: "

    for h_spec in H:
        if isinstance(h_spec, Qobj):
            h = h_spec

            if not isinstance(h, Qobj):
                raise TypeError(msg + "expected Qobj")

            if h.isoper:
                Lconst += -1j * (spre(h) - spost(h))
            elif h.issuper:
                Lconst += h
            else:
                raise TypeError(msg + "expected operator or superoperator")

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

            if not isinstance(h, Qobj):
                raise TypeError(msg + "expected Qobj")

            if h.isoper:
                L = -1j * (spre(h) - spost(h))
            elif h.issuper:
                L = h
            else:
                raise TypeError(msg + "expected operator or superoperator")

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

        else:
            raise TypeError(msg + "expected string format")

    # loop over all collapse operators
    for c_spec in c_ops:
        if isinstance(c_spec, Qobj):
            c = c_spec

            if not isinstance(c, Qobj):
                raise TypeError(msg + "expected Qobj")

            if c.isoper:
                cdc = c.dag() * c
                Lconst += spre(c) * spost(c.dag()) - 0.5 * spre(cdc) \
#.........这里部分代码省略.........

def _mesolve_list_str_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: must be a density matrix
    #
    if isket(rho0):
        rho0 = rho0 * rho0.dag()

    #
    # construct liouvillian
    #
    Lconst = 0

    Ldata = []
    Linds = []
    Lptrs = []
    Lcoeff = []
    Lobj = []
    me_cops_coeff = []
    me_cops_obj = []
    me_cops_obj_flags = []

    # loop over all hamiltonian terms, convert to superoperator form and
    # add the data of sparse matrix representation to
    n_not_const_terms = 0
    for h_spec in H_list:
        if isinstance(h_spec, Qobj):
            h = h_spec

            if isoper(h):
                Lconst += -1j * (spre(h) - spost(h))
            elif issuper(h):
                Lconst += h
            else:
                raise TypeError("Incorrect specification of time-dependent " +
                                "Hamiltonian (expected operator or " +
                                "superoperator)")

        elif isinstance(h_spec, list):
            n_not_const_terms +=1
            h = h_spec[0]
            h_coeff = h_spec[1]

            if isoper(h):
                L = -1j * (spre(h) - spost(h))
            elif issuper(h):
                L = h
            else:
                raise TypeError("Incorrect specification of time-dependent " +
                                "Hamiltonian (expected operator or " +
                                "superoperator)")

            Ldata.append(L.data.data)
            Linds.append(L.data.indices)
            Lptrs.append(L.data.indptr)
            if isinstance(h_coeff, Cubic_Spline):
                Lobj.append(h_coeff.coeffs)
            Lcoeff.append(h_coeff)

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


    
    # loop over all collapse operators
    for c_spec in c_list:
        if isinstance(c_spec, Qobj):
            c = c_spec

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

        elif isinstance(c_spec, list):
            n_not_const_terms +=1
            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)
                if isinstance(c_coeff, Cubic_Spline):
                    me_cops_obj.append(c_coeff.coeffs)
                    me_cops_obj_flags.append(n_not_const_terms)
#.........这里部分代码省略.........

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 configure(self, H_sys, coup_op, coup_strength, temperature,
                     N_cut, N_exp, cut_freq, planck=None, boltzmann=None,
                     renorm=None, bnd_cut_approx=None,
                     options=None, progress_bar=None, stats=None):
        """
        Calls configure from :class:`HEOMSolver` and sets any attributes
        that are specific to this subclass
        """
        start_config = timeit.default_timer()

        HEOMSolver.configure(self, H_sys, coup_op, coup_strength,
                    temperature, N_cut, N_exp,
                    planck=planck, boltzmann=boltzmann,
                    options=options, progress_bar=progress_bar, stats=stats)
        self.cut_freq = cut_freq
        if renorm is not None: self.renorm = renorm
        if bnd_cut_approx is not None: self.bnd_cut_approx = bnd_cut_approx

        # Load local values for optional parameters
        # Constants and Hamiltonian.
        hbar = self.planck
        options = self.options
        progress_bar = self.progress_bar
        stats = self.stats


        if stats:
            ss_conf = stats.sections.get('config')
            if ss_conf is None:
                ss_conf = stats.add_section('config')

        c, nu = self._calc_matsubara_params()

        if renorm:
            norm_plus, norm_minus = self._calc_renorm_factors()
            if stats:
                stats.add_message('options', 'renormalisation', ss_conf)
        # Dimensions et by system
        sup_dim = H_sys.dims[0][0]**2
        unit_sys = qeye(H_sys.dims[0])

        # Use shorthands (mainly as in referenced PRL)
        lam0 = self.coup_strength
        gam = self.cut_freq
        N_c = self.N_cut
        N_m = self.N_exp
        Q = coup_op # Q as shorthand for coupling operator
        beta = 1.0/(self.boltzmann*self.temperature)

        # Ntot is the total number of ancillary elements in the hierarchy
        # Ntot = factorial(N_c + N_m) / (factorial(N_c)*factorial(N_m))
        # Turns out to be the same as nstates from state_number_enumerate
        N_he, he2idx, idx2he = enr_state_dictionaries([N_c + 1]*N_m , N_c)

        unit_helems = sp.identity(N_he, format='csr')
        if self.bnd_cut_approx:
            # the Tanimura boundary cut off operator
            if stats:
                stats.add_message('options', 'boundary cutoff approx', ss_conf)
            op = -2*spre(Q)*spost(Q.dag()) + spre(Q.dag()*Q) + spost(Q.dag()*Q)

            approx_factr = ((2*lam0 / (beta*gam*hbar)) - 1j*lam0) / hbar
            for k in range(N_m):
                approx_factr -= (c[k] / nu[k])
            L_bnd = -approx_factr*op.data
            L_helems = sp.kron(unit_helems, L_bnd)
        else:
            L_helems = sp.csr_matrix((N_he*sup_dim, N_he*sup_dim),
                                     dtype=complex)

        # Build the hierarchy element interaction matrix
        if stats: start_helem_constr = timeit.default_timer()

        unit_sup = spre(unit_sys).data
        spreQ = spre(Q).data
        spostQ = spost(Q).data
        commQ = (spre(Q) - spost(Q)).data
        N_he_interact = 0

        for he_idx in range(N_he):
            he_state = list(idx2he[he_idx])
            n_excite = sum(he_state)

            # The diagonal elements for the hierarchy operator
            # coeff for diagonal elements
            sum_n_m_freq = 0.0
            for k in range(N_m):
                sum_n_m_freq += he_state[k]*nu[k]

            op = -sum_n_m_freq*unit_sup
            L_he = _pad_csr(op, N_he, N_he, he_idx, he_idx)
            L_helems += L_he

            # Add the neighour interations
            he_state_neigh = copy(he_state)
            for k in range(N_m):

                n_k = he_state[k]
                if n_k >= 1:
                    # find the hierarchy element index of the neighbour before
#.........这里部分代码省略.........

def bloch_redfield_tensor(H, a_ops, spectra_cb, use_secular=True):
    """
    Calculate the Bloch-Redfield tensor for a system given a set of operators
    and corresponding spectral functions that describes the system's coupling
    to its environment.

    Parameters
    ----------

    H : :class:`qutip.qobj`
        System Hamiltonian.

    a_ops : list of :class:`qutip.qobj`
        List of system operators that couple to the environment.

    spectra_cb : list of callback functions
        List of callback functions that evaluate the noise power spectrum
        at a given frequency.

    use_secular : bool
        Flag (True of False) that indicates if the secular approximation should
        be used.

    Returns
    -------

    R, kets: :class:`qutip.qobj`, list of :class:`qutip.qobj`

        R is the Bloch-Redfield tensor and kets is a list eigenstates of the
        Hamiltonian.

    """

    # Sanity checks for input parameters
    if not isinstance(H, Qobj):
        raise TypeError("H must be an instance of Qobj")

    for a in a_ops:
        if not isinstance(a, Qobj) or not a.isherm:
            raise TypeError("Operators in a_ops must be Hermitian Qobj.")

    # default spectrum
    if not spectra_cb:
        spectra_cb = [lambda w: 1.0 for _ in a_ops]

    # use the eigenbasis
    evals, ekets = H.eigenstates()

    N = len(evals)
    K = len(a_ops)
    A = np.zeros((K, N, N), dtype=complex)  # TODO: use sparse here
    W = np.zeros((N, N))

    # pre-calculate matrix elements
    for n in range(N):
        for m in range(N):
            W[m, n] = np.real(evals[m] - evals[n])

    for k in range(K):
        # A[k,n,m] = a_ops[k].matrix_element(ekets[n], ekets[m])
        A[k, :, :] = a_ops[k].transform(ekets).full()

    dw_min = abs(W[W.nonzero()]).min()

    # unitary part
    Heb = H.transform(ekets)
    R = -1.0j * (spre(Heb) - spost(Heb))
    R.data = R.data.tolil()
    for I in range(N * N):
        a, b = vec2mat_index(N, I)
        for J in range(N * N):
            c, d = vec2mat_index(N, J)

            # unitary part: use spre and spost above, same as this:
            # R.data[I,J] = -1j * W[a,b] * (a == c) * (b == d)

            if use_secular is False or abs(W[a, b] - W[c, d]) < dw_min / 10.0:

                # dissipative part:
                for k in range(K):
                    # for each operator coupling the system to the environment