From 7e40fb44ae98c226eefd824d3212d5f74ce6328f Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 11 Jul 2026 20:43:27 +0200 Subject: [PATCH 1/4] feat: update to SQA v0.9 --- Project.toml | 2 +- TODO.md | 31 +- benchmarks/correlations.jl | 8 +- benchmarks/interaction_picture.jl | 4 +- benchmarks/translation.jl | 12 +- docs/src/api.md | 4 +- docs/src/implementation.md | 8 +- docs/src/tutorial.md | 14 +- ...cattering__PRL2019_123-123604_fig2-fig3.jl | 20 +- ...lated-emission__PRL2019_123-123604_fig4.jl | 12 +- ...y-phase-noise__PRA2020_102- 023717_fig2.jl | 4 +- ...e-entanglement__PRA2020_102-023717_fig4.jl | 8 +- ...-combiner__PRA2023_107-023715_fig2-fig3.jl | 14 +- ...-1_two-sided-cavity_with-atom_coh-drive.jl | 16 +- ..._bidirectional-waveguide_coherent-pulse.jl | 8 +- ...irectional-waveguide_feedback-reduction.jl | 8 +- ...action-picture__PRA2023_107-013706_fig2.jl | 8 +- .../07-1_beamsplitter_loss__quantum-pulse.jl | 4 +- .../07-2_hong-ou-mandel__quantum-pulse.jl | 4 +- examples/08-1_pulse-delay__simple.jl | 22 +- examples/10-1_SUPER_excitation.jl | 10 +- ...ing-cat-state__PRA2020_102-023717_fig3b.jl | 8 +- ...pagation-delay__PRA2026_113-013730_fig4.jl | 20 +- examples/drafts/10-1_SUPER_excitation.jl | 6 +- src/QuantumInputOutput.jl | 15 +- src/SLH.jl | 20 + src/translate.jl | 358 ------------------ src/utils.jl | 39 -- test/test_compare_example_05_1_05_2.jl | 8 +- test/test_example_cavity_scattering.jl | 14 +- test/test_interaction_picture.jl | 16 +- test/test_translate.jl | 74 ++-- 32 files changed, 197 insertions(+), 602 deletions(-) delete mode 100644 src/translate.jl delete mode 100644 src/utils.jl diff --git a/Project.toml b/Project.toml index 191405a..60bab9f 100644 --- a/Project.toml +++ b/Project.toml @@ -32,7 +32,7 @@ OrdinaryDiffEq = "6" QuantumCumulants = "0.5" QuantumOptics = "1" QuantumOpticsBase = "0.5" -SecondQuantizedAlgebra = "0.8.2" +SecondQuantizedAlgebra = "0.9" SpecialFunctions = "2" StaticArrays = "1" SymbolicUtils = "4" diff --git a/TODO.md b/TODO.md index f8c64b1..6ba13c9 100644 --- a/TODO.md +++ b/TODO.md @@ -1,7 +1,7 @@ # TODOs - [ ] examples - - [ ] SUPER example + - [x] SUPER example - [ ] pulse delay advanced example - [x] feedback reduction (see SLH paper) - [ ] cumulants correlation matrix modes @@ -24,36 +24,9 @@ ## Additional TODOs - [x] u_eff and v_eff equation (PRA 2020, appendix) theory, API, implementation -- [ ] rename translate_qo (to to_numerics ?) +- [x] rename translate_qo (now `to_numeric`, from SQA v0.9; `substitute_operators` now `substitute`) - -## Architecture & maintainability (review after SQA 0.6 / QC 0.5 migration) - -### Layering: SQA vs QC vs QIO -Clean intended split: SQA owns symbolic operator algebra plus the numeric backend; QC owns dynamics derivation (meanfield/cumulants) plus the ModelingToolkit ODE bridge; QIO owns SLH network composition plus the pulse/input-output physics. The split is mostly real. Every leak is in one spot: the "symbolic to numeric coefficient" bridge. - -Every hard bug in the migration was the same bug rediscovered at a different layer ("reduce or compile a `Complex{Num}` coefficient under Symbolics 7"): -- `to_numeric` not constant-folding `sqrt(1.0)` (worked around in QIO `_translate_numeric`). -- `build_function` emitting `complex(re, im)` which throws on complex inputs (QC already solved this with `_im_form` / `_qc_maketerm`; QIO re-solved it with `_compile_coeff`). -- Substituting a `Complex{Num}` parameter key not matching (QIO `_expand_parameter`). -- `substitute_operators`: SQA `substitute` stack-overflows on self-referential operator replacement (`a => 2a`); QIO hand-rolled dict-lookup; QC already has robust `rewrite`. - -Conclusion: QC and QIO each carry a private copy of the coefficient-codegen machinery and each independently ate the Symbolics 7 upgrade pain. That duplication is the core smell. - -### translate_qo really belongs in SQA (completing the QO bridge) -SQA already translates to QuantumOptics (`to_numeric` in `SQA/src/numeric.jl`, already depends on QuantumOpticsBase). `translate_qo` is not a new capability; it is `to_numeric` completed with the three things `to_numeric` punts on: parameter substitution, constant-folding, and time-dependence. Symmetry: SQA to QuantumOptics is one numeric backend (currently shipped half in SQA, half stranded in QIO); QC to ModelingToolkit is the other (fully owned by QC). The complete QO backend is SQA's to own. - -- [ ] Move the generic bridge into SQA (ideally a QuantumOptics package extension so SQA core stays algebra-only): `to_numeric(op, basis; parameters, time_parameters)` returning an operator or a `t -> op(t)`; parameter substitution incl. `Complex{Num}` keys; constant-folding; the `im`-recombination coefficient compiler (one robust copy). -- [ ] Keep in QIO only the domain-specific bits: the `SLH` overload, the `operators` kwarg (substitute symbolic ops with user QO ops on another basis), the `adjoint_ops` convenience, and pulse/interaction-picture/correlation physics. QIO `translate_qo` then becomes a thin SLH wrapper over the SQA primitive. -- [ ] This collapses QC `_im_form` and QIO `_compile_coeff` into one shared path; fix the Symbolics-version pain once, not N times. -- Caveats: `t -> op(t)` is solver-flavored (be deliberate that the abstraction is "time-dependent operator", not "thing for master_dynamic"); time path pulls `FunctionWrappers` + `build_function` (small, QOBase already a dep); this is a qojulia-wide scope decision for the SQA maintainers. -- Related: the existing "rename translate_qo (to to_numerics ?)" item above. - -### Things to fix in SQA (or lift from QC) -- [ ] `to_numeric` should accept a parameter dict and constant-fold residual symbolic constants (`sqrt`, `exp`), instead of erroring in `_to_complex`. -- [ ] Robust operator substitution: fix SQA `substitute` for self-referential replacement, or lift QC `rewrite`/`_qc_maketerm` into SQA; then QIO `substitute_operators` is a thin wrapper. - ### QIO maintainability wins - [ ] Examples in CI as smoke tests (highest leverage). Examples were untested and rotted: stale API, several did not run. Real bugs (04-2 MTK flow, `get_solution` new signature, complex-coupling crash, `QAdd + BasicSymbolic` scalar drive) were only caught by running them by hand. A CI job running each example's symbolic setup plus a tiny ODE (skip plots/long solves) would catch API drift. - [ ] Document the `::Real` vs `::Complex` contract for couplings. The `conj` in `g'` is load-bearing: a coupling that is physically complex must be `::Complex` or the physics is silently wrong. Old `rnumbers` made couplings real; that was a latent modeling bug. Consider a runtime check ("real parameter received a complex value"). diff --git a/benchmarks/correlations.jl b/benchmarks/correlations.jl index 99b7524..3e31a3b 100644 --- a/benchmarks/correlations.jl +++ b/benchmarks/correlations.jl @@ -39,8 +39,8 @@ function benchmark_correlations!(SUITE) dict_p = Dict([γ, Δ, g_v] .=> [γ_, 0.0, 0]) dict_p_t = Dict(g_u => gu_t) - H_QO = translate_qo(H_sym, b; parameter = dict_p, time_parameter = dict_p_t) - L_QO = translate_qo(L_sym, b; parameter = dict_p, time_parameter = dict_p_t) + 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) @@ -51,8 +51,8 @@ function benchmark_correlations!(SUITE) ψ0 = fockstate(bu1, 1) ⊗ fockstate(bc1, 0) ⊗ fockstate(bv1, 0) _, ρt = timeevolution.master_dynamic(T, ψ0, input_output) - au_qo = translate_qo(au, b) - c_qo = translate_qo(c, b) + au_qo = to_numeric(au, b) + c_qo = to_numeric(c, b) Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo ## --- Two-time correlation --- diff --git a/benchmarks/interaction_picture.jl b/benchmarks/interaction_picture.jl index 755b924..cb743be 100644 --- a/benchmarks/interaction_picture.jl +++ b/benchmarks/interaction_picture.jl @@ -81,8 +81,8 @@ function benchmark_interaction_picture!(SUITE) SUITE["Interaction Picture"]["operator substitution"]["TLS cascade"] = @benchmarkable begin - simplify(substitute_operators($H_int_, $int_dict)) - simplify(substitute_operators($L, $int_dict)) + simplify(substitute($H_int_, $int_dict)) + simplify(substitute($L, $int_dict)) end return nothing diff --git a/benchmarks/translation.jl b/benchmarks/translation.jl index 8b92946..7195870 100644 --- a/benchmarks/translation.jl +++ b/benchmarks/translation.jl @@ -45,7 +45,7 @@ function benchmark_translation!(SUITE) expr_composite = a * 3 + Δ_sym * σ(2, 2) SUITE["Translation"]["static"]["atom-cavity"] = - @benchmarkable translate_qo($expr_composite, $b_multi; parameter = $dict_p_static) + @benchmarkable to_numeric($expr_composite, $b_multi; parameter = $dict_p_static) ## --- Time-dependent translation --- @@ -56,7 +56,7 @@ function benchmark_translation!(SUITE) expr_td = a * 3 * conj(E) + Δ_sym * σ(2, 2) - SUITE["Translation"]["time-dependent"]["atom-cavity"] = @benchmarkable translate_qo( + SUITE["Translation"]["time-dependent"]["atom-cavity"] = @benchmarkable to_numeric( $expr_td, $b_multi; parameter = $dict_p_static, @@ -80,8 +80,8 @@ function benchmark_translation!(SUITE) dict_p_t_cav = Dict([gu_sym, gv_sym] .=> [gu_t, gv_t]) SUITE["Translation"]["time-dependent"]["3-cavity H+L"] = @benchmarkable begin - translate_qo($H_sym, $b_cav; parameter = $dict_p_cav, time_parameter = $dict_p_t_cav) - translate_qo( + to_numeric($H_sym, $b_cav; parameter = $dict_p_cav, time_parameter = $dict_p_t_cav) + to_numeric( $L_sym, $b_cav; parameter = $dict_p_cav, @@ -93,8 +93,8 @@ function benchmark_translation!(SUITE) SUITE["Translation"]["closure evaluation"] = BenchmarkGroup() - H_QO = translate_qo(H_sym, b_cav; parameter = dict_p_cav, time_parameter = dict_p_t_cav) - L_QO = translate_qo(L_sym, b_cav; parameter = dict_p_cav, time_parameter = dict_p_t_cav) + H_QO = to_numeric(H_sym, b_cav; parameter = dict_p_cav, time_parameter = dict_p_t_cav) + L_QO = to_numeric(L_sym, b_cav; parameter = dict_p_cav, time_parameter = dict_p_t_cav) t_mid = T[length(T)÷2] diff --git a/docs/src/api.md b/docs/src/api.md index 1b8eca1..c8e295c 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -41,7 +41,7 @@ feedback ## [Translation](@id API: Translation) ```@docs -translate_qo +to_numeric ``` ## [Pulses](@id API: Pulses) @@ -97,5 +97,5 @@ correlation_matrix ## [Utilities](@id API: Utilities) ```@docs -substitute_operators +substitute ``` diff --git a/docs/src/implementation.md b/docs/src/implementation.md index 2516c39..0f2f89d 100644 --- a/docs/src/implementation.md +++ b/docs/src/implementation.md @@ -37,12 +37,12 @@ If you directly want to use [QuantumOptics.jl](https://github.com/qojulia/Quantu ## Translate to numerics To solve the dynamics of the master equation we first need to create the corresponding numeric operators for the Hamilton and the Lindblad terms. -This can be done with [`translate_qo`](@ref), which converts symbolic operators into [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators on a chosen basis. It accepts two parameter substitution dictionaries: +This can be done with [`to_numeric`](@ref), which converts symbolic operators into [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators on a chosen basis. It accepts two parameter substitution dictionaries: - `parameter`: numeric parameters used in algebraic substitution - `time_parameter`: time-dependent parameters, given as functions of `t` -If `time_parameter` is non-empty, [`translate_qo`](@ref) returns a callable `t -> op(t)` so that the Hamiltonian and jump operators can be supplied to `timeevolution.master_dynamic`. +If `time_parameter` is non-empty, [`to_numeric`](@ref) returns a callable `t -> op(t)` so that the Hamiltonian and jump operators can be supplied to `timeevolution.master_dynamic`. ```julia bu = FockBasis(2) @@ -54,8 +54,8 @@ dict_p = Dict(γ => 1.0, g_v => 0.0) gu_t = coupling_input(t -> exp(-t^2), 0:0.01:5) dict_p_t = Dict(g_u => gu_t) -H_QO = translate_qo(H, b; parameter=dict_p, time_parameter=dict_p_t) -L_QO = translate_qo(L, b; parameter=dict_p, time_parameter=dict_p_t) +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) ``` In some cases it can be useful to define your own set of numeric operators which should replace the symbolic expressions, e.g. to reduce the Hilbert space if the output cavities are not analyzed but they are already included in the symbolic derivation. Such a list of operators can be provide with the dictionary `operators`. diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index b2c22c6..94f64f0 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -87,8 +87,8 @@ bc1 = FockBasis(2) bv1 = FockBasis(2) b = bu1 ⊗ bc1 ⊗ bv1 -H_QO = translate_qo(H, b; parameter=dict_p, time_parameter=dict_p_t) -L_QO = translate_qo(L, b; parameter=dict_p, time_parameter=dict_p_t) +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 ``` @@ -114,9 +114,9 @@ To extract the dominant output mode, we compute the two-time correlation matrix ``g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle`` and diagonalize it. To this end, we first define the desired numerical operators. ```@example tutorial -au_qo = translate_qo(au, b) -c_qo = translate_qo(c, b) -av_qo = translate_qo(av, b) +au_qo = to_numeric(au, b) +c_qo = to_numeric(c, b) +av_qo = to_numeric(av, b) Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo g1_m = correlation_matrix(T, ρt, input_output_1, Ls) @@ -160,8 +160,8 @@ gv_t = coupling_output(v_mode, T) dict_p_t_2 = Dict([g_u, g_v] .=> [gu_t, gv_t]) dict_p_2 = Dict([γ, Δ] .=> [γ_, Δ_]) -H_QO_2 = translate_qo(H, b; parameter=dict_p_2, time_parameter=dict_p_t_2) -L_QO_2 = translate_qo(L, b; parameter=dict_p_2, time_parameter=dict_p_t_2) +H_QO_2 = to_numeric(H, b; parameter=dict_p_2, time_parameter=dict_p_t_2) +L_QO_2 = to_numeric(L, b; parameter=dict_p_2, time_parameter=dict_p_t_2) function input_output_2(t, ρ) Ht = H_QO_2(t) 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 88b1dda..5edf9b4 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 @@ -74,10 +74,10 @@ bv1 = FockBasis(1) b = bu1 ⊗ bc1 ⊗ bv1 nothing # hide -# We now use the function [`translate_qo`](@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 = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = translate_qo(L, b; parameter = dict_p, time_parameter = dict_p_t) +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 @@ -99,9 +99,9 @@ nothing # hide # We create the desired numerical operators to calculate expectation values. -au_qo = translate_qo(au, b) -c_qo = translate_qo(c, b) -av_qo = translate_qo(av, b) +au_qo = to_numeric(au, b) +c_qo = to_numeric(c, b) +av_qo = to_numeric(av, b) n_c_t = real.(expect(c_qo'*c_qo, ρt)) n_u1_t = real.(expect(au_qo'*au_qo, ρt)) @@ -147,8 +147,8 @@ gv_t = coupling_output(v_mode, T) dict_p_t_2 = Dict([g_u, g_v] .=> [gu_t, gv_t]) -H_QO_2 = translate_qo(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) -L_QO_2 = translate_qo(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) +H_QO_2 = to_numeric(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) +L_QO_2 = to_numeric(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) function input_output_2(t, ρ) H = H_QO_2(t) J = [L_QO_2(t)] @@ -213,8 +213,8 @@ 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 -H_QO_3 = translate_qo(H, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2) -L_QO_3 = translate_qo(L, b_3; parameter = dict_p_2, time_parameter = dict_p_t_2) +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, ρ) H = H_QO_3(t) J = [L_QO_3(t), √(γ_)*cdc_3] diff --git a/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl b/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl index e786f37..ca752a4 100644 --- a/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl +++ b/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl @@ -62,7 +62,7 @@ T = [0.001:0.001:1;]*4.0 ΔT = T[2] - T[1] nothing # hide -# We use the function [`translate_qo`](@ref) to create the numeric operators and solve the dynamics with the function `timeevolution.master_dynamic` in QuantumOptics.jl. +# We use the function [`to_numeric`](@ref) to create the numeric operators and solve the dynamics with the function `timeevolution.master_dynamic` in QuantumOptics.jl. ## numeric bases bu1 = FockBasis(1) @@ -70,8 +70,8 @@ ba1 = NLevelBasis(2) bv1 = FockBasis(2) b = bu1 ⊗ ba1 ⊗ bv1 -H_QO = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = translate_qo(L, b; parameter = dict_p, time_parameter = dict_p_t) +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) function input_output(t, ρ) H = H_QO(t) @@ -88,9 +88,9 @@ nothing # hide # To calculate expectation values we create the desired numerical operators. -au_qo = translate_qo(au, b) -σ_qo(i, j) = translate_qo(σ(i, j), b) -av_qo = translate_qo(av, b) +au_qo = to_numeric(au, b) +σ_qo(i, j) = to_numeric(σ(i, j), b) +av_qo = to_numeric(av, b) proj_u(n) = embed(b, 1, projector(fockstate(bu1, n))) proj_v(n) = embed(b, 3, projector(fockstate(bv1, n))) diff --git a/examples/02-1_cavity-phase-noise__PRA2020_102- 023717_fig2.jl b/examples/02-1_cavity-phase-noise__PRA2020_102- 023717_fig2.jl index 2b7a14e..686f0d0 100644 --- a/examples/02-1_cavity-phase-noise__PRA2020_102- 023717_fig2.jl +++ b/examples/02-1_cavity-phase-noise__PRA2020_102- 023717_fig2.jl @@ -74,8 +74,8 @@ c_qo = to_numeric(c, b) cdc_qo = c_qo'c_qo ## translate to numeric Hamiltonian and Lindblad -H_QO = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = translate_qo(L, b; parameter = dict_p, time_parameter = dict_p_t) +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 # We additionally include a cavity dephasing term and solve the dynamics. diff --git a/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl b/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl index 024d481..08eec90 100644 --- a/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl +++ b/examples/02-3_mode-entanglement__PRA2020_102-023717_fig4.jl @@ -80,8 +80,8 @@ nothing # hide # determine the two dominant temporal output modes from the first-order # correlation matrix g⁽¹⁾(t₁, t₂). -H_QO_1 = translate_qo(H, b; parameter = dict_p_1) -L_QO_1 = translate_qo(L, b; parameter = dict_p_1) +H_QO_1 = to_numeric(H, b; parameter = dict_p_1) +L_QO_1 = to_numeric(L, b; parameter = dict_p_1) function input_output_1(t, ρ) J = [L_QO_1] @@ -115,8 +115,8 @@ gv2_t = coupling_output(v2_eff, T) dict_p_2 = Dict([γ, g, ω12] .=> [γ_, g_, ω12_]) dict_p_t_2 = Dict(g_v1 => gv1_t, g_v2 => gv2_t) -H_QO_2 = translate_qo(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) -L_QO_2 = translate_qo(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) +H_QO_2 = to_numeric(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) +L_QO_2 = to_numeric(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) function input_output_2(t, ρ) J = [L_QO_2(t)] diff --git a/examples/03-1_beam-combiner__PRA2023_107-023715_fig2-fig3.jl b/examples/03-1_beam-combiner__PRA2023_107-023715_fig2-fig3.jl index 5879483..2f5c476 100644 --- a/examples/03-1_beam-combiner__PRA2023_107-023715_fig2-fig3.jl +++ b/examples/03-1_beam-combiner__PRA2023_107-023715_fig2-fig3.jl @@ -76,7 +76,7 @@ nothing # hide # To obtain the output modes we do not use the second input mode and the output mode cavity. # However, to keep the example short we include them already from the beginning since they are needed later. # To perform time consuming parameter scans one should merely use the necessary Hilbert spaces. In this case, this would correspond to one input cavity and the two-level system. -# The kwarg `operators` of the function [translate_qo](@ref) provides a convenient way to use predefined numerical operators, see the example `Two-sided Cavity with Atom`. +# The kwarg `operators` of the function [to_numeric](@ref) provides a convenient way to use predefined numerical operators, see the example `Two-sided Cavity with Atom`. ## numeric bases bu2 = FockBasis(2) @@ -85,8 +85,8 @@ bs1 = NLevelBasis(2) bv1 = FockBasis(2) b = bu2 ⊗ bu1 ⊗ bs1 ⊗ bv1; -H_QO = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = translate_qo(L, b; parameter = dict_p, time_parameter = dict_p_t) +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) function input_output(t, ρ) H = H_QO(t) J = [L_QO(t)] @@ -103,8 +103,8 @@ nothing # hide # Now we analyze the output modes with the two-time autocorrelation function $g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle$. -au1_qo = translate_qo(au1, b) -σ_qo(i, j) = translate_qo(σ(i, j), b) +au1_qo = to_numeric(au1, b) +σ_qo(i, j) = to_numeric(σ(i, j), b) Ls(t) = (gu_(t))'*au1_qo + √(γ_)*σ_qo(1, 2) g1_m = correlation_matrix(T, ρt, input_output, Ls); @@ -182,8 +182,8 @@ nothing # hide # The time-dependent couplings are used to define the numeric Hamiltonian and Lindblad term, and then solve the dynamics of the system. -H_QO_2 = translate_qo(H, b; parameter = dict_p_out, time_parameter = dict_p_t_out) -L_QO_2 = translate_qo(L, b; parameter = dict_p_out, time_parameter = dict_p_t_out) +H_QO_2 = to_numeric(H, b; parameter = dict_p_out, time_parameter = dict_p_t_out) +L_QO_2 = to_numeric(L, b; parameter = dict_p_out, time_parameter = dict_p_t_out) function input_output_2(t, ρ) H = H_QO_2(t) J = [L_QO_2(t)] diff --git a/examples/04-1_two-sided-cavity_with-atom_coh-drive.jl b/examples/04-1_two-sided-cavity_with-atom_coh-drive.jl index b835e08..9d34434 100644 --- a/examples/04-1_two-sided-cavity_with-atom_coh-drive.jl +++ b/examples/04-1_two-sided-cavity_with-atom_coh-drive.jl @@ -53,7 +53,7 @@ L1_L = lindblad(G_cav_L_R_drive)[1] L1_R = lindblad(G_cav_L_R_drive)[2] # Here, the usual classical cavity drive-term $\sqrt{\kappa_L} E (a^\dagger + a)$ appears as a combination of Hamiltonian and Lindblad term. -# To solve the dynamics of the system we translate the symbolic expressions into numeric operators (matrices) of [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl). Since we do not want to include the basis of the atoms, we provide a dictionary of operators with the kwarg `operators` in the function [`translate_qo`](@ref). +# To solve the dynamics of the system we translate the symbolic expressions into numeric operators (matrices) of [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl). Since we do not want to include the basis of the atoms, we provide a dictionary of operators with the kwarg `operators` in the function [`to_numeric`](@ref). ## numerical parameters En = 0.5 @@ -72,9 +72,9 @@ bc1 = FockBasis(4) a_QO = destroy(bc1) ops_dict = Dict([a, a'] .=> [a_QO, dagger(a_QO)]) -H1_QO = translate_qo(H1, bc1; parameter = dict_p1, operators = ops_dict) -L1_L_QO = translate_qo(L1_L, bc1; parameter = dict_p1, operators = ops_dict) -L1_R_QO = translate_qo(L1_R, bc1; parameter = dict_p1, operators = ops_dict) +H1_QO = to_numeric(H1, bc1; parameter = dict_p1, operators = ops_dict) +L1_L_QO = to_numeric(L1_L, bc1; parameter = dict_p1, operators = ops_dict) +L1_R_QO = to_numeric(L1_R, bc1; parameter = dict_p1, operators = ops_dict) J1_QO = [L1_L_QO, L1_R_QO] nothing # hide @@ -104,7 +104,7 @@ plot(p1, p2; layout = (2, 1), size = (600, 500)) # Now we scan the laser-cavity detuning $\Delta$ to plot the transmission and reflection spectrum. dict_p_Δ(Δn) = Dict(p_sym .=> [En, κ_Rn, κ_Ln, Δn]) -H1_QO_Δ_(Δn) = translate_qo(H1, bc1; parameter = dict_p_Δ(Δn), operators = ops_dict) +H1_QO_Δ_(Δn) = to_numeric(H1, bc1; parameter = dict_p_Δ(Δn), operators = ops_dict) n_ref_Δ = zeros(lΔ) n_trans_Δ = zeros(lΔ) @@ -185,9 +185,9 @@ a_QO2 = to_numeric(a, b) σ_QO(α, i, j) = to_numeric(σ(α, i, j), b) ## translate to numeric Hamiltonian and Lindblad -H_QO = translate_qo(H2, b; parameter = dict_p2, time_parameter = dict_p_t2) -L2_L_QO = translate_qo(L2_L, b; parameter = dict_p2, time_parameter = dict_p_t2) -L2_R_QO = translate_qo(L2_R, b; parameter = dict_p2) +H_QO = to_numeric(H2, b; parameter = dict_p2, time_parameter = dict_p_t2) +L2_L_QO = to_numeric(L2_L, b; parameter = dict_p2, time_parameter = dict_p_t2) +L2_R_QO = to_numeric(L2_R, b; parameter = dict_p2) ## additional atomic decay into free space J_add = [√(γn)*σ_QO(α, 1, 2) for α = 1:Natoms] 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 d54be23..ab91ea5 100644 --- a/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl +++ b/examples/05-1_N-QDs_bidirectional-waveguide_coherent-pulse.jl @@ -99,11 +99,11 @@ nothing # hide ba = NLevelBasis(2) b = tensor([ba for i = 1:N]...) -H_QO = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_R_QO = translate_qo(L_R, b; parameter = dict_p, time_parameter = dict_p_t) -L_L_QO = translate_qo(L_L, b; parameter = dict_p, time_parameter = dict_p_t) +H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) +L_R_QO = to_numeric(L_R, b; parameter = dict_p, time_parameter = dict_p_t) +L_L_QO = to_numeric(L_L, b; parameter = dict_p, time_parameter = dict_p_t) -σ_qo(α, i, j) = translate_qo(σ(α, i, j), b) +σ_qo(α, i, j) = to_numeric(σ(α, i, j), b) J_add = [√(γ_add[i])*σ_qo(i, 1, 2) for i = 1:N] function input_output(t, ρ) diff --git a/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl b/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl index 6a7c05f..f872698 100644 --- a/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl +++ b/examples/05-3_N-QDs_bidirectional-waveguide_feedback-reduction.jl @@ -89,11 +89,11 @@ nothing # hide ba = NLevelBasis(2) b = tensor([ba for _ = 1:N]...) -H_QO = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_R_QO = translate_qo(L_R, b; parameter = dict_p, time_parameter = dict_p_t) -L_L_QO = translate_qo(L_L, b; parameter = dict_p, time_parameter = dict_p_t) +H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) +L_R_QO = to_numeric(L_R, b; parameter = dict_p, time_parameter = dict_p_t) +L_L_QO = to_numeric(L_L, b; parameter = dict_p, time_parameter = dict_p_t) -σ_qo(α, i, j) = translate_qo(σ(α, i, j), b) +σ_qo(α, i, j) = to_numeric(σ(α, i, j), b) J_add = [√(γ_add[i]) * σ_qo(i, 1, 2) for i = 1:N] function input_output(t, ρ) diff --git a/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.jl b/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.jl index ac2cd7e..63a3713 100644 --- a/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.jl +++ b/examples/06-1_interaction-picture__PRA2023_107-013706_fig2.jl @@ -64,11 +64,11 @@ nothing # hide # -H_int_sym = simplify(substitute_operators(H_int_sym_, int_dict)) +H_int_sym = simplify(substitute(H_int_sym_, int_dict)) # -L_int_sym = simplify(substitute_operators(L, int_dict)) +L_int_sym = simplify(substitute(L, int_dict)) # The above Hamiltonian and Lindblad operator are the ones in the interaction picture of the virtual cavity interaction. # We define the numerical parameters of the system, calculate the solution for the coefficient matrix $M(t)$ and solve the time evolution of the system. @@ -109,8 +109,8 @@ ba = NLevelBasis(2) bv = FockBasis(5) b = bu ⊗ ba ⊗ bv -H_int_QO = translate_qo(H_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = translate_qo(L_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) +H_int_QO = to_numeric(H_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) +L_QO = to_numeric(L_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) function input_output(t, ρ) Ht = H_int_QO(t) diff --git a/examples/07-1_beamsplitter_loss__quantum-pulse.jl b/examples/07-1_beamsplitter_loss__quantum-pulse.jl index 20221a8..b87131d 100644 --- a/examples/07-1_beamsplitter_loss__quantum-pulse.jl +++ b/examples/07-1_beamsplitter_loss__quantum-pulse.jl @@ -92,8 +92,8 @@ au_qo = destroy(bu) ⊗ one(bv) av_qo = one(bu) ⊗ destroy(bv) ## translate to numeric operators -H_QO = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = [translate_qo(Li, b; parameter = dict_p, time_parameter = dict_p_t) for Li in L] +H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) +L_QO = [to_numeric(Li, b; parameter = dict_p, time_parameter = dict_p_t) for Li in L] function input_output(t, ρ) Ht = H_QO(t) diff --git a/examples/07-2_hong-ou-mandel__quantum-pulse.jl b/examples/07-2_hong-ou-mandel__quantum-pulse.jl index d52ae46..013ffa1 100644 --- a/examples/07-2_hong-ou-mandel__quantum-pulse.jl +++ b/examples/07-2_hong-ou-mandel__quantum-pulse.jl @@ -105,8 +105,8 @@ av1_qo = one(bu1) ⊗ one(bu2) ⊗ destroy(bv1) ⊗ one(bv2) av2_qo = one(bu1) ⊗ one(bu2) ⊗ one(bv1) ⊗ destroy(bv2) ## translate to numeric operators -H_QO = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) -L_QO = [translate_qo(Li, b; parameter = dict_p, time_parameter = dict_p_t) for Li in L] +H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) +L_QO = [to_numeric(Li, b; parameter = dict_p, time_parameter = dict_p_t) for Li in L] function input_output(t, ρ) Ht = H_QO(t) diff --git a/examples/08-1_pulse-delay__simple.jl b/examples/08-1_pulse-delay__simple.jl index 07b4cd7..78a1a59 100644 --- a/examples/08-1_pulse-delay__simple.jl +++ b/examples/08-1_pulse-delay__simple.jl @@ -82,8 +82,8 @@ bd = FockBasis(n) bv = FockBasis(n) b = bu ⊗ bd ⊗ bv -H_QO = translate_qo(H, b; time_parameter = dict_p_t) -L_QO = [translate_qo(L[i], b; time_parameter = dict_p_t) for i = 1:length(L)] +H_QO = to_numeric(H, b; time_parameter = dict_p_t) +L_QO = [to_numeric(L[i], b; time_parameter = dict_p_t) for i = 1:length(L)] function input_output(t, ρ) Ht = H_QO(t) @@ -98,9 +98,9 @@ nothing # hide ψ0 = fockstate(bu, n) ⊗ fockstate(bd, 0) ⊗ fockstate(bv, 0) t_, ρt = timeevolution.master_dynamic(T, ψ0, input_output) -au_qo = translate_qo(au, b) -ad_qo = translate_qo(ad, b) -av_qo = translate_qo(av, b) +au_qo = to_numeric(au, b) +ad_qo = to_numeric(ad, b) +av_qo = to_numeric(av, b) nu = real.(expect(dagger(au_qo)*au_qo, ρt)) nd = real.(expect(dagger(ad_qo)*ad_qo, ρt)) @@ -146,11 +146,11 @@ nothing # hide # -H_int_sym = simplify(substitute_operators(H_int_sym_, int_dict)) +H_int_sym = simplify(substitute(H_int_sym_, int_dict)) # -L_int_sym = simplify.(substitute_operators.(L, Ref(int_dict))) +L_int_sym = simplify.(substitute.(L, Ref(int_dict))) L_int_sym[1] # @@ -198,13 +198,13 @@ operators = Dict( [au, au', ad, ad', av, av'] .=> [au_int, au_int', ad_int, ad_int', av_int, av_int'], ) -H_int_QO = translate_qo(H_int_sym, b; time_parameter = dict_p_t_int, operators) +H_int_QO = to_numeric(H_int_sym, b; time_parameter = dict_p_t_int, operators) L_int_QO = [ - translate_qo(L_int_sym[i], b; time_parameter = dict_p_t_int, operators) for + to_numeric(L_int_sym[i], b; time_parameter = dict_p_t_int, operators) for i = 1:length(L_int_sym) ] -## H_int_QO = translate_qo(H_int_sym, b; time_parameter=dict_p_t_int) # hide -## L_int_QO = [translate_qo(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; 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 nothing # hide # diff --git a/examples/10-1_SUPER_excitation.jl b/examples/10-1_SUPER_excitation.jl index 398421a..3a11a44 100644 --- a/examples/10-1_SUPER_excitation.jl +++ b/examples/10-1_SUPER_excitation.jl @@ -210,9 +210,9 @@ la = length(a0_ls) a_int_ls = [sum(M(i, j)*a0_ls[j] for j = 1:la) for i = 1:la] int_dict = Dict(a0_ls .=> a_int_ls) -H_int = substitute_operators(H_int_, int_dict) -L_int = simplify(substitute_operators(Lcas, int_dict)) -Ld_int = simplify(substitute_operators(Lcasd, int_dict)) +H_int = substitute(H_int_, int_dict) +L_int = simplify(substitute(Lcas, int_dict)) +Ld_int = simplify(substitute(Lcasd, int_dict)) nothing # hide ## Coefficient matrix M @@ -300,8 +300,8 @@ g_t_ls = [gu2_t_, gu1_t_, gv1_t_, gv2_t_] M_t_ls = [t -> M_t(t)[i, j] for i = 1:la for j = 1:la] dict_fock = Dict([g_ls; M_ls] .=> [g_t_ls; M_t_ls]) -H_int_fock = translate_qo(H_int, b; parameter = Dict(γ=>γ_), time_parameter = dict_fock) -L_int_fock = translate_qo(L_int, b; parameter = Dict(γ=>γ_), time_parameter = dict_fock) +H_int_fock = to_numeric(H_int, b; parameter = Dict(γ=>γ_), time_parameter = dict_fock) +L_int_fock = to_numeric(L_int, b; parameter = Dict(γ=>γ_), time_parameter = dict_fock) ## To solve the dynamics, we create the time-dependent function for the open quantum system and define the initial state. 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 ce3f080..6d43b19 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 @@ -87,8 +87,8 @@ b1 = bc dict_p = Dict([γ, K, Δ] .=> [γ_, K_, Δ_]) dict_p_t = Dict(p => p_t) -H_QO = translate_qo(H, b1; parameter = dict_p, time_parameter = dict_p_t) -L_QO = translate_qo(L, b1; parameter = dict_p, time_parameter = dict_p_t) +H_QO = to_numeric(H, b1; parameter = dict_p, time_parameter = dict_p_t) +L_QO = to_numeric(L, b1; parameter = dict_p, time_parameter = dict_p_t) function input_output_1(t, ρ) Ht = H_QO(t) @@ -134,8 +134,8 @@ dict_p_t_2 = Dict(p => p_t, g_v => gv_t) bv = FockBasis(n_cut) b2 = bc ⊗ bv -H_QO_2 = translate_qo(H_2, b2; parameter = dict_p, time_parameter = dict_p_t_2) -L_QO_2 = translate_qo(L_2, b2; parameter = dict_p, time_parameter = dict_p_t_2) +H_QO_2 = to_numeric(H_2, b2; parameter = dict_p, time_parameter = dict_p_t_2) +L_QO_2 = to_numeric(L_2, b2; parameter = dict_p, time_parameter = dict_p_t_2) function input_output_2(t, ρ) Ht = H_QO_2(t) 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 3236bb0..849a918 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 @@ -95,8 +95,8 @@ int_dict = Dict(a_ls .=> a_int_ls) # int_dict = Dict([a1_ls; a2_ls] .=> [a1_int_ls; a2_int_ls]) ## substitute interaction picture operators -H_int = simplify(substitute_operators(H_int_, int_dict)) -L_int = [simplify(substitute_operators(L_, int_dict)) for L_ in L] +H_int = simplify(substitute(H_int_, int_dict)) +L_int = [simplify(substitute(L_, int_dict)) for L_ in L] ## Pulse parameters n = 5#9 @@ -148,13 +148,13 @@ b = bu ⊗ bd1 ⊗ bd2 ⊗ bs dict_p_Δ(Δn) = Dict([γ, Δ, r, t] .=> [γ_, Δn, rn, tn]) -H_QO_Δ(Δn) = translate_qo(H_int, b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) +H_QO_Δ(Δn) = to_numeric(H_int, b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) L_QO_Δ(Δn) = [ - translate_qo(L_int[i], b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) + to_numeric(L_int[i], b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) for i = 1:length(L) ] -# H_QO_Δ(Δn) = translate_qo(H, b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) -# L_QO_Δ(Δn) = [translate_qo(L[i], b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) for i=1:length(L)] +# H_QO_Δ(Δn) = to_numeric(H, b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) +# L_QO_Δ(Δn) = [to_numeric(L[i], b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) for i=1:length(L)] nothing # hide # @@ -162,10 +162,10 @@ nothing # hide # -σ22_QO = translate_qo(σ(2, 2), b) -au_QO = translate_qo(au, b) -ad1_QO = translate_qo(ad1, b) -ad2_QO = translate_qo(ad2, b) +σ22_QO = to_numeric(σ(2, 2), b) +au_QO = to_numeric(au, b) +ad1_QO = to_numeric(ad1, b) +ad2_QO = to_numeric(ad2, b) Δn = 0 H_QO = H_QO_Δ(Δn) diff --git a/examples/drafts/10-1_SUPER_excitation.jl b/examples/drafts/10-1_SUPER_excitation.jl index 02b5ec0..6e4a7c6 100644 --- a/examples/drafts/10-1_SUPER_excitation.jl +++ b/examples/drafts/10-1_SUPER_excitation.jl @@ -305,9 +305,9 @@ a_c_int_ls = [sum(conj(M(i, j))*a0_ls[j]' for j = 1:la) for i = 1:la] # # TODO int_dict = Dict([a0_ls; adjoint.(a0_ls)] .=> [a_int_ls; a_c_int_ls]) -H_int_sym = substitute_operators(H_int_sym_, int_dict) -L_int_sym = simplify(substitute_operators(Lcas, int_dict)) -Ld_int_sym = simplify(substitute_operators(Lcasd, int_dict)) +H_int_sym = substitute(H_int_sym_, int_dict) +L_int_sym = simplify(substitute(Lcas, int_dict)) +Ld_int_sym = simplify(substitute(Lcasd, int_dict)) # M-matrix Mat = Matrix{Any}(undef, la, la) diff --git a/src/QuantumInputOutput.jl b/src/QuantumInputOutput.jl index 9ee3a8f..600b1c0 100644 --- a/src/QuantumInputOutput.jl +++ b/src/QuantumInputOutput.jl @@ -1,10 +1,10 @@ module QuantumInputOutput -using SecondQuantizedAlgebra: SecondQuantizedAlgebra, QSym, to_numeric -using QuantumOpticsBase: QuantumOpticsBase, expect, basis, dagger, sparse +using SecondQuantizedAlgebra: SecondQuantizedAlgebra, to_numeric +using QuantumOpticsBase: QuantumOpticsBase, expect, dagger using QuantumOptics: QuantumOptics, timeevolution using SymbolicUtils: SymbolicUtils, substitute, BasicSymbolic, simplify -using Symbolics: Symbolics, build_function, Num +using Symbolics: Symbolics using SpecialFunctions: erf using DataInterpolations: LinearInterpolation, ExtrapolationType using NumericalIntegration: cumul_integrate @@ -27,7 +27,8 @@ export SLH, concatenate, feedback, # Translation - translate_qo, + to_numeric, + substitute, # Pulse coupling coupling_input, coupling_output, @@ -41,13 +42,9 @@ export SLH, solve_mode_evolution, solve_mode_evolution_symmetric, # Correlations - correlation_matrix, - # Operators - substitute_operators + correlation_matrix include("SLH.jl") -include("translate.jl") -include("utils.jl") include("pulses.jl") include("correlations.jl") include("interaction_picture.jl") diff --git a/src/SLH.jl b/src/SLH.jl index 9aaec7d..210c4fd 100644 --- a/src/SLH.jl +++ b/src/SLH.jl @@ -445,3 +445,23 @@ function _feedback_maps(n::Int, connections) end return mapped end + +# ────────────────────────────────────────────── +# Numeric translation of an SLH object +# ────────────────────────────────────────────── + +""" + to_numeric(G::SLH, b::QuantumOpticsBase.Basis; kwargs...) + +Translate the Hamiltonian and Lindblad operators of an SLH object `G` into numeric +[QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators on the basis `b`. +Returns the tuple `(H_QO, L_QO)`, where `L_QO` is a vector holding one translated operator +per jump operator in `lindblad(G)`. All keyword arguments (`parameter`, `time_parameter`, +`operators`, `adjoint_ops`, `op_type`) are forwarded to +[`SecondQuantizedAlgebra.to_numeric`](@ref). +""" +function SQA.to_numeric(G::SLH, b::QuantumOpticsBase.Basis; kwargs...) + H_QO = SQA.to_numeric(hamiltonian(G), b; kwargs...) + L_QO = [SQA.to_numeric(L_, b; kwargs...) for L_ in lindblad(G)] + return H_QO, L_QO +end diff --git a/src/translate.jl b/src/translate.jl deleted file mode 100644 index 12c1de4..0000000 --- a/src/translate.jl +++ /dev/null @@ -1,358 +0,0 @@ -# ── Numeric conversion helpers ── -# `to_numeric` requires the substitution dict to be `AbstractDict{<:QSym}`; the -# empty default `Dict()` is not, so route through the dict-less method instead. -_to_numeric(op, b, operators) = - isempty(operators) ? to_numeric(op, b) : to_numeric(op, b, operators) - -_translate_numeric( - op, - b; - parameter = Dict(), - operators = Dict{QSym,Any}(), - op_type = sparse, -) = op_type(_to_numeric(substitute(op, parameter), b, operators)) - -# Static numeric translation of a full operator sum. SQA's `to_numeric` reduces -# coefficients via `_to_complex`, which fails on constant-but-symbolic values such -# as `sqrt(1.0)` (the result of substituting a number into `√(symbol)`) because -# Symbolics 7 does not constant-fold them. Fold each term's coefficient through -# `_const_coeff` (compiled evaluation) before assembling the operator. -function _translate_numeric( - op::SQA.QAdd, - b; - parameter = Dict(), - operators = Dict{QSym,Any}(), - op_type = sparse, -) - op_ = substitute(op, parameter) - iszero(op_) && return op_type(0 * _one_op(b, operators)) - result = nothing - for (term, c_) in op_.arguments - c = _coeff_num(c_) - _coeff_is_const(c) || throw( - ArgumentError( - "cannot translate `$op` to a static operator: coefficient `$c` still " * - "depends on a free variable; supply it via `parameter` or `time_parameter`.", - ), - ) - contrib = _const_coeff(c) * _numeric_product(term.ops, b, operators, op_type) - result = result === nothing ? contrib : result + contrib - end - return op_type(result) -end - -_translate_numeric_raw(op, b; operators = Dict{QSym,Any}(), op_type = sparse) = - op_type(_to_numeric(op, b, operators)) - -# Numeric operator product of a `QTerm.ops` vector. -function _numeric_product(ops, b, operators, op_type) - isempty(ops) && return op_type(_one_op(b, operators)) - return op_type(prod(_to_numeric(o, b, operators) for o in ops)) -end - -_one_op(b, operators) = isempty(operators) ? one(b) : one(basis(first(values(operators)))) - -_translate_one(b, operators, op_type) = op_type(_one_op(b, operators)) - -# Build the operator-substitution dict, optionally adding adjoint entries. -# Always returns a `Dict{QSym,Any}` so it dispatches to `to_numeric(_, _, ::AbstractDict{<:QSym})`. -function _translate_operator_dict(operators, adjoint_ops::Bool) - out = Dict{QSym,Any}() - for (k, v) in operators - out[k] = v - end - if adjoint_ops - for (k, v) in operators - k_adj = Base.adjoint(k) - haskey(out, k_adj) || (out[k_adj] = Base.adjoint(v)) - end - end - return out -end - -# Complex-valued parameters (`Complex{Num}`) are stored in coefficients as separate -# `real`/`imag` parts, so a whole-variable substitution `gv => 0.0` never matches. -# Expand each complex key into independent real- and imaginary-part substitutions; -# the imaginary leg is dropped when it is a literal (e.g. a real variable that was -# promoted to `Complex{Num}` inside a mixed key vector). -function _expand_parameter(parameter) - isempty(parameter) && return parameter - out = Dict{Any,Any}() - for (k, v) in parameter - if k isa Complex - out[real(k)] = real(v) - ik = imag(k) - SymbolicUtils.unwrap(ik) isa Number || (out[ik] = imag(v)) - else - out[k] = v - end - end - return out -end - -# ── Time-parameter handling ── - -# Wrap plain-number values as constant functions; pass functions through. -function _normalize_time_parameter(time_parameter) - if isempty(time_parameter) - return time_parameter - end - KT = keytype(time_parameter) - out = Dict{KT,Any}() - for (k, v) in time_parameter - out[k] = v isa Number ? (t -> v) : v - end - return out -end - -# Reduce the time-parameter keys to the underlying symbolic base variables and -# the runtime value of each. Complex coefficients are stored in `real`/`imag` -# form in v0.6, so `build_function` only accepts bare variables as arguments; -# an expression key such as `conj(E)` is inverted to its base variable `E` with -# the value function conjugated accordingly. -function _time_basis(time_parameter) - basevars = Any[] - valuefuncs = Any[] - for (k, f) in time_parameter - uk = SymbolicUtils.unwrap(k) - if SymbolicUtils.issym(uk) - push!(basevars, k) - push!(valuefuncs, f) - continue - end - vs = collect(Symbolics.get_variables(k)) - length(vs) == 1 || error( - "time_parameter key `$k` must depend on exactly one variable, got $(length(vs)).", - ) - wv = Symbolics.wrap(vs[1]) - if isequal(uk, SymbolicUtils.unwrap(conj(wv))) - push!(basevars, wv) - push!(valuefuncs, t -> conj(f(t))) - else - error( - "unsupported time_parameter key `$k`; only a bare variable `v` or `conj(v)` are supported.", - ) - end - end - return basevars, Tuple(valuefuncs) -end - -# SQA v0.7 stores each `QTerm` prefactor as a `Coeff` (native/poly/symbolic forms); -# lower it to a `Complex{Num}` at the symbolic boundary so the helpers below operate -# on plain Symbolics expressions. A `Complex{Num}` passes through unchanged. -_coeff_num(c) = SQA.to_num(c) -_coeff_num(c::Complex{Num}) = c - -# A coefficient with no free variables is concrete. -_coeff_is_const(c::Complex{Num}) = - isempty(Symbolics.get_variables(real(c))) && isempty(Symbolics.get_variables(imag(c))) - -# Reduce a concrete `Complex{Num}` to a plain Julia number. Coefficients that are -# already numeric take a fast path; constant symbolic expressions (e.g. `exp(0.5im)` -# produced by substituting a numeric value into `exp(im*ϕ)`) are compiled and evaluated. -function _const_coeff(c::Complex{Num}) - re = Symbolics.value(SymbolicUtils.unwrap(real(c))) - im = Symbolics.value(SymbolicUtils.unwrap(imag(c))) - if re isa Number && im isa Number - return iszero(im) ? re : complex(re, im) - end - # `build_function` returns either a single function or an (out-of-place, in-place) - # tuple; take the out-of-place callable and evaluate the (variable-free) expression. - f = build_function(c; expression = Val(false)) - g = f isa Tuple ? first(f) : f - return g() -end - -_as_cnum(x::Complex) = Complex{Num}(Num(real(x)), Num(imag(x))) -_as_cnum(x) = Complex{Num}(Num(x), Num(false)) - -# Compile a `Complex{Num}` coefficient into a runtime function of the time-parameter -# base variables. A naive `build_function(c, vars...)` emits a `complex(re, im)` call, -# which requires both arguments to be `<:Real` and therefore throws `MethodError: -# complex(::Int64, ::ComplexF64)` when a base variable is fed a complex-valued time -# function. Following QuantumCumulants (`_im_form`/`_qc_maketerm`), we instead compile -# the real and imaginary parts separately and recombine them with `im` in Julia, so the -# generated code uses ordinary `+`/`*` and stays valid for real- and complex-valued -# arguments alike. This is behaviour-preserving for the previously-working cases. -function _compile_coeff(c::Complex{Num}, vars...) - f_re = build_function(real(c), vars...; expression = Val(false)) - f_im = build_function(imag(c), vars...; expression = Val(false)) - g_re = f_re isa Tuple ? first(f_re) : f_re - g_im = f_im isa Tuple ? first(f_im) : f_im - return (vals...) -> g_re(vals...) + im * g_im(vals...) -end - -# ── Per-term translation ── -# Translate a single `(ops, coefficient)` term into a time-dependent function -# `t -> op`. A concrete coefficient yields a constant function; a symbolic -# coefficient is compiled against the time-parameter base variables. -function _translate_term(ops, c::Complex{Num}, b, time_parameter, operators, op_type) - prodop = _numeric_product(ops, b, operators, op_type) - if _coeff_is_const(c) - op = _const_coeff(c) * prodop - return t -> op - end - basevars, valuefuncs = _time_basis(time_parameter) - pref = _compile_coeff(c, basevars...) - return t -> pref(map(g -> g(t), valuefuncs)...) * prodop -end - -""" - translate_qo(op, b::QuantumOpticsBase.Basis; parameter=Dict(), time_parameter=Dict(), - operators=Dict(), adjoint_ops=true, op_type=sparse) - -Translate a symbolic operator `op` into a numeric QuantumOptics.jl operator with the corresponding basis `b`. -The dictionary `parameter` substitutes symbolic parameters with numbers. Time-dependent functions can be provided -with the dictionary `time_parameter`. -If `time_parameter` is non-empty, the result is a time-dependent function `t -> op(t)`. -The operator type which should be returned can be set with the kwarg `op_type=sparse` and -a list of user-defined operators (e.g. on a different basis than `b`) can be provided with the -dictionary `operators=Dict()`. These operators will then be used to replace the symbolic expressions. -If `adjoint_ops=true`, adjoint entries are added automatically for operators missing from the -dictionary, e.g. `a => A` also provides `a' => A'`. -""" -function translate_qo( - op, - b::QuantumOpticsBase.Basis; - parameter = Dict(), - time_parameter = Dict(), - operators = Dict(), - adjoint_ops = true, - op_type = sparse, -) - param = _expand_parameter(parameter) - tp = _normalize_time_parameter(time_parameter) - operators_ = _translate_operator_dict(operators, adjoint_ops) - return _translate_qo( - op, - b; - parameter = param, - time_parameter = tp, - operators = operators_, - op_type, - ) -end - -# ── Internal dispatch: time_parameter already normalized ── - -function _translate_qo( - op::SQA.QAdd, - b::QuantumOpticsBase.Basis; - parameter = Dict(), - time_parameter = Dict(), - operators = Dict{QSym,Any}(), - op_type = sparse, -) - if isempty(time_parameter) - return _translate_numeric(op, b; parameter, operators, op_type) - end - - op_ = substitute(op, parameter) - iszero(op_) && return t -> op_type(0 * one(b)) - - pairs = collect(op_.arguments) - if length(pairs) == 1 - (term, c) = pairs[1] - return _translate_term( - term.ops, - _coeff_num(c), - b, - time_parameter, - operators, - op_type, - ) - end - - # Multi-term: combine via FunctionWrappers with a common concrete type so the - # returned closure is type-stable. Inner products use `sparse` for a uniform type. - first_res = _translate_term( - pairs[1].first.ops, - _coeff_num(pairs[1].second), - b, - time_parameter, - operators, - sparse, - ) - OpType = typeof(first_res(0.0)) - FW = FunctionWrapper{OpType,Tuple{Float64}} - - args_wrapped = ntuple(length(pairs)) do k - res = - k == 1 ? first_res : - _translate_term( - pairs[k].first.ops, - _coeff_num(pairs[k].second), - b, - time_parameter, - operators, - sparse, - ) - FW(res) - end - - return t -> begin - result = args_wrapped[1](t) - for i = 2:length(args_wrapped) - result = result + args_wrapped[i](t) - end - result - end -end - -function _translate_qo( - op::QSym, - b::QuantumOpticsBase.Basis; - parameter = Dict(), - time_parameter = Dict(), - operators = Dict{QSym,Any}(), - op_type = sparse, -) - if isempty(time_parameter) - return _translate_numeric(op, b; parameter, operators, op_type) - end - output_op = _translate_numeric_raw(op, b; operators, op_type) - return t -> output_op -end - -function _translate_qo( - arg_c_, - b::QuantumOpticsBase.Basis; - parameter = Dict(), - time_parameter = Dict(), - operators = Dict{QSym,Any}(), - op_type = sparse, -) - arg_c = substitute(arg_c_, parameter) - one_b = _translate_one(b, operators, op_type) - c = _as_cnum(arg_c) - - if isempty(time_parameter) - _coeff_is_const(c) || throw( - ArgumentError( - "cannot translate symbolic scalar `$arg_c` without a value: supply it via " * - "`parameter` (numeric) or `time_parameter` (time-dependent).", - ), - ) - return _const_coeff(c) * one_b - end - - if _coeff_is_const(c) - val = _const_coeff(c) - return t -> val * one_b - end - basevars, valuefuncs = _time_basis(time_parameter) - pref = _compile_coeff(c, basevars...) - return t -> pref(map(g -> g(t), valuefuncs)...) * one_b -end - -function translate_qo(ops::Vector, b::QuantumOpticsBase.Basis; kwargs...) - return [translate_qo(op, b; kwargs...) for op in ops] -end - -function translate_qo(G::SLH, b::QuantumOpticsBase.Basis; kwargs...) - L = lindblad(G) - H = hamiltonian(G) - H_QO = translate_qo(H, b; kwargs...) - L_QO = [translate_qo(L_, b; kwargs...) for L_ in L] - return H_QO, L_QO -end diff --git a/src/utils.jl b/src/utils.jl deleted file mode 100644 index 9295f52..0000000 --- a/src/utils.jl +++ /dev/null @@ -1,39 +0,0 @@ -""" - substitute_operators(op, dict::Dict; replace_adjoint=true) - -Like `substitute(op, dict::Dict)` but with special handling for operator products and sums. -This is needed if an operator is substituted by a product or sum of operators, e.g. - ``a_1 \\rightarrow g_2 a_2 + g_3 a_3`` - -Plain `substitute` recurses into the replacement, so substituting an operator by an -expression that still contains that operator loops forever. This routine instead replaces -each leaf operator with a single dictionary lookup and rebuilds the product/sum. - -If `replace_adjoint=true`, the dictionary is extended with adjoint substitutions -for all key/value pairs, i.e. `adjoint(key) => adjoint(value)`. -""" -function substitute_operators(op, dict::Dict; replace_adjoint = true) - dict_ = replace_adjoint ? _extend_with_adjoint(dict) : dict - return substitute(op, dict_) -end - -function substitute_operators(op::SQA.QSym, dict::Dict; replace_adjoint = true) - dict_ = replace_adjoint ? _extend_with_adjoint(dict) : dict - return get(dict_, op, op) -end - -function substitute_operators(op::SQA.QAdd, dict::Dict; replace_adjoint = true) - dict_ = replace_adjoint ? _extend_with_adjoint(dict) : dict - iszero(op) && return op - return sum(_coeff_num(c) * _replace_ops(term.ops, dict_) for (term, c) in op) -end - -# Single-level operator replacement: each leaf is replaced by a direct dictionary -# lookup (never re-substituted), avoiding the infinite recursion of `substitute`. -_replace_ops(ops, dict_) = isempty(ops) ? 1 : prod(get(dict_, o, o) for o in ops) - -function _extend_with_adjoint(dict::Dict) - pairs_ = collect(dict) - adj_pairs = [Base.adjoint(k) => Base.adjoint(v) for (k, v) in pairs_] - return Dict(vcat(pairs_, adj_pairs)) -end diff --git a/test/test_compare_example_05_1_05_2.jl b/test/test_compare_example_05_1_05_2.jl index f554fb0..bc141d5 100644 --- a/test/test_compare_example_05_1_05_2.jl +++ b/test/test_compare_example_05_1_05_2.jl @@ -58,11 +58,11 @@ using Test ba = NLevelBasis(2) b = tensor([ba for _ = 1:N]...) - H_QO = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) - L_R_QO = translate_qo(L_R, b; parameter = dict_p, time_parameter = dict_p_t) - L_L_QO = translate_qo(L_L, b; parameter = dict_p, time_parameter = dict_p_t) + H_QO = to_numeric(H, b; parameter = dict_p, time_parameter = dict_p_t) + L_R_QO = to_numeric(L_R, b; parameter = dict_p, time_parameter = dict_p_t) + L_L_QO = to_numeric(L_L, b; parameter = dict_p, time_parameter = dict_p_t) - σ_qo(α, i, j) = translate_qo(σ(α, i, j), b) + σ_qo(α, i, j) = to_numeric(σ(α, i, j), b) J_add = [√(γ_add[i]) * σ_qo(i, 1, 2) for i = 1:N] function input_output_1(t, ρ) diff --git a/test/test_example_cavity_scattering.jl b/test/test_example_cavity_scattering.jl index 4bcc30e..d68037a 100644 --- a/test/test_example_cavity_scattering.jl +++ b/test/test_example_cavity_scattering.jl @@ -47,8 +47,8 @@ using Test bv1 = FockBasis(2) b = bu1 ⊗ bc1 ⊗ bv1 - H_QO = translate_qo(H, b; parameter = dict_p, time_parameter = dict_p_t) - L_QO = translate_qo(L, b; parameter = dict_p, time_parameter = dict_p_t) + 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) function input_output_1(t, ρ) Ht = H_QO(t) @@ -59,9 +59,9 @@ using Test ψ0 = fockstate(bu1, 1) ⊗ fockstate(bc1, 0) ⊗ fockstate(bv1, 0) t_, ρt = timeevolution.master_dynamic(T, ψ0, input_output_1) - au_qo = translate_qo(au, b) - c_qo = translate_qo(c, b) - av_qo = translate_qo(av, b) + au_qo = to_numeric(au, b) + c_qo = to_numeric(c, b) + av_qo = to_numeric(av, b) n_u_t = real.(expect(au_qo' * au_qo, ρt)) n_c_t = real.(expect(c_qo' * c_qo, ρt)) @@ -90,8 +90,8 @@ using Test dict_p_t_2 = Dict([gu, gv] .=> [gu_t, gv_t]); - H_QO_2 = translate_qo(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) - L_QO_2 = translate_qo(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) + H_QO_2 = to_numeric(H, b; parameter = dict_p_2, time_parameter = dict_p_t_2) + L_QO_2 = to_numeric(L, b; parameter = dict_p_2, time_parameter = dict_p_t_2) function input_output_2(t, ρ) H = H_QO_2(t) J = [L_QO_2(t)] diff --git a/test/test_interaction_picture.jl b/test/test_interaction_picture.jl index ac375d4..39f88ff 100644 --- a/test/test_interaction_picture.jl +++ b/test/test_interaction_picture.jl @@ -50,8 +50,8 @@ using Test a_int_ls = [sum(M(i, j) * a0_ls[j] for j = 1:la) for i = 1:la] int_dict = Dict(a0_ls .=> a_int_ls) - H_int_sym = simplify(substitute_operators(H_int_sym_, int_dict)) - L_int_sym = simplify(substitute_operators(L, int_dict)) + H_int_sym = simplify(substitute(H_int_sym_, int_dict)) + L_int_sym = simplify(substitute(L, int_dict)) # Virtual-cavity couplings gu_t = coupling_input(u, T) @@ -90,8 +90,8 @@ using Test p_t_num = [gu_t, gv_t, M_t_ls...] dict_p_t = Dict(p_t_sym .=> p_t_num) - H_int_QO = translate_qo(H_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) - L_QO = translate_qo(L_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) + H_int_QO = to_numeric(H_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) + L_QO = to_numeric(L_int_sym, b; parameter = dict_p, time_parameter = dict_p_t) function input_output_I(t, ρ) Ht = H_int_QO(t) @@ -114,8 +114,8 @@ using Test dict_p_s = Dict(γ_sym => γ) dict_p_t_s = Dict(gu_sym => gu_t, gv_sym => gv_t) - H_QO = translate_qo(H, b; parameter = dict_p_s, time_parameter = dict_p_t_s) - L_QO_S = translate_qo(L, b; parameter = dict_p_s, time_parameter = dict_p_t_s) + H_QO = to_numeric(H, b; parameter = dict_p_s, time_parameter = dict_p_t_s) + L_QO_S = to_numeric(L, b; parameter = dict_p_s, time_parameter = dict_p_t_s) function input_output_S(t, ρ) Ht = H_QO(t) @@ -132,9 +132,9 @@ using Test γ_small = 1e-4 dict_p_small = Dict(γ_sym => γ_small) H_int_QO_small = - translate_qo(H_int_sym, b; parameter = dict_p_small, time_parameter = dict_p_t) + to_numeric(H_int_sym, b; parameter = dict_p_small, time_parameter = dict_p_t) L_QO_small = - translate_qo(L_int_sym, b; parameter = dict_p_small, time_parameter = dict_p_t) + to_numeric(L_int_sym, b; parameter = dict_p_small, time_parameter = dict_p_t) function input_output_I_small(t, ρ) Ht = H_int_QO_small(t) diff --git a/test/test_translate.jl b/test/test_translate.jl index 9173b78..3a5a5b7 100644 --- a/test/test_translate.jl +++ b/test/test_translate.jl @@ -37,19 +37,19 @@ using Test @testset "kwarg_operators" begin @test isequal( - translate_qo(Δ, bc1; parameter = dict_p1, operators = ops_dict), + to_numeric(Δ, bc1; parameter = dict_p1, operators = ops_dict), one(bc1)*Δn, ) @test isequal( - translate_qo(2, bc1; parameter = dict_p1, operators = ops_dict), + to_numeric(2, bc1; parameter = dict_p1, operators = ops_dict), one(bc1)*2, ) - @test isequal(translate_qo(a, bc1; parameter = dict_p1, operators = ops_dict), a_QO) + @test isequal(to_numeric(a, bc1; parameter = dict_p1, operators = ops_dict), a_QO) @test isequal( - translate_qo(a*3, bc1; parameter = dict_p1, operators = ops_dict), + to_numeric(a*3, bc1; parameter = dict_p1, operators = ops_dict), a_QO*3, ) - F1 = translate_qo( + F1 = to_numeric( a*3, bc1; parameter = dict_p1, @@ -58,7 +58,7 @@ using Test ) @test isa(F1, Function) @test isequal(F1(0.1), a_QO*3) - F2 = translate_qo( + F2 = to_numeric( a*E, bc1; parameter = dict_p1, @@ -71,7 +71,7 @@ using Test abs.( ( dense( - translate_qo(Δ*a'a, bc1; parameter = dict_p1, operators = ops_dict), + to_numeric(Δ*a'a, bc1; parameter = dict_p1, operators = ops_dict), ) - dense(Δn*dagger(a_QO)*a_QO) ).data, ), @@ -88,7 +88,7 @@ using Test @test isequal( a_QO2'σ_QO(1, 2), dense( - translate_qo( + to_numeric( a'σ(1, 2), b; operators = Dict([a, σ(1, 2)] .=> [a_QO2, σ_QO(1, 2)]), @@ -97,17 +97,19 @@ using Test ), ) @test isequal(a_QO2, dense(to_numeric(a, b))) - @test isequal(translate_qo(Δ, b; parameter = dict_p1), one(b)*Δn) - F3 = translate_qo(Δ, b; parameter = dict_p1, time_parameter = dict_p_t2) - @test sum(abs.((F3(4) - one(b)*Δn).data)) < 1e-8 - F4 = translate_qo( + @test sum( + 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 + 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 - F5 = translate_qo( + F5 = to_numeric( a*conj(E) + Δ*σ(2, 2), b; parameter = dict_p1, @@ -115,11 +117,11 @@ using Test ) @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 - F5_ = translate_qo(a*conj(E), b; parameter = dict_p1, time_parameter = dict_p_t2) - @test sum(abs.((F5_(0.2) - dense(a_QO2*E_t_c(0.2))).data)) < 1e-8 - F6 = translate_qo(conj(E), b; parameter = dict_p1, time_parameter = dict_p_t2) - @test sum(abs.((F6(0.2) - dense(E_t_c(0.2)*one(b))).data)) < 1e-8 - F7 = translate_qo( + 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 + 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 + F7 = to_numeric( conj(E) + Δ*σ(2, 2), b; parameter = dict_p1, @@ -127,15 +129,15 @@ using Test ) @test sum(abs.((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 translate_qo(conj(E), b; parameter = dict_p1) - F8 = translate_qo(E^2, b; parameter = dict_p1, time_parameter = dict_p_t2) - @test sum(abs.((F8(0.2) - dense(E_t(0.2)^2*one(b))).data)) < 1e-8 + @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 @testset "time_parameter_normalization" begin dict_p_t_num = Dict([E] .=> [2.5]) - F_num = translate_qo(a*E, b; parameter = dict_p1, time_parameter = dict_p_t_num) - @test sum(abs.((F_num(0.2) - dense(a_QO2 * 2.5)).data)) < 1e-8 + 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 end @testset "multiple_time_prefactors" begin @@ -143,19 +145,19 @@ using Test E2_t(t) = 0.7 - 0.1im + 0.2t E2_t_c(t) = conj(E2_t(t)) dict_p_t_multi = Dict([E1, conj(E2)] .=> [E1_t, E2_t_c]) - F_multi = translate_qo( + F_multi = to_numeric( a * E1 * conj(E2), b; parameter = dict_p1, time_parameter = dict_p_t_multi, ) expected = dense(a_QO2 * E1_t(0.4) * E2_t_c(0.4)) - @test sum(abs.((F_multi(0.4) - expected).data)) < 1e-8 + @test sum(abs.((dense(F_multi(0.4)) - expected).data)) < 1e-8 end - @testset "translate_qo vector overload" begin + @testset "to_numeric vector overload" begin ops = [a, Δ, a * E] - translated_ops = translate_qo( + translated_ops = to_numeric( ops, bc1; parameter = dict_p1, @@ -169,9 +171,9 @@ using Test @test isequal(translated_ops[3](0.3), a_QO * E_t(0.3)) end - @testset "translate_qo SLH overload" begin + @testset "to_numeric SLH overload" begin G_sym = SLH(1, [sqrt(κ_R) * a * E], Δ * a' * a) - H_QO, L_QO = translate_qo( + H_QO, L_QO = to_numeric( G_sym, bc1; parameter = dict_p1, @@ -186,7 +188,7 @@ using Test @test sum(abs.((L_QO[1](0.4) - dense(sqrt(κ_Rn) * a_QO * E_t(0.4))).data)) < 1e-8 end - @testset "substitute_operators_qmul" begin + @testset "substitute operators qmul" begin h2 = FockSpace(:h2) a = Destroy(h2, :a) @variables c1::Complex c2::Complex c3::Complex @@ -196,16 +198,16 @@ using Test dict_sub2 = Dict(a => a_2) x = a*a*c1 + a*c3 - y = substitute_operators(x, dict_sub) - y2 = substitute_operators(x, dict_sub2) + y = substitute(x, dict_sub) + y2 = substitute(x, dict_sub2) @test iszero(y - (a_1*a_1*c1 + a_1*c3)) @test iszero(y2 - (a_2*a_2*c1 + a_2*c3)) - @test substitute_operators(5, dict_sub) == 5 + @test substitute(5, dict_sub) == 5 end - @testset "substitute_operators adjoint in args_nc" begin + @testset "substitute adjoint in args_nc" begin h2 = FockSpace(:h2) a = Destroy(h2, :a) b = Destroy(h2, :b) @@ -217,7 +219,7 @@ using Test op = g * ad * a dict_sub = Dict(a => b) - result = substitute_operators(op, dict_sub) + result = substitute(op, dict_sub) expected = g * bd * b @test iszero(result - expected) end @@ -233,7 +235,7 @@ using Test a3_QO = destroy(b3) @variables gR::Real gR_t(t) = 1.0 + 2.0im - F = translate_qo(im * gR * a3, b3; time_parameter = Dict(gR => gR_t)) + 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 end From 3d9ceffe33bfce2080d295ad80bf66eefe42a3dc Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sat, 11 Jul 2026 20:50:16 +0200 Subject: [PATCH 2/4] format --- benchmarks/translation.jl | 7 +------ ...8-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl | 4 ++-- test/test_translate.jl | 7 +------ 3 files changed, 4 insertions(+), 14 deletions(-) diff --git a/benchmarks/translation.jl b/benchmarks/translation.jl index 7195870..4e5db47 100644 --- a/benchmarks/translation.jl +++ b/benchmarks/translation.jl @@ -81,12 +81,7 @@ function benchmark_translation!(SUITE) SUITE["Translation"]["time-dependent"]["3-cavity H+L"] = @benchmarkable begin to_numeric($H_sym, $b_cav; parameter = $dict_p_cav, time_parameter = $dict_p_t_cav) - to_numeric( - $L_sym, - $b_cav; - parameter = $dict_p_cav, - time_parameter = $dict_p_t_cav, - ) + to_numeric($L_sym, $b_cav; parameter = $dict_p_cav, time_parameter = $dict_p_t_cav) end ## --- Closure evaluation (the ODE hot loop) --- 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 849a918..38ddc01 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 @@ -150,8 +150,8 @@ b = bu ⊗ bd1 ⊗ bd2 ⊗ bs dict_p_Δ(Δn) = Dict([γ, Δ, r, t] .=> [γ_, Δn, rn, tn]) H_QO_Δ(Δn) = to_numeric(H_int, b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) L_QO_Δ(Δn) = [ - to_numeric(L_int[i], b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) - for i = 1:length(L) + to_numeric(L_int[i], b; parameter = dict_p_Δ(Δn), time_parameter = dict_p_t_int) for + i = 1:length(L) ] # H_QO_Δ(Δn) = to_numeric(H, b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) # L_QO_Δ(Δn) = [to_numeric(L[i], b; parameter=dict_p_Δ(Δn), time_parameter=dict_p_t) for i=1:length(L)] diff --git a/test/test_translate.jl b/test/test_translate.jl index 3a5a5b7..6b16ddc 100644 --- a/test/test_translate.jl +++ b/test/test_translate.jl @@ -121,12 +121,7 @@ using Test @test sum(abs.((dense(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 - F7 = to_numeric( - conj(E) + Δ*σ(2, 2), - b; - parameter = dict_p1, - time_parameter = dict_p_t2, - ) + 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 # a bare symbolic scalar without a numeric/time value cannot be translated @test_throws ArgumentError to_numeric(conj(E), b; parameter = dict_p1) From 6ab076cfdea477c1fac9604dcedb173a707453fd Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 12 Jul 2026 17:33:43 +0200 Subject: [PATCH 3/4] fix docs --- docs/src/api.md | 6 ------ 1 file changed, 6 deletions(-) diff --git a/docs/src/api.md b/docs/src/api.md index c8e295c..85650a1 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -93,9 +93,3 @@ solve_mode_evolution_symmetric ```@docs correlation_matrix ``` - -## [Utilities](@id API: Utilities) - -```@docs -substitute -``` From cfac0114ec971413036f85ded7a1197f58b8c4e3 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 12 Jul 2026 18:12:15 +0200 Subject: [PATCH 4/4] bump minor --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 60bab9f..22b3834 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "QuantumInputOutput" uuid = "18f9eda6-924c-47c7-a881-996695cfd7c6" -version = "0.2.0" +version = "0.3.0" [deps] DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"