feat: SciMLOperators-backed lazy operator extension and benchmark suite

Project.toml

@@ -26,6 +26,7 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 
 [extensions]
 QuantumOpticsBaseMakieExt = ["Makie"]
+QuantumOpticsBaseSciMLOperatorsExt = "SciMLOperators"
 
 [compat]
 Adapt = "4.1"

benchmark/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"

benchmark/sciml_lazy_benchmark.jl

@@ -0,0 +1,122 @@
+# benchmark/sciml_lazy_benchmark.jl
+#
+# Comparative benchmark: QuantumOpticsBase native lazy operators
+# vs SciMLOperators-backed prototype.
+#
+# Usage:
+#   cd QuantumOpticsBase.jl
+#   julia --project=benchmark benchmark/sciml_lazy_benchmark.jl
+#
+# Output: a Markdown table of minimum times (BenchmarkTools, 200-sample budget).
+
+using BenchmarkTools
+using QuantumOpticsBase
+using SciMLOperators
+using LinearAlgebra
+using Printf
+
+BenchmarkTools.DEFAULT_PARAMETERS.samples  = 200
+BenchmarkTools.DEFAULT_PARAMETERS.seconds  = 2.0
+BenchmarkTools.DEFAULT_PARAMETERS.evals    = 1
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+fmt(t) = @sprintf("%.2f μs", t * 1e6)
+
+function run_triple(label, lazy_op, psi)
+    w       = sciml_lazy_operator(lazy_op)
+    w_cache = cache_sciml_lazy_operator(w, psi.data)
+
+    t_lazy   = @belapsed $lazy_op * $psi
+    t_sciml  = @belapsed $w      * $psi
+    t_cached = @belapsed $w_cache * $psi
+
+    @printf("| %-45s | %10s | %14s | %13s |\n",
+        label, fmt(t_lazy), fmt(t_sciml), fmt(t_cached))
+end
+
+# ---------------------------------------------------------------------------
+# Print header
+# ---------------------------------------------------------------------------
+println()
+println("## SciMLOperators prototype benchmark")
+println()
+println("Julia $(VERSION), SciMLOperators $(pkgversion(SciMLOperators))")
+println()
+@printf("| %-45s | %10s | %14s | %13s |\n",
+    "case", "native lazy", "SciML uncached", "SciML cached")
+println("|" * "-"^47 * "|" * "-"^12 * "|" * "-"^16 * "|" * "-"^15 * "|")
+
+b0 = SpinBasis(1//2)
+b1 = SpinBasis(1)
+sx, sy, sz = sigmax(b0), sigmay(b0), sigmaz(b0)
+sp, sm     = sigmap(b0), sigmam(b0)
+
+# ---------------------------------------------------------------------------
+# 1. LazySum: local transverse-field terms on spin chains (n = 4, 6, 8)
+# ---------------------------------------------------------------------------
+for n in (4, 6, 8)
+    b  = tensor(fill(b0, n)...)
+    H  = LazySum([LazyTensor(b, k, sx) for k in 1:n]...)
+    psi = Ket(b, normalize!(randn(ComplexF64, length(b))))
+    run_triple("LazySum  n=$n (transverse field)", H, psi)
+end
+
+# ---------------------------------------------------------------------------
+# 2. LazyProduct: depth-2, -4, -6 chains
+# ---------------------------------------------------------------------------
+for depth in (2, 4, 6)
+    n  = max(depth, 4)
+    b  = tensor(fill(b0, n)...)
+    psi = Ket(b, normalize!(randn(ComplexF64, length(b))))
+    ops = [LazyTensor(b, mod1(k, n), sx) for k in 1:depth]
+    P  = LazyProduct(ops...)
+    run_triple("LazyProduct  depth=$depth n=$n", P, psi)
+end
+
+# ---------------------------------------------------------------------------
+# 3. LazyTensor: edge vs mid, dense vs sparse, spin-½ n=6
+# ---------------------------------------------------------------------------
+let n = 6, b = tensor(fill(b0, n)...)
+    psi = Ket(b, normalize!(randn(ComplexF64, length(b))))
+
+    run_triple("LazyTensor  edge-1  dense  n=$n", LazyTensor(b, 1, sx), psi)
+    run_triple("LazyTensor  edge-$n dense  n=$n", LazyTensor(b, n, sz), psi)
+    run_triple("LazyTensor  mid-3   dense  n=$n", LazyTensor(b, 3, sy), psi)
+    run_triple("LazyTensor  mid-3   sparse n=$n", LazyTensor(b, 3, SparseOperator(sy)), psi)
+    run_triple("LazyTensor  edge-1  sparse n=$n", LazyTensor(b, 1, SparseOperator(sx)), psi)
+end
+
+# ---------------------------------------------------------------------------
+# 4. LazyTensor: spin-1 local dimension (d=3)
+# ---------------------------------------------------------------------------
+let n = 4, b = tensor(fill(b1, n)...),
+    sx1 = sigmax(b1), sz1 = sigmaz(b1)
+    psi = Ket(b, normalize!(randn(ComplexF64, length(b))))
+
+    run_triple("LazyTensor  mid-2   spin-1 n=$n", LazyTensor(b, 2, sx1), psi)
+    run_triple("LazyTensor  mid-3   spin-1 n=$n", LazyTensor(b, 3, sz1), psi)
+end
+
+# ---------------------------------------------------------------------------
+# 5. Mixed Heisenberg-like Hamiltonian (nearest-neighbour XX+YY+ZZ)
+# ---------------------------------------------------------------------------
+for n in (4, 6)
+    b  = tensor(fill(b0, n)...)
+    psi = Ket(b, normalize!(randn(ComplexF64, length(b))))
+    terms = AbstractOperator[]
+    for k in 1:(n-1)
+        push!(terms, LazyProduct(LazyTensor(b, k, sx), LazyTensor(b, k+1, sx)))
+        push!(terms, LazyProduct(LazyTensor(b, k, sy), LazyTensor(b, k+1, sy)))
+        push!(terms, LazyProduct(LazyTensor(b, k, sz), LazyTensor(b, k+1, sz)))
+    end
+    H = LazySum(terms...)
+    run_triple("Mixed Heisenberg  n=$n ($(length(terms)) terms)", H, psi)
+end
+
+println()
+println("> Timings are minimum over $(BenchmarkTools.DEFAULT_PARAMETERS.samples) samples.")
+println("> 'SciML uncached' allocates intermediate buffers on every call.")
+println("> 'SciML cached'  pre-allocates via `cache_sciml_lazy_operator`.")

ext/QuantumOpticsBaseSciMLOperatorsExt.jl

@@ -0,0 +1,181 @@
+module QuantumOpticsBaseSciMLOperatorsExt
+
+using QuantumOpticsBase
+using SciMLOperators
+import LinearAlgebra
+import LinearAlgebra: mul!, I
+import QuantumOpticsBase: AbstractOperator, Operator, Ket, Bra,
+    LazySum, LazyProduct, LazyTensor, CompositeBasis,
+    dense, dagger, check_samebases
+
+export SciMLOperatorWrapper
+
+# ---------------------------------------------------------------------------
+# Core wrapper
+# ---------------------------------------------------------------------------
+
+"""
+    SciMLOperatorWrapper{BL,BR,T} <: AbstractOperator{BL,BR}
+
+Wraps an `AbstractSciMLOperator` while preserving QuantumOpticsBase basis
+information. Construct via [`sciml_lazy_operator`](@ref).
+"""
+struct SciMLOperatorWrapper{BL,BR,T} <: AbstractOperator{BL,BR}
+    basis_l::BL
+    basis_r::BR
+    sciml_op::T
+end
+
+Base.eltype(::SciMLOperatorWrapper{BL,BR,T}) where {BL,BR,T} = eltype(T)
+Base.size(w::SciMLOperatorWrapper) = (length(w.basis_l), length(w.basis_r))
+
+# ---------------------------------------------------------------------------
+# Conversion helpers: QO types → SciML operator trees
+# ---------------------------------------------------------------------------
+
+function _qo_to_sciml(op::Operator)
+    MatrixOperator(op.data)
+end
+
+function _qo_to_sciml(op::AbstractOperator)
+    MatrixOperator(dense(op).data)
+end
+
+function _qo_to_sciml(op::LazySum)
+    if isempty(op.operators)
+        T = eltype(op.factors)
+        return MatrixOperator(zeros(T, length(op.basis_l), length(op.basis_r)))
+    end
+    result = op.factors[1] * _qo_to_sciml(op.operators[1])
+    for i in 2:length(op.operators)
+        result = result + op.factors[i] * _qo_to_sciml(op.operators[i])
+    end
+    result
+end
+
+function _qo_to_sciml(op::LazyProduct)
+    isempty(op.operators) &&
+        error("Cannot convert an empty LazyProduct to SciMLOperators.")
+    result = _qo_to_sciml(op.operators[1])
+    for i in 2:length(op.operators)
+        result = result * _qo_to_sciml(op.operators[i])
+    end
+    isone(op.factor) ? result : op.factor * result
+end
+
+function _qo_to_sciml(op::LazyTensor)
+    cb     = op.basis_l
+    nsites = length(cb.bases)
+    T      = eltype(op)
+    site_ops = map(1:nsites) do k
+        idx = findfirst(==(k), op.indices)
+        if idx !== nothing
+            _qo_to_sciml(op.operators[idx])
+        else
+            d = length(cb.bases[k])
+            MatrixOperator(Matrix{T}(I, d, d))
+        end
+    end
+    # QuantumOpticsBase uses Fortran/column-major ordering: site 1 is the
+    # fastest-varying index (stride 1). SciMLOperators' TensorProductOperator
+    # uses C/row-major ordering: the first argument is slowest-varying.
+    # Reversing aligns the two conventions.
+    tp = TensorProductOperator(reverse(site_ops)...)
+    isone(op.factor) ? tp : op.factor * tp
+end
+
+# ---------------------------------------------------------------------------
+# Public API
+# ---------------------------------------------------------------------------
+
+"""
+    sciml_lazy_operator(op::AbstractOperator) -> SciMLOperatorWrapper
+
+Convert a QuantumOpticsBase lazy operator to a SciMLOperators-backed wrapper,
+preserving `basis_l` and `basis_r`. Requires `SciMLOperators` to be loaded.
+
+```julia
+using QuantumOpticsBase, SciMLOperators
+b0 = SpinBasis(1//2)
+b  = tensor(b0, b0, b0, b0)
+H  = LazySum(LazyTensor(b,1,sigmax(b0)), LazyTensor(b,2,sigmaz(b0)))
+H_sciml = sciml_lazy_operator(H)
+psi = Ket(b, randn(ComplexF64, length(b)))
+@assert dense(H_sciml) ≈ dense(H)
+@assert H_sciml * psi  ≈ H * psi
+```
+"""
+function QuantumOpticsBase.sciml_lazy_operator(op::AbstractOperator)
+    SciMLOperatorWrapper(op.basis_l, op.basis_r, _qo_to_sciml(op))
+end
+
+"""
+    cache_sciml_lazy_operator(w::SciMLOperatorWrapper, u::AbstractVector)
+
+Pre-allocate intermediate buffers via `SciMLOperators.cache_operator`.
+Repeated `mul!` calls on the returned wrapper avoid per-call heap allocation.
+"""
+function QuantumOpticsBase.cache_sciml_lazy_operator(w::SciMLOperatorWrapper, u::AbstractVector)
+    SciMLOperatorWrapper(w.basis_l, w.basis_r, cache_operator(w.sciml_op, u))
+end
+
+# ---------------------------------------------------------------------------
+# QuantumOpticsBase interface
+# ---------------------------------------------------------------------------
+
+function QuantumOpticsBase.dense(w::SciMLOperatorWrapper)
+    n   = length(w.basis_l)
+    m   = length(w.basis_r)
+    dat = zeros(ComplexF64, n, m)
+    e_j = zeros(ComplexF64, m)
+    col = zeros(ComplexF64, n)
+    # TensorProductOperator (and any tree containing it) requires cache before mul!.
+    # dense() is O(n^2) allocation anyway, so caching here is fine.
+    cached_op = cache_operator(w.sciml_op, e_j)
+    for j in 1:m
+        e_j[j] = one(ComplexF64)
+        mul!(col, cached_op, e_j)
+        dat[:, j] .= col
+        e_j[j] = zero(ComplexF64)
+    end
+    Operator(w.basis_l, w.basis_r, dat)
+end
+
+function mul!(result::Ket{B1}, op::SciMLOperatorWrapper{B1,B2},
+              psi::Ket{B2}, alpha::Number, beta::Number) where {B1,B2}
+    cached_op = cache_operator(op.sciml_op, psi.data)
+    mul!(result.data, cached_op, psi.data, alpha, beta)
+    return result
+end
+
+function Base.:*(op::SciMLOperatorWrapper{B1,B2}, psi::Ket{B2}) where {B1,B2}
+    out = zeros(eltype(psi), length(op.basis_l))
+    cached_op = cache_operator(op.sciml_op, psi.data)
+    mul!(out, cached_op, psi.data)
+    return Ket(op.basis_l, out)
+end
+
+Base.:*(α::Number, op::SciMLOperatorWrapper) =
+    SciMLOperatorWrapper(op.basis_l, op.basis_r, α * op.sciml_op)
+Base.:*(op::SciMLOperatorWrapper, α::Number) =
+    SciMLOperatorWrapper(op.basis_l, op.basis_r, op.sciml_op * α)
+Base.:/(op::SciMLOperatorWrapper, α::Number) =
+    SciMLOperatorWrapper(op.basis_l, op.basis_r, op.sciml_op / α)
+
+function Base.:+(a::SciMLOperatorWrapper{B1,B2},
+                  b::SciMLOperatorWrapper{B1,B2}) where {B1,B2}
+    check_samebases(a, b)
+    SciMLOperatorWrapper(a.basis_l, a.basis_r, a.sciml_op + b.sciml_op)
+end
+
+function Base.:*(a::SciMLOperatorWrapper{B1,B2},
+                  b::SciMLOperatorWrapper{B2,B3}) where {B1,B2,B3}
+    SciMLOperatorWrapper(a.basis_l, b.basis_r, a.sciml_op * b.sciml_op)
+end
+
+function QuantumOpticsBase.dagger(w::SciMLOperatorWrapper)
+    d = dense(w)
+    SciMLOperatorWrapper(w.basis_r, w.basis_l, MatrixOperator(adjoint(d.data)))
+end
+
+end  # module QuantumOpticsBaseSciMLOperatorsExt

src/QuantumOpticsBase.jl

@@ -107,4 +107,11 @@ include("printing.jl")
 include("apply.jl")
 include("visualization.jl")
 
+# SciMLOperators extension — stub declarations so these symbols live in the
+# QuantumOpticsBase namespace and are accessible via `using QuantumOpticsBase`.
+# Methods are added by QuantumOpticsBaseSciMLOperatorsExt when SciMLOperators loads.
+function sciml_lazy_operator end
+function cache_sciml_lazy_operator end
+export sciml_lazy_operator, cache_sciml_lazy_operator
+
 end # module

test/Project.toml

@@ -11,8 +11,10 @@ OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 QuantumInterface = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
 QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
+QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"