From 8c41b87b075400505545d064cce3ea1076622901 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Sun, 16 Feb 2025 20:13:07 +0100 Subject: [PATCH 1/5] [no ci] Support for single `sc_ops` for faster specific method --- src/QuantumToolbox.jl | 4 +-- src/time_evolution/smesolve.jl | 50 +++++++++++++++------------- src/time_evolution/ssesolve.jl | 50 +++++++++++++++------------- src/time_evolution/time_evolution.jl | 36 +++++++++++++++----- test/core-test/time_evolution.jl | 23 +++++++++++-- 5 files changed, 103 insertions(+), 60 deletions(-) diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl index f1bb81b95..b88f8c519 100644 --- a/src/QuantumToolbox.jl +++ b/src/QuantumToolbox.jl @@ -40,7 +40,7 @@ import SciMLBase: AbstractSciMLProblem, AbstractODEIntegrator, AbstractODESolution -import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA1 +import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA1, SRIW1 import SciMLOperators: SciMLOperators, AbstractSciMLOperator, @@ -56,7 +56,7 @@ import DiffEqBase: get_tstops import DiffEqCallbacks: PeriodicCallback, PresetTimeCallback, TerminateSteadyState import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm import OrdinaryDiffEqTsit5: Tsit5 -import DiffEqNoiseProcess: RealWienerProcess! +import DiffEqNoiseProcess: RealWienerProcess!, RealWienerProcess # other dependencies (in alphabetical order) import ArrayInterface: allowed_getindex, allowed_setindex! diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl index 32f02be70..f5f52d024 100644 --- a/src/time_evolution/smesolve.jl +++ b/src/time_evolution/smesolve.jl @@ -15,7 +15,7 @@ _smesolve_ScalarOperator(op_vec) = ψ0::QuantumObject, tlist::AbstractVector, c_ops::Union{Nothing,AbstractVector,Tuple} = nothing, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -50,7 +50,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref) or a [`Operator`](@ref). - `tlist`: List of times at which to save either the state or the expectation values of the system. - `c_ops`: List of collapse operators ``\{\hat{C}_i\}_i``. It can be either a `Vector` or a `Tuple`. -- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector` or a `Tuple`. +- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided. - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`. - `params`: `NullParameters` of parameters to pass to the solver. - `rng`: Random number generator for reproducibility. @@ -74,7 +74,7 @@ function smesolveProblem( ψ0::QuantumObject{StateOpType}, tlist::AbstractVector, c_ops::Union{Nothing,AbstractVector,Tuple} = nothing, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -87,10 +87,12 @@ function smesolveProblem( isnothing(sc_ops) && throw(ArgumentError("The list of stochastic collapse operators must be provided. Use mesolveProblem instead.")) + sc_ops_list = _make_c_ops_list(sc_ops) # If it is an AbstractQuantumObject but we need to iterate + sc_ops_isa_Qobj = sc_ops isa AbstractQuantumObject # We can avoid using non-diagonal noise if sc_ops is just an AbstractQuantumObject tlist = _check_tlist(tlist, _FType(ψ0)) - L_evo = _mesolve_make_L_QobjEvo(H, c_ops) + _mesolve_make_L_QobjEvo(nothing, sc_ops) + L_evo = _mesolve_make_L_QobjEvo(H, c_ops) + _mesolve_make_L_QobjEvo(nothing, sc_ops_list) check_dimensions(L_evo, ψ0) dims = L_evo.dimensions @@ -99,7 +101,7 @@ function smesolveProblem( progr = ProgressBar(length(tlist), enable = getVal(progress_bar)) - sc_ops_evo_data = Tuple(map(get_data ∘ QobjEvo, sc_ops)) + sc_ops_evo_data = Tuple(map(get_data ∘ QobjEvo, sc_ops_list)) K = get_data(L_evo) @@ -116,7 +118,7 @@ function smesolveProblem( kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_SOLVER_OPTIONS; kwargs...) kwargs3 = _generate_stochastic_kwargs( e_ops, - sc_ops, + sc_ops_list, makeVal(progress_bar), tlist, makeVal(store_measurement), @@ -125,14 +127,8 @@ function smesolveProblem( ) tspan = (tlist[1], tlist[end]) - noise = RealWienerProcess!( - tlist[1], - zeros(length(sc_ops)), - zeros(length(sc_ops)), - save_everystep = getVal(store_measurement), - rng = rng, - ) - noise_rate_prototype = similar(ρ0, length(ρ0), length(sc_ops)) + noise = _make_noise(tspan[1], sc_ops, store_measurement, rng) + noise_rate_prototype = sc_ops_isa_Qobj ? nothing : similar(ρ0, length(ρ0), length(sc_ops_list)) prob = SDEProblem{true}( K, D, @@ -153,7 +149,7 @@ end ψ0::QuantumObject, tlist::AbstractVector, c_ops::Union{Nothing,AbstractVector,Tuple} = nothing, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -192,7 +188,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref) or a [`Operator`](@ref). - `tlist`: List of times at which to save either the state or the expectation values of the system. - `c_ops`: List of collapse operators ``\{\hat{C}_i\}_i``. It can be either a `Vector` or a `Tuple`. -- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector` or a `Tuple`. +- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided. - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`. - `params`: `NullParameters` of parameters to pass to the solver. - `rng`: Random number generator for reproducibility. @@ -220,7 +216,7 @@ function smesolveEnsembleProblem( ψ0::QuantumObject{StateOpType}, tlist::AbstractVector, c_ops::Union{Nothing,AbstractVector,Tuple} = nothing, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -239,7 +235,7 @@ function smesolveEnsembleProblem( ntraj, tlist, _stochastic_prob_func; - n_sc_ops = length(sc_ops), + sc_ops = sc_ops, store_measurement = makeVal(store_measurement), ) : prob_func _output_func = @@ -276,8 +272,8 @@ end ψ0::QuantumObject, tlist::AbstractVector, c_ops::Union{Nothing,AbstractVector,Tuple} = nothing, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; - alg::StochasticDiffEqAlgorithm = SRA1(), + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; + alg::Union{Nothing,StochasticDiffEqAlgorithm} = nothing, e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -316,8 +312,8 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref) or a [`Operator`](@ref). - `tlist`: List of times at which to save either the state or the expectation values of the system. - `c_ops`: List of collapse operators ``\{\hat{C}_i\}_i``. It can be either a `Vector` or a `Tuple`. -- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector` or a `Tuple`. -- `alg`: The algorithm to use for the stochastic differential equation. Default is `SRA1()`. +- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided. +- `alg`: The algorithm to use for the stochastic differential equation. Default is `SRIW1()` if `sc_ops` is an [`AbstractQuantumObject`](@ref) (diagonal noise), and `SRA1()` otherwise (non-diagonal noise). - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`. - `params`: `NullParameters` of parameters to pass to the solver. - `rng`: Random number generator for reproducibility. @@ -345,8 +341,8 @@ function smesolve( ψ0::QuantumObject{StateOpType}, tlist::AbstractVector, c_ops::Union{Nothing,AbstractVector,Tuple} = nothing, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; - alg::StochasticDiffEqAlgorithm = SRA1(), + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; + alg::Union{Nothing,StochasticDiffEqAlgorithm} = nothing, e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -376,6 +372,12 @@ function smesolve( kwargs..., ) + sc_ops_isa_Qobj = sc_ops isa AbstractQuantumObject # We can avoid using non-diagonal noise if sc_ops is just an AbstractQuantumObject + + if isnothing(alg) + alg = sc_ops_isa_Qobj ? SRIW1() : SRA1() + end + return smesolve(ensemble_prob, alg, ntraj, ensemblealg) end diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl index c428b1547..264ff32f6 100644 --- a/src/time_evolution/ssesolve.jl +++ b/src/time_evolution/ssesolve.jl @@ -17,7 +17,7 @@ _ScalarOperator_e2_2(op, f = +) = H::Union{AbstractQuantumObject{OperatorQuantumObject},Tuple}, ψ0::QuantumObject{KetQuantumObject}, tlist::AbstractVector, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -52,7 +52,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th - `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. -- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector` or a `Tuple`. +- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided. - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`. - `params`: `NullParameters` of parameters to pass to the solver. - `rng`: Random number generator for reproducibility. @@ -75,7 +75,7 @@ function ssesolveProblem( H::Union{AbstractQuantumObject{OperatorQuantumObject},Tuple}, ψ0::QuantumObject{KetQuantumObject}, tlist::AbstractVector, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -88,10 +88,12 @@ function ssesolveProblem( sc_ops isa Nothing && throw(ArgumentError("The list of stochastic collapse operators must be provided. Use sesolveProblem instead.")) + sc_ops_list = _make_c_ops_list(sc_ops) # If it is an AbstractQuantumObject but we need to iterate + sc_ops_isa_Qobj = sc_ops isa AbstractQuantumObject # We can avoid using non-diagonal noise if sc_ops is just an AbstractQuantumObject tlist = _check_tlist(tlist, _FType(ψ0)) - H_eff_evo = _mcsolve_make_Heff_QobjEvo(H, sc_ops) + H_eff_evo = _mcsolve_make_Heff_QobjEvo(H, sc_ops_list) isoper(H_eff_evo) || throw(ArgumentError("The Hamiltonian must be an Operator.")) check_dimensions(H_eff_evo, ψ0) dims = H_eff_evo.dimensions @@ -100,7 +102,7 @@ function ssesolveProblem( progr = ProgressBar(length(tlist), enable = getVal(progress_bar)) - sc_ops_evo_data = Tuple(map(get_data ∘ QobjEvo, sc_ops)) + sc_ops_evo_data = Tuple(map(get_data ∘ QobjEvo, sc_ops_list)) # Here the coefficients depend on the state, so this is a non-linear operator, which should be implemented with FunctionOperator instead. However, the nonlinearity is only on the coefficients, and it should be safe. K_l = sum( @@ -116,7 +118,7 @@ function ssesolveProblem( kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_SOLVER_OPTIONS; kwargs...) kwargs3 = _generate_stochastic_kwargs( e_ops, - sc_ops, + sc_ops_list, makeVal(progress_bar), tlist, makeVal(store_measurement), @@ -125,14 +127,8 @@ function ssesolveProblem( ) tspan = (tlist[1], tlist[end]) - noise = RealWienerProcess!( - tlist[1], - zeros(length(sc_ops)), - zeros(length(sc_ops)), - save_everystep = getVal(store_measurement), - rng = rng, - ) - noise_rate_prototype = similar(ψ0, length(ψ0), length(sc_ops)) + noise = _make_noise(tspan[1], sc_ops, store_measurement, rng) + noise_rate_prototype = sc_ops_isa_Qobj ? nothing : similar(ψ0, length(ψ0), length(sc_ops_list)) prob = SDEProblem{true}( K, D, @@ -152,7 +148,7 @@ end H::Union{AbstractQuantumObject{OperatorQuantumObject},Tuple}, ψ0::QuantumObject{KetQuantumObject}, tlist::AbstractVector, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -191,7 +187,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is t - `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. -- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector` or a `Tuple`. +- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided. - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`. - `params`: `NullParameters` of parameters to pass to the solver. - `rng`: Random number generator for reproducibility. @@ -219,7 +215,7 @@ function ssesolveEnsembleProblem( H::Union{AbstractQuantumObject{OperatorQuantumObject},Tuple}, ψ0::QuantumObject{KetQuantumObject}, tlist::AbstractVector, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -238,7 +234,7 @@ function ssesolveEnsembleProblem( ntraj, tlist, _stochastic_prob_func; - n_sc_ops = length(sc_ops), + sc_ops = sc_ops, store_measurement = makeVal(store_measurement), ) : prob_func _output_func = @@ -273,8 +269,8 @@ end H::Union{AbstractQuantumObject{OperatorQuantumObject},Tuple}, ψ0::QuantumObject{KetQuantumObject}, tlist::AbstractVector, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; - alg::StochasticDiffEqAlgorithm = SRA1(), + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; + alg::Union{Nothing,StochasticDiffEqAlgorithm} = nothing, e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -316,8 +312,8 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th - `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. -- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector` or a `Tuple`. -- `alg`: The algorithm to use for the stochastic differential equation. Default is `SRA1()`. +- `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided. +- `alg`: The algorithm to use for the stochastic differential equation. Default is `SRIW1()` if `sc_ops` is an [`AbstractQuantumObject`](@ref) (diagonal noise), and `SRA1()` otherwise (non-diagonal noise). - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`. - `params`: `NullParameters` of parameters to pass to the solver. - `rng`: Random number generator for reproducibility. @@ -345,8 +341,8 @@ function ssesolve( H::Union{AbstractQuantumObject{OperatorQuantumObject},Tuple}, ψ0::QuantumObject{KetQuantumObject}, tlist::AbstractVector, - sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing; - alg::StochasticDiffEqAlgorithm = SRA1(), + sc_ops::Union{Nothing,AbstractVector,Tuple,AbstractQuantumObject} = nothing; + alg::Union{Nothing,StochasticDiffEqAlgorithm} = nothing, e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), @@ -375,6 +371,12 @@ function ssesolve( kwargs..., ) + sc_ops_isa_Qobj = sc_ops isa AbstractQuantumObject # We can avoid using non-diagonal noise if sc_ops is just an AbstractQuantumObject + + if isnothing(alg) + alg = sc_ops_isa_Qobj ? SRIW1() : SRA1() + end + return ssesolve(ens_prob, alg, ntraj, ensemblealg) end diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl index 737837358..aad79c00f 100644 --- a/src/time_evolution/time_evolution.jl +++ b/src/time_evolution/time_evolution.jl @@ -238,6 +238,9 @@ end ####################################### +_make_c_ops_list(c_ops::Union{AbstractVector,Tuple}) = c_ops +_make_c_ops_list(c_ops::AbstractQuantumObject) = (c_ops,) + function _merge_saveat(tlist, e_ops, default_options; kwargs...) is_empty_e_ops = isnothing(e_ops) ? true : isempty(e_ops) saveat = is_empty_e_ops ? tlist : [tlist[end]] @@ -347,13 +350,9 @@ function _stochastic_prob_func(prob, i, repeat, rng, seeds, tlist; kwargs...) traj_rng = typeof(rng)() seed!(traj_rng, seed) - noise = RealWienerProcess!( - prob.prob.tspan[1], - zeros(kwargs[:n_sc_ops]), - zeros(kwargs[:n_sc_ops]), - save_everystep = getVal(kwargs[:store_measurement]), - rng = traj_rng, - ) + sc_ops = kwargs[:sc_ops] + store_measurement = kwargs[:store_measurement] + noise = _make_noise(prob.prob.tspan[1], sc_ops, store_measurement, traj_rng) return remake(prob.prob, noise = noise, seed = seed) end @@ -361,6 +360,27 @@ end # Standard output function _stochastic_output_func(sol, i) = (sol, false) +#= + Define diagonal or non-diagonal noise depending on the type of `sc_ops`. + If `sc_ops` is a `AbstractQuantumObject`, we avoid using the non-diagonal noise. +=# +function _make_noise(t0::Real, sc_ops::Union{AbstractVector,Tuple}, store_measurement::Val, rng) + noise = RealWienerProcess!( + t0, + zeros(length(sc_ops)), + zeros(length(sc_ops)), + save_everystep = getVal(store_measurement), + rng = rng, + ) + + return noise +end +function _make_noise(t0::Real, sc_ops::AbstractQuantumObject, store_measurement::Val, rng) + noise = RealWienerProcess(t0, 0.0, 0.0, save_everystep = getVal(store_measurement), rng = rng) + + return noise +end + #= struct DiffusionOperator @@ -391,7 +411,7 @@ end N = length(ops_types) quote M = length(u) - S = size(v) + S = (size(v, 1), size(v, 2)) # This supports also `v` as a `Vector` (S[1] == M && S[2] == $N) || throw(DimensionMismatch("The size of the output vector is incorrect.")) Base.@nexprs $N i -> begin mul!(@view(v[:, i]), L.ops[i], u) diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl index 419c17677..55be066a3 100644 --- a/test/core-test/time_evolution.jl +++ b/test/core-test/time_evolution.jl @@ -21,6 +21,9 @@ sme_η = 0.7 # Efficiency of the homodyne detector for smesolve c_ops_sme = [sqrt(1 - sme_η) * op for op in c_ops] sc_ops_sme = [sqrt(sme_η) * op for op in c_ops] + # The following definition is to test the case of `sc_ops` as an `AbstractQuantumObject` + c_ops_sme2 = c_ops[2:end] + sc_ops_sme2 = c_ops[1] ψ0_int = Qobj(round.(Int, real.(ψ0.data)), dims = ψ0.dims) # Used for testing the type inference @@ -153,6 +156,7 @@ progress_bar = Val(false), store_measurement = Val(true), ) + sol_sme3 = smesolve(H, ψ0, tlist, c_ops_sme2, sc_ops_sme2, e_ops = e_ops, progress_bar = Val(false)) ρt_mc = [ket2dm.(normalize.(states)) for states in sol_mc_states.states] expect_mc_states = mapreduce(states -> expect.(Ref(e_ops[1]), states), hcat, ρt_mc) @@ -175,6 +179,7 @@ @test sum(abs, vec(expect_mc_states_mean2) .- vec(sol_me.expect[1, saveat_idxs])) / length(tlist) < 0.1 @test sum(abs, sol_sse.expect .- sol_me.expect) / length(tlist) < 0.1 @test sum(abs, sol_sme.expect .- sol_me.expect) / length(tlist) < 0.1 + @test sum(abs, sol_sme3.expect .- sol_me.expect) / length(tlist) < 0.1 @test length(sol_me.times) == length(tlist) @test length(sol_me.states) == 1 @test size(sol_me.expect) == (length(e_ops), length(tlist)) @@ -544,8 +549,11 @@ end @testset "Type Inference smesolve" begin - c_ops_sme_tuple = Tuple(c_ops_sme) # To avoid type instability, we must have a Tuple instead of a Vector - sc_ops_sme_tuple = Tuple(sc_ops_sme) # To avoid type instability, we must have a Tuple instead of a Vector + # To avoid type instability, we must have a Tuple instead of a Vector + c_ops_sme_tuple = Tuple(c_ops_sme) + sc_ops_sme_tuple = Tuple(sc_ops_sme) + c_ops_sme2_tuple = Tuple(c_ops_sme2) + sc_ops_sme2_tuple = sc_ops_sme2 # This is an `AbstractQuantumObject` @inferred smesolveEnsembleProblem( H, ψ0, @@ -568,6 +576,17 @@ progress_bar = Val(false), rng = rng, ) + @inferred smesolve( + H, + ψ0, + tlist, + c_ops_sme2_tuple, + sc_ops_sme2_tuple, + ntraj = 5, + e_ops = e_ops, + progress_bar = Val(false), + rng = rng, + ) @inferred smesolve( H, ψ0, From 9b20aeabc1d69bb470c4005ffd9bf9f3032b6b9a Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Sun, 16 Feb 2025 20:39:28 +0100 Subject: [PATCH 2/5] Make changelog and change default non-diagonal solver to `SRA2()` --- CHANGELOG.md | 3 +++ src/QuantumToolbox.jl | 2 +- src/time_evolution/smesolve.jl | 6 +++--- src/time_evolution/ssesolve.jl | 6 +++--- 4 files changed, 10 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b59a7bf4f..ce0837a02 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main) +- Support for single `AbstractQuantumObject` in `sc_ops` for faster specific method in `ssesolve` and `smesolve`. ([#408]) + ## [v0.27.0] Release date: 2025-02-14 @@ -140,3 +142,4 @@ Release date: 2024-11-13 [#403]: https://github.com/qutip/QuantumToolbox.jl/issues/403 [#404]: https://github.com/qutip/QuantumToolbox.jl/issues/404 [#405]: https://github.com/qutip/QuantumToolbox.jl/issues/405 +[#408]: https://github.com/qutip/QuantumToolbox.jl/issues/408 diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl index b88f8c519..e4835221b 100644 --- a/src/QuantumToolbox.jl +++ b/src/QuantumToolbox.jl @@ -40,7 +40,7 @@ import SciMLBase: AbstractSciMLProblem, AbstractODEIntegrator, AbstractODESolution -import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA1, SRIW1 +import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA2, SRIW1 import SciMLOperators: SciMLOperators, AbstractSciMLOperator, diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl index f5f52d024..60700b66f 100644 --- a/src/time_evolution/smesolve.jl +++ b/src/time_evolution/smesolve.jl @@ -313,7 +313,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio - `tlist`: List of times at which to save either the state or the expectation values of the system. - `c_ops`: List of collapse operators ``\{\hat{C}_i\}_i``. It can be either a `Vector` or a `Tuple`. - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided. -- `alg`: The algorithm to use for the stochastic differential equation. Default is `SRIW1()` if `sc_ops` is an [`AbstractQuantumObject`](@ref) (diagonal noise), and `SRA1()` otherwise (non-diagonal noise). +- `alg`: The algorithm to use for the stochastic differential equation. Default is `SRIW1()` if `sc_ops` is an [`AbstractQuantumObject`](@ref) (diagonal noise), and `SRA2()` otherwise (non-diagonal noise). - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`. - `params`: `NullParameters` of parameters to pass to the solver. - `rng`: Random number generator for reproducibility. @@ -375,7 +375,7 @@ function smesolve( sc_ops_isa_Qobj = sc_ops isa AbstractQuantumObject # We can avoid using non-diagonal noise if sc_ops is just an AbstractQuantumObject if isnothing(alg) - alg = sc_ops_isa_Qobj ? SRIW1() : SRA1() + alg = sc_ops_isa_Qobj ? SRIW1() : SRA2() end return smesolve(ensemble_prob, alg, ntraj, ensemblealg) @@ -383,7 +383,7 @@ end function smesolve( ens_prob::TimeEvolutionProblem, - alg::StochasticDiffEqAlgorithm = SRA1(), + alg::StochasticDiffEqAlgorithm = SRA2(), ntraj::Int = 500, ensemblealg::EnsembleAlgorithm = EnsembleThreads(), ) diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl index 264ff32f6..9127c63e4 100644 --- a/src/time_evolution/ssesolve.jl +++ b/src/time_evolution/ssesolve.jl @@ -313,7 +313,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th - `ψ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. - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided. -- `alg`: The algorithm to use for the stochastic differential equation. Default is `SRIW1()` if `sc_ops` is an [`AbstractQuantumObject`](@ref) (diagonal noise), and `SRA1()` otherwise (non-diagonal noise). +- `alg`: The algorithm to use for the stochastic differential equation. Default is `SRIW1()` if `sc_ops` is an [`AbstractQuantumObject`](@ref) (diagonal noise), and `SRA2()` otherwise (non-diagonal noise). - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`. - `params`: `NullParameters` of parameters to pass to the solver. - `rng`: Random number generator for reproducibility. @@ -374,7 +374,7 @@ function ssesolve( sc_ops_isa_Qobj = sc_ops isa AbstractQuantumObject # We can avoid using non-diagonal noise if sc_ops is just an AbstractQuantumObject if isnothing(alg) - alg = sc_ops_isa_Qobj ? SRIW1() : SRA1() + alg = sc_ops_isa_Qobj ? SRIW1() : SRA2() end return ssesolve(ens_prob, alg, ntraj, ensemblealg) @@ -382,7 +382,7 @@ end function ssesolve( ens_prob::TimeEvolutionProblem, - alg::StochasticDiffEqAlgorithm = SRA1(), + alg::StochasticDiffEqAlgorithm = SRA2(), ntraj::Int = 500, ensemblealg::EnsembleAlgorithm = EnsembleThreads(), ) From 7980f916de61345bf4c1d78f69ecdeb672f72a0e Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Sun, 16 Feb 2025 21:36:42 +0100 Subject: [PATCH 3/5] Fix Code Quality --- src/time_evolution/smesolve.jl | 2 +- src/time_evolution/ssesolve.jl | 2 +- src/time_evolution/time_evolution.jl | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl index 60700b66f..a9d151d7f 100644 --- a/src/time_evolution/smesolve.jl +++ b/src/time_evolution/smesolve.jl @@ -127,7 +127,7 @@ function smesolveProblem( ) tspan = (tlist[1], tlist[end]) - noise = _make_noise(tspan[1], sc_ops, store_measurement, rng) + noise = _make_noise(tspan[1], sc_ops, makeVal(store_measurement), rng) noise_rate_prototype = sc_ops_isa_Qobj ? nothing : similar(ρ0, length(ρ0), length(sc_ops_list)) prob = SDEProblem{true}( K, diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl index 9127c63e4..5b622335f 100644 --- a/src/time_evolution/ssesolve.jl +++ b/src/time_evolution/ssesolve.jl @@ -127,7 +127,7 @@ function ssesolveProblem( ) tspan = (tlist[1], tlist[end]) - noise = _make_noise(tspan[1], sc_ops, store_measurement, rng) + noise = _make_noise(tspan[1], sc_ops, makeVal(store_measurement), rng) noise_rate_prototype = sc_ops_isa_Qobj ? nothing : similar(ψ0, length(ψ0), length(sc_ops_list)) prob = SDEProblem{true}( K, diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl index aad79c00f..b487d6511 100644 --- a/src/time_evolution/time_evolution.jl +++ b/src/time_evolution/time_evolution.jl @@ -238,7 +238,7 @@ end ####################################### -_make_c_ops_list(c_ops::Union{AbstractVector,Tuple}) = c_ops +_make_c_ops_list(c_ops) = c_ops _make_c_ops_list(c_ops::AbstractQuantumObject) = (c_ops,) function _merge_saveat(tlist, e_ops, default_options; kwargs...) @@ -364,7 +364,7 @@ _stochastic_output_func(sol, i) = (sol, false) Define diagonal or non-diagonal noise depending on the type of `sc_ops`. If `sc_ops` is a `AbstractQuantumObject`, we avoid using the non-diagonal noise. =# -function _make_noise(t0::Real, sc_ops::Union{AbstractVector,Tuple}, store_measurement::Val, rng) +function _make_noise(t0, sc_ops, store_measurement::Val, rng) noise = RealWienerProcess!( t0, zeros(length(sc_ops)), @@ -375,7 +375,7 @@ function _make_noise(t0::Real, sc_ops::Union{AbstractVector,Tuple}, store_measur return noise end -function _make_noise(t0::Real, sc_ops::AbstractQuantumObject, store_measurement::Val, rng) +function _make_noise(t0, sc_ops::AbstractQuantumObject, store_measurement::Val, rng) noise = RealWienerProcess(t0, 0.0, 0.0, save_everystep = getVal(store_measurement), rng = rng) return noise From e083c31900e8a6e7b43008aedfdfdea5a45b07f9 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Mon, 17 Feb 2025 10:01:46 +0100 Subject: [PATCH 4/5] Add a note in the docstrings --- src/time_evolution/smesolve.jl | 9 +++++++++ src/time_evolution/ssesolve.jl | 9 +++++++++ 2 files changed, 18 insertions(+) diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl index a9d151d7f..157b0b42f 100644 --- a/src/time_evolution/smesolve.jl +++ b/src/time_evolution/smesolve.jl @@ -65,6 +65,9 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio - The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`. - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) +!!! tip "Performance Tip" + When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + # Returns - `prob`: The [`TimeEvolutionProblem`](@ref) containing the `SDEProblem` for the Stochastic Master Equation time evolution. @@ -207,6 +210,9 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio - The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`. - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) +!!! tip "Performance Tip" + When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + # Returns - `prob`: The [`TimeEvolutionProblem`](@ref) containing the Ensemble `SDEProblem` for the Stochastic Master Equation time evolution. @@ -332,6 +338,9 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio - The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`. - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) +!!! tip "Performance Tip" + When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + # Returns - `sol::TimeEvolutionStochasticSol`: The solution of the time evolution. See [`TimeEvolutionStochasticSol`](@ref). diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl index 5b622335f..21b65c4fa 100644 --- a/src/time_evolution/ssesolve.jl +++ b/src/time_evolution/ssesolve.jl @@ -67,6 +67,9 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th - The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`. - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) +!!! tip "Performance Tip" + When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + # Returns - `prob`: The `SDEProblem` for the Stochastic Schrödinger time evolution of the system. @@ -207,6 +210,9 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is t - The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`. - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) +!!! tip "Performance Tip" + When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + # Returns - `prob::EnsembleProblem with SDEProblem`: The Ensemble SDEProblem for the Stochastic Shrödinger time evolution. @@ -333,6 +339,9 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th - For more details about `alg` please refer to [`DifferentialEquations.jl` (SDE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/) - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) +!!! tip "Performance Tip" + When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + # Returns - `sol::TimeEvolutionStochasticSol`: The solution of the time evolution. See [`TimeEvolutionStochasticSol`](@ref). From 4029a5bfdf66568973a875306954845645df5011 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Mon, 17 Feb 2025 10:07:17 +0100 Subject: [PATCH 5/5] Make the note more clear --- src/time_evolution/smesolve.jl | 6 +++--- src/time_evolution/ssesolve.jl | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl index 157b0b42f..4890030d2 100644 --- a/src/time_evolution/smesolve.jl +++ b/src/time_evolution/smesolve.jl @@ -66,7 +66,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) !!! tip "Performance Tip" - When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + When `sc_ops` contains only a single operator, it is recommended to pass only that operator as the argument. This ensures that the stochastic noise is diagonal, making the simulation faster. # Returns @@ -211,7 +211,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) !!! tip "Performance Tip" - When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + When `sc_ops` contains only a single operator, it is recommended to pass only that operator as the argument. This ensures that the stochastic noise is diagonal, making the simulation faster. # Returns @@ -339,7 +339,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) !!! tip "Performance Tip" - When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + When `sc_ops` contains only a single operator, it is recommended to pass only that operator as the argument. This ensures that the stochastic noise is diagonal, making the simulation faster. # Returns diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl index 21b65c4fa..5d7e89d0b 100644 --- a/src/time_evolution/ssesolve.jl +++ b/src/time_evolution/ssesolve.jl @@ -68,7 +68,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) !!! tip "Performance Tip" - When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + When `sc_ops` contains only a single operator, it is recommended to pass only that operator as the argument. This ensures that the stochastic noise is diagonal, making the simulation faster. # Returns @@ -211,7 +211,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is t - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) !!! tip "Performance Tip" - When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + When `sc_ops` contains only a single operator, it is recommended to pass only that operator as the argument. This ensures that the stochastic noise is diagonal, making the simulation faster. # Returns @@ -340,7 +340,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) !!! tip "Performance Tip" - When `sc_ops` contains only one operator, it is highly recommended to put only the operator as the argument. This will ensure the stochastic noise to be diagonal, allowing for a faster simulation. + When `sc_ops` contains only a single operator, it is recommended to pass only that operator as the argument. This ensures that the stochastic noise is diagonal, making the simulation faster. # Returns