Store both times and times_states in time evolution solutions (#506)

Author: ytdHuang

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Implement `EnrSpace` and corresponding functionality. ([#500])
 - Check for orthogonality breakdown in `Lanczos` solver for `spectrum`. ([#501])
+- Store both `times` and `times_states` in time evolution solutions. ([#506], [#504])
 - Fix errors in `Julia v1.12`. ([#507])
 
 ## [v0.32.1]
@@ -259,4 +260,6 @@ Release date: 2024-11-13
 [#494]: https://github.com/qutip/QuantumToolbox.jl/issues/494
 [#500]: https://github.com/qutip/QuantumToolbox.jl/issues/500
 [#501]: https://github.com/qutip/QuantumToolbox.jl/issues/501
+[#504]: https://github.com/qutip/QuantumToolbox.jl/issues/504
+[#506]: https://github.com/qutip/QuantumToolbox.jl/issues/506
 [#507]: https://github.com/qutip/QuantumToolbox.jl/issues/507

docs/src/users_guide/time_evolution/mcsolve.md

@@ -59,7 +59,7 @@ CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
 ```
 
 ```@example mcsolve
-times = LinRange(0.0, 10.0, 200)
+tlist = LinRange(0.0, 10.0, 200)
 
 ψ0 = tensor(fock(2, 0), fock(10, 8))
 a  = tensor(qeye(2), destroy(10))
@@ -69,7 +69,7 @@ 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)
+sol_500 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops)
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
@@ -78,8 +78,8 @@ ax = Axis(fig[1, 1],
     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)
+lines!(ax, tlist, real(sol_500.expect[1,:]), label = "cavity photon number", linestyle = :solid)
+lines!(ax, tlist, real(sol_500.expect[2,:]), label = "atom excitation probability", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 
@@ -93,9 +93,9 @@ We can also change the number of trajectories (`ntraj`). This can be used to exp
 ```@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)
+sol_1   = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 1)
+sol_10  = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 10)
+sol_100 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 100)
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
@@ -104,9 +104,9 @@ ax = Axis(fig[1, 1],
     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)
+lines!(ax, tlist, real(sol_1.expect[1,:]), label = "1 trajectory", linestyle = :dashdot)
+lines!(ax, tlist, real(sol_10.expect[1,:]), label = "10 trajectories", linestyle = :dash)
+lines!(ax, tlist, real(sol_100.expect[1,:]), label = "100 trajectories", linestyle = :solid)
 
 axislegend(ax, position = :rt)
 
@@ -126,8 +126,8 @@ Monte Carlo evolutions often need hundreds of trajectories to obtain sufficient
     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());
+sol_serial   = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ensemblealg=EnsembleSerial())
+sol_parallel = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ensemblealg=EnsembleThreads());
 ```
 
 !!! tip "Parallelization on a Cluster"

docs/src/users_guide/time_evolution/mesolve.md

@@ -128,15 +128,14 @@ sol = mesolve(H, ψ0, tlist, c_ops, e_ops = [sigmaz(), sigmay()])
 We can therefore plot the expectation values:
 
 ```@example mesolve
-times = sol.times
 expt_z = real(sol.expect[1,:])
 expt_y = real(sol.expect[2,:])
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
-lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
-lines!(ax, times, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
+lines!(ax, tlist, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
+lines!(ax, tlist, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 
@@ -210,17 +209,15 @@ tlist = LinRange(0.0, 10.0, 200)
 L = liouvillian(H_a + H_c + H_I, c_ops)
 sol = mesolve(L, ψ0, tlist, e_ops=[σm' * σm, a' * a])
 
-times = sol.times
-
 # expectation value of Number operator
 N_atom   = real(sol.expect[1,:])
 N_cavity = real(sol.expect[2,:])
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
-lines!(ax, times, N_atom, label = "atom excitation probability", linestyle = :solid)
-lines!(ax, times, N_cavity, label = "cavity photon number", linestyle = :dash)
+lines!(ax, tlist, N_atom, label = "atom excitation probability", linestyle = :solid)
+lines!(ax, tlist, N_cavity, label = "cavity photon number", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 

docs/src/users_guide/time_evolution/sesolve.md

@@ -59,8 +59,8 @@ sol = sesolve(H, ψ0, tlist, e_ops = [sigmaz(), sigmay()])
     Here, we call [`sesolve`](@ref) directly instead of pre-defining [`sesolveProblem`](@ref) first (as shown previously).
 
 ```@example sesolve
-times = sol.times
-print(size(times))
+println(size(sol.times)) # time points corresponds to stored expectation values 
+println(size(sol.times_states)) # time points corresponds to stored states
 ```
 
 ```@example sesolve
@@ -77,8 +77,8 @@ expt_y = real(expt[2,:])
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
-lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
-lines!(ax, times, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
+lines!(ax, tlist, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
+lines!(ax, tlist, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
 
 axislegend(ax, position = :rb)
 
@@ -91,6 +91,9 @@ If the keyword argument `e_ops` is not specified (or given as an empty `Vector`)
 tlist = [0, 10]
 sol = sesolve(H, ψ0, tlist) # or specify: e_ops = []
 
+println(size(sol.times))
+println(size(sol.times_states))
+
 sol.states
 ```
 
@@ -104,9 +107,11 @@ sol = sesolve(H, ψ0, tlist, e_ops = [sigmay()], saveat = tlist)
 ```
 
 ```@example sesolve
+println(size(sol.times))
 sol.expect
 ```
 
 ```@example sesolve
+println(size(sol.times_states))
 sol.states
 ```

docs/src/users_guide/time_evolution/solution.md

@@ -12,8 +12,9 @@ CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
 
 | **Fields (Attributes)** | **Description** |
 |:------------------------|:----------------|
-| `sol.times` | The time list of the evolution. |
-| `sol.states` | The list of result states. |
+| `sol.times` | The list of time points at which the expectation values are calculated during the evolution. |
+| `sol.times_states` | The list of time points at which the states are stored during the evolution. |
+| `sol.states` | The list of result states corresponding to each time point in `sol.times_states`. |
 | `sol.expect` | The expectation values corresponding to each time point in `sol.times`. |
 | `sol.alg` | The algorithm which is used during the solving process. |
 | `sol.abstol` | The absolute tolerance which is used during the solving process. |
@@ -54,7 +55,7 @@ nothing # hide
 
 Recall that `Julia` uses `Fortran`-style indexing that begins with one (i.e., `[1,:]` represents the 1-st observable, where `:` represents all values corresponding to `tlist`).
 
-Together with the array of times at which these expectation values are calculated:
+Together with the list of time points at which these expectation values are calculated:
 
 ```@example TE-solution
 times = sol.times
@@ -83,6 +84,13 @@ State vectors, or density matrices, are accessed in a similar manner:
 sol.states
 ```
 
+Together with the list of time points at which these states are stored:
+
+```@example TE-solution
+times = sol.times_states
+nothing # hide
+```
+
 Here, the solution contains only one (final) state. Because the `states` will be saved depend on the keyword argument `saveat` in `kwargs`. If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). One can also specify `e_ops` and `saveat` separately.
 
 Some other solvers can have other output.

src/correlations.jl

@@ -19,7 +19,7 @@ end
         C::QuantumObject;
         kwargs...)
 
-Returns the two-times correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
+Returns the two-time correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
 
 If the initial state `ψ0` is given as `nothing`, then the [`steadystate`](@ref) will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
 """
@@ -96,7 +96,7 @@ end
         reverse::Bool=false,
         kwargs...)
 
-Returns the two-times correlation function of two operators ``\hat{A}`` and ``\hat{B}`` : ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
+Returns the two-time correlation function of two operators ``\hat{A}`` and ``\hat{B}`` : ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
 
 If the initial state `ψ0` is given as `nothing`, then the [`steadystate`](@ref) will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
 

src/spectrum.jl

@@ -313,7 +313,7 @@ end
 Calculate the power spectrum corresponding to a two-time correlation function using fast Fourier transform (FFT).
 
 # Parameters
-- `tlist::AbstractVector`: List of times at which the two-time correlation function is given.
+- `tlist::AbstractVector`: List of time points at which the two-time correlation function is given.
 - `corr::AbstractVector`: List of two-time correlations corresponding to the given time point in `tlist`.
 - `inverse::Bool`: Whether to use the inverse Fourier transform or not. Default to `false`.
 

src/time_evolution/brmesolve.jl

@@ -158,7 +158,7 @@ Solves for the dynamics of a system using the Bloch-Redfield master equation, gi
 
 - `H`: The system Hamiltonian. Must be an [`Operator`](@ref)
 - `ψ0`:  Initial state of the system $|\psi(0)\rangle$. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `a_ops`: Nested list with each element is a `Tuple` of operator-function pairs `(a_op, spectra)`, and the coupling [`Operator`](@ref) `a_op` must be hermitian with corresponding `spectra` being a `Function` of transition energy
 - `c_ops`: List of collapse operators corresponding to Lindblad dissipators
 - `sec_cutoff`: Cutoff for secular approximation. Use `-1` if secular approximation is not used when evaluating bath-coupling terms.

src/time_evolution/lr_mesolve.jl

@@ -7,10 +7,11 @@ A structure storing the results and some information from solving low-rank maste
 
 # Fields (Attributes)
 
-- `times::AbstractVector`: The time list of the evolution.
-- `states::Vector{QuantumObject}`: The list of result states.
+- `times::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
+- `times_states::AbstractVector`: The list of time points at which the states are stored during the evolution.
+- `states::Vector{QuantumObject}`: The list of result states corresponding to each time point in `times_states`.
 - `expect::Matrix`: The expectation values corresponding to each time point in `times`.
-- `fexpect::Matrix`: The function values at each time point.
+- `fexpect::Matrix`: The function values corresponding to each time point in `times`.
 - `retcode`: The return code from the solver.
 - `alg`: The algorithm which is used during the solving process.
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
@@ -19,7 +20,8 @@ A structure storing the results and some information from solving low-rank maste
 - `B::Vector{QuantumObject}`: The `B` matrix of the low-rank algorithm at each time point.
 """
 struct TimeEvolutionLRSol{
-    TT<:AbstractVector{<:Real},
+    TT1<:AbstractVector{<:Real},
+    TT2<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Matrix{ComplexF64},
     RetT<:Enum,
@@ -29,7 +31,8 @@ struct TimeEvolutionLRSol{
     TSZB<:AbstractVector,
     TM<:Vector{<:Integer},
 }
-    times::TT
+    times::TT1
+    times_states::TT2
     states::TS
     expect::TE
     fexpect::TE
@@ -549,6 +552,7 @@ function lr_mesolve(prob::ODEProblem; kwargs...)
     ρt = map(x -> Qobj(x[1] * x[2] * x[1]', type = Operator(), dims = prob.p.Hdims), zip(zt, Bt))
 
     return TimeEvolutionLRSol(
+        prob.p.times,
         sol.t,
         ρt,
         prob.p.expvals,

src/time_evolution/mcsolve.jl

@@ -89,7 +89,7 @@ If the environmental measurements register a quantum jump, the wave function und
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -195,7 +195,7 @@ If the environmental measurements register a quantum jump, the wave function und
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -320,7 +320,7 @@ If the environmental measurements register a quantum jump, the wave function und
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `alg`: The algorithm to use for the ODE solver. Default to `Tsit5()`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -412,6 +412,7 @@ function mcsolve(
     return TimeEvolutionMCSol(
         ntraj,
         ens_prob_mc.times,
+        _sol_1.t,
         states,
         expvals,
         expvals, # This is average_expect