Set default trajectories to 500 and rename the keyword argument ensemble_method to ensemblealg

CHANGELOG.md

@@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Fix time evolution output when using `saveat` keyword argument. ([#398])
 - Align some attributes of `mcsolve`, `ssesolve` and `smesolve` results with `QuTiP`. ([#402])
 - Improve ensemble generation of `ssesolve` and change parameters handling on stochastic processes. ([#403])
+- Set default trajectories to 500 and rename the keyword argument `ensemble_method` to `ensemblealg`. ([#405])
 
 ## [v0.26.0]
 Release date: 2025-02-09
@@ -132,3 +133,4 @@ Release date: 2024-11-13
 [#398]: https://github.com/qutip/QuantumToolbox.jl/issues/398
 [#402]: https://github.com/qutip/QuantumToolbox.jl/issues/402
 [#403]: https://github.com/qutip/QuantumToolbox.jl/issues/403
+[#405]: https://github.com/qutip/QuantumToolbox.jl/issues/405

benchmarks/timeevolution.jl

@@ -50,7 +50,7 @@ function benchmark_timeevolution!(SUITE)
         ntraj = 100,
         e_ops = $e_ops,
         progress_bar = Val(false),
-        ensemble_method = EnsembleSerial(),
+        ensemblealg = EnsembleSerial(),
     )
     SUITE["Time Evolution"]["time-independent"]["mcsolve"]["Multithreaded"] = @benchmarkable mcsolve(
         $H,
@@ -60,7 +60,7 @@ function benchmark_timeevolution!(SUITE)
         ntraj = 100,
         e_ops = $e_ops,
         progress_bar = Val(false),
-        ensemble_method = EnsembleThreads(),
+        ensemblealg = EnsembleThreads(),
     )
 
     return nothing

docs/make.jl

@@ -34,7 +34,7 @@ const PAGES = [
     "Home" => "index.md",
     "Getting Started" => [
         "Brief Example" => "getting_started/brief_example.md",
-        "Key differences from QuTiP" => "getting_started/qutip_differences.md",
+        # "Key differences from QuTiP" => "getting_started/qutip_differences.md",
         "The Importance of Type-Stability" => "getting_started/type_stability.md",
         "Example: Create QuantumToolbox.jl Logo" => "getting_started/logo.md",
         # "Cite QuantumToolbox.jl" => "getting_started/cite.md",
@@ -51,7 +51,7 @@ const PAGES = [
             "Time Evolution Solutions" => "users_guide/time_evolution/solution.md",
             "Schrödinger Equation Solver" => "users_guide/time_evolution/sesolve.md",
             "Lindblad Master Equation Solver" => "users_guide/time_evolution/mesolve.md",
-            "Monte-Carlo Solver" => "users_guide/time_evolution/mcsolve.md",
+            "Monte Carlo Solver" => "users_guide/time_evolution/mcsolve.md",
             "Stochastic Solver" => "users_guide/time_evolution/stochastic.md",
             "Solving Problems with Time-dependent Hamiltonians" => "users_guide/time_evolution/time_dependent.md",
         ],

docs/src/users_guide/cluster.md

@@ -112,7 +112,7 @@ e_ops = [Sx, Sy, Sz]
 
 tlist = range(0, 10, 100)
 
-sol_mc = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=5000, ensemble_method=EnsembleSplitThreads())
+sol_mc = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=5000, ensemblealg=EnsembleSplitThreads())
 
 ##
 
@@ -131,7 +131,7 @@ println("Hello! You have $(nworkers()) workers with $(remotecall_fetch(Threads.n
 
 command, we test that the distributed network is correctly initialized. The `remotecall_fetch(Threads.nthreads, 2)` command returns the number of threads of the worker with ID `2`.
 
-We then write the main part of the script, where we define the lattice through the [`Lattice`](@ref) function. We set the parameters and define the Hamiltonian and collapse operators with the [`DissipativeIsing`](@ref) function. We also define the expectation operators `e_ops` and the initial state `ψ0`. Finally, we perform the Monte Carlo quantum trajectories with the [`mcsolve`](@ref) function. The `ensemble_method=EnsembleSplitThreads()` argument is used to parallelize the Monte Carlo quantum trajectories, by splitting the ensemble of trajectories among the workers. For a more detailed explanation of the different ensemble methods, you can check the [official documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/) of the [**DifferentialEquations.jl**](https://github.com/SciML/DifferentialEquations.jl/) package. Finally, the `rmprocs(workers())` command is used to remove the workers after the computation is finished.
+We then write the main part of the script, where we define the lattice through the [`Lattice`](@ref) function. We set the parameters and define the Hamiltonian and collapse operators with the [`DissipativeIsing`](@ref) function. We also define the expectation operators `e_ops` and the initial state `ψ0`. Finally, we perform the Monte Carlo quantum trajectories with the [`mcsolve`](@ref) function. The `ensemblealg=EnsembleSplitThreads()` argument is used to parallelize the Monte Carlo quantum trajectories, by splitting the ensemble of trajectories among the workers. For a more detailed explanation of the different ensemble methods, you can check the [official documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/) of the [**DifferentialEquations.jl**](https://github.com/SciML/DifferentialEquations.jl/) package. Finally, the `rmprocs(workers())` command is used to remove the workers after the computation is finished.
 
 The output of the script will be printed in the `output.out` file, which contains an output similar to the following:
 

docs/src/users_guide/time_evolution/intro.md

@@ -17,6 +17,9 @@
     - [Example: Harmonic oscillator in thermal bath](@ref doc-TE:Example:Harmonic-oscillator-in-thermal-bath)
     - [Example: Two-level atom coupled to dissipative single-mode cavity](@ref doc-TE:Example:Two-level-atom-coupled-to-dissipative-single-mode-cavity)
 - [Monte Carlo Solver](@ref doc-TE:Monte-Carlo-Solver)
+    - [Monte Carlo wave-function](@ref doc-TE:Monte-Carlo-wave-function)
+    - [Example: Two-level atom coupled to dissipative single-mode cavity (MC)](@ref doc-TE:Example:Two-level-atom-coupled-to-dissipative-single-mode-cavity-(MC))
+    - [Running trajectories in parallel](@ref doc-TE:Running-trajectories-in-parallel)
 - [Stochastic Solver](@ref doc-TE:Stochastic-Solver)
     - [Stochastic Schrödinger equation](@ref doc-TE:Stochastic-Schrödinger-equation)
     - [Stochastic master equation](@ref doc-TE:Stochastic-master-equation)

docs/src/users_guide/time_evolution/mcsolve.md

@@ -1,3 +1,134 @@
 # [Monte Carlo Solver](@id doc-TE:Monte-Carlo-Solver)
 
-This page is still under construction, please visit [API](@ref doc-API) first.
+## [Monte Carlo wave-function](@id doc-TE:Monte-Carlo-wave-function)
+
+Where as the density matrix formalism describes the ensemble average over many identical realizations of a quantum system, the Monte Carlo (MC), or quantum-jump approach to wave function evolution, allows for simulating an individual realization of the system dynamics. Here, the environment is continuously monitored, resulting in a series of quantum jumps in the system wave function, conditioned on the increase in information gained about the state of the system via the environmental measurements. In general, this evolution is governed by the Schrödinger equation with a non-Hermitian effective Hamiltonian
+
+```math
+\hat{H}_{\textrm{eff}} = \hat{H} - \frac{i}{2} \sum_{n=1}^N \hat{C}_n^\dagger \hat{C}_n.
+```
+
+where ``\hat{H}`` is the system Hamiltonian and ``\hat{C}_n`` are collapse (jump) operators (assume ``N`` is the total number of collapse operators). Each collapse operator corresponds to a separate irreversible process with rate ``\gamma_n``. Here, the strictly negative non-Hermitian portion of the above equation gives rise to a reduction in the norm of the wave function, that to first-order in a small time ``\delta t``, is given by
+
+```math
+\langle \psi(t + \delta t) | \psi(t + \delta t) \rangle = 1 - \delta p,
+```
+
+where
+
+```math
+\delta p = \delta t \sum_{n=1}^N \langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle,
+```
+
+and ``\delta t`` is such that ``\delta p \ll 1``. With a probability of remaining in the state ``| \psi(t + \delta t) \rangle`` given by ``1 - \delta p``, the corresponding quantum jump probability is thus ``\delta p``. If the environmental measurements register a quantum jump, say via the emission of a photon into the environment, or a change in the spin of a quantum dot, the wave function undergoes a jump into a state defined by projecting ``| \psi(t) \rangle`` using the collapse operator ``\hat{C}_n`` corresponding to the measurement
+
+```math
+| \psi(t+\delta t) \rangle = \frac{\hat{C}_n |\psi(t)\rangle}{ \sqrt{\langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle} }.
+```
+
+If more than a single collapse operator is present in ``\hat{H}_{\textrm{eff}}``, the probability of collapse due to the ``n``-th operator ``\hat{C}_n`` is given by
+
+```math
+P_n(t) = \frac{1}{\delta p}\langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle.
+```
+
+Note that the probability of all collapses should be normalized to unity for all time ``t``, namely
+
+```math
+\sum_{n=1}^N P_n(t) = 1 ~~~\forall~~t.
+```
+
+Evaluating the MC evolution to first-order in time is quite tedious. Instead, `QuantumToolbox.jl` provides the function [`mcsolve`](@ref) which uses the following algorithm to simulate a single realization of a quantum system. Starting from a pure state ``| \psi(0) \rangle``:
+
+1. Choose two random numbers (``r_1`` and ``r_2``) between 0 and 1, where ``r_1`` represents the probability that a quantum jump occurs and  ``r_2`` is used to select which collapse operator was responsible for the jump.
+1. Integrate the Schrödinger equation with respect to the effective Hamiltonian ``\hat{H}_{\textrm{eff}}`` until a time ``\tau`` such that the norm of the wave function satisfies ``\langle \psi(\tau) | \psi(\tau) \rangle = r_1``, at which point a jump occurs
+1. The resultant jump projects the system at time ``\tau`` into one of the renormalized states ``| \psi(\tau + \delta t) \rangle``. The corresponding collapse operator ``\hat{C}_n`` is chosen such that ``\tilde{n} \leq N`` is the smallest integer satisfying ``\sum_{n=1}^{\tilde{n}} P_n(\tau) \geq r_2``.
+1. Using the renormalized state from previous step as the new initial condition at time ``\tau`` and repeat the above procedure until the final simulation time is reached.
+
+## [Example: Two-level atom coupled to dissipative single-mode cavity (MC)](@id doc-TE:Example:Two-level-atom-coupled-to-dissipative-single-mode-cavity-(MC))
+
+In `QuantumToolbox.jl`, Monte Carlo evolution is implemented with the [`mcsolve`](@ref) function. It takes nearly the same arguments as the [`mesolve`](@ref) function for [Lindblad master equation evolution](@ref doc-TE:Lindblad-Master-Equation-Solver), except that the initial state must be a [`Ket`](@ref) vector, as oppose to a density matrix, and there is an optional keyword argument `ntraj` that defines the number of stochastic trajectories to be simulated. By default, `ntraj=500` indicating that `500` Monte Carlo trajectories will be performed.
+
+To illustrate the use of the Monte Carlo evolution of quantum systems in `QuantumToolbox.jl`, let’s again consider the case of a two-level atom coupled to a leaky cavity. The only differences to the master equation treatment is that in this case we invoke the [`mcsolve`](@ref) function instead of [`mesolve`](@ref)
+
+```@setup mcsolve
+using QuantumToolbox
+
+using CairoMakie
+CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
+```
+
+```@example mcsolve
+times = LinRange(0.0, 10.0, 200)
+
+ψ0 = tensor(fock(2, 0), fock(10, 8))
+a  = tensor(qeye(2), destroy(10))
+σm = tensor(destroy(2), qeye(10))
+H = 2 * π * a' * a + 2 * π * σm' * σm + 2 * π * 0.25 * (σm * a' + σm' * a)
+
+c_ops = [sqrt(0.1) * a]
+e_ops = [a' * a, σm' * σm]
+
+sol_500 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops)
+
+# plot by CairoMakie.jl
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1],
+    xlabel = "Time",
+    ylabel = "Expectation values",
+    title = "Monte Carlo time evolution (500 trajectories)",
+)
+lines!(ax, times, real(sol_500.expect[1,:]), label = "cavity photon number", linestyle = :solid)
+lines!(ax, times, real(sol_500.expect[2,:]), label = "atom excitation probability", linestyle = :dash)
+
+axislegend(ax, position = :rt)
+
+fig
+```
+
+The advantage of the Monte Carlo method over the master equation approach is that only the state vector ([`Ket`](@ref)) 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. 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.
+
+We can also change the number of trajectories (`ntraj`). This can be used to explore the convergence of the Monte Carlo solver. For example, the following code plots the expectation values for `1`, `10` and `100` trajectories:
+
+```@example mcsolve
+e_ops = [a' * a]
+
+sol_1   = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 1)
+sol_10  = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 10)
+sol_100 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 100)
+
+# plot by CairoMakie.jl
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1],
+    xlabel = "Time",
+    ylabel = "Expectation values",
+    title = "Monte Carlo time evolution",
+)
+lines!(ax, times, real(sol_1.expect[1,:]), label = "1 trajectory", linestyle = :dashdot)
+lines!(ax, times, real(sol_10.expect[1,:]), label = "10 trajectories", linestyle = :dash)
+lines!(ax, times, real(sol_100.expect[1,:]), label = "100 trajectories", linestyle = :solid)
+
+axislegend(ax, position = :rt)
+
+fig
+```
+
+## [Running trajectories in parallel](@id doc-TE:Running-trajectories-in-parallel)
+
+Monte Carlo evolutions often need hundreds of trajectories to obtain sufficient statistics. Since all trajectories are independent of each other, they can be computed in parallel. The keyword argument `ensemblealg` can specify how the multiple trajectories are handled. The common ensemble methods are:
+
+- `EnsembleSerial()`: No parallelism
+- `EnsembleThreads()`: **The default.** This uses multithreading.
+- `EnsembleDistributed()`: This uses as many processors as you have Julia processes.
+- `EnsembleSplitThreads()`: This uses multithreading on each process.
+
+!!! note "Other Ensemble Algorithms"
+    See the [documentation of `DifferentialEquations.jl`](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/) for more details. Also, see Julia's documentation for more details about multithreading and adding more processes.
+
+```julia
+sol_serial   = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ensemblealg=EnsembleSerial())
+sol_parallel = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ensemblealg=EnsembleThreads());
+```
+
+!!! tip "Parallelization on a Cluster"
+    See the section [Intensive parallelization on a Cluster](@ref doc:Intensive-parallelization-on-a-Cluster) for more details.

src/QuantumToolbox.jl

@@ -28,6 +28,7 @@ import SciMLBase:
     ODEProblem,
     SDEProblem,
     EnsembleProblem,
+    EnsembleAlgorithm,
     EnsembleSerial,
     EnsembleThreads,
     EnsembleSplitThreads,