From 2f9600c75c48ec9e269eec8a1c45bc089cf263d4 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Thu, 23 Jul 2026 11:18:06 +0200 Subject: [PATCH 1/7] build: move to SQA v0.10 --- Project.toml | 4 +- docs/Project.toml | 3 +- ..._bidirectional-waveguide_coherent-pulse.jl | 9 ++-- examples/08-1_pulse-delay__simple.jl | 11 ++-- ...ing-cat-state__PRA2020_102-023717_fig3b.jl | 3 +- ...pagation-delay__PRA2026_113-013730_fig4.jl | 10 ++-- src/pulses.jl | 23 +++++--- test/test_translate.jl | 53 ++++++++++--------- 8 files changed, 69 insertions(+), 47 deletions(-) diff --git a/Project.toml b/Project.toml index 22b38344..e0b55941 100644 --- a/Project.toml +++ b/Project.toml @@ -29,10 +29,10 @@ LinearAlgebra = "1.10" ModelingToolkitBase = "1.42.2" NumericalIntegration = "0.3" OrdinaryDiffEq = "6" -QuantumCumulants = "0.5" +QuantumCumulants = "0.7" QuantumOptics = "1" QuantumOpticsBase = "0.5" -SecondQuantizedAlgebra = "0.9" +SecondQuantizedAlgebra = "0.10" SpecialFunctions = "2" StaticArrays = "1" SymbolicUtils = "4" diff --git a/docs/Project.toml b/docs/Project.toml index 99b50570..b9f60894 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -24,8 +24,9 @@ MacroTools = "0.5" ModelingToolkitBase = "1" OrdinaryDiffEq = "6" Plots = "1" -QuantumCumulants = "0.5" +QuantumCumulants = "0.7" QuantumOptics = "1" QuantumOpticsBase = "0.5" +SecondQuantizedAlgebra = "0.10" SymbolicUtils = "4" Symbolics = "7" diff --git a/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl b/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl index ab91ea51..ad7f9f4c 100644 --- a/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl +++ b/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl @@ -8,7 +8,7 @@ using QuantumInputOutput using SecondQuantizedAlgebra using Symbolics: Symbolics using QuantumOptics -using QuantumOpticsBase: dagger +using QuantumOpticsBase: dagger, static_operator using Plots using LaTeXStrings @@ -159,9 +159,12 @@ lT = length(T) G2 = zeros(lT, lT) # transmission G2_ref = zeros(lT, lT) # reflection -L0(t) = L_R_QO(t) +# `L_R_QO(t)`/`L_L_QO(t)` return a lazy `TimeDependentSum`; materialize it to a concrete +# operator at each time so the quantum-regression products below yield a plain operator +# (a `TimeDependentSum` product cannot serve as the solver's initial state). +L0(t) = dense(static_operator(L_R_QO(t))) L0_dag(t) = dagger(L0(t)) -L0_ref(t) = L_L_QO(t) +L0_ref(t) = dense(static_operator(L_L_QO(t))) L0_ref_dag(t) = dagger(L0_ref(t)) for it1 = 1:(lT-1) diff --git a/examples/08-1_pulse-delay__simple.jl b/examples/08-1_pulse-delay__simple.jl index 78a1a592..16715be0 100644 --- a/examples/08-1_pulse-delay__simple.jl +++ b/examples/08-1_pulse-delay__simple.jl @@ -190,6 +190,9 @@ p_t_sym = [g_u, g_in, g_out, g_v, M_ls...] p_t_num = [gu_, gin_, gout_, gv_, M_t_ls...] dict_p_t_int = Dict(p_t_sym .=> p_t_num) +# The interaction picture eliminates the delay cavity `d`, so the numeric operators live +# on the two-mode basis `bu ⊗ bv` (with `ad` mapped to the identity). Build on that basis. +b_int = bu ⊗ bv au_int = destroy(bu) ⊗ one(bv) ## This is amazing! # hide ad_int = one(bu ⊗ bv) @@ -198,13 +201,13 @@ operators = Dict( [au, au', ad, ad', av, av'] .=> [au_int, au_int', ad_int, ad_int', av_int, av_int'], ) -H_int_QO = to_numeric(H_int_sym, b; time_parameter = dict_p_t_int, operators) +H_int_QO = to_numeric(H_int_sym, b_int; time_parameter = dict_p_t_int, operators) L_int_QO = [ - to_numeric(L_int_sym[i], b; time_parameter = dict_p_t_int, operators) for + to_numeric(L_int_sym[i], b_int; time_parameter = dict_p_t_int, operators) for i = 1:length(L_int_sym) ] -## H_int_QO = to_numeric(H_int_sym, b; time_parameter=dict_p_t_int) # hide -## L_int_QO = [to_numeric(L_int_sym[i], b; time_parameter=dict_p_t_int) for i=1:length(L_int_sym)] # hide +## H_int_QO = to_numeric(H_int_sym, b_int; time_parameter=dict_p_t_int) # hide +## L_int_QO = [to_numeric(L_int_sym[i], b_int; time_parameter=dict_p_t_int) for i=1:length(L_int_sym)] # hide nothing # hide # diff --git a/examples/drafts/02-2_traveling-cat-state__PRA2020_102-023717_fig3b.jl b/examples/drafts/02-2_traveling-cat-state__PRA2020_102-023717_fig3b.jl index 6d43b192..f16f70b2 100644 --- a/examples/drafts/02-2_traveling-cat-state__PRA2020_102-023717_fig3b.jl +++ b/examples/drafts/02-2_traveling-cat-state__PRA2020_102-023717_fig3b.jl @@ -13,6 +13,7 @@ using QuantumInputOutput using SecondQuantizedAlgebra using QuantumOptics +using QuantumOpticsBase: dagger using Plots using LinearAlgebra using DataInterpolations @@ -23,7 +24,7 @@ using DataInterpolations hc = FockSpace(:c) ## symbolic operator -a = Destroy(hc, :a, 1) +a = Destroy(hc, :a) ## symbolic parameters @variables γ::Real K::Real Δ::Real p::Number diff --git a/examples/drafts/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl b/examples/drafts/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl index 38ddc011..40dcca0f 100644 --- a/examples/drafts/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl +++ b/examples/drafts/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl @@ -11,6 +11,7 @@ using QuantumInputOutput using SecondQuantizedAlgebra using Symbolics: Symbolics using QuantumOptics +using QuantumOpticsBase: dagger using SymbolicUtils using LinearAlgebra using Plots @@ -181,9 +182,9 @@ t_, ρt = timeevolution.master_dynamic(T, ψ0_fock, input_output) expect(au_QO'au_QO, ρt) -p = plot(T, expect(au_QO' * au_QO, ρt); label = "") -plot!(p, T, expect(ad1_QO' * ad1_QO, ρt); label = "") -plot!(p, T, expect(ad2_QO' * ad2_QO, ρt); label = "") +p = plot(T, real.(expect(au_QO' * au_QO, ρt)); label = "") +plot!(p, T, real.(expect(ad1_QO' * ad1_QO, ρt)); label = "") +plot!(p, T, real.(expect(ad2_QO' * ad2_QO, ρt)); label = "") p @@ -197,8 +198,6 @@ expect(σ22_QO, ρt[end]) pop_fock = zeros(length(Δ_ls)) pop_coh = zeros(length(Δ_ls)) idx_t1 = findmin(abs.(T .- t1))[2] -using ProgressMeter -prog = Progress(length(Δ_ls)) for (it, Δn) in enumerate(Δ_ls) Δn = Δ_ls[it] @@ -219,7 +218,6 @@ for (it, Δn) in enumerate(Δ_ls) # t_, ρt = timeevolution.master_dynamic(T, ψ0_coh, input_output) # pop_coh[it] = real(expect(σ22_QO, ρt[idx_t1])) - next!(prog) end nothing # hide diff --git a/src/pulses.jl b/src/pulses.jl index 0b9b616c..60c3e92d 100644 --- a/src/pulses.jl +++ b/src/pulses.jl @@ -6,6 +6,15 @@ const _tol_div = 1e-10 const _extrapolate = ExtrapolationType.Extension const _ϵ = 1e-10 +# Wrap a sampled interpolation in an explicit single-argument closure `t -> g(t)`. +# SecondQuantizedAlgebra v0.10 rejects `time_parameter` values whose callable has +# ambiguous arity (a `LinearInterpolation` responds to both `g(t)` and `g(t, order)`), +# so the coupling functions hand back an unambiguous `t -> value` closure that plugs +# straight into `to_numeric`'s `time_parameter` while still being callable as `g(t)`. +_as_time_fn(itp) = let g = itp + t -> g(t) +end + _mode_interp(mode::AbstractVector, T::AbstractVector) = LinearInterpolation(mode, T; extrapolation = _extrapolate) _mode_interp(mode, T::AbstractVector) = mode @@ -54,7 +63,7 @@ function _compute_coupling(mode::Vector, T::Vector, denom_fn) g[i] = mode[i]' / sqrt(d) end end - return LinearInterpolation(g, T; extrapolation = _extrapolate) + return _as_time_fn(LinearInterpolation(g, T; extrapolation = _extrapolate)) end # ────────────────────────────────────────────── @@ -65,7 +74,8 @@ end coupling_input(u, T) Compute the virtual-cavity input coupling ``g_u(t)`` from an input mode `u(t)` -sampled on time grid `T`. Returns the `LinearInterpolation` directly (callable as `g(t)`). +sampled on time grid `T`. Returns a single-argument closure `g(t)` (an interpolation of +the sampled coupling) that plugs directly into a `to_numeric` `time_parameter`. """ coupling_input(u::Vector, T::Vector) = _compute_coupling(u, T, x -> abs(1 - x) + _ϵ) coupling_input(u::Function, T::Vector) = coupling_input(u.(T), T) @@ -85,7 +95,8 @@ end coupling_output(v, T) Compute the virtual-cavity output coupling ``g_v(t)`` from an output mode `v(t)` -sampled on time grid `T`. Returns the `LinearInterpolation` directly (callable as `g(t)`). +sampled on time grid `T`. Returns a single-argument closure `g(t)` (an interpolation of +the sampled coupling) that plugs directly into a `to_numeric` `time_parameter`. """ coupling_output(v::Vector, T::Vector) = _compute_coupling(-v, T, x -> x + _ϵ) coupling_output(v::Function, T::Vector) = coupling_output(v.(T), T) @@ -290,14 +301,14 @@ function _compute_coupling_delay(num_mode::Vector, u::Vector, v::Vector, T::Vect g[i] = num_mode[i]' / sqrt(d + _ϵ) end end - return LinearInterpolation(g, T; extrapolation = _extrapolate) + return _as_time_fn(LinearInterpolation(g, T; extrapolation = _extrapolate)) end """ coupling_delay_out(u, v, T) Compute the out-coupling strength for a delay cavity. -Returns `LinearInterpolation` directly. +Returns a single-argument closure `g(t)`. """ coupling_delay_out(u::Vector, v::Vector, T::Vector) = _compute_coupling_delay(u, u, v, T) coupling_delay_out(u::Function, v::Function, T::Vector) = @@ -309,7 +320,7 @@ coupling_delay_out(u::LinearInterpolation, v::LinearInterpolation, T::Vector) = coupling_delay_in(u, v, T) Compute the in-coupling strength for a delay cavity. -Returns `LinearInterpolation` directly. +Returns a single-argument closure `g(t)`. """ coupling_delay_in(u::Vector, v::Vector, T::Vector) = _compute_coupling_delay(-v, u, v, T) coupling_delay_in(u::Function, v::Function, T::Vector) = coupling_delay_in(u.(T), v.(T), T) diff --git a/test/test_translate.jl b/test/test_translate.jl index 6b16ddcf..51e7dda6 100644 --- a/test/test_translate.jl +++ b/test/test_translate.jl @@ -2,9 +2,15 @@ using QuantumInputOutput using SecondQuantizedAlgebra using SymbolicUtils using QuantumOpticsBase -using QuantumOpticsBase: dagger +using QuantumOpticsBase: dagger, static_operator, TimeDependentSum using Test +# Since SecondQuantizedAlgebra v0.10, `to_numeric` returns a `TimeDependentSum` for +# time-dependent inputs instead of a plain callable. Calling it at a time `t` sets the +# current time and still returns a `TimeDependentSum`; `mat` materializes it to a dense +# `Operator` for numeric comparison. +mat(F, t) = dense(static_operator(F(t))) + @testset "translate" begin @variables κ_L::Real κ_R::Real Δ::Real g::Real γ::Real @variables E::Complex E1::Complex E2::Complex @@ -56,8 +62,8 @@ using Test time_parameter = dict_p_t2, operators = ops_dict, ) - @test isa(F1, Function) - @test isequal(F1(0.1), a_QO*3) + @test F1 isa TimeDependentSum + @test sum(abs.((mat(F1, 0.1) - dense(a_QO*3)).data)) < 1e-12 F2 = to_numeric( a*E, bc1; @@ -65,7 +71,7 @@ using Test time_parameter = dict_p_t2, operators = ops_dict, ) - @test isequal(F2(0.1), a_QO*E_t(0.1)) + @test sum(abs.((mat(F2, 0.1) - dense(a_QO*E_t(0.1))).data)) < 1e-12 # association of the scalar prefactor differs from `Δn*A'*A`, compare numerically @test sum( abs.( @@ -101,38 +107,37 @@ using Test abs.((dense(to_numeric(Δ, b; parameter = dict_p1)) - dense(one(b)*Δn)).data), ) < 1e-12 F3 = to_numeric(Δ, b; parameter = dict_p1, time_parameter = dict_p_t2) - @test sum(abs.((dense(F3(4)) - dense(one(b)*Δn)).data)) < 1e-8 + @test sum(abs.((mat(F3, 4) - dense(one(b)*Δn)).data)) < 1e-8 F4 = to_numeric( a*3*conj(E) + Δ*σ(2, 2), b; parameter = dict_p1, time_parameter = dict_p_t2, ) - @test sum(abs.((F4(0.2) - dense(a_QO2*3*E_t_c(0.2) + Δn*σ_QO(2, 2))).data)) < 1e-8 + @test sum(abs.((mat(F4, 0.2) - dense(a_QO2*3*E_t_c(0.2) + Δn*σ_QO(2, 2))).data)) < 1e-8 F5 = to_numeric( a*conj(E) + Δ*σ(2, 2), b; parameter = dict_p1, time_parameter = dict_p_t2, ) - @test sum(abs.((F5(0.2) - dense(a_QO2*E_t_c(0.2) + Δn*σ_QO(2, 2))).data)) < 1e-8 - @inferred F5(0.2) # QAdd time-dependent path should return concrete type + @test sum(abs.((mat(F5, 0.2) - dense(a_QO2*E_t_c(0.2) + Δn*σ_QO(2, 2))).data)) < 1e-8 F5_ = to_numeric(a*conj(E), b; parameter = dict_p1, time_parameter = dict_p_t2) - @test sum(abs.((dense(F5_(0.2)) - dense(a_QO2*E_t_c(0.2))).data)) < 1e-8 + @test sum(abs.((mat(F5_, 0.2) - dense(a_QO2*E_t_c(0.2))).data)) < 1e-8 F6 = to_numeric(conj(E), b; parameter = dict_p1, time_parameter = dict_p_t2) - @test sum(abs.((dense(F6(0.2)) - dense(E_t_c(0.2)*one(b))).data)) < 1e-8 + @test sum(abs.((mat(F6, 0.2) - dense(E_t_c(0.2)*one(b))).data)) < 1e-8 F7 = to_numeric(conj(E) + Δ*σ(2, 2), b; parameter = dict_p1, time_parameter = dict_p_t2) - @test sum(abs.((F7(0.2) - dense(E_t_c(0.2)*one(b) + Δn*σ_QO(2, 2))).data)) < 1e-8 + @test sum(abs.((mat(F7, 0.2) - dense(E_t_c(0.2)*one(b) + Δn*σ_QO(2, 2))).data)) < 1e-8 # a bare symbolic scalar without a numeric/time value cannot be translated @test_throws ArgumentError to_numeric(conj(E), b; parameter = dict_p1) F8 = to_numeric(E^2, b; parameter = dict_p1, time_parameter = dict_p_t2) - @test sum(abs.((dense(F8(0.2)) - dense(E_t(0.2)^2*one(b))).data)) < 1e-8 + @test sum(abs.((mat(F8, 0.2) - dense(E_t(0.2)^2*one(b))).data)) < 1e-8 @testset "time_parameter_normalization" begin dict_p_t_num = Dict([E] .=> [2.5]) F_num = to_numeric(a*E, b; parameter = dict_p1, time_parameter = dict_p_t_num) - @test sum(abs.((dense(F_num(0.2)) - dense(a_QO2 * 2.5)).data)) < 1e-8 + @test sum(abs.((mat(F_num, 0.2) - dense(a_QO2 * 2.5)).data)) < 1e-8 end @testset "multiple_time_prefactors" begin @@ -147,7 +152,7 @@ using Test time_parameter = dict_p_t_multi, ) expected = dense(a_QO2 * E1_t(0.4) * E2_t_c(0.4)) - @test sum(abs.((dense(F_multi(0.4)) - expected).data)) < 1e-8 + @test sum(abs.((mat(F_multi, 0.4) - expected).data)) < 1e-8 end @testset "to_numeric vector overload" begin @@ -160,10 +165,10 @@ using Test operators = ops_dict, ) @test length(translated_ops) == length(ops) - @test all(op -> op isa Function, translated_ops) - @test isequal(translated_ops[1](0.3), a_QO) - @test isequal(translated_ops[2](0.3), one(bc1) * Δn) - @test isequal(translated_ops[3](0.3), a_QO * E_t(0.3)) + @test all(op -> op isa TimeDependentSum, translated_ops) + @test sum(abs.((mat(translated_ops[1], 0.3) - dense(a_QO)).data)) < 1e-12 + @test sum(abs.((mat(translated_ops[2], 0.3) - dense(one(bc1) * Δn)).data)) < 1e-12 + @test sum(abs.((mat(translated_ops[3], 0.3) - dense(a_QO * E_t(0.3))).data)) < 1e-12 end @testset "to_numeric SLH overload" begin @@ -176,11 +181,11 @@ using Test operators = ops_dict, ) - @test H_QO isa Function + @test H_QO isa TimeDependentSum @test length(L_QO) == 1 - @test L_QO[1] isa Function - @test sum(abs.((H_QO(0.4) - dense(Δn * dagger(a_QO) * a_QO)).data)) < 1e-8 - @test sum(abs.((L_QO[1](0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) < 1e-8 + @test L_QO[1] isa TimeDependentSum + @test sum(abs.((mat(H_QO, 0.4) - dense(Δn * dagger(a_QO) * a_QO)).data)) < 1e-8 + @test sum(abs.((mat(L_QO[1], 0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) < 1e-8 end @testset "substitute operators qmul" begin @@ -231,7 +236,7 @@ using Test @variables gR::Real gR_t(t) = 1.0 + 2.0im F = to_numeric(im * gR * a3, b3; time_parameter = Dict(gR => gR_t)) - @test F isa Function - @test sum(abs.((F(0.0) - dense((im * (1.0 + 2.0im)) * a3_QO)).data)) < 1e-8 + @test F isa TimeDependentSum + @test sum(abs.((mat(F, 0.0) - dense((im * (1.0 + 2.0im)) * a3_QO)).data)) < 1e-8 end end From 718e7ea631ed2b5fd67464461779f32d58fe6327 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Thu, 23 Jul 2026 12:11:21 +0200 Subject: [PATCH 2/7] format --- src/pulses.jl | 7 ++++--- test/test_translate.jl | 3 ++- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/pulses.jl b/src/pulses.jl index 60c3e92d..bb99b0ae 100644 --- a/src/pulses.jl +++ b/src/pulses.jl @@ -11,9 +11,10 @@ const _ϵ = 1e-10 # ambiguous arity (a `LinearInterpolation` responds to both `g(t)` and `g(t, order)`), # so the coupling functions hand back an unambiguous `t -> value` closure that plugs # straight into `to_numeric`'s `time_parameter` while still being callable as `g(t)`. -_as_time_fn(itp) = let g = itp - t -> g(t) -end +_as_time_fn(itp) = + let g = itp + t -> g(t) + end _mode_interp(mode::AbstractVector, T::AbstractVector) = LinearInterpolation(mode, T; extrapolation = _extrapolate) diff --git a/test/test_translate.jl b/test/test_translate.jl index 51e7dda6..471499db 100644 --- a/test/test_translate.jl +++ b/test/test_translate.jl @@ -185,7 +185,8 @@ mat(F, t) = dense(static_operator(F(t))) @test length(L_QO) == 1 @test L_QO[1] isa TimeDependentSum @test sum(abs.((mat(H_QO, 0.4) - dense(Δn * dagger(a_QO) * a_QO)).data)) < 1e-8 - @test sum(abs.((mat(L_QO[1], 0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) < 1e-8 + @test sum(abs.((mat(L_QO[1], 0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) < + 1e-8 end @testset "substitute operators qmul" begin From 2630ae5310a3c21e9391663ad7742845d48fb45e Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Thu, 23 Jul 2026 12:28:25 +0200 Subject: [PATCH 3/7] fix perf --- benchmarks/correlations.jl | 14 +++-- ...cattering__PRL2019_123-123604_fig2-fig3.jl | 5 +- src/correlations.jl | 51 ++++++++++++++----- 3 files changed, 47 insertions(+), 23 deletions(-) diff --git a/benchmarks/correlations.jl b/benchmarks/correlations.jl index 3e31a3bf..8a6c5979 100644 --- a/benchmarks/correlations.jl +++ b/benchmarks/correlations.jl @@ -42,24 +42,22 @@ function benchmark_correlations!(SUITE) H_QO = to_numeric(H_sym, b; parameter = dict_p, time_parameter = dict_p_t) L_QO = to_numeric(L_sym, b; parameter = dict_p, time_parameter = dict_p_t) - function input_output(t, ρ) - Ht = H_QO(t) - J = [L_QO(t)] - return Ht, J, QuantumOpticsBase.dagger.(J) - end - ψ0 = fockstate(bu1, 1) ⊗ fockstate(bc1, 0) ⊗ fockstate(bv1, 0) - _, ρt = timeevolution.master_dynamic(T, ψ0, input_output) + _, ρt = timeevolution.master_dynamic(T, ψ0, H_QO, [L_QO]) au_qo = to_numeric(au, b) c_qo = to_numeric(c, b) Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo ## --- Two-time correlation --- + # Pass the time-dependent operators (`H_QO`, `[L_QO]`) directly to the solver rather + # than wrapping them in a `(t, ρ)` closure: since v0.10 `to_numeric` returns a lazy + # `TimeDependentSum`, the direct path lets the solver build its integrator once and is + # markedly faster than re-reading the operator from a function at every step. SUITE["Correlations"]["two-time"] = BenchmarkGroup() SUITE["Correlations"]["two-time"]["single photon cavity"] = - @benchmarkable correlation_matrix($T, $ρt, $input_output, $Ls) + @benchmarkable correlation_matrix($T, $ρt, $H_QO, [$L_QO], $Ls) return nothing end diff --git a/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl b/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl index 5edf9b4b..85b32172 100644 --- a/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl +++ b/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl @@ -110,7 +110,10 @@ nothing # hide # In order to determine suitable temporal output modes we calculate the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle$ and diagonalize the matrix to obtain the eigenvalues with the corresponding eigenvectors. The eigenvalues correspond to the mean photon number $n_i$ in the corresponding temporal eigenvector mode $v_i$. Ls(t) = gu_t(t)*au_qo + √(γ_)*c_qo -g1_m = correlation_matrix(T, ρt, input_output_1, Ls) +## Pass the time-dependent operators (`H_QO`, `[L_QO]`) straight to the solver instead of +## the `input_output_1` closure: the direct path lets the solver reuse its integrator and +## is much faster for the many independent solves of the quantum-regression loop. +g1_m = correlation_matrix(T, ρt, H_QO, [L_QO], Ls) nothing # hide # diff --git a/src/correlations.jl b/src/correlations.jl index 177dbc58..627d6908 100644 --- a/src/correlations.jl +++ b/src/correlations.jl @@ -5,35 +5,58 @@ Compute the two-time correlation matrix ``g^{(1)}(t_1, t_2) = \\langle L_s^\\dagger(t_1) L_s(t_2) \\rangle`` on the time grid `T`. Writes directly into output matrix. + +The dynamics can be supplied either as a `master_dynamic`-style function `f(t, ρ)` or, +much more efficiently for time-dependent problems, as operators passed straight to the +solver: a time-dependent `H` (an `AbstractTimeDependentOperator`, e.g. the +`TimeDependentSum` returned by [`to_numeric`](@ref)) together with its jump operators `J`, +or a constant `H` with constant `J`. Passing the operators directly lets the solver build +the integrator once instead of rebuilding it from `f`'s return value at every step. + +`Ls` may be a constant operator or a function `Ls(t)` returning the (concrete) operator at +time `t`. """ -function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls::Function; kwargs...) - Ls_vec = Ls.(T) - Ls_dag_vec = dagger.(Ls_vec) +function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls; kwargs...) + Ls_vec, Ls_dag_vec = _ls_vectors(T, Ls) _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 timeevolution.master_dynamic(T_slice, ρ0, f; kwargs...) end end -function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls; kwargs...) - Ls_dag = dagger(Ls) - l_T = length(T) - Ls_vec = fill(Ls, l_T) - Ls_dag_vec = fill(Ls_dag, l_T) +function correlation_matrix( + T::Vector, + ρt::Vector, + H::QuantumOpticsBase.AbstractTimeDependentOperator, + J::AbstractVector, + Ls; + kwargs..., +) + Ls_vec, Ls_dag_vec = _ls_vectors(T, Ls) + # `master_dynamic` mutates the operator's current time via `set_time!`, so each parallel + # iteration solves with its own copy to avoid a data race on the shared `H`/`J`. _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 - timeevolution.master_dynamic(T_slice, ρ0, f; kwargs...) + timeevolution.master_dynamic(T_slice, ρ0, copy(H), [copy(j) for j in J]; kwargs...) end end -function correlation_matrix(T::Vector, ρt::Vector, H, J::Vector, Ls; kwargs...) - Ls_dag = dagger(Ls) - l_T = length(T) - Ls_vec = fill(Ls, l_T) - Ls_dag_vec = fill(Ls_dag, l_T) +function correlation_matrix(T::Vector, ρt::Vector, H, J::AbstractVector, Ls; kwargs...) + Ls_vec, Ls_dag_vec = _ls_vectors(T, Ls) _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 timeevolution.master(T_slice, ρ0, H, J; kwargs...) end end +# Build the per-time `Ls`/`Ls†` vectors. A function `Ls(t)` is sampled on `T`; a constant +# operator is broadcast to every time. +function _ls_vectors(T::Vector, Ls::Function) + Ls_vec = Ls.(T) + return Ls_vec, dagger.(Ls_vec) +end +function _ls_vectors(T::Vector, Ls) + l_T = length(T) + return fill(Ls, l_T), fill(dagger(Ls), l_T) +end + function _correlation_loop(solve_fn, T, ρt, Ls_vec, Ls_dag_vec) l_T = length(T) @assert l_T == length(ρt) From 7f7dd975bcb57ff34292f9aa61767ce135c1c66d Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Thu, 23 Jul 2026 16:08:21 +0200 Subject: [PATCH 4/7] PulseCoupling --- docs/src/api.md | 4 + ...26-07-23-qio-v0.10-numeric-layer-design.md | 181 ++++++++++++++++++ ...cattering__PRL2019_123-123604_fig2-fig3.jl | 33 ++-- ..._bidirectional-waveguide_coherent-pulse.jl | 5 +- src/QuantumInputOutput.jl | 1 + src/correlations.jl | 47 ++--- src/pulses.jl | 124 ++++++------ test/test_code_quality.jl | 1 + test/test_translate.jl | 26 ++- test/test_utils.jl | 8 + 10 files changed, 322 insertions(+), 108 deletions(-) create mode 100644 docs/superpowers/specs/2026-07-23-qio-v0.10-numeric-layer-design.md diff --git a/docs/src/api.md b/docs/src/api.md index 85650a1a..d6475011 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -50,6 +50,10 @@ to_numeric Gaussian ``` +```@docs +PulseCoupling +``` + ```@docs coupling_input ``` diff --git a/docs/superpowers/specs/2026-07-23-qio-v0.10-numeric-layer-design.md b/docs/superpowers/specs/2026-07-23-qio-v0.10-numeric-layer-design.md new file mode 100644 index 00000000..063fa0c2 --- /dev/null +++ b/docs/superpowers/specs/2026-07-23-qio-v0.10-numeric-layer-design.md @@ -0,0 +1,181 @@ +# QIO v0.10-native numeric layer + +Date: 2026-07-23 +Status: Draft (awaiting review) + +## Context + +QuantumInputOutput.jl (QIO) was written against SecondQuantizedAlgebra (SQA) v0.9, whose +numeric model was: a time-dependent numeric operator **is a callable** `t -> operator`, and a +`time_parameter` value **is any callable** `t -> value`. Two QIO subsystems were shaped by +that model and now read as mismatches against SQA v0.10. + +1. **Couplings return a `LinearInterpolation`.** SQA v0.10 introduced parameter-aware + coefficients (`(p, t) -> value`) and now **introspects the arity** of every + `time_parameter` value (`_tp_reads_p` in `SQA/src/numeric/coeff.jl`): it requires a single, + unambiguous arity (1 for `t`, 2 for `(p, t)`) and rejects anything else. A + `DataInterpolations.LinearInterpolation`'s call operator carries methods of more than one + arity (a 1-arg scalar/vector form plus 2-arg vector forms; note there is no `itp(t, deriv)` + scalar method), so `_tp_reads_p` cannot resolve a single arity and rejects it with + `ArgumentError: ... has methods of conflicting arity [1, 2]; write it as an explicit + t -> value or (p, t) -> value closure`. QIO papered over this with `_as_time_fn` (an + anonymous `t -> itp(t)` wrapper, exactly the closure SQA's error prescribes). The wrapper is + v0.10-correct, but it signals that QIO leaks an interpolation *object* where SQA's currency + is a *time function*. + +2. **SLH composition uses `FunctionWrapper`/callables.** `SLH.jl` carries a callable-composition + layer (`_Callable = Union{Function, FunctionWrapper}`, the `_add`/`_mul`/`_adj` closure + branches, `_detect_operator_type`, `_maybe_wrap_lindblad`, `_maybe_wrap_hamiltonian`, + `_fw_return_type`, `_op_type`, `_to_func`). It exists to compose SLH triples whose numeric + time-dependent entries are closures `t -> operator`, and to pin their return type via + `FunctionWrapper{OpType, Tuple{Float64}}` for type stability. This looked like a v0.9 relic, + but measurement (see "Investigation: is `FunctionWrapper` still needed?" below) shows it + provides concreteness that native `TimeDependentSum` composition does not replace. It is + therefore **kept** in this PR and revisited in the future backend-agnostic PR. + +### How SQA v0.10 consumes a `time_parameter` value (verified by reading the pipeline) + +`api.jl` calls `core.jl:_to_numeric_td`, which calls `coeff.jl:_time_basis`/`_td_coeff`. + +- `_time_basis` collects the value functions into a **`Tuple`** and runs `_tp_reads_p` on each. +- `_td_coeff` compiles the symbolic prefactor with `build_function` and returns the per-term + coefficient closure `t -> pref(map(f -> f(t), valuefuncs)...)`. +- The QuantumOptics backend stores these as `TimeDependentSum(ComplexF64, b, b, coeffs, ops)`. + +SQA's contract is therefore exact: **a `time_parameter` value must be a single-arity, concrete +callable returning a `Number`.** For `map(f -> f(t), valuefuncs)` to stay type-stable, each +value function must be a concrete callable type. + +## Goals + +Make QIO's numeric layer speak SQA v0.10 natively, while staying representation-agnostic so a +later backend-neutral PR needs no rework here. + +## Non-goals (explicitly deferred to a future PR) + +- Moving `QuantumOptics`/`QuantumOpticsBase` to weakdeps behind package extensions. +- QuantumToolbox support, or any second numeric backend. +- Any backend-hook abstraction for the solver/correlation layer. + +This PR stays **QuantumOptics-only**, but every new type and rewrite here is **backend-agnostic +by construction** (it names no backend-specific numeric type), so the future PR inherits it +unchanged. + +## Deliverable 1: `PulseCoupling`, a first-class time-function type + +### The type + +```julia +""" + PulseCoupling(f) + +A single-argument time function `g(t)` used as a `to_numeric` `time_parameter` value. +Wraps any evaluator `f` (a sampled interpolation, an analytic closure, ...) behind exactly +one call method, so it satisfies SQA's single-arity `time_parameter` contract by +construction. Callable as `g(t)` and broadcastable as `g.(T)`. +""" +struct PulseCoupling{F} + f::F +end +(g::PulseCoupling)(t::Real) = g.f(t) +``` + +- **Single-arity by construction**, so it satisfies `_tp_reads_p` with no wrapper. + `_as_time_fn` is deleted. +- **Concrete/static**: `PulseCoupling{LinearInterpolation{...}}` is concrete; it lands in SQA's + `valuefuncs` `Tuple` concretely, keeping `map(f -> f(t), ...)` inferrable, and it returns a + concrete `ComplexF64`. +- **Performant**: one direct call into `f`, no extra allocation, identical hot-path cost to the + raw interpolation minus the arity ambiguity. +- **Backend-agnostic**: a scalar time function that names no QuantumOptics/QuantumToolbox type. +- **Public/introspectable**: documented and exported like `Gaussian`, with `g.f` reachable. + +### Routing + +Every coupling constructor returns a `PulseCoupling`. + +- `coupling_input(u::Vector, T)`, `coupling_output`, `coupling_delay_in`, `coupling_delay_out`: + wrap the computed `LinearInterpolation` in `PulseCoupling` (replacing `_as_time_fn`). +- `coupling_input(::Gaussian)`, `coupling_output(::Gaussian)`: wrap the analytic closure in + `PulseCoupling` too, for a uniform return type. These were already single-arity, so this is + for consistency, not correctness. +- Internal uses that relied on the return being callable or broadcastable (for example + `coupling_output(v_, T).(T)` and the `effective_*_mode` collections in `pulses.jl`) keep + working, because `PulseCoupling` is callable and `g.(T)` broadcasts. + +`_as_time_fn` is removed. The existing `coupling_*(u::LinearInterpolation, T)` input-mode +dispatch methods stay, since they accept an interpolation as the *input mode*, unrelated to the +return type. + +### Interaction picture (light touch) + +`interaction_picture.jl` builds matrix-valued time functions (`M(t)`), and example 08-1 extracts +scalar elements with `t -> M_t(t)[i, j]` before putting them in `time_parameter`. Those +per-element accessors are the same "scalar time function" currency, so they can be wrapped as +`PulseCoupling(t -> M_t(t)[i, j])` at the point they enter `time_parameter`. The matrix-valued +API itself does not change; this is an example-level consistency update. + +## Investigation: is `FunctionWrapper` still needed? (yes) + +The `FunctionWrapper` layer looked like a v0.9 relic to delete in favour of native +`TimeDependentSum` arithmetic. Direct measurement (QuantumOptics backend, `FockBasis`) shows +otherwise: + +| Probe | Result | +|:--|:--| +| Two SQA `TimeDependentSum` on the same basis share one concrete type | yes | +| `SVector(tds, tds)` eltype concrete | yes | +| `SVector(static_op, tds)` eltype | `AbstractOperator` (**not** concrete) | +| A single `tds`'s internal containers (coeffs, `LazySum` ops) | abstract (`Vector{Function}`, `Vector{AbstractOperator}`) even though the outer wrapper type is concrete | +| `tds + tds`, `scalar * tds`, `adjoint(tds)` stay `<: TimeDependentSum` | yes, but each returns a **differently-parameterized** concrete type (`typeof(x) != typeof(tds)`) | +| `tds * tds` preserves the **exact** concrete type | yes | +| `SLH(tds) ▷ SLH(tds)` (cascade) | **errors** (`MethodError: zero(::TimeDependentSum)`) | +| `iszero(tds)` | **errors** (`MethodError`) | + +Conclusions: + +- `FunctionWrapper` provided a **single uniform concrete type** for mixed and composed numeric + SLH entries. Native `TimeDependentSum` composition does not: a mixed static/time-dependent + Lindblad vector degrades to `SVector{N, AbstractOperator}`; the outer `TimeDependentSum` type + is concrete but wraps abstract containers; and `+`/`scalar *`/`adjoint` each return a + differently-parameterized `TimeDependentSum` (only `*` preserves the exact concrete type), so a + heterogeneous collection of composed entries still widens to an abstract eltype. +- Restoring that concreteness the v0.10 way would mean promoting every entry to a + `TimeDependentSum`. But **QuantumToolbox has no `TimeDependentSum`** (it uses `QobjEvo`), so a + promotion step is backend-specific and must live in the future backend-agnostic PR, alongside + the per-backend abstraction. +- Separately, numeric SLH composition on `TimeDependentSum` entries currently **errors** (the + `iszero`-on-operator bug), and no example exercises time-dependent numeric SLH composition, so + nothing in the current workflows is blocked. + +**Decision:** keep `FunctionWrapper` unchanged in this PR. Its removal, the concreteness-restoring +promotion, and the `iszero`-on-operator fix are all entangled with the per-backend abstraction and +are deferred to the future backend-agnostic PR. + +## Testing + +- **`PulseCoupling`**: unit tests that a coupling `isa PulseCoupling`, is callable, broadcasts, + returns `ComplexF64`, and is accepted by `to_numeric`'s `time_parameter` without error (the + case that previously required `_as_time_fn`). Keep the existing pulse-coupling numeric + comparisons. +- Full `Pkg.test()` green, and all examples still run headless. The `SLH.jl` and `test_SLH.jl` + `FunctionWrapper` tests are untouched. + +## Migration and back-compat + +- Examples using `coupling_*` are unaffected (still callable and broadcastable). +- Example 08-1 interaction-picture element accessors optionally wrapped in `PulseCoupling`. +- No breaking changes to SLH construction or composition. + +## Risks + +- **`PulseCoupling` broadcast/dispatch**: confirm `g.(T)` still broadcasts and that the internal + `effective_*_mode` collections (which call and broadcast couplings) are unaffected. Low risk, + covered by the existing pulse tests. + +## Scope note for review + +This PR is scoped to Deliverable 1 (`PulseCoupling`) only. The `FunctionWrapper` question is +resolved above: it stays, and the SLH-composition modernization is deferred to the future +backend-agnostic PR where it can be done correctly (per-backend promotion, not +`TimeDependentSum`-specific code). diff --git a/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl b/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl index 85b32172..982c46a6 100644 --- a/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl +++ b/examples/01-1_cavity-scattering__PRL2019_123-123604_fig2-fig3.jl @@ -1,8 +1,8 @@ # # Cavity Scattering of a Single Photon -# In this example, we simulate the scattering of a resonant single photon on an empty one-sided cavity. The temporal mode of the light pulse is a Gaussian with width $\sigma$ and the cavity has a decay rate $\gamma$. This system has been studied in [A. Kiilerich, et al., Phys. Rev. Lett. 123, 123604 (2019)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.123604). +# In this example, we simulate the scattering of a resonant single photon on an empty one-sided cavity. The temporal mode of the light pulse is a Gaussian with width $\sigma$ and the cavity has a decay rate $\gamma$. This system has been studied in [A. Kiilerich, et al., Phys. Rev. Lett. 123, 123604 (2019)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.123604). -# We start by loading the needed packages and specifying the model. +# We start by loading the needed packages and specifying the model. using QuantumInputOutput using SecondQuantizedAlgebra @@ -29,9 +29,9 @@ av = Destroy(h, :a_v, 3) @variables g_u::Number Δ::Real γ::Real g_v::Number nothing # hide -# We use the symbolic operators and parameters to define the SLH triples and cascade them to obtain the Hamiltonian and Lindblad for the system. +# We use the symbolic operators and parameters to define the SLH triples and cascade them to obtain the Hamiltonian and Lindblad for the system. -G_u = SLH(1, g_u'*au, 0) # input cavity +G_u = SLH(1, g_u'*au, 0) # input cavity G_c = SLH(1, √(γ)*c, Δ*c'c) # system cavity G_v = SLH(1, g_v'*av, 0) # output cavity @@ -46,9 +46,9 @@ H = hamiltonian(G_cas) L = lindblad(G_cas)[1] # only one Lindblad term in this example -# To solve the dynamics of the system we translate the symbolic expressions into numeric operators (matrices) of QuantumOptics.jl. To do so, we define the numerical parameters and operator basis. +# To solve the dynamics of the system we translate the symbolic expressions into numeric operators (matrices) of QuantumOptics.jl. To do so, we define the numerical parameters and operator basis. -## numerical parameters +## numerical parameters γ_ = 1.0 Δ_ = 0.0 @@ -67,21 +67,21 @@ T = [0:0.002:1;]*T_end gu_t = coupling_input(u1, T) dict_p_t = Dict(g_u => gu_t) -## numeric bases +## numeric bases bu1 = FockBasis(1) bc1 = FockBasis(1) bv1 = FockBasis(1) b = bu1 ⊗ bc1 ⊗ bv1 nothing # hide -# We now use the function [`to_numeric`](@ref) to create the numeric operators. If the kwarg `time_parameter` is provided the created operator is a time-dependent function. +# We now use the function [`to_numeric`](@ref) to create the numeric operators. If the kwarg `time_parameter` is provided the created operator is a time-dependent function. H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) L_QO = to_numeric(L, b; parameter = dict_p, time_parameter = dict_p_t) nothing # hide -# To solve the dynamics we use the QuantumOptics.jl function `timeevolution.master_dynamic`. +# To solve the dynamics we use the QuantumOptics.jl function `timeevolution.master_dynamic`. ## time-dependent function for timeevolution.master_dynamic that returns H(t), J(t) and Jd(t) function input_output_1(t, ρ) @@ -107,12 +107,9 @@ n_c_t = real.(expect(c_qo'*c_qo, ρt)) n_u1_t = real.(expect(au_qo'*au_qo, ρt)) nothing # hide -# In order to determine suitable temporal output modes we calculate the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle$ and diagonalize the matrix to obtain the eigenvalues with the corresponding eigenvectors. The eigenvalues correspond to the mean photon number $n_i$ in the corresponding temporal eigenvector mode $v_i$. +# In order to determine suitable temporal output modes we calculate the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle$ and diagonalize the matrix to obtain the eigenvalues with the corresponding eigenvectors. The eigenvalues correspond to the mean photon number $n_i$ in the corresponding temporal eigenvector mode $v_i$. Ls(t) = gu_t(t)*au_qo + √(γ_)*c_qo -## Pass the time-dependent operators (`H_QO`, `[L_QO]`) straight to the solver instead of -## the `input_output_1` closure: the direct path lets the solver reuse its integrator and -## is much faster for the many independent solves of the quantum-regression loop. g1_m = correlation_matrix(T, ρt, H_QO, [L_QO], Ls) nothing # hide @@ -130,7 +127,7 @@ p = heatmap( ) p -# The eigenvalues and corresponding eigenvectors are sorted in ascending order, which means the last eigenvalue corresponds to the highest populated temporal mode. +# The eigenvalues and corresponding eigenvectors are sorted in ascending order, which means the last eigenvalue corresponds to the highest populated temporal mode. F = eigen(g1_m) n_avg = round.(real.(F.values)*ΔT; digits = 3) @@ -140,7 +137,7 @@ v_mode = (modes[:, end]) / sqrt(ΔT) @show n_avg[(end-1):end] nothing # hide -# We use now the mode with the highest mean photon number as our out-mode to determine its quantum state. +# We use now the mode with the highest mean photon number as our out-mode to determine its quantum state. p_sym_2 = [γ, Δ] p_num_2 = [γ_, Δ_] dict_p_2 = Dict(p_sym_2 .=> p_num_2) @@ -201,7 +198,7 @@ plot(p1, p2; layout = (2, 1), size = (600, 550)) # ## Cavity with phase noise -# We slightly adapt the above example by assuming the initial pulse to be in a coherent state and adding phase noise to the cavity. This results in scattering into multiple modes. +# We slightly adapt the above example by assuming the initial pulse to be in a coherent state and adding phase noise to the cavity. This results in scattering into multiple modes. ## new basis of the system bu1_3 = FockBasis(12) @@ -215,7 +212,7 @@ c_3 = embed(b_3, 2, destroy(bc1_3)) av_3 = embed(b_3, 3, destroy(bv1_3)) cdc_3 = c_3'c_3 -## we use the same Hamiltonian as before but add a depasing term to the dissipation +## we use the same Hamiltonian as before but add a depasing term to the dissipation H_QO_3 = to_numeric(H, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2) L_QO_3 = to_numeric(L, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2) function input_output_3(t, ρ) @@ -226,7 +223,7 @@ end; # -# Due to the larger Hilbert space the time evolution takes a few seconds. +# Due to the larger Hilbert space the time evolution takes a few seconds. ψ0_3 = coherentstate(bu1_3, 2) ⊗ fockstate(bc1_3, 0) ⊗ fockstate(bv1_3, 0) t_3, ρt_3 = timeevolution.master_dynamic(T, ψ0_3, input_output_3) diff --git a/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl b/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl index ad7f9f4c..910fc103 100644 --- a/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl +++ b/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl @@ -159,9 +159,8 @@ lT = length(T) G2 = zeros(lT, lT) # transmission G2_ref = zeros(lT, lT) # reflection -# `L_R_QO(t)`/`L_L_QO(t)` return a lazy `TimeDependentSum`; materialize it to a concrete -# operator at each time so the quantum-regression products below yield a plain operator -# (a `TimeDependentSum` product cannot serve as the solver's initial state). +# Materialize the lazy `TimeDependentSum` to a concrete operator at each time, so the +# quantum-regression products below give a plain operator usable as the solver's initial state. L0(t) = dense(static_operator(L_R_QO(t))) L0_dag(t) = dagger(L0(t)) L0_ref(t) = dense(static_operator(L_L_QO(t))) diff --git a/src/QuantumInputOutput.jl b/src/QuantumInputOutput.jl index 600b1c03..3f6e7b05 100644 --- a/src/QuantumInputOutput.jl +++ b/src/QuantumInputOutput.jl @@ -30,6 +30,7 @@ export SLH, to_numeric, substitute, # Pulse coupling + PulseCoupling, coupling_input, coupling_output, coupling_delay_in, diff --git a/src/correlations.jl b/src/correlations.jl index 627d6908..3ea313a7 100644 --- a/src/correlations.jl +++ b/src/correlations.jl @@ -6,18 +6,14 @@ Compute the two-time correlation matrix ``g^{(1)}(t_1, t_2) = \\langle L_s^\\dagger(t_1) L_s(t_2) \\rangle`` on the time grid `T`. Writes directly into output matrix. -The dynamics can be supplied either as a `master_dynamic`-style function `f(t, ρ)` or, -much more efficiently for time-dependent problems, as operators passed straight to the -solver: a time-dependent `H` (an `AbstractTimeDependentOperator`, e.g. the -`TimeDependentSum` returned by [`to_numeric`](@ref)) together with its jump operators `J`, -or a constant `H` with constant `J`. Passing the operators directly lets the solver build -the integrator once instead of rebuilding it from `f`'s return value at every step. - -`Ls` may be a constant operator or a function `Ls(t)` returning the (concrete) operator at -time `t`. +Supply the dynamics either as a `master_dynamic`-style function `f(t, ρ)`, or as operators +passed straight to the solver: a time-dependent `H` (e.g. the `TimeDependentSum` from +[`to_numeric`](@ref)) with jump operators `J`, or a constant `H` with constant `J`. The +operator form is much faster for time-dependent problems (the integrator is built once). +`Ls` is either a constant operator or a function `Ls(t)` returning the operator at `t`. """ function correlation_matrix(T::Vector, ρt::Vector, f::Function, Ls; kwargs...) - Ls_vec, Ls_dag_vec = _ls_vectors(T, Ls) + Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls) _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 timeevolution.master_dynamic(T_slice, ρ0, f; kwargs...) end @@ -31,30 +27,37 @@ function correlation_matrix( Ls; kwargs..., ) - Ls_vec, Ls_dag_vec = _ls_vectors(T, Ls) - # `master_dynamic` mutates the operator's current time via `set_time!`, so each parallel - # iteration solves with its own copy to avoid a data race on the shared `H`/`J`. + Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls) _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 timeevolution.master_dynamic(T_slice, ρ0, copy(H), [copy(j) for j in J]; kwargs...) end end function correlation_matrix(T::Vector, ρt::Vector, H, J::AbstractVector, Ls; kwargs...) - Ls_vec, Ls_dag_vec = _ls_vectors(T, Ls) + Ls_vec, Ls_dag_vec = _sample_operator_and_adjoint(T, Ls) _correlation_loop(T, ρt, Ls_vec, Ls_dag_vec) do T_slice, ρ0 timeevolution.master(T_slice, ρ0, H, J; kwargs...) end end -# Build the per-time `Ls`/`Ls†` vectors. A function `Ls(t)` is sampled on `T`; a constant -# operator is broadcast to every time. -function _ls_vectors(T::Vector, Ls::Function) - Ls_vec = Ls.(T) - return Ls_vec, dagger.(Ls_vec) +function _sample_operator_and_adjoint(T::Vector, op::Function) + vals = op.(T) + return vals, dagger.(vals) end -function _ls_vectors(T::Vector, Ls) - l_T = length(T) - return fill(Ls, l_T), fill(dagger(Ls), l_T) +function _sample_operator_and_adjoint( + ::Vector, + op::QuantumOpticsBase.AbstractTimeDependentOperator, +) + throw( + ArgumentError( + "`Ls` is a time-dependent operator ($(nameof(typeof(op)))); pass it as a " * + "function `Ls(t)` returning the concrete operator at time `t`.", + ), + ) +end +function _sample_operator_and_adjoint(T::Vector, op) + nt = length(T) + return fill(op, nt), fill(dagger(op), nt) end function _correlation_loop(solve_fn, T, ρt, Ls_vec, Ls_dag_vec) diff --git a/src/pulses.jl b/src/pulses.jl index bb99b0ae..c5efc871 100644 --- a/src/pulses.jl +++ b/src/pulses.jl @@ -6,19 +6,23 @@ const _tol_div = 1e-10 const _extrapolate = ExtrapolationType.Extension const _ϵ = 1e-10 -# Wrap a sampled interpolation in an explicit single-argument closure `t -> g(t)`. -# SecondQuantizedAlgebra v0.10 rejects `time_parameter` values whose callable has -# ambiguous arity (a `LinearInterpolation` responds to both `g(t)` and `g(t, order)`), -# so the coupling functions hand back an unambiguous `t -> value` closure that plugs -# straight into `to_numeric`'s `time_parameter` while still being callable as `g(t)`. -_as_time_fn(itp) = - let g = itp - t -> g(t) - end +""" + PulseCoupling(f) + +A single-argument time function `g(t)` used as a `to_numeric` `time_parameter` value. +Wraps any evaluator `f` (a sampled interpolation, an analytic closure, ...) behind exactly +one call method, so it satisfies SecondQuantizedAlgebra's single-arity `time_parameter` +contract by construction (a raw interpolation is rejected for having methods of conflicting +arity). Callable as `g(t)` and broadcastable as `g.(T)`. +""" +struct PulseCoupling{F} + f::F +end +(g::PulseCoupling)(t::Real) = g.f(t) -_mode_interp(mode::AbstractVector, T::AbstractVector) = +_interpolate_mode(mode::AbstractVector, T::AbstractVector) = LinearInterpolation(mode, T; extrapolation = _extrapolate) -_mode_interp(mode, T::AbstractVector) = mode +_interpolate_mode(mode, T::AbstractVector) = mode # ────────────────────────────────────────────── # Gaussian pulse type @@ -52,19 +56,19 @@ end # Shared coupling core # ────────────────────────────────────────────── -function _compute_coupling(mode::Vector, T::Vector, denom_fn) - l_T = length(T) - buf = Vector{Float64}(undef, l_T) - map!(abs2, buf, mode) - ∫m2 = cumul_integrate(T, buf) - g = zeros(ComplexF64, l_T) - @inbounds for i = 2:l_T - d = denom_fn(∫m2[i]) - if sqrt(abs(d)) > _tol_div - g[i] = mode[i]' / sqrt(d) +function _coupling_from_mode(mode::Vector, T::Vector, denominator) + nt = length(T) + mode_sq = Vector{Float64}(undef, nt) + map!(abs2, mode_sq, mode) + ∫mode2 = cumul_integrate(T, mode_sq) + coupling = zeros(ComplexF64, nt) + @inbounds for i = 2:nt + denom = denominator(∫mode2[i]) + if sqrt(abs(denom)) > _tol_div + coupling[i] = mode[i]' / sqrt(denom) end end - return _as_time_fn(LinearInterpolation(g, T; extrapolation = _extrapolate)) + return PulseCoupling(LinearInterpolation(coupling, T; extrapolation = _extrapolate)) end # ────────────────────────────────────────────── @@ -75,10 +79,10 @@ end coupling_input(u, T) Compute the virtual-cavity input coupling ``g_u(t)`` from an input mode `u(t)` -sampled on time grid `T`. Returns a single-argument closure `g(t)` (an interpolation of -the sampled coupling) that plugs directly into a `to_numeric` `time_parameter`. +sampled on time grid `T`. Returns a [`PulseCoupling`](@ref) (callable as `g(t)`, +broadcastable as `g.(T)`) that plugs directly into a `to_numeric` `time_parameter`. """ -coupling_input(u::Vector, T::Vector) = _compute_coupling(u, T, x -> abs(1 - x) + _ϵ) +coupling_input(u::Vector, T::Vector) = _coupling_from_mode(u, T, x -> abs(1 - x) + _ϵ) coupling_input(u::Function, T::Vector) = coupling_input(u.(T), T) coupling_input(u::LinearInterpolation, T::Vector) = coupling_input(u.(T), T) @@ -89,17 +93,17 @@ Analytical input coupling ``g_u(t)`` for a Gaussian pulse. """ function coupling_input(g::Gaussian) mode, ∫m2 = _gaussian_mode(g) - return t -> mode(t)' / √(abs(1 - ∫m2(t)) + _ϵ) + return PulseCoupling(t -> mode(t)' / √(abs(1 - ∫m2(t)) + _ϵ)) end """ coupling_output(v, T) Compute the virtual-cavity output coupling ``g_v(t)`` from an output mode `v(t)` -sampled on time grid `T`. Returns a single-argument closure `g(t)` (an interpolation of -the sampled coupling) that plugs directly into a `to_numeric` `time_parameter`. +sampled on time grid `T`. Returns a [`PulseCoupling`](@ref) (callable as `g(t)`, +broadcastable as `g.(T)`) that plugs directly into a `to_numeric` `time_parameter`. """ -coupling_output(v::Vector, T::Vector) = _compute_coupling(-v, T, x -> x + _ϵ) +coupling_output(v::Vector, T::Vector) = _coupling_from_mode(-v, T, x -> x + _ϵ) coupling_output(v::Function, T::Vector) = coupling_output(v.(T), T) coupling_output(v::LinearInterpolation, T::Vector) = coupling_output(v.(T), T) @@ -110,7 +114,7 @@ Analytical output coupling ``g_v(t)`` for a Gaussian pulse. """ function coupling_output(g::Gaussian) mode, ∫m2 = _gaussian_mode(g) - return t -> -mode(t)' / √(∫m2(t) + _ϵ) + return PulseCoupling(t -> -mode(t)' / √(∫m2(t) + _ϵ)) end # ────────────────────────────────────────────── @@ -130,14 +134,13 @@ All kwargs are passed on to the ODE solver. function effective_output_mode(v_fcts, gv_fcts, T, i; alg = Tsit5(), kwargs...) @assert i > 1 n = i - 1 - # Capture as tuples for concrete closure types - gv_t = ntuple(k -> gv_fcts[k], n) + # Capture as a tuple for concrete closure types + gv = ntuple(k -> gv_fcts[k], n) v_i = v_fcts[i] - gv_all = ntuple(k -> gv_fcts[k], n) function multiple_outputs_α!(dα, α, p, t) gv_buf = p @inbounds for k = 1:n - gv_buf[k] = gv_t[k](t) + gv_buf[k] = gv[k](t) end vi_t = v_i(t) @inbounds for j = 1:n @@ -157,7 +160,7 @@ function effective_output_mode(v_fcts, gv_fcts, T, i; alg = Tsit5(), kwargs...) α_t = sol_α(t) result = v_i(t) @inbounds for k = 1:n - result += gv_all[k](t)' * α_t[k] + result += gv[k](t)' * α_t[k] end return result end @@ -180,8 +183,8 @@ function effective_output_mode( alg = Tsit5(), kwargs..., ) - v_fcts = [_mode_interp(v_, T) for v_ in v_data] - gv_fcts = [_mode_interp(gv_, T) for gv_ in gv_data] + v_fcts = [_interpolate_mode(v_, T) for v_ in v_data] + gv_fcts = [_interpolate_mode(gv_, T) for gv_ in gv_data] return effective_output_mode(v_fcts, gv_fcts, T, i; alg, kwargs...) end @@ -217,15 +220,14 @@ All kwargs are passed on to the ODE solver. function effective_input_mode(u_fcts, gu_fcts, T, i; alg = Tsit5(), kwargs...) @assert i > 1 n = i - 1 - # Capture as tuples for concrete closure types - gu_t = ntuple(k -> gu_fcts[k], n) + # Capture as a tuple for concrete closure types + gu = ntuple(k -> gu_fcts[k], n) u_i = u_fcts[i] gu_i = gu_fcts[i] - gu_all = ntuple(k -> gu_fcts[k], n) function multiple_inputs_α!(dα, α, p, t) gu_buf = p @inbounds for k = 1:n - gu_buf[k] = gu_t[k](t) + gu_buf[k] = gu[k](t) end ui_t = u_i(t) gui_t = gu_i(t) @@ -246,7 +248,7 @@ function effective_input_mode(u_fcts, gu_fcts, T, i; alg = Tsit5(), kwargs...) α_t = sol_α(t) result = u_i(t) @inbounds for k = 1:n - result -= gu_all[k](t)' * α_t[k] + result -= gu[k](t)' * α_t[k] end return result end @@ -264,8 +266,8 @@ function effective_input_mode( alg = Tsit5(), kwargs..., ) - u_fcts = [_mode_interp(u_, T) for u_ in u_data] - gu_fcts = [_mode_interp(gu_, T) for gu_ in gu_data] + u_fcts = [_interpolate_mode(u_, T) for u_ in u_data] + gu_fcts = [_interpolate_mode(gu_, T) for gu_ in gu_data] return effective_input_mode(u_fcts, gu_fcts, T, i; alg, kwargs...) end @@ -288,30 +290,30 @@ effective_input_mode( # coupling_delay_out / coupling_delay_in # ────────────────────────────────────────────── -function _compute_coupling_delay(num_mode::Vector, u::Vector, v::Vector, T::Vector) - l_T = length(T) - buf = Vector{Float64}(undef, l_T) - map!(abs2, buf, u) - ∫u2 = cumul_integrate(T, buf) - map!(abs2, buf, v) - ∫v2 = cumul_integrate(T, buf) - g = zeros(ComplexF64, l_T) - @inbounds for i = 2:l_T - d = abs(∫v2[i] - ∫u2[i]) - if sqrt(abs(d)) > _tol_div - g[i] = num_mode[i]' / sqrt(d + _ϵ) +function _delay_coupling_from_modes(num_mode::Vector, u::Vector, v::Vector, T::Vector) + nt = length(T) + mode_sq = Vector{Float64}(undef, nt) + map!(abs2, mode_sq, u) + ∫u2 = cumul_integrate(T, mode_sq) + map!(abs2, mode_sq, v) + ∫v2 = cumul_integrate(T, mode_sq) + coupling = zeros(ComplexF64, nt) + @inbounds for i = 2:nt + denom = abs(∫v2[i] - ∫u2[i]) + if sqrt(abs(denom)) > _tol_div + coupling[i] = num_mode[i]' / sqrt(denom + _ϵ) end end - return _as_time_fn(LinearInterpolation(g, T; extrapolation = _extrapolate)) + return PulseCoupling(LinearInterpolation(coupling, T; extrapolation = _extrapolate)) end """ coupling_delay_out(u, v, T) Compute the out-coupling strength for a delay cavity. -Returns a single-argument closure `g(t)`. +Returns a [`PulseCoupling`](@ref) (callable as `g(t)`, broadcastable as `g.(T)`). """ -coupling_delay_out(u::Vector, v::Vector, T::Vector) = _compute_coupling_delay(u, u, v, T) +coupling_delay_out(u::Vector, v::Vector, T::Vector) = _delay_coupling_from_modes(u, u, v, T) coupling_delay_out(u::Function, v::Function, T::Vector) = coupling_delay_out(u.(T), v.(T), T) coupling_delay_out(u::LinearInterpolation, v::LinearInterpolation, T::Vector) = @@ -321,9 +323,9 @@ coupling_delay_out(u::LinearInterpolation, v::LinearInterpolation, T::Vector) = coupling_delay_in(u, v, T) Compute the in-coupling strength for a delay cavity. -Returns a single-argument closure `g(t)`. +Returns a [`PulseCoupling`](@ref) (callable as `g(t)`, broadcastable as `g.(T)`). """ -coupling_delay_in(u::Vector, v::Vector, T::Vector) = _compute_coupling_delay(-v, u, v, T) +coupling_delay_in(u::Vector, v::Vector, T::Vector) = _delay_coupling_from_modes(-v, u, v, T) coupling_delay_in(u::Function, v::Function, T::Vector) = coupling_delay_in(u.(T), v.(T), T) coupling_delay_in(u::LinearInterpolation, v::LinearInterpolation, T::Vector) = coupling_delay_in(u.(T), v.(T), T) diff --git a/test/test_code_quality.jl b/test/test_code_quality.jl index 7e8e89ce..8c8b30a9 100644 --- a/test/test_code_quality.jl +++ b/test/test_code_quality.jl @@ -38,4 +38,5 @@ end all_concrete(typeof(QIO.SLH(1, 0, 0))) # TODO all_concrete(QIO.Gaussian) + all_concrete(QIO.PulseCoupling) end diff --git a/test/test_translate.jl b/test/test_translate.jl index 471499db..5fe1168c 100644 --- a/test/test_translate.jl +++ b/test/test_translate.jl @@ -5,10 +5,6 @@ using QuantumOpticsBase using QuantumOpticsBase: dagger, static_operator, TimeDependentSum using Test -# Since SecondQuantizedAlgebra v0.10, `to_numeric` returns a `TimeDependentSum` for -# time-dependent inputs instead of a plain callable. Calling it at a time `t` sets the -# current time and still returns a `TimeDependentSum`; `mat` materializes it to a dense -# `Operator` for numeric comparison. mat(F, t) = dense(static_operator(F(t))) @testset "translate" begin @@ -240,4 +236,26 @@ mat(F, t) = dense(static_operator(F(t))) @test F isa TimeDependentSum @test sum(abs.((mat(F, 0.0) - dense((im * (1.0 + 2.0im)) * a3_QO)).data)) < 1e-8 end + + @testset "PulseCoupling as time_parameter" begin + # A `PulseCoupling` from a coupling constructor must plug into `to_numeric`'s + # `time_parameter` without the conflicting-arity error a raw interpolation triggers. + hcp = FockSpace(:cp) + acp = Destroy(hcp, :a) + bcp = FockBasis(3) + acp_QO = destroy(bcp) + @variables gu::Complex + Tg = collect(0.0:0.05:1.0) + umode(t) = exp(-(t - 0.5)^2) + gu_t = coupling_input(umode, Tg) + @test gu_t isa PulseCoupling + Fpc = to_numeric( + acp * gu, + bcp; + time_parameter = Dict(gu => gu_t), + operators = Dict([acp, acp'] .=> [acp_QO, dagger(acp_QO)]), + ) + @test Fpc isa TimeDependentSum + @test sum(abs.((mat(Fpc, 0.3) - dense(gu_t(0.3) * acp_QO)).data)) < 1e-10 + end end diff --git a/test/test_utils.jl b/test/test_utils.jl index 1ee3b631..e7637800 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -23,6 +23,14 @@ using Test gu_ana = coupling_input(Gaussian(τ, σ; δ = δ)) gv_ana = coupling_output(Gaussian(τ, σ; δ = δ)) + # Every coupling constructor returns a `PulseCoupling`: callable and broadcastable. + @test gu_num isa PulseCoupling + @test gv_num isa PulseCoupling + @test gu_ana isa PulseCoupling + @test gv_ana isa PulseCoupling + @test gu_num(T[5]) isa ComplexF64 + @test gu_num.(T[1:3]) isa AbstractVector{ComplexF64} + gv_err = maximum(abs.(gv_num.(T[2:end]) .- gv_ana.(T[2:end]))) gu_err = maximum(abs.(gu_num.(T[2:end]) .- gu_ana.(T[2:end]))) From 411ac715391d232e36a9786b12eebeedc39aa7a0 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Thu, 23 Jul 2026 17:35:37 +0200 Subject: [PATCH 5/7] test coverage --- ...26-07-23-qio-v0.10-numeric-layer-design.md | 181 ------------------ test/test_correlations.jl | 15 +- test/test_utils.jl | 12 ++ 3 files changed, 26 insertions(+), 182 deletions(-) delete mode 100644 docs/superpowers/specs/2026-07-23-qio-v0.10-numeric-layer-design.md diff --git a/docs/superpowers/specs/2026-07-23-qio-v0.10-numeric-layer-design.md b/docs/superpowers/specs/2026-07-23-qio-v0.10-numeric-layer-design.md deleted file mode 100644 index 063fa0c2..00000000 --- a/docs/superpowers/specs/2026-07-23-qio-v0.10-numeric-layer-design.md +++ /dev/null @@ -1,181 +0,0 @@ -# QIO v0.10-native numeric layer - -Date: 2026-07-23 -Status: Draft (awaiting review) - -## Context - -QuantumInputOutput.jl (QIO) was written against SecondQuantizedAlgebra (SQA) v0.9, whose -numeric model was: a time-dependent numeric operator **is a callable** `t -> operator`, and a -`time_parameter` value **is any callable** `t -> value`. Two QIO subsystems were shaped by -that model and now read as mismatches against SQA v0.10. - -1. **Couplings return a `LinearInterpolation`.** SQA v0.10 introduced parameter-aware - coefficients (`(p, t) -> value`) and now **introspects the arity** of every - `time_parameter` value (`_tp_reads_p` in `SQA/src/numeric/coeff.jl`): it requires a single, - unambiguous arity (1 for `t`, 2 for `(p, t)`) and rejects anything else. A - `DataInterpolations.LinearInterpolation`'s call operator carries methods of more than one - arity (a 1-arg scalar/vector form plus 2-arg vector forms; note there is no `itp(t, deriv)` - scalar method), so `_tp_reads_p` cannot resolve a single arity and rejects it with - `ArgumentError: ... has methods of conflicting arity [1, 2]; write it as an explicit - t -> value or (p, t) -> value closure`. QIO papered over this with `_as_time_fn` (an - anonymous `t -> itp(t)` wrapper, exactly the closure SQA's error prescribes). The wrapper is - v0.10-correct, but it signals that QIO leaks an interpolation *object* where SQA's currency - is a *time function*. - -2. **SLH composition uses `FunctionWrapper`/callables.** `SLH.jl` carries a callable-composition - layer (`_Callable = Union{Function, FunctionWrapper}`, the `_add`/`_mul`/`_adj` closure - branches, `_detect_operator_type`, `_maybe_wrap_lindblad`, `_maybe_wrap_hamiltonian`, - `_fw_return_type`, `_op_type`, `_to_func`). It exists to compose SLH triples whose numeric - time-dependent entries are closures `t -> operator`, and to pin their return type via - `FunctionWrapper{OpType, Tuple{Float64}}` for type stability. This looked like a v0.9 relic, - but measurement (see "Investigation: is `FunctionWrapper` still needed?" below) shows it - provides concreteness that native `TimeDependentSum` composition does not replace. It is - therefore **kept** in this PR and revisited in the future backend-agnostic PR. - -### How SQA v0.10 consumes a `time_parameter` value (verified by reading the pipeline) - -`api.jl` calls `core.jl:_to_numeric_td`, which calls `coeff.jl:_time_basis`/`_td_coeff`. - -- `_time_basis` collects the value functions into a **`Tuple`** and runs `_tp_reads_p` on each. -- `_td_coeff` compiles the symbolic prefactor with `build_function` and returns the per-term - coefficient closure `t -> pref(map(f -> f(t), valuefuncs)...)`. -- The QuantumOptics backend stores these as `TimeDependentSum(ComplexF64, b, b, coeffs, ops)`. - -SQA's contract is therefore exact: **a `time_parameter` value must be a single-arity, concrete -callable returning a `Number`.** For `map(f -> f(t), valuefuncs)` to stay type-stable, each -value function must be a concrete callable type. - -## Goals - -Make QIO's numeric layer speak SQA v0.10 natively, while staying representation-agnostic so a -later backend-neutral PR needs no rework here. - -## Non-goals (explicitly deferred to a future PR) - -- Moving `QuantumOptics`/`QuantumOpticsBase` to weakdeps behind package extensions. -- QuantumToolbox support, or any second numeric backend. -- Any backend-hook abstraction for the solver/correlation layer. - -This PR stays **QuantumOptics-only**, but every new type and rewrite here is **backend-agnostic -by construction** (it names no backend-specific numeric type), so the future PR inherits it -unchanged. - -## Deliverable 1: `PulseCoupling`, a first-class time-function type - -### The type - -```julia -""" - PulseCoupling(f) - -A single-argument time function `g(t)` used as a `to_numeric` `time_parameter` value. -Wraps any evaluator `f` (a sampled interpolation, an analytic closure, ...) behind exactly -one call method, so it satisfies SQA's single-arity `time_parameter` contract by -construction. Callable as `g(t)` and broadcastable as `g.(T)`. -""" -struct PulseCoupling{F} - f::F -end -(g::PulseCoupling)(t::Real) = g.f(t) -``` - -- **Single-arity by construction**, so it satisfies `_tp_reads_p` with no wrapper. - `_as_time_fn` is deleted. -- **Concrete/static**: `PulseCoupling{LinearInterpolation{...}}` is concrete; it lands in SQA's - `valuefuncs` `Tuple` concretely, keeping `map(f -> f(t), ...)` inferrable, and it returns a - concrete `ComplexF64`. -- **Performant**: one direct call into `f`, no extra allocation, identical hot-path cost to the - raw interpolation minus the arity ambiguity. -- **Backend-agnostic**: a scalar time function that names no QuantumOptics/QuantumToolbox type. -- **Public/introspectable**: documented and exported like `Gaussian`, with `g.f` reachable. - -### Routing - -Every coupling constructor returns a `PulseCoupling`. - -- `coupling_input(u::Vector, T)`, `coupling_output`, `coupling_delay_in`, `coupling_delay_out`: - wrap the computed `LinearInterpolation` in `PulseCoupling` (replacing `_as_time_fn`). -- `coupling_input(::Gaussian)`, `coupling_output(::Gaussian)`: wrap the analytic closure in - `PulseCoupling` too, for a uniform return type. These were already single-arity, so this is - for consistency, not correctness. -- Internal uses that relied on the return being callable or broadcastable (for example - `coupling_output(v_, T).(T)` and the `effective_*_mode` collections in `pulses.jl`) keep - working, because `PulseCoupling` is callable and `g.(T)` broadcasts. - -`_as_time_fn` is removed. The existing `coupling_*(u::LinearInterpolation, T)` input-mode -dispatch methods stay, since they accept an interpolation as the *input mode*, unrelated to the -return type. - -### Interaction picture (light touch) - -`interaction_picture.jl` builds matrix-valued time functions (`M(t)`), and example 08-1 extracts -scalar elements with `t -> M_t(t)[i, j]` before putting them in `time_parameter`. Those -per-element accessors are the same "scalar time function" currency, so they can be wrapped as -`PulseCoupling(t -> M_t(t)[i, j])` at the point they enter `time_parameter`. The matrix-valued -API itself does not change; this is an example-level consistency update. - -## Investigation: is `FunctionWrapper` still needed? (yes) - -The `FunctionWrapper` layer looked like a v0.9 relic to delete in favour of native -`TimeDependentSum` arithmetic. Direct measurement (QuantumOptics backend, `FockBasis`) shows -otherwise: - -| Probe | Result | -|:--|:--| -| Two SQA `TimeDependentSum` on the same basis share one concrete type | yes | -| `SVector(tds, tds)` eltype concrete | yes | -| `SVector(static_op, tds)` eltype | `AbstractOperator` (**not** concrete) | -| A single `tds`'s internal containers (coeffs, `LazySum` ops) | abstract (`Vector{Function}`, `Vector{AbstractOperator}`) even though the outer wrapper type is concrete | -| `tds + tds`, `scalar * tds`, `adjoint(tds)` stay `<: TimeDependentSum` | yes, but each returns a **differently-parameterized** concrete type (`typeof(x) != typeof(tds)`) | -| `tds * tds` preserves the **exact** concrete type | yes | -| `SLH(tds) ▷ SLH(tds)` (cascade) | **errors** (`MethodError: zero(::TimeDependentSum)`) | -| `iszero(tds)` | **errors** (`MethodError`) | - -Conclusions: - -- `FunctionWrapper` provided a **single uniform concrete type** for mixed and composed numeric - SLH entries. Native `TimeDependentSum` composition does not: a mixed static/time-dependent - Lindblad vector degrades to `SVector{N, AbstractOperator}`; the outer `TimeDependentSum` type - is concrete but wraps abstract containers; and `+`/`scalar *`/`adjoint` each return a - differently-parameterized `TimeDependentSum` (only `*` preserves the exact concrete type), so a - heterogeneous collection of composed entries still widens to an abstract eltype. -- Restoring that concreteness the v0.10 way would mean promoting every entry to a - `TimeDependentSum`. But **QuantumToolbox has no `TimeDependentSum`** (it uses `QobjEvo`), so a - promotion step is backend-specific and must live in the future backend-agnostic PR, alongside - the per-backend abstraction. -- Separately, numeric SLH composition on `TimeDependentSum` entries currently **errors** (the - `iszero`-on-operator bug), and no example exercises time-dependent numeric SLH composition, so - nothing in the current workflows is blocked. - -**Decision:** keep `FunctionWrapper` unchanged in this PR. Its removal, the concreteness-restoring -promotion, and the `iszero`-on-operator fix are all entangled with the per-backend abstraction and -are deferred to the future backend-agnostic PR. - -## Testing - -- **`PulseCoupling`**: unit tests that a coupling `isa PulseCoupling`, is callable, broadcasts, - returns `ComplexF64`, and is accepted by `to_numeric`'s `time_parameter` without error (the - case that previously required `_as_time_fn`). Keep the existing pulse-coupling numeric - comparisons. -- Full `Pkg.test()` green, and all examples still run headless. The `SLH.jl` and `test_SLH.jl` - `FunctionWrapper` tests are untouched. - -## Migration and back-compat - -- Examples using `coupling_*` are unaffected (still callable and broadcastable). -- Example 08-1 interaction-picture element accessors optionally wrapped in `PulseCoupling`. -- No breaking changes to SLH construction or composition. - -## Risks - -- **`PulseCoupling` broadcast/dispatch**: confirm `g.(T)` still broadcasts and that the internal - `effective_*_mode` collections (which call and broadcast couplings) are unaffected. Low risk, - covered by the existing pulse tests. - -## Scope note for review - -This PR is scoped to Deliverable 1 (`PulseCoupling`) only. The `FunctionWrapper` question is -resolved above: it stays, and the SLH-composition modernization is deferred to the future -backend-agnostic PR where it can be done correctly (per-backend promotion, not -`TimeDependentSum`-specific code). diff --git a/test/test_correlations.jl b/test/test_correlations.jl index 718aa86c..eb4ca207 100644 --- a/test/test_correlations.jl +++ b/test/test_correlations.jl @@ -1,6 +1,6 @@ using QuantumInputOutput using QuantumOptics -using QuantumOpticsBase: dagger +using QuantumOpticsBase: dagger, TimeDependentSum using LinearAlgebra using Test @@ -45,4 +45,17 @@ using Test mode_num_aligned = mode_num * exp(-1im * angle(overlap)) overlap_aligned = abs(sum(conj.(mode_exp) .* mode_num_aligned) * ΔT) @test overlap_aligned > 0.99 + + @testset "time-dependent Ls rejected" begin + td_Ls = TimeDependentSum([t -> 1.0 + 0im], J) + @test_throws ArgumentError correlation_matrix(T, ρt_static, H, J, td_Ls) + end + + @testset "time-dependent operator path matches static" begin + H_td = TimeDependentSum([t -> 0.0 + 0im], [one(ba)]) + J_td = [TimeDependentSum([t -> sqrt(γ) + 0im], [σm])] + _, ρt_td = timeevolution.master_dynamic(T, ψ0, H_td, J_td) + g1_td = correlation_matrix(T, ρt_td, H_td, J_td, Ls) + @test maximum(abs.(g1_td .- g1_static)) < 1e-8 + end end diff --git a/test/test_utils.jl b/test/test_utils.jl index e7637800..7d2f8c70 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -61,4 +61,16 @@ using Test u_eff_data2 = effective_input_mode([u1.(T), u2.(T)], [gu1.(T), gu2.(T)], T, 2) @test maximum(abs.(u_eff_f.(T) .- u_eff_data2.(T))) / maximum(abs.(u_eff_f.(T) .+ u_eff_data2.(T))) < 1e-2 + + v3(t) = 1 / (sqrt(σ) * π^(1/4)) * exp(-0.5 * ((t - τ) / σ)^2) + v_eff3_f = effective_output_mode([v1, v2, v3], T, 3) + v_eff3_data = effective_output_mode([v1.(T), v2.(T), v3.(T)], T, 3) + @test maximum(abs.(v_eff3_f.(T) .- v_eff3_data.(T))) / + maximum(abs.(v_eff3_f.(T) .+ v_eff3_data.(T))) < 1e-2 + + u3(t) = 1 / (sqrt(σ) * π^(1/4)) * exp(-0.5 * ((t - τ) / σ)^2) + u_eff3_f = effective_input_mode([u1, u2, u3], T, 3) + u_eff3_data = effective_input_mode([u1.(T), u2.(T), u3.(T)], T, 3) + @test maximum(abs.(u_eff3_f.(T) .- u_eff3_data.(T))) / + maximum(abs.(u_eff3_f.(T) .+ u_eff3_data.(T))) < 1e-2 end From 98d9f1c34b424124627510c597912686917982b8 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Thu, 23 Jul 2026 18:33:00 +0200 Subject: [PATCH 6/7] test: strengthen coupling/effective-mode coverage into correctness checks --- test/test_utils.jl | 43 +++++++++++++++++++++++++++++++++---------- 1 file changed, 33 insertions(+), 10 deletions(-) diff --git a/test/test_utils.jl b/test/test_utils.jl index 7d2f8c70..15f2448b 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -1,4 +1,5 @@ using QuantumInputOutput +using DataInterpolations: LinearInterpolation using Test @testset "utils" begin @@ -62,15 +63,37 @@ using Test @test maximum(abs.(u_eff_f.(T) .- u_eff_data2.(T))) / maximum(abs.(u_eff_f.(T) .+ u_eff_data2.(T))) < 1e-2 + # zero coupling ⇒ no distortion ⇒ effective mode equals the bare mode (exercises i = 3) + zero_c(t) = 0.0 + 0.0im v3(t) = 1 / (sqrt(σ) * π^(1/4)) * exp(-0.5 * ((t - τ) / σ)^2) - v_eff3_f = effective_output_mode([v1, v2, v3], T, 3) - v_eff3_data = effective_output_mode([v1.(T), v2.(T), v3.(T)], T, 3) - @test maximum(abs.(v_eff3_f.(T) .- v_eff3_data.(T))) / - maximum(abs.(v_eff3_f.(T) .+ v_eff3_data.(T))) < 1e-2 - - u3(t) = 1 / (sqrt(σ) * π^(1/4)) * exp(-0.5 * ((t - τ) / σ)^2) - u_eff3_f = effective_input_mode([u1, u2, u3], T, 3) - u_eff3_data = effective_input_mode([u1.(T), u2.(T), u3.(T)], T, 3) - @test maximum(abs.(u_eff3_f.(T) .- u_eff3_data.(T))) / - maximum(abs.(u_eff3_f.(T) .+ u_eff3_data.(T))) < 1e-2 + v_eff3 = effective_output_mode([v1, v2, v3], [zero_c, zero_c, zero_c], T, 3) + u_eff3 = effective_input_mode([u1, u2, v3], [zero_c, zero_c, zero_c], T, 3) + @test maximum(abs.(v_eff3.(T) .- v3.(T))) < 1e-8 + @test maximum(abs.(u_eff3.(T) .- v3.(T))) < 1e-8 + + # constant modes give the closed form g_out(tᵢ) = u / sqrt((|v|² - |u|²) i Δt) + a_c = 0.5 + b_c = 0.8 + Td = collect(0.0:0.01:2.0) + dt_c = Td[2] - Td[1] + uc(t) = a_c + vc(t) = b_c + g_out = coupling_delay_out(uc, vc, Td) + g_in = coupling_delay_in(uc, vc, Td) + exp_out = [a_c / sqrt((b_c^2 - a_c^2) * (i * dt_c)) for i = 2:length(Td)] + exp_in = [-b_c / sqrt((b_c^2 - a_c^2) * (i * dt_c)) for i = 2:length(Td)] + @test g_out isa PulseCoupling + @test maximum(abs.(g_out.(Td[2:end]) .- exp_out)) < 1e-6 + @test maximum(abs.(g_in.(Td[2:end]) .- exp_in)) < 1e-6 + # the shared denominator cancels, so in/out = -v/u exactly + @test maximum(abs.(g_in.(Td[2:end]) ./ g_out.(Td[2:end]) .- (-b_c / a_c))) < 1e-12 + + # LinearInterpolation dispatch must reproduce the analytic references + u_itp = LinearInterpolation(u.(T), T) + @test maximum(abs.(coupling_input(u_itp, T).(T[2:end]) .- gu_ana.(T[2:end]))) < 5e-4 + @test maximum(abs.(coupling_output(u_itp, T).(T[2:end]) .- gv_ana.(T[2:end]))) < 5e-4 + uc_itp = LinearInterpolation(uc.(Td), Td) + vc_itp = LinearInterpolation(vc.(Td), Td) + @test maximum(abs.(coupling_delay_out(uc_itp, vc_itp, Td).(Td[2:end]) .- exp_out)) < 1e-6 + @test maximum(abs.(coupling_delay_in(uc_itp, vc_itp, Td).(Td[2:end]) .- exp_in)) < 1e-6 end From ae175688a6df3f226537f2c5f8817fc8bfe18f35 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Thu, 23 Jul 2026 18:33:13 +0200 Subject: [PATCH 7/7] format --- test/test_utils.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_utils.jl b/test/test_utils.jl index 15f2448b..e92fadf0 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -94,6 +94,7 @@ using Test @test maximum(abs.(coupling_output(u_itp, T).(T[2:end]) .- gv_ana.(T[2:end]))) < 5e-4 uc_itp = LinearInterpolation(uc.(Td), Td) vc_itp = LinearInterpolation(vc.(Td), Td) - @test maximum(abs.(coupling_delay_out(uc_itp, vc_itp, Td).(Td[2:end]) .- exp_out)) < 1e-6 + @test maximum(abs.(coupling_delay_out(uc_itp, vc_itp, Td).(Td[2:end]) .- exp_out)) < + 1e-6 @test maximum(abs.(coupling_delay_in(uc_itp, vc_itp, Td).(Td[2:end]) .- exp_in)) < 1e-6 end