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 add_annotation(self, state_or_vector, text, **kwargs):
        """Add a text or LaTeX annotation to Bloch sphere,
        parametrized by a qubit state or a vector.

        Parameters
        ----------
        state_or_vector : Qobj/array/list/tuple
            Position for the annotaion.
            Qobj of a qubit or a vector of 3 elements.

        text : str/unicode
            Annotation text.
            You can use LaTeX, but remember to use raw string
            e.g. r"$\\langle x \\rangle$"
            or escape backslashes
            e.g. "$\\\\langle x \\\\rangle$".

        **kwargs :
            Options as for mplot3d.axes3d.text, including:
            fontsize, color, horizontalalignment, verticalalignment.
        """
        if isinstance(state_or_vector, Qobj):
            vec = [expect(sigmax(), state_or_vector),
                   expect(sigmay(), state_or_vector),
                   expect(sigmaz(), state_or_vector)]
        elif isinstance(state_or_vector, (list, ndarray, tuple)) \
                and len(state_or_vector) == 3:
            vec = state_or_vector
        else:
            raise Exception("Position needs to be specified by a qubit " +
                            "state or a 3D vector.")
        self.annotations.append({'position': vec,
                                 'text': text,
                                 'opts': kwargs})

def spectrum_ss(H, wlist, c_op_list, 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_op_list : 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`.

    """

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

    # find the steady state density matrix and a_op and b_op expecation values
    rho0 = steady(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_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 _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 _smepdpsolve_single_trajectory(data, L, dt, times, N_store, N_substeps, rho_t, dims, c_ops, e_ops):
    """
    Internal function. See smepdpsolve.
    """
    states_list = []

    rho_t = np.copy(rho_t)
    sigma_t = np.copy(rho_t)

    prng = RandomState()  # todo: seed it
    r_jump, r_op = prng.rand(2)

    jump_times = []
    jump_op_idx = []

    for t_idx, t in enumerate(times):

        if e_ops:
            for e_idx, e in enumerate(e_ops):
                data.expect[e_idx, t_idx] += expect_rho_vec(e, rho_t)
        else:
            states_list.append(Qobj(vec2mat(rho_t), dims=dims))

        for j in range(N_substeps):

            if sigma_t.norm() < r_jump:
                # jump occurs
                p = np.array([expect(c.dag() * c, rho_t) for c in c_ops])
                p = np.cumsum(p / np.sum(p))
                n = np.where(p >= r_op)[0][0]

                # apply jump
                rho_t = c_ops[n] * rho_t * c_ops[n].dag()
                rho_t /= expect(c_ops[n].dag() * c_ops[n], rho_t)
                sigma_t = np.copy(rho_t)

                # store info about jump
                jump_times.append(times[t_idx] + dt * j)
                jump_op_idx.append(n)

                # get new random numbers for next jump
                r_jump, r_op = prng.rand(2)

            # deterministic evolution wihtout correction for norm decay
            dsigma_t = spmv(L.data, sigma_t) * dt

            # deterministic evolution with correction for norm decay
            drho_t = spmv(L.data, rho_t) * dt

            rho_t += drho_t

            # increment density matrices
            sigma_t += dsigma_t
            rho_t += drho_t

    return states_list, jump_times, jump_op_idx

def covariance_matrix(basis, rho):
    """
    The covariance matrix given a basis of operators.

    .. note::

        Experimental.
    """

    return np.array([[expect(op1*op2+op2*op1, rho)-expect(op1, rho)*expect(op2, rho) for op1 in basis] for op2 in basis])

    def testOperatorListState(self):
        """
        expect: operator list and state
        """
        res = expect([sigmax(), sigmay(), sigmaz()], fock(2, 0))
        assert_(len(res) == 3)
        assert_(all(abs(res - [0, 0, 1]) < 1e-12))

        res = expect([sigmax(), sigmay(), sigmaz()], fock_dm(2, 1))
        assert_(len(res) == 3)
        assert_(all(abs(res - [0, 0, -1]) < 1e-12))

 def testOperatorDensityMatrix(self):
     """
     expect: operator and density matrix
     """
     N = 10
     op_N = num(N)
     op_a = destroy(N)
     for n in range(N):
         e = expect(op_N, fock_dm(N, n))
         assert_(e == n)
         assert_(type(e) == float)
         e = expect(op_a, fock_dm(N, n))
         assert_(e == 0)
         assert_(type(e) == complex)

 def testOperatorKet(self):
     """
     expect: operator and ket
     """
     N = 10
     op_N = num(N)
     op_a = destroy(N)
     for n in range(N):
         e = expect(op_N, fock(N, n))
         assert_(e == n)
         assert_(type(e) == float)
         e = expect(op_a, fock(N, n))
         assert_(e == 0)
         assert_(type(e) == complex)

def wigner_covariance_matrix(a1=None, a2=None, R=None, rho=None):
    """
    calculate the wigner covariance matrix given the quadrature correlation
    matrix (R) and a state.

    .. note::

        Experimental.

    """
    if R != None:

        if rho is None:
            return np.array([[np.real(R[i,j]+R[j,i]) for i in range(4)] for j in range(4)])
        else:          
            return np.array([[np.real(expect(R[i,j]+R[j,i], rho)) for i in range(4)] for j in range(4)])

    elif a1 != None and a2 != None:

        if rho != None:
            x1 =      (a1 + a1.dag()) / np.sqrt(2)
            p1 = 1j * (a1 - a1.dag()) / np.sqrt(2)
            x2 =      (a2 + a2.dag()) / np.sqrt(2)
            p2 = 1j * (a2 - a2.dag()) / np.sqrt(2)
            return covariance_matrix([x1, p1, x2, p2], rho)
        else:
            raise ValueError("Must give rho if using field operators (a1 and a2)")
    
    else:
        raise ValueError("Must give either field operators (a1 and a2) or a precomputed correlation matrix (R)")

def _smesolve_single_trajectory(L, dt, tlist, N_store, N_substeps, rho_t, A_ops, e_ops, data, rhs, d1, d2):
    """
    Internal function. See smesolve.
    """

    dW = np.sqrt(dt) * scipy.randn(len(A_ops), N_store, N_substeps)

    states_list = []

    for t_idx, t in enumerate(tlist):

        if e_ops:
            for e_idx, e in enumerate(e_ops):
                # XXX: need to keep hilbert space structure
                data.expect[e_idx, t_idx] += expect(e, Qobj(vec2mat(rho_t)))
        else:
            states_list.append(Qobj(rho_t))  # dito

        for j in range(N_substeps):

            drho_t = spmv(L.data.data, L.data.indices, L.data.indptr, rho_t) * dt

            for a_idx, A in enumerate(A_ops):

                drho_t += rhs(L.data, rho_t, A, dt, dW[a_idx, t_idx, j], d1, d2)

            rho_t += drho_t

    return states_list

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 _ssesolve_single_trajectory(H, dt, tlist, N_store, N_substeps, psi_t, A_ops, e_ops, data, rhs, d1, d2):
    """
    Internal function. See ssesolve.
    """

    dW = np.sqrt(dt) * scipy.randn(len(A_ops), N_store, N_substeps)

    states_list = []

    for t_idx, t in enumerate(tlist):

        if e_ops:
            for e_idx, e in enumerate(e_ops):
                data.expect[e_idx, t_idx] += expect(e, Qobj(psi_t))
        else:
            states_list.append(Qobj(psi_t))

        for j in range(N_substeps):

            dpsi_t = (-1.0j * dt) * (H.data * psi_t)

            for a_idx, A in enumerate(A_ops):

                dpsi_t += rhs(H.data, psi_t, A, dt, dW[a_idx, t_idx, j], d1, d2)

            # increment and renormalize the wave function
            psi_t += dpsi_t
            psi_t /= norm(psi_t, 2)

    return states_list

def correlation_matrix(basis, rho=None):
    """
    Given a basis set of operators :math:`\\{a\\}_n`, calculate the correlation
    matrix:

    .. math::

        C_{mn} = \\langle a_m a_n \\rangle

    Parameters
    ----------

    basis : list of :class:`qutip.qobj.Qobj`
        List of operators that defines the basis for the correlation matrix.

    rho : :class:`qutip.qobj.Qobj`
        Density matrix for which to calculate the correlation matrix. If
        `rho` is `None`, then a matrix of correlation matrix operators is
        returned instead of expectation values of those operators.

    Returns
    -------

    corr_mat: *array*
        A 2-dimensional *array* of correlation values or operators.


    """

    if rho is None:
        # return array of operators
        return np.array([[op1 * op2 for op1 in basis] for op2 in basis], dtype=object)
    else:
        # return array of expectation values
        return np.array([[expect(op1 * op2, rho) for op1 in basis] for op2 in basis], dtype=object)

def correlation_es(H, rho0, tlist, taulist, c_op_list, a_op, b_op):
    """
    Internal function for calculating correlation functions using the
    exponential series solver. See :func:`correlation` usage.
    """

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

    if rho0 is None:
        rho0 = steady(L)

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

    solES_t = ode2es(L, rho0)

    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 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 wigner_covariance_matrix(a1=None, a2=None, R=None, rho=None):
    """
    Calculate the Wigner covariance matrix
    :math:`V_{ij} = \\frac{1}{2}(R_{ij} + R_{ji})`, given
    the quadrature correlation matrix
    :math:`R_{ij} = \\langle R_{i} R_{j}\\rangle -
    \\langle R_{i}\\rangle \\langle R_{j}\\rangle`, where
    :math:`R = (q_1, p_1, q_2, p_2)^T` is the vector with quadrature operators
    for the two modes.

    Alternatively, if `R = None`, and if annihilation operators `a1` and `a2`
    for the two modes are supplied instead, the quadrature correlation matrix
    is constructed from the annihilation operators before then the covariance
    matrix is calculated.

    Parameters
    ----------

    a1 : :class:`qutip.qobj.Qobj`
        Field operator for mode 1.

    a2 : :class:`qutip.qobj.Qobj`
        Field operator for mode 2.

    R : *array*
        The quadrature correlation matrix.

    rho : :class:`qutip.qobj.Qobj`
        Density matrix for which to calculate the covariance matrix.

    Returns
    -------

    cov_mat: *array*
        A 2-dimensional *array* of covariance values.

    """
    if R is not None:

        if rho is None:
            return np.array([[0.5 * np.real(R[i, j] + R[j, i]) for i in range(4)] for j in range(4)], dtype=object)
        else:
            return np.array(
                [[0.5 * np.real(expect(R[i, j] + R[j, i], rho)) for i in range(4)] for j in range(4)], dtype=object
            )

    elif a1 is not None and a2 is not None:

        if rho is not None:
            x1 = (a1 + a1.dag()) / np.sqrt(2)
            p1 = -1j * (a1 - a1.dag()) / np.sqrt(2)
            x2 = (a2 + a2.dag()) / np.sqrt(2)
            p2 = -1j * (a2 - a2.dag()) / np.sqrt(2)
            return covariance_matrix([x1, p1, x2, p2], rho)
        else:
            raise ValueError("Must give rho if using field operators " + "(a1 and a2)")

    else:
        raise ValueError("Must give either field operators (a1 and a2) " + "or a precomputed correlation matrix (R)")

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 essolve(H, rho0, tlist, c_op_list, expt_op_list):
    """
    Evolution of a state vector or density matrix (`rho0`) for a given
    Hamiltonian (`H`) and set of collapse operators (`c_op_list`), by
    expressing the ODE as an exponential series. The output is either
    the state vector at arbitrary points in time (`tlist`), or the
    expectation values of the supplied operators (`expt_op_list`).

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

    rho0 : :class:`qutip.qobj`
        Initial state density matrix.

    tlist : list/array
        ``list`` of times for :math:`t`.

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

    expt_op_list : list of :class:`qutip.qobj`
        ``list`` of :class:`qutip.qobj` operators for which to evaluate
        expectation values.


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


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

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

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

    es = ode2es(L, rho0)

    # evaluate the expectation values
    if n_expt_op == 0:
        result_list = [Qobj()] * n_tsteps
    else:
        result_list = np.zeros([n_expt_op, n_tsteps], dtype=complex)

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

    return result_list