feat: update to SQA v0.9

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"
@@ -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"

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").

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 ---

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

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,21 +80,16 @@ 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(
-            $L_sym,
-            $b_cav;
-            parameter = $dict_p_cav,
-            time_parameter = $dict_p_t_cav,
-        )
+        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)
     end
 
     ## --- Closure evaluation (the ODE hot loop) ---
 
     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]
 

docs/src/api.md

@@ -41,7 +41,7 @@ feedback
 ## [Translation](@id API: Translation)
 
 ```@docs
-translate_qo
+to_numeric
 ```
 
 ## [Pulses](@id API: Pulses)
@@ -93,9 +93,3 @@ solve_mode_evolution_symmetric
 ```@docs
 correlation_matrix
 ```
-
-## [Utilities](@id API: Utilities)
-
-```@docs
-substitute_operators
-```

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`.

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)

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]

examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl

@@ -62,16 +62,16 @@ 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)
 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)))

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.

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)]

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)]