From 1b5c14e0947e9a26e4e1af09454cd17b716da2b9 Mon Sep 17 00:00:00 2001 From: Gavin Date: Sun, 3 Aug 2025 23:38:46 -0400 Subject: [PATCH 1/9] playing with the ensemble_sesolve function to make it more efficient. Currently cannot do multiple parameters, but can do multiple initial states along with being type stable. --- Project.toml | 2 + src/QuantumToolbox.jl | 1 + src/time_evolution/ensemble_sesolve.jl | 92 ++++++++++++++++++++++++++ 3 files changed, 95 insertions(+) create mode 100644 src/time_evolution/ensemble_sesolve.jl diff --git a/Project.toml b/Project.toml index 60f71f7c3..5ce8f99bb 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +ProtoStructs = "437b6fc4-8e8e-11e9-3fa1-ad391e66c018" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" @@ -60,6 +61,7 @@ Makie = "0.24" OrdinaryDiffEqCore = "1" OrdinaryDiffEqTsit5 = "1" Pkg = "1" +ProtoStructs = "1.2.1" Random = "1" SciMLBase = "2.105" SciMLOperators = "1.4" diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl index c3d9be754..3b86d897c 100644 --- a/src/QuantumToolbox.jl +++ b/src/QuantumToolbox.jl @@ -117,6 +117,7 @@ include("time_evolution/mcsolve.jl") include("time_evolution/ssesolve.jl") include("time_evolution/smesolve.jl") include("time_evolution/time_evolution_dynamical.jl") +include("time_evolution/ensemble_sesolve.jl") # Others include("correlations.jl") diff --git a/src/time_evolution/ensemble_sesolve.jl b/src/time_evolution/ensemble_sesolve.jl new file mode 100644 index 000000000..36b20f32b --- /dev/null +++ b/src/time_evolution/ensemble_sesolve.jl @@ -0,0 +1,92 @@ +export EnsembleTimeEvolutionProblem +""" + EnsembleTimeEvolutionProblem{PT<:TimeEvolutionProblem, PF<:Function} + +A structure representing an ensemble time evolution problem for quantum systems. + +# Fields +- `prob::PT`: The base time evolution problem. +- `func::PF`: A function used to modify or sample parameters for each trajectory in the ensemble. +- `iterate_params::Bool`: If `true`, parameters are iterated for each trajectory; otherwise, the same parameters are used. +- `full_iterator::AbstractArray`: An array containing all parameter sets or states to be used in the ensemble. +- `n_states::Int`: The number of initial states. +- `trajectories::Int`: The total number of trajectories to simulate. + +# Usage +This is used when setting up ensemble sesolve problems, useful for simulating multiple quantum states or parameter sets in parallel. + +Example: +```julia + H = 2 * π * 0.1 * sigmax() + ψ0 = basis(2, 0) # spin-up + tlist = LinRange(0.0, 100.0, 100) + + ψs = [ψ0, basis(2, 1)] # spin-up and spin-down + + params = collect(Iterators.product([0,1,2,3,4,5], [0,1,2,3,4,5], [0,1,2,3,4,5])) + res = sesolve(H, ψs, tlist; params = params, iterate_params = true, alg = Tsit5(), progress_bar=false); +``` +""" +struct EnsembleTimeEvolutionProblem{PT<:TimeEvolutionProblem,PF<:Function} + prob::PT + func::PF + full_iterator::AbstractArray + trajectories::Int +end + + +function sesolve(prob::EnsembleTimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5(); backend = EnsembleThreads()) + ensemble_prob = EnsembleProblem(prob.prob.prob, prob_func = prob.func) + sols = solve(ensemble_prob, alg, backend, trajectories = prob.trajectories) + + to_return = Array{TimeEvolutionSol}(undef, size(prob.full_iterator)) + for i in 1:length(sols) + ψt = map(ϕ -> QuantumObject(ϕ, type = Ket(), dims = prob.prob.dimensions), sols[i].u) + sol = TimeEvolutionSol( + prob.prob.times, + sols[i].t, + ψt, + _get_expvals(sols[i], SaveFuncSESolve), + sols[i].retcode, + sols[i].alg, + sols[i].prob.kwargs[:abstol], + sols[i].prob.kwargs[:reltol], + ) + to_return[CartesianIndices(to_return)[i]] = sol + end + + return to_return +end + +function sesolve( + H::Union{AbstractQuantumObject{Operator},Tuple}, + ψ0s::Vector{T}, + tlist::AbstractVector; + alg::OrdinaryDiffEqAlgorithm = Tsit5(), + e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, + params = NullParameters(), + progress_bar::Union{Val,Bool} = Val(false), + inplace::Union{Val,Bool} = Val(true), + backend = EnsembleThreads(), + kwargs...,) where T<:QuantumObject{Ket} + + prob_init = sesolveProblem( + H, + ψ0s[1], + tlist; + e_ops = e_ops, + params = params, + progress_bar = progress_bar, + inplace = inplace, + kwargs..., + ) + + trajectories = length(ψ0s) + + function ensemble_func(prob, i, repeat) + return remake(prob, u0 = ψ0s[i].data) + end + + ensemble_prob = EnsembleTimeEvolutionProblem(prob_init, ensemble_func, ψ0s, trajectories) + return sesolve(ensemble_prob, alg; backend = backend) +end \ No newline at end of file From 1746ffce49d4e409cf6aaea61b5576d2ad81020e Mon Sep 17 00:00:00 2001 From: Gavin Date: Mon, 4 Aug 2025 13:10:02 -0400 Subject: [PATCH 2/9] hmm --- src/time_evolution/ensemble_sesolve.jl | 37 +++++++++++++++++++++----- 1 file changed, 30 insertions(+), 7 deletions(-) diff --git a/src/time_evolution/ensemble_sesolve.jl b/src/time_evolution/ensemble_sesolve.jl index 36b20f32b..7a7a97209 100644 --- a/src/time_evolution/ensemble_sesolve.jl +++ b/src/time_evolution/ensemble_sesolve.jl @@ -27,19 +27,37 @@ Example: res = sesolve(H, ψs, tlist; params = params, iterate_params = true, alg = Tsit5(), progress_bar=false); ``` """ -struct EnsembleTimeEvolutionProblem{PT<:TimeEvolutionProblem,PF<:Function} +struct EnsembleTimeEvolutionProblem{PT<:TimeEvolutionProblem,PF<:Function, X<:Vector{T} where T<:QuantumObject{Ket}, Y<:AbstractArray} prob::PT func::PF - full_iterator::AbstractArray + states::X + params::Y + problem_dims::Tuple trajectories::Int end +function EnsembleTimeEvolutionProblem( + prob::PT, + states::Vector{T}, + params::AbstractArray = [NullParameters()] + ) where {PT<:TimeEvolutionProblem, T<:QuantumObject{Ket}} + + problem_dims = (length(states), length(params)) + + function ensemble_func(prob, i, repeat) + state_id = mod1(i, problem_dims[1]) + param_id = div(i - 1, problem_dims[1]) + 1 + return remake(prob, u0 = states[state_id].data, p = params[param_id]) + end + trajectories = prod(problem_dims) + return EnsembleTimeEvolutionProblem(prob, ensemble_func, states, params, problem_dims, trajectories) +end function sesolve(prob::EnsembleTimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5(); backend = EnsembleThreads()) ensemble_prob = EnsembleProblem(prob.prob.prob, prob_func = prob.func) sols = solve(ensemble_prob, alg, backend, trajectories = prob.trajectories) - to_return = Array{TimeEvolutionSol}(undef, size(prob.full_iterator)) + to_return = Array{TimeEvolutionSol}(undef, prob.problem_dims) for i in 1:length(sols) ψt = map(ϕ -> QuantumObject(ϕ, type = Ket(), dims = prob.prob.dimensions), sols[i].u) sol = TimeEvolutionSol( @@ -83,10 +101,15 @@ function sesolve( trajectories = length(ψ0s) - function ensemble_func(prob, i, repeat) - return remake(prob, u0 = ψ0s[i].data) - end + # function ensemble_func(prob, i, repeat) + # return remake(prob, u0 = ψ0s[i].data) + # end - ensemble_prob = EnsembleTimeEvolutionProblem(prob_init, ensemble_func, ψ0s, trajectories) + # ensemble_prob = EnsembleTimeEvolutionProblem(prob_init, ensemble_func, ψ0s, trajectories) + ensemble_prob = EnsembleTimeEvolutionProblem( + prob_init, + ψ0s, + [params] + ) return sesolve(ensemble_prob, alg; backend = backend) end \ No newline at end of file From ad467ce8940ed89b19689fd83ebddf5febfb3518 Mon Sep 17 00:00:00 2001 From: Gavin Date: Mon, 1 Sep 2025 16:13:41 -0500 Subject: [PATCH 3/9] Added propagator.jl and qeye_like.jl --- other_missing_things.md | 48 ++++++++++ propagator_test.jl | 14 +++ src/QuantumToolbox.jl | 3 +- src/qobj/qeye_like.jl | 35 ++++++++ src/time_evolution/propagator.jl | 147 +++++++++++++++++++++++++++++++ 5 files changed, 246 insertions(+), 1 deletion(-) create mode 100644 other_missing_things.md create mode 100644 propagator_test.jl create mode 100644 src/qobj/qeye_like.jl create mode 100644 src/time_evolution/propagator.jl diff --git a/other_missing_things.md b/other_missing_things.md new file mode 100644 index 000000000..cf334ef06 --- /dev/null +++ b/other_missing_things.md @@ -0,0 +1,48 @@ +# Things that are missing from the QuantumToolbox.jl package. + +1) Powers are not supported for `QobjEvo` objects. + - There is a bug with the one I did, as the following returns a stackoverflow error. + ```julia + import Base.^ + import Base.^ + function ^(x::QobjEvo, y::Integer) + if y< 0 + error("Exponent must be non-negative integer") + end + if y == 0 + return QobjEvo(qeye(size(x)[1], dims = x.dims)) + end + to_return = x + for i in 2:y + to_return = to_return * x + end + return to_return + end + + 1*(sum([(QobjEvo((qeye(2), (p,t) -> 1)))^n*1 for n in 0:1:2])) + 1*(sum([1*(QobjEvo((qeye(2), (p,t) -> 1)))^n for n in 0:1:2])) + ``` + - However, the following do not + ```julia + 1*(sum([(QobjEvo((qeye(2), (p,t) -> 1)))^n for n in 0:1:2])) + 1*(sum([(1*QobjEvo((qeye(2), (p,t) -> 1)))^n for n in 0:1:2])) + 1*(sum([(QobjEvo((qeye(2), (p,t) -> 1))*1)^n for n in 0:1:2])) + ``` + + The error comes from expressions like + ```julia + a = QobjEvo((qeye(2), (p,t) -> 1)) + 1*(a+a*a*1) + ``` + + The error is: + StackOverflowError: + ```julia + StackOverflowError: + + Stacktrace: + [1] *(α::ScalarOperator{Int64, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, x::UniformScaling{Int64}) + @ SciMLOperators ~/.julia/packages/SciMLOperators/DKrpP/src/scalar.jl:366 + [2] *(α::ScalarOperator{Int64, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, x::UniformScaling{Int64}) + @ SciMLOperators ~/.julia/packages/SciMLOperators/DKrpP/src/scalar.jl:369 + ``` \ No newline at end of file diff --git a/propagator_test.jl b/propagator_test.jl new file mode 100644 index 000000000..41fa96cd1 --- /dev/null +++ b/propagator_test.jl @@ -0,0 +1,14 @@ +Pkg.activate(".") +include("src/QuantumToolbox.jl") +using .QuantumToolbox +using LinearAlgebra + +using Revise + + +p_se = propagator(sigmax()) +p_me = propagator(sigmax(), [sigmaz()]) + +println(p_se(1.0)) +println("-------------------") +println(p_me(2, t0 = 1)) \ No newline at end of file diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl index 3b86d897c..01d725b02 100644 --- a/src/QuantumToolbox.jl +++ b/src/QuantumToolbox.jl @@ -100,7 +100,7 @@ include("qobj/operators.jl") include("qobj/superoperators.jl") include("qobj/synonyms.jl") include("qobj/block_diagonal_form.jl") - +include("qobj/qeye_like.jl") # time evolution include("time_evolution/time_evolution.jl") include("time_evolution/callback_helpers/callback_helpers.jl") @@ -118,6 +118,7 @@ include("time_evolution/ssesolve.jl") include("time_evolution/smesolve.jl") include("time_evolution/time_evolution_dynamical.jl") include("time_evolution/ensemble_sesolve.jl") +include("time_evolution/propagator.jl") # Others include("correlations.jl") diff --git a/src/qobj/qeye_like.jl b/src/qobj/qeye_like.jl new file mode 100644 index 000000000..36765f27d --- /dev/null +++ b/src/qobj/qeye_like.jl @@ -0,0 +1,35 @@ +export qeye_like + +function eye(x :: T) where T<:Union{AbstractQuantumObject{Operator}, AbstractQuantumObject{SuperOperator}} + qeye_like(x) +end +function qeye_like(x :: T) where T<:QuantumObject{Operator} + y = 0*x + for i in 1:size(x)[1] + y[i,i] = 1 + end + return y +end + +function qeye_like(x :: T) where T<:QobjEvo + y = 0*x(0) + for i in 1:size(x)[1] + y[i,i] = 1 + end + return y +end + +function qeye_like(x :: T) where T<:QuantumObject{Ket} + qeye_like(x*x') +end +function qeye_like(x :: T) where T<:QuantumObject{Bra} + qeye_like(x'*x) +end + +function qeye_like(x::T) where T<:QuantumObject{SuperOperator} + y = 0*x + for i in 1:size(x)[1] + y[i,i] = 1 + end + return y +end \ No newline at end of file diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl new file mode 100644 index 000000000..35154b65b --- /dev/null +++ b/src/time_evolution/propagator.jl @@ -0,0 +1,147 @@ +export Propagator +export propagator + +function propagatorProblem( + H::AbstractQuantumObject{Operator}, + tspan::AbstractVector; + params = NullParameters(), + progress_bar::Union{Val, Bool} = Val(true), + inplace::Union{Val, Bool} = Val(true), + kwargs... +) + H_evo = QobjEvo(H) + + p0 = qeye(H).data + U = -1im*H_evo.data + + kwargs2 = _merge_saveat(tspan, nothing, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...) + kwargs3 = _generate_se_me_kwargs(nothing, makeVal(progress_bar), tspan, kwargs2, SaveFuncSESolve) + + return ODEProblem{getVal(inplace),FullSpecialize}(U, p0, tspan, params; kwargs3...) +end + +function propagatorProblem( + H::AbstractQuantumObject{SuperOperator}, + tspan::AbstractVector; + params = NullParameters(), + progress_bar::Union{Val, Bool} = Val(true), + inplace::Union{Val, Bool} = Val(true), + kwargs... +) + H_evo = QobjEvo(H) + + p0 = qeye(H).data + U = H_evo.data + + kwargs2 = _merge_saveat(tspan, nothing, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...) + kwargs3 = _generate_se_me_kwargs(nothing, makeVal(progress_bar), tspan, kwargs2, SaveFuncMESolve) + + return ODEProblem{getVal(inplace),FullSpecialize}(U, p0, tspan, params; kwargs3...) +end + +struct Propagator{T<:Union{AbstractQuantumObject{Operator}, AbstractQuantumObject{SuperOperator}}} + H::T + times::AbstractArray{Float64} + props::AbstractArray{A} where A <: AbstractQuantumObject + tol::Float64 + solver_kwargs::Base.Pairs + dims::Union{AbstractArray, Tuple} + solver::OrdinaryDiffEqAlgorithm + type +end + +function Base.show(io::IO, p::Propagator) + println("Propagator: ") + println(" dims: $(p.dims)") + println(" issuper: $(issuper(p.H))") + println(" size: $(size(p.H))") + println(" Number of Saved Times: $(length(p.times))") +end + +function propagator( + H::Union{QobjEvo, QuantumObject}; + tol = 1e-6, + solver = Tsit5(), + kwargs... + ) + return Propagator(H, [0.0], [qeye_like(H)], tol, kwargs, H.dims, solver, Operator()) +end + +function propagator( + H::Union{QobjEvo, QuantumObject}, + c_ops :: Union{AbstractArray, Tuple}; + tol = 1e-6, + solver = Tsit5(), + kwargs... + ) + L = liouvillian(H, c_ops) + return Propagator(L, [0.0], [L], tol, kwargs, L.dims, solver, SuperOperator()) +end + +function _lookup_or_compute(p::Propagator, t::Real) + """ + Get U(t) from cache or compute it. + """ + idx = searchsorted(p.times, t) + + # Check if we found an exact match or close enough within tolerance + if (idx.start <= length(p.times)) && (abs(t - p.times[idx.start]) <= p.tol) + U = p.props[idx.start] + elseif (idx.start > 1) && (abs(t - p.times[idx.start - 1]) <= p.tol) + U = p.props[idx.start - 1] + else + t0 = (idx.start > 1) ? p.times[idx.start - 1] : 0.0 + U = _compute_propagator(p, t) + insert!(p.props, idx.start, U) + insert!(p.times, idx.start, t) + end + + return U +end + +function _compute_propagator(p::Propagator, + t::Real; + t0 = 0.0 + ) + prob = propagatorProblem(p.H, [t0, t]; p.solver_kwargs...) + res = solve(prob, p.solver) + U = QuantumObject(res.u[end], dims = p.dims, type = p.type) + return U +end + + +""" + (p::Propagator)(t; t0=0.0, just_interval=false) + +Evaluate the propagator at time t. + +- t: target time. +- t0: initial time. If nonzero: + - when just_interval == true, integrate directly on [t0, t] to obtain U(t, t0). + - otherwise, compute U(t) * inv(U(t0)) using cached values when available. +- just_interval: force a fresh integration on [t0, t] instead of composing from cache. This is useful if both t0 and t are large but t-t0 is small. + +Returns a QuantumObject with the same dims/type as p.H. Results are cached at requested +times (within tolerance p.tol) to avoid recomputation. +""" +function (p::Propagator)(t::Real; t0 = 0.0, just_interval = false) + if just_interval + return _compute_propagator(p, t; t0 = t0) + end + # We start by computing U_t0 if needed. This means that the computation for U will, at worst, start at t0 instead of 0.0. + if t0 != 0.0 + U_t0 = (_lookup_or_compute(p, t0)) + end + U = _lookup_or_compute(p, t) + + # We only do the multiplication if needed + if t0 != 0.0 + U = U*inv(U_t0) + end + return U +end + +function Base.size(p::Propagator) + return size(p.H) +end + From f86e1ef9b4f22b7f5eb5a46fd1ab5e5a34d8dbf6 Mon Sep 17 00:00:00 2001 From: Gavin Date: Sat, 6 Sep 2025 19:03:55 +0000 Subject: [PATCH 4/9] made a few changes, including removing qeye_like because I discovered that was already implimented as `one`. --- .gitignore | 3 +- Project.toml | 2 + docs/make.jl | 2 +- propagator_test.jl | 3 +- src/QuantumToolbox.jl | 1 - src/qobj/qeye_like.jl | 35 ---------- src/time_evolution/propagator.jl | 113 +++++++++++++++++++++++++------ 7 files changed, 101 insertions(+), 58 deletions(-) delete mode 100644 src/qobj/qeye_like.jl diff --git a/.gitignore b/.gitignore index b07886a94..0d3f129f9 100644 --- a/.gitignore +++ b/.gitignore @@ -10,4 +10,5 @@ Manifest.toml benchmarks/benchmarks_output.json .ipynb_checkpoints -*.ipynb \ No newline at end of file +*.ipynb +.devcontainer \ No newline at end of file diff --git a/Project.toml b/Project.toml index 5ce8f99bb..f9b61278c 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,7 @@ OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" ProtoStructs = "437b6fc4-8e8e-11e9-3fa1-ad391e66c018" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -63,6 +64,7 @@ OrdinaryDiffEqTsit5 = "1" Pkg = "1" ProtoStructs = "1.2.1" Random = "1" +Revise = "3.9.0" SciMLBase = "2.105" SciMLOperators = "1.4" SparseArrays = "1" diff --git a/docs/make.jl b/docs/make.jl index d9ac8be98..028d99cac 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -19,7 +19,7 @@ DocMeta.setdocmeta!(QuantumToolbox, :DocTestSetup, doctest_setup; recursive = tr # some options for `makedocs` const DRAFT = get(ENV, "DRAFT", false) == "true" # `DRAFT = true` disables cell evaluation -const DOCTEST = get(ENV, "DOCTEST", true) == true # `DOCTEST = false` skips doc tests +const DOCTEST = get(ENV, "DOCTEST", true) == false # `DOCTEST = false` skips doc tests # generate bibliography bib = CitationBibliography( diff --git a/propagator_test.jl b/propagator_test.jl index 41fa96cd1..486fb6fbc 100644 --- a/propagator_test.jl +++ b/propagator_test.jl @@ -1,3 +1,4 @@ +using Pkg Pkg.activate(".") include("src/QuantumToolbox.jl") using .QuantumToolbox @@ -11,4 +12,4 @@ p_me = propagator(sigmax(), [sigmaz()]) println(p_se(1.0)) println("-------------------") -println(p_me(2, t0 = 1)) \ No newline at end of file +println(p_me(2, 0.1)) \ No newline at end of file diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl index 01d725b02..d23c3242e 100644 --- a/src/QuantumToolbox.jl +++ b/src/QuantumToolbox.jl @@ -100,7 +100,6 @@ include("qobj/operators.jl") include("qobj/superoperators.jl") include("qobj/synonyms.jl") include("qobj/block_diagonal_form.jl") -include("qobj/qeye_like.jl") # time evolution include("time_evolution/time_evolution.jl") include("time_evolution/callback_helpers/callback_helpers.jl") diff --git a/src/qobj/qeye_like.jl b/src/qobj/qeye_like.jl deleted file mode 100644 index 36765f27d..000000000 --- a/src/qobj/qeye_like.jl +++ /dev/null @@ -1,35 +0,0 @@ -export qeye_like - -function eye(x :: T) where T<:Union{AbstractQuantumObject{Operator}, AbstractQuantumObject{SuperOperator}} - qeye_like(x) -end -function qeye_like(x :: T) where T<:QuantumObject{Operator} - y = 0*x - for i in 1:size(x)[1] - y[i,i] = 1 - end - return y -end - -function qeye_like(x :: T) where T<:QobjEvo - y = 0*x(0) - for i in 1:size(x)[1] - y[i,i] = 1 - end - return y -end - -function qeye_like(x :: T) where T<:QuantumObject{Ket} - qeye_like(x*x') -end -function qeye_like(x :: T) where T<:QuantumObject{Bra} - qeye_like(x'*x) -end - -function qeye_like(x::T) where T<:QuantumObject{SuperOperator} - y = 0*x - for i in 1:size(x)[1] - y[i,i] = 1 - end - return y -end \ No newline at end of file diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl index 35154b65b..af2c0ee72 100644 --- a/src/time_evolution/propagator.jl +++ b/src/time_evolution/propagator.jl @@ -11,7 +11,7 @@ function propagatorProblem( ) H_evo = QobjEvo(H) - p0 = qeye(H).data + p0 = one(H).data U = -1im*H_evo.data kwargs2 = _merge_saveat(tspan, nothing, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...) @@ -30,7 +30,7 @@ function propagatorProblem( ) H_evo = QobjEvo(H) - p0 = qeye(H).data + p0 = one(H).data U = H_evo.data kwargs2 = _merge_saveat(tspan, nothing, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...) @@ -39,6 +39,54 @@ function propagatorProblem( return ODEProblem{getVal(inplace),FullSpecialize}(U, p0, tspan, params; kwargs3...) end +@doc raw""" +Propagator{T} + +Container for time-evolution propagators built from `H` (a Hamiltonian or a Louivillian superoperator). +This immutable struct groups the generator, time list and solver settings together with any cached +propagators so that previous results can be reused where applicable. + +# Parameters +- `H::T` + The generator of dynamics. `T` is typically an `AbstractQuantumObject{Operator}` for closed, + unitary dynamics (Hamiltonian) or an `AbstractQuantumObject{SuperOperator}` for open-system + dynamics (Liouvillian or other superoperator). +- `times::AbstractArray{Float64}` + A monotonic array of time points at which propagators are requested or saved. Values are + interpreted in the same time units used to build `H`. +- `props::AbstractArray{A} where A <: AbstractQuantumObject` + An array holding computed propagators corresponding to `times`. Elements are quantum objects + (operators or superoperators) matching the expected output of propagating with `H`. Both `times` and `props` are initialized with the identity operator/superoperator at t=0. +- `tol::Float64` + This is the time tolerance used to determine if a previous result should be reused. If \$t - t_{\text{cached}} < \text{tol}\$, where $t_\text{cached}\$ is the closed cached + time to \$t\$, then the propagator at \$t_\text{cached}\$ will be reused. +- `solver_kwargs::Base.Pairs` + Additional keyword arguments forwarded to the ODE integrator (for example `abstol`, `reltol`, + `saveat`, or solver-specific options). Stored as a `Base.Pairs` collection for convenient + dispatch into solver APIs. +- `dims::Union{AbstractArray, Tuple}` + Dimension metadata of `H`. +- `solver::OrdinaryDiffEqAlgorithm` + The chosen ODE algorithm (a concrete algorithm type from OrdinaryDiffEq, the default is `Tsit5`) that will be used to integrate the generator to obtain + time-ordered evolution. +- `type` + Whether the propagator is an operator or superoperator. + +# Initialization +A propagator is initialized by calling the `propagator` function. + +# Calling +If `U` is a propagator object, to get the propagator between t0 and t, call `U(t, t0=t0; just_interval::Bool=false)`. + - t0 is optional, if left out then t0=0.0. + - `just_interval` determines how the propagator from t0 to t is calculated. + - if `just_interval = false`, then the propagator is calculated by getting U1(t0,0.0) and U2(t, 0.0) and setting + ```math + U = U2U1^{-1} + ``` + This is useful if U1 is alread cached or one plans to reuse U1 or U2. + - if `just_interval=true` then the propagator is calulated by directly integrating from t0 to t. This approach does't save + anything but is useful if t-t0 >> t0 as it avoids the long integration required to get t0. +""" struct Propagator{T<:Union{AbstractQuantumObject{Operator}, AbstractQuantumObject{SuperOperator}}} H::T times::AbstractArray{Float64} @@ -58,13 +106,54 @@ function Base.show(io::IO, p::Propagator) println(" Number of Saved Times: $(length(p.times))") end +@doc raw""" +propagator(H::Union{QobjEvo, QuantumObject}; tol=1e-6, solver=Tsit5(), kwargs...) + +Construct and return a Propagator object prepared for time evolution using the +given Hamiltonian or time-dependent operator. + +# Arguments +- H::Union{QobjEvo, QuantumObject} + The Hamiltonian or generator for the dynamics. Can be a time-independent + QuantumObject or a time-dependent QobjEvo describing H(t). +- c_ops::Union{AbstractArray, Tuple} (Optional) + - if no c_ops are given, then this is treated as a lossless propagator and the propagator itself will be of type `AbstractQuantumObject{Operator}`. If c_ops are given, + then the propagator will be of type `AbstractQuantumObject{SuperOperator}` and the save `H` will be the Louivillian. +- tol::Real=1e-6 + Absolute/relative tolerance used by the ODE integrator. This value is + forwarded to the integrator/Propagator to control numerical accuracy. +- solver + The ODE solver algorithm to use (default: `Tsit5()`); any solver object + compatible with the underlying ODE interface is accepted. +- kwargs... + Additional keyword arguments are forwarded to the Propagator constructor + (e.g., options for step control, callback functions, or integration + settings supported by the backend). + +# Notes +- The function does not perform propagation itself; it only constructs and + configures the Propagator object which can then be used to evolve states or + operators. +- For time-dependent `H` provide a `QobjEvo` with an appropriate time + dependence. For time-independent dynamics provide a `QuantumObject`. + +# Returns +- `Propagator` + An initialized Propagator instance set up with + - initial time [0.0], + - initial operator equal to the identity on the Hilbert space of `H` + (created with `one(H)`), + - the requested tolerance and solver, + - dimensions inferred from `H`, and + - any extra options passed via `kwargs`. +""" function propagator( H::Union{QobjEvo, QuantumObject}; tol = 1e-6, solver = Tsit5(), kwargs... ) - return Propagator(H, [0.0], [qeye_like(H)], tol, kwargs, H.dims, solver, Operator()) + return Propagator(H, [0.0], [one(H)], tol, kwargs, H.dims, solver, Operator()) end function propagator( @@ -75,7 +164,7 @@ function propagator( kwargs... ) L = liouvillian(H, c_ops) - return Propagator(L, [0.0], [L], tol, kwargs, L.dims, solver, SuperOperator()) + return Propagator(L, [0.0], [one(L)], tol, kwargs, L.dims, solver, SuperOperator()) end function _lookup_or_compute(p::Propagator, t::Real) @@ -110,21 +199,7 @@ function _compute_propagator(p::Propagator, end -""" - (p::Propagator)(t; t0=0.0, just_interval=false) - -Evaluate the propagator at time t. - -- t: target time. -- t0: initial time. If nonzero: - - when just_interval == true, integrate directly on [t0, t] to obtain U(t, t0). - - otherwise, compute U(t) * inv(U(t0)) using cached values when available. -- just_interval: force a fresh integration on [t0, t] instead of composing from cache. This is useful if both t0 and t are large but t-t0 is small. - -Returns a QuantumObject with the same dims/type as p.H. Results are cached at requested -times (within tolerance p.tol) to avoid recomputation. -""" -function (p::Propagator)(t::Real; t0 = 0.0, just_interval = false) +function (p::Propagator)(t::Real, t0 = 0.0; just_interval = false) if just_interval return _compute_propagator(p, t; t0 = t0) end From d0c2dfebfe4bfdddd4f23f83d6de5f7d6f5a2916 Mon Sep 17 00:00:00 2001 From: Gavin Date: Sat, 6 Sep 2025 21:36:59 +0000 Subject: [PATCH 5/9] added tests in core_test/propagators.jl and docstrings. --- propagator_test.jl | 54 ++++++++++++++++++++++++++------ src/time_evolution/propagator.jl | 18 ++++++----- test/core-test/propagator.jl | 41 ++++++++++++++++++++++++ 3 files changed, 96 insertions(+), 17 deletions(-) create mode 100644 test/core-test/propagator.jl diff --git a/propagator_test.jl b/propagator_test.jl index 486fb6fbc..0023d96b5 100644 --- a/propagator_test.jl +++ b/propagator_test.jl @@ -1,15 +1,51 @@ -using Pkg -Pkg.activate(".") -include("src/QuantumToolbox.jl") -using .QuantumToolbox +using QuantumToolbox using LinearAlgebra using Revise -p_se = propagator(sigmax()) -p_me = propagator(sigmax(), [sigmaz()]) +# p_se = propagator(QobjEvo(sigmax())) +# p_me = propagator(sigmax(), [sigmaz()]) -println(p_se(1.0)) -println("-------------------") -println(p_me(2, 0.1)) \ No newline at end of file +# println(p_se(1.0)) +# println("-------------------") +# println(p_me(2, 0.1)) + +sx = sigmax() +sy = sigmay() +sz = sigmaz() +I = one(sz) + +H = sz + QobjEvo((0.1*sx, (p,t) -> sin(2*t))) + +p_se = propagator(H, progress_bar = false) + +c_ops = [0.1*sz] +p_me = propagator(H, c_ops; progress_bar = false) + +sesolve_res = [] +mesolve_res = [] + +t0 = 0.0 +t = 10.0 +for i in 0:1 + psi = fock(2, i) + push!(sesolve_res, sesolve(H, psi, [t0, t], progress_bar = false)) +end + +sup_psi_base = mat2vec(fock(2,0)*fock(2,0)') +for i in 1:4 + sup_psi = 0 * sup_psi_base + sup_psi[i] = 1 + push!(mesolve_res, mesolve(H, vec2mat(sup_psi), [t0, t], c_ops, progress_bar = false)) +end + +U_se = p_se(t) +U_me = p_me(t) + + +U_sesolve = QuantumObject(hcat([sesolve_res[i].states[end].data for i in 1:2]...)) +U_mesolve = QuantumObject(hcat([mat2vec(mesolve_res[i].states[end]).data for i in 1:4]...), type = SuperOperator()) + +println(norm(U_se-U_sesolve)<1e-8) +println(norm(U_me-U_mesolve)<1e-8) \ No newline at end of file diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl index af2c0ee72..c56b16cfe 100644 --- a/src/time_evolution/propagator.jl +++ b/src/time_evolution/propagator.jl @@ -9,10 +9,10 @@ function propagatorProblem( inplace::Union{Val, Bool} = Val(true), kwargs... ) - H_evo = QobjEvo(H) + H_evo = _sesolve_make_U_QobjEvo(H) + U = H_evo.data - p0 = one(H).data - U = -1im*H_evo.data + p0 = one(H_evo(0)).data kwargs2 = _merge_saveat(tspan, nothing, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...) kwargs3 = _generate_se_me_kwargs(nothing, makeVal(progress_bar), tspan, kwargs2, SaveFuncSESolve) @@ -28,10 +28,10 @@ function propagatorProblem( inplace::Union{Val, Bool} = Val(true), kwargs... ) - H_evo = QobjEvo(H) + L_evo = _mesolve_make_L_QobjEvo(H, []) + U = L_evo.data - p0 = one(H).data - U = H_evo.data + p0 = one(L_evo(0)).data kwargs2 = _merge_saveat(tspan, nothing, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...) kwargs3 = _generate_se_me_kwargs(nothing, makeVal(progress_bar), tspan, kwargs2, SaveFuncMESolve) @@ -153,7 +153,8 @@ function propagator( solver = Tsit5(), kwargs... ) - return Propagator(H, [0.0], [one(H)], tol, kwargs, H.dims, solver, Operator()) + H_evo = QobjEvo(H) + return Propagator(H, [0.0], [one(H_evo(0))], tol, kwargs, H.dims, solver, Operator()) end function propagator( @@ -164,7 +165,8 @@ function propagator( kwargs... ) L = liouvillian(H, c_ops) - return Propagator(L, [0.0], [one(L)], tol, kwargs, L.dims, solver, SuperOperator()) + L_evo = QobjEvo(L) + return Propagator(L, [0.0], [one(L_evo(0))], tol, kwargs, L.dims, solver, SuperOperator()) end function _lookup_or_compute(p::Propagator, t::Real) diff --git a/test/core-test/propagator.jl b/test/core-test/propagator.jl new file mode 100644 index 000000000..940c61927 --- /dev/null +++ b/test/core-test/propagator.jl @@ -0,0 +1,41 @@ +@testitem "Propagators" begin + sx = sigmax() + sy = sigmay() + sz = sigmaz() + I = one(sz) + + H = sz + QobjEvo((0.1*sx, (p,t) -> sin(2*t))) + + p_se = propagator(H) + + c_ops = [0.1*sz] + p_me = propagator(H, c_ops) + + sesolve_res = [] + mesolve_res = [] + + t0 = 0.0 + t = 10.0 + for i in 0:1 + psi = fock(2, i) + push!(sesolve_res, sesolve(H, psi, [t0, t])) + end + + sup_psi_base = mat2vec(fock(2,0)*fock(2,0)') + for i in 1:4 + sup_psi = 0 * sup_psi_base + sup_psi[i] = 1 + push!(mesolve_res, mesolve(H, vec2mat(sup_psi), [t0, t], c_ops)) + end + + U_se = p_se(t) + U_me = p_me(t) + + + U_sesolve = QuantumObject(hcat([sesolve_res[i].states[end].data for i in 1:2]...)) + U_mesolve = QuantumObject(hcat([mat2vec(mesolve_res[i].states[end]).data for i in 1:4]...), type = SuperOperator()) + + @test norm(U_se-U_sesolve) < 1e-8 + @test norm(U_me-U_mesolve) < 1e-8 + +end \ No newline at end of file From 386e0086dfaf405253502d9261860fabf1158845 Mon Sep 17 00:00:00 2001 From: Gavin Date: Sat, 6 Sep 2025 21:40:26 +0000 Subject: [PATCH 6/9] ran formatter --- propagator_test.jl | 8 +-- src/time_evolution/ensemble_sesolve.jl | 90 ++++++++++++++------------ src/time_evolution/propagator.jl | 58 +++++++---------- test/core-test/propagator.jl | 18 +++--- 4 files changed, 82 insertions(+), 92 deletions(-) diff --git a/propagator_test.jl b/propagator_test.jl index 0023d96b5..7702aa94e 100644 --- a/propagator_test.jl +++ b/propagator_test.jl @@ -3,7 +3,6 @@ using LinearAlgebra using Revise - # p_se = propagator(QobjEvo(sigmax())) # p_me = propagator(sigmax(), [sigmaz()]) @@ -16,7 +15,7 @@ sy = sigmay() sz = sigmaz() I = one(sz) -H = sz + QobjEvo((0.1*sx, (p,t) -> sin(2*t))) +H = sz + QobjEvo((0.1*sx, (p, t) -> sin(2*t))) p_se = propagator(H, progress_bar = false) @@ -33,7 +32,7 @@ for i in 0:1 push!(sesolve_res, sesolve(H, psi, [t0, t], progress_bar = false)) end -sup_psi_base = mat2vec(fock(2,0)*fock(2,0)') +sup_psi_base = mat2vec(fock(2, 0)*fock(2, 0)') for i in 1:4 sup_psi = 0 * sup_psi_base sup_psi[i] = 1 @@ -43,9 +42,8 @@ end U_se = p_se(t) U_me = p_me(t) - U_sesolve = QuantumObject(hcat([sesolve_res[i].states[end].data for i in 1:2]...)) U_mesolve = QuantumObject(hcat([mat2vec(mesolve_res[i].states[end]).data for i in 1:4]...), type = SuperOperator()) println(norm(U_se-U_sesolve)<1e-8) -println(norm(U_me-U_mesolve)<1e-8) \ No newline at end of file +println(norm(U_me-U_mesolve)<1e-8) diff --git a/src/time_evolution/ensemble_sesolve.jl b/src/time_evolution/ensemble_sesolve.jl index 7a7a97209..2324f732a 100644 --- a/src/time_evolution/ensemble_sesolve.jl +++ b/src/time_evolution/ensemble_sesolve.jl @@ -5,29 +5,37 @@ export EnsembleTimeEvolutionProblem A structure representing an ensemble time evolution problem for quantum systems. # Fields -- `prob::PT`: The base time evolution problem. -- `func::PF`: A function used to modify or sample parameters for each trajectory in the ensemble. -- `iterate_params::Bool`: If `true`, parameters are iterated for each trajectory; otherwise, the same parameters are used. -- `full_iterator::AbstractArray`: An array containing all parameter sets or states to be used in the ensemble. -- `n_states::Int`: The number of initial states. -- `trajectories::Int`: The total number of trajectories to simulate. + + - `prob::PT`: The base time evolution problem. + - `func::PF`: A function used to modify or sample parameters for each trajectory in the ensemble. + - `iterate_params::Bool`: If `true`, parameters are iterated for each trajectory; otherwise, the same parameters are used. + - `full_iterator::AbstractArray`: An array containing all parameter sets or states to be used in the ensemble. + - `n_states::Int`: The number of initial states. + - `trajectories::Int`: The total number of trajectories to simulate. # Usage -This is used when setting up ensemble sesolve problems, useful for simulating multiple quantum states or parameter sets in parallel. + +This is used when setting up ensemble sesolve problems, useful for simulating multiple quantum states or parameter sets in parallel. Example: + ```julia - H = 2 * π * 0.1 * sigmax() - ψ0 = basis(2, 0) # spin-up - tlist = LinRange(0.0, 100.0, 100) +H = 2 * π * 0.1 * sigmax() +ψ0 = basis(2, 0) # spin-up +tlist = LinRange(0.0, 100.0, 100) - ψs = [ψ0, basis(2, 1)] # spin-up and spin-down +ψs = [ψ0, basis(2, 1)] # spin-up and spin-down - params = collect(Iterators.product([0,1,2,3,4,5], [0,1,2,3,4,5], [0,1,2,3,4,5])) - res = sesolve(H, ψs, tlist; params = params, iterate_params = true, alg = Tsit5(), progress_bar=false); +params = collect(Iterators.product([0, 1, 2, 3, 4, 5], [0, 1, 2, 3, 4, 5], [0, 1, 2, 3, 4, 5])) +res = sesolve(H, ψs, tlist; params = params, iterate_params = true, alg = Tsit5(), progress_bar = false); ``` """ -struct EnsembleTimeEvolutionProblem{PT<:TimeEvolutionProblem,PF<:Function, X<:Vector{T} where T<:QuantumObject{Ket}, Y<:AbstractArray} +struct EnsembleTimeEvolutionProblem{ + PT<:TimeEvolutionProblem, + PF<:Function, + X<:Vector{T} where T<:QuantumObject{Ket}, + Y<:AbstractArray, +} prob::PT func::PF states::X @@ -38,9 +46,8 @@ end function EnsembleTimeEvolutionProblem( prob::PT, states::Vector{T}, - params::AbstractArray = [NullParameters()] - ) where {PT<:TimeEvolutionProblem, T<:QuantumObject{Ket}} - + params::AbstractArray = [NullParameters()], +) where {PT<:TimeEvolutionProblem,T<:QuantumObject{Ket}} problem_dims = (length(states), length(params)) function ensemble_func(prob, i, repeat) @@ -52,27 +59,30 @@ function EnsembleTimeEvolutionProblem( return EnsembleTimeEvolutionProblem(prob, ensemble_func, states, params, problem_dims, trajectories) end +function sesolve( + prob::EnsembleTimeEvolutionProblem, + alg::OrdinaryDiffEqAlgorithm = Tsit5(); + backend = EnsembleThreads(), +) + ensemble_prob = EnsembleProblem(prob.prob.prob, prob_func = prob.func) + sols = solve(ensemble_prob, alg, backend, trajectories = prob.trajectories) -function sesolve(prob::EnsembleTimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5(); backend = EnsembleThreads()) - ensemble_prob = EnsembleProblem(prob.prob.prob, prob_func = prob.func) - sols = solve(ensemble_prob, alg, backend, trajectories = prob.trajectories) - - to_return = Array{TimeEvolutionSol}(undef, prob.problem_dims) + to_return = Array{TimeEvolutionSol}(undef, prob.problem_dims) for i in 1:length(sols) ψt = map(ϕ -> QuantumObject(ϕ, type = Ket(), dims = prob.prob.dimensions), sols[i].u) sol = TimeEvolutionSol( - prob.prob.times, - sols[i].t, - ψt, - _get_expvals(sols[i], SaveFuncSESolve), - sols[i].retcode, - sols[i].alg, - sols[i].prob.kwargs[:abstol], - sols[i].prob.kwargs[:reltol], - ) + prob.prob.times, + sols[i].t, + ψt, + _get_expvals(sols[i], SaveFuncSESolve), + sols[i].retcode, + sols[i].alg, + sols[i].prob.kwargs[:abstol], + sols[i].prob.kwargs[:reltol], + ) to_return[CartesianIndices(to_return)[i]] = sol end - + return to_return end @@ -86,8 +96,8 @@ function sesolve( progress_bar::Union{Val,Bool} = Val(false), inplace::Union{Val,Bool} = Val(true), backend = EnsembleThreads(), - kwargs...,) where T<:QuantumObject{Ket} - + kwargs..., +) where {T<:QuantumObject{Ket}} prob_init = sesolveProblem( H, ψ0s[1], @@ -97,8 +107,8 @@ function sesolve( progress_bar = progress_bar, inplace = inplace, kwargs..., - ) - + ) + trajectories = length(ψ0s) # function ensemble_func(prob, i, repeat) @@ -106,10 +116,6 @@ function sesolve( # end # ensemble_prob = EnsembleTimeEvolutionProblem(prob_init, ensemble_func, ψ0s, trajectories) - ensemble_prob = EnsembleTimeEvolutionProblem( - prob_init, - ψ0s, - [params] - ) + ensemble_prob = EnsembleTimeEvolutionProblem(prob_init, ψ0s, [params]) return sesolve(ensemble_prob, alg; backend = backend) -end \ No newline at end of file +end diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl index c56b16cfe..f55f9be12 100644 --- a/src/time_evolution/propagator.jl +++ b/src/time_evolution/propagator.jl @@ -5,9 +5,9 @@ function propagatorProblem( H::AbstractQuantumObject{Operator}, tspan::AbstractVector; params = NullParameters(), - progress_bar::Union{Val, Bool} = Val(true), - inplace::Union{Val, Bool} = Val(true), - kwargs... + progress_bar::Union{Val,Bool} = Val(true), + inplace::Union{Val,Bool} = Val(true), + kwargs..., ) H_evo = _sesolve_make_U_QobjEvo(H) U = H_evo.data @@ -16,7 +16,7 @@ function propagatorProblem( kwargs2 = _merge_saveat(tspan, nothing, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...) kwargs3 = _generate_se_me_kwargs(nothing, makeVal(progress_bar), tspan, kwargs2, SaveFuncSESolve) - + return ODEProblem{getVal(inplace),FullSpecialize}(U, p0, tspan, params; kwargs3...) end @@ -24,9 +24,9 @@ function propagatorProblem( H::AbstractQuantumObject{SuperOperator}, tspan::AbstractVector; params = NullParameters(), - progress_bar::Union{Val, Bool} = Val(true), - inplace::Union{Val, Bool} = Val(true), - kwargs... + progress_bar::Union{Val,Bool} = Val(true), + inplace::Union{Val,Bool} = Val(true), + kwargs..., ) L_evo = _mesolve_make_L_QobjEvo(H, []) U = L_evo.data @@ -35,7 +35,7 @@ function propagatorProblem( kwargs2 = _merge_saveat(tspan, nothing, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...) kwargs3 = _generate_se_me_kwargs(nothing, makeVal(progress_bar), tspan, kwargs2, SaveFuncMESolve) - + return ODEProblem{getVal(inplace),FullSpecialize}(U, p0, tspan, params; kwargs3...) end @@ -87,15 +87,15 @@ If `U` is a propagator object, to get the propagator between t0 and t, call `U(t - if `just_interval=true` then the propagator is calulated by directly integrating from t0 to t. This approach does't save anything but is useful if t-t0 >> t0 as it avoids the long integration required to get t0. """ -struct Propagator{T<:Union{AbstractQuantumObject{Operator}, AbstractQuantumObject{SuperOperator}}} +struct Propagator{T<:Union{AbstractQuantumObject{Operator},AbstractQuantumObject{SuperOperator}}} H::T times::AbstractArray{Float64} - props::AbstractArray{A} where A <: AbstractQuantumObject + props::AbstractArray{A} where {A<:AbstractQuantumObject} tol::Float64 solver_kwargs::Base.Pairs - dims::Union{AbstractArray, Tuple} + dims::Union{AbstractArray,Tuple} solver::OrdinaryDiffEqAlgorithm - type + type::Any end function Base.show(io::IO, p::Propagator) @@ -103,7 +103,7 @@ function Base.show(io::IO, p::Propagator) println(" dims: $(p.dims)") println(" issuper: $(issuper(p.H))") println(" size: $(size(p.H))") - println(" Number of Saved Times: $(length(p.times))") + return println(" Number of Saved Times: $(length(p.times))") end @doc raw""" @@ -147,23 +147,18 @@ given Hamiltonian or time-dependent operator. - dimensions inferred from `H`, and - any extra options passed via `kwargs`. """ -function propagator( - H::Union{QobjEvo, QuantumObject}; - tol = 1e-6, - solver = Tsit5(), - kwargs... - ) +function propagator(H::Union{QobjEvo,QuantumObject}; tol = 1e-6, solver = Tsit5(), kwargs...) H_evo = QobjEvo(H) return Propagator(H, [0.0], [one(H_evo(0))], tol, kwargs, H.dims, solver, Operator()) end function propagator( - H::Union{QobjEvo, QuantumObject}, - c_ops :: Union{AbstractArray, Tuple}; + H::Union{QobjEvo,QuantumObject}, + c_ops::Union{AbstractArray,Tuple}; tol = 1e-6, solver = Tsit5(), - kwargs... - ) + kwargs..., +) L = liouvillian(H, c_ops) L_evo = QobjEvo(L) return Propagator(L, [0.0], [one(L_evo(0))], tol, kwargs, L.dims, solver, SuperOperator()) @@ -178,10 +173,10 @@ function _lookup_or_compute(p::Propagator, t::Real) # Check if we found an exact match or close enough within tolerance if (idx.start <= length(p.times)) && (abs(t - p.times[idx.start]) <= p.tol) U = p.props[idx.start] - elseif (idx.start > 1) && (abs(t - p.times[idx.start - 1]) <= p.tol) - U = p.props[idx.start - 1] + elseif (idx.start > 1) && (abs(t - p.times[idx.start-1]) <= p.tol) + U = p.props[idx.start-1] else - t0 = (idx.start > 1) ? p.times[idx.start - 1] : 0.0 + t0 = (idx.start > 1) ? p.times[idx.start-1] : 0.0 U = _compute_propagator(p, t) insert!(p.props, idx.start, U) insert!(p.times, idx.start, t) @@ -190,17 +185,13 @@ function _lookup_or_compute(p::Propagator, t::Real) return U end -function _compute_propagator(p::Propagator, - t::Real; - t0 = 0.0 - ) +function _compute_propagator(p::Propagator, t::Real; t0 = 0.0) prob = propagatorProblem(p.H, [t0, t]; p.solver_kwargs...) res = solve(prob, p.solver) U = QuantumObject(res.u[end], dims = p.dims, type = p.type) return U end - function (p::Propagator)(t::Real, t0 = 0.0; just_interval = false) if just_interval return _compute_propagator(p, t; t0 = t0) @@ -218,7 +209,4 @@ function (p::Propagator)(t::Real, t0 = 0.0; just_interval = false) return U end -function Base.size(p::Propagator) - return size(p.H) -end - +Base.size(p::Propagator) = size(p.H) diff --git a/test/core-test/propagator.jl b/test/core-test/propagator.jl index 940c61927..5cc6cf944 100644 --- a/test/core-test/propagator.jl +++ b/test/core-test/propagator.jl @@ -4,7 +4,7 @@ sz = sigmaz() I = one(sz) - H = sz + QobjEvo((0.1*sx, (p,t) -> sin(2*t))) + H = sz + QobjEvo((0.1*sx, (p, t) -> sin(2*t))) p_se = propagator(H) @@ -17,25 +17,23 @@ t0 = 0.0 t = 10.0 for i in 0:1 - psi = fock(2, i) - push!(sesolve_res, sesolve(H, psi, [t0, t])) + psi = fock(2, i) + push!(sesolve_res, sesolve(H, psi, [t0, t])) end - sup_psi_base = mat2vec(fock(2,0)*fock(2,0)') + sup_psi_base = mat2vec(fock(2, 0)*fock(2, 0)') for i in 1:4 - sup_psi = 0 * sup_psi_base - sup_psi[i] = 1 - push!(mesolve_res, mesolve(H, vec2mat(sup_psi), [t0, t], c_ops)) + sup_psi = 0 * sup_psi_base + sup_psi[i] = 1 + push!(mesolve_res, mesolve(H, vec2mat(sup_psi), [t0, t], c_ops)) end U_se = p_se(t) U_me = p_me(t) - U_sesolve = QuantumObject(hcat([sesolve_res[i].states[end].data for i in 1:2]...)) U_mesolve = QuantumObject(hcat([mat2vec(mesolve_res[i].states[end]).data for i in 1:4]...), type = SuperOperator()) @test norm(U_se-U_sesolve) < 1e-8 @test norm(U_me-U_mesolve) < 1e-8 - -end \ No newline at end of file +end From 2a361eb352073ab2a08e10dfcfd187db58570a21 Mon Sep 17 00:00:00 2001 From: Gavin Date: Wed, 10 Sep 2025 17:19:15 +0000 Subject: [PATCH 7/9] fixing .gitignore --- .gitignore | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 0d3f129f9..014c79889 100644 --- a/.gitignore +++ b/.gitignore @@ -11,4 +11,7 @@ benchmarks/benchmarks_output.json .ipynb_checkpoints *.ipynb -.devcontainer \ No newline at end of file +.devcontainer +package-lock.json +package.json +/node_modules \ No newline at end of file From a002306a3af2da73d088afbe4fe8bed30b6bd81c Mon Sep 17 00:00:00 2001 From: Gavin Date: Wed, 10 Sep 2025 18:46:04 +0000 Subject: [PATCH 8/9] cleaned up extra files and Project.toml --- Project.toml | 4 --- propagator_test.jl | 49 ------------------------------------ test/core-test/propagator.jl | 9 +++++++ 3 files changed, 9 insertions(+), 53 deletions(-) delete mode 100644 propagator_test.jl diff --git a/Project.toml b/Project.toml index 9b460d1d2..5fe3b1b63 100644 --- a/Project.toml +++ b/Project.toml @@ -18,9 +18,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -ProtoStructs = "437b6fc4-8e8e-11e9-3fa1-ad391e66c018" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -62,9 +60,7 @@ Makie = "0.24" OrdinaryDiffEqCore = "1" OrdinaryDiffEqTsit5 = "1" Pkg = "1" -ProtoStructs = "1.2.1" Random = "1" -Revise = "3.9.0" SciMLBase = "2.105" SciMLOperators = "1.4" SparseArrays = "1" diff --git a/propagator_test.jl b/propagator_test.jl deleted file mode 100644 index 7702aa94e..000000000 --- a/propagator_test.jl +++ /dev/null @@ -1,49 +0,0 @@ -using QuantumToolbox -using LinearAlgebra - -using Revise - -# p_se = propagator(QobjEvo(sigmax())) -# p_me = propagator(sigmax(), [sigmaz()]) - -# println(p_se(1.0)) -# println("-------------------") -# println(p_me(2, 0.1)) - -sx = sigmax() -sy = sigmay() -sz = sigmaz() -I = one(sz) - -H = sz + QobjEvo((0.1*sx, (p, t) -> sin(2*t))) - -p_se = propagator(H, progress_bar = false) - -c_ops = [0.1*sz] -p_me = propagator(H, c_ops; progress_bar = false) - -sesolve_res = [] -mesolve_res = [] - -t0 = 0.0 -t = 10.0 -for i in 0:1 - psi = fock(2, i) - push!(sesolve_res, sesolve(H, psi, [t0, t], progress_bar = false)) -end - -sup_psi_base = mat2vec(fock(2, 0)*fock(2, 0)') -for i in 1:4 - sup_psi = 0 * sup_psi_base - sup_psi[i] = 1 - push!(mesolve_res, mesolve(H, vec2mat(sup_psi), [t0, t], c_ops, progress_bar = false)) -end - -U_se = p_se(t) -U_me = p_me(t) - -U_sesolve = QuantumObject(hcat([sesolve_res[i].states[end].data for i in 1:2]...)) -U_mesolve = QuantumObject(hcat([mat2vec(mesolve_res[i].states[end]).data for i in 1:4]...), type = SuperOperator()) - -println(norm(U_se-U_sesolve)<1e-8) -println(norm(U_me-U_mesolve)<1e-8) diff --git a/test/core-test/propagator.jl b/test/core-test/propagator.jl index 5cc6cf944..1f07b7ba8 100644 --- a/test/core-test/propagator.jl +++ b/test/core-test/propagator.jl @@ -36,4 +36,13 @@ @test norm(U_se-U_sesolve) < 1e-8 @test norm(U_me-U_mesolve) < 1e-8 + + U_se = p_se(t, 3.0) + U_se_JI = p_se(t, 3.0; just_interval = true) + U_me = p_me(t, 3.0) + U_me_JI = p_me(t, 3.0; just_interval = true) + + @test norm(U_se-U_se_JI)<1e-6 + @test norm(U_me-U_me_JI) < 1e-6 + end From 40eda7e365fee97dd257cb241558eb95725026f9 Mon Sep 17 00:00:00 2001 From: Gavin Date: Wed, 10 Sep 2025 18:47:58 +0000 Subject: [PATCH 9/9] missed a few files. --- other_missing_things.md | 48 ---------- src/QuantumToolbox.jl | 1 - src/time_evolution/ensemble_sesolve.jl | 121 ------------------------- 3 files changed, 170 deletions(-) delete mode 100644 other_missing_things.md delete mode 100644 src/time_evolution/ensemble_sesolve.jl diff --git a/other_missing_things.md b/other_missing_things.md deleted file mode 100644 index cf334ef06..000000000 --- a/other_missing_things.md +++ /dev/null @@ -1,48 +0,0 @@ -# Things that are missing from the QuantumToolbox.jl package. - -1) Powers are not supported for `QobjEvo` objects. - - There is a bug with the one I did, as the following returns a stackoverflow error. - ```julia - import Base.^ - import Base.^ - function ^(x::QobjEvo, y::Integer) - if y< 0 - error("Exponent must be non-negative integer") - end - if y == 0 - return QobjEvo(qeye(size(x)[1], dims = x.dims)) - end - to_return = x - for i in 2:y - to_return = to_return * x - end - return to_return - end - - 1*(sum([(QobjEvo((qeye(2), (p,t) -> 1)))^n*1 for n in 0:1:2])) - 1*(sum([1*(QobjEvo((qeye(2), (p,t) -> 1)))^n for n in 0:1:2])) - ``` - - However, the following do not - ```julia - 1*(sum([(QobjEvo((qeye(2), (p,t) -> 1)))^n for n in 0:1:2])) - 1*(sum([(1*QobjEvo((qeye(2), (p,t) -> 1)))^n for n in 0:1:2])) - 1*(sum([(QobjEvo((qeye(2), (p,t) -> 1))*1)^n for n in 0:1:2])) - ``` - - The error comes from expressions like - ```julia - a = QobjEvo((qeye(2), (p,t) -> 1)) - 1*(a+a*a*1) - ``` - - The error is: - StackOverflowError: - ```julia - StackOverflowError: - - Stacktrace: - [1] *(α::ScalarOperator{Int64, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, x::UniformScaling{Int64}) - @ SciMLOperators ~/.julia/packages/SciMLOperators/DKrpP/src/scalar.jl:366 - [2] *(α::ScalarOperator{Int64, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, x::UniformScaling{Int64}) - @ SciMLOperators ~/.julia/packages/SciMLOperators/DKrpP/src/scalar.jl:369 - ``` \ No newline at end of file diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl index d23c3242e..73165da06 100644 --- a/src/QuantumToolbox.jl +++ b/src/QuantumToolbox.jl @@ -116,7 +116,6 @@ include("time_evolution/mcsolve.jl") include("time_evolution/ssesolve.jl") include("time_evolution/smesolve.jl") include("time_evolution/time_evolution_dynamical.jl") -include("time_evolution/ensemble_sesolve.jl") include("time_evolution/propagator.jl") # Others diff --git a/src/time_evolution/ensemble_sesolve.jl b/src/time_evolution/ensemble_sesolve.jl deleted file mode 100644 index 2324f732a..000000000 --- a/src/time_evolution/ensemble_sesolve.jl +++ /dev/null @@ -1,121 +0,0 @@ -export EnsembleTimeEvolutionProblem -""" - EnsembleTimeEvolutionProblem{PT<:TimeEvolutionProblem, PF<:Function} - -A structure representing an ensemble time evolution problem for quantum systems. - -# Fields - - - `prob::PT`: The base time evolution problem. - - `func::PF`: A function used to modify or sample parameters for each trajectory in the ensemble. - - `iterate_params::Bool`: If `true`, parameters are iterated for each trajectory; otherwise, the same parameters are used. - - `full_iterator::AbstractArray`: An array containing all parameter sets or states to be used in the ensemble. - - `n_states::Int`: The number of initial states. - - `trajectories::Int`: The total number of trajectories to simulate. - -# Usage - -This is used when setting up ensemble sesolve problems, useful for simulating multiple quantum states or parameter sets in parallel. - -Example: - -```julia -H = 2 * π * 0.1 * sigmax() -ψ0 = basis(2, 0) # spin-up -tlist = LinRange(0.0, 100.0, 100) - -ψs = [ψ0, basis(2, 1)] # spin-up and spin-down - -params = collect(Iterators.product([0, 1, 2, 3, 4, 5], [0, 1, 2, 3, 4, 5], [0, 1, 2, 3, 4, 5])) -res = sesolve(H, ψs, tlist; params = params, iterate_params = true, alg = Tsit5(), progress_bar = false); -``` -""" -struct EnsembleTimeEvolutionProblem{ - PT<:TimeEvolutionProblem, - PF<:Function, - X<:Vector{T} where T<:QuantumObject{Ket}, - Y<:AbstractArray, -} - prob::PT - func::PF - states::X - params::Y - problem_dims::Tuple - trajectories::Int -end -function EnsembleTimeEvolutionProblem( - prob::PT, - states::Vector{T}, - params::AbstractArray = [NullParameters()], -) where {PT<:TimeEvolutionProblem,T<:QuantumObject{Ket}} - problem_dims = (length(states), length(params)) - - function ensemble_func(prob, i, repeat) - state_id = mod1(i, problem_dims[1]) - param_id = div(i - 1, problem_dims[1]) + 1 - return remake(prob, u0 = states[state_id].data, p = params[param_id]) - end - trajectories = prod(problem_dims) - return EnsembleTimeEvolutionProblem(prob, ensemble_func, states, params, problem_dims, trajectories) -end - -function sesolve( - prob::EnsembleTimeEvolutionProblem, - alg::OrdinaryDiffEqAlgorithm = Tsit5(); - backend = EnsembleThreads(), -) - ensemble_prob = EnsembleProblem(prob.prob.prob, prob_func = prob.func) - sols = solve(ensemble_prob, alg, backend, trajectories = prob.trajectories) - - to_return = Array{TimeEvolutionSol}(undef, prob.problem_dims) - for i in 1:length(sols) - ψt = map(ϕ -> QuantumObject(ϕ, type = Ket(), dims = prob.prob.dimensions), sols[i].u) - sol = TimeEvolutionSol( - prob.prob.times, - sols[i].t, - ψt, - _get_expvals(sols[i], SaveFuncSESolve), - sols[i].retcode, - sols[i].alg, - sols[i].prob.kwargs[:abstol], - sols[i].prob.kwargs[:reltol], - ) - to_return[CartesianIndices(to_return)[i]] = sol - end - - return to_return -end - -function sesolve( - H::Union{AbstractQuantumObject{Operator},Tuple}, - ψ0s::Vector{T}, - tlist::AbstractVector; - alg::OrdinaryDiffEqAlgorithm = Tsit5(), - e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, - params = NullParameters(), - progress_bar::Union{Val,Bool} = Val(false), - inplace::Union{Val,Bool} = Val(true), - backend = EnsembleThreads(), - kwargs..., -) where {T<:QuantumObject{Ket}} - prob_init = sesolveProblem( - H, - ψ0s[1], - tlist; - e_ops = e_ops, - params = params, - progress_bar = progress_bar, - inplace = inplace, - kwargs..., - ) - - trajectories = length(ψ0s) - - # function ensemble_func(prob, i, repeat) - # return remake(prob, u0 = ψ0s[i].data) - # end - - # ensemble_prob = EnsembleTimeEvolutionProblem(prob_init, ensemble_func, ψ0s, trajectories) - ensemble_prob = EnsembleTimeEvolutionProblem(prob_init, ψ0s, [params]) - return sesolve(ensemble_prob, alg; backend = backend) -end