Merge pull request #2042 from Ericgig/solver.doc

Documentation for v5 mcsolve

doc/apidoc/classes.rst

@@ -44,10 +44,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

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
 ----------------------

doc/guide/dynamics/dynamics-master.rst

@@ -35,7 +35,7 @@ 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()])
 
 
 See the next section for examples on how dissipation is included by defining a list of collapse operators and using :func:`qutip.mesolve` instead.
@@ -49,7 +49,7 @@ Adding operators to this list results in a larger output list returned by the fu
 .. 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,
@@ -78,7 +78,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.Result` instance that contains a list of state vectors for the times specified in ``times``
+If an empty list of operators is passed to the ``e_ops`` parameter, the :func:`qutip.sesolve` and :func:`qutip.mesolve` functions return a :class:`qutip.Result` instance that contains a list of state vectors for the times specified in ``times``
 
 .. plot::
     :context: close-figs
@@ -176,7 +176,7 @@ operators ``[sigmaz(), sigmay()]`` to the fifth argument.
     :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
@@ -199,7 +199,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

doc/guide/dynamics/dynamics-monte.rst

@@ -4,6 +4,7 @@
 Monte Carlo Solver
 *******************************************
 
+
 .. _monte-intro:
 
 Introduction
@@ -70,12 +71,14 @@ 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))
+    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], [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,44 +92,79 @@ 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.
 
+Monte Carlo Solver Result
+-------------------------
 
-.. _monte-ntraj:
-
-Changing the Number of Trajectories
------------------------------------
+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 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 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.
 
-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)
+    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])
 
-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:
+    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:
+
+Changing the Number of 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, 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:
+    :context: close-figs
 
-    ntraj = [1, 10, 100, 1000]
+    data = mcsolve(H, psi0, times, [np.sqrt(0.1) * a], e_ops=[a.dag() * a, sm.dag() * sm], ntraj=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.
+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:
 
-    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=100, seeds=1)
 
 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)
+
+    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.title('Monte Carlo time evolution')
+    plt.xlabel('Time')
+    plt.ylabel('Expectation values')
+    plt.legend()
+    plt.show()
 
 
 .. _monte-reuse:
@@ -136,27 +174,31 @@ 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))
     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)
-    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=100)
     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=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)
@@ -165,4 +207,28 @@ 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.
+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
+
+    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=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()
+    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()

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)

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):