add field rng for multi-trajectory solutions

Author: ytdHuang

src/time_evolution/mcsolve.jl

@@ -256,7 +256,7 @@ function mcsolveEnsembleProblem(
         EnsembleProblem(prob_mc.prob, prob_func = _prob_func, output_func = _output_func[1], safetycopy = false),
         prob_mc.times,
         prob_mc.dimensions,
-        (progr = _output_func[2], channel = _output_func[3]),
+        (progr = _output_func[2], channel = _output_func[3], rng = rng),
     )
 
     return ensemble_prob
@@ -412,11 +412,12 @@ function mcsolve(
         ens_prob_mc.times,
         states,
         expvals_all,
+        ens_prob_mc.kwargs.rng,
         col_times,
         col_which,
         sol.converged,
         _sol_1.alg,
-        NamedTuple(_sol_1.prob.kwargs).abstol,
-        NamedTuple(_sol_1.prob.kwargs).reltol,
+        _sol_1.prob.kwargs[:abstol],
+        _sol_1.prob.kwargs[:reltol],
     )
 end

src/time_evolution/smesolve.jl

@@ -271,7 +271,7 @@ function smesolveEnsembleProblem(
         EnsembleProblem(prob_sme, prob_func = _prob_func, output_func = _output_func[1], safetycopy = true),
         prob_sme.times,
         prob_sme.dimensions,
-        merge(prob_sme.kwargs, (progr = _output_func[2], channel = _output_func[3])),
+        merge(prob_sme.kwargs, (progr = _output_func[2], channel = _output_func[3], rng = rng)),
     )
 
     return ensemble_prob
@@ -423,6 +423,7 @@ function smesolve(
         ens_prob.times,
         states,
         expvals_all,
+        ens_prob.kwargs.rng,
         m_expvals, # Measurement expectation values
         sol.converged,
         _sol_1.alg,

src/time_evolution/ssesolve.jl

@@ -264,7 +264,7 @@ function ssesolveEnsembleProblem(
         EnsembleProblem(prob_sme, prob_func = _prob_func, output_func = _output_func[1], safetycopy = true),
         prob_sme.times,
         prob_sme.dimensions,
-        (progr = _output_func[2], channel = _output_func[3]),
+        (progr = _output_func[2], channel = _output_func[3], rng = rng),
     )
 
     return ensemble_prob
@@ -417,6 +417,7 @@ function ssesolve(
         ens_prob.times,
         states,
         expvals_all,
+        ens_prob.kwargs.rng,
         m_expvals, # Measurement expectation values
         sol.converged,
         _sol_1.alg,

src/time_evolution/time_evolution.jl

@@ -103,6 +103,7 @@ A structure storing the results and some information from solving quantum trajec
 - `times::AbstractVector`: The time list of the evolution.
 - `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
 - `expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`.
+- `rng::AbstractRNG`: Random number generator for reproducibility.
 - `col_times::Vector{Vector{Real}}`: The time records of every quantum jump occurred in each trajectory.
 - `col_which::Vector{Vector{Int}}`: The indices of which collapse operator was responsible for each quantum jump in `col_times`.
 - `converged::Bool`: Whether the solution is converged or not.
@@ -127,6 +128,7 @@ struct TimeEvolutionMCSol{
     TT<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Union{AbstractArray,Nothing},
+    TR<:AbstractRNG,
     TJT<:Vector{<:Vector{<:Real}},
     TJW<:Vector{<:Vector{<:Integer}},
     AlgT<:OrdinaryDiffEqAlgorithm,
@@ -137,6 +139,7 @@ struct TimeEvolutionMCSol{
     times::TT
     states::TS
     expect::TE
+    rng::TR
     col_times::TJT
     col_which::TJW
     converged::Bool
@@ -173,6 +176,8 @@ A structure storing the results and some information from solving trajectories o
 - `times::AbstractVector`: The time list of the evolution.
 - `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times`.
 - `expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`.
+- `rng::AbstractRNG`: Random number generator for reproducibility.
+- `measurement::Union{AbstractArray,Nothing}`: Measurements for each trajectories and stochastic collapse operators (`sc_ops`).
 - `converged::Bool`: Whether the solution is converged or not.
 - `alg`: The algorithm which is used during the solving process.
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
@@ -195,6 +200,7 @@ struct TimeEvolutionStochasticSol{
     TT<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Union{AbstractArray,Nothing},
+    TR<:AbstractRNG,
     TEM<:Union{AbstractArray,Nothing},
     AlgT<:StochasticDiffEqAlgorithm,
     AT<:Real,
@@ -204,6 +210,7 @@ struct TimeEvolutionStochasticSol{
     times::TT
     states::TS
     expect::TE
+    rng::TR
     measurement::TEM
     converged::Bool
     alg::AlgT

test/core-test/time_evolution.jl

@@ -751,11 +751,9 @@
                 rng = rng,
             )
 
-            rng = MersenneTwister(1234)
-            sol_mc2 = mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = rng)
-            rng = MersenneTwister(1234)
-            sol_sse2 = ssesolve(H, ψ0, tlist, c_ops, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = rng)
-            rng = MersenneTwister(1234)
+            sol_mc2 = mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = sol_mc1.rng)
+            sol_sse2 =
+                ssesolve(H, ψ0, tlist, c_ops, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = sol_sse1.rng)
             sol_sme2 = smesolve(
                 H,
                 ψ0,
@@ -765,14 +763,13 @@
                 ntraj = 50,
                 e_ops = e_ops,
                 progress_bar = Val(false),
-                rng = rng,
+                rng = sol_sme1.rng,
             )
 
-            rng = MersenneTwister(1234)
-            sol_mc3 = mcsolve(H, ψ0, tlist, c_ops, ntraj = 510, e_ops = e_ops, progress_bar = Val(false), rng = rng)
-            rng = MersenneTwister(1234)
-            sol_sse3 = ssesolve(H, ψ0, tlist, c_ops, ntraj = 60, e_ops = e_ops, progress_bar = Val(false), rng = rng)
-            rng = MersenneTwister(1234)
+            sol_mc3 =
+                mcsolve(H, ψ0, tlist, c_ops, ntraj = 510, e_ops = e_ops, progress_bar = Val(false), rng = sol_mc1.rng)
+            sol_sse3 =
+                ssesolve(H, ψ0, tlist, c_ops, ntraj = 60, e_ops = e_ops, progress_bar = Val(false), rng = sol_sse1.rng)
             sol_sme3 = smesolve(
                 H,
                 ψ0,
@@ -782,7 +779,7 @@
                 ntraj = 60,
                 e_ops = e_ops,
                 progress_bar = Val(false),
-                rng = rng,
+                rng = sol_sme1.rng,
             )
 
             @test sol_mc1.expect ≈ sol_mc2.expect atol = 1e-10