refactor: More type-stable and concrete implementation

.other/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4.jl

@@ -59,8 +59,8 @@ L_atom = [0, 0, √(γ/2)*σ(1, 2), √(γ/2)*σ(1, 2)]
 G_atom = SLH(I4, L_atom, Δ*σ(2, 2))
 
 G_u_bs_d1_d2_atom = G_u_bs_d1_d2 ▷ G_atom
-H = get_hamiltonian(G_u_bs_d1_d2_atom)
-L = get_lindblad(G_u_bs_d1_d2_atom)
+H = hamiltonian(G_u_bs_d1_d2_atom)
+L = lindblad(G_u_bs_d1_d2_atom)
 
 #
 

.other/backups/08-2_pulse-propagation-delay__PRA2026_113-013730_fig4__bu.jl

@@ -60,8 +60,8 @@ L_atom = [0, 0, √(γ/2)*σ(1, 2), √(γ/2)*σ(1, 2)]
 G_atom = SLH(I4, L_atom, Δ*σ(2, 2))
 
 G_u_bs_d1_d2_atom = G_u_bs_d1_d2 ▷ G_atom
-H = get_hamiltonian(G_u_bs_d1_d2_atom)
-L = get_lindblad(G_u_bs_d1_d2_atom)
+H = hamiltonian(G_u_bs_d1_d2_atom)
+L = lindblad(G_u_bs_d1_d2_atom)
 
 #
 

Project.toml

@@ -5,20 +5,23 @@ version = "0.1.0"
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
+FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
 QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
 SecondQuantizedAlgebra = "f7aa4685-e143-4cb6-a7f3-073579757907"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [compat]
 # Aqua = "0.8"
 # CheckConcreteStructs = "0.1"
 DataInterpolations = "8"
+FunctionWrappers = "1"
 DiffEqNoiseProcess = "5.27.0"
 # ExplicitImports = "1.12"
 LinearAlgebra = "1.10"
@@ -30,6 +33,7 @@ QuantumOptics = "1"
 QuantumOpticsBase = "0.4, 0.5"
 SecondQuantizedAlgebra = "0.4.4"
 SpecialFunctions = "2"
+StaticArrays = "1"
 Symbolics = "6"
 SymbolicUtils = "3.6"
 Test = "1.10"

benchmarks/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuantumInputOutput = "18f9eda6-924c-47c7-a881-996695cfd7c6"
 QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"

benchmarks/correlations.jl

@@ -23,14 +23,14 @@ function benchmark_correlations!(SUITE)
     G_v = SLH(1, gv * av, 0)
     G_cas = ▷(G_u, G_c, G_v)
 
-    H_sym = get_hamiltonian(G_cas)
-    L_sym = get_lindblad(G_cas)[1]
+    H_sym = hamiltonian(G_cas)
+    L_sym = lindblad(G_cas)[1]
 
     γ_ = 1.0
     σ_pulse = 1 / γ_
     T = [0:0.002:1;] * 12σ_pulse
     u1(t) = 1 / (sqrt(σ_pulse) * π^(1 / 4)) * exp(-(t - 4σ_pulse)^2 / (2 * σ_pulse^2))
-    gu_t = u_to_gu(u1, T)
+    gu_t = coupling_input(u1, T)
 
     bu1 = FockBasis(1)
     bc1 = FockBasis(1)
@@ -61,7 +61,7 @@ function benchmark_correlations!(SUITE)
     SUITE["Correlations"]["two-time"] = BenchmarkGroup()
 
     SUITE["Correlations"]["two-time"]["single photon cavity"] =
-        @benchmarkable two_time_corr_matrix($T, $ρt, $input_output, $Ls)
+        @benchmarkable correlation_matrix($T, $ρt, $input_output, $Ls)
 
     return nothing
 end

benchmarks/interaction_picture.jl

@@ -1,6 +1,5 @@
 """
 Interaction picture benchmarks — computing coefficient matrices for mode transformations.
-Based on Christiansen et al., PRA 107, 013706 (2023).
 """
 function benchmark_interaction_picture!(SUITE)
     SUITE["Interaction Picture"] = BenchmarkGroup()
@@ -14,9 +13,9 @@ function benchmark_interaction_picture!(SUITE)
     T_end = 12.0
     T = [0:0.005:1;] * T_end
 
-    gu_t = u_to_gu(u, T)
-    gv_t = v_to_gv(u, T)
-    A_uv = interaction_picture_A_2modes(gu_t, gv_t)
+    gu_t = coupling_input(u, T)
+    gv_t = coupling_output(u, T)
+    A_uv = coupling_matrix((gu_t, gv_t))
 
     ## --- Coupling matrix A(t) evaluation (called every ODE step) ---
 
@@ -28,11 +27,11 @@ function benchmark_interaction_picture!(SUITE)
     SUITE["Interaction Picture"]["coupling matrix evaluation"]["2 modes"] =
         @benchmarkable $A_uv($t_mid)
 
-    # 4-mode A(t) — more realistic for multi-port systems
+    # 4-mode A(t)
     u2(t) = 1 / (sqrt(τ) * π^(1 / 4)) * exp(-0.5 * ((t - t_p * 1.5) / τ)^2)
-    g3_t = u_to_gu(u2, T)
-    g4_t = v_to_gv(u2, T)
-    A_4m = interaction_picture_A_4modes(gu_t, gv_t, g3_t, g4_t)
+    g3_t = coupling_input(u2, T)
+    g4_t = coupling_output(u2, T)
+    A_4m = coupling_matrix((gu_t, gv_t, g3_t, g4_t))
 
     SUITE["Interaction Picture"]["coupling matrix evaluation"]["4 modes"] =
         @benchmarkable $A_4m($t_mid)
@@ -42,10 +41,10 @@ function benchmark_interaction_picture!(SUITE)
     SUITE["Interaction Picture"]["coefficient matrix M"] = BenchmarkGroup()
 
     SUITE["Interaction Picture"]["coefficient matrix M"]["numerical (ODE)"] =
-        @benchmarkable interaction_picture_M($A_uv, $T)
+        @benchmarkable solve_mode_evolution($A_uv, $T)
 
     SUITE["Interaction Picture"]["coefficient matrix M"]["analytical (2 equal modes)"] =
-        @benchmarkable interaction_picture_M_2modes_equal($u, $T)
+        @benchmarkable solve_mode_evolution_symmetric($u, $T)
 
     ## --- Symbolic operator substitution ---
 
@@ -65,9 +64,9 @@ function benchmark_interaction_picture!(SUITE)
     G_v = SLH(1, gv_sym' * av_sym, 0)
     G_cas = ▷(G_u, G_s, G_v)
 
-    H = get_hamiltonian(G_cas)
-    L = get_lindblad(G_cas)[1]
-    H_uv = get_hamiltonian(▷(G_u, G_v))
+    H = hamiltonian(G_cas)
+    L = lindblad(G_cas)[1]
+    H_uv = hamiltonian(▷(G_u, G_v))
     H_int_ = simplify(H - H_uv)
 
     M_sym(i, j) = cnumber("M_{$(i)$(j)}")

benchmarks/pulse_couplings.jl

@@ -17,9 +17,18 @@ function benchmark_pulse_couplings!(SUITE)
 
     SUITE["Pulse Couplings"]["single-pulse"] = BenchmarkGroup()
 
-    SUITE["Pulse Couplings"]["single-pulse"]["input"] = @benchmarkable u_to_gu($u, $T)
+    SUITE["Pulse Couplings"]["single-pulse"]["input"] =
+        @benchmarkable coupling_input($u, $T)
 
-    SUITE["Pulse Couplings"]["single-pulse"]["output"] = @benchmarkable v_to_gv($v, $T)
+    SUITE["Pulse Couplings"]["single-pulse"]["output"] =
+        @benchmarkable coupling_output($v, $T)
+
+    # Gaussian analytical
+    SUITE["Pulse Couplings"]["single-pulse"]["input Gaussian"] =
+        @benchmarkable coupling_input(Gaussian($τ, $σ; δ = $δ))
+
+    SUITE["Pulse Couplings"]["single-pulse"]["output Gaussian"] =
+        @benchmarkable coupling_output(Gaussian($τ, $σ; δ = $δ))
 
     ## --- Effective multi-pulse couplings ---
 
@@ -31,10 +40,10 @@ function benchmark_pulse_couplings!(SUITE)
     u2(t) = 1 / (sqrt(σ) * π^(1 / 4)) * exp(-0.5 * ((t - (τ + 2σ)) / σ)^2)
 
     SUITE["Pulse Couplings"]["multi-pulse"]["output 2 modes"] =
-        @benchmarkable v_eff([$v1, $v2], $T, 2)
+        @benchmarkable effective_output_mode([$v1, $v2], $T, 2)
 
     SUITE["Pulse Couplings"]["multi-pulse"]["input 2 modes"] =
-        @benchmarkable u_eff([$u1, $u2], $T, 2)
+        @benchmarkable effective_input_mode([$u1, $u2], $T, 2)
 
     return nothing
 end

benchmarks/runbenchmarks.jl

@@ -5,6 +5,7 @@ using QuantumOptics
 using QuantumOpticsBase
 using SymbolicUtils
 using LinearAlgebra
+using FunctionWrappers
 
 const SUITE = BenchmarkGroup()
 

benchmarks/slh_algebra.jl

@@ -8,7 +8,6 @@ function benchmark_slh_algebra!(SUITE)
 
     SUITE["SLH Algebra"]["symbolic"] = BenchmarkGroup()
 
-    # Setup: 3-cavity cascade (example 01-1)
     hu1 = FockSpace(:u1)
     hc1 = FockSpace(:c1)
     hv1 = FockSpace(:v1)
@@ -27,8 +26,8 @@ function benchmark_slh_algebra!(SUITE)
 
     SUITE["SLH Algebra"]["symbolic"]["3-cavity cascade"] = @benchmarkable begin
         G_cas = ▷($G_u, $G_c, $G_v)
-        get_hamiltonian(G_cas)
-        get_lindblad(G_cas)
+        hamiltonian(G_cas)
+        lindblad(G_cas)
     end
 
     SUITE["SLH Algebra"]["symbolic"]["concatenate + cascade"] = @benchmarkable begin
@@ -37,7 +36,7 @@ function benchmark_slh_algebra!(SUITE)
         G_R ⊞ G_L
     end
 
-    # Feedback: coherent-feedback OPO loop (from test_feedback.jl)
+    # Feedback: coherent-feedback OPO loop
     hs = FockSpace(:s)
     a_fb = Destroy(hs, :a, 1)
     κ_fb = 0.7
@@ -52,11 +51,10 @@ function benchmark_slh_algebra!(SUITE)
     SUITE["SLH Algebra"]["symbolic"]["feedback OPO loop"] =
         @benchmarkable feedback($G_fb_unconnected, 1 => 2, 2 => 1)
 
-    ## --- Numeric SLH (SLHqo with QuantumOptics operators) ---
+    ## --- Numeric SLH (unified SLH with QuantumOptics operators) ---
 
-    SUITE["SLH Algebra"]["numeric (SLHqo)"] = BenchmarkGroup()
+    SUITE["SLH Algebra"]["numeric"] = BenchmarkGroup()
 
-    # Setup: 2-QD bidirectional waveguide (example 05-2, realistic basis)
     N = 2
     bu = FockBasis(20)
     ba = NLevelBasis(2)
@@ -76,43 +74,45 @@ function benchmark_slh_algebra!(SUITE)
     t0 = 4σt
     T = [0:0.005:1;] * 3t0
     u1(t) = 1 / (sqrt(σt) * π^(1 / 4)) * exp(-(t - t0)^2 / (2 * σt^2))
-    gu_t = u_to_gu(u1, T)
-
-    G_u_qo = SLHqo(1, t -> gu_t(t) * a_u, 0 * one(b))
-    G_ϕ_12 = SLHqo(exp(1im * ϕn[1]), 0 * one(b), 0 * one(b))
-    G_R1 = SLHqo(1, √(γRn[1]) * σ(1, 1, 2), 0 * one(b))
-    G_R2 = SLHqo(1, √(γRn[2]) * σ(2, 1, 2), 0 * one(b))
-    G_L1 = SLHqo(1, √(γLn[1]) * σ(1, 1, 2), 0 * one(b))
-    G_L2 = SLHqo(1, √(γLn[2]) * σ(2, 1, 2), 0 * one(b))
-
-    SUITE["SLH Algebra"]["numeric (SLHqo)"]["2-QD waveguide composition"] =
-        @benchmarkable begin
-            G_R_t = $G_u_qo ▷ $G_R1 ▷ $G_ϕ_12 ▷ $G_R2
-            G_L_t = $G_L2 ▷ $G_ϕ_12 ▷ $G_L1
-            G_t = G_R_t ⊞ G_L_t
-            get_hamiltonian(G_t)
-            get_lindblad(G_t)
-        end
+    gu_t = coupling_input(u1, T)
+
+    G_u_qo = SLH(1, t -> gu_t(t) * a_u, 0 * one(b))
+    G_ϕ_12 = SLH(exp(1im * ϕn[1]), 0 * one(b), 0 * one(b))
+    G_R1 = SLH(1, √(γRn[1]) * σ(1, 1, 2), 0 * one(b))
+    G_R2 = SLH(1, √(γRn[2]) * σ(2, 1, 2), 0 * one(b))
+    G_L1 = SLH(1, √(γLn[1]) * σ(1, 1, 2), 0 * one(b))
+    G_L2 = SLH(1, √(γLn[2]) * σ(2, 1, 2), 0 * one(b))
+
+    SUITE["SLH Algebra"]["numeric"]["2-QD waveguide composition"] = @benchmarkable begin
+        G_R_t = $G_u_qo ▷ $G_R1 ▷ $G_ϕ_12 ▷ $G_R2
+        G_L_t = $G_L2 ▷ $G_ϕ_12 ▷ $G_L1
+        G_t = G_R_t ⊞ G_L_t
+        hamiltonian(G_t)
+        lindblad(G_t)
+    end
 
     ## --- Time-dependent closure evaluation (the ODE hot loop) ---
 
     SUITE["SLH Algebra"]["closure evaluation"] = BenchmarkGroup()
 
-    # Build the composed network once, then benchmark calling H(t) and L(t)
     G_R_t = G_u_qo ▷ G_R1 ▷ G_ϕ_12 ▷ G_R2
     G_L_t = G_L2 ▷ G_ϕ_12 ▷ G_L1
     G_t = G_R_t ⊞ G_L_t
 
-    H_f = get_hamiltonian(G_t)
-    L_f = get_lindblad(G_t)
+    H_f = hamiltonian(G_t)
+    L_f = lindblad(G_t)
 
     t_mid = T[length(T)÷2]
 
+    _callable(x) = x isa Union{Function,FunctionWrappers.FunctionWrapper}
+    Hf = _callable(H_f) ? H_f : (t -> H_f)
+    L_callables = [_callable(Li) ? Li : (t -> Li) for Li in L_f]
+
     SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide H(t)"] =
-        @benchmarkable $H_f($t_mid)
+        @benchmarkable $Hf($t_mid)
 
     SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide L(t)"] = @benchmarkable begin
-        for Li in $L_f
+        for Li in $L_callables
             Li($t_mid)
         end
     end

benchmarks/translation.jl

@@ -29,8 +29,8 @@ function benchmark_translation!(SUITE)
     G_c = SLH(1, √(γ_sym) * c_, Δ_sym * c_'c_)
     G_v = SLH(1, gv_sym * av_, 0)
     G_cas = ▷(G_u, G_c, G_v)
-    H_sym = get_hamiltonian(G_cas)
-    L_sym = get_lindblad(G_cas)[1]
+    H_sym = hamiltonian(G_cas)
+    L_sym = lindblad(G_cas)[1]
 
     ## --- Static translation (no time dependence) ---
 
@@ -42,7 +42,6 @@ function benchmark_translation!(SUITE)
 
     dict_p_static = Dict([κ_R, κ_L, Δ_sym] .=> [1.5, 1.0, 0.2])
 
-    # Composite expression: operator + parameter + transition
     expr_composite = a * 3 + Δ_sym * σ(2, 2)
 
     SUITE["Translation"]["static"]["atom-cavity"] =
@@ -55,7 +54,6 @@ function benchmark_translation!(SUITE)
     E_t(t) = 2 * t + 1im
     dict_p_t = Dict(E => E_t)
 
-    # Expression with time-dependent + static parameters and operators
     expr_td = a * 3 * conj(E) + Δ_sym * σ(2, 2)
 
     SUITE["Translation"]["time-dependent"]["atom-cavity"] = @benchmarkable translate_qo(
@@ -65,15 +63,14 @@ function benchmark_translation!(SUITE)
         time_parameter = $dict_p_t,
     )
 
-    # Full cascade H and L translation (realistic workflow from example 06-1)
+    # Full cascade H and L translation
     γ_ = 1.0
     σ_pulse = 1 / γ_
     T = [0:0.002:1;] * 12σ_pulse
     u_pulse(t) = 1 / (sqrt(σ_pulse) * π^(1 / 4)) * exp(-(t - 4σ_pulse)^2 / (2 * σ_pulse^2))
-    gu_t = u_to_gu(u_pulse, T)
-    gv_t = v_to_gv(u_pulse, T)
+    gu_t = coupling_input(u_pulse, T)
+    gv_t = coupling_output(u_pulse, T)
 
-    # Realistic basis sizes from interaction picture example (n_ph=20)
     bu = FockBasis(20)
     bc3 = FockBasis(6)
     bv = FockBasis(6)

docs/Project.toml

@@ -9,6 +9,7 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuantumCumulants = "35bcea6d-e19f-57db-af74-8011de6c7255"
+QuantumInputOutput = "18f9eda6-924c-47c7-a881-996695cfd7c6"
 QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
 QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
 SecondQuantizedAlgebra = "f7aa4685-e143-4cb6-a7f3-073579757907"