Stochastic solver using v5 class interface

.coveragerc

@@ -3,8 +3,12 @@ source = qutip
 omit =
     # QuTiP test files
     */qutip/tests/*
+    # Tool for tests
+    */qutip/solver/sode/_noise.py
 
 [report]
 exclude_lines =
     # Skip Python wrappers which help load in C extension modules.
     __bootstrap__()
+    # Skip empty method that are never meant to be used.
+    raise NotImplementedError

doc/apidoc/classes.rst

@@ -53,6 +53,13 @@ Solver
 .. autoclass:: qutip.solver.brmesolve.BRSolver
     :members:
 
+.. autoclass:: qutip.solver.stochastic.SMESolver
+    :members:
+
+.. autoclass:: qutip.solver.stochastic.SSESolver
+    :members:
+
+
 
 .. _classes-non_markov_heom:
 
@@ -132,6 +139,39 @@ Integrator
     :members: options
 
 
+.. _classes-sode:
+
+Stochastic Integrator
+---------------------
+
+.. autoclass:: qutip.solver.sode.rouchon.RouchonSODE
+    :members: options
+
+.. autoclass:: qutip.solver.sode.itotaylor.EulerSODE
+    :members: options
+
+.. autoclass:: qutip.solver.sode.itotaylor.Milstein_SODE
+    :members: options
+
+.. autoclass:: qutip.solver.sode.itotaylor.Taylor1_5_SODE
+    :members: options
+
+.. autoclass:: qutip.solver.sode.itotaylor.Implicit_Milstein_SODE
+    :members: options
+
+.. autoclass:: qutip.solver.sode.itotaylor.Implicit_Taylor1_5_SODE
+    :members: options
+
+.. autoclass:: qutip.solver.sode.sode.PlatenSODE
+    :members: options
+
+.. autoclass:: qutip.solver.sode.itotaylor.Explicit1_5_SODE
+    :members: options
+
+.. autoclass:: qutip.solver.sode.sode.PredCorr_SODE
+    :members: options
+
+
 .. _classes-non_markov_mc_and_tt:
 
 Non-Markovian Memory Cascade and Transfer Tensor Solvers

doc/apidoc/functions.rst

@@ -169,8 +169,8 @@ Floquet States and Floquet-Markov Master Equation
 Stochastic Schrödinger Equation and Master Equation
 ---------------------------------------------------
 
-.. automodule:: qutip.solve.stochastic
-    :members: ssesolve, photocurrent_sesolve, smepdpsolve, smesolve, photocurrent_mesolve, ssepdpsolve, stochastic_solvers, general_stochastic
+.. automodule:: qutip.solver.stochastic
+    :members: ssesolve, smesolve
 
 
 Hierarchical Equations of Motion

doc/changes/2131.feature

@@ -0,0 +1 @@
+Rewrite the stochastic solver to use the v5 solver interface.

doc/guide/dynamics/dynamics-monte.rst

@@ -102,7 +102,12 @@ When trajectories are stored, ``result.runs_expect`` is a list over the expectat
 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.
+
+
+Photocurrent
+------------
+
+The photocurrent, previously computed using the ``photocurrent_sesolve`` and ``photocurrent_sesolve`` functions, are now included in the output of :func:`qutip.solver.mcsolve` as ``result.photocurrent``.
 
 
 .. plot::
@@ -236,6 +241,34 @@ For example, the following code block plots expectation values for 1, 10 and 100
     plt.legend()
     plt.show()
 
+
+.. openmcsolve:
+
+Open Systems
+------------
+
+``mcsolve`` can be used to study system with have measured and dissipative interaction with the bath.
+This is done by using a liouvillian including the dissipative interaction instead of an Hamiltonian.
+
+.. 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)
+    L = liouvillian(H, [0.01 * sm, np.sqrt(0.1) * a])
+    data = mcsolve(L, 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()
+
+
 .. plot::
     :context: reset
     :include-source: false