Python Result.expect方法代碼示例

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def _smepdpsolve_generic(sso, options, progress_bar):
    """
    For internal use. See smepdpsolve.
    """
    if debug:
        logger.debug(inspect.stack()[0][3])

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

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

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

    progress_bar.start(sso.ntraj)

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

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

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

    progress_bar.finished()

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

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

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

    return data

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def _gather(sols):
    # gather list of Result objects, sols, into one.
    sol = Result()
    # 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 (config.e_num == 0):
            if (config.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 (config.options.average_expect):
                # collect expectation values, averaged
                for i in range(config.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 (config.options.average_states or
                    config.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 (config.options.average_states or config.options.average_expect):
        if (config.e_num == 0):
            sol.states = sol.states / len(sols)
        else:
            sol.expect = list(sol.expect / len(sols))
            inds = np.where(config.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 config.options.average_expect) and config.e_num != 0:
        temp = [list(sol.expect[ii]) for ii in range(ntraj)]
        for ii in range(ntraj):
            for jj in np.where(config.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

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
    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 (config.c_num != 0):
            _init_c_ops()
        if (config.e_num != 0):
            _init_e_ops()
        # set options
        qtf90.qutraj_run.n_c_ops = config.c_num
        qtf90.qutraj_run.n_e_ops = config.e_num
        qtf90.qutraj_run.ntraj = ntraj
        qtf90.qutraj_run.unravel_type = self.unravel_type
        qtf90.qutraj_run.average_states = config.options.average_states
        qtf90.qutraj_run.average_expect = config.options.average_expect
        qtf90.qutraj_run.init_result(config.psi0_shape[0],
                                     config.options.atol,
                                     config.options.rtol, mf=self.mf,
                                     norm_steps=config.norm_steps,
                                     norm_tol=config.norm_tol)
        # set optional arguments
        qtf90.qutraj_run.order = config.options.order
        qtf90.qutraj_run.nsteps = config.options.nsteps
        qtf90.qutraj_run.first_step = config.options.first_step
        qtf90.qutraj_run.min_step = config.options.min_step
        qtf90.qutraj_run.max_step = config.options.max_step
        qtf90.qutraj_run.norm_steps = config.options.norm_steps
        qtf90.qutraj_run.norm_tol = config.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 Result instance
        sol = Result()
        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 (config.e_num == 0):
            sol.states = self.get_states(len(config.tlist), ntraj)
        else:
            sol.expect = self.get_expect(len(config.tlist), ntraj)
        if (self.calc_entropy):
            sol.entropy = self.get_entropy(len(config.tlist))

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

        # deallocate stuff
        # finalize()
        return sol

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def floquet_markov_mesolve(R, ekets, rho0, tlist, e_ops, f_modes_table=None,
                           options=None, floquet_basis=True):
    """
    Solve the dynamics for the system using the Floquet-Markov master equation.
    """

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

    if opt.tidy:
        R.tidyup()

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

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

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

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

    elif isinstance(e_ops, list):

        n_expt_op = len(e_ops)
        expt_callback = False

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

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

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

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

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

    #
    # start evolution
    #
    rho = Qobj(rho0)

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

        rho.data = vec2mat(r.y)

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

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
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.solver.Result`

        An instance of the class :class:`qutip.solver.Result`, 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 Result for storing the results
    output = Result()
    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)
    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, True)

        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

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar):
    """
    Internal function for solving ME. Solve an ODE which solver parameters
    already setup (r). Calculate the required expectation values or invoke
    callback function at each time step.
    """

    #
    # prepare output array
    #
    n_tsteps = len(tlist)
    e_sops_data = []

    output = Result()
    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)

    dt = np.diff(tlist)
    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, 0)
            else:
                output.expect[m][t_idx] = expect_rho_vec(e_sops_data[m], r.y, 1)

        if t_idx < n_tsteps - 1:
            r.integrate(r.t + dt[t_idx])

    progress_bar.finished()

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

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

    return output

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def brmesolve(H, psi0, tlist, a_ops, e_ops=[], spectra_cb=[], c_ops=[],
              args={}, options=Options()):
    """
    Solve the dynamics for a system using the Bloch-Redfield master equation.

    .. note::

        This solver does not currently support time-dependent Hamiltonians.

    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.

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

    args : *dictionary*
        Placeholder for future implementation, kept for API consistency.

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

    Returns
    -------

    result: :class:`qutip.solver.Result`

        An instance of the class :class:`qutip.solver.Result`, which contains
        either an array 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, c_ops)

    output = Result()
    output.solver = "brmesolve"
    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

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar):
    """
    Internal function for solving ME. Solve an ODE which solver parameters
    already setup (r). Calculate the required expectation values or invoke
    callback function at each time step.
    """

    #
    # prepare output array
    #
    n_tsteps = len(tlist)
    e_sops_data = []

    output = Result()
    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)

    dt = np.diff(tlist)
    for t_idx, t in enumerate(tlist):
        progress_bar.update(t_idx)

        if not r.successful():
            raise Exception("ODE integration error: Try to increase "
                            "the allowed number of substeps by increasing "
                            "the nsteps parameter in the Options class.")

        if opt.store_states or expt_callback:
            rho.data = dense2D_to_fastcsr_fmode(vec2mat(r.y), rho.shape[0], rho.shape[1])

            if opt.store_states:
                output.states.append(Qobj(rho, isherm=True))

            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, 0)
            else:
                output.expect[m][t_idx] = expect_rho_vec(e_sops_data[m],
                                                         r.y, 1)

        if t_idx < n_tsteps - 1:
            r.integrate(r.t + dt[t_idx])

    progress_bar.finished()

    if (not opt.rhs_reuse) and (config.tdname is not None):
        _cython_build_cleanup(config.tdname)

    if opt.store_final_state:
        rho.data = dense2D_to_fastcsr_fmode(vec2mat(r.y), rho.shape[0], rho.shape[1])
        output.final_state = Qobj(rho, dims=rho0.dims, isherm=True)

    return output

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]

#.........這裏部分代碼省略.........
    _ode.set_integrator('zvode', method=options.method, 
                    order=options.order, atol=options.atol, 
                    rtol=options.rtol, nsteps=options.nsteps,
                    first_step=options.first_step, 
                    min_step=options.min_step,
                    max_step=options.max_step)
    _ode.set_initial_value(initial_vector, tlist[0])
    exec(code, locals())
    
    #
    # prepare output array
    #
    n_tsteps = len(tlist)
    e_sops_data = []

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

    if options.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 = []
            options.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:
                    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
    #
    if type(progress_bar)==BaseProgressBar and verbose:
        _run_time = time.time()
    
    progress_bar.start(n_tsteps)

    rho = Qobj(rho0)

    dt = np.diff(tlist)
    for t_idx, t in enumerate(tlist):
        progress_bar.update(t_idx)

        if not _ode.successful():
            raise Exception("ODE integration error: Try to increase "
                            "the allowed number of substeps by increasing "
                            "the nsteps parameter in the Options class.")

        if options.store_states or expt_callback:
            rho.data = dense2D_to_fastcsr_fmode(vec2mat(_ode.y), rho.shape[0], rho.shape[1])

            if options.store_states:
                output.states.append(Qobj(rho, isherm=True))

            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],
                                                         _ode.y, 0)
            else:
                output.expect[m][t_idx] = expect_rho_vec(e_sops_data[m],
                                                         _ode.y, 1)

        if t_idx < n_tsteps - 1:
            _ode.integrate(_ode.t + dt[t_idx])

    progress_bar.finished()
    
    if type(progress_bar)==BaseProgressBar and verbose:
        print('BR runtime:', time.time()-_run_time)

    if (not options.rhs_reuse) and (config.tdname is not None):
        _cython_build_cleanup(config.tdname)
    
    if options.store_final_state:
        rho.data = dense2D_to_fastcsr_fmode(vec2mat(_ode.y), rho.shape[0], rho.shape[1])
        output.final_state = Qobj(rho, dims=rho0.dims, isherm=True)

    return output

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def brmesolve(H, psi0, tlist, a_ops=[], e_ops=[], c_ops=[],
              args={}, use_secular=True, sec_cutoff = 0.1,
              tol=qset.atol,
              spectra_cb=None, options=None,
              progress_bar=None, _safe_mode=True, verbose=False):
    """
    Solves for the dynamics of a system using the Bloch-Redfield master equation,
    given an input Hamiltonian, Hermitian bath-coupling terms and their associated 
    spectrum functions, as well as possible Lindblad collapse operators.
              
    For time-independent systems, the Hamiltonian must be given as a Qobj,
    whereas the bath-coupling terms (a_ops), must be written as a nested list
    of operator - spectrum function pairs, where the frequency is specified by
    the `w` variable.
              
    *Example*

        a_ops = [[a+a.dag(),lambda w: 0.2*(w>=0)]] 
              
    For time-dependent systems, the Hamiltonian, a_ops, and Lindblad collapse
    operators (c_ops), can be specified in the QuTiP string-based time-dependent
    format.  For the a_op spectra, the frequency variable must be `w`, and the 
    string cannot contain any other variables other than the possibility of having
    a time-dependence through the time variable `t`:
                            
    *Example*

        a_ops = [[a+a.dag(), '0.2*exp(-t)*(w>=0)']]
              
    It is also possible to use Cubic_Spline objects for time-dependence.  In
    the case of a_ops, Cubic_Splines must be passed as a tuple:
              
    *Example*
              
        a_ops = [ [a+a.dag(), ( f(w), g(t)] ]
              
    where f(w) and g(t) are strings or Cubic_spline objects for the bath
    spectrum and time-dependence, respectively.
              
    Finally, if one has bath-couplimg terms of the form
    H = f(t)*a + conj[f(t)]*a.dag(), then the correct input format is
              
    *Example*
    
              a_ops = [ [(a,a.dag()), (f(w), g1(t), g2(t))],... ]

    where f(w) is the spectrum of the operators while g1(t) and g2(t)
    are the time-dependence of the operators `a` and `a.dag()`, respectively 
    
    Parameters
    ----------
    H : Qobj / list
        System Hamiltonian given as a Qobj or
        nested list in string-based format.

    psi0: Qobj
        Initial density matrix or state vector (ket).

    tlist : array_like
        List of times for evaluating evolution

    a_ops : list
        Nested list of Hermitian system operators that couple to 
        the bath degrees of freedom, along with their associated
        spectra.

    e_ops : list
        List of operators for which to evaluate expectation values.

    c_ops : list
        List of system collapse operators, or nested list in
        string-based format.

    args : dict 
        Placeholder for future implementation, kept for API consistency.

    use_secular : bool {True}
        Use secular approximation when evaluating bath-coupling terms.
    
    sec_cutoff : float {0.1}
        Cutoff for secular approximation.
    
    tol : float {qutip.setttings.atol}
        Tolerance used for removing small values after 
        basis transformation.
              
    spectra_cb : list
        DEPRECIATED. Do not use.
    
    options : :class:`qutip.solver.Options`
        Options for the solver.
              
    progress_bar : BaseProgressBar
        Optional instance of BaseProgressBar, or a subclass thereof, for
        showing the progress of the simulation.

    Returns
    -------
    result: :class:`qutip.solver.Result`

#.........這裏部分代碼省略.........

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def _generic_ode_solve(r, psi0, tlist, e_ops, opt, progress_bar, dims=None):
    """
    Internal function for solving ODEs.
    """

    #
    # prepare output array
    #
    n_tsteps = len(tlist)
    output = Result()
    output.solver = "sesolve"
    output.times = tlist

    if psi0.isunitary:
        oper_evo = True
        oper_n = dims[0][0]
        norm_dim_factor = np.sqrt(oper_n)
    else:
        oper_evo = False
        norm_dim_factor = 1.0

    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:
            # fallback on storing states
            output.states = []
            opt.store_states = True
        else:
            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")

    def get_curr_state_data():
        if oper_evo:
            return vec2mat(r.y)
        else:
            return r.y

    #
    # start evolution
    #
    progress_bar.start(n_tsteps)

    dt = np.diff(tlist)
    for t_idx, t in enumerate(tlist):
        progress_bar.update(t_idx)

        if not r.successful():
            raise Exception("ODE integration error: Try to increase "
                            "the allowed number of substeps by increasing "
                            "the nsteps parameter in the Options class.")

        # get the current state / oper data if needed
        cdata = None
        if opt.store_states or opt.normalize_output or n_expt_op > 0:
            cdata = get_curr_state_data()

        if opt.normalize_output:
            # cdata *= _get_norm_factor(cdata, oper_evo)
            cdata *= norm_dim_factor / la_norm(cdata)
            if oper_evo:
                r.set_initial_value(cdata.ravel('F'), r.t)
            else:
                r.set_initial_value(cdata, r.t)

        if opt.store_states:
            output.states.append(Qobj(cdata, dims=dims))

        if expt_callback:
            # use callback method
            e_ops(t, Qobj(cdata, dims=dims))

        for m in range(n_expt_op):
            output.expect[m][t_idx] = cy_expect_psi(e_ops[m].data,
                                                    cdata, e_ops[m].isherm)

        if t_idx < n_tsteps - 1:
            r.integrate(r.t + dt[t_idx])

    progress_bar.finished()

    if not opt.rhs_reuse and config.tdname is not None:
        try:
            os.remove(config.tdname + ".pyx")
        except:
#.........這裏部分代碼省略.........

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def _generic_ode_solve(r, psi0, tlist, e_ops, opt, progress_bar,
                       state_norm_func=None, dims=None):
    """
    Internal function for solving ODEs.
    """

    #
    # prepare output array
    #
    n_tsteps = len(tlist)
    output = Result()
    output.solver = "sesolve"
    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:
            # fallback on storing states
            output.states = []
            opt.store_states = True
        else:
            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")

    #
    # start evolution
    #
    progress_bar.start(n_tsteps)

    dt = np.diff(tlist)
    for t_idx, t in enumerate(tlist):
        progress_bar.update(t_idx)

        if not r.successful():
            break

        if state_norm_func:
            data = r.y / state_norm_func(r.y)
            r.set_initial_value(data, r.t)

        if opt.store_states:
            output.states.append(Qobj(r.y, dims=dims))

        if expt_callback:
            # use callback method
            e_ops(t, Qobj(r.y, dims=psi0.dims))

        for m in range(n_expt_op):
            output.expect[m][t_idx] = cy_expect_psi(e_ops[m].data,
                                                    r.y, e_ops[m].isherm)

        if t_idx < n_tsteps - 1:
            r.integrate(r.t + dt[t_idx])

    progress_bar.finished()

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

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

    return output

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def mcsolve(H, psi0, tlist, c_ops, e_ops, ntraj=None, args={}, options=Options()):
    """Monte-Carlo evolution of a state vector :math:`|\psi \\rangle` for a
    given Hamiltonian and sets of collapse operators, and possibly, operators
    for calculating expectation values. Options for the underlying ODE solver
    are given by the Options class.

    mcsolve supports time-dependent Hamiltonians and collapse operators using
    either Python functions of strings to represent time-dependent
    coefficients. Note that, the system Hamiltonian MUST have at least one
    constant term.

    As an example of a time-dependent problem, consider a Hamiltonian with two
    terms ``H0`` and ``H1``, where ``H1`` is time-dependent with coefficient
    ``sin(w*t)``, and collapse operators ``C0`` and ``C1``, where ``C1`` is
    time-dependent with coeffcient ``exp(-a*t)``.  Here, w and a are constant
    arguments with values ``W`` and ``A``.

    Using the Python function time-dependent format requires two Python
    functions, one for each collapse coefficient. Therefore, this problem could
    be expressed as::

        def H1_coeff(t,args):
            return sin(args['w']*t)

        def C1_coeff(t,args):
            return exp(-args['a']*t)

        H=[H0,[H1,H1_coeff]]

        c_op_list=[C0,[C1,C1_coeff]]

        args={'a':A,'w':W}

    or in String (Cython) format we could write::

        H=[H0,[H1,'sin(w*t)']]

        c_op_list=[C0,[C1,'exp(-a*t)']]

        args={'a':A,'w':W}

    Constant terms are preferably placed first in the Hamiltonian and collapse
    operator lists.

    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
        single collapse operator or ``list`` or ``array`` of collapse
        operators.
    e_ops : array_like
        single operator or ``list`` or ``array`` of operators for calculating
        expectation values.
    args : dict
        Arguments for time-dependent Hamiltonian and collapse operator terms.
    options : Options
        Instance of ODE solver options.

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

    """

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

    if ntraj is None:
        ntraj = options.ntraj

    if not psi0.isket:
        raise Exception("Initial state must be a state vector.")

    if isinstance(c_ops, Qobj):
        c_ops = [c_ops]

    if isinstance(e_ops, Qobj):
        e_ops = [e_ops]

    if isinstance(e_ops, dict):
        e_ops_dict = e_ops
        e_ops = [e for e in e_ops.values()]
    else:
        e_ops_dict = None

    config.options = options
    if isinstance(ntraj, list):
        config.progress_bar = TextProgressBar(max(ntraj))
    else:
        config.progress_bar = TextProgressBar(ntraj)

    # set num_cpus to the value given in qutip.settings if none in Options
#.........這裏部分代碼省略.........

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def mcsolve_f90(H, psi0, tlist, c_ops, e_ops, ntraj=None,
                options=Options(), 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 : Options
        Instance of 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 : Result
        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.")
    config.options = options
    # set num_cpus to the value given in qutip.settings
    # if none in Options
    if not config.options.num_cpus:
        config.options.num_cpus = qutip.settings.num_cpus
    # set initial value data
    if options.tidy:
        config.psi0 = psi0.tidyup(options.atol).full()
    else:
        config.psi0 = psi0.full()
    config.psi0_dims = psi0.dims
    config.psi0_shape = psi0.shape
    # set general items
    config.tlist = tlist
    if isinstance(ntraj, (list, np.ndarray)):
        raise Exception("ntraj as list argument is not supported.")
    else:
        config.ntraj = ntraj
        # ntraj_list = [ntraj]
    # set norm finding constants
    config.norm_tol = options.norm_tol
    config.norm_steps = options.norm_steps

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

    # 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
#.........這裏部分代碼省略.........

# 需要導入模塊: from qutip.solver import Result [as 別名]
# 或者: from qutip.solver.Result import expect [as 別名]
def mcsolve(H, psi0, tlist, c_ops, e_ops, ntraj=None,
            args={}, options=None, progress_bar=True,
            map_func=None, map_kwargs=None):
    """Monte Carlo evolution of a state vector :math:`|\psi \\rangle` for a
    given Hamiltonian and sets of collapse operators, and possibly, operators
    for calculating expectation values. Options for the underlying ODE solver
    are given by the Options class.

    mcsolve supports time-dependent Hamiltonians and collapse operators using
    either Python functions of strings to represent time-dependent
    coefficients. Note that, the system Hamiltonian MUST have at least one
    constant term.

    As an example of a time-dependent problem, consider a Hamiltonian with two
    terms ``H0`` and ``H1``, where ``H1`` is time-dependent with coefficient
    ``sin(w*t)``, and collapse operators ``C0`` and ``C1``, where ``C1`` is
    time-dependent with coeffcient ``exp(-a*t)``.  Here, w and a are constant
    arguments with values ``W`` and ``A``.

    Using the Python function time-dependent format requires two Python
    functions, one for each collapse coefficient. Therefore, this problem could
    be expressed as::

        def H1_coeff(t,args):
            return sin(args['w']*t)

        def C1_coeff(t,args):
            return exp(-args['a']*t)

        H = [H0, [H1, H1_coeff]]

        c_ops = [C0, [C1, C1_coeff]]

        args={'a': A, 'w': W}

    or in String (Cython) format we could write::

        H = [H0, [H1, 'sin(w*t)']]

        c_ops = [C0, [C1, 'exp(-a*t)']]

        args={'a': A, 'w': W}

    Constant terms are preferably placed first in the Hamiltonian and collapse
    operator lists.

    Parameters
    ----------
    H : :class:`qutip.Qobj`
        System Hamiltonian.

    psi0 : :class:`qutip.Qobj`
        Initial state vector

    tlist : array_like
        Times at which results are recorded.

    ntraj : int
        Number of trajectories to run.

    c_ops : array_like
        single collapse operator or ``list`` or ``array`` of collapse
        operators.

    e_ops : array_like
        single operator or ``list`` or ``array`` of operators for calculating
        expectation values.

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

    options : Options
        Instance of ODE solver options.

    progress_bar: BaseProgressBar
        Optional instance of BaseProgressBar, or a subclass thereof, for
        showing the progress of the simulation. Set to None to disable the
        progress bar.

    map_func: function
        A map function for managing the calls to the single-trajactory solver.

    map_kwargs: dictionary
        Optional keyword arguments to the map_func function.

    Returns
    -------
    results : :class:`qutip.solver.Result`
        Object storing all results from the simulation.

    .. note::

        It is possible to reuse the random number seeds from a previous run
        of the mcsolver by passing the output Result object seeds via the
        Options class, i.e. Options(seeds=prev_result.seeds).
    """

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

#.........這裏部分代碼省略.........