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, 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 sepdpsolve_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 = "spdpsolve"
    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.jump_times = []
    data.jump_op_idx = []

    # effective hamiltonian for deterministic part
    Heff = ssdata.H
    for c in ssdata.c_ops:
        Heff += -0.5j * c.dag() * c
        
    progress_bar.start(ssdata.ntraj)

    for n in range(ssdata.ntraj):
        progress_bar.update(n)
        psi_t = ssdata.psi0.full()

        states_list, jump_times, jump_op_idx = \
            _sepdpsolve_single_trajectory(Heff, dt, ssdata.tlist,
                                          N_store, N_substeps,
                                          psi_t, ssdata.c_ops, ssdata.e_ops, 
                                          data)

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

        data.jump_times.append(jump_times)
        data.jump_op_idx.append(jump_op_idx)

    progress_bar.finished()

    # average
    data.expect = data.expect / NT

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

    return data

def ssesolve_generic(H, psi0, 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 Schrodinger equations
    A_ops = []
    for c_idx, c in enumerate(c_ops):
        A_ops.append([c.data, (c + c.dag()).data, (c.dag() * c).data])

    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

        psi_t = psi0.full()

        states_list = _ssesolve_single_trajectory(
            H, dt, tlist, N_store, N_substeps, psi_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 ssesolve_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 = "ssesolve"
    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)

    # pre-compute collapse operator combinations that are commonly needed
    # when evaluating the RHS of stochastic Schrodinger equations
    A_ops = []
    for c_idx, c in enumerate(ssdata.c_ops):
        A_ops.append([c.data,
                      (c + c.dag()).data,
                      (c - c.dag()).data,
                      (c.dag() * c).data])

    progress_bar.start(ssdata.ntraj)

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

        psi_t = ssdata.state0.full()

        states_list = _ssesolve_single_trajectory(ssdata.H, dt, ssdata.tlist, N_store,
                                                  N_substeps, psi_t, A_ops,
                                                  ssdata.e_ops, data, ssdata.rhs_func,
                                                  ssdata.d1, ssdata.d2, ssdata.d2_len,
                                                  ssdata.homogeneous, ssdata)

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

    progress_bar.finished()

    # average
    data.expect = data.expect / NT

    # standard error
    data.ss = (data.ss - NT * (data.expect ** 2)) / (NT * (NT - 1))

    return data

    def evolve_serial(self, args):

        if debug:
            print(inspect.stack()[0][3] + ":" + str(os.getpid()))

        # run ntraj trajectories for one process via fortran
        # get args
        queue, ntraj, instanceno, rngseed = args
        # initialize the problem in fortran
        _init_tlist()
        _init_psi0()
        if (self.ptrace_sel != []):
            _init_ptrace_stuff(self.ptrace_sel)
        _init_hamilt()
        if (odeconfig.c_num != 0):
            _init_c_ops()
        if (odeconfig.e_num != 0):
            _init_e_ops()
        # set options
        qtf90.qutraj_run.n_c_ops = odeconfig.c_num
        qtf90.qutraj_run.n_e_ops = odeconfig.e_num
        qtf90.qutraj_run.ntraj = ntraj
        qtf90.qutraj_run.unravel_type = self.unravel_type
        qtf90.qutraj_run.average_states = odeconfig.options.average_states 
        qtf90.qutraj_run.average_expect = odeconfig.options.average_expect
        qtf90.qutraj_run.init_odedata(odeconfig.psi0_shape[0],
                                      odeconfig.options.atol,
                                      odeconfig.options.rtol, mf=self.mf,
                                      norm_steps=odeconfig.norm_steps,
                                      norm_tol=odeconfig.norm_tol)
        # set optional arguments
        qtf90.qutraj_run.order = odeconfig.options.order
        qtf90.qutraj_run.nsteps = odeconfig.options.nsteps
        qtf90.qutraj_run.first_step = odeconfig.options.first_step
        qtf90.qutraj_run.min_step = odeconfig.options.min_step
        qtf90.qutraj_run.max_step = odeconfig.options.max_step
        qtf90.qutraj_run.norm_steps = odeconfig.options.norm_steps
        qtf90.qutraj_run.norm_tol = odeconfig.options.norm_tol
        # use sparse density matrices during computation?
        qtf90.qutraj_run.rho_return_sparse = self.sparse_dms
        # calculate entropy of reduced density matrice?
        qtf90.qutraj_run.calc_entropy = self.calc_entropy
        # run
        show_progress = 1 if debug else 0
        qtf90.qutraj_run.evolve(instanceno, rngseed, show_progress)
    

        # construct Odedata instance
        sol = Odedata()
        sol.ntraj = ntraj
        # sol.col_times = qtf90.qutraj_run.col_times
        # sol.col_which = qtf90.qutraj_run.col_which-1
        sol.col_times, sol.col_which = self.get_collapses(ntraj)
        if (odeconfig.e_num == 0):
            sol.states = self.get_states(len(odeconfig.tlist), ntraj)
        else:
            sol.expect = self.get_expect(len(odeconfig.tlist), ntraj)
        if (self.calc_entropy):
            sol.entropy = self.get_entropy(len(odeconfig.tlist))

        if (not self.serial_run):            
            # put to queue
            queue.put(sol)
            queue.join()

        # deallocate stuff
        # finalize()
        return sol

def brmesolve(H, psi0, tlist, a_ops, e_ops=[], spectra_cb=[],
              args={}, options=Odeoptions()):
    """
    Solve the dynamics for the system using the Bloch-Redfeild master equation.

    .. note::

        This solver does not currently support time-dependent Hamiltonian or
        collapse operators.

    Parameters
    ----------

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

    rho0 / psi0: :class:`qutip.qobj`
        Initial density matrix or state vector (ket).

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

    a_ops : list of :class:`qutip.qobj`
        List of system operators that couple to bath degrees of freedom.

    e_ops : list of :class:`qutip.qobj` / callback function
        List of operators for which to evaluate expectation values.

    args : *dictionary*
        Dictionary of parameters for time-dependent Hamiltonians and collapse
        operators.

    options : :class:`qutip.Qdeoptions`
        Options for the ODE solver.

    Returns
    -------

    output: :class:`qutip.odedata`

        An instance of the class :class:`qutip.odedata`, which contains either
        a list of expectation values, for operators given in e_ops, or a list
        of states for the times specified by `tlist`.
    """

    if not spectra_cb:
        # default to infinite temperature white noise
        spectra_cb = [lambda w: 1.0 for _ in a_ops]

    R, ekets = bloch_redfield_tensor(H, a_ops, spectra_cb)

    output = Odedata()
    output.times = tlist
    
    results = bloch_redfield_solve(R, ekets, psi0, tlist, e_ops, options)

    if e_ops:
        output.expect = results
    else:
        output.states = results

    return output

def _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar):
    """
    Internal function for solving ME.
    """

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

    output = Odedata()
    output.solver = "mesolve"
    output.times = tlist

    if opt.store_states:
        output.states = []

    if isinstance(e_ops, types.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:
            # fall back on storing states
            output.states = []
            opt.store_states = True
        else:
            output.expect = []
            output.num_expect = n_expt_op
            for op in e_ops:
                e_sops_data.append(spre(op).data)
                if op.isherm and rho0.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")

    #
    # start evolution
    #
    progress_bar.start(n_tsteps)

    rho = Qobj(rho0)

    for t_idx, t in enumerate(tlist):
        progress_bar.update(t_idx)

        if not r.successful():
            break

        if opt.store_states or expt_callback:
            rho.data = vec2mat(r.y)

            if opt.store_states:
                output.states.append(Qobj(rho))

            if expt_callback:
                # use callback method
                e_ops(t, rho)

        for m in range(n_expt_op):
            if output.expect[m].dtype == complex:
                output.expect[m][t_idx] = expect_rho_vec(e_sops_data[m], r.y)
            else:
                output.expect[m][t_idx] = np.real(
                    expect_rho_vec(e_sops_data[m], r.y))

        r.integrate(r.t + dt)

    progress_bar.finished()

    if not opt.rhs_reuse and odeconfig.tdname is not None:
        try:
            os.remove(odeconfig.tdname + ".pyx")
        except:
            pass

    if opt.store_final_state:
        rho.data = vec2mat(r.y)
        output.final_state = Qobj(rho)

    return output

def fsesolve(H, psi0, tlist, e_ops=[], T=None, args={}, Tsteps=100):
    """
    Solve the Schrodinger equation using the Floquet formalism.

    Parameters
    ----------

    H : :class:`qutip.qobj.Qobj`
        System Hamiltonian, time-dependent with period `T`.

    psi0 : :class:`qutip.qobj`
        Initial state vector (ket).

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

    e_ops : list of :class:`qutip.qobj` / callback function
        list of operators for which to evaluate expectation values. If this
        list is empty, the state vectors for each time in `tlist` will be
        returned instead of expectation values.

    T : float
        The period of the time-dependence of the hamiltonian.

    args : dictionary
        Dictionary with variables required to evaluate H.

    Tsteps : integer
        The number of time steps in one driving period for which to
        precalculate the Floquet modes. `Tsteps` should be an even number.

    Returns
    -------

    output : :class:`qutip.odedata.Odedata`

        An instance of the class :class:`qutip.odedata.Odedata`, which
        contains either an *array* of expectation values or an array of
        state vectors, for the times specified by `tlist`.
    """

    if not T:
        # assume that tlist span exactly one period of the driving
        T = tlist[-1]

    # find the floquet modes for the time-dependent hamiltonian
    f_modes_0, f_energies = floquet_modes(H, T, args)

    # calculate the wavefunctions using the from the floquet modes
    f_modes_table_t = floquet_modes_table(f_modes_0, f_energies,
                                          np.linspace(0, T, Tsteps + 1),
                                          H, T, args)

    # setup Odedata for storing the results
    output = Odedata()
    output.times = tlist
    output.solver = "fsesolve"

    if isinstance(e_ops, FunctionType):
        output.num_expect = 0
        expt_callback = True

    elif isinstance(e_ops, list):

        output.num_expect = len(e_ops)
        expt_callback = False

        if output.num_expect == 0:
            output.states = []
        else:
            output.expect = []
            for op in e_ops:
                if op.isherm:
                    output.expect.append(np.zeros(len(tlist)))
                else:
                    output.expect.append(np.zeros(len(tlist), dtype=complex))

    else:
        raise TypeError("e_ops must be a list Qobj or a callback function")

    psi0_fb = psi0.transform(f_modes_0, True)
    for t_idx, t in enumerate(tlist):
        f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)
        f_states_t = floquet_states(f_modes_t, f_energies, t)
        psi_t = psi0_fb.transform(f_states_t, False)

        if expt_callback:
            # use callback method
            e_ops(t, psi_t)
        else:
            # calculate all the expectation values, or output psi if
            # no expectation value operators where defined
            if output.num_expect == 0:
                output.states.append(Qobj(psi_t))
            else:
                for e_idx, e in enumerate(e_ops):
                    output.expect[e_idx][t_idx] = expect(e, psi_t)

    return output

def essolve(H, rho0, tlist, c_op_list, e_ops):
    """
    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 (`e_ops`).

    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.

    e_ops : 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(e_ops)
    n_tsteps = len(tlist)

    # Calculate the Liouvillian
    if (c_op_list is None or len(c_op_list) == 0) and isket(rho0):
        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 # XXX
    else:
        result_list = np.zeros([n_expt_op, n_tsteps], dtype=complex)

    for n, e in enumerate(e_ops):
        result_list[n, :] = expect(e, esval(es, tlist))

    data = Odedata()
    data.solver = "essolve"
    data.times = tlist
    data.expect = [np.real(result_list[n, :]) if e.isherm else result_list[n, :]
                   for n, e in enumerate(e_ops)]

    return data

def brmesolve(H, psi0, tlist, c_ops, e_ops=[], spectra_cb=[],
              args={}, options=Odeoptions()):
    """
    Solve the dynamics for the system using the Bloch-Redfeild master equation.

    .. note::

        This solver does not currently support time-dependent Hamiltonian or
        collapse operators.

    Parameters
    ----------

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

    rho0 / psi0: :class:`qutip.qobj`
        Initial density matrix or state vector (ket).

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

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

    expt_ops : list of :class:`qutip.qobj` / callback function
        List of operators for which to evaluate expectation values.

    args : *dictionary*
        Dictionary of parameters for time-dependent Hamiltonians and collapse
        operators.

    options : :class:`qutip.Qdeoptions`
        Options for the ODE solver.

    Returns
    -------

    output: :class:`qutip.odedata`

        An instance of the class :class:`qutip.odedata`, which contains either
        an *array* of expectation values for the times specified by `tlist`.
    """

    if len(spectra_cb) == 0:
        for n in range(len(c_ops)):
            # add white noise callbacks if absent
            spectra_cb.append(lambda w: 1.0)

    R, ekets = bloch_redfield_tensor(H, c_ops, spectra_cb)

    output = Odedata()
    output.times = tlist

    results = bloch_redfield_solve(R, ekets, psi0, tlist, e_ops, options)

    if len(e_ops):
        output.expect = results
    else:
        output.states = results

    return output

def _generic_ode_solve(r, psi0, tlist, expt_ops, opt,
                       state_vectorize, state_norm_func=None):
    """
    Internal function for solving ODEs.
    """

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

    output = Odedata()
    output.solver = "mesolve"
    output.times = tlist

    if isinstance(expt_ops, types.FunctionType):
        n_expt_op = 0
        expt_callback = True

    elif isinstance(expt_ops, list):

        n_expt_op = len(expt_ops)
        expt_callback = False

        if n_expt_op == 0:
            output.states = []
        else:
            output.expect = []
            output.num_expect = n_expt_op
            for op in expt_ops:
                if op.isherm and psi0.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")

    #
    # start evolution
    #
    psi = Qobj(psi0)

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

        if state_norm_func:
            psi.data = state_vectorize(r.y)
            state_norm = state_norm_func(psi.data)
            psi.data = psi.data / state_norm
            r.set_initial_value(r.y / state_norm, r.t)
        else:
            psi.data = state_vectorize(r.y)

        if expt_callback:
            # use callback method
            expt_ops(t, Qobj(psi))
        else:
            # calculate all the expectation values,
            # or output rho if no operators
            if n_expt_op == 0:
                output.states.append(Qobj(psi))  # copy psi/rho
            else:
                for m in range(0, n_expt_op):
                    output.expect[m][t_idx] = expect(expt_ops[m], psi)

        r.integrate(r.t + dt)
        t_idx += 1

    if not opt.rhs_reuse and odeconfig.tdname is not None:
        try:
            os.remove(odeconfig.tdname + ".pyx")
        except:
            pass

    return output

def _gather(sols):
    # gather list of Odedata objects, sols, into one.
    sol = Odedata()
    # sol = sols[0]
    ntraj = sum([a.ntraj for a in sols])
    sol.col_times = np.zeros((ntraj), dtype=np.ndarray)
    sol.col_which = np.zeros((ntraj), dtype=np.ndarray)
    sol.col_times[0:sols[0].ntraj] = sols[0].col_times
    sol.col_which[0:sols[0].ntraj] = sols[0].col_which
    sol.states = np.array(sols[0].states)
    sol.expect = np.array(sols[0].expect)
    if (hasattr(sols[0], 'entropy')):
        sol.entropy = np.array(sols[0].entropy)
    sofar = 0
    for j in range(1, len(sols)):
        sofar = sofar + sols[j - 1].ntraj
        sol.col_times[sofar:sofar + sols[j].ntraj] = (
            sols[j].col_times)
        sol.col_which[sofar:sofar + sols[j].ntraj] = (
            sols[j].col_which)
        if (odeconfig.e_num == 0):
            if (odeconfig.options.average_states):
                # collect states, averaged over trajectories
                sol.states += np.array(sols[j].states)
            else:
                # collect states, all trajectories
                sol.states = np.vstack((sol.states,
                                        np.array(sols[j].states)))
        else:
            if (odeconfig.options.average_expect):
                # collect expectation values, averaged
                for i in range(odeconfig.e_num):
                    sol.expect[i] += np.array(sols[j].expect[i])
            else:
                # collect expectation values, all trajectories
                sol.expect = np.vstack((sol.expect,
                                        np.array(sols[j].expect)))
        if (hasattr(sols[j], 'entropy')):
            if (odeconfig.options.average_states or odeconfig.options.average_expect):
                # collect entropy values, averaged
                sol.entropy += np.array(sols[j].entropy)
            else:
                # collect entropy values, all trajectories
                sol.entropy = np.vstack((sol.entropy,
                                         np.array(sols[j].entropy)))
    if (odeconfig.options.average_states or odeconfig.options.average_expect):
        if (odeconfig.e_num == 0):
            sol.states = sol.states / len(sols)
        else:
            sol.expect = list(sol.expect / len(sols))
            inds=np.where(odeconfig.e_ops_isherm)[0]
            for jj in inds:
                sol.expect[jj]=np.real(sol.expect[jj])
        if (hasattr(sols[0], 'entropy')):
            sol.entropy = sol.entropy / len(sols)
    
    #convert sol.expect array to list and fix dtypes of arrays
    if (not odeconfig.options.average_expect) and odeconfig.e_num!=0:
        temp=[list(sol.expect[ii]) for ii in range(ntraj)]
        for ii in range(ntraj):
            for jj in np.where(odeconfig.e_ops_isherm)[0]:
                temp[ii][jj]=np.real(temp[ii][jj])
        sol.expect=temp
    # convert to list/array to be consistent with qutip mcsolve
    sol.states = list(sol.states)
    return sol

def ssesolve_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 = "ssesolve"
    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 collapse operator combinations that are commonly needed
    # when evaluating the RHS of stochastic Schrodinger equations
    A_ops = []
    for c_idx, c in enumerate(ssdata.sc_ops):
        A_ops.append([c.data,
                      (c + c.dag()).data,
                      (c - c.dag()).data,
                      (c.dag() * c).data])

    progress_bar.start(ssdata.ntraj)

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

        psi_t = ssdata.state0.full().ravel()

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

        states_list, dW, m = _ssesolve_single_trajectory(
            ssdata.H, dt, ssdata.tlist, N_store, N_substeps, psi_t, A_ops,
            ssdata.e_ops, data, ssdata.rhs_func, ssdata.d1, ssdata.d2,
            ssdata.d2_len, ssdata.homogeneous, ssdata.distribution, ssdata.args,
            store_measurement=ssdata.store_measurement, 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 mcsolve_f90(H, psi0, tlist, c_ops, e_ops, ntraj=None,
                options=Odeoptions(), sparse_dms=True, serial=False,
                ptrace_sel=[], calc_entropy=False):
    """
    Monte-Carlo wave function solver with fortran 90 backend.
    Usage is identical to qutip.mcsolve, for problems without explicit
    time-dependence, and with some optional input:

    Parameters
    ----------
    H : qobj
        System Hamiltonian.
    psi0 : qobj
        Initial state vector
    tlist : array_like
        Times at which results are recorded.
    ntraj : int
        Number of trajectories to run.
    c_ops : array_like
        ``list`` or ``array`` of collapse operators.
    e_ops : array_like
        ``list`` or ``array`` of operators for calculating expectation values.
    options : Odeoptions
        Instance of ODE solver options.
    sparse_dms : boolean
        If averaged density matrices are returned, they will be stored as
        sparse (Compressed Row Format) matrices during computation if
        sparse_dms = True (default), and dense matrices otherwise. Dense
        matrices might be preferable for smaller systems.
    serial : boolean
        If True (default is False) the solver will not make use of the
        multiprocessing module, and simply run in serial.
    ptrace_sel: list
        This optional argument specifies a list of components to keep when
        returning a partially traced density matrix. This can be convenient for
        large systems where memory becomes a problem, but you are only
        interested in parts of the density matrix.
    calc_entropy : boolean
        If ptrace_sel is specified, calc_entropy=True will have the solver
        return the averaged entropy over trajectories in results.entropy. This
        can be interpreted as a measure of entanglement. See Phys. Rev. Lett.
        93, 120408 (2004), Phys. Rev. A 86, 022310 (2012).

    Returns
    -------
    results : Odedata
        Object storing all results from simulation.

    """
    if ntraj is None:
        ntraj = options.ntraj

    if psi0.type != 'ket':
        raise Exception("Initial state must be a state vector.")
    odeconfig.options = options
    # set num_cpus to the value given in qutip.settings
    # if none in Odeoptions
    if not odeconfig.options.num_cpus:
        odeconfig.options.num_cpus = qutip.settings.num_cpus
    # set initial value data
    if options.tidy:
        odeconfig.psi0 = psi0.tidyup(options.atol).full()
    else:
        odeconfig.psi0 = psi0.full()
    odeconfig.psi0_dims = psi0.dims
    odeconfig.psi0_shape = psi0.shape
    # set general items
    odeconfig.tlist = tlist
    if isinstance(ntraj, (list, np.ndarray)):
        raise Exception("ntraj as list argument is not supported.")
    else:
        odeconfig.ntraj = ntraj
        # ntraj_list = [ntraj]
    # set norm finding constants
    odeconfig.norm_tol = options.norm_tol
    odeconfig.norm_steps = options.norm_steps

    if not options.rhs_reuse:
        odeconfig.soft_reset()
        # no time dependence
        odeconfig.tflag = 0
        # check for collapse operators
        if len(c_ops) > 0:
            odeconfig.cflag = 1
        else:
            odeconfig.cflag = 0
        # Configure data
        _mc_data_config(H, psi0, [], c_ops, [], [], e_ops, options, odeconfig)

    # Load Monte Carlo class
    mc = _MC_class()
    # Set solver type
    if (options.method == 'adams'):
        mc.mf = 10
    elif (options.method == 'bdf'):
        mc.mf = 22
    else:
        if debug:
            print('Unrecognized method for ode solver, using "adams".')
        mc.mf = 10
#.........这里部分代码省略.........

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 = Odeoptions()
    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 = Odedata()
    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, True)

    #
    # 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.o