def _correlation_me_2op_2t(H, rho0, tlist, taulist, c_ops, a_op, b_op,
                           reverse=False, args=None, options=Odeoptions()):
    """
    Internal function for calculating correlation functions using the master
    equation solver. See :func:`correlation` for usage.
    """

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

    if rho0 is None:
        rho0 = steadystate(H, c_ops)
    elif rho0 and isket(rho0):
        rho0 = ket2dm(rho0)

    C_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)

    rho_t_list = mesolve(
        H, rho0, tlist, c_ops, [], args=args, options=options).states

    if reverse:
        # <A(t)B(t+tau)>
        for t_idx, rho_t in enumerate(rho_t_list):
            C_mat[t_idx, :] = mesolve(H, rho_t * a_op, taulist,
                                      c_ops, [b_op], args=args,
                                      options=options).expect[0]
    else:
        # <A(t+tau)B(t)>
        for t_idx, rho_t in enumerate(rho_t_list):
            C_mat[t_idx, :] = mesolve(H, b_op * rho_t, taulist,
                                      c_ops, [a_op], args=args,
                                      options=options).expect[0]

    return C_mat

def _correlation_me_4op_2t(H, rho0, tlist, taulist, c_ops,
                           a_op, b_op, c_op, d_op, reverse=False,
                           args=None, options=Odeoptions()):
    """
    Calculate the four-operator two-time correlation function on the form
    <A(t)B(t+tau)C(t+tau)D(t)>.

    See, Gardiner, Quantum Noise, Section 5.2.1
    """

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

    if rho0 is None:
        rho0 = steadystate(H, c_ops)
    elif rho0 and isket(rho0):
        rho0 = ket2dm(rho0)

    C_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)

    rho_t = mesolve(
        H, rho0, tlist, c_ops, [], args=args, options=options).states

    for t_idx, rho in enumerate(rho_t):
        C_mat[t_idx, :] = mesolve(H, d_op * rho * a_op, taulist,
                                  c_ops, [b_op * c_op],
                                  args=args, options=options).expect[0]

    return C_mat

def _correlation_es_2op_1t(H, rho0, tlist, c_ops, a_op, b_op, reverse=False,
                           args=None, options=Odeoptions()):
    """
    Internal function for calculating correlation functions using the
    exponential series solver. See :func:`correlation_ss` usage.
    """

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

    # contruct the Liouvillian
    L = liouvillian(H, c_ops)

    # find the steady state
    if rho0 is None:
        rho0 = steadystate(L)
    elif rho0 and isket(rho0):
        rho0 = ket2dm(rho0)

    # evaluate the correlation function
    if reverse:
        # <A(t)B(t+tau)>
        solC_tau = ode2es(L, rho0 * a_op)
        return esval(expect(b_op, solC_tau), tlist)
    else:
        # default: <A(t+tau)B(t)>
        solC_tau = ode2es(L, b_op * rho0)
        return esval(expect(a_op, solC_tau), tlist)

def _spectrum_es(H, wlist, c_ops, a_op, b_op):
    """
    Internal function for calculating the spectrum of the correlation function
    :math:`\left<A(\\tau)B(0)\\right>`.
    """
    if debug:
        print(inspect.stack()[0][3])

    # construct the Liouvillian
    L = liouvillian(H, c_ops)

    # find the steady state density matrix and a_op and b_op expecation values
    rho0 = steadystate(L)

    a_op_ss = expect(a_op, rho0)
    b_op_ss = expect(b_op, rho0)

    # eseries solution for (b * rho0)(t)
    es = ode2es(L, b_op * rho0)

    # correlation
    corr_es = expect(a_op, es)

    # covariance
    cov_es = corr_es - np.real(np.conjugate(a_op_ss) * b_op_ss)

    # spectrum
    spectrum = esspec(cov_es, wlist)

    return spectrum

def _correlation_es_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op):
    """
    Internal function for calculating the three-operator two-time
    correlation function:
    <A(t)B(t+tau)C(t)>
    using an exponential series solver.
    """

    # the solvers only work for positive time differences and the correlators
    # require positive tau
    if state0 is None:
        rho0 = steadystate(H, c_ops)
        tlist = [0]
    elif isket(state0):
        rho0 = ket2dm(state0)
    else:
        rho0 = state0

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

    # contruct the Liouvillian
    L = liouvillian(H, c_ops)

    corr_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
    solES_t = ode2es(L, rho0)

    # evaluate the correlation function
    for t_idx in range(len(tlist)):
        rho_t = esval(solES_t, [tlist[t_idx]])
        solES_tau = ode2es(L, c_op * rho_t * a_op)
        corr_mat[t_idx, :] = esval(expect(b_op, solES_tau), taulist)

    return corr_mat

def _spectrum_es(H, wlist, c_ops, a_op, b_op):
    """
    Internal function for calculating the spectrum of the correlation function
    :math:`\left<A(\\tau)B(0)\\right>`.
    """
    if debug:
        print(inspect.stack()[0][3])

    # construct the Liouvillian
    L = liouvillian(H, c_ops)

    # find the steady state density matrix and a_op and b_op expecation values
    rho0 = steadystate(L)

    a_op_ss = expect(a_op, rho0)
    b_op_ss = expect(b_op, rho0)

    # eseries solution for (b * rho0)(t)
    es = ode2es(L, b_op * rho0)

    # correlation
    corr_es = expect(a_op, es)

    # covariance
    cov_es = corr_es - a_op_ss * b_op_ss
    # tidy up covariance (to combine, e.g., zero-frequency components that cancel)
    cov_es.tidyup()

    # spectrum
    spectrum = esspec(cov_es, wlist)

    return spectrum

def coherence_function_g1(H, taulist, c_ops, a_op, solver="me", args=None,
                          options=Options(ntraj=[20, 100])):
    """
    Calculate the normalized first-order quantum coherence function:

    .. math::

        g^{(1)}(\\tau) = \lim_{t \to \infty}
        \\frac{\\langle a^\\dagger(t+\\tau)a(t)\\rangle}
        {\\langle a^\\dagger(t)a(t)\\rangle}

    using the quantum regression theorem and the evolution solver indicated by
    the `solver` parameter. Note: g1 is only defined for stationary
    statistics (uses steady state).

    Parameters
    ----------

    H : :class:`qutip.qobj.Qobj`
        system Hamiltonian.

    taulist : *list* / *array*
        list of times for :math:`\\tau`. taulist must be positive and contain
        the element `0`.

    c_ops : list of :class:`qutip.qobj.Qobj`
        list of collapse operators.

    a_op : :class:`qutip.qobj.Qobj`
        The annihilation operator of the mode.

    solver : str
        choice of solver (`me` for master-equation and
        `es` for exponential series)

    options : :class:`qutip.solver.Options`
        solver options class. `ntraj` is taken as a two-element list because
        the `mc` correlator calls `mcsolve()` recursively; by default,
        `ntraj=[20, 100]`. `mc_corr_eps` prevents divide-by-zero errors in
        the `mc` correlator; by default, `mc_corr_eps=1e-10`.

    Returns
    -------

    g1: *array*
        The normalized first-order coherence function.

    """

    # first calculate the steady state photon number
    rho0 = steadystate(H, c_ops)
    n = np.array([expect(rho0, a_op.dag() * a_op)])

    # calculate the correlation function G1 and normalize with n to obtain g1
    G1 = correlation_2op_1t(H, None, taulist, c_ops, a_op.dag(), a_op,
                            args=args, solver=solver, options=options)
    g1 = G1 / n

    return g1

def floquet_master_equation_steadystate(H, A):
    """
    Returns the steadystate density matrix (in the floquet basis!) for the
    Floquet-Markov master equation.
    """
    c_ops = floquet_collapse_operators(A)
    rho_ss = steadystate(H, c_ops)
    return rho_ss

def coherence_function_g2(H, rho0, taulist, c_ops, a_op, solver="me",
                          args=None, options=Odeoptions()):
    """
    Calculate the second-order quantum coherence function:

    .. math::

        g^{(2)}(\\tau) =
        \\frac{\\langle a^\\dagger(0)a^\\dagger(\\tau)a(\\tau)a(0)\\rangle}
        {\\langle a^\\dagger(\\tau)a(\\tau)\\rangle
         \\langle a^\\dagger(0)a(0)\\rangle}

    Parameters
    ----------

    H : :class:`qutip.qobj.Qobj`
        system Hamiltonian.

    rho0 : :class:`qutip.qobj.Qobj`
        Initial state density matrix (or state vector). If 'rho0' is
        'None', then the steady state will be used as initial state.

    taulist : *list* / *array*
        list of times for :math:`\\tau`.

    c_ops : list of :class:`qutip.qobj.Qobj`
        list of collapse operators.

    a_op : :class:`qutip.qobj.Qobj`
        The annihilation operator of the mode.

    solver : str
        choice of solver (currently only 'me')

    Returns
    -------

    g2, G2: tuble of *array*
        The normalized and unnormalized second-order coherence function.

    """

    # first calculate the photon number
    if rho0 is None:
        rho0 = steadystate(H, c_ops)
        n = np.array([expect(rho0, a_op.dag() * a_op)])
    else:
        n = mesolve(
            H, rho0, taulist, c_ops, [a_op.dag() * a_op], 
            args=args, options=options).expect[0]

    # calculate the correlation function G2 and normalize with n to obtain g2
    G2 = correlation_4op_1t(H, rho0, taulist, c_ops,
                            a_op.dag(), a_op.dag(), a_op, a_op,
                            solver=solver, args=args, options=options)
    g2 = G2 / (n[0] * n)

    return g2, G2

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 _correlation_me_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op,
                       args={}, options=Options()):
    """
    Internal function for calculating the three-operator two-time
    correlation function:
    <A(t)B(t+tau)C(t)>
    using a master equation solver.
    """

    # the solvers only work for positive time differences and the correlators
    # require positive tau
    if state0 is None:
        rho0 = steadystate(H, c_ops)
        tlist = [0]
    elif isket(state0):
        rho0 = ket2dm(state0)
    else:
        rho0 = state0

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

    rho_t = mesolve(H, rho0, tlist, c_ops, [],
                    args=args, options=options).states
    corr_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
    H_shifted, c_ops_shifted, _args = _transform_L_t_shift(H, c_ops, args)
    if config.tdname:
        _cython_build_cleanup(config.tdname)
    rhs_clear()

    for t_idx, rho in enumerate(rho_t):
        if not isinstance(H, Qobj):
            _args["_t0"] = tlist[t_idx]

        corr_mat[t_idx, :] = mesolve(
            H_shifted, c_op * rho * a_op, taulist, c_ops_shifted,
            [b_op], args=_args, options=options
        ).expect[0]

        if t_idx == 1:
            options.rhs_reuse = True

    if config.tdname:
        _cython_build_cleanup(config.tdname)
    rhs_clear()

    return corr_mat

def _correlation_me_4op_1t(H, rho0, tlist, c_ops, a_op, b_op, c_op, d_op, args=None, options=Options()):
    """
    Calculate the four-operator two-time correlation function on the form
    <A(0)B(tau)C(tau)D(0)>.

    See, Gardiner, Quantum Noise, Section 5.2.1
    """

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

    if rho0 is None:
        rho0 = steadystate(H, c_ops)
    elif rho0 and isket(rho0):
        rho0 = ket2dm(rho0)

    return mesolve(H, d_op * rho0 * a_op, tlist, c_ops, [b_op * c_op], args=args, options=options).expect[0]

def _correlation_me_2op_1t(H, rho0, tlist, c_ops, a_op, b_op, reverse=False, args=None, options=Options()):
    """
    Internal function for calculating correlation functions using the master
    equation solver. See :func:`correlation_ss` for usage.
    """

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

    if rho0 is None:
        rho0 = steadystate(H, c_ops)
    elif rho0 and isket(rho0):
        rho0 = ket2dm(rho0)

    if reverse:
        # <A(t)B(t+tau)>
        return mesolve(H, rho0 * a_op, tlist, c_ops, [b_op], args=args, options=options).expect[0]
    else:
        # <A(t+tau)B(t)>
        return mesolve(H, b_op * rho0, tlist, c_ops, [a_op], args=args, options=options).expect[0]

def _correlation_es_2op_2t(H, rho0, tlist, taulist, c_ops, a_op, b_op,
                           reverse=False, args=None, options=Odeoptions()):
    """
    Internal function for calculating correlation functions using the
    exponential series solver. See :func:`correlation` usage.
    """

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

    # contruct the Liouvillian
    L = liouvillian(H, c_ops)

    if rho0 is None:
        rho0 = steadystate(L)
    elif rho0 and isket(rho0):
        rho0 = ket2dm(rho0)

    C_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)

    solES_t = ode2es(L, rho0)

    # evaluate the correlation function
    if reverse:
        # <A(t)B(t+tau)>
        for t_idx in range(len(tlist)):
            rho_t = esval(solES_t, [tlist[t_idx]])
            solES_tau = ode2es(L, rho_t * a_op)
            C_mat[t_idx, :] = esval(expect(b_op, solES_tau), taulist)

    else:
        # default: <A(t+tau)B(t)>
        for t_idx in range(len(tlist)):
            rho_t = esval(solES_t, [tlist[t_idx]])
            solES_tau = ode2es(L, b_op * rho_t)
            C_mat[t_idx, :] = esval(expect(a_op, solES_tau), taulist)

    return C_mat

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 countstat_current_noise(L, c_ops, rhoss=None, J_ops=None, R=False):
    """
    Compute the cross-current noise spectrum for a list of collapse operators
    `c_ops` corresponding to monitored currents, given the system
    Liouvillian `L`. The current collapse operators `c_ops` should be part
    of the dissipative processes in `L`, but the `c_ops` given here does not
    necessarily need to be all collapse operators contributing to dissipation
    in the Liouvillian. Optionally, the steadystate density matrix `rhoss`
    and/or the pseudo inverse `R` of the Liouvillian `L`, and the current
    operators `J_ops` correpsonding to the current collapse operators `c_ops`
    can also be specified. If `R` is not given, the cross-current correlations
    will be computed directly without computing `R` explicitly. If either of
    `rhoss` and `J_ops` are omitted, they will be computed internally.

    Parameters
    ----------

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

    c_ops : array / list
        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.

    R : :class:`qutip.Qobj` (optional)
        Qobj representing the pseudo inverse of the system Liouvillian `L`.

    Returns
    --------
    I, S : tuple of arrays
        The currents `I` corresponding to each current collapse operator
        `c_ops` (or, equivalently, each current superopeator `J_ops`) and the
        zero-frequency cross-current correlation `S`.
    """

    if rhoss is None:
        rhoss = steadystate(L, c_ops)

    if J_ops is None:
        J_ops = [sprepost(c, c.dag()) for c in c_ops]

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

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

    if R:
        if R is True:
            R = pseudo_inverse(L, rhoss)

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

                S[i, j] -= expect_rho_vec((Ji * R * Jj + Jj * R * Ji).data,
                                          rhoss_vec, 1)
    else:
        N = np.prod(L.dims[0][0])

        rhoss_vec = operator_to_vector(rhoss)

        tr_op = tensor([identity(n) for n in L.dims[0][0]])
        tr_op_vec = operator_to_vector(tr_op)

        Pop = sp.kron(rhoss_vec.data, tr_op_vec.data.T, format='csc')
        Iop = sp.eye(N*N, N*N, format='csc')
        Q = Iop - Pop
        A = L.data.tocsc()
        rhoss_vec = mat2vec(rhoss.full()).ravel()

        for j, Jj in enumerate(J_ops):
            Qj = Q * Jj.data * rhoss_vec
            X_rho_vec = sp.linalg.splu(A, permc_spec='COLAMD').solve(Qj)
            for i, Ji in enumerate(J_ops):
                if i == j:
                    S[i, i] = I[i] = expect_rho_vec(Ji.data, rhoss_vec, 1)

                S[i, j] -= expect_rho_vec(Ji.data * Q, X_rho_vec, 1)

    return I, S

def coherence_function_g2(H, state0, taulist, c_ops, a_op, solver="me", args={},
                          options=Options(ntraj=[20, 100])):
    """
    Calculate the normalized second-order quantum coherence function:

    .. math::

         g^{(2)}(\\tau) =
        \\frac{\\langle A^\\dagger(0)A^\\dagger(\\tau)A(\\tau)A(0)\\rangle}
        {\\langle A^\\dagger(\\tau)A(\\tau)\\rangle
         \\langle A^\\dagger(0)A(0)\\rangle}

    using the quantum regression theorem and the evolution solver indicated by
    the `solver` parameter.

    Parameters
    ----------
    H : Qobj
        system Hamiltonian, may be time-dependent for solver choice of `me` or
        `mc`.
    state0 : Qobj
        Initial state density matrix :math:`\\rho(t_0)` or state vector
        :math:`\\psi(t_0)`. If 'state0' is 'None', then the steady state will
        be used as the initial state. The 'steady-state' is only implemented
        for the `me` and `es` solvers.
    taulist : array_like
        list of times for :math:`\\tau`. taulist must be positive and contain
        the element `0`.
    c_ops : list
        list of collapse operators, may be time-dependent for solver choice of
        `me` or `mc`.
    a_op : Qobj
        operator A.
    solver : str
        choice of solver (`me` for master-equation and
        `es` for exponential series).
    options : Options
        solver options class. `ntraj` is taken as a two-element list because
        the `mc` correlator calls `mcsolve()` recursively; by default,
        `ntraj=[20, 100]`. `mc_corr_eps` prevents divide-by-zero errors in
        the `mc` correlator; by default, `mc_corr_eps=1e-10`.

    Returns
    -------
    g2, G2 : tuple
        The normalized and unnormalized second-order coherence function.

    """

    # first calculate the photon number
    if state0 is None:
        state0 = steadystate(H, c_ops)
        n = np.array([expect(state0, a_op.dag() * a_op)])
    else:
        n = mesolve(H, state0, taulist, c_ops, [a_op.dag() * a_op]).expect[0]

    # calculate the correlation function G2 and normalize with n to obtain g2
    G2 = correlation_3op_1t(H, state0, taulist, c_ops,
                            a_op.dag(), a_op.dag()*a_op, a_op,
                            solver=solver, args=args, options=options)
    g2 = G2 / (n[0] * n)

    return g2, G2

def countstat_current_noise(L, c_ops, wlist=None, rhoss=None, J_ops=None, 
                            sparse=True, method='direct'):
    """
    Compute the cross-current noise spectrum for a list of collapse operators
    `c_ops` corresponding to monitored currents, given the system
    Liouvillian `L`. The current collapse operators `c_ops` should be part
    of the dissipative processes in `L`, but the `c_ops` given here does not
    necessarily need to be all collapse operators contributing to dissipation
    in the Liouvillian. Optionally, the steadystate density matrix `rhoss`
    and the current operators `J_ops` correpsonding to the current collapse 
    operators `c_ops` can also be specified. If either of
    `rhoss` and `J_ops` are omitted, they will be computed internally.
    'wlist' is an optional list of frequencies at which to evaluate the noise 
    spectrum.  
    
    Note:
    The default method is a direct solution using dense matrices, as sparse 
    matrix methods fail for some examples of small systems.
    For larger systems it is reccomended to use the sparse solver
    with the direct method, as it avoids explicit calculation of the
    pseudo-inverse, as described in page 67 of "Electrons in nanostructures"
    C. Flindt, PhD Thesis, available online:
    http://orbit.dtu.dk/fedora/objects/orbit:82314/datastreams/file_4732600/content
    
    Parameters
    ----------

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

    c_ops : array / list
        List of current collapse operators.

    rhoss : :class:`qutip.Qobj` (optional)
        The steadystate density matrix corresponding the system Liouvillian
        `L`.
        
    wlist : array / list (optional)
        List of frequencies at which to evaluate (if none are given, evaluates 
        at zero frequency)

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

    sparse : bool
        Flag that indicates whether to use sparse or dense matrix methods when
        computing the pseudo inverse. Default is false, as sparse solvers
        can fail for small systems. For larger systems the sparse solvers
        are reccomended. 
        
        
    Returns
    --------
    I, S : tuple of arrays
        The currents `I` corresponding to each current collapse operator
        `c_ops` (or, equivalently, each current superopeator `J_ops`) and the
        zero-frequency cross-current correlation `S`.
    """

    if rhoss is None:
        rhoss = steadystate(L, c_ops)

    if J_ops is None:
        J_ops = [sprepost(c, c.dag()) for c in c_ops]

    

    N = len(J_ops)
    I = np.zeros(N)
    
    if wlist is None:
        S = np.zeros((N, N,1))
        wlist=[0.]
    else:
        S = np.zeros((N, N,len(wlist)))
        
    if sparse == False: 
        rhoss_vec = mat2vec(rhoss.full()).ravel()
        for k,w in enumerate(wlist):
            R = pseudo_inverse(L, rhoss=rhoss, w= w, sparse = sparse, method=method)
            for i, Ji in enumerate(J_ops):
                for j, Jj in enumerate(J_ops):
                    if i == j:
                        I[i] = expect_rho_vec(Ji.data, rhoss_vec, 1)
                        S[i, j,k] = I[i]
                    S[i, j,k] -= expect_rho_vec((Ji * R * Jj 
                                                + Jj * R * Ji).data,
                                                rhoss_vec, 1)
    else:
        if method == "direct":
            N = np.prod(L.dims[0][0])
            
            rhoss_vec = operator_to_vector(rhoss)
            
            tr_op = tensor([identity(n) for n in L.dims[0][0]])
            tr_op_vec = operator_to_vector(tr_op)
            
            Pop = sp.kron(rhoss_vec.data, tr_op_vec.data.T, format='csr')
            Iop = sp.eye(N*N, N*N, format='csr')
            Q = Iop - Pop
#.........这里部分代码省略.........

def spectrum_ss(H, wlist, c_ops, a_op, b_op):
    """
    Calculate the spectrum corresponding to a correlation function
    :math:`\left<A(\\tau)B(0)\\right>`, i.e., the Fourier transform of the
    correlation function:

    .. math::

        S(\omega) = \int_{-\infty}^{\infty} \left<A(\\tau)B(0)\\right>
        e^{-i\omega\\tau} d\\tau.

    Parameters
    ----------

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

    wlist : *list* / *array*
        list of frequencies for :math:`\\omega`.

    c_ops : list of :class:`qutip.qobj`
        list of collapse operators.

    a_op : :class:`qutip.qobj`
        operator A.

    b_op : :class:`qutip.qobj`
        operator B.

    Returns
    -------

    spectrum: *array*
        An *array* with spectrum :math:`S(\omega)` for the frequencies
        specified in `wlist`.

    """

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

    # contruct the Liouvillian
    L = liouvillian(H, c_ops)

    # find the steady state density matrix and a_op and b_op expecation values
    rho0 = steadystate(L)

    a_op_ss = expect(a_op, rho0)
    b_op_ss = expect(b_op, rho0)

    # eseries solution for (b * rho0)(t)
    es = ode2es(L, b_op * rho0)

    # correlation
    corr_es = expect(a_op, es)

    # covarience
    cov_es = corr_es - np.real(np.conjugate(a_op_ss) * b_op_ss)

    # spectrum
    spectrum = esspec(cov_es, wlist)

    return spectrum

def spectrum_pi(H, wlist, c_ops, a_op, b_op, use_pinv=False):
    """
    Calculate the spectrum corresponding to a correlation function
    :math:`\left<A(\\tau)B(0)\\right>`, i.e., the Fourier transform of the
    correlation function:

    .. math::

        S(\omega) = \int_{-\infty}^{\infty} \left<A(\\tau)B(0)\\right>
        e^{-i\omega\\tau} d\\tau.

    Parameters
    ----------

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

    wlist : *list* / *array*
        list of frequencies for :math:`\\omega`.

    c_ops : list of :class:`qutip.qobj`
        list of collapse operators.


    a_op : :class:`qutip.qobj`
        operator A.

    b_op : :class:`qutip.qobj`
        operator B.

    Returns
    -------

    s_vec: *array*
        An *array* with spectrum :math:`S(\omega)` for the frequencies
        specified in `wlist`.

    """

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

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

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

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

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

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

    s_vec = np.zeros(len(wlist))

    for idx, w in enumerate(wlist):

        if use_pinv:
            MMR = numpy.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, transpose(rho)))))
        s_vec[idx] = -2 * np.real(s[0, 0])

    return s_vec