From 26ac589e613d8b69792bbc71e1d4cb12d2c39d6c Mon Sep 17 00:00:00 2001 From: Eric Giguere Date: Thu, 8 Dec 2022 15:15:15 -0500 Subject: [PATCH 1/8] Constant capitalization solver --- doc/apidoc/classes.rst | 4 +- doc/guide/dynamics/dynamics-master.rst | 10 +-- doc/guide/dynamics/dynamics-monte.rst | 88 ++++++++++++++++++------- doc/guide/dynamics/dynamics-options.rst | 8 +-- qutip/solver/correlation.py | 14 ++-- qutip/solver/mcsolve.py | 10 +-- qutip/solver/mesolve.py | 10 +-- qutip/solver/propagator.py | 13 ++-- qutip/solver/result.py | 20 ++++-- qutip/solver/sesolve.py | 8 +-- qutip/tests/solver/test_correlation.py | 4 +- qutip/tests/solver/test_integrator.py | 25 +++---- qutip/tests/solver/test_mcsolve.py | 8 +-- qutip/tests/solver/test_mesolve.py | 16 ++--- qutip/tests/solver/test_options.py | 4 +- qutip/tests/solver/test_propagator.py | 18 ++--- qutip/tests/solver/test_sesolve.py | 18 ++--- 17 files changed, 167 insertions(+), 111 deletions(-) diff --git a/doc/apidoc/classes.rst b/doc/apidoc/classes.rst index 1e85f82f1e..1fd31c4764 100644 --- a/doc/apidoc/classes.rst +++ b/doc/apidoc/classes.rst @@ -42,10 +42,10 @@ Distributions Solver ------ -.. autoclass:: qutip.solver.sesolve.SeSolver +.. autoclass:: qutip.solver.sesolve.SESolver :members: -.. autoclass:: qutip.solver.mesolve.MeSolver +.. autoclass:: qutip.solver.mesolve.MESolver :members: .. autoclass:: qutip.solver.brmesolve.BRSolver diff --git a/doc/guide/dynamics/dynamics-master.rst b/doc/guide/dynamics/dynamics-master.rst index 49720ea763..5384ee4d6f 100644 --- a/doc/guide/dynamics/dynamics-master.rst +++ b/doc/guide/dynamics/dynamics-master.rst @@ -33,17 +33,17 @@ For example, the time evolution of a quantum spin-1/2 system with tunneling rate >>> H = 2*np.pi * 0.1 * sigmax() >>> psi0 = basis(2, 0) >>> times = np.linspace(0.0, 10.0, 20) - >>> result = sesolve(H, psi0, times, [sigmaz()]) + >>> result = sesolve(H, psi0, times, e_ops=[sigmaz()]) The brackets in the fourth argument is an empty list of collapse operators, since we consider unitary evolution in this example. See the next section for examples on how dissipation is included by defining a list of collapse operators. -The function returns an instance of :class:`qutip.solve.solver.Result`, as described in the previous section :ref:`solver_result`. The attribute ``expect`` in ``result`` is a list of expectation values for the operators that are included in the list in the fifth argument. Adding operators to this list results in a larger output list returned by the function (one array of numbers, corresponding to the times in times, for each operator) +The function returns an instance of :class:`qutip.Result`, as described in the previous section :ref:`solver_result`. The attribute ``expect`` in ``result`` is a list of expectation values for the operators that are included in the list in the fifth argument. Adding operators to this list results in a larger output list returned by the function (one array of numbers, corresponding to the times in times, for each operator) .. plot:: :context: - >>> result = sesolve(H, psi0, times, [sigmaz(), sigmay()]) + >>> result = sesolve(H, psi0, times, e_ops=[sigmaz(), sigmay()]) >>> result.expect # doctest: +NORMALIZE_WHITESPACE [array([ 1. , 0.78914057, 0.24548559, -0.40169513, -0.8794735 , -0.98636142, -0.67728219, -0.08258023, 0.54694721, 0.94581685, @@ -169,7 +169,7 @@ the previously empty list in the fourth parameter to the :func:`qutip.mesolve` f :context: >>> times = np.linspace(0.0, 10.0, 100) - >>> result = mesolve(H, psi0, times, [np.sqrt(0.05) * sigmax()], [sigmaz(), sigmay()]) + >>> result = mesolve(H, psi0, times, [np.sqrt(0.05) * sigmax()], e_ops=[sigmaz(), sigmay()]) >>> fig, ax = plt.subplots() >>> ax.plot(times, result.expect[0]) # doctest: +SKIP >>> ax.plot(times, result.expect[1]) # doctest: +SKIP @@ -192,7 +192,7 @@ Now a slightly more complex example: Consider a two-level atom coupled to a leak >>> a = tensor(qeye(2), destroy(10)) >>> sm = tensor(destroy(2), qeye(10)) >>> H = 2 * np.pi * a.dag() * a + 2 * np.pi * sm.dag() * sm + 2 * np.pi * 0.25 * (sm * a.dag() + sm.dag() * a) - >>> result = mesolve(H, psi0, times, [np.sqrt(0.1)*a], [a.dag()*a, sm.dag()*sm]) + >>> result = mesolve(H, psi0, times, [np.sqrt(0.1)*a], e_ops=[a.dag()*a, sm.dag()*sm]) >>> plt.figure() # doctest: +SKIP >>> plt.plot(times, result.expect[0]) # doctest: +SKIP >>> plt.plot(times, result.expect[1]) # doctest: +SKIP diff --git a/doc/guide/dynamics/dynamics-monte.rst b/doc/guide/dynamics/dynamics-monte.rst index 797892b824..c1f528e3e5 100644 --- a/doc/guide/dynamics/dynamics-monte.rst +++ b/doc/guide/dynamics/dynamics-monte.rst @@ -4,6 +4,11 @@ Monte Carlo Solver ******************************************* +.. plot:: + :include-source: False + + from qutip.solver.mcsolve import * + .. _monte-intro: Introduction @@ -75,7 +80,7 @@ To illustrate the use of the Monte Carlo evolution of quantum systems in QuTiP, a = tensor(qeye(2), destroy(10)) sm = tensor(destroy(2), qeye(10)) H = 2*np.pi*a.dag()*a + 2*np.pi*sm.dag()*sm + 2*np.pi*0.25*(sm*a.dag() + sm.dag()*a) - data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], [a.dag() * a, sm.dag() * sm]) + data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], e_ops=[a.dag() * a, sm.dag() * sm]) plt.figure() plt.plot(times, data.expect[0], times, data.expect[1]) @@ -89,7 +94,7 @@ To illustrate the use of the Monte Carlo evolution of quantum systems in QuTiP, The advantage of the Monte Carlo method over the master equation approach is that only the state vector is required to be kept in the computers memory, as opposed to the entire density matrix. For large quantum system this becomes a significant advantage, and the Monte Carlo solver is therefore generally recommended for such systems. For example, simulating a Heisenberg spin-chain consisting of 10 spins with random parameters and initial states takes almost 7 times longer using the master equation rather than Monte Carlo approach with the default number of trajectories running on a quad-CPU machine. Furthermore, it takes about 7 times the memory as well. However, for small systems, the added overhead of averaging a large number of stochastic trajectories to obtain the open system dynamics, as well as starting the multiprocessing functionality, outweighs the benefit of the minor (in this case) memory saving. Master equation methods are therefore generally more efficient when Hilbert space sizes are on the order of a couple of hundred states or smaller. -Like the master equation solver :func:`qutip.mesolve`, the Monte Carlo solver returns a :class:`qutip.solve.solver.Result` object consisting of expectation values, if the user has defined expectation value operators in the 5th argument to ``mcsolve``, or state vectors if no expectation value operators are given. If state vectors are returned, then the :class:`qutip.solve.solver.Result` returned by :func:`qutip.mcsolve` will be an array of length ``ntraj``, with each element containing an array of ket-type qobjs with the same number of elements as ``times``. Furthermore, the output :class:`qutip.solve.solver.Result` object will also contain a list of times at which collapse occurred, and which collapse operators did the collapse, in the ``col_times`` and ``col_which`` properties, respectively. +Like the master equation solver :func:`qutip.mesolve`, the Monte Carlo solver returns a :class:`qutip.Result` object consisting of expectation values, if the user has defined expectation value operators in the 5th argument to ``mcsolve``, or state vectors if no expectation value operators are given. If state vectors are returned, then the :class:`qutip.Result` returned by :func:`qutip.mcsolve` will be an array of length ``ntraj``, with each element containing an array of ket-type qobjs with the same number of elements as ``times``. Furthermore, the output :class:`qutip.Result` object will also contain a list of times at which collapse occurred, and which collapse operators did the collapse, in the ``col_times`` and ``col_which`` properties, respectively. .. _monte-ntraj: @@ -97,36 +102,46 @@ Like the master equation solver :func:`qutip.mesolve`, the Monte Carlo solver re Changing the Number of Trajectories ----------------------------------- -As mentioned earlier, by default, the ``mcsolve`` function runs 500 trajectories. This value was chosen because it gives good accuracy, Monte Carlo errors scale as :math:`1/n` where :math:`n` is the number of trajectories, and simultaneously does not take an excessive amount of time to run. However, like many other options in QuTiP you are free to change the number of trajectories to fit your needs. If we want to run 1000 trajectories in the above example, we can simply modify the call to ``mcsolve`` like: +As mentioned earlier, by default, the ``mcsolve`` function runs 500 trajectories. +This value was chosen because it gives good accuracy, Monte Carlo errors scale as :math:`1/n` where :math:`n` is the number of trajectories, and simultaneously does not take an excessive amount of time to run. +However, like many other options in QuTiP you are free to change the number of trajectories to fit your needs. +If we want to run 1000 trajectories in the above example, we can simply modify the call to ``mcsolve`` like: .. plot:: :context: close-figs data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], [a.dag() * a, sm.dag() * sm], ntraj=1000) -where we have added the keyword argument ``ntraj=1000`` at the end of the inputs. Now, the Monte Carlo solver will calculate expectation values for both operators, ``a.dag() * a, sm.dag() * sm`` averaging over 1000 trajectories. Sometimes one is also interested in seeing how the Monte Carlo trajectories converge to the master equation solution by calculating expectation values over a range of trajectory numbers. If, for example, we want to average over 1, 10, 100, and 1000 trajectories, then we can input this into the solver using: - -.. plot:: - :context: - - ntraj = [1, 10, 100, 1000] - -Keep in mind that the input list must be in ascending order since the total number of trajectories run by ``mcsolve`` will be calculated using the last element of ``ntraj``. In this case, we need to use an extra index when getting the expectation values from the :class:`qutip.solve.solver.Result` object returned by ``mcsolve``. In the above example using: +where we have added the keyword argument ``ntraj=1000`` at the end of the inputs. +Now, the Monte Carlo solver will calculate expectation values for both operators, ``a.dag() * a, sm.dag() * sm`` averaging over 1000 trajectories. +Sometimes one is also interested in seeing how the Monte Carlo trajectories converge to the master equation solution by calculating expectation values over a range of trajectory numbers. +If, for example, we want to average over 1, 10, 100, and 1000 trajectories, we must first run the solver keeping expectation values for each trajectories: .. plot:: :context: - data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], [a.dag() * a, sm.dag() * sm], ntraj=[1, 10, 100, 1000]) + data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], e_ops=[a.dag() * a, sm.dag() * sm], options={"keep_runs_results": True}, ntraj=1000) we can extract the relevant expectation values using: .. plot:: :context: - expt1 = data.expect[0] - expt10 = data.expect[1] - expt100 = data.expect[2] - expt1000 = data.expect[3] + expt1 = data.expect_traj_avg(1) + expt10 = data.expect_traj_avg(10) + expt100 = data.expect_traj_avg(100) + expt1000 = data.expect_traj_avg(1000) + + plt.figure() + plt.plot(times, expt1, label="ntraj=1") + plt.plot(times, expt10, label="ntraj=10") + plt.plot(times, expt100, label="ntraj=100") + plt.plot(times, expt1000, label="ntraj=1000") + plt.title('Monte Carlo time evolution') + plt.xlabel('Time') + plt.ylabel('Expectation values') + plt.legend() + plt.show() .. _monte-reuse: @@ -136,12 +151,16 @@ Reusing Hamiltonian Data .. note:: This section covers a specialized topic and may be skipped if you are new to QuTiP. -In order to solve a given simulation as fast as possible, the solvers in QuTiP take the given input operators and break them down into simpler components before passing them on to the ODE solvers. Although these operations are reasonably fast, the time spent organizing data can become appreciable when repeatedly solving a system over, for example, many different initial conditions. In cases such as this, the Hamiltonian and other operators may be reused after the initial configuration, thus speeding up calculations. Note that, unless you are planning to reuse the data many times, this functionality will not be very useful. + +In order to solve a given simulation as fast as possible, the solvers in QuTiP take the given input operators and break them down into simpler components before passing them on to the ODE solvers. +Although these operations are reasonably fast, the time spent organizing data can become appreciable when repeatedly solving a system over, for example, many different initial conditions. +In cases such as this, the Monte Carlo Solver may be reused after the initial configuration, thus speeding up calculations. + Using the previous example, we will calculate the dynamics for two different initial states, with the Hamiltonian data being reused on the second call .. plot:: - :context: + :context: close-figs times = np.linspace(0.0, 10.0, 200) psi0 = tensor(fock(2, 0), fock(10, 5)) @@ -149,10 +168,10 @@ Using the previous example, we will calculate the dynamics for two different ini sm = tensor(destroy(2), qeye(10)) H = 2*np.pi*a.dag()*a + 2*np.pi*sm.dag()*sm + 2*np.pi*0.25*(sm*a.dag() + sm.dag()*a) - data1 = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], [a.dag() * a, sm.dag() * sm]) + solver = MCSolver(H, c_ops=[np.sqrt(0.1) * a]) + data1 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=1000) psi1 = tensor(fock(2, 0), coherent(10, 2 - 1j)) - opts = SolverOptions() # Run a second time, reusing RHS - data2 = mcsolve(H, psi1, times, [np.sqrt(0.1) * a], [a.dag() * a, sm.dag() * sm], options=opts) + data2 = solver.run(psi1, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=1000) plt.figure() plt.plot(times, data1.expect[0], times, data1.expect[1], lw=2) @@ -165,4 +184,29 @@ Using the previous example, we will calculate the dynamics for two different ini .. guide-dynamics-mc2: -In addition to the initial state, one may reuse the Hamiltonian data when changing the number of trajectories ``ntraj`` or simulation times ``times``. The reusing of Hamiltonian data is also supported for time-dependent Hamiltonians. See :ref:`time` for further details. +In addition to the initial state, one may reuse the Hamiltonian data when adding more trajectories to the result or changing simulation times ``times``. + + +.. plot:: + :context: close-figs + + times = np.linspace(0.0, 10.0, 200) + psi0 = tensor(fock(2, 0), fock(10, 5)) + a = tensor(qeye(2), destroy(10)) + sm = tensor(destroy(2), qeye(10)) + + H = 2*np.pi*a.dag()*a + 2*np.pi*sm.dag()*sm + 2*np.pi*0.25*(sm*a.dag() + sm.dag()*a) + solver = MCSolver(H, c_ops=[np.sqrt(0.1) * a]) + data1 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=10) + data2 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=10) + data_merged = data1 + data2 + + plt.figure() + plt.plot(times, data1.expect[0], times, data1.expect[1], lw=2) + plt.plot(times, data2.expect[0], '--', times, data2.expect[1], '--', lw=2) + plt.plot(times, data_merged.expect[0], ':', times, data_merged.expect[1], ':', lw=2) + plt.title('Monte Carlo time evolution') + plt.xlabel('Time', fontsize=14) + plt.ylabel('Expectation values', fontsize=14) + plt.legend(("cavity photon number", "atom excitation probability")) + plt.show() diff --git a/doc/guide/dynamics/dynamics-options.rst b/doc/guide/dynamics/dynamics-options.rst index 4a7e29979b..587000ca84 100644 --- a/doc/guide/dynamics/dynamics-options.rst +++ b/doc/guide/dynamics/dynamics-options.rst @@ -6,7 +6,7 @@ Setting Options for the Dynamics Solvers .. testsetup:: [dynamics_options] - from qutip.solver.mesolve import MeSolver, mesolve + from qutip.solver.mesolve import MESolver, mesolve import numpy as np Occasionally it is necessary to change the built in parameters of the dynamics solvers used by for example the :func:`qutip.mesolve` and :func:`qutip.mcsolve` functions. @@ -21,13 +21,13 @@ Supported solver options and their default can be seen using the class interface .. testcode:: [dynamics_options] - help(MeSolver.options) + help(MESolver.options) Options supported by the ODE integration depend on the "method" options of the solver, they can be listed through the integrator method of the solvers: .. testcode:: [dynamics_options] - help(MeSolver.integrator("adams").options) + help(MESolver.integrator("adams").options) See `Integrator <../../apidoc/classes.html#classes-ode>`_ for a list of supported methods. @@ -44,4 +44,4 @@ To use these new settings we can use the keyword argument ``options`` in either or:: - >>> McSolver(H0, c_op_list, options=options) + >>> MCSolver(H0, c_op_list, options=options) diff --git a/qutip/solver/correlation.py b/qutip/solver/correlation.py index 239aa80cc4..a2d956065e 100644 --- a/qutip/solver/correlation.py +++ b/qutip/solver/correlation.py @@ -11,8 +11,8 @@ qeye, Qobj, QobjEvo, liouvillian, spre, unstack_columns, stack_columns, tensor, qzero, expect ) -from .mesolve import MeSolver -from .mcsolve import McSolver +from .mesolve import MESolver +from .mcsolve import MCSolver from .brmesolve import BRSolver from .heom.bofin_solvers import HEOMSolver @@ -419,10 +419,10 @@ def _make_solver(H, c_ops, args, options, solver): H = QobjEvo(H, args=args) c_ops = [QobjEvo(c_op, args=args) for c_op in c_ops] if solver == "me": - solver_instance = MeSolver(H, c_ops, options=options) + solver_instance = MESolver(H, c_ops, options=options) elif solver == "es": options = {"method": "diag"} - solver_instance = MeSolver(H, c_ops, options=options) + solver_instance = MESolver(H, c_ops, options=options) elif solver == "mc": raise ValueError("MC solver for correlation has been removed") return solver_instance @@ -441,7 +441,7 @@ def correlation_3op(solver, state0, tlist, taulist, A=None, B=None, C=None): Parameters ---------- - solver : :class:`MeSolver`, :class:`BRSolver` + solver : :class:`MESolver`, :class:`BRSolver` Qutip solver for an open system. state0 : :class:`Qobj` Initial state density matrix :math:`\rho(t_0)` or state vector @@ -472,9 +472,9 @@ def correlation_3op(solver, state0, tlist, taulist, A=None, B=None, C=None): B = QobjEvo(qeye(dims) if B in [None, 1] else B) C = QobjEvo(qeye(dims) if C in [None, 1] else C) - if isinstance(solver, (MeSolver, BRSolver)): + if isinstance(solver, (MESolver, BRSolver)): out = _correlation_3op_dm(solver, state0, tlist, taulist, A, B, C) - elif isinstance(solver, McSolver): + elif isinstance(solver, MCSolver): raise TypeError("Monte Carlo support for correlation was removed. " "Please, tell us on GitHub issues if you need it!") elif isinstance(solver, HEOMSolver): diff --git a/qutip/solver/mcsolve.py b/qutip/solver/mcsolve.py index fb38f93cb9..c9669e2ee1 100644 --- a/qutip/solver/mcsolve.py +++ b/qutip/solver/mcsolve.py @@ -1,4 +1,4 @@ -__all__ = ['mcsolve', "McSolver"] +__all__ = ['mcsolve', "MCSolver"] import warnings @@ -8,7 +8,7 @@ from .multitraj import MultiTrajSolver from .solver_base import Solver from .result import McResult, Result -from .mesolve import mesolve, MeSolver +from .mesolve import mesolve, MESolver import qutip.core.data as _data from time import time @@ -138,7 +138,7 @@ def mcsolve(H, state, tlist, c_ops=(), e_ops=None, ntraj=500, *, options = { key: options[key] for key in options - if key in MeSolver.solver_options + if key in MESolver.solver_options } return mesolve(H, state, tlist, e_ops=e_ops, args=args, options=options) @@ -149,7 +149,7 @@ def mcsolve(H, state, tlist, c_ops=(), e_ops=None, ntraj=500, *, "with the options `keep_runs_results=True`." ) - mc = McSolver(H, c_ops, options=options) + mc = MCSolver(H, c_ops, options=options) result = mc.run(state, tlist=tlist, ntraj=ntraj, e_ops=e_ops, seed=seeds, target_tol=target_tol, timeout=timeout) return result @@ -317,7 +317,7 @@ def integrator_options(self): # ----------------------------------------------------------------------------- # MONTE CARLO CLASS # ----------------------------------------------------------------------------- -class McSolver(MultiTrajSolver): +class MCSolver(MultiTrajSolver): r""" Monte Carlo Solver of a state vector :math:`|\psi \rangle` for a given Hamiltonian and sets of collapse operators. Options for the diff --git a/qutip/solver/mesolve.py b/qutip/solver/mesolve.py index ed1cd23e42..ed865dcabd 100644 --- a/qutip/solver/mesolve.py +++ b/qutip/solver/mesolve.py @@ -3,7 +3,7 @@ equation. """ -__all__ = ['mesolve', 'MeSolver'] +__all__ = ['mesolve', 'MESolver'] import numpy as np from time import time @@ -11,7 +11,7 @@ from ..core import stack_columns, unstack_columns from ..core.data import to from .solver_base import Solver -from .sesolve import sesolve, SeSolver +from .sesolve import sesolve, SESolver def mesolve(H, rho0, tlist, c_ops=None, e_ops=None, args=None, options=None): @@ -138,12 +138,12 @@ def mesolve(H, rho0, tlist, c_ops=None, e_ops=None, args=None, options=None): return sesolve(H, rho0, tlist, e_ops=e_ops, args=args, options=options) - solver = MeSolver(H, c_ops, options=options) + solver = MESolver(H, c_ops, options=options) return solver.run(rho0, tlist, e_ops=e_ops) -class MeSolver(SeSolver): +class MESolver(SESolver): """ Master equation evolution of a density matrix for a given Hamiltonian and set of collapse operators, or a Liouvillian. @@ -170,7 +170,7 @@ class MeSolver(SeSolver): of Liouvillian superoperators. None is equivalent to an empty list. options : dict, optional - Options for the solver, see :obj:`SeSolver.options` and + Options for the solver, see :obj:`SESolver.options` and `Integrator <./classes.html#classes-ode>`_ for a list of all options. attributes diff --git a/qutip/solver/propagator.py b/qutip/solver/propagator.py index 53253abe83..fbc83698ad 100644 --- a/qutip/solver/propagator.py +++ b/qutip/solver/propagator.py @@ -5,9 +5,8 @@ from .. import Qobj, qeye, unstack_columns, QobjEvo from ..core import data as _data -from .mesolve import mesolve, MeSolver -from .sesolve import sesolve, SeSolver -from .mcsolve import McSolver +from .mesolve import mesolve, MESolver +from .sesolve import sesolve, SESolver from .heom.bofin_solvers import HEOMSolver from .solver_base import Solver from .multitraj import MultiTrajSolver @@ -116,12 +115,12 @@ class Propagator: ---------- system : :class:`Qobj`, :class:`QobjEvo`, :class:`Solver` Possibly time-dependent system driving the evolution, either already - packaged in a solver, such as :class:`SeSolver` or :class:`BrSolver`, + packaged in a solver, such as :class:`SESolver` or :class:`BRSolver`, or the Liouvillian or Hamiltonian as a :class:`Qobj`, :class:`QobjEvo`. ``list`` of [:class:`Qobj`, :class:`Coefficient`] or callable that can be made into :class:`QobjEvo` are also accepted. - Solvers that run non-deterministacilly, such as :class:`McSolver`, are + Solvers that run non-deterministacilly, such as :class:`MCSolver`, are not supported. c_ops : list, optional @@ -163,9 +162,9 @@ def __init__(self, system, *, c_ops=(), args=None, options=None, Hevo = QobjEvo(system, args=args) c_ops = [QobjEvo(op, args=args) for op in c_ops] if Hevo.issuper or c_ops: - self.solver = MeSolver(Hevo, c_ops=c_ops, options=options) + self.solver = MESolver(Hevo, c_ops=c_ops, options=options) else: - self.solver = SeSolver(Hevo, options=options) + self.solver = SESolver(Hevo, options=options) self.times = [0] self.invs = [None] diff --git a/qutip/solver/result.py b/qutip/solver/result.py index 8c8ffd3719..24c5a50562 100644 --- a/qutip/solver/result.py +++ b/qutip/solver/result.py @@ -785,6 +785,11 @@ def e_data_traj_avg(self, ntraj=-1): """ if not self.trajectories: return None + if ntraj > len(self.trajectories): + raise ValueError( + f"Cannot compute statistic for {ntraj} trajectories. " + f"Only {len(self.trajectories)} trajectories are stored." + ) return { k: np.mean(np.stack([ traj.e_data[k] for traj in self.trajectories[:ntraj] @@ -820,6 +825,11 @@ def e_data_traj_std(self, ntraj=-1): """ if not self.trajectories: return None + if ntraj > len(self.trajectories): + raise ValueError( + f"Cannot compute statistic for {ntraj} trajectories. " + f"Only {len(self.trajectories)} trajectories are stored." + ) return { k: np.std(np.stack([ traj.e_data[k] for traj in self.trajectories[:ntraj] @@ -916,10 +926,12 @@ def __add__(self, other): new.expect = new.average_expect new.e_data = new.average_e_data - new.std_e_data = { - k: np.sqrt(avg_expect2 - abs(avg_expect**2)) - for k, avg_expect, avg_expect2 in zip(self._raw_ops, avg, avg2) - } + new.std_e_data = {} + for i, key in enumerate(self._raw_ops): + std2 = avg2[i] - abs(avg[i]**2) + std2[std2 < 0] = 0. + new.std_e_data[key] = np.sqrt(std2) + new.std_expect = list(new.std_e_data.values()) if new.trajectories: diff --git a/qutip/solver/sesolve.py b/qutip/solver/sesolve.py index 2450fc652f..0914624232 100644 --- a/qutip/solver/sesolve.py +++ b/qutip/solver/sesolve.py @@ -2,7 +2,7 @@ This module provides solvers for the unitary Schrodinger equation. """ -__all__ = ['sesolve', 'SeSolver'] +__all__ = ['sesolve', 'SESolver'] import numpy as np from time import time @@ -99,11 +99,11 @@ def sesolve(H, psi0, tlist, e_ops=None, args=None, options=None): is an empty list of `store_states=True` in options]. """ H = QobjEvo(H, args=args, tlist=tlist) - solver = SeSolver(H, options=options) + solver = SESolver(H, options=options) return solver.run(psi0, tlist, e_ops=e_ops) -class SeSolver(Solver): +class SESolver(Solver): """ Schrodinger equation evolution of a state vector or unitary matrix for a given Hamiltonian. @@ -116,7 +116,7 @@ class SeSolver(Solver): that can be made into :class:`QobjEvo` are also accepted. options : dict, optional - Options for the solver, see :obj:`SeSolver.options` and + Options for the solver, see :obj:`SESolver.options` and `Integrator <./classes.html#classes-ode>`_ for a list of all options. attributes diff --git a/qutip/tests/solver/test_correlation.py b/qutip/tests/solver/test_correlation.py index 3e7f0f0820..d17c6e526b 100644 --- a/qutip/tests/solver/test_correlation.py +++ b/qutip/tests/solver/test_correlation.py @@ -286,7 +286,7 @@ def test_correlation_timedependant_op(): def test_alternative_solver(): - from qutip.solver.mesolve import MeSolver + from qutip.solver.mesolve import MESolver from qutip.solver.brmesolve import BRSolver H = qutip.num(5) @@ -294,7 +294,7 @@ def test_alternative_solver(): a_ops = [(a+a.dag(), qutip.coefficient(lambda _, w: w>0, args={"w":0}))] br = BRSolver(H, a_ops) - me = MeSolver(H, [a]) + me = MESolver(H, [a]) times = np.arange(4) br_corr = qutip.correlation_3op(br, qutip.basis(5), [0], times, a, a.dag()) diff --git a/qutip/tests/solver/test_integrator.py b/qutip/tests/solver/test_integrator.py index 3ba60add1c..424c39bd0f 100644 --- a/qutip/tests/solver/test_integrator.py +++ b/qutip/tests/solver/test_integrator.py @@ -1,5 +1,6 @@ -from qutip.solver.sesolve import SeSolver -from qutip.solver.mesolve import MeSolver +from qutip.solver.sesolve import SESolver +from qutip.solver.mesolve import MESolver +from qutip.solver.mcsolve import MCSolver from qutip.solver.solver_base import Solver from qutip.solver.ode.scipy_integrator import * import qutip @@ -15,21 +16,21 @@ class TestIntegratorCte(): me_system = qutip.liouvillian(qutip.QobjEvo(qutip.qeye(2)), c_ops=[qutip.destroy(2)]) - @pytest.fixture(params=list(SeSolver.avail_integrators().keys())) + @pytest.fixture(params=list(SESolver.avail_integrators().keys())) def se_method(self, request): return request.param - @pytest.fixture(params=list(MeSolver.avail_integrators().keys())) + @pytest.fixture(params=list(MESolver.avail_integrators().keys())) def me_method(self, request): return request.param - # TODO: Change when the McSolver is added - @pytest.fixture(params=list(Solver.avail_integrators().keys())) + # TODO: Change when the MCSolver is added + @pytest.fixture(params=list(MCSolver.avail_integrators().keys())) def mc_method(self, request): return request.param def test_se_integration(self, se_method): - evol = SeSolver.avail_integrators()[se_method](self.se_system, {}) + evol = SESolver.avail_integrators()[se_method](self.se_system, {}) state0 = qutip.core.unstack_columns(qutip.basis(6,0).data, (2, 3)) evol.set_state(0, state0) for t, state in evol.run(np.linspace(0, 2, 21)): @@ -38,7 +39,7 @@ def test_se_integration(self, se_method): assert state.shape == (2, 3) def test_me_integration(self, me_method): - evol = MeSolver.avail_integrators()[me_method](self.me_system, {}) + evol = MESolver.avail_integrators()[me_method](self.me_system, {}) state0 = qutip.operator_to_vector(qutip.fock_dm(2,1)).data evol.set_state(0, state0) for t in np.linspace(0, 2, 21): @@ -48,7 +49,7 @@ def test_me_integration(self, me_method): state.to_array()[0, 0], atol=2e-5) def test_mc_integration(self, mc_method): - evol = Solver.avail_integrators()[mc_method](self.se_system, {}) + evol = MCSolver.avail_integrators()[mc_method](self.se_system, {}) state = qutip.basis(2,0).data evol.set_state(0, state) t = 0 @@ -99,21 +100,21 @@ class TestIntegrator(TestIntegratorCte): ) @pytest.fixture( - params=[key for key, integrator in SeSolver.avail_integrators().items() + params=[key for key, integrator in SESolver.avail_integrators().items() if integrator.support_time_dependant] ) def se_method(self, request): return request.param @pytest.fixture( - params=[key for key, integrator in MeSolver.avail_integrators().items() + params=[key for key, integrator in MESolver.avail_integrators().items() if integrator.support_time_dependant] ) def me_method(self, request): return request.param @pytest.fixture( - params=[key for key, integrator in Solver.avail_integrators().items() + params=[key for key, integrator in MCSolver.avail_integrators().items() if integrator.support_time_dependant] ) def mc_method(self, request): diff --git a/qutip/tests/solver/test_mcsolve.py b/qutip/tests/solver/test_mcsolve.py index ec2be2e8b8..fdac5eb55b 100644 --- a/qutip/tests/solver/test_mcsolve.py +++ b/qutip/tests/solver/test_mcsolve.py @@ -2,7 +2,7 @@ import numpy as np import qutip from copy import copy -from qutip.solver.mcsolve import mcsolve, McSolver +from qutip.solver.mcsolve import mcsolve, MCSolver from qutip.solver.solver_base import Solver def _return_constant(t, args): @@ -341,7 +341,7 @@ def test_stepping(self): size = 10 a = qutip.QobjEvo([qutip.destroy(size), 'alpha'], args={'alpha': 0}) H = qutip.num(size) - mcsolver = McSolver(H, a, options={'map': 'serial'}) + mcsolver = MCSolver(H, a, options={'map': 'serial'}) mcsolver.start(qutip.basis(size, size-1), 0, seed=5) state_1 = mcsolver.step(1, args={'alpha':1}) @@ -389,7 +389,7 @@ def test_McSolver_run(): size = 10 a = qutip.QobjEvo([qutip.destroy(size), 'coupling'], args={'coupling':0}) H = qutip.num(size) - solver = McSolver(H, a) + solver = MCSolver(H, a) solver.options = {'store_final_state': True} res = solver.run(qutip.basis(size, size-1), np.linspace(0, 5.0, 11), e_ops=[qutip.qeye(size)], args={'coupling': 1}) @@ -409,7 +409,7 @@ def test_McSolver_stepping(): size = 10 a = qutip.QobjEvo([qutip.destroy(size), 'coupling'], args={'coupling':0}) H = qutip.num(size) - solver = McSolver(H, a) + solver = MCSolver(H, a) solver.start(qutip.basis(size, size-1), 0, seed=0) solver.options = {'method': 'lsoda'} state = solver.step(1) diff --git a/qutip/tests/solver/test_mesolve.py b/qutip/tests/solver/test_mesolve.py index 90bfe615d7..6f6bbab6c7 100644 --- a/qutip/tests/solver/test_mesolve.py +++ b/qutip/tests/solver/test_mesolve.py @@ -2,13 +2,13 @@ from types import FunctionType import qutip -from qutip.solver.mesolve import mesolve, MeSolver +from qutip.solver.mesolve import mesolve, MESolver from qutip.solver.solver_base import Solver import pickle import pytest all_ode_method = [ - method for method, integrator in MeSolver.avail_integrators().items() + method for method, integrator in MESolver.avail_integrators().items() if integrator.support_time_dependant ] @@ -222,7 +222,7 @@ def test_mesolve_normalization(self, state_type): def test_mesolver_pickling(self): options = {"progress_bar": None} - solver_obj = MeSolver(self.ada, c_ops=[self.a], options=options) + solver_obj = MESolver(self.ada, c_ops=[self.a], options=options) solver_copy = pickle.loads(pickle.dumps(solver_obj)) e1 = solver_obj.run(qutip.basis(self.N, 9), [0, 1, 2, 3], e_ops=[self.ada]).expect[0] @@ -234,7 +234,7 @@ def test_mesolver_pickling(self): all_ode_method, ids=all_ode_method) def test_mesolver_stepping(self, method): options = {"method": method, "progress_bar": None} - solver_obj = MeSolver( + solver_obj = MESolver( self.ada, c_ops=qutip.QobjEvo( [self.a, lambda t, kappa: np.sqrt(kappa * np.exp(-t))], @@ -661,17 +661,17 @@ def test_num_collapse_set(): def test_mesolve_bad_H(): with pytest.raises(TypeError): - MeSolver([qutip.qeye(3), 't'], []) + MESolver([qutip.qeye(3), 't'], []) with pytest.raises(TypeError): - MeSolver(qutip.qeye(3), [[qutip.qeye(3), 't']]) + MESolver(qutip.qeye(3), [[qutip.qeye(3), 't']]) def test_mesolve_bad_state(): - solver = MeSolver(qutip.qeye(4), []) + solver = MESolver(qutip.qeye(4), []) with pytest.raises(TypeError): solver.start(qutip.basis(2,1) & qutip.basis(2,0), 0) def test_mesolve_bad_options(): with pytest.raises(TypeError): - MeSolver(qutip.qeye(4), [], options=False) + MESolver(qutip.qeye(4), [], options=False) diff --git a/qutip/tests/solver/test_options.py b/qutip/tests/solver/test_options.py index 79d16e1d0b..d7758a1a71 100644 --- a/qutip/tests/solver/test_options.py +++ b/qutip/tests/solver/test_options.py @@ -89,7 +89,7 @@ def test_print(): def test_in_solver(): opt = {"method": "adams", "store_states": True, "atol": 1} - solver = qutip.solver.sesolve.SeSolver(qutip.qeye(1), options=opt) + solver = qutip.solver.sesolve.SESolver(qutip.qeye(1), options=opt) adams = qutip.solver.ode.scipy_integrator.IntegratorScipyAdams lsoda = qutip.solver.ode.scipy_integrator.IntegratorScipylsoda bdf = qutip.solver.ode.scipy_integrator.IntegratorScipyBDF @@ -119,7 +119,7 @@ def test_in_solver(): def test_options_update_solver(): opt = {"method": "adams", "normalize_output": False} - solver = qutip.solver.sesolve.SeSolver(1j * qutip.qeye(1), options=opt) + solver = qutip.solver.sesolve.SESolver(1j * qutip.qeye(1), options=opt) solver.start(qutip.basis(1), 0) err_atol_def = (solver.step(1) - np.exp(1)).norm() diff --git a/qutip/tests/solver/test_propagator.py b/qutip/tests/solver/test_propagator.py index 72f668cdc8..8970022629 100644 --- a/qutip/tests/solver/test_propagator.py +++ b/qutip/tests/solver/test_propagator.py @@ -4,9 +4,9 @@ import qutip import pytest from qutip.solver.brmesolve import BRSolver -from qutip.solver.mesolve import MeSolver -from qutip.solver.sesolve import SeSolver -from qutip.solver.mcsolve import McSolver +from qutip.solver.mesolve import MESolver +from qutip.solver.sesolve import SESolver +from qutip.solver.mcsolve import MCSolver def testPropHOB(): @@ -71,7 +71,7 @@ def testPropHOSteady(): kappa = 0.1 n_th = 2 rate = kappa * (1 + n_th) - c_op_list.append(np.sqrt(rate) * a) + c_op_list.append(np.sqrt(rate) * a)e rate = kappa * n_th c_op_list.append(np.sqrt(rate) * a.dag()) U = propagator(H, 2*np.pi, c_op_list) @@ -101,11 +101,11 @@ def testPropEvo(): def _make_se(H, a): - return SeSolver(H) + return SESolver(H) def _make_me(H, a): - return MeSolver(H, [a]) + return MESolver(H, [a]) def _make_br(H, a): @@ -114,8 +114,8 @@ def _make_br(H, a): @pytest.mark.parametrize('solver', [ - pytest.param(_make_se, id='SeSolver'), - pytest.param(_make_me, id='MeSolver'), + pytest.param(_make_se, id='SESolver'), + pytest.param(_make_me, id='MESolver'), pytest.param(_make_br, id='BRSolver'), ]) def testPropSolver(solver): @@ -134,7 +134,7 @@ def testPropSolver(solver): def testPropMcSolver(): a = destroy(5) H = a.dag()*a - solver = McSolver(H, [a]) + solver = MCSolver(H, [a]) with pytest.raises(TypeError) as err: Propagator(solver) assert str(err.value).startswith("Non-deterministic") diff --git a/qutip/tests/solver/test_sesolve.py b/qutip/tests/solver/test_sesolve.py index 135be222bf..b5a44e58ba 100644 --- a/qutip/tests/solver/test_sesolve.py +++ b/qutip/tests/solver/test_sesolve.py @@ -2,11 +2,11 @@ import pickle import qutip import numpy as np -from qutip.solver.sesolve import sesolve, SeSolver +from qutip.solver.sesolve import sesolve, SESolver from qutip.solver.solver_base import Solver all_ode_method = [ - method for method, integrator in SeSolver.avail_integrators().items() + method for method, integrator in SESolver.avail_integrators().items() if integrator.support_time_dependant ] @@ -200,7 +200,7 @@ def test_compare_evolution(self, H, normalize, args, tol=5e-5): def test_sesolver_args(self): options = {"progress_bar": None} - solver_obj = SeSolver(qutip.QobjEvo([self.H0, [self.H1,'a']], + solver_obj = SESolver(qutip.QobjEvo([self.H0, [self.H1,'a']], args={'a': 1}), options=options) res = solver_obj.run(qutip.basis(2,1), [0, 1, 2, 3], @@ -210,7 +210,7 @@ def test_sesolver_args(self): def test_sesolver_pickling(self): e_ops = [qutip.sigmax(), qutip.sigmay(), qutip.sigmaz()] options = {"progress_bar": None} - solver_obj = SeSolver(self.H0 + self.H1, + solver_obj = SESolver(self.H0 + self.H1, options=options) solver_copy = pickle.loads(pickle.dumps(solver_obj)) sx, sy, sz = solver_obj.run(qutip.basis(2,1), [0, 1, 2, 3], @@ -229,7 +229,7 @@ def test_sesolver_stepping(self, method): "rtol": 1e-8, "progress_bar": None } - solver_obj = SeSolver( + solver_obj = SESolver( qutip.QobjEvo([self.H1, lambda t, a: a], args={"a":0.25}), options=options ) @@ -268,13 +268,13 @@ def test_sesolver_stepping(self, method): def test_sesolve_bad_H(): with pytest.raises(TypeError): - SeSolver(np.eye(3)) + SESolver(np.eye(3)) with pytest.raises(ValueError): - SeSolver(qutip.basis(3,1)) + SESolver(qutip.basis(3,1)) def test_sesolve_bad_state(): - solver = SeSolver(qutip.qeye(4)) + solver = SESolver(qutip.qeye(4)) with pytest.raises(TypeError): solver.start(qutip.basis(4,1).dag(), 0) with pytest.raises(TypeError): @@ -282,6 +282,6 @@ def test_sesolve_bad_state(): def test_sesolve_step_no_start(): - solver = SeSolver(qutip.qeye(4)) + solver = SESolver(qutip.qeye(4)) with pytest.raises(RuntimeError): solver.step(1) From 0cfab4e1f8f35fbbb2d747d0ac6e913d7f00dd68 Mon Sep 17 00:00:00 2001 From: Eric Giguere Date: Thu, 8 Dec 2022 16:59:50 -0500 Subject: [PATCH 2/8] Fix doc build --- doc/apidoc/functions.rst | 6 +----- doc/guide/dynamics/dynamics-monte.rst | 25 ++++++++++++++----------- qutip/solver/mcsolve.py | 1 + 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/doc/apidoc/functions.rst b/doc/apidoc/functions.rst index fd359553b5..13d2377376 100644 --- a/doc/apidoc/functions.rst +++ b/doc/apidoc/functions.rst @@ -141,13 +141,9 @@ Master Equation Monte Carlo Evolution --------------------- -.. automodule:: qutip.solve.mcsolve +.. automodule:: qutip.solver.mcsolve :members: mcsolve -.. ignore f90 stuff for now - .. automodule:: qutip.fortran.mcsolve_f90 - :members: mcsolve_f90 - Krylov Subspace Solver ---------------------- diff --git a/doc/guide/dynamics/dynamics-monte.rst b/doc/guide/dynamics/dynamics-monte.rst index c1f528e3e5..2e88aa0de8 100644 --- a/doc/guide/dynamics/dynamics-monte.rst +++ b/doc/guide/dynamics/dynamics-monte.rst @@ -4,10 +4,6 @@ Monte Carlo Solver ******************************************* -.. plot:: - :include-source: False - - from qutip.solver.mcsolve import * .. _monte-intro: @@ -75,6 +71,8 @@ To illustrate the use of the Monte Carlo evolution of quantum systems in QuTiP, .. plot:: :context: + from qutip.solver.mcsolve import MCSolver, mcsolve + times = np.linspace(0.0, 10.0, 200) psi0 = tensor(fock(2, 0), fock(10, 5)) a = tensor(qeye(2), destroy(10)) @@ -110,7 +108,7 @@ If we want to run 1000 trajectories in the above example, we can simply modify t .. plot:: :context: close-figs - data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], [a.dag() * a, sm.dag() * sm], ntraj=1000) + data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], e_ops=[a.dag() * a, sm.dag() * sm], ntraj=1000) where we have added the keyword argument ``ntraj=1000`` at the end of the inputs. Now, the Monte Carlo solver will calculate expectation values for both operators, ``a.dag() * a, sm.dag() * sm`` averaging over 1000 trajectories. @@ -132,11 +130,11 @@ we can extract the relevant expectation values using: expt100 = data.expect_traj_avg(100) expt1000 = data.expect_traj_avg(1000) - plt.figure() - plt.plot(times, expt1, label="ntraj=1") - plt.plot(times, expt10, label="ntraj=10") - plt.plot(times, expt100, label="ntraj=100") - plt.plot(times, expt1000, label="ntraj=1000") + plt.figure() + plt.plot(times, expt1[0], label="ntraj=1") + plt.plot(times, expt10[0], label="ntraj=10") + plt.plot(times, expt100[0], label="ntraj=100") + plt.plot(times, expt1000[0], label="ntraj=1000") plt.title('Monte Carlo time evolution') plt.xlabel('Time') plt.ylabel('Expectation values') @@ -144,6 +142,11 @@ we can extract the relevant expectation values using: plt.show() +Monte Carlo Solver Result +------------------------- + +... + .. _monte-reuse: Reusing Hamiltonian Data @@ -196,7 +199,7 @@ In addition to the initial state, one may reuse the Hamiltonian data when adding sm = tensor(destroy(2), qeye(10)) H = 2*np.pi*a.dag()*a + 2*np.pi*sm.dag()*sm + 2*np.pi*0.25*(sm*a.dag() + sm.dag()*a) - solver = MCSolver(H, c_ops=[np.sqrt(0.1) * a]) + # solver = MCSolver(H, c_ops=[np.sqrt(0.1) * a]) data1 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=10) data2 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=10) data_merged = data1 + data2 diff --git a/qutip/solver/mcsolve.py b/qutip/solver/mcsolve.py index c9669e2ee1..24458ce5dc 100644 --- a/qutip/solver/mcsolve.py +++ b/qutip/solver/mcsolve.py @@ -103,6 +103,7 @@ def mcsolve(H, state, tlist, c_ops=(), e_ops=None, ntraj=500, *, Seed for the random number generator. It can be a single seed used to spawn seeds for each trajectory or a list of seeds, one for each trajectory. Seeds are saved in the result and they can be reused with:: + seeds=prev_result.seeds target_tol : {float, tuple, list}, optional From 377ba3de98dfc118d7dde7ac2f8cdf57256c46f6 Mon Sep 17 00:00:00 2001 From: Eric Date: Fri, 9 Dec 2022 10:39:46 -0500 Subject: [PATCH 3/8] Fix typo --- qutip/tests/solver/test_propagator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qutip/tests/solver/test_propagator.py b/qutip/tests/solver/test_propagator.py index 8970022629..6cb3ac598b 100644 --- a/qutip/tests/solver/test_propagator.py +++ b/qutip/tests/solver/test_propagator.py @@ -71,7 +71,7 @@ def testPropHOSteady(): kappa = 0.1 n_th = 2 rate = kappa * (1 + n_th) - c_op_list.append(np.sqrt(rate) * a)e + c_op_list.append(np.sqrt(rate) * a) rate = kappa * n_th c_op_list.append(np.sqrt(rate) * a.dag()) U = propagator(H, 2*np.pi, c_op_list) From c75ab06709842b47dc01239b5ccc2fc10b0c8b40 Mon Sep 17 00:00:00 2001 From: Eric Date: Fri, 9 Dec 2022 14:36:38 -0500 Subject: [PATCH 4/8] Add mcsolve's result explanations --- doc/guide/dynamics/dynamics-monte.rst | 57 ++++++++++++++++++--------- qutip/solver/result.py | 2 +- 2 files changed, 40 insertions(+), 19 deletions(-) diff --git a/doc/guide/dynamics/dynamics-monte.rst b/doc/guide/dynamics/dynamics-monte.rst index 2e88aa0de8..47f299db8f 100644 --- a/doc/guide/dynamics/dynamics-monte.rst +++ b/doc/guide/dynamics/dynamics-monte.rst @@ -74,7 +74,7 @@ To illustrate the use of the Monte Carlo evolution of quantum systems in QuTiP, from qutip.solver.mcsolve import MCSolver, mcsolve times = np.linspace(0.0, 10.0, 200) - psi0 = tensor(fock(2, 0), fock(10, 5)) + psi0 = tensor(fock(2, 0), fock(10, 8)) a = tensor(qeye(2), destroy(10)) sm = tensor(destroy(2), qeye(10)) H = 2*np.pi*a.dag()*a + 2*np.pi*sm.dag()*sm + 2*np.pi*0.25*(sm*a.dag() + sm.dag()*a) @@ -92,7 +92,35 @@ To illustrate the use of the Monte Carlo evolution of quantum systems in QuTiP, The advantage of the Monte Carlo method over the master equation approach is that only the state vector is required to be kept in the computers memory, as opposed to the entire density matrix. For large quantum system this becomes a significant advantage, and the Monte Carlo solver is therefore generally recommended for such systems. For example, simulating a Heisenberg spin-chain consisting of 10 spins with random parameters and initial states takes almost 7 times longer using the master equation rather than Monte Carlo approach with the default number of trajectories running on a quad-CPU machine. Furthermore, it takes about 7 times the memory as well. However, for small systems, the added overhead of averaging a large number of stochastic trajectories to obtain the open system dynamics, as well as starting the multiprocessing functionality, outweighs the benefit of the minor (in this case) memory saving. Master equation methods are therefore generally more efficient when Hilbert space sizes are on the order of a couple of hundred states or smaller. -Like the master equation solver :func:`qutip.mesolve`, the Monte Carlo solver returns a :class:`qutip.Result` object consisting of expectation values, if the user has defined expectation value operators in the 5th argument to ``mcsolve``, or state vectors if no expectation value operators are given. If state vectors are returned, then the :class:`qutip.Result` returned by :func:`qutip.mcsolve` will be an array of length ``ntraj``, with each element containing an array of ket-type qobjs with the same number of elements as ``times``. Furthermore, the output :class:`qutip.Result` object will also contain a list of times at which collapse occurred, and which collapse operators did the collapse, in the ``col_times`` and ``col_which`` properties, respectively. + +Monte Carlo Solver Result +------------------------- + +The Monte Carlo solver returns a :class:`qutip.MultitrajResult` object consisting of expectation values and/or states. +The main difference with :func:`qutip.mesolve`'s :class:`qutip.Result` is that it optionally store the result of each trajectories with the averages. +When trajectories are stored, ``result.runs_expect`` is a list over the expectation operators, trajectories and times in that order. +The averages are stored in ``result.average_expect`` and the standard derivation of the expectation values in ``result.std_expect``. +When states are returned, ``result.runs_states`` will be an array of length ``ntraj``, with each element containing an array of ket-type qobjs with the same number of elements as ``times`` and ``result.average_states`` is a list of density matrices for each times. +Furthermore, the output will also contain a list of times at which collapse occurred, and which collapse operators did the collapse, in the ``col_times`` and ``col_which`` properties, respectively. +Lastly ``result.photocurrent`` contain the measurement of the evolution. + + +.. plot:: + :context: close-figs + + times = np.linspace(0.0, 10.0, 200) + psi0 = tensor(fock(2, 0), fock(10, 8)) + a = tensor(qeye(2), destroy(10)) + sm = tensor(destroy(2), qeye(10)) + H = 2*np.pi*a.dag()*a + 2*np.pi*sm.dag()*sm + 2*np.pi*0.25*(sm*a.dag() + sm.dag()*a) + data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], e_ops=[a.dag() * a, sm.dag() * sm]) + + plt.figure() + plt.plot((times[:-1] + times[1:])/2, data.photocurrent[0]) + plt.title('Monte Carlo Photocurrent') + plt.xlabel('Time') + plt.ylabel('Photon detections') + plt.show() .. _monte-ntraj: @@ -113,12 +141,12 @@ If we want to run 1000 trajectories in the above example, we can simply modify t where we have added the keyword argument ``ntraj=1000`` at the end of the inputs. Now, the Monte Carlo solver will calculate expectation values for both operators, ``a.dag() * a, sm.dag() * sm`` averaging over 1000 trajectories. Sometimes one is also interested in seeing how the Monte Carlo trajectories converge to the master equation solution by calculating expectation values over a range of trajectory numbers. -If, for example, we want to average over 1, 10, 100, and 1000 trajectories, we must first run the solver keeping expectation values for each trajectories: +If, for example, we want to average over 1, 10 and 100 trajectories, we must first run the solver keeping expectation values for each trajectories: .. plot:: :context: - data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], e_ops=[a.dag() * a, sm.dag() * sm], options={"keep_runs_results": True}, ntraj=1000) + data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], e_ops=[a.dag() * a, sm.dag() * sm], options={"keep_runs_results": True}, ntraj=100, seeds=1) we can extract the relevant expectation values using: @@ -128,13 +156,11 @@ we can extract the relevant expectation values using: expt1 = data.expect_traj_avg(1) expt10 = data.expect_traj_avg(10) expt100 = data.expect_traj_avg(100) - expt1000 = data.expect_traj_avg(1000) plt.figure() plt.plot(times, expt1[0], label="ntraj=1") plt.plot(times, expt10[0], label="ntraj=10") plt.plot(times, expt100[0], label="ntraj=100") - plt.plot(times, expt1000[0], label="ntraj=1000") plt.title('Monte Carlo time evolution') plt.xlabel('Time') plt.ylabel('Expectation values') @@ -142,11 +168,6 @@ we can extract the relevant expectation values using: plt.show() -Monte Carlo Solver Result -------------------------- - -... - .. _monte-reuse: Reusing Hamiltonian Data @@ -172,13 +193,13 @@ Using the previous example, we will calculate the dynamics for two different ini H = 2*np.pi*a.dag()*a + 2*np.pi*sm.dag()*sm + 2*np.pi*0.25*(sm*a.dag() + sm.dag()*a) solver = MCSolver(H, c_ops=[np.sqrt(0.1) * a]) - data1 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=1000) + data1 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=100) psi1 = tensor(fock(2, 0), coherent(10, 2 - 1j)) - data2 = solver.run(psi1, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=1000) + data2 = solver.run(psi1, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=100) plt.figure() - plt.plot(times, data1.expect[0], times, data1.expect[1], lw=2) - plt.plot(times, data2.expect[0], '--', times, data2.expect[1], '--', lw=2) + plt.plot(times, data1.expect[0], "b", times, data1.expect[1], "r", lw=2) + plt.plot(times, data2.expect[0], 'b--', times, data2.expect[1], 'r--', lw=2) plt.title('Monte Carlo time evolution') plt.xlabel('Time', fontsize=14) plt.ylabel('Expectation values', fontsize=14) @@ -199,9 +220,9 @@ In addition to the initial state, one may reuse the Hamiltonian data when adding sm = tensor(destroy(2), qeye(10)) H = 2*np.pi*a.dag()*a + 2*np.pi*sm.dag()*sm + 2*np.pi*0.25*(sm*a.dag() + sm.dag()*a) - # solver = MCSolver(H, c_ops=[np.sqrt(0.1) * a]) - data1 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=10) - data2 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=10) + solver = MCSolver(H, c_ops=[np.sqrt(0.1) * a]) + data1 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=1, seed=1) + data2 = solver.run(psi0, times, e_ops=[a.dag() * a, sm.dag() * sm], ntraj=1, seed=3) data_merged = data1 + data2 plt.figure() diff --git a/qutip/solver/result.py b/qutip/solver/result.py index 24c5a50562..48d7dd12c1 100644 --- a/qutip/solver/result.py +++ b/qutip/solver/result.py @@ -376,7 +376,7 @@ class MultiTrajResult(_BaseResult): The final state (if the recording of the final state was requested) averaged over all trajectories as a density matrix. - runs_state : list of :obj:`~Qobj` + runs_final_state : list of :obj:`~Qobj` The final state for each trajectory (if the recording of the final state and trajectories was requested). From c9d66d6892823e8676d88f9bf9ab6cc31ab3e39d Mon Sep 17 00:00:00 2001 From: Eric Date: Fri, 9 Dec 2022 16:30:58 -0500 Subject: [PATCH 5/8] Change link location --- doc/guide/dynamics/dynamics-master.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/guide/dynamics/dynamics-master.rst b/doc/guide/dynamics/dynamics-master.rst index 5384ee4d6f..932398efe0 100644 --- a/doc/guide/dynamics/dynamics-master.rst +++ b/doc/guide/dynamics/dynamics-master.rst @@ -72,7 +72,7 @@ The resulting list of expectation values can easily be visualized using matplotl >>> ax.legend(("Sigma-Z", "Sigma-Y")) # doctest: +SKIP >>> plt.show() # doctest: +SKIP -If an empty list of operators is passed as fifth parameter, the :func:`qutip.mesolve` function returns a :class:`qutip.solve.solver.Result` instance that contains a list of state vectors for the times specified in ``times`` +If an empty list of operators is passed as fifth parameter, the :func:`qutip.mesolve` function returns a :class:`qutip.Result` instance that contains a list of state vectors for the times specified in ``times`` .. plot:: :context: close-figs From e93975168770cf81cad8caed0ea3dad98e2a6be9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Eric=20Gigu=C3=A8re?= Date: Tue, 13 Dec 2022 09:17:42 -0500 Subject: [PATCH 6/8] Apply suggestions from code review Co-authored-by: Asier Galicia <57414022+AGaliciaMartinez@users.noreply.github.com> --- doc/guide/dynamics/dynamics-monte.rst | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/doc/guide/dynamics/dynamics-monte.rst b/doc/guide/dynamics/dynamics-monte.rst index 47f299db8f..7ebf030cd3 100644 --- a/doc/guide/dynamics/dynamics-monte.rst +++ b/doc/guide/dynamics/dynamics-monte.rst @@ -97,11 +97,11 @@ Monte Carlo Solver Result ------------------------- The Monte Carlo solver returns a :class:`qutip.MultitrajResult` object consisting of expectation values and/or states. -The main difference with :func:`qutip.mesolve`'s :class:`qutip.Result` is that it optionally store the result of each trajectories with the averages. +The main difference with :func:`qutip.mesolve`'s :class:`qutip.Result` is that it optionally stores the result of each trajectory together with their averages. When trajectories are stored, ``result.runs_expect`` is a list over the expectation operators, trajectories and times in that order. The averages are stored in ``result.average_expect`` and the standard derivation of the expectation values in ``result.std_expect``. -When states are returned, ``result.runs_states`` will be an array of length ``ntraj``, with each element containing an array of ket-type qobjs with the same number of elements as ``times`` and ``result.average_states`` is a list of density matrices for each times. -Furthermore, the output will also contain a list of times at which collapse occurred, and which collapse operators did the collapse, in the ``col_times`` and ``col_which`` properties, respectively. +When the states are returned, ``result.runs_states`` will be an array of length ``ntraj``. Each element contains an array of "Qobj" type ket with the same number of elements as ``times``. ``result.average_states`` is a list of density matrices computed as the average of the states at each time step. +Furthermore, the output will also contain a list of times at which the collapse occurred, and which collapse operators did the collapse. These can be obtained in ``result.col_times`` and ``result.col_which`` respectively. Lastly ``result.photocurrent`` contain the measurement of the evolution. @@ -128,10 +128,10 @@ Lastly ``result.photocurrent`` contain the measurement of the evolution. Changing the Number of Trajectories ----------------------------------- -As mentioned earlier, by default, the ``mcsolve`` function runs 500 trajectories. +By default, the ``mcsolve`` function runs 500 trajectories. This value was chosen because it gives good accuracy, Monte Carlo errors scale as :math:`1/n` where :math:`n` is the number of trajectories, and simultaneously does not take an excessive amount of time to run. -However, like many other options in QuTiP you are free to change the number of trajectories to fit your needs. -If we want to run 1000 trajectories in the above example, we can simply modify the call to ``mcsolve`` like: +However, you can change the number of trajectories to fit your needs. +In order to run 1000 trajectories in the above example, we can simply modify the call to ``mcsolve`` like: .. plot:: :context: close-figs @@ -140,8 +140,7 @@ If we want to run 1000 trajectories in the above example, we can simply modify t where we have added the keyword argument ``ntraj=1000`` at the end of the inputs. Now, the Monte Carlo solver will calculate expectation values for both operators, ``a.dag() * a, sm.dag() * sm`` averaging over 1000 trajectories. -Sometimes one is also interested in seeing how the Monte Carlo trajectories converge to the master equation solution by calculating expectation values over a range of trajectory numbers. -If, for example, we want to average over 1, 10 and 100 trajectories, we must first run the solver keeping expectation values for each trajectories: +In order to explore the convergence of the Monte Carlo solver, it is possible to retrieve expectation values as a function of the number of trajectories. For example, the following code block plots expectation values for 1, 10 and 100 trajectories by first running the solver for 100 trajectories but keeping the expectation values for each trajectory: .. plot:: :context: @@ -208,8 +207,7 @@ Using the previous example, we will calculate the dynamics for two different ini .. guide-dynamics-mc2: -In addition to the initial state, one may reuse the Hamiltonian data when adding more trajectories to the result or changing simulation times ``times``. - +The ``MCSolver`` also allows adding new trajectories after the first computation. This is shown in the next example where the results of two separated runs with identical conditions are merged into a single ``result`` object. .. plot:: :context: close-figs From 6b2479cfe11a255064717991d451ed8ec6ab3c16 Mon Sep 17 00:00:00 2001 From: Eric Giguere Date: Tue, 13 Dec 2022 09:52:50 -0500 Subject: [PATCH 7/8] Computing statistic when unavailable raise error --- qutip/solver/result.py | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/qutip/solver/result.py b/qutip/solver/result.py index 48d7dd12c1..42383501f5 100644 --- a/qutip/solver/result.py +++ b/qutip/solver/result.py @@ -784,7 +784,10 @@ def e_data_traj_avg(self, ntraj=-1): Default: all trajectories. """ if not self.trajectories: - return None + raise ValueError( + f"Trajectories information are not available. " + f"Use the options 'keep_runs_results' to store trajectories." + ) if ntraj > len(self.trajectories): raise ValueError( f"Cannot compute statistic for {ntraj} trajectories. " @@ -809,7 +812,10 @@ def expect_traj_avg(self, ntraj=-1): Default: all trajectories. """ if not self.trajectories: - return None + raise ValueError( + f"Trajectories information are not available. " + f"Use the options 'keep_runs_results' to store trajectories." + ) return list(self.e_data_traj_avg(ntraj).values()) def e_data_traj_std(self, ntraj=-1): @@ -824,7 +830,10 @@ def e_data_traj_std(self, ntraj=-1): Default: all trajectories. """ if not self.trajectories: - return None + raise ValueError( + f"Trajectories information are not available. " + f"Use the options 'keep_runs_results' to store trajectories." + ) if ntraj > len(self.trajectories): raise ValueError( f"Cannot compute statistic for {ntraj} trajectories. " @@ -849,7 +858,10 @@ def expect_traj_std(self, ntraj=-1): Default: all trajectories. """ if not self.trajectories: - return None + raise ValueError( + f"Trajectories information are not available. " + f"Use the options 'keep_runs_results' to store trajectories." + ) return list(self.e_data_traj_std(ntraj).values()) def __repr__(self): From 551768aa9677df6fd9cbf1d2f086ff678eef1c50 Mon Sep 17 00:00:00 2001 From: Eric Giguere Date: Tue, 13 Dec 2022 10:50:47 -0500 Subject: [PATCH 8/8] Fix tests --- qutip/tests/solver/test_results.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/qutip/tests/solver/test_results.py b/qutip/tests/solver/test_results.py index ea5de3f801..b11b19fd47 100644 --- a/qutip/tests/solver/test_results.py +++ b/qutip/tests/solver/test_results.py @@ -202,10 +202,14 @@ def _expect_check_types(self, multiresult): else: assert multiresult.runs_expect == [] assert multiresult.runs_e_data == {} - assert multiresult.expect_traj_avg() is None - assert multiresult.expect_traj_std() is None - assert multiresult.e_data_traj_avg() is None - assert multiresult.e_data_traj_std() is None + with pytest.raises(ValueError): + multiresult.expect_traj_avg() + with pytest.raises(ValueError): + multiresult.expect_traj_std() + with pytest.raises(ValueError): + multiresult.e_data_traj_avg() + with pytest.raises(ValueError): + multiresult.e_data_traj_std() @pytest.mark.parametrize('keep_runs_results', [True, False]) @pytest.mark.parametrize('dm', [True, False])