Add benchmark suite

.github/workflows/Benchmarks.yaml

@@ -0,0 +1,67 @@
+name: Benchmark Tracking
+on:
+  push:
+    branches:
+      - 'main'
+    paths:
+      - '.github/workflows/Benchmarks.yml'
+      - 'src/**'
+      - 'ext/**'
+      - 'test/runtests.jl'
+      - 'Project.toml'
+  pull_request:
+    branches:
+      - 'main'
+    paths:
+      - '.github/workflows/Benchmarks.yml'
+      - 'src/**'
+      - 'ext/**'
+      - 'benchmarks/**'
+      - 'Project.toml'
+    types:
+      - opened
+      - reopened
+      - synchronize
+      - ready_for_review
+  workflow_dispatch:
+
+permissions:
+    actions: write
+    contents: write
+    deployments: write
+
+jobs:
+  benchmark:
+    runs-on: ubuntu-latest
+    if: ${{ !github.event.pull_request.draft }}
+    steps:
+      - uses: actions/checkout@v6
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: '1.12' # should run on same julia version 
+          arch: x64
+      - uses: julia-actions/cache@v3
+      - name: Run benchmark
+        run: |
+          cd benchmarks
+          julia --project --threads=2 --color=yes -e '
+            using Pkg;
+            Pkg.develop(PackageSpec(path=joinpath(pwd(), "..")));
+            Pkg.instantiate();
+            include("runbenchmarks.jl")'
+
+      # this will update benchmarks/data.js in gh-pages branch
+      - name: Parse & Upload Benchmark Results
+        uses: benchmark-action/github-action-benchmark@v1
+        with:
+          name: Benchmark Results
+          tool: "julia"
+          output-file-path: benchmarks/benchmarks_output.json
+          github-token: ${{ secrets.GITHUB_TOKEN }}
+          alert-threshold: "130%"
+          fail-threshold: "170%"
+          comment-on-alert: true
+          fail-on-alert: true
+          benchmark-data-dir-path: benchmarks
+          max-items-in-chart: 100
+          auto-push: ${{ github.event_name != 'pull_request' }}

.gitignore

@@ -30,4 +30,5 @@ JuliaLocalPreferences.toml
 # additional
 /.vscode/
 docs/src/examples/jupyter_notebooks/.ipynb_checkpoints
-docs/Manifest.toml
+docs/Manifest.toml
+benchmarks_output.json

Makefile

@@ -20,7 +20,7 @@ docs:
 
 bench:
 	${JULIA} --project=benchmarks -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
-	${JULIA} --project=docs benchmarks/runbenchmarks.jl
+	${JULIA} --project=benchmarks benchmarks/runbenchmarks.jl
 
 
 all: setup format test docs 

benchmarks/Project.toml

@@ -0,0 +1,8 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+QuantumInputOutput = "18f9eda6-924c-47c7-a881-996695cfd7c6"
+QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
+QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
+SecondQuantizedAlgebra = "f7aa4685-e143-4cb6-a7f3-073579757907"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"

benchmarks/correlations.jl

@@ -0,0 +1,67 @@
+"""
+Correlation function benchmarks — two-time correlation matrices.
+Based on example 01-1: single photon cavity scattering.
+"""
+function benchmark_correlations!(SUITE)
+    SUITE["Correlations"] = BenchmarkGroup()
+
+    ## Setup: single photon scattering through a cavity
+    hu1 = FockSpace(:u1)
+    hc1 = FockSpace(:c1)
+    hv1 = FockSpace(:v1)
+    h = hu1 ⊗ hc1 ⊗ hv1
+
+    au = Destroy(h, :a_u, 1)
+    c = Destroy(h, :c, 2)
+    av = Destroy(h, :a_v, 3)
+
+    gu, Δ, γ = rnumbers("g_u Δ γ")
+    gv = cnumber("g_v")
+
+    G_u = SLH(1, gu * au, 0)
+    G_c = SLH(1, √(γ) * c, Δ * c'c)
+    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]
+
+    γ_ = 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)
+
+    bu1 = FockBasis(1)
+    bc1 = FockBasis(1)
+    bv1 = FockBasis(1)
+    b = bu1 ⊗ bc1 ⊗ bv1
+
+    dict_p = Dict([γ, Δ, gv] .=> [γ_, 0.0, 0])
+    dict_p_t = Dict(gu => 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)
+
+    function input_output(t, ρ)
+        Ht = H_QO(t)
+        J = [L_QO(t)]
+        return Ht, J, dagger.(J)
+    end
+
+    ψ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)
+    Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo
+
+    ## --- Two-time correlation ---
+
+    SUITE["Correlations"]["two-time"] = BenchmarkGroup()
+
+    SUITE["Correlations"]["two-time"]["single photon cavity"] =
+        @benchmarkable two_time_corr_matrix($T, $ρt, $input_output, $Ls)
+
+    return nothing
+end

benchmarks/interaction_picture.jl

@@ -0,0 +1,88 @@
+"""
+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()
+
+    ## Setup: u -> TLS -> v system
+    γ_ = 1.0
+    τ = 1 / γ_
+    t_p = 4 / γ_
+    u(t) = 1 / (sqrt(τ) * π^(1 / 4)) * exp(-0.5 * ((t - t_p) / τ)^2)
+
+    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)
+
+    ## --- Coupling matrix A(t) evaluation (called every ODE step) ---
+
+    SUITE["Interaction Picture"]["coupling matrix evaluation"] = BenchmarkGroup()
+
+    t_mid = T[length(T)÷2]
+
+    # 2-mode A(t)
+    SUITE["Interaction Picture"]["coupling matrix evaluation"]["2 modes"] =
+        @benchmarkable $A_uv($t_mid)
+
+    # 4-mode A(t) — more realistic for multi-port systems
+    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)
+
+    SUITE["Interaction Picture"]["coupling matrix evaluation"]["4 modes"] =
+        @benchmarkable $A_4m($t_mid)
+
+    ## --- Coefficient matrix M(t) ---
+
+    SUITE["Interaction Picture"]["coefficient matrix M"] = BenchmarkGroup()
+
+    SUITE["Interaction Picture"]["coefficient matrix M"]["numerical (ODE)"] =
+        @benchmarkable interaction_picture_M($A_uv, $T)
+
+    SUITE["Interaction Picture"]["coefficient matrix M"]["analytical (2 equal modes)"] =
+        @benchmarkable interaction_picture_M_2modes_equal($u, $T)
+
+    ## --- Symbolic operator substitution ---
+
+    hu = FockSpace(:u)
+    hs = NLevelSpace(:s, 2)
+    hv = FockSpace(:v)
+    h = hu ⊗ hs ⊗ hv
+
+    au_sym = Destroy(h, :a_u, 1)
+    av_sym = Destroy(h, :a_v, 3)
+    σ_sym = Transition(h, :σ, 1, 2, 2)
+
+    gu_sym, γ_sym, gv_sym = rnumbers("gu γ gv")
+
+    G_u = SLH(1, gu_sym' * au_sym, 0)
+    G_s = SLH(1, sqrt(γ_sym) * σ_sym, 0)
+    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_int_ = simplify(H - H_uv)
+
+    M_sym(i, j) = cnumber("M_{$(i)$(j)}")
+    a0_ls = [au_sym, av_sym]
+    la = length(a0_ls)
+    a_int_ls = [sum(M_sym(i, j) * a0_ls[j] for j = 1:la) for i = 1:la]
+    int_dict = Dict(a0_ls .=> a_int_ls)
+
+    SUITE["Interaction Picture"]["operator substitution"] = BenchmarkGroup()
+
+    SUITE["Interaction Picture"]["operator substitution"]["TLS cascade"] =
+        @benchmarkable begin
+            simplify(substitute_operators($H_int_, $int_dict))
+            simplify(substitute_operators($L, $int_dict))
+        end
+
+    return nothing
+end

benchmarks/pulse_couplings.jl

@@ -0,0 +1,40 @@
+"""
+Pulse coupling benchmarks — converting temporal mode shapes to virtual cavity couplings.
+"""
+function benchmark_pulse_couplings!(SUITE)
+    SUITE["Pulse Couplings"] = BenchmarkGroup()
+
+    ## Setup: Gaussian pulse on a realistic time grid
+    τ = 3.0
+    σ = 0.5
+    δ = 0.3
+    T = [0.0:0.0001:1;] * 6  # 60001 points
+
+    u(t) = 1 / (sqrt(σ) * π^(1 / 4)) * exp(-0.5 * ((t - τ) / σ)^2) * exp(1im * δ * t)
+    v(t) = u(t)
+
+    ## --- Single-pulse coupling computation ---
+
+    SUITE["Pulse Couplings"]["single-pulse"] = BenchmarkGroup()
+
+    SUITE["Pulse Couplings"]["single-pulse"]["input"] = @benchmarkable u_to_gu($u, $T)
+
+    SUITE["Pulse Couplings"]["single-pulse"]["output"] = @benchmarkable v_to_gv($v, $T)
+
+    ## --- Effective multi-pulse couplings ---
+
+    SUITE["Pulse Couplings"]["multi-pulse"] = BenchmarkGroup()
+
+    v1(t) = 1 / (sqrt(σ) * π^(1 / 4)) * exp(-0.5 * ((t - (τ - 2σ)) / σ)^2)
+    v2(t) = 1 / (sqrt(σ) * π^(1 / 4)) * exp(-0.5 * ((t - (τ + 2σ)) / σ)^2)
+    u1(t) = 1 / (sqrt(σ) * π^(1 / 4)) * exp(-0.5 * ((t - (τ - 2σ)) / σ)^2)
+    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)
+
+    SUITE["Pulse Couplings"]["multi-pulse"]["input 2 modes"] =
+        @benchmarkable u_eff([$u1, $u2], $T, 2)
+
+    return nothing
+end

benchmarks/runbenchmarks.jl

@@ -0,0 +1,27 @@
+using BenchmarkTools
+using QuantumInputOutput
+using SecondQuantizedAlgebra
+using QuantumOptics
+using QuantumOpticsBase
+using SymbolicUtils
+using LinearAlgebra
+
+const SUITE = BenchmarkGroup()
+
+include("slh_algebra.jl")
+include("translation.jl")
+include("pulse_couplings.jl")
+include("interaction_picture.jl")
+include("correlations.jl")
+
+benchmark_slh_algebra!(SUITE)
+benchmark_translation!(SUITE)
+benchmark_pulse_couplings!(SUITE)
+benchmark_interaction_picture!(SUITE)
+benchmark_correlations!(SUITE)
+
+BenchmarkTools.tune!(SUITE)
+results = BenchmarkTools.run(SUITE; verbose = true)
+display(median(results))
+
+BenchmarkTools.save("benchmarks_output.json", median(results))