feat: backend agnostic numerical conversion

CLAUDE.md

@@ -15,7 +15,15 @@ src/SecondQuantizedAlgebra.jl       # Main module, exports, imports
 src/types.jl                        # QField, QSym abstract hierarchy
 src/precompile.jl                   # PrecompileTools workload
 src/average.jl                      # AvgFunc/AvgSym, average(), undo_average()
-src/numeric.jl                      # to_numeric(), numeric_average() via QuantumOpticsBase
+
+src/numeric/backend.jl              # NumericBackend + singletons; open hooks; _default_backend
+src/numeric/coeff.jl                # _to_complex family, constant folder, parameter/time splitting
+src/numeric/core.jl                 # NumericContext, _numeric_leaf, _to_numeric_static/_td
+src/numeric/indexed.jl              # backend-neutral indexed unroll (sites, ne, sub_op/sub_coef)
+src/numeric/api.jl                  # public to_numeric / numeric_average / expect (QI types)
+
+ext/SecondQuantizedAlgebraQuantumOpticsBaseExt.jl  # QuantumOpticsBase backend (vector LazySum)
+ext/SecondQuantizedAlgebraQuantumToolboxExt.jl     # QuantumToolbox backend (VecSum over QobjEvo)
 
 src/expressions/cnum.jl             # CNum = Complex{Num} arithmetic, fast paths, constants
 src/expressions/qterm.jl            # QTerm struct (ops, ne) — dict key for QAdd
@@ -72,7 +80,7 @@ HilbertSpace (abstract)
 - **Dict-based term storage**: `QAdd` stores `Dict{QTerm, CNum}` where `QTerm` bundles `ops::Vector{QSym}` with `ne::Vector{NonEqualPair}` index-inequality scope, plus a cached `hash::UInt` (computed once at construction; the key is hashed repeatedly per dict insert/probe/rehash). Like terms are collected on construction.
 - **CNum prefactors**: prefactors are `Complex{Num}` (from Symbolics.jl), not parameterized. Dedicated fast paths in `cnum.jl` short-circuit for numeric (non-symbolic) cases.
 - **Site-indexed operators**: each `Op` carries `space_index` and `index::Index`. Operators interact only if `_same_site(a, b)`.
-- **Five operator hooks**: `_site_compare`, `_can_commute`, `_commute_pair`, `_reduce_pair`, `_ground_state_expand` are each a single `(::Op, ::Op)` method (in `operators.jl`) that branches on `kind`. The algebra talks to operators exclusively through these. A sixth, defaulted hook `_may_reduce(a, b)::Bool` gates the reduce pass (`true` only for `Transition`/`Pauli` pairs). Adding a new operator role means adding an `OP_*` enum arm, a constructor, an `is_*` predicate, and a `kind` branch in each hook plus `adjoint`/`order_key`/`to_numeric`/printing. With the concrete `Op` eltype the hooks now infer concrete return types (`_commute_pair`/`_reduce_pair` return `Tuple{Op, …}`), so `_may_reduce`'s original boxing-avoidance role is moot; it remains as a cheap same-site skip.
+- **Five operator hooks**: `_site_compare`, `_can_commute`, `_commute_pair`, `_reduce_pair`, `_ground_state_expand` are each a single `(::Op, ::Op)` method (in `operators.jl`) that branches on `kind`. The algebra talks to operators exclusively through these. A sixth, defaulted hook `_may_reduce(a, b)::Bool` gates the reduce pass (`true` only for `Transition`/`Pauli` pairs). Adding a new operator role means adding an `OP_*` enum arm, a constructor, an `is_*` predicate, and a `kind` branch in each hook plus `adjoint`/`order_key`/`numeric_operator` (per backend extension)/printing. With the concrete `Op` eltype the hooks now infer concrete return types (`_commute_pair`/`_reduce_pair` return `Tuple{Op, …}`), so `_may_reduce`'s original boxing-avoidance role is moot; it remains as a cheap same-site skip.
 - **Concrete struct fields**: all struct fields are concretely typed (enforced by CheckConcreteStructs in tests).
 - **Index tracking**: `QAdd` carries `indices::Vector{Index}` for summation scope; per-term inequality constraints live on `QTerm.ne`, not on `QAdd`.
 
@@ -160,11 +168,13 @@ Before merging any PR:
 | Package | Purpose |
 |---------|---------|
 | Combinatorics | Product enumeration for operator generation |
-| QuantumOpticsBase | Numeric basis types (FockBasis, NLevelBasis, SpinBasis), `⊗`, LazyTensor |
+| QuantumInterface | Lightweight owner of `⊗`/`tensor`/`expect`/`basis` and `Basis`/`AbstractOperator`/`StateVector` types (hard dep) |
 | SymbolicUtils | Symbolic tree traversal interface |
 | Symbolics | Symbolic variables (`@variables`), `Num` type for CNum prefactors |
 | TermInterface | `iscall`, `operation`, `arguments` protocol |
 | Latexify | LaTeX rendering recipes |
 | PrecompileTools | `@setup_workload`/`@compile_workload` in `precompile.jl` |
 
-**Test dependencies:** Aqua, JET, CheckConcreteStructs, ExplicitImports, LaTeXStrings
+**Weak dependencies (numeric extensions):** `QuantumOpticsBase` (FockBasis/NLevelBasis/SpinBasis, vector `LazySum`/`TimeDependentSum`), `QuantumToolbox` + `SciMLOperators` (QuantumObject builders, the `VecSum` over `QobjEvo`). The numeric API errors until one backend is loaded.
+
+**Test dependencies:** Aqua, JET, CheckConcreteStructs, ExplicitImports, LaTeXStrings, QuantumOpticsBase, QuantumToolbox, SciMLOperators, LinearAlgebra

Changelog.md

@@ -4,6 +4,33 @@ All notable changes to [`SecondQuantizedAlgebra.jl`](https://github.com/qojulia/
 
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+## [v0.10.0]
+
+Numeric conversion (`to_numeric`/`numeric_average`/`expect`) was redesigned to be extensible, type-stable, lazy, and multi-backend. This is a breaking release.
+
+### Added
+
+- A second numeric backend, [QuantumToolbox.jl](https://github.com/qutip/QuantumToolbox.jl), alongside QuantumOpticsBase.jl. Both are wired in through Julia package extensions and selected with the backend singletons `QuantumOpticsBackend()` / `QuantumToolboxBackend()`.
+- A uniform Hilbert-space entry point `to_numeric(op, h::HilbertSpace, dims; backend, parameter, time_parameter, operators, adjoint_ops)`. It builds the backend basis from `dims` (Fock cutoff / spin number) and is the only form that works for both backends. The backend defaults to the single loaded one.
+- Open extension hooks `numeric_operator`, `numeric_basis`, `numeric_subbasis`, `numeric_embed`, `numeric_identity`, `numeric_assemble`, `numeric_assemble_td`, `numeric_expect` are exported. Downstream packages add methods to support custom operator roles or backends without editing this package.
+
+### Changed (breaking)
+
+- `to_numeric` of a `QAdd` (any sum or product) now returns a **lazy** operator: a vector-backed `LazySum` on QuantumOptics, a `QobjEvo` over a vector-backed `VecSum` on QuantumToolbox. The type is one concrete type for any number of terms, so `to_numeric`/`numeric_average` stay `@inferred`-stable. A single operator still returns the bare eager matrix. Materialize a lazy result with `dense`/`sparse` (QuantumOptics) or `QuantumToolbox.sparse`/`SciMLOperators.concretize` (QuantumToolbox). Tests that compared `to_numeric(a'*a, b)` against an eager matrix must materialize first.
+- The time-dependent form (`time_parameter` non-empty) returns the backend's **native** time-dependent operator: a `TimeDependentSum` (QuantumOptics) or a `QobjEvo` (QuantumToolbox), both callable at a time (`H(t)` / `H(p, t)`) and directly consumable by `mesolve`/`master_dynamic`/`sesolve`. It is no longer a `t -> op(t)` closure.
+- `QuantumOpticsBase` moved from a hard dependency to a weak dependency. Using the numeric API now requires loading a backend: add `using QuantumOpticsBase` (or `using QuantumToolbox`) next to `using SecondQuantizedAlgebra`. The lightweight `QuantumInterface.jl` is a new hard dependency (it owns the `⊗`/`tensor`/`expect`/`basis` generics that the algebra extends).
+
+### Removed (breaking)
+
+- The `op_type` keyword of `to_numeric` is gone. The result is lazy by default; call `dense`/`sparse` (or `concretize`) on it instead of passing `op_type=dense`/`sparse`.
+- The per-basis leaf methods `to_numeric(op::Op, b::FockBasis)` (and the `NLevelBasis`/`SpinBasis` variants) are gone. Custom backends implement `numeric_operator` instead.
+- The `FunctionWrappers.jl` dependency was dropped (the time-dependent path no longer wraps closures).
+
+### Notes
+
+- `expect`/`numeric_average` and `sesolve`/`schroedinger` are unaffected by laziness: they apply `H` to the state directly. Only `mesolve` and other superoperator solvers pay a cost, because QuantumToolbox builds the Liouvillian by tensoring `H` with identity, which a lazy `H` makes lazy (the "lazy tensor can hurt performance" warning). For fast `mesolve` on small or medium systems, concretize `H` first (`QuantumToolbox.sparse(H)` / `dense(H)`); keep it lazy only when materializing the Liouvillian is itself the problem.
+- A user doing `using SecondQuantizedAlgebra, QuantumToolbox` has two `⊗`/`tensor`/`expect` in scope (QuantumInterface's and QuantumToolbox's) and must qualify them; importing QuantumToolbox as `import QuantumToolbox as QTB` avoids the clash.
+
 ## [v0.9.0]
 
 ### Added

Makefile

@@ -7,7 +7,7 @@ setup:
 	${JULIA} -e 'using Pkg; Pkg.Apps.add("Runic")'
 
 format: ## Format all Julia files with Runic
-	runic --inplace src/ test/ benchmark/ examples/ docs/
+	runic --inplace src/ ext/ test/ benchmark/ examples/ docs/
 
 servedocs:
 	${JULIA} --project=docs -e 'using LiveServer; LiveServer.servedocs(skip_files=[joinpath("docs", "src", "changelog.md")])'

Project.toml

@@ -1,26 +1,36 @@
 name = "SecondQuantizedAlgebra"
 uuid = "f7aa4685-e143-4cb6-a7f3-073579757907"
-version = "0.9.0"
+version = "0.10.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
-FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
+QuantumInterface = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
 SciMLPublic = "431bcebd-1456-4ced-9d72-93c2757fff0b"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TermInterface = "8ea1fca8-c5ef-4a55-8b96-4e9afe9c9a3c"
 
+[weakdeps]
+QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
+QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
+
+[extensions]
+SecondQuantizedAlgebraQuantumOpticsBaseExt = "QuantumOpticsBase"
+SecondQuantizedAlgebraQuantumToolboxExt = ["QuantumToolbox", "SciMLOperators"]
+
 [compat]
 Combinatorics = "1.0.2"
-FunctionWrappers = "1"
 Latexify = "0.13, 0.14, 0.15, 0.16"
 MutableArithmetics = "1"
 PrecompileTools = "1"
-QuantumOpticsBase = "0.4, 0.5"
+QuantumInterface = "0.4"
+QuantumOpticsBase = "0.5"
+QuantumToolbox = "0.47"
+SciMLOperators = "1.22"
 SciMLPublic = "1"
 SymbolicUtils = "4.30"
 Symbolics = "7.26"

docs/src/devdocs.md

@@ -107,7 +107,7 @@ A sixth, defaulted hook supports the reduce pass: `_may_reduce(a, b)::Bool` answ
 
 **Site families inside `_site_compare`.** Cross-role comparison is family-scoped: Fock `{Destroy, Create}` and PhaseSpace `{Position, Momentum}` compare within their family (PhaseSpace ignores the operator name, treating x and p as conjugate variables on one site); the other roles are singleton families. Pairs in different families fall back to the `kind` integer order (`OP_DESTROY=0 < … < OP_MOMENTUM=6`), which is exactly the old `_type_order`.
 
-**Adding an operator role.** Add an `OP_*` enum arm, a constructor function and an `is_*` predicate, then a branch in each of the six hooks plus `adjoint`, `order_key`, `to_numeric`, and the printing/Latexify methods. The hooks are written so a future open-extension escape hatch (an `OP_CUSTOM` arm carrying a payload) could be slotted in without reshaping them; a `Union{Nothing,QSym}` payload field is deliberately *not* added now because it would fail the `CheckConcreteStructs` (`all_concrete`) gate.
+**Adding an operator role.** Add an `OP_*` enum arm, a constructor function and an `is_*` predicate, then a branch in each of the six hooks plus `adjoint`, `order_key`, the `numeric_operator` arm in each backend extension (`ext/`), and the printing/Latexify methods. The hooks are written so a future open-extension escape hatch (an `OP_CUSTOM` arm carrying a payload) could be slotted in without reshaping them; a `Union{Nothing,QSym}` payload field is deliberately *not* added now because it would fail the `CheckConcreteStructs` (`all_concrete`) gate.
 
 
 ## Hilbert spaces
@@ -378,20 +378,73 @@ The scope rides as a `SumScope` *argument* of the `Term`, not as metadata, becau
 
 ## Numeric conversion
 
-`to_numeric` maps symbolic operators to QuantumOpticsBase matrices:
-
-- **Simple basis:** Direct dispatch — `Destroy → destroy(b)`, `Transition → transition(b, i, j)`, etc.
-- **Composite basis:** Embeds the single-site matrix as a `LazyTensor`:
-  ```julia
-  op_num = to_numeric(op, b.bases[idx])
-  LazyTensor(b, [idx], (op_num,))
-  ```
-
-**`_to_number`** extracts plain Julia numbers from `Num`/`CNum` wrappers for numeric evaluation. Falls back to the symbolic value if it can't be unwrapped (for symbolic prefactors that haven't been substituted yet).
-
-**`_reduce_const`** reduces a fully-substituted coefficient part to a number. `Symbolics.value` handles numeric constants directly, but a part it leaves symbolic (for example `conj` of a complex literal, which SymbolicUtils does not fold) is compiled with `build_function` and evaluated. Because a `Number`-symtype parameter is held opaquely in the real slot, that real part can itself reduce to a `Complex`, so `_to_complex` recombines as `_reduce_const(real(x)) + im * _reduce_const(imag(x))` rather than `Complex(re, im)`.
-
-**`_lazy_one`** creates the identity operator. For simple bases it returns `one(b)` (dense identity). For composite bases it returns a `LazyTensor` identity rather than materializing the full Kronecker-product identity matrix.
+`to_numeric`/`numeric_average` are a layered, open-via-dispatch, multi-backend design. The
+backend-neutral core lives in `src/numeric/` and never names a concrete numeric type; the
+two backends (QuantumOpticsBase, QuantumToolbox) are Julia package extensions in `ext/`. The
+algebra itself depends only on the lightweight `QuantumInterface.jl` (it owns the
+`⊗`/`tensor`/`expect`/`basis` generics the symbolic side extends), so building Hilbert
+spaces and operators needs no numeric backend loaded.
+
+**Backends and hooks.** A backend is a zero-field singleton (`QuantumOpticsBackend`,
+`QuantumToolboxBackend`) so dispatch is static. The core talks to a backend only through the
+open generic hooks declared in `src/numeric/backend.jl` and defined per backend in `ext/`:
+`numeric_operator` (operator role to matrix), `numeric_basis` (Hilbert space + user `dims`
+to a backend basis), `numeric_subbasis`, `numeric_embed`, `numeric_identity`,
+`numeric_assemble` (static), `numeric_assemble_td` (time-dependent), and `numeric_expect`.
+`numeric_operator` is a closed `k === OP_*` value-branch on `op.kind` per backend with an
+open `Val(op.kind)` fallthrough (a separate throwing `numeric_operator(be, ::Val{K}, …)`
+method) so a downstream custom role can add a method without editing the package. The hooks
+are exported, which is the supported extension surface.
+
+**One backend-neutral core.** `NumericContext` bundles the backend singleton, the opaque
+basis, the operator/scalar substitution dicts, and the indexed `sites` map, with all fields
+concrete via type parameters. `_to_numeric_static` walks the canonical `QAdd` into a
+backend-neutral term list (`(c::ComplexF64, factors)`) and hands it to `numeric_assemble`;
+the core does no operator arithmetic. The kept coefficient machinery (`_to_complex` and its
+`_reduce_const`/`_fold_const` constant folder, the `_expand_parameter`/`_time_basis`
+splitting) is in `src/numeric/coeff.jl`. The `_to_complex(::Any)::ComplexF64` keystone stays
+exactly two methods (a third trips Julia's union-split budget; a test pins this).
+
+**Lazy and type-stable.** A single operator returns the bare eager matrix (so per-leaf uses
+stay simple), but a `QAdd` (any sum or product) returns a lazy operator. The key constraint
+is that a lazy sum's type must not depend on the runtime term count, and must be
+inference-stable so `@inferred(to_numeric(...))` holds. QuantumOptics uses the native
+vector-backed `LazySum`, constructed with the **5-argument** form
+`LazySum(ComplexF64, b, b, coeffs, ops)`: passing the concrete basis pins the basis type
+parameters, which the 2-argument `LazySum(coeffs, ops)` leaves abstract (runtime-`n`-stable
+but not `@inferred`-stable). QuantumToolbox's built-in `AddedOperator` is type-locked to a
+`Tuple`, so the QTB extension defines a small vector-backed `VecSum <:
+SciMLOperators.AbstractSciMLOperator` (abstract-eltype `Vector` fields, exactly like
+`LazySum`'s `Vector{AbstractOperator}`) wrapped in a `QobjEvo`; the result is one concrete
+`QobjEvo{Operator, Dims, VecSum{ComplexF64}}` for any term count and for static or
+time-dependent coefficients. Within a term, QuantumOptics embeds lazily (`LazyTensor`/
+`LazyProduct`) while QuantumToolbox embeds eagerly (a plain concrete `tensor`, because it
+warns on lazy tensors); the laziness lives only at the sum level on both backends.
+Materialize on demand with `dense`/`sparse` (QuantumOptics) or `sparse`/`concretize`
+(QuantumToolbox).
+
+**Time-dependent.** When `time_parameter` is non-empty the per-term coefficient becomes a
+`t -> ComplexF64` closure and `numeric_assemble_td` builds the backend's native
+time-dependent operator (`TimeDependentSum` / a `QobjEvo` whose `VecSum` coefficients are
+`ScalarOperator`s carrying the update function). These are callable at a time and feed
+ODE/SciML solvers directly, replacing the old `t -> op(t)` closure return.
+
+**Indexed path.** `src/numeric/indexed.jl` keeps the bound-index unroll (dependent-index
+detection, the `ne` diagonal filter, the `sub_op`/`sub_coef` dual substitution) verbatim,
+but routes each unrolled combination through the same hooks at the resolved site, so it emits
+entries of the same term list the static path assembles and gains the QTB backend for free.
+
+**QuantumObject 0-indexing.** QuantumToolbox is 0-indexed, so `Transition(i, j)` (SQA levels
+are 1-indexed, matching QuantumOptics) maps to `projection(N, i-1, j-1)` in the QTB
+extension.
+
+**Solver cost of laziness.** `expect`/`numeric_average` and `sesolve`/`schroedinger` apply
+`H` to the state directly and are unaffected. Only `mesolve` and other superoperator solvers
+pay a cost: QuantumToolbox builds the Liouvillian by tensoring `H` with identity, which a
+lazy `H` makes lazy (the "lazy tensor can hurt performance" warning). For fast `mesolve` on
+small or medium systems, concretize `H` first; keep it lazy only when materializing the
+Liouvillian is itself the bottleneck. The custom `VecSum` therefore implements `transpose`/
+`adjoint` returning a `VecSum` (not a wrapper), which the Liouvillian construction needs.
 
 
 ## Hermitian conjugation (operators.jl)