diff --git a/CLAUDE.md b/CLAUDE.md index ffed21aa..4881676e 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -14,8 +14,16 @@ A Julia package for symbolic manipulation of second-quantized operators used in src/SecondQuantizedAlgebra.jl # Main module, exports, imports src/types.jl # QField, QSym abstract hierarchy src/precompile.jl # PrecompileTools workload -src/average.jl # AvgFunc, average(), undo_average(), Hermitian-Real symtype -src/numeric.jl # to_numeric(), numeric_average() via QuantumOpticsBase +src/average.jl # AvgFunc/AvgSym, average(), undo_average(), Hermitian-Real symtype + +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 diff --git a/Changelog.md b/Changelog.md index 661b402c..482c60f8 100644 --- a/Changelog.md +++ b/Changelog.md @@ -4,14 +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, and multi-backend. This is a breaking release. -## [v0.10.0] +### 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, op_type)`. 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 backend hooks `numeric_operator`, `numeric_basis`, `numeric_subbasis`, `numeric_embed`, `numeric_identity`, `numeric_num_subsystems`, `numeric_assemble`, `numeric_assemble_td`, `numeric_materialize`, `numeric_expect`, and `numeric_backend` are exported. Downstream packages can implement another backend and add numeric support for existing operator roles; `OpKind` remains a closed symbolic-role set. + +### Changed (breaking) + +- `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). +- The time-dependent form (`time_parameter` non-empty) returns the backend's **native** time-dependent operator: a `TimeDependentSum` (QuantumOptics) or a `QobjEvo` (QuantumToolbox), both directly consumable by `mesolve`/`master_dynamic`/`sesolve`. It is no longer a `t -> op(t)` closure. +- Time-dependent conversion accepts only `op_type=nothing` or `identity`; eager `op_type` materializers are rejected because a time-varying operator cannot be materialized once during conversion. +- Product-space `dims` must have exactly one entry per symbolic subspace. Indexed numeric layouts validate that every physical slot is unique and in range instead of silently collapsing multiple sites onto a simple basis. ### Changed - Averages of provably Hermitian operators (`adjoint(A) == A`) now carry `symtype === Real` instead of `Number`. This gives a faster `simplify` path and makes `conj(⟨a'a⟩)` fold to `⟨a'a⟩` rather than an inert `conj(...)` wrapper; indexed sums and lifted time-dependent variables inherit the typing, which survives `substitute`. Resolves [#171](https://github.com/qojulia/SecondQuantizedAlgebra.jl/issues/171). +### Notes + +- For static conversion, the `op_type` contract from v0.9.2 is now backend-neutral: `to_numeric` assembles the operator lazily and materializes it once via `op_type`, so the return type depends only on `op_type`, not on the expression shape. The default consistently returns an eager backend operator (sparse on both bundled backends). Pass an explicit `op_type` for another eager representation (`dense` on QuantumOptics; `QuantumToolbox.to_sparse`/`to_dense` or `SciMLOperators.concretize` on QuantumToolbox), or `op_type=identity` for the natural lazy representation (`LazyTensor`/`LazyProduct`/`LazySum` / `QobjEvo` over `VecSum`). The lazy form is the internal assembly primitive and is what `numeric_average`/`expect` consume directly, so `LazyKet` states work without materializing. +- 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.4] ### Fixed diff --git a/Makefile b/Makefile index 7512c743..1e219c16 100644 --- a/Makefile +++ b/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")])' diff --git a/Project.toml b/Project.toml index da5cda54..eef8702e 100644 --- a/Project.toml +++ b/Project.toml @@ -4,23 +4,35 @@ 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] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678" +QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab" +SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" + +[extensions] +SecondQuantizedAlgebraQuantumOpticsBaseExt = "QuantumOpticsBase" +SecondQuantizedAlgebraQuantumToolboxExt = ["QuantumToolbox", "SciMLOperators", "LinearAlgebra"] + [compat] Combinatorics = "1.0.2" -FunctionWrappers = "1" Latexify = "0.13, 0.14, 0.15, 0.16" +LinearAlgebra = "1" 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" diff --git a/README.md b/README.md index 7718ebf3..176d9484 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,7 @@ The package provides: - Automatic commutation relations and canonical-form arithmetic - Normal ordering, simplification, and completeness expansion - Symbolic summations with automatic diagonal splitting -- Averaging and numeric conversion via QuantumOpticsBase +- Averaging and extensible numeric conversion via [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) or [QuantumToolbox.jl](https://github.com/qutip/QuantumToolbox.jl) - Extensible for custom operator types The code was refactored out of [QuantumCumulants.jl](https://github.com/qojulia/QuantumCumulants.jl). diff --git a/docs/src/API.md b/docs/src/API.md index f9bfd770..b24cd3bb 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -280,6 +280,29 @@ to_numeric numeric_average ``` +### Numeric backends + +Use `to_numeric(op, h, dims; backend=...)` for backend-portable code. The complete extension +contract and a minimal implementation are in [Adding a numeric backend](@ref +numeric-backend-interface). + +```@docs +NumericBackend +QuantumOpticsBackend +QuantumToolboxBackend +numeric_operator +numeric_basis +numeric_subbasis +numeric_embed +numeric_identity +numeric_num_subsystems +numeric_assemble +numeric_assemble_td +numeric_materialize +numeric_expect +numeric_backend +``` + ## [Symbolic Summations](@id API: Sums) diff --git a/docs/src/devdocs.md b/docs/src/devdocs.md index 63c5b795..f49f29f5 100644 --- a/docs/src/devdocs.md +++ b/docs/src/devdocs.md @@ -116,7 +116,7 @@ leaf fast path sums both slots. **Site families inside `_site_compare`.** The comparator first orders by `space_index`, so operators are grouped by the subspace they act on before anything else is consulted (operators on different subspaces commute, so this reordering is always safe and gives the more natural per-subspace canonical form). Within a subspace, 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. A subspace carries a single Hilbert-space type, so same-`space_index` operators always share a family; the `kind` integer fallback preserves the existing values `OP_DESTROY=0 < … < OP_MOMENTUM=6` and appends `OP_COLLECTIVE_TRANSITION=7`. It only distinguishes the pathological case of two operators sharing a `space_index` across unrelated simple spaces. -**Adding an operator role.** Add an `OP_*` enum arm, a constructor function and an `is_*` predicate, then add any non-default branches required in 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 add any non-default branches required in 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 @@ -395,28 +395,89 @@ The scope rides as a `SumScope` *argument* of the `Term`, not as metadata, becau ## Numeric conversion -`to_numeric` maps symbolic operators to QuantumOpticsBase matrices. A term is -assembled in a natural lazy representation and materialized to a concrete -operator exactly once, at the top, via `op_type` (default `sparse`). This keeps -the return type independent of the shape of the expression (a bare operator, a -product, and a sum all return the same `op_type`), while still avoiding -per-factor conversions during assembly. +The backend-neutral core in `src/numeric/` never names a concrete numeric type. +QuantumOpticsBase and QuantumToolbox support lives in package extensions under `ext/`, so +the symbolic algebra can be loaded without either package. -- **Simple basis:** Direct dispatch, e.g. `Destroy → destroy(b)`, `Transition → transition(b, i, j)`. -- **Composite basis:** The internal `_embed` embeds the single-site matrix as a `LazyTensor`: - ```julia - op_num = _embed(op, b.bases[idx]) - LazyTensor(b, [idx], (op_num,)) - ``` - `_embed` is the composable building block for products and sums. The public - `to_numeric` methods wrap the assembled result in `op_type`, so the positional - forms return `sparse` and `op_type=identity` returns the lazy form as-is. +### [Adding a numeric backend](@id numeric-backend-interface) -**`_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). +Define a concrete singleton subtype of `NumericBackend` and extend the exported hooks. A +static-only backend needs the following methods: -**`_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)`. +| Hook | Responsibility | +|:--|:--| +| `numeric_basis(be, h, dims)` | Build the full backend basis or dimension descriptor. | +| `numeric_num_subsystems(be, basis)` | Return the number of tensor factors. | +| `numeric_subbasis(be, basis, slot)` | Select one tensor factor. | +| `numeric_operator(be, op, subbasis)` | Convert a symbolic leaf. | +| `numeric_embed(be, basis, slot, leaf)` | Place a leaf in the full space. | +| `numeric_identity(be, basis)` | Build the full-space identity. | +| `numeric_assemble(be, basis, terms)` | Combine `(ComplexF64, factors)` terms. | +| `numeric_materialize(be, assembled, op_type)` | Select the public representation. | -**`_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. +`numeric_assemble` may return a lazy object. `numeric_materialize` must give `op_type` these +backend-neutral meanings: + +- `nothing`: return the backend's ordinary eager operator; +- `identity`: return `assembled` unchanged; +- any other callable: a backend-defined explicit conversion. + +The ordinary eager type need not be sparse; that is only the choice made by both bundled +backends. Materialization happens once, after the whole expression is assembled, so a leaf, +product, and sum have the same representation for a given `op_type`. + +A minimal static backend looks like this: + +```julia +import SecondQuantizedAlgebra as SQA + +struct MyBackend <: SQA.NumericBackend end + +SQA.numeric_basis(::MyBackend, h::SQA.HilbertSpace, dims) = error("build the basis") +SQA.numeric_num_subsystems(::MyBackend, basis) = error("count subsystems") +SQA.numeric_subbasis(::MyBackend, basis, slot::Int) = error("select a subsystem") +SQA.numeric_operator(::MyBackend, op::SQA.Op, subbasis) = error("construct a leaf") +SQA.numeric_embed(::MyBackend, basis, slot::Int, leaf) = error("embed a leaf") +SQA.numeric_identity(::MyBackend, basis) = error("construct the identity") +SQA.numeric_assemble(::MyBackend, basis, terms) = error("assemble terms") +SQA.numeric_materialize(::MyBackend, assembled, ::Nothing) = + error("make an eager operator") +SQA.numeric_materialize(::MyBackend, assembled, ::typeof(identity)) = assembled +``` + +Optional capabilities add three small groups of methods: + +| Capability | Additional hooks | +|:--|:--| +| time-dependent conversion | `numeric_assemble_td` | +| conversion from backend states | `numeric_backend(state)`, `numeric_basis(state)` | +| `numeric_average` / `expect` | the state hooks plus `numeric_expect` | + +Third-party backends are passed explicitly as `backend=MyBackend()`; automatic discovery is +limited to the bundled extensions. `numeric_operator` covers the existing closed `OpKind` +roles. Packages can support a new backend/basis combination, but cannot currently define a +new symbolic operator role through this interface. + +### Implementation notes + +- `NumericContext` holds the concrete backend singleton, opaque basis, substitutions, and + indexed-site map. The core emits `(ComplexF64, factors)` terms and performs no backend + operator arithmetic. +- The internal assembled type must not depend on the number of terms. QuantumOptics uses + the five-argument, vector-backed `LazySum`; QuantumToolbox uses the extension's + vector-backed `VecSum` inside `QobjEvo`. Tests enforce inference stability. +- `numeric_average` calls `numeric_expect` on the lazy assembly directly. It therefore does + not build an eager matrix merely to compute an expectation value. +- A non-empty `time_parameter` uses `numeric_assemble_td` and returns the native + `TimeDependentSum` or `QobjEvo`. Only `op_type=nothing` and `identity` are accepted because + a time-varying operator cannot be materialized once during conversion. +- Indexed conversion validates sites, unrolls sums, and emits the same term format as the + non-indexed path. A site must be unique and within the backend basis. +- QuantumToolbox levels are zero-based internally, so symbolic `Transition(i, j)` maps to + `projection(N, i-1, j-1)`. +- Lazy Hamiltonians are useful when materializing a large tensor product is prohibitive, + but can slow superoperator construction. Prefer the eager default for small and medium + solver problems; request `op_type=identity` deliberately. ## Hermitian conjugation (operators.jl) diff --git a/docs/src/implementation.md b/docs/src/implementation.md index 62c3a812..ce0e07d6 100644 --- a/docs/src/implementation.md +++ b/docs/src/implementation.md @@ -202,7 +202,17 @@ Averaged expressions also preserve summation metadata when working with indexed ## Numeric conversion -Operators can be converted to numeric representations using [QuantumOpticsBase.jl](https://github.com/qojulia/QuantumOpticsBase.jl). +Numeric conversion supports [QuantumOpticsBase.jl](https://github.com/qojulia/QuantumOpticsBase.jl) +and [QuantumToolbox.jl](https://github.com/qutip/QuantumToolbox.jl). With one bundled backend +loaded it is inferred; if both are loaded, select one explicitly: + +```julia +to_numeric(op, h, dims; backend = QuantumOpticsBackend()) +to_numeric(op, h, dims; backend = QuantumToolboxBackend()) +``` + +This Hilbert-space form is the portable API. Backend-native basis or dimension forms remain +available as conveniences. ### Direct conversion @@ -230,8 +240,9 @@ H_num = to_numeric(H, b; parameter = Dict(Δ => 2.0)) nothing # hide ``` -For time-dependent parameters, pass numbers or functions in `time_parameter`. -The result is a callable `t -> op(t)`. +For time-dependent parameters, pass numbers or functions in `time_parameter`. The result is +the backend's native time-dependent operator (`TimeDependentSum` or `QobjEvo`), ready for +that backend's solvers. ### Numeric averages @@ -265,7 +276,10 @@ nothing # hide ### Composite systems -For product spaces, [`to_numeric`](@ref) returns a concrete `sparse` operator by default, independent of the shape of the expression. Pass `op_type=identity` to get the `LazyTensor`/`LazySum` representation instead, which avoids materializing the full Kronecker product and is preferable for large tensor products where an operator is local to a few subsystems: +For static expressions, [`to_numeric`](@ref) returns an eager backend operator by default; +both bundled backends use sparse storage. This is independent of whether the expression is +a single operator, product, or sum. Pass `op_type=identity` to retain the backend's lazy +assembly, which can avoid a large materialized tensor product: ```@example numeric-composite using SecondQuantizedAlgebra, QuantumOpticsBase @@ -280,8 +294,8 @@ bf = FockBasis(10) bn = NLevelBasis(3) bc = bf ⊗ bn -a_num = to_numeric(a, bc) # sparse Operator (default) -a_lazy = to_numeric(a, bc; op_type = identity) # LazyTensor +a_num = to_numeric(a, bc) # eager sparse Operator +a_lazy = to_numeric(a, bc; op_type = identity) # lazy assembly s_num = to_numeric(s(:a, :c), bc) nothing # hide ``` diff --git a/docs/src/index.md b/docs/src/index.md index 8189b769..0e7bad07 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -15,7 +15,7 @@ The package provides: - Eager canonical-form arithmetic: every `*` applies (anti-)commutation, local algebraic identities, and `NLevelSpace` completeness in one pass, so the return type is always a canonical `QAdd`. - Explicit pipeline functions when you want piecewise control: `normal_order`, `simplify`, `commutator`, `anticommutator`, `expand`, `expand_completeness`. - Symbolic summations via `Index` and `Σ` for indexed families, with automatic diagonal splitting and `assume_distinct_index` for free-index constraints. -- Averaging to symbolic scalars via `average` / `undo_average`, and numeric conversion to QuantumOpticsBase matrices via `to_numeric` / `numeric_average`. +- Averaging to symbolic scalars via `average` / `undo_average`, and extensible numeric conversion via QuantumOpticsBase or QuantumToolbox with `to_numeric` / `numeric_average`. - Hermitian conjugation across mixed operator + symbolic expressions via `qadjoint` (aliased as `qconj`, `dagger`) and the average-aware `inner_adjoint`. - Extensible for custom operator types via five small hooks — see the developer docs. diff --git a/ext/SecondQuantizedAlgebraQuantumOpticsBaseExt.jl b/ext/SecondQuantizedAlgebraQuantumOpticsBaseExt.jl new file mode 100644 index 00000000..a7cc0f8c --- /dev/null +++ b/ext/SecondQuantizedAlgebraQuantumOpticsBaseExt.jl @@ -0,0 +1,164 @@ +module SecondQuantizedAlgebraQuantumOpticsBaseExt + +# QuantumOpticsBase backend for `to_numeric`/`numeric_average`. Provides the operator-to- +# matrix map, basis construction, lazy embedding, and vector-backed lazy assembly. Loaded +# automatically by `using QuantumOpticsBase`. + +import SecondQuantizedAlgebra as SQA +using SecondQuantizedAlgebra: QuantumOpticsBackend, Op +import QuantumOpticsBase as QOB + +const QOBState = Union{QOB.StateVector, QOB.AbstractOperator} + +# --- operator -> matrix (closed value-branch ladder + open Val fallthrough) ------------- + +@inline function SQA.numeric_operator(be::QuantumOpticsBackend, op::Op, b::QOB.Basis) + k = op.kind + k === SQA.OP_DESTROY && return QOB.destroy(b) + k === SQA.OP_CREATE && return QOB.create(b) + k === SQA.OP_POSITION && return (QOB.destroy(b) + QOB.create(b)) / sqrt(2) + k === SQA.OP_MOMENTUM && return im * (QOB.create(b) - QOB.destroy(b)) / sqrt(2) + k === SQA.OP_TRANSITION && return QOB.transition(b, Int(op.l1), Int(op.l2)) + k === SQA.OP_PAULI && return _pauli_matrix(op, b) + k === SQA.OP_SPIN && return _spin_matrix(op, b) + return SQA.numeric_operator(be, Val(k), op, b) +end + +# Separate, throwing extension point for custom kinds (`::Val` is `::Val{K} where K`, so it +# must NOT share a signature with the main `op::Op` method above). +SQA.numeric_operator(::QuantumOpticsBackend, ::Val{K}, op::Op, b::QOB.Basis) where {K} = + throw(ArgumentError("no QuantumOpticsBase numeric_operator for role $K on $(typeof(b))")) + +# CollectiveTransition lives on a `ManyBodyBasis` over an `NLevelBasis`: the single-body +# transition is promoted to the symmetric many-body operator via `manybodyoperator`. This +# method wins over the generic `op::Op, b::QOB.Basis` ladder for a `ManyBodyBasis`, so a +# non-collective role on such a basis (or a collective role elsewhere) throws cleanly. +function SQA.numeric_operator(::QuantumOpticsBackend, op::Op, b::QOB.ManyBodyBasis) + SQA.is_collective_transition(op) || + throw(ArgumentError("Op kind $(op.kind) does not act on a ManyBodyBasis")) + onebody = b.onebodybasis + onebody isa QOB.NLevelBasis || throw( + ArgumentError( + "CollectiveTransition requires a ManyBodyBasis with an NLevelBasis one-body basis; got $(typeof(onebody))", + ), + ) + onebody_op = QOB.transition(onebody, Int(op.l1), Int(op.l2)) + return QOB.manybodyoperator(b, onebody_op) +end + +function _pauli_matrix(op::Op, b::QOB.SpinBasis) + b.spinnumber == 1 // 2 || + throw(ArgumentError("Pauli operators require SpinBasis(1//2), got SpinBasis($(b.spinnumber))")) + op.l1 == 1 && return QOB.sigmax(b) + op.l1 == 2 && return QOB.sigmay(b) + return QOB.sigmaz(b) +end + +function _spin_matrix(op::Op, b::QOB.SpinBasis) + op.l1 == 1 && return 0.5 * QOB.sigmax(b) + op.l1 == 2 && return 0.5 * QOB.sigmay(b) + return 0.5 * QOB.sigmaz(b) +end + +# --- basis construction (HilbertSpace form) -------------------------------------------- + +SQA.numeric_basis(::QuantumOpticsBackend, ::SQA.FockSpace, N) = QOB.FockBasis(Int(N)) +SQA.numeric_basis(::QuantumOpticsBackend, h::SQA.NLevelSpace, _) = QOB.NLevelBasis(h.n) +SQA.numeric_basis(::QuantumOpticsBackend, ::SQA.PauliSpace, _) = QOB.SpinBasis(1 // 2) +SQA.numeric_basis(::QuantumOpticsBackend, ::SQA.SpinSpace, s) = QOB.SpinBasis(s) +SQA.numeric_basis(::QuantumOpticsBackend, ::SQA.PhaseSpace, N) = QOB.FockBasis(Int(N)) +function SQA.numeric_basis(be::QuantumOpticsBackend, h::SQA.ProductSpace, dims) + SQA._check_product_dims(h, dims) + subs = h.spaces + return QOB.tensor(ntuple(i -> SQA.numeric_basis(be, subs[i], dims[i]), length(subs))...) +end + +# --- subsystem basis / embedding / identity -------------------------------------------- + +function SQA.numeric_subbasis(::QuantumOpticsBackend, b::QOB.Basis, slot::Int) + slot == 1 || throw(ArgumentError("simple QuantumOptics basis has no subsystem slot $slot")) + return b +end +function SQA.numeric_subbasis(::QuantumOpticsBackend, b::QOB.CompositeBasis, slot::Int) + 1 <= slot <= length(b.bases) || throw( + ArgumentError( + "QuantumOptics composite basis has $(length(b.bases)) subsystems, not slot $slot", + ), + ) + return b.bases[slot] +end + +SQA.numeric_num_subsystems(::QuantumOpticsBackend, ::QOB.Basis) = 1 +SQA.numeric_num_subsystems(::QuantumOpticsBackend, b::QOB.CompositeBasis) = length(b.bases) + +SQA.numeric_embed(::QuantumOpticsBackend, b::QOB.Basis, slot::Int, m) = m +SQA.numeric_embed(::QuantumOpticsBackend, b::QOB.CompositeBasis, slot::Int, m) = + QOB.LazyTensor(b, [slot], (m,)) + +SQA.numeric_identity(::QuantumOpticsBackend, b::QOB.Basis) = one(b) +SQA.numeric_identity(::QuantumOpticsBackend, b::QOB.CompositeBasis) = + QOB.LazyTensor(b, collect(1:length(b.bases)), Tuple(one(bi) for bi in b.bases)) + +# --- vector-backed lazy assembly ------------------------------------------------------- + +# Static: one concrete `LazySum` type for any term count. Pass the concrete basis to the +# 5-arg constructor so `BL`/`BR` are pinned and the result is `@inferred`-stable (the 2-arg +# `LazySum(coeffs, ops)` is only runtime-n-stable, not inference-stable). +function SQA.numeric_assemble(be::QuantumOpticsBackend, b, terms) + coeffs = ComplexF64[] + ops = QOB.AbstractOperator[] + for (c, factors) in terms + push!(coeffs, c) + push!(ops, _fuse(be, b, factors)) + end + if isempty(ops) + push!(coeffs, 0.0im) + push!(ops, SQA.numeric_identity(be, b)) + end + return QOB.LazySum(ComplexF64, b, b, coeffs, ops) +end + +# Time-dependent: native TimeDependentSum, vector-backed (same n-stability story). The +# per-term operators are materialised `sparse` so the TD operator is solver-friendly in the +# hot ODE loop (the lazy form is kept only for the static `op_type = identity` opt-in). +function SQA.numeric_assemble_td(be::QuantumOpticsBackend, b, td_terms) + coeffs = Function[] + ops = QOB.AbstractOperator[] + for (c, factors) in td_terms + push!( + coeffs, c isa Function ? c : ( + let c = c + t -> c + end + ) + ) + push!(ops, QOB.sparse(_fuse(be, b, factors))) + end + return QOB.TimeDependentSum(ComplexF64, b, b, coeffs, ops) +end + +# Per-term operator: identity for an empty product, the bare leaf for one factor, else the +# product of the embedded factors. `*` combines disjoint-slot `LazyTensor`s into a single +# `LazyTensor` (rather than a `LazyProduct`), which keeps the lazy form usable with a +# `LazyKet` state and matches the eager product of the same operators. +function _fuse(be::QuantumOpticsBackend, b, factors) + isempty(factors) && return SQA.numeric_identity(be, b) + length(factors) == 1 && return factors[1] + return foldl(*, factors) +end + +# --- materialization (op_type applied once at the top of `to_numeric`) ------------------ + +# Default (`nothing`) materialises `sparse`; any explicit `op_type` (`sparse`, `dense`, +# `identity`) is applied as a plain callable. +SQA.numeric_materialize(::QuantumOpticsBackend, op, ::Nothing) = QOB.sparse(op) +SQA.numeric_materialize(::QuantumOpticsBackend, op, op_type) = op_type(op) + +# --- expectation ----------------------------------------------------------------------- + +SQA.numeric_expect(::QuantumOpticsBackend, numop, state::QOBState) = + ComplexF64(QOB.expect(numop, state)) +SQA.numeric_expect(::QuantumOpticsBackend, numop, states::AbstractVector) = + QOB.expect(QOB.sparse(numop), states) + +end # module diff --git a/ext/SecondQuantizedAlgebraQuantumToolboxExt.jl b/ext/SecondQuantizedAlgebraQuantumToolboxExt.jl new file mode 100644 index 00000000..86ec29a0 --- /dev/null +++ b/ext/SecondQuantizedAlgebraQuantumToolboxExt.jl @@ -0,0 +1,304 @@ +module SecondQuantizedAlgebraQuantumToolboxExt + +# QuantumToolbox backend for `to_numeric`/`numeric_average`. Per-term operators are eager +# concrete `QuantumObject`s; the laziness lives in a vector-backed `VecSum` (a small custom +# `SciMLOperators.AbstractSciMLOperator`) wrapped in a `QobjEvo`, so the result is one +# concrete type for any term count and for static or time-dependent coefficients. Loaded by +# `using QuantumToolbox` (which also loads SciMLOperators, the second extension trigger). + +import SecondQuantizedAlgebra as SQA +using SecondQuantizedAlgebra: QuantumToolboxBackend, Op +import QuantumToolbox as QTB +import SciMLOperators as SO +import SymbolicUtils: BasicSymbolic +# `mul!` is imported from its owner `LinearAlgebra` (a weakdep trigger of this extension); +# `transpose`/`adjoint` come through Base. +import LinearAlgebra: mul! + +# QuantumToolbox represents a "basis" as integer dims: an `Int` for a simple space, a +# `Vector{Int}` for a composite one. +const QTBDims = Union{Integer, AbstractVector{<:Integer}} + +# ======================================================================================= +# Vector-backed lazy sum +# +# SciMLOperators' built-in `AddedOperator` is type-locked to a `Tuple`, so a lazy sum built +# from it is type-unstable for a runtime term count. `VecSum` stores the per-term coefficient +# operators and operators in `Vector`s instead, so `VecSum{ComplexF64}` is ONE concrete type +# for any term count (the heterogeneous elements live behind the abstract `Vector` eltype, +# one dynamic dispatch per term in the apply loop, exactly like QuantumOptics' own +# `Vector{AbstractOperator}`-backed `LazySum`). The field eltypes are abstract on purpose: a +# `ScalarOperator`'s full type encodes its update function, so constant and each distinct +# time-dependent coefficient are different concrete types. +# ======================================================================================= + +struct VecSum{T} <: SO.AbstractSciMLOperator{T} + coeffs::Vector{SO.ScalarOperator} + ops::Vector{SO.AbstractSciMLOperator} +end +VecSum(coeffs, ops) = VecSum{ComplexF64}(coeffs, ops) + +_cv(c) = convert(Number, c) + +Base.size(L::VecSum) = size(L.ops[1]) +Base.size(L::VecSum, i::Int) = size(L.ops[1], i) +Base.eltype(::VecSum{T}) where {T} = T +SO.issquare(::VecSum) = true +SO.isconstant(L::VecSum) = all(SO.isconstant, L.coeffs) +SO.update_coefficients!(L::VecSum, u, p, t) = + (foreach(c -> SO.update_coefficients!(c, u, p, t), L.coeffs); L) + +function mul!(w::AbstractVecOrMat, L::VecSum, v::AbstractVecOrMat) + mul!(w, L.ops[1], v, _cv(L.coeffs[1]), false) + for i in 2:length(L.ops) + mul!(w, L.ops[i], v, _cv(L.coeffs[i]), true) + end + return w +end +function mul!(w::AbstractVecOrMat, L::VecSum, v::AbstractVecOrMat, α, β) + mul!(w, L.ops[1], v, _cv(L.coeffs[1]) * α, β) + for i in 2:length(L.ops) + mul!(w, L.ops[i], v, _cv(L.coeffs[i]) * α, true) + end + return w +end + +# Allocating apply (needed by `expect`, which routes through `dot(ψ, L, ψ) = dot(ψ, L*ψ)`). +function Base.:*(L::VecSum, v::AbstractVector) + w = similar(v, promote_type(eltype(L), eltype(v)), size(L, 1)) + return mul!(w, L, v) +end +function Base.:*(L::VecSum, v::AbstractMatrix) + w = similar(v, promote_type(eltype(L), eltype(v)), size(L, 1), size(v, 2)) + return mul!(w, L, v) +end + +# In-place call forms used by the SciML/QuantumToolbox solver stepping. +(L::VecSum)(w, v, u, p, t) = (SO.update_coefficients!(L, u, p, t); mul!(w, L, v)) +(L::VecSum)(w, v, u, p, t, α, β) = (SO.update_coefficients!(L, u, p, t); mul!(w, L, v, α, β)) + +SO.concretize(L::VecSum) = sum(_cv(L.coeffs[i]) * SO.concretize(L.ops[i]) for i in eachindex(L.ops)) +Base.convert(::Type{AbstractMatrix}, L::VecSum) = SO.concretize(L) +# `::AbstractVecOrMat` avoids an ambiguity with the SciMLOperators fallback. +SO.cache_operator(L::VecSum, ::AbstractVecOrMat) = L + +# transpose/adjoint MUST return a VecSum (not be wrapped in a TransposedOperator), or the +# mesolve Liouvillian (which tensors transpose(H) with identity) fails its cache reshape. +Base.transpose(L::VecSum{T}) where {T} = + VecSum{T}(L.coeffs, SO.AbstractSciMLOperator[transpose(o) for o in L.ops]) +Base.adjoint(L::VecSum{T}) where {T} = + VecSum{T}( + SO.AbstractSciMLOperator[adjoint(c) for c in L.coeffs], + SO.AbstractSciMLOperator[adjoint(o) for o in L.ops], +) + +# ======================================================================================= +# operator -> QuantumObject (closed value-branch ladder + open Val fallthrough) +# ======================================================================================= + +@inline function SQA.numeric_operator(be::QuantumToolboxBackend, op::Op, N::Integer) + k = op.kind + k === SQA.OP_DESTROY && return QTB.destroy(N) + k === SQA.OP_CREATE && return QTB.create(N) + k === SQA.OP_POSITION && return QTB.position(N) + k === SQA.OP_MOMENTUM && return QTB.momentum(N) + # QuantumToolbox is 0-indexed; SQA `Transition` levels are 1-indexed (as in QuantumOptics). + k === SQA.OP_TRANSITION && return QTB.projection(N, Int(op.l1) - 1, Int(op.l2) - 1) + k === SQA.OP_PAULI && return _qtb_pauli(op, N) + k === SQA.OP_SPIN && return QTB.jmat((N - 1) // 2, _axis(op)) + return SQA.numeric_operator(be, Val(k), op, N) +end + +SQA.numeric_operator(::QuantumToolboxBackend, ::Val{K}, op::Op, N::Integer) where {K} = + throw(ArgumentError("no QuantumToolbox numeric_operator for role $K on dim $N")) + +function _qtb_pauli(op::Op, N::Integer) + N == 2 || throw(ArgumentError("Pauli operators require dim 2, got $N")) + op.l1 == 1 && return QTB.sigmax() + op.l1 == 2 && return QTB.sigmay() + return QTB.sigmaz() +end + +_axis(op::Op) = op.l1 == 1 ? :x : op.l1 == 2 ? :y : :z + +# ======================================================================================= +# basis (integer dims) / subsystem / embedding / identity +# ======================================================================================= + +SQA.numeric_basis(::QuantumToolboxBackend, ::SQA.FockSpace, N) = Int(N) + 1 +SQA.numeric_basis(::QuantumToolboxBackend, h::SQA.NLevelSpace, _) = Int(h.n) +SQA.numeric_basis(::QuantumToolboxBackend, ::SQA.PauliSpace, _) = 2 +SQA.numeric_basis(::QuantumToolboxBackend, ::SQA.SpinSpace, s) = Int(2s + 1) +SQA.numeric_basis(::QuantumToolboxBackend, ::SQA.PhaseSpace, N) = Int(N) + 1 +function SQA.numeric_basis(be::QuantumToolboxBackend, h::SQA.ProductSpace, dims) + SQA._check_product_dims(h, dims) + return Int[SQA.numeric_basis(be, h.spaces[i], dims[i]) for i in eachindex(h.spaces)] +end + +function SQA.numeric_subbasis(::QuantumToolboxBackend, N::Integer, slot::Int) + slot == 1 || throw(ArgumentError("simple QuantumToolbox dimension has no subsystem slot $slot")) + return N +end +function SQA.numeric_subbasis(::QuantumToolboxBackend, ds::AbstractVector, slot::Int) + 1 <= slot <= length(ds) || throw( + ArgumentError( + "QuantumToolbox dimensions have $(length(ds)) subsystems, not slot $slot", + ), + ) + return ds[slot] +end + +SQA.numeric_num_subsystems(::QuantumToolboxBackend, ::Integer) = 1 +SQA.numeric_num_subsystems(::QuantumToolboxBackend, ds::AbstractVector) = length(ds) + +# Within-term embedding is EAGER concrete (QuantumToolbox warns on lazy tensor); the lazy sum +# is the only lazy layer. +SQA.numeric_embed(::QuantumToolboxBackend, N::Integer, slot::Int, m) = m +function SQA.numeric_embed(::QuantumToolboxBackend, ds::AbstractVector, slot::Int, m) + length(ds) == 1 && return m + return QTB.tensor((i == slot ? m : QTB.qeye(ds[i]) for i in eachindex(ds))...) +end + +SQA.numeric_identity(::QuantumToolboxBackend, N::Integer) = QTB.qeye(N) +function SQA.numeric_identity(::QuantumToolboxBackend, ds::AbstractVector) + length(ds) == 1 && return QTB.qeye(ds[1]) + return QTB.tensor((QTB.qeye(d) for d in ds)...) +end + +# ======================================================================================= +# vector-backed lazy assembly (static + time-dependent share the one VecSum type) +# ======================================================================================= + +function SQA.numeric_assemble(be::QuantumToolboxBackend, basis, terms) + refid = SQA.numeric_identity(be, basis) + coeffs = SO.ScalarOperator[] + ops = SO.AbstractSciMLOperator[] + for (c, factors) in terms + push!(coeffs, SO.ScalarOperator(ComplexF64(c))) + push!(ops, SO.MatrixOperator(_fuse(be, basis, factors).data)) + end + if isempty(ops) + push!(coeffs, SO.ScalarOperator(0.0im)) + push!(ops, SO.MatrixOperator(refid.data)) + end + return QTB.QobjEvo(VecSum(coeffs, ops), QTB.Operator(), refid.dimensions) +end + +function SQA.numeric_assemble_td(be::QuantumToolboxBackend, basis, td_terms) + refid = SQA.numeric_identity(be, basis) + coeffs = SO.ScalarOperator[] + ops = SO.AbstractSciMLOperator[] + for (c, factors) in td_terms + push!(coeffs, _td_scalar(c)) + push!(ops, SO.MatrixOperator(_fuse(be, basis, factors).data)) + end + return QTB.QobjEvo(VecSum(coeffs, ops), QTB.Operator(), refid.dimensions) +end + +# SciML update signature is the four-arg `(a, u, p, t) -> c`, not the two-arg `(p, t)`. +_td_scalar(c::Function) = SO.ScalarOperator(0.0im, (a, u, p, t) -> ComplexF64(c(t))) +_td_scalar(c) = SO.ScalarOperator(ComplexF64(c)) + +function _fuse(be::QuantumToolboxBackend, basis, factors) + isempty(factors) && return SQA.numeric_identity(be, basis) + length(factors) == 1 && return factors[1] + return foldl(*, factors) +end + +# ======================================================================================= +# materialization (op_type applied once at the top of `to_numeric`) +# ======================================================================================= + +# The default (`nothing`) has the same semantics as QuantumOptics: return an eager sparse +# backend operator. `op_type = identity` is the explicit opt-in to the lazily assembled +# `QobjEvo`/`VecSum`. Other callables are applied to an eager `QuantumObject`; +# `SciMLOperators.concretize` is special-cased to return its raw matrix as advertised. +SQA.numeric_materialize(::QuantumToolboxBackend, op, ::Nothing) = QTB.to_sparse(_qtb_eager(op)) +function SQA.numeric_materialize(::QuantumToolboxBackend, op, op_type) + op_type === identity && return op + eager = _qtb_eager(op) + op_type === SO.concretize && return eager.data + return op_type(eager) +end + +# Reduce the assembled object to an eager `QuantumObject`: concretize the lazy `VecSum` +# behind a `QobjEvo`; a bare scalar path is already an eager `QuantumObject`. +_qtb_eager(op::QTB.QuantumObjectEvolution) = QTB.Qobj(SO.concretize(op.data); dims = op.dimensions) +_qtb_eager(op::QTB.QuantumObject) = op + +# ======================================================================================= +# expectation + backend resolution + state/dims convenience +# ======================================================================================= + +SQA.numeric_expect(::QuantumToolboxBackend, numop, state::QTB.QuantumObject) = + ComplexF64(QTB.expect(numop, state)) +SQA.numeric_expect(::QuantumToolboxBackend, numop, states::AbstractVector) = + ComplexF64[ComplexF64(QTB.expect(numop, s)) for s in states] + +SQA.numeric_backend(::QTB.AbstractQuantumObject) = QuantumToolboxBackend() +function SQA.numeric_basis(o::QTB.AbstractQuantumObject) + ds = collect(Int, o.dims[1]) + return length(ds) == 1 ? ds[1] : ds +end + +# Positional dims form (the QuantumToolbox analog of `to_numeric(op, basis, d)`). +SQA.to_numeric(op::Op, dims::QTBDims, d::AbstractDict{<:SQA.QSym} = SQA._NO_SUBS) = + SQA.numeric_materialize( + QuantumToolboxBackend(), + SQA._to_numeric_static( + SQA._single_qadd(SQA._CNUM_ONE, Op[op]), + SQA.NumericContext(QuantumToolboxBackend(), dims, d), + ), + nothing, +) +SQA.to_numeric(q::SQA.QAdd, dims::QTBDims, d::AbstractDict{<:SQA.QSym} = SQA._NO_SUBS) = + SQA.numeric_materialize( + QuantumToolboxBackend(), + SQA._to_numeric_static(q, SQA.NumericContext(QuantumToolboxBackend(), dims, d)), + nothing, +) +SQA.to_numeric(x::Number, dims::QTBDims, ::AbstractDict{<:SQA.QSym} = SQA._NO_SUBS) = + SQA.numeric_materialize( + QuantumToolboxBackend(), + SQA._to_complex(x) * SQA.numeric_identity(QuantumToolboxBackend(), dims), + nothing, +) + +# Indexed positional dims form (the QTB analog of `to_numeric(q, basis, sites, …)`). +SQA.to_numeric( + q::SQA.QAdd, dims::QTBDims, + sites::AbstractDict{Int, Vector{Int}}, + d::AbstractDict{<:SQA.QSym} = SQA._NO_SUBS, + scalar_subs::AbstractDict = SQA._NO_SCALAR_SUBS, +) = SQA._to_numeric_indexed(QuantumToolboxBackend(), dims, q, sites, d, scalar_subs) + +# Keyword dims form. +function SQA.to_numeric( + op::SQA.QField, dims::QTBDims; + parameter = Dict(), time_parameter = Dict(), + operators = Dict{SQA.QSym, Any}(), adjoint_ops = true, op_type = nothing, + ) + return SQA._to_numeric_kw(QuantumToolboxBackend(), op, dims; parameter, time_parameter, operators, adjoint_ops, op_type) +end +function SQA.to_numeric( + x::Union{Number, BasicSymbolic}, dims::QTBDims; + parameter = Dict(), time_parameter = Dict(), + operators = Dict{SQA.QSym, Any}(), adjoint_ops = true, op_type = nothing, + ) + return SQA._to_numeric_kw(QuantumToolboxBackend(), x, dims; parameter, time_parameter, operators, adjoint_ops, op_type) +end + +# Vector of operators on the direct QTB-dims form (mirrors the `Basis` method in api.jl). +SQA.to_numeric(ops::AbstractVector, dims::QTBDims; kwargs...) = + [SQA.to_numeric(op, dims; kwargs...) for op in ops] + +# State convenience (dims derived from the state). +SQA.to_numeric(op, state::QTB.QuantumObject; kwargs...) = SQA.to_numeric(op, SQA.numeric_basis(state); kwargs...) +SQA.to_numeric(op::SQA.QField, state::QTB.QuantumObject, d::AbstractDict{<:SQA.QSym} = SQA._NO_SUBS) = + SQA.to_numeric(op, SQA.numeric_basis(state), d) +SQA.to_numeric(x::Number, state::QTB.QuantumObject) = SQA.to_numeric(x, SQA.numeric_basis(state)) + +SQA.numeric_average(op, state::QTB.QuantumObject, d::AbstractDict{<:SQA.QSym} = SQA._NO_SUBS) = + SQA._numeric_average(op, state, d) + +end # module diff --git a/src/SecondQuantizedAlgebra.jl b/src/SecondQuantizedAlgebra.jl index c5dc3db3..542e4896 100644 --- a/src/SecondQuantizedAlgebra.jl +++ b/src/SecondQuantizedAlgebra.jl @@ -4,13 +4,10 @@ using SymbolicUtils: SymbolicUtils, simplify, substitute, add_worker using Symbolics: Symbolics, Num, expand, @variables, build_function, symbolic_to_float using TermInterface: TermInterface -# `CNum` (the coefficient type `Coeff`) is defined in `expressions/cnum.jl`. - -using QuantumOpticsBase: QuantumOpticsBase -import QuantumOpticsBase: ⊗, tensor, expect +import QuantumInterface: ⊗, tensor, expect, basis +using QuantumInterface: AbstractOperator, StateVector, Basis using Combinatorics: with_replacement_combinations -using FunctionWrappers: FunctionWrapper using Latexify: Latexify, latexify, @latexrecipe using PrecompileTools: @setup_workload, @compile_workload using SciMLPublic: @public @@ -43,7 +40,12 @@ include("algebra/mutable_arithmetics.jl") include("algebra/weyl.jl") include("average.jl") -include("numeric.jl") + +include("numeric/backend.jl") +include("numeric/coeff.jl") +include("numeric/core.jl") +include("numeric/indexed.jl") +include("numeric/api.jl") include("printing/printing.jl") include("printing/latexify_recipes.jl") @@ -108,6 +110,11 @@ export FockSpace, ProductSpace, normal_order, normal_to_symmetric, symmetric_to_normal, simplify, expand, expand_completeness, assume_distinct_index, commutator, anticommutator, to_numeric, numeric_average, + NumericBackend, QuantumOpticsBackend, QuantumToolboxBackend, + numeric_operator, numeric_basis, numeric_subbasis, numeric_embed, + numeric_identity, numeric_num_subsystems, + numeric_assemble, numeric_assemble_td, numeric_materialize, numeric_expect, + numeric_backend, qadjoint, qconj, dagger, inner_adjoint, Op, operator_name, operator_index, is_destroy, is_create, is_transition, is_collective_transition, is_pauli, is_spin, is_position, is_momentum, optype @@ -116,6 +123,7 @@ export FockSpace, ProductSpace, # Public API that is intentionally NOT exported — accessed as # `SecondQuantizedAlgebra.symbol`. @public HilbertSpace, QField, QSym, OpKind, + NumericContext, expect, OP_DESTROY, OP_CREATE, OP_TRANSITION, OP_PAULI, OP_SPIN, OP_POSITION, OP_MOMENTUM, OP_COLLECTIVE_TRANSITION, QAdd, QTerm, QTermDict, Coeff, CNum, diff --git a/src/numeric.jl b/src/numeric.jl deleted file mode 100644 index cf30235c..00000000 --- a/src/numeric.jl +++ /dev/null @@ -1,787 +0,0 @@ -const QuantumState = Union{QuantumOpticsBase.StateVector, QuantumOpticsBase.AbstractOperator} - -# `Union{}` value keeps the haskey-fallback branch concrete-typed. -const _NO_SUBS = Dict{QSym, Union{}}() - -""" - to_numeric(op, basis [, d::AbstractDict{<:QSym}]) - to_numeric(op, state [, d::AbstractDict{<:QSym}]) - to_numeric(op, basis; parameter=Dict(), time_parameter=Dict(), operators=Dict(), - adjoint_ops=true, op_type=sparse) - to_numeric(ops::AbstractVector, basis; kwargs...) - -Convert a symbolic operator expression to a numeric QuantumOpticsBase operator. -`d` substitutes individual `QSym`s with custom numeric operators. Throws -`ArgumentError` if any prefactor cannot be reduced to a concrete `ComplexF64`. - -The keyword form first substitutes scalar `parameter`s, then translates using -custom numeric `operators`. Missing adjoint entries are added to `operators` -when `adjoint_ops=true`. If `time_parameter` is non-empty, values may be numbers -or functions of time and the result is a closure `t -> op(t)`. - -The result representation is controlled by `op_type` and does not depend on the -shape of the expression: a bare operator, a product, and a sum all return the -same operator type. `op_type` is applied once to the assembled operator: - - - `op_type=sparse` (default): a sparse `Operator`. A good general default: it - is memory-scalable, type-stable, and the form the QuantumOptics solvers - consume most efficiently. - - `op_type=dense`: a dense `Operator`. Prefer for small Hilbert spaces. - - `op_type=identity`: the natural lazy representation (`LazyTensor` / - `LazyProduct` / `LazySum`). Opt in for large tensor-product spaces where an - operator is local to a few subsystems, so the full Kronecker product is - never materialised. (In the time-dependent path with more than one term, - terms are still materialised to keep the returned closure type-stable.) - -The positional forms `to_numeric(op, basis[, d])` always return `sparse`; use -the keyword form to select `dense` or `identity`. - -# Examples - -```jldoctest -julia> using QuantumOpticsBase - -julia> h = FockSpace(:f); - -julia> @qnumbers a::Destroy(h); - -julia> b = FockBasis(5); ψ = fockstate(b, 2); - -julia> real(QuantumOpticsBase.expect(to_numeric(a' * a, ψ), ψ)) ≈ 2 -true -``` - -See also [`numeric_average`](@ref). -""" -function to_numeric end - -function to_numeric(op::QField, b::QuantumOpticsBase.Basis; kwargs...) - return _to_numeric_kw(op, b; kwargs...) -end - -function to_numeric(x::Union{Number, SymbolicUtils.BasicSymbolic}, b::QuantumOpticsBase.Basis; kwargs...) - return _to_numeric_kw(x, b; kwargs...) -end - -function to_numeric(op, state::QuantumState; kwargs...) - return to_numeric(op, QuantumOpticsBase.basis(state); kwargs...) -end - -function to_numeric(ops::AbstractVector, b::QuantumOpticsBase.Basis; kwargs...) - return [to_numeric(op, b; kwargs...) for op in ops] -end - -# Per-basis kind-chain on the single concrete `Op`. -function to_numeric(op::Op, b::QuantumOpticsBase.FockBasis) - is_destroy(op) && return QuantumOpticsBase.destroy(b) - is_create(op) && return QuantumOpticsBase.create(b) - is_position(op) && return (QuantumOpticsBase.destroy(b) + QuantumOpticsBase.create(b)) / sqrt(2) - is_momentum(op) && return im * (QuantumOpticsBase.create(b) - QuantumOpticsBase.destroy(b)) / sqrt(2) - throw(ArgumentError("Op kind $(op.kind) does not act on a FockBasis")) -end -function to_numeric(op::Op, b::QuantumOpticsBase.NLevelBasis) - is_transition(op) && return QuantumOpticsBase.transition(b, Int(op.l1), Int(op.l2)) - throw(ArgumentError("Op kind $(op.kind) does not act on an NLevelBasis")) -end -function to_numeric(op::Op, b::QuantumOpticsBase.ManyBodyBasis) - is_collective_transition(op) || throw(ArgumentError("Op kind $(op.kind) does not act on a ManyBodyBasis")) - onebody = b.onebodybasis - onebody isa QuantumOpticsBase.NLevelBasis || - throw( - ArgumentError( - "CollectiveTransition requires a ManyBodyBasis with NLevelBasis one-body basis; got $(typeof(onebody))", - ), - ) - onebody_op = QuantumOpticsBase.transition(onebody, Int(op.l1), Int(op.l2)) - return QuantumOpticsBase.manybodyoperator(b, onebody_op) -end -function to_numeric(op::Op, b::QuantumOpticsBase.SpinBasis) - if is_pauli(op) - b.spinnumber == 1 // 2 || throw(ArgumentError("Pauli operators require SpinBasis(1//2), got SpinBasis($(b.spinnumber))")) - op.l1 == 1 && return QuantumOpticsBase.sigmax(b) - op.l1 == 2 && return QuantumOpticsBase.sigmay(b) - return QuantumOpticsBase.sigmaz(b) - end - if is_spin(op) - op.l1 == 1 && return 0.5 * QuantumOpticsBase.sigmax(b) - op.l1 == 2 && return 0.5 * QuantumOpticsBase.sigmay(b) - return 0.5 * QuantumOpticsBase.sigmaz(b) - end - throw(ArgumentError("Op kind $(op.kind) does not act on a SpinBasis")) -end -# Lazy embedding of a single operator on `b`; building block for products/sums. -_embed(op::Op, b::QuantumOpticsBase.Basis) = to_numeric(op, b) -function _embed(op::Op, b::QuantumOpticsBase.CompositeBasis) - si = Int(op.space_index) - return QuantumOpticsBase.LazyTensor(b, [si], (_embed(op, b.bases[si]),)) -end -function _embed(op::Op, b::QuantumOpticsBase.Basis, d::AbstractDict{<:QSym}) - haskey(d, op) && return d[op] - return _embed(op, b) -end - -# Positional forms materialise sparse; use the keyword form for other `op_type`. -to_numeric(op::Op, b::QuantumOpticsBase.CompositeBasis) = - QuantumOpticsBase.sparse(_embed(op, b)) -to_numeric(op::Op, b::QuantumOpticsBase.Basis, d::AbstractDict{<:QSym}) = - QuantumOpticsBase.sparse(_embed(op, b, d)) -to_numeric(s::QAdd, b::QuantumOpticsBase.Basis, d::AbstractDict{<:QSym} = _NO_SUBS) = - QuantumOpticsBase.sparse(_assemble(s, b, d)) - -# Lazy assembly of an operator sum; the caller materialises it. -function _assemble(s::QAdd, b::QuantumOpticsBase.Basis, d::AbstractDict{<:QSym} = _NO_SUBS) - iter = iterate(s.arguments) - iter === nothing && return _to_complex(_CNUM_ZERO) * _lazy_one(b) - (term, c), st = iter - result = isempty(term.ops) ? _to_complex(c) * _lazy_one(b) : _to_complex(c) * _product(term.ops, b, d) - while true - next = iterate(s.arguments, st) - next === nothing && return result - (term, c), st = next - result += isempty(term.ops) ? _to_complex(c) * _lazy_one(b) : _to_complex(c) * _product(term.ops, b, d) - end - return -end - -function _product(ops::Vector{Op}, b::QuantumOpticsBase.Basis, d::AbstractDict{<:QSym} = _NO_SUBS) - acc = _embed(first(ops), b, d) - for i in 2:length(ops) - acc *= _embed(ops[i], b, d) - end - return acc -end - -# Lazy form for consumers that contract immediately (e.g. `expect` on a `LazyKet`). -_to_numeric_lazy(op::Op, b::QuantumOpticsBase.Basis, d::AbstractDict{<:QSym}) = _embed(op, b, d) -_to_numeric_lazy(s::QAdd, b::QuantumOpticsBase.Basis, d::AbstractDict{<:QSym}) = _assemble(s, b, d) - -function _to_numeric_kw( - op, - b::QuantumOpticsBase.Basis; - parameter = Dict(), - time_parameter = Dict(), - operators = Dict{QSym, Any}(), - adjoint_ops = true, - op_type = QuantumOpticsBase.sparse, - ) - param = _expand_parameter(parameter) - tp = _normalize_time_parameter(time_parameter) - ops = _numeric_operator_dict(operators, adjoint_ops) - return _to_numeric_translated( - op, b; parameter = param, time_parameter = tp, operators = ops, op_type, - ) -end - -function _to_numeric_translated( - op::QSym, - b::QuantumOpticsBase.Basis; - parameter, - time_parameter, - operators, - op_type, - ) - return _to_numeric_translated( - _single_qadd(_CNUM_ONE, Op[op]), b; - parameter, time_parameter, operators, op_type, - ) -end - -function _to_numeric_translated( - op::QAdd, - b::QuantumOpticsBase.Basis; - parameter, - time_parameter, - operators, - op_type, - ) - if isempty(time_parameter) - return _to_numeric_static(op, b; parameter, operators, op_type) - end - - op_ = substitute(op, parameter) - if iszero(op_) - z = op_type(_to_complex(_CNUM_ZERO) * _lazy_one(b)) - return t -> z - end - - pairs = collect(op_) - if length(pairs) == 1 - term, c = pairs[1] - return _translate_term( - term.ops, to_num(c), b, time_parameter, operators, op_type, - ) - end - - # Use a concrete operator representation for every term so the tuple of - # wrapped closures has one return type, even for scalar-plus-operator sums. - term_op_type = op_type === identity ? QuantumOpticsBase.sparse : op_type - first_res = _translate_term( - pairs[1].first.ops, to_num(pairs[1].second), b, - time_parameter, operators, term_op_type, - ) - OpType = typeof(first_res(0.0)) - FW = FunctionWrapper{OpType, Tuple{Float64}} - wrapped = ntuple(length(pairs)) do k - res = k == 1 ? first_res : _translate_term( - pairs[k].first.ops, to_num(pairs[k].second), b, - time_parameter, operators, term_op_type, - ) - FW(res) - end - - return t -> begin - tt = Float64(t) - result = wrapped[1](tt) - for i in 2:length(wrapped) - result = result + wrapped[i](tt) - end - result - end -end - -function _to_numeric_translated( - arg, - b::QuantumOpticsBase.Basis; - parameter, - time_parameter, - operators, - op_type, - ) - arg_sub = substitute(arg, parameter) - one_b = _lazy_one(b) - c = _as_cnum(arg_sub) - - if isempty(time_parameter) - _coeff_is_const(c) || throw( - ArgumentError( - "cannot translate symbolic scalar `$arg_sub` without a value: supply it via " * - "`parameter` or `time_parameter`.", - ), - ) - return op_type(_const_coeff(c) * one_b) - end - - if _coeff_is_const(c) - val = _const_coeff(c) - op = op_type(val * one_b) - return t -> op - end - basevars, valuefuncs = _time_basis(time_parameter) - _check_time_variables(c, basevars) - pref = _compile_coeff(c, basevars...) - return t -> op_type(pref(map(f -> f(t), valuefuncs)...) * one_b) -end - -function _to_numeric_static( - op::QAdd, - b::QuantumOpticsBase.Basis; - parameter, - operators, - op_type, - ) - op_ = substitute(op, parameter) - iszero(op_) && return op_type(_to_complex(_CNUM_ZERO) * _lazy_one(b)) - - result = nothing - for (term, c_) in op_ - c = to_num(c_) - _coeff_is_const(c) || throw( - ArgumentError( - "cannot translate `$op` to a static operator: coefficient `$c` still " * - "depends on a free variable; supply it via `parameter` or `time_parameter`.", - ), - ) - contrib = _const_coeff(c) * _numeric_product(term.ops, b, operators, identity) - result = result === nothing ? contrib : result + contrib - end - return op_type(result) -end - -_to_numeric_with_ops(op, b::QuantumOpticsBase.Basis, operators::AbstractDict{<:QSym}) = - isempty(operators) ? _embed(op, b) : _embed(op, b, operators) - -function _numeric_product( - ops::Vector{Op}, - b::QuantumOpticsBase.Basis, - operators::AbstractDict{<:QSym}, - op_type, - ) - if isempty(ops) - return op_type(_lazy_one(b)) - end - acc = _to_numeric_with_ops(ops[1], b, operators) - for i in 2:length(ops) - acc *= _to_numeric_with_ops(ops[i], b, operators) - end - return op_type(acc) -end - -function _numeric_operator_dict(operators, adjoint_ops::Bool) - out = Dict{QSym, Any}() - for (k, v) in operators - k isa QSym || throw(ArgumentError("operator substitution key `$k` is not a QSym")) - out[k] = v - end - if adjoint_ops - for (k, v) in operators - k isa QSym || continue - k_adj = Base.adjoint(k) - haskey(out, k_adj) || (out[k_adj] = Base.adjoint(v)) - end - end - return out -end - -function _expand_parameter(parameter) - isempty(parameter) && return parameter - out = Dict{Any, Any}() - for (k, v) in parameter - if k isa Complex - out[real(k)] = real(v) - ik = imag(k) - SymbolicUtils.unwrap(ik) isa Number || (out[ik] = imag(v)) - else - out[k] = v - end - end - return out -end - -function _normalize_time_parameter(time_parameter) - isempty(time_parameter) && return time_parameter - out = Dict{Any, Any}() - for (k, v) in time_parameter - out[k] = v isa Number ? (t -> v) : v - end - return out -end - -function _time_basis(time_parameter) - basevars = Any[] - valuefuncs = Any[] - for (k, f) in time_parameter - uk = SymbolicUtils.unwrap(k) - if SymbolicUtils.issym(uk) - push!(basevars, k) - push!(valuefuncs, f) - continue - end - vs = collect(Symbolics.get_variables(k)) - length(vs) == 1 || throw( - ArgumentError( - "time_parameter key `$k` must depend on exactly one variable, got $(length(vs)).", - ), - ) - wv = Symbolics.wrap(vs[1]) - if isequal(uk, SymbolicUtils.unwrap(conj(wv))) - push!(basevars, wv) - push!(valuefuncs, t -> conj(f(t))) - else - throw( - ArgumentError( - "unsupported time_parameter key `$k`; only a bare variable `v` or `conj(v)` are supported.", - ), - ) - end - end - return basevars, Tuple(valuefuncs) -end - -_coeff_is_const(c::Complex{Num}) = - isempty(Symbolics.get_variables(real(c))) && isempty(Symbolics.get_variables(imag(c))) - -_const_coeff(c::Complex{Num}) = _to_complex(c) - -_as_cnum(x::Complex) = Complex{Num}(Num(real(x)), Num(imag(x))) -_as_cnum(x) = Complex{Num}(Num(x), Num(false)) - -function _coefficient_variables(c::Complex{Num}) - vars = Any[] - append!(vars, Symbolics.get_variables(real(c))) - append!(vars, Symbolics.get_variables(imag(c))) - unique!(vars) - return vars -end - -function _check_time_variables(c::Complex{Num}, basevars) - vars = _coefficient_variables(c) - base_unwrapped = SymbolicUtils.unwrap.(basevars) - missing = Any[] - for v in vars - any(b -> isequal(v, b), base_unwrapped) || push!(missing, v) - end - isempty(missing) && return nothing - throw( - ArgumentError( - "time-dependent coefficient `$c` depends on variables without time values: " * - join(string.(missing), ", "), - ), - ) -end - -function _compile_coeff(c::Complex{Num}, vars...) - f_re = build_function(real(c), vars...; expression = Val(false)) - f_im = build_function(imag(c), vars...; expression = Val(false)) - g_re = f_re isa Tuple ? first(f_re) : f_re - g_im = f_im isa Tuple ? first(f_im) : f_im - return (vals...) -> g_re(vals...) + im * g_im(vals...) -end - -function _translate_term( - ops::Vector{Op}, - c::Complex{Num}, - b::QuantumOpticsBase.Basis, - time_parameter, - operators, - op_type, - ) - prodop = _numeric_product(ops, b, operators, op_type) - if _coeff_is_const(c) - op = _const_coeff(c) * prodop - return t -> op - end - basevars, valuefuncs = _time_basis(time_parameter) - _check_time_variables(c, basevars) - pref = _compile_coeff(c, basevars...) - return t -> pref(map(f -> f(t), valuefuncs)...) * prodop -end - -const _NO_SCALAR_SUBS = Dict{Num, Any}() - -""" - to_numeric(q::QAdd, b::CompositeBasis, sites::Dict{Int, Vector{Int}}[, d[, scalar_subs]]) - -Convert a `QAdd` that carries bound summation indices (`q.indices`) to a numeric -operator on a `CompositeBasis` whose layout replicates one or more abstract -subspaces. `sites[orig_space_index]` is the list of slots into `b.bases` that -realize that subspace; non-indexed subspaces map to a single-slot entry. - -For each term, every bound index that the term depends on is unrolled over the -length of its `sites` vector, with `term.ne` constraints filtering out diagonal -combinations. Concrete-site operators are routed to `sites[op.space_index][k]`, -where `k` is the substituted integer site. Terms that do not depend on any -bound index are emitted once as written, matching the `Σ` convention that -i-independent residuals already carry their range factor. - -`d` substitutes individual `QSym` operators with custom numeric operators (same -semantics as the single-basis-slot path). `scalar_subs` substitutes symbolic -scalar parameters (e.g. `Δ` or `g(k)` for the integer `k` produced by index -unrolling) with numeric values; supply this to resolve indexed parameters that -only become fully concrete after the symbolic sum is unrolled. -""" -function to_numeric( - q::QAdd, b::QuantumOpticsBase.CompositeBasis, - sites::AbstractDict{Int, Vector{Int}}, - d::AbstractDict{<:QSym} = _NO_SUBS, - scalar_subs::AbstractDict = _NO_SCALAR_SUBS, - ) - if isempty(q.indices) - return to_numeric(q, b, d) - end - sub_re, sub_im, has_imag = _split_scalar_subs(scalar_subs) - result = _to_complex(_CNUM_ZERO) * _lazy_one(b) - for (term, c) in q.arguments - result = _accumulate_indexed_term!( - result, term, c, q.indices, b, sites, d, sub_re, sub_im, has_imag, - ) - end - return QuantumOpticsBase.sparse(result) -end - -function _accumulate_indexed_term!( - acc, term::QTerm, c::CNum, indices::Vector{Index}, - b::QuantumOpticsBase.CompositeBasis, - sites::AbstractDict{Int, Vector{Int}}, - d::AbstractDict{<:QSym}, - sub_re::Dict, sub_im::Dict, has_imag::Bool, - ) - dep_indices = Index[idx for idx in indices if _depends_on_index_term(c, term.ops, idx)] - if isempty(dep_indices) - c_resolved = _apply_scalar_subs(c, sub_re, sub_im, has_imag) - return acc + _emit_indexed_combo(term.ops, c_resolved, b, sites, d) - end - lens = Int[length(sites[Int(idx.space_index)]) for idx in dep_indices] - total = prod(lens) - sub_op = Dict{Index, Index}() - sub_coef = Dict{Index, Index}() - for combo in 1:total - empty!(sub_op) - empty!(sub_coef) - rem = combo - 1 - for k in 1:length(dep_indices) - kpos = (rem % lens[k]) + 1 - rem ÷= lens[k] - idx = dep_indices[k] - # Operators keep the index name (slot kpos drives routing / ne checks and - # lets a resolved op still match a user `d` key); coefficients use the - # anonymous name_id-0 form so `index_sym` is Num(kpos), resolving g(i)→g(k). - sub_op[idx] = Index(idx.name_id, idx.range_id, idx.space_index, Int32(kpos)) - sub_coef[idx] = Index(Int32(0), idx.range_id, idx.space_index, Int32(kpos)) - end - if _violates_ne(term.ne, sub_op) - continue - end - new_ops = Op[change_index(op, sub_op) for op in term.ops] - new_c = change_index(c, sub_coef) - new_c = _apply_scalar_subs(new_c, sub_re, sub_im, has_imag) - acc = acc + _emit_indexed_combo(new_ops, new_c, b, sites, d) - end - return acc -end - -# Split user-supplied scalar substitutions into real and imag legs once per -# `to_numeric` call. A complex RHS `re + im*ii` propagates as separate real -# and imaginary substitutions, preserving the `Complex{Num}` invariant. -function _split_scalar_subs(scalar_subs::AbstractDict) - sub_re = Dict{Any, Any}() - sub_im = Dict{Any, Any}() - has_imag = false - for (k, v) in scalar_subs - kraw = SymbolicUtils.unwrap(k) - if v isa Complex - sub_re[kraw] = real(v) - sub_im[kraw] = imag(v) - iszero(imag(v)) || (has_imag = true) - else - sub_re[kraw] = v - sub_im[kraw] = 0 - end - end - return sub_re, sub_im, has_imag -end - -function _apply_scalar_subs(c::CNum, sub_re::Dict, sub_im::Dict, has_imag::Bool) - (isempty(sub_re) || _is_native(c)) && return c - re_part = SymbolicUtils.unwrap(real(c)) - im_part = SymbolicUtils.unwrap(imag(c)) - if !has_imag - return _cnum( - Num(Symbolics.substitute(re_part, sub_re)), - Num(Symbolics.substitute(im_part, sub_re)), - ) - end - # (re + i*im) -> (re|sub_re - im|sub_im) + i*(re|sub_im + im|sub_re) - new_re = Symbolics.substitute(re_part, sub_re) - - Symbolics.substitute(im_part, sub_im) - new_im = Symbolics.substitute(re_part, sub_im) + - Symbolics.substitute(im_part, sub_re) - return _cnum(Num(new_re), Num(new_im)) -end - -function _violates_ne(ne::Vector{NonEqualPair}, sub::Dict{Index, Index}) - for (a, b) in ne - ra = get(sub, a, a) - rb = get(sub, b, b) - va = index_slot(ra) - vb = index_slot(rb) - va === nothing && continue - vb === nothing && continue - ra.space_index == rb.space_index || continue - va == vb && return true - end - return false -end - -function _emit_indexed_combo( - ops::Vector{Op}, c::CNum, - b::QuantumOpticsBase.CompositeBasis, - sites::AbstractDict{Int, Vector{Int}}, - d::AbstractDict{<:QSym}, - ) - if isempty(ops) - return _to_complex(c) * _lazy_one(b) - end - acc = _site_routed_op(ops[1], b, sites, d) - for k in 2:length(ops) - acc = acc * _site_routed_op(ops[k], b, sites, d) - end - return _to_complex(c) * acc -end - -function _site_routed_op( - op::QSym, b::QuantumOpticsBase.CompositeBasis, - sites::AbstractDict{Int, Vector{Int}}, - d::AbstractDict{<:QSym}, - ) - slot = _resolve_slot(op, sites) - haskey(d, op) && return d[op] - return QuantumOpticsBase.LazyTensor(b, [slot], (to_numeric(op, b.bases[slot]),)) -end - -function _resolve_slot(op::QSym, sites::AbstractDict{Int, Vector{Int}}) - si = Int(op.space_index) - slots = get(sites, si, Int[]) - if isempty(slots) - return si - end - if has_index(op.index) - v = index_slot(op.index) - if v !== nothing && 1 <= v <= length(slots) - return slots[v] - end - end - length(slots) == 1 && return slots[1] - throw( - ArgumentError( - "cannot resolve slot for operator $(op): index is not a concrete integer site, " * - "and sites[$si] has $(length(slots)) candidates", - ) - ) -end - -to_numeric(x::Number, b::QuantumOpticsBase.Basis, ::AbstractDict{<:QSym} = _NO_SUBS) = QuantumOpticsBase.sparse(_to_complex(x) * _lazy_one(b)) -to_numeric(op::QField, state::QuantumState, d::AbstractDict{<:QSym} = _NO_SUBS) = to_numeric(op, QuantumOpticsBase.basis(state), d) -to_numeric(x::Number, state::QuantumState) = to_numeric(x, QuantumOpticsBase.basis(state)) - -_lazy_one(b::QuantumOpticsBase.Basis) = one(b) -_lazy_one(b::QuantumOpticsBase.CompositeBasis) = QuantumOpticsBase.LazyTensor(b, collect(1:length(b.bases)), Tuple(one(bi) for bi in b.bases)) - -function _reduce_const(n::Num)::ComplexF64 - x = Symbolics.value(n) - try - return _fold_const(x) - catch err - err isa ArgumentError || rethrow() - isempty(Symbolics.get_variables(n)) || rethrow() - return _compile_const(n) - end -end - -_compile_const(n::Num)::ComplexF64 = ComplexF64(symbolic_to_float(n)) - -function _fold_const(x)::ComplexF64 - x isa Number && return x - if x isa SymbolicUtils.BasicSymbolic - if SymbolicUtils.iscall(x) - op = SymbolicUtils.operation(x) - args = SymbolicUtils.arguments(x) - if op === (+) - acc = zero(ComplexF64) - for a in args - acc += _fold_const(a) - end - return acc - elseif op === (*) - acc = one(ComplexF64) - for a in args - acc *= _fold_const(a) - end - return acc - elseif op === conj - return conj(_fold_const(first(args))) - elseif op === real - return real(_fold_const(first(args))) - elseif op === imag - return imag(_fold_const(first(args))) - elseif op === (/) - return _fold_const(first(args)) / _fold_const(last(args)) - elseif op === (^) - return _fold_const(first(args))^_fold_const(last(args)) - elseif op === (-) - return length(args) == 1 ? -_fold_const(first(args)) : - _fold_const(first(args)) - _fold_const(last(args)) - end - elseif SymbolicUtils.isconst(x) - return x.val - end - end - throw(ArgumentError("cannot reduce symbolic expression $x to a concrete number")) -end - -_to_complex(c::Coeff) = _is_native(c) ? c.z : _to_complex(to_num(c)) - -# One method (union-split budget) routing every input through `convert ∘ Complex`, -# the only pattern that infers to ComplexF64 from `Any` after `Symbolics.value`. -function _to_complex(x) - x isa ComplexF64 && return x - x isa Complex{Num} && return _reduce_const(real(x)) + im * _reduce_const(imag(x)) - x isa Complex && return convert(ComplexF64, x) - x isa Num && return _reduce_const(x) - if x isa SymbolicUtils.BasicSymbolic - SymbolicUtils.isconst(x) || throw(ArgumentError("cannot reduce symbolic expression $x to a concrete number")) - return convert(ComplexF64, Complex(x.val, false)) - end - return convert(ComplexF64, Complex(x, false)) -end - -""" - numeric_average(op, state[, d::AbstractDict{<:QSym}]) -> ComplexF64 - numeric_average(op, states::AbstractVector[, d::AbstractDict{<:QSym}]) -> Vector - -Expectation value ``\\langle \\psi | \\hat{O} | \\psi \\rangle`` of a symbolic -operator expression. Averaged `BasicSymbolic` expressions (from -[`average`](@ref)) are unwrapped automatically. `d` substitutes symbolic -operators with custom numeric operators before evaluation. - -# Examples - -```jldoctest -julia> using QuantumOpticsBase - -julia> h = FockSpace(:f); - -julia> @qnumbers a::Destroy(h); - -julia> b = FockBasis(5); ψ = fockstate(b, 2); - -julia> real(numeric_average(a' * a, ψ)) ≈ 2 -true -``` - -See also [`to_numeric`](@ref), [`average`](@ref). -""" -function numeric_average end - -numeric_average(op, state::QuantumState, d::AbstractDict{<:QSym} = _NO_SUBS) = _numeric_average(op, state, d) -function numeric_average(op, states::AbstractVector, d::AbstractDict{<:QSym} = _NO_SUBS) - isempty(states) && throw(ArgumentError("numeric_average: states vector is empty")) - return _numeric_average_vec(op, states, d) -end - -_numeric_average(op::QField, state::QuantumState, d::AbstractDict{<:QSym}) = ComplexF64(QuantumOpticsBase.expect(_to_numeric_lazy(op, QuantumOpticsBase.basis(state), d), state)) -_numeric_average(x::Number, ::QuantumState, ::AbstractDict{<:QSym}) = _to_complex(x) -_numeric_average(x::Num, state::QuantumState, d::AbstractDict{<:QSym}) = _numeric_average(SymbolicUtils.unwrap(x), state, d) - -# `arguments(::BasicSymbolic)` is statically `Any`; the `_to_complex` wraps -# reseal each recursive result to `ComplexF64`. -function _numeric_average(avg::SymbolicUtils.BasicSymbolic, state::QuantumState, d::AbstractDict{<:QSym}) - if SymbolicUtils.hasmetadata(avg, AverageOperator) || - (SymbolicUtils.iscall(avg) && SymbolicUtils.operation(avg) isa AvgFunc) || - is_indexed_sum(avg) - return _to_complex(_numeric_average(undo_average(avg), state, d)) - end - SymbolicUtils.isconst(avg) && return _to_complex(_numeric_average(avg.val, state, d)) - SymbolicUtils.iscall(avg) || throw(ArgumentError("numeric_average not implemented for $(typeof(avg))")) - f, args = SymbolicUtils.operation(avg), SymbolicUtils.arguments(avg) - f === (^) && return _to_complex(_numeric_average(args[1], state, d))^_to_complex(args[2]) - f === (+) && return _sum_avg(args, state, d) - f === (*) && return _prod_avg(args, state, d) - throw(ArgumentError("numeric_average not implemented for operation $f")) -end - -function _sum_avg(args, state::QuantumState, d::AbstractDict{<:QSym}) - r = ComplexF64(0) - for a in args - r += _to_complex(_numeric_average(a, state, d)) - end - return r -end -function _prod_avg(args, state::QuantumState, d::AbstractDict{<:QSym}) - r = ComplexF64(1) - for a in args - r *= _to_complex(_numeric_average(a, state, d)) - end - return r -end - -_numeric_average_vec(op::QField, states::AbstractVector, d::AbstractDict{<:QSym}) = QuantumOpticsBase.expect(QuantumOpticsBase.sparse(to_numeric(op, first(states), d)), states) -_numeric_average_vec(op, states::AbstractVector, d::AbstractDict{<:QSym}) = [_numeric_average(op, ψ, d) for ψ in states] - -""" - expect(op::QField, state[, d::AbstractDict{<:QSym}]) - -Alias for [`numeric_average`](@ref) on operator expressions. For symbolic -scalar expressions such as `average(op)`, call [`numeric_average`](@ref) -directly. -""" -expect(op::QField, state, d::AbstractDict{<:QSym} = _NO_SUBS) = numeric_average(op, state, d) diff --git a/src/numeric/api.jl b/src/numeric/api.jl new file mode 100644 index 00000000..adaaa25e --- /dev/null +++ b/src/numeric/api.jl @@ -0,0 +1,216 @@ +# Public numeric API. Signatures use only QuantumInterface types (`Basis`, `StateVector`, +# `AbstractOperator`) plus the backend singletons, so `src/` never loads a numeric backend. +# The QuantumToolbox extension adds the parallel `QuantumObject` state methods. + +const QuantumState = Union{StateVector, AbstractOperator} + +""" + to_numeric(op, h::HilbertSpace, dims; backend, kwargs...) + to_numeric(op, basis [, d]; kwargs...) + to_numeric(op, state [, d]; kwargs...) + to_numeric(ops::AbstractVector, ...; kwargs...) + +Convert a symbolic operator expression to a numeric operator on the chosen backend. + +The first form is backend-neutral and is recommended for portable code. Load +`QuantumOpticsBase` or `QuantumToolbox` and pass `QuantumOpticsBackend()` or +`QuantumToolboxBackend()`; `backend` may be omitted when exactly one bundled backend is +loaded. Backend-native basis, state, and dimension forms are convenience overloads. + +The representation contract is independent of the expression shape: + +| Input | `op_type` | Result | +|:--|:--|:--| +| static | `nothing` (default) | ordinary eager backend operator | +| static | `identity` | backend's lazy assembly | +| static | another callable | backend-defined representation | +| time-dependent | `nothing` or `identity` | native time-dependent operator | + +Both bundled backends use sparse storage for the default static result. Examples of explicit +conversions are `dense` for QuantumOptics and `QuantumToolbox.to_dense` for QuantumToolbox. +Other backends may choose a different eager representation. Time-dependent conversion +cannot apply a one-time materializer, so other `op_type` values are rejected. + +Keyword arguments: + +- `parameter=Dict()`: substitutions for scalar symbolic coefficients. +- `time_parameter=Dict()`: scalar values or functions of time. A non-empty mapping selects + time-dependent conversion (`TimeDependentSum` for QuantumOptics, `QobjEvo` for + QuantumToolbox). +- `operators=Dict()`: replacements for individual symbolic operators. Missing adjoints are + inferred when `adjoint_ops=true` (the default). Positional `d` supplies exact replacements + without this inference. +- `op_type=nothing`: representation selection as described above. + +`dims` follows the symbolic space convention: a Fock value is a cutoff, so `5` creates a +six-dimensional space. In contrast, QuantumToolbox's native `to_numeric(op, dims)` +convenience form accepts raw matrix dimensions. Conversion throws `ArgumentError` when a +coefficient cannot be reduced to `ComplexF64`. + +# Examples + +```jldoctest +julia> using QuantumOpticsBase + +julia> h = FockSpace(:f); + +julia> @qnumbers a::Destroy(h); + +julia> b = FockBasis(5); ψ = fockstate(b, 2); + +julia> real(expect(to_numeric(a' * a, b), ψ)) ≈ 2 +true +``` + +See also [`numeric_average`](@ref). +""" +function to_numeric end + +# --- leaf / QAdd / scalar on an explicit basis (positional, no kwargs) ----------------- +# The positional forms always materialize `sparse` (`op_type` is only on the keyword form). + +to_numeric(op::Op, b::Basis, d::AbstractDict{<:QSym} = _NO_SUBS) = + numeric_materialize(QuantumOpticsBackend(), _numeric_leaf(op, NumericContext(QuantumOpticsBackend(), b, d)), nothing) + +to_numeric(q::QAdd, b::Basis, d::AbstractDict{<:QSym} = _NO_SUBS) = + numeric_materialize(QuantumOpticsBackend(), _to_numeric_static(q, NumericContext(QuantumOpticsBackend(), b, d)), nothing) + +to_numeric(x::Number, b::Basis, ::AbstractDict{<:QSym} = _NO_SUBS) = + numeric_materialize(QuantumOpticsBackend(), _to_complex(x) * numeric_identity(QuantumOpticsBackend(), b), nothing) + +# --- keyword form on an explicit basis ------------------------------------------------- + +function to_numeric( + op::QField, b::Basis; + parameter = Dict(), time_parameter = Dict(), + operators = Dict{QSym, Any}(), adjoint_ops = true, op_type = nothing, + ) + return _to_numeric_kw(QuantumOpticsBackend(), op, b; parameter, time_parameter, operators, adjoint_ops, op_type) +end + +function to_numeric( + x::Union{Number, SymbolicUtils.BasicSymbolic}, b::Basis; + parameter = Dict(), time_parameter = Dict(), + operators = Dict{QSym, Any}(), adjoint_ops = true, op_type = nothing, + ) + return _to_numeric_kw(QuantumOpticsBackend(), x, b; parameter, time_parameter, operators, adjoint_ops, op_type) +end + +# --- uniform HilbertSpace form (the only one that works for both backends) ------------- + +function to_numeric( + op, h::HilbertSpace, dims; + backend::NumericBackend = _default_backend(), + parameter = Dict(), time_parameter = Dict(), + operators = Dict{QSym, Any}(), adjoint_ops = true, op_type = nothing, + ) + b = numeric_basis(backend, h, dims) + return _to_numeric_kw(backend, op, b; parameter, time_parameter, operators, adjoint_ops, op_type) +end + +# Vector of operators on the uniform HilbertSpace form (both backends). +to_numeric(ops::AbstractVector, h::HilbertSpace, dims; backend::NumericBackend = _default_backend(), kwargs...) = + [to_numeric(op, h, dims; backend, kwargs...) for op in ops] + +# --- state convenience (basis derived from the state) ---------------------------------- + +to_numeric(op, state::QuantumState; kwargs...) = to_numeric(op, numeric_basis(state); kwargs...) +to_numeric(op::QField, state::QuantumState, d::AbstractDict{<:QSym} = _NO_SUBS) = + to_numeric(op, numeric_basis(state), d) +to_numeric(x::Number, state::QuantumState) = to_numeric(x, numeric_basis(state)) + +# --- vector of operators --------------------------------------------------------------- + +to_numeric(ops::AbstractVector, b::Basis; kwargs...) = [to_numeric(op, b; kwargs...) for op in ops] + +# ======================================================================================= + +""" + numeric_average(op, state[, d::AbstractDict{<:QSym}]) -> ComplexF64 + numeric_average(op, states::AbstractVector[, d::AbstractDict{<:QSym}]) -> Vector + +Expectation value ``\\langle \\psi | \\hat{O} | \\psi \\rangle`` of a symbolic operator +expression. Averaged `BasicSymbolic` expressions (from [`average`](@ref)) are unwrapped +automatically. `d` substitutes symbolic operators with custom numeric operators before +evaluation. + +# Examples + +```jldoctest +julia> using QuantumOpticsBase + +julia> h = FockSpace(:f); + +julia> @qnumbers a::Destroy(h); + +julia> b = FockBasis(5); ψ = fockstate(b, 2); + +julia> real(numeric_average(a' * a, ψ)) ≈ 2 +true +``` + +See also [`to_numeric`](@ref), [`average`](@ref). +""" +function numeric_average end + +numeric_average(op, state::QuantumState, d::AbstractDict{<:QSym} = _NO_SUBS) = + _numeric_average(op, state, d) +function numeric_average(op, states::AbstractVector, d::AbstractDict{<:QSym} = _NO_SUBS) + isempty(states) && throw(ArgumentError("numeric_average: states vector is empty")) + return _numeric_average_vec(op, states, d) +end + +# Leaf: assemble the lazy numeric operator on the state's basis (never materialized, so a +# `LazyKet` state works) and take the backend expectation. +_numeric_average(op::QField, state, d::AbstractDict{<:QSym}) = + numeric_expect(numeric_backend(state), _to_numeric_lazy(op, NumericContext(numeric_backend(state), numeric_basis(state), d)), state) +_numeric_average(x::Number, ::Any, ::AbstractDict{<:QSym}) = _to_complex(x) +_numeric_average(x::Num, state, d::AbstractDict{<:QSym}) = + _numeric_average(SymbolicUtils.unwrap(x), state, d) + +# `arguments(::BasicSymbolic)` is statically `Any`; the `_to_complex` wraps reseal each +# recursive result to `ComplexF64`. +function _numeric_average(avg::SymbolicUtils.BasicSymbolic, state, d::AbstractDict{<:QSym}) + if SymbolicUtils.hasmetadata(avg, AverageOperator) || + (SymbolicUtils.iscall(avg) && SymbolicUtils.operation(avg) isa AvgFunc) || + is_indexed_sum(avg) + return _to_complex(_numeric_average(undo_average(avg), state, d)) + end + SymbolicUtils.isconst(avg) && return _to_complex(_numeric_average(avg.val, state, d)) + SymbolicUtils.iscall(avg) || throw(ArgumentError("numeric_average not implemented for $(typeof(avg))")) + f, args = SymbolicUtils.operation(avg), SymbolicUtils.arguments(avg) + f === (^) && return _to_complex(_numeric_average(args[1], state, d))^_to_complex(args[2]) + f === (+) && return _sum_avg(args, state, d) + f === (*) && return _prod_avg(args, state, d) + throw(ArgumentError("numeric_average not implemented for operation $f")) +end + +function _sum_avg(args, state, d::AbstractDict{<:QSym}) + r = ComplexF64(0) + for a in args + r += _to_complex(_numeric_average(a, state, d)) + end + return r +end +function _prod_avg(args, state, d::AbstractDict{<:QSym}) + r = ComplexF64(1) + for a in args + r *= _to_complex(_numeric_average(a, state, d)) + end + return r +end + +# Operator expressions assemble lazily once and expect over all states (the vector +# `numeric_expect` materialises internally); symbolic averages map. +_numeric_average_vec(op::QField, states, d::AbstractDict{<:QSym}) = + numeric_expect(numeric_backend(first(states)), _to_numeric_lazy(op, NumericContext(numeric_backend(first(states)), numeric_basis(first(states)), d)), states) +_numeric_average_vec(op, states, d::AbstractDict{<:QSym}) = + [_numeric_average(op, ψ, d) for ψ in states] + +""" + expect(op::QField, state[, d::AbstractDict{<:QSym}]) + +Alias for [`numeric_average`](@ref) on operator expressions. For symbolic scalar +expressions such as `average(op)`, call [`numeric_average`](@ref) directly. +""" +expect(op::QField, state, d::AbstractDict{<:QSym} = _NO_SUBS) = numeric_average(op, state, d) diff --git a/src/numeric/backend.jl b/src/numeric/backend.jl new file mode 100644 index 00000000..ed3df132 --- /dev/null +++ b/src/numeric/backend.jl @@ -0,0 +1,196 @@ +# Numeric backends and the open extension surface. +# +# `to_numeric`/`numeric_average` are backend-neutral in `src/`; every place that would +# name a concrete numeric type (a `Basis`, a lazy operator, an `expect`) goes through a +# generic hook whose methods live in a package extension (`ext/`). A backend is selected +# by a zero-field singleton so dispatch is static and the backend axis is free for +# inference. Loading `QuantumOpticsBase` or `QuantumToolbox` activates the matching +# extension, which adds the methods below. + +""" + NumericBackend + +Abstract supertype for numeric-conversion backends. A backend is normally a concrete +singleton whose methods implement the exported numeric hooks. See [Adding a numeric +backend](@ref numeric-backend-interface) for the interface. + +Bundled implementations are selected with [`QuantumOpticsBackend`](@ref) or +[`QuantumToolboxBackend`](@ref) after loading the corresponding package. +""" +abstract type NumericBackend end + +""" + QuantumOpticsBackend() + +Select [QuantumOpticsBase.jl](https://github.com/qojulia/QuantumOpticsBase.jl). Active after +`using QuantumOpticsBase`. +""" +struct QuantumOpticsBackend <: NumericBackend end + +""" + QuantumToolboxBackend() + +Select [QuantumToolbox.jl](https://github.com/qutip/QuantumToolbox.jl). Active after +`using QuantumToolbox`. +""" +struct QuantumToolboxBackend <: NumericBackend end + +# --- Open extension hooks (bare generic functions; methods defined in ext/) ------------ +# +# `be` is always a backend singleton. `basis` is opaque to the core: a `Basis` for +# QuantumOptics, integer dims (`Int`/`Vector{Int}`) for QuantumToolbox. + +""" + numeric_operator(be, op::Op, subbasis) + +Map a single-site operator `op` to its numeric matrix on `subbasis` (the subsystem basis +for `op`'s slot). The main extension method is a closed `k === OP_*` value-branch on +`op.kind`; the separate extension point `numeric_operator(be, ::Val{kind}, op, subbasis)` +is the fallthrough for adding backend support for an existing role and throws by default. +""" +function numeric_operator end + +""" + numeric_basis(be, h::HilbertSpace, dims) + +Build the full numeric basis for Hilbert space `h` with user dimensions `dims` (Fock +cutoff / spin number; ignored for spaces that carry their own level count). Returns a +`Basis` (QuantumOptics) or integer dims (QuantumToolbox). +""" +function numeric_basis end + +""" + numeric_subbasis(be, basis, slot::Int) + +The subsystem basis for `slot` of a composite `basis` (the basis itself for a simple one). +""" +function numeric_subbasis end + +""" + numeric_embed(be, basis, slot::Int, leaf) + +Place the single-site numeric operator `leaf` at `slot` of `basis`, returning an operator +on the full space (the leaf unchanged for a simple basis, a lazy/concrete tensor for a +composite one). +""" +function numeric_embed end + +""" + numeric_identity(be, basis) + +The identity operator on `basis` (lazy on a composite basis). +""" +function numeric_identity end + +""" + numeric_num_subsystems(be, basis) -> Int + +Number of tensor-product subsystems represented by `basis`. Backends use this hook to +validate indexed-site layouts before emitting any numeric operators. +""" +function numeric_num_subsystems end + +""" + numeric_assemble(be, basis, terms) + +Assemble a static lazy operator from a backend-neutral `terms` list. Each entry is +`(c::ComplexF64, factors::Vector{})` with `factors` already embedded at their slots. +The result is one concrete vector-backed lazy type for any term count. +""" +function numeric_assemble end + +""" + numeric_assemble_td(be, basis, td_terms) + +Like [`numeric_assemble`](@ref) but each coefficient is `Union{ComplexF64, Function}` +(a constant or a `t -> ComplexF64`), producing the backend's native time-dependent +operator. +""" +function numeric_assemble_td end + +""" + numeric_materialize(be, op, op_type) + +Materialize a static lazily assembled operator `op` according to `op_type`, applied exactly +once at the public boundary. Implementations must treat `nothing` as a request for their +ordinary eager operator and `identity` as a request for `op` unchanged. Other callables are +backend-defined. The eager representation need not be sparse, although both bundled +backends choose sparse storage. + +This hook is not used by [`numeric_average`](@ref), which consumes the lazy assembly, or by +time-dependent conversion, which returns the backend's native time-dependent operator. +""" +function numeric_materialize end + +""" + numeric_expect(be, numop, state) + +Expectation value `⟨state| numop |state⟩` (a `ComplexF64`), or a `Vector{ComplexF64}` when +`state` is a vector of states. +""" +function numeric_expect end + +""" + numeric_backend(state) -> NumericBackend + +Return the numeric backend owning `state`. Third-party backends must implement this hook +for their state types to support [`numeric_average`](@ref) and state-based [`to_numeric`](@ref). +""" +function numeric_backend end + +# The one-argument `numeric_basis(state)` is the public state-to-basis counterpart of the +# three-argument Hilbert-space builder declared above. +""" + numeric_basis(state) + +Return the backend basis/dimensions carried by `state`. Third-party backends must implement +this together with [`numeric_backend`](@ref) for state-based numeric conversion. +""" +numeric_basis(s::Union{StateVector, AbstractOperator}) = basis(s) + +function _check_product_dims(h::ProductSpace, dims) + nspaces = length(h.spaces) + ndims = try + length(dims) + catch + throw( + ArgumentError( + "dims for a ProductSpace must be an indexable collection with $nspaces entries", + ), + ) + end + ndims == nspaces || throw( + ArgumentError( + "dims has $ndims entries, but the ProductSpace has $nspaces subspaces", + ), + ) + return dims +end + +""" + _default_backend() + +The backend to use when none is given and no backend object pins one. Resolves to the sole +loaded backend extension, erroring helpfully when zero or both are loaded. +""" +function _default_backend() + m = @__MODULE__ + qob = Base.get_extension(m, :SecondQuantizedAlgebraQuantumOpticsBaseExt) !== nothing + qtb = Base.get_extension(m, :SecondQuantizedAlgebraQuantumToolboxExt) !== nothing + qob && !qtb && return QuantumOpticsBackend() + qtb && !qob && return QuantumToolboxBackend() + if qob && qtb + throw( + ArgumentError( + "both numeric backends are loaded; pass `backend = QuantumOpticsBackend()` " * + "or `backend = QuantumToolboxBackend()` explicitly.", + ), + ) + end + throw( + ArgumentError( + "no numeric backend loaded; run `using QuantumOpticsBase` or `using QuantumToolbox` " * + "to enable `to_numeric`/`numeric_average`.", + ), + ) +end diff --git a/src/numeric/coeff.jl b/src/numeric/coeff.jl new file mode 100644 index 00000000..d81b54b7 --- /dev/null +++ b/src/numeric/coeff.jl @@ -0,0 +1,213 @@ +# Coefficient lowering: reduce a (possibly symbolic) `CNum`/`Num` prefactor to a concrete +# `ComplexF64`, split user parameter/time substitutions, and compile time-dependent +# coefficient functions. Entirely backend-independent. + +const _NO_SCALAR_SUBS = Dict{Num, Any}() + +# --- parameter / time-parameter normalisation ----------------------------------------- + +function _expand_parameter(parameter) + isempty(parameter) && return parameter + out = Dict{Any, Any}() + for (k, v) in parameter + if k isa Complex + out[real(k)] = real(v) + ik = imag(k) + SymbolicUtils.unwrap(ik) isa Number || (out[ik] = imag(v)) + else + out[k] = v + end + end + return out +end + +function _normalize_time_parameter(time_parameter) + isempty(time_parameter) && return time_parameter + out = Dict{Any, Any}() + for (k, v) in time_parameter + out[k] = v isa Number ? (t -> v) : v + end + return out +end + +function _time_basis(time_parameter) + basevars = Any[] + valuefuncs = Any[] + for (k, f) in time_parameter + uk = SymbolicUtils.unwrap(k) + if SymbolicUtils.issym(uk) + push!(basevars, k) + push!(valuefuncs, f) + continue + end + vs = collect(Symbolics.get_variables(k)) + length(vs) == 1 || throw( + ArgumentError( + "time_parameter key `$k` must depend on exactly one variable, got $(length(vs)).", + ), + ) + wv = Symbolics.wrap(vs[1]) + if isequal(uk, SymbolicUtils.unwrap(conj(wv))) + push!(basevars, wv) + push!(valuefuncs, t -> conj(f(t))) + else + throw( + ArgumentError( + "unsupported time_parameter key `$k`; only a bare variable `v` or `conj(v)` are supported.", + ), + ) + end + end + return basevars, Tuple(valuefuncs) +end + +_coeff_is_const(c::Complex{Num}) = + isempty(Symbolics.get_variables(real(c))) && isempty(Symbolics.get_variables(imag(c))) + +_const_coeff(c::Complex{Num}) = _to_complex(c) + +_as_cnum(x::Complex) = Complex{Num}(Num(real(x)), Num(imag(x))) +_as_cnum(x) = Complex{Num}(Num(x), Num(false)) + +function _coefficient_variables(c::Complex{Num}) + vars = Any[] + append!(vars, Symbolics.get_variables(real(c))) + append!(vars, Symbolics.get_variables(imag(c))) + unique!(vars) + return vars +end + +function _check_time_variables(c::Complex{Num}, basevars) + vars = _coefficient_variables(c) + base_unwrapped = SymbolicUtils.unwrap.(basevars) + missing = Any[] + for v in vars + any(b -> isequal(v, b), base_unwrapped) || push!(missing, v) + end + isempty(missing) && return nothing + throw( + ArgumentError( + "time-dependent coefficient `$c` depends on variables without time values: " * + join(string.(missing), ", "), + ), + ) +end + +function _compile_coeff(c::Complex{Num}, vars...) + f_re = build_function(real(c), vars...; expression = Val(false)) + f_im = build_function(imag(c), vars...; expression = Val(false)) + g_re = f_re isa Tuple ? first(f_re) : f_re + g_im = f_im isa Tuple ? first(f_im) : f_im + return (vals...) -> g_re(vals...) + im * g_im(vals...) +end + +# --- indexed scalar substitutions (used by the indexed unroll) ------------------------ + +# Split user-supplied scalar substitutions into real and imag legs once per +# `to_numeric` call. A complex RHS `re + im*ii` propagates as separate real +# and imaginary substitutions, preserving the `Complex{Num}` invariant. +function _split_scalar_subs(scalar_subs::AbstractDict) + sub_re = Dict{Any, Any}() + sub_im = Dict{Any, Any}() + has_imag = false + for (k, v) in scalar_subs + kraw = SymbolicUtils.unwrap(k) + if v isa Complex + sub_re[kraw] = real(v) + sub_im[kraw] = imag(v) + iszero(imag(v)) || (has_imag = true) + else + sub_re[kraw] = v + sub_im[kraw] = 0 + end + end + return sub_re, sub_im, has_imag +end + +function _apply_scalar_subs(c::CNum, sub_re::Dict, sub_im::Dict, has_imag::Bool) + (isempty(sub_re) || _is_native(c)) && return c + re_part = SymbolicUtils.unwrap(real(c)) + im_part = SymbolicUtils.unwrap(imag(c)) + if !has_imag + return _cnum( + Num(Symbolics.substitute(re_part, sub_re)), + Num(Symbolics.substitute(im_part, sub_re)), + ) + end + # (re + i*im) -> (re|sub_re - im|sub_im) + i*(re|sub_im + im|sub_re) + new_re = Symbolics.substitute(re_part, sub_re) - + Symbolics.substitute(im_part, sub_im) + new_im = Symbolics.substitute(re_part, sub_im) + + Symbolics.substitute(im_part, sub_re) + return _cnum(Num(new_re), Num(new_im)) +end + +# --- concrete reduction ---------------------------------------------------------------- + +function _reduce_const(n::Num)::ComplexF64 + x = Symbolics.value(n) + try + return _fold_const(x) + catch err + err isa ArgumentError || rethrow() + isempty(Symbolics.get_variables(n)) || rethrow() + return _compile_const(n) + end +end + +_compile_const(n::Num)::ComplexF64 = ComplexF64(symbolic_to_float(n)) + +function _fold_const(x)::ComplexF64 + x isa Number && return x + if x isa SymbolicUtils.BasicSymbolic + if SymbolicUtils.iscall(x) + op = SymbolicUtils.operation(x) + args = SymbolicUtils.arguments(x) + if op === (+) + acc = zero(ComplexF64) + for a in args + acc += _fold_const(a) + end + return acc + elseif op === (*) + acc = one(ComplexF64) + for a in args + acc *= _fold_const(a) + end + return acc + elseif op === conj + return conj(_fold_const(first(args))) + elseif op === real + return real(_fold_const(first(args))) + elseif op === imag + return imag(_fold_const(first(args))) + elseif op === (/) + return _fold_const(first(args)) / _fold_const(last(args)) + elseif op === (^) + return _fold_const(first(args))^_fold_const(last(args)) + elseif op === (-) + return length(args) == 1 ? -_fold_const(first(args)) : + _fold_const(first(args)) - _fold_const(last(args)) + end + elseif SymbolicUtils.isconst(x) + return x.val + end + end + throw(ArgumentError("cannot reduce symbolic expression $x to a concrete number")) +end + +_to_complex(c::Coeff) = _is_native(c) ? c.z : _to_complex(to_num(c)) + +# One method (union-split budget) routing every input through `convert ∘ Complex`, +# the only pattern that infers to ComplexF64 from `Any` after `Symbolics.value`. +function _to_complex(x) + x isa ComplexF64 && return x + x isa Complex{Num} && return _reduce_const(real(x)) + im * _reduce_const(imag(x)) + x isa Complex && return convert(ComplexF64, x) + x isa Num && return _reduce_const(x) + if x isa SymbolicUtils.BasicSymbolic + SymbolicUtils.isconst(x) || throw(ArgumentError("cannot reduce symbolic expression $x to a concrete number")) + return convert(ComplexF64, Complex(x.val, false)) + end + return convert(ComplexF64, Complex(x, false)) +end diff --git a/src/numeric/core.jl b/src/numeric/core.jl new file mode 100644 index 00000000..a3f1070d --- /dev/null +++ b/src/numeric/core.jl @@ -0,0 +1,161 @@ +# Backend-neutral core: walk a canonical `QAdd` into a backend-neutral term list and hand +# it to the backend's lazy assembler. The core never does operator arithmetic and never +# names a concrete numeric type; everything numeric goes through the hooks in `backend.jl`. + +# `Union{}` value keeps the haskey-fallback branch concrete-typed (a real op_subs dict +# widens the leaf to `Any`, which only the kwargs/dict paths use and which are not on the +# `@inferred` contract). +const _NO_SUBS = Dict{QSym, Union{}}() +const _EMPTY_SITES = Dict{Int, Vector{Int}}() + +""" + NumericContext(backend, basis, op_subs, scalar_subs, sites) + +Immutable, all-concrete-fields bundle threaded through the numeric core: the backend +singleton, the opaque per-backend `basis`, the operator-substitution dict `op_subs`, the +indexed `scalar_subs`, and the indexed `sites` map (empty for the non-indexed path). All +fields are concrete via the type parameters, so the context never widens inference. +""" +struct NumericContext{BE <: NumericBackend, B, SUBS <: AbstractDict, SS <: AbstractDict} + backend::BE + basis::B + op_subs::SUBS + scalar_subs::SS + sites::Dict{Int, Vector{Int}} +end + +NumericContext(be::NumericBackend, basis, op_subs::AbstractDict = _NO_SUBS) = + NumericContext(be, basis, op_subs, _NO_SCALAR_SUBS, _EMPTY_SITES) + +# --- backend resolution for QuantumOptics states (QI types, so it lives here) ----------- +# A `Basis` is inherently QuantumOptics, so the Basis-form `to_numeric` hardcodes +# `QuantumOpticsBackend()`. The QuantumToolbox extension adds its state method. +numeric_backend(::Union{StateVector, AbstractOperator}) = QuantumOpticsBackend() + +# --- leaf / term construction ---------------------------------------------------------- + +# The numeric leaf for one operator: a custom override from `op_subs`, else the backend +# operator embedded at its slot. With `op_subs = _NO_SUBS` the override branch has value +# type `Union{}` and contributes nothing, so the inferred type is exactly the embed type. +function _numeric_leaf(op::Op, ctx::NumericContext) + haskey(ctx.op_subs, op) && return ctx.op_subs[op] + slot = Int(op.space_index) + sub = numeric_subbasis(ctx.backend, ctx.basis, slot) + leaf = numeric_operator(ctx.backend, op, sub) + return numeric_embed(ctx.backend, ctx.basis, slot, leaf) +end + +_term_factors(ops::Vector{Op}, ctx::NumericContext) = [_numeric_leaf(op, ctx) for op in ops] + +# --- static assembly ------------------------------------------------------------------- + +# Build the (ComplexF64, factors) term list and assemble. Coefficients must already be +# concrete (the kwargs path substitutes `parameter` before calling this); `_to_complex` +# throws `ArgumentError` for a still-symbolic coefficient. +function _to_numeric_static(q::QAdd, ctx::NumericContext) + terms = [(_to_complex(c), _term_factors(term.ops, ctx)) for (term, c) in q.arguments] + return numeric_assemble(ctx.backend, ctx.basis, terms) +end + +# Lazy builders used by `numeric_average`/`expect`: assemble but never materialize, so a +# `LazyKet` state works and no concrete matrix is built just to take an expectation value. +_to_numeric_lazy(op::Op, ctx::NumericContext) = _numeric_leaf(op, ctx) +_to_numeric_lazy(q::QAdd, ctx::NumericContext) = _to_numeric_static(q, ctx) + +# --- time-dependent assembly ----------------------------------------------------------- + +# Each coefficient becomes a `ComplexF64` (constant term) or a `t -> ComplexF64` closure +# (genuinely time-dependent), reusing the kept coefficient machinery. Factors are stored +# behind `Vector{Any}` because the TD path is not on the `@inferred` contract and the +# backend assembler only needs to iterate them. +const _TDTerm = Tuple{Union{ComplexF64, Function}, Vector{Any}} + +function _td_coeff(c::Complex{Num}, basevars, valuefuncs) + if _coeff_is_const(c) + return _const_coeff(c) + end + _check_time_variables(c, basevars) + pref = _compile_coeff(c, basevars...) + return let pref = pref, valuefuncs = valuefuncs + t -> pref(map(f -> f(t), valuefuncs)...) + end +end + +function _to_numeric_td(q::QAdd, ctx::NumericContext, time_parameter) + basevars, valuefuncs = _time_basis(time_parameter) + td_terms = _TDTerm[] + for (term, c_) in q.arguments + coef = _td_coeff(to_num(c_), basevars, valuefuncs) + push!(td_terms, (coef, collect(Any, _term_factors(term.ops, ctx)))) + end + isempty(td_terms) && push!(td_terms, (0.0im, Any[])) + return numeric_assemble_td(ctx.backend, ctx.basis, td_terms) +end + +# --- keyword translation boundary ------------------------------------------------------ + +function _numeric_operator_dict(operators, adjoint_ops::Bool) + out = Dict{QSym, Any}() + for (k, v) in operators + k isa QSym || throw(ArgumentError("operator substitution key `$k` is not a QSym")) + out[k] = v + end + if adjoint_ops + for (k, v) in operators + k isa QSym || continue + k_adj = Base.adjoint(k) + haskey(out, k_adj) || (out[k_adj] = Base.adjoint(v)) + end + end + return out +end + +function _to_numeric_kw(be::NumericBackend, op, b; parameter, time_parameter, operators, adjoint_ops, op_type) + param = _expand_parameter(parameter) + tp = _normalize_time_parameter(time_parameter) + ops = _numeric_operator_dict(operators, adjoint_ops) + ctx = NumericContext(be, b, ops) + return _to_numeric_translated(op, ctx, param, tp, op_type) +end + +# QSym: wrap as a one-term QAdd and reuse the QAdd path. +_to_numeric_translated(op::QSym, ctx::NumericContext, parameter, time_parameter, op_type) = + _to_numeric_translated(_single_qadd(_CNUM_ONE, Op[op]), ctx, parameter, time_parameter, op_type) + +# Static path materializes via `op_type` (`nothing` requests the backend's ordinary eager +# operator); the time-dependent path returns the backend's native TD operator. +function _to_numeric_translated(op::QAdd, ctx::NumericContext, parameter, time_parameter, op_type) + op_ = substitute(op, parameter) + isempty(time_parameter) && return numeric_materialize(ctx.backend, _to_numeric_static(op_, ctx), op_type) + _check_td_op_type(op_type) + return _to_numeric_td(op_, ctx, time_parameter) +end + +# Bare scalar: static -> scaled identity; time-dependent -> native TD over identity. +function _to_numeric_translated(arg, ctx::NumericContext, parameter, time_parameter, op_type) + arg_sub = substitute(arg, parameter) + c = _as_cnum(arg_sub) + if isempty(time_parameter) + _coeff_is_const(c) || throw( + ArgumentError( + "cannot translate symbolic scalar `$arg_sub` without a value: supply it via " * + "`parameter` or `time_parameter`.", + ), + ) + return numeric_materialize(ctx.backend, _const_coeff(c) * numeric_identity(ctx.backend, ctx.basis), op_type) + end + _check_td_op_type(op_type) + basevars, valuefuncs = _time_basis(time_parameter) + coef = _td_coeff(c, basevars, valuefuncs) + return numeric_assemble_td(ctx.backend, ctx.basis, _TDTerm[(coef, Any[])]) +end + +function _check_td_op_type(op_type) + (op_type === nothing || op_type === identity) && return nothing + throw( + ArgumentError( + "op_type is not supported when time_parameter is non-empty; native " * + "time-dependent operators cannot be materialized once at conversion time", + ), + ) +end diff --git a/src/numeric/indexed.jl b/src/numeric/indexed.jl new file mode 100644 index 00000000..89562590 --- /dev/null +++ b/src/numeric/indexed.jl @@ -0,0 +1,158 @@ +# Indexed numeric conversion. The unroll logic (which bound indices a term depends on, +# the `ne` diagonal filter, the `sub_op`/`sub_coef` dual substitution) is backend-neutral +# and kept verbatim; only the leaf emission is routed through the backend hooks at the +# resolved site, so every combination becomes one entry of the same term list the static +# path assembles. + +""" + to_numeric(q::QAdd, basis, sites::Dict{Int, Vector{Int}}[, d[, scalar_subs]]) + +Convert a `QAdd` that carries bound summation indices (`q.indices`) to a numeric operator +on a `basis` whose layout replicates one or more abstract subspaces. `sites[orig_space_index]` +is the list of slots that realise that subspace; non-indexed subspaces map to a single-slot +entry. + +For each term, every bound index the term depends on is unrolled over the length of its +`sites` vector, with `term.ne` constraints filtering out diagonal combinations. Concrete-site +operators are routed to `sites[op.space_index][k]`, where `k` is the substituted integer +site. Terms that do not depend on any bound index are emitted once as written, matching the +`Σ` convention that i-independent residuals already carry their range factor. + +`d` substitutes individual `QSym` operators with custom numeric operators. `scalar_subs` +substitutes symbolic scalar parameters (e.g. `Δ` or `g(k)` for the integer `k` produced by +index unrolling) with numeric values. +""" +function to_numeric( + q::QAdd, b::Basis, + sites::AbstractDict{Int, Vector{Int}}, + d::AbstractDict{<:QSym} = _NO_SUBS, + scalar_subs::AbstractDict = _NO_SCALAR_SUBS, + ) + return _to_numeric_indexed(QuantumOpticsBackend(), b, q, sites, d, scalar_subs) +end + +# Backend-neutral indexed conversion. Every place that would name a numeric type goes +# through the hooks, so this drives both backends; the `Basis`/QTB-dims entry points only +# pick the backend singleton and hand off here. +function _to_numeric_indexed( + be::NumericBackend, basis, q::QAdd, + sites::AbstractDict{Int, Vector{Int}}, + d::AbstractDict{<:QSym}, scalar_subs::AbstractDict, + ) + isempty(q.indices) && return to_numeric(q, basis, d) + _validate_sites(be, basis, sites) + ctx = NumericContext(be, basis, d, scalar_subs, Dict{Int, Vector{Int}}(sites)) + sub_re, sub_im, has_imag = _split_scalar_subs(scalar_subs) + terms = Tuple{ComplexF64, Vector{Any}}[] + for (term, c) in q.arguments + _accumulate_indexed_term!(terms, term, c, q.indices, ctx, sub_re, sub_im, has_imag) + end + return numeric_materialize(be, numeric_assemble(be, basis, terms), nothing) +end + +function _validate_sites(be::NumericBackend, basis, sites) + nsubsystems = numeric_num_subsystems(be, basis) + seen = Set{Int}() + for (space, slots) in sites + for slot in slots + 1 <= slot <= nsubsystems || throw( + ArgumentError( + "sites[$space] contains slot $slot, but the numeric basis has " * + "$nsubsystems subsystem$(nsubsystems == 1 ? "" : "s")", + ), + ) + slot in seen && throw( + ArgumentError( + "numeric subsystem slot $slot occurs more than once in the sites layout", + ), + ) + push!(seen, slot) + end + end + return nothing +end + +function _numeric_leaf_indexed(op::Op, ctx::NumericContext) + haskey(ctx.op_subs, op) && return ctx.op_subs[op] + slot = _resolve_slot(op, ctx.sites) + sub = numeric_subbasis(ctx.backend, ctx.basis, slot) + leaf = numeric_operator(ctx.backend, op, sub) + return numeric_embed(ctx.backend, ctx.basis, slot, leaf) +end + +function _push_indexed_combo!(terms, ops::Vector{Op}, c::CNum, ctx::NumericContext) + factors = Any[_numeric_leaf_indexed(op, ctx) for op in ops] + push!(terms, (_to_complex(c), factors)) + return terms +end + +function _accumulate_indexed_term!( + terms, term::QTerm, c::CNum, indices::Vector{Index}, + ctx::NumericContext, sub_re::Dict, sub_im::Dict, has_imag::Bool, + ) + dep_indices = Index[idx for idx in indices if _depends_on_index_term(c, term.ops, idx)] + if isempty(dep_indices) + c_resolved = _apply_scalar_subs(c, sub_re, sub_im, has_imag) + return _push_indexed_combo!(terms, term.ops, c_resolved, ctx) + end + lens = Int[length(ctx.sites[Int(idx.space_index)]) for idx in dep_indices] + total = prod(lens) + sub_op = Dict{Index, Index}() + sub_coef = Dict{Index, Index}() + for combo in 1:total + empty!(sub_op) + empty!(sub_coef) + rem = combo - 1 + for k in 1:length(dep_indices) + kpos = (rem % lens[k]) + 1 + rem ÷= lens[k] + idx = dep_indices[k] + # Operators keep the index name (slot kpos drives routing / ne checks and + # lets a resolved op still match a user `d` key); coefficients use the + # anonymous name_id-0 form so `index_sym` is Num(kpos), resolving g(i)→g(k). + sub_op[idx] = Index(idx.name_id, idx.range_id, idx.space_index, Int32(kpos)) + sub_coef[idx] = Index(Int32(0), idx.range_id, idx.space_index, Int32(kpos)) + end + _violates_ne(term.ne, sub_op) && continue + new_ops = Op[change_index(op, sub_op) for op in term.ops] + new_c = change_index(c, sub_coef) + new_c = _apply_scalar_subs(new_c, sub_re, sub_im, has_imag) + _push_indexed_combo!(terms, new_ops, new_c, ctx) + end + return terms +end + +function _violates_ne(ne::Vector{NonEqualPair}, sub::Dict{Index, Index}) + for (a, b) in ne + ra = get(sub, a, a) + rb = get(sub, b, b) + va = index_slot(ra) + vb = index_slot(rb) + va === nothing && continue + vb === nothing && continue + ra.space_index == rb.space_index || continue + va == vb && return true + end + return false +end + +function _resolve_slot(op::QSym, sites::AbstractDict{Int, Vector{Int}}) + si = Int(op.space_index) + slots = get(sites, si, Int[]) + if isempty(slots) + return si + end + if has_index(op.index) + v = index_slot(op.index) + if v !== nothing && 1 <= v <= length(slots) + return slots[v] + end + end + length(slots) == 1 && return slots[1] + throw( + ArgumentError( + "cannot resolve slot for operator $(op): index is not a concrete integer site, " * + "and sites[$si] has $(length(slots)) candidates", + ), + ) +end diff --git a/test/Project.toml b/test/Project.toml index 543f6ea6..7778b301 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,11 +4,14 @@ CheckConcreteStructs = "73c92db5-9da6-4938-911a-6443a7e94a58" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc" QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678" +QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" SecondQuantizedAlgebra = "f7aa4685-e143-4cb6-a7f3-073579757907" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" @@ -24,7 +27,9 @@ LaTeXStrings = "1" Latexify = "0.13, 0.14, 0.15, 0.16" MutableArithmetics = "1" ParallelTestRunner = "2" -QuantumOpticsBase = "0.4, 0.5" +QuantumOpticsBase = "0.5" +QuantumToolbox = "0.47" +SciMLOperators = "1.22" SymbolicUtils = "4" Symbolics = "7" TermInterface = "2" diff --git a/test/average_test.jl b/test/average_test.jl index 0bf6c222..4283c61f 100644 --- a/test/average_test.jl +++ b/test/average_test.jl @@ -851,6 +851,38 @@ import SecondQuantizedAlgebra: simplify, QAdd, QSym, QField, CNum, _to_cnum, _si @test symtype(SymbolicUtils.unwrap(lifted_h)) === Real @test symtype(SymbolicUtils.unwrap(lifted_n)) === Number end + + @testset "lifted QAdd variable name encodes NE constraints" begin + # A QAdd term carrying an index-inequality pair must contribute an `_ne…` + # token to the auto-generated time-dependent variable name, so two terms + # differing only in their NE scope lift to distinct unknowns. + ha = NLevelSpace(:atom, 2, 1) + hf = FockSpace(:f) + hp = hf ⊗ ha + @variables Nsites t + iv = SymbolicUtils.unwrap(t) + i = Index(hp, :i, Nsites, ha) + j = Index(hp, :j, Nsites, ha) + σ_i = IndexedOperator(Transition(hp, :σ, 2, 2, 2), i) + σ_j = IndexedOperator(Transition(hp, :σ, 2, 2, 2), j) + + with_ne = QAdd( + QTermDict(QTerm(Op[σ_i, σ_j], NonEqualPair[(i, j)]) => _to_cnum(1)), + Index[], + ) + without_ne = QAdd( + QTermDict(QTerm(Op[σ_i, σ_j], NonEqualPair[]) => _to_cnum(1)), + Index[], + ) + + lifted_ne = make_time_dependent(average(with_ne), iv) + lifted_plain = make_time_dependent(average(without_ne), iv) + name_ne = string(SymbolicUtils.operation(SymbolicUtils.unwrap(lifted_ne))) + name_plain = string(SymbolicUtils.operation(SymbolicUtils.unwrap(lifted_plain))) + @test occursin("_ne", name_ne) + @test !occursin("_ne", name_plain) + @test name_ne != name_plain + end end end diff --git a/test/numeric_indexed_test.jl b/test/numeric_indexed_test.jl index 1faaeca2..f1760a85 100644 --- a/test/numeric_indexed_test.jl +++ b/test/numeric_indexed_test.jl @@ -141,3 +141,102 @@ end M_over = to_numeric(H, b, sites, Dict{S.QSym, Any}(sig(1, 2, i) => custom)) @test sparse(M_plain).data != sparse(M_over).data end + +@testset "Indexed numeric: ne-constrained double sum" begin + # A double sum `Σ_{i,j} σ⁺_i σ⁻_j` splits into a diagonal part (`i == j`) and an + # off-diagonal part carrying an `ne = (i, j)` constraint that `_violates_ne` uses to + # skip the `i == j` combinations during the unroll. The result must equal the full + # explicit sum over all site pairs. + ha = NLevelSpace(:atom, 2) + σ(α, β, k) = IndexedOperator(Transition(ha, :σ, α, β), k) + i = Index(ha, :i, 2, ha) + j = Index(ha, :j, 2, ha) + H = Σ(σ(1, 2, i) * σ(2, 1, j), i, j) + + bn = NLevelBasis(2) + b = bn ⊗ bn + sites = Dict{Int, Vector{Int}}(1 => [1, 2]) + M = to_numeric(H, b, sites, Dict{S.QSym, Any}()) + + he = NLevelSpace(:a1, 2) ⊗ NLevelSpace(:a2, 2) + t(α, β, k) = Transition(he, Symbol(:s, k), α, β, k) + Hex = t(1, 2, 1) * t(2, 1, 1) + t(1, 2, 1) * t(2, 1, 2) + + t(1, 2, 2) * t(2, 1, 1) + t(1, 2, 2) * t(2, 1, 2) + Mex = to_numeric(Hex, b) + @test _matnorm(sparse(M).data, sparse(Mex).data) < 1.0e-10 +end + +@testset "Indexed numeric: real-valued scalar_subs" begin + # Real (non-`Complex`-typed) `scalar_subs` values exercise the real-only substitution + # legs of `_split_scalar_subs`/`_apply_scalar_subs` (the `has_imag == false` path). + hc = FockSpace(:cavity) + ha = NLevelSpace(:atom, 2) + h = hc ⊗ ha + a = Destroy(h, :a, 1) + σ(α, β, k) = IndexedOperator(Transition(h, :σ, α, β, 2), k) + gv(k) = IndexedVariable(:g, k) + i = Index(h, :i, 2, ha) + H = Σ(gv(i) * (a' * σ(1, 2, i) + a * σ(2, 1, i)), i) + + bc = FockBasis(2) + bn = NLevelBasis(2) + b = bc ⊗ bn ⊗ bn + sites = Dict{Int, Vector{Int}}(1 => [1], 2 => [2, 3]) + + ks_sym = SymbolicUtils.Sym{SymbolicUtils.SymReal}( + :g; + type = SymbolicUtils.FnType{Tuple{Int}, Real, Nothing}, + shape = UnitRange{Int}[], + ) + scalar_subs = Dict{Num, Float64}(Num(ks_sym(1)) => 0.3, Num(ks_sym(2)) => 0.7) + M = to_numeric(H, b, sites, Dict{S.QSym, Any}(), scalar_subs) + + # Explicit reference on the same layout. + he = FockSpace(:c) ⊗ NLevelSpace(:a1, 2) ⊗ NLevelSpace(:a2, 2) + ae = Destroy(he, :a, 1) + te(α, β, k) = Transition(he, Symbol(:s, k), α, β, k + 1) + Hex = 0.3 * (ae' * te(1, 2, 1) + ae * te(2, 1, 1)) + + 0.7 * (ae' * te(1, 2, 2) + ae * te(2, 1, 2)) + Mex = to_numeric(Hex, b) + @test _matnorm(sparse(M).data, sparse(Mex).data) < 1.0e-10 +end + +@testset "Indexed numeric: slot resolution edge cases" begin + hc = FockSpace(:cavity) + ha = NLevelSpace(:atom, 2) + h = hc ⊗ ha + a = Destroy(h, :a, 1) + σ(α, β, k) = IndexedOperator(Transition(h, :σ, α, β, 2), k) + i = Index(h, :i, 2, ha) + H = Σ(a' * σ(1, 2, i) + a * σ(2, 1, i), i) + + bc = FockBasis(2) + bn = NLevelBasis(2) + b = bc ⊗ bn ⊗ bn + sites = Dict{Int, Vector{Int}}(1 => [1], 2 => [2, 3]) + + # A non-indexed subspace (the cavity) may be omitted from `sites`: its operator falls + # back to its own `space_index` slot, so the result is unchanged. + sites_no_cavity = Dict{Int, Vector{Int}}(2 => [2, 3]) + M_full = to_numeric(H, b, sites, Dict{S.QSym, Any}()) + M_omit = to_numeric(H, b, sites_no_cavity, Dict{S.QSym, Any}()) + @test _matnorm(sparse(M_full).data, sparse(M_omit).data) < 1.0e-10 + + # An operator carrying an unresolved (non-integer) index on a multi-slot subspace is + # ambiguous and must error. `σ_j` never binds to the summed index `i`. + j = Index(h, :j, 2, ha) + Hbad = Σ(σ(1, 2, i) * σ(1, 2, j), i) + @test_throws ArgumentError to_numeric(Hbad, b, sites, Dict{S.QSym, Any}()) +end + +@testset "Indexed numeric: invalid site layouts are rejected" begin + h = FockSpace(:site) + a = Destroy(h, :a) + i = Index(h, :i, 2, h) + H = Σ(IndexedOperator(a, i), i) + + @test_throws ArgumentError to_numeric(H, FockBasis(2), Dict(1 => [1, 2])) + b = FockBasis(2) ⊗ FockBasis(2) + @test_throws ArgumentError to_numeric(H, b, Dict(1 => [1, 1])) + @test_throws ArgumentError to_numeric(H, b, Dict(1 => [1, 3])) +end diff --git a/test/numeric_qtb_test.jl b/test/numeric_qtb_test.jl new file mode 100644 index 00000000..5f0fde87 --- /dev/null +++ b/test/numeric_qtb_test.jl @@ -0,0 +1,312 @@ +# QuantumToolbox backend for `to_numeric`/`numeric_average`. +# +# `QuantumToolbox` is imported (not `using`) so its `⊗`/`tensor`/`expect` do not collide +# with the QuantumInterface ones that `SecondQuantizedAlgebra` extends; `⊗` here is the +# shared QuantumInterface operator (works on Hilbert spaces and QuantumOptics bases). Parity +# is checked through `expect`/`numeric_average` scalars (never raw matrices, whose kron +# convention differs across backends), with states built per-backend in slot order. + +using SecondQuantizedAlgebra +import QuantumToolbox as QTB +import QuantumOpticsBase as QOB +import SciMLOperators as SO +using Symbolics: @variables +import LinearAlgebra: tr, norm, diag +using Test + +const QTBB = QuantumToolboxBackend() + +@testset "QuantumToolbox backend" begin + @testset "eager default; lazy assembly is n-stable" begin + h = FockSpace(:f) + @qnumbers a::Destroy(h) + H2 = to_numeric(a' * a, h, 5; backend = QTBB, op_type = identity) + H3 = to_numeric(a' * a + a + a', h, 5; backend = QTBB, op_type = identity) + @variables E::Number + Htd = to_numeric(E * a + conj(E) * a', h, 5; backend = QTBB, time_parameter = Dict(E => t -> 1 + im * t)) + + # The vector-backed VecSum gives ONE concrete type for any term count and for static + # vs time-dependent (the property the custom VecSum buys over SciMLOperators' + # Tuple-locked AddedOperator). + @test H2 isa QTB.QuantumObjectEvolution + @test typeof(H2) === typeof(H3) + @test typeof(H2) === typeof(Htd) + # Static defaults are eager and sparse on both direct and uniform forms; lazy is + # an explicit representation choice. + @test to_numeric(a, 5) isa QTB.QuantumObject + @test typeof(to_numeric(a, 5).data) === typeof(QTB.destroy(5).data) + @test to_numeric(a' * a, h, 5; backend = QTBB) isa QTB.QuantumObject + @test typeof(to_numeric(a' * a, h, 5; backend = QTBB).data) === + typeof(QTB.destroy(6).data) + end + + @testset "Fock parity vs QuantumOptics" begin + h = FockSpace(:f) + @qnumbers a::Destroy(h) + ψq = QOB.coherentstate(QOB.FockBasis(7), 0.4 + 0.2im) + ψt = QTB.coherent(8, 0.4 + 0.2im) + for op in (a, a', a' * a, a + a', 2 * a' * a + a) + @test numeric_average(op, ψq) ≈ numeric_average(op, ψt) atol = 1.0e-6 + end + end + + @testset "Transition parity (1-indexed vs 0-indexed)" begin + h = NLevelSpace(:atom, 3, 1) + σ(i, j) = Transition(h, :σ, i, j) + ψq = QOB.nlevelstate(QOB.NLevelBasis(3), 1) + + 2 * QOB.nlevelstate(QOB.NLevelBasis(3), 2) + + 3 * QOB.nlevelstate(QOB.NLevelBasis(3), 3) + ψq = ψq / norm(ψq.data) + ψt = QTB.basis(3, 0) + 2 * QTB.basis(3, 1) + 3 * QTB.basis(3, 2) + ψt = ψt / norm(ψt) + for i in 1:3, j in 1:3 + @test numeric_average(σ(i, j), ψq) ≈ numeric_average(σ(i, j), ψt) atol = 1.0e-8 + end + end + + @testset "Pauli / Spin parity" begin + hp = PauliSpace(:p) + ψq = (QOB.spindown(QOB.SpinBasis(1 // 2)) + QOB.spinup(QOB.SpinBasis(1 // 2))) / sqrt(2) + ψt = (QTB.basis(2, 1) + QTB.basis(2, 0)) / sqrt(2) + for ax in 1:3 + @test numeric_average(Pauli(hp, :σ, ax), ψq) ≈ numeric_average(Pauli(hp, :σ, ax), ψt) atol = 1.0e-8 + end + + hs = SpinSpace(:s) + for s in (1 // 2, 1) + d = Int(2s + 1) + bq = QOB.SpinBasis(s) + ψq = QOB.spinup(bq) + ψt = QTB.basis(d, 0) + for ax in 1:3 + @test numeric_average(Spin(hs, :S, ax), ψq) ≈ numeric_average(Spin(hs, :S, ax), ψt) atol = 1.0e-8 + end + end + end + + @testset "Composite parity" begin + h = FockSpace(:c) ⊗ NLevelSpace(:atom, 3, 1) + a = Destroy(h, :a, 1) + σ(i, j) = Transition(h, :σ, i, j, 2) + ba = QOB.NLevelBasis(3) + ψq = QOB.coherentstate(QOB.FockBasis(4), 0.3) ⊗ + ((QOB.nlevelstate(ba, 1) + QOB.nlevelstate(ba, 2)) / sqrt(2)) + ψt = QTB.tensor(QTB.coherent(5, 0.3), (QTB.basis(3, 0) + QTB.basis(3, 1)) / sqrt(2)) + for op in (a' * a, a' * σ(1, 2) + a * σ(2, 1), a' * a * σ(2, 2)) + @test numeric_average(op, ψq) ≈ numeric_average(op, ψt) atol = 1.0e-6 + end + end + + @testset "Time-dependent (native QobjEvo)" begin + h = FockSpace(:f) + @qnumbers a::Destroy(h) + @variables E::Number + N = 6 + f = to_numeric(E * a + conj(E) * a', h, N; backend = QTBB, time_parameter = Dict(E => t -> 1 + im * t)) + @test f isa QTB.QuantumObjectEvolution + ψt = QTB.coherent(N + 1, 0.3 - 0.1im) + aT = QTB.destroy(N + 1) + @test QTB.expect(f(nothing, 0.5), ψt) ≈ QTB.expect((1 + 0.5im) * aT + (1 - 0.5im) * aT', ψt) atol = 1.0e-8 + + @test to_numeric( + E * a, h, N; + backend = QTBB, time_parameter = Dict(E => t -> 1), op_type = identity, + ) isa QTB.QuantumObjectEvolution + @test_throws ArgumentError to_numeric( + E * a, h, N; + backend = QTBB, time_parameter = Dict(E => t -> 1), op_type = QTB.to_dense, + ) + end + + @testset "sesolve / mesolve consume the lazy VecSum" begin + h = FockSpace(:f) + @qnumbers a::Destroy(h) + N = 6 + H = to_numeric( + 2.0 * a' * a + 0.5 * (a + a'), h, N; + backend = QTBB, op_type = identity, + ) + Hc = QTB.Qobj(SO.concretize(H.data); dims = N + 1) # eager reference + ψ = QTB.coherent(N + 1, 0.3 + 0.1im) + tl = collect(range(0, 1.0, length = 11)) + eop = QTB.create(N + 1) * QTB.destroy(N + 1) + + rs = QTB.sesolve(H, ψ, tl; e_ops = [eop], progress_bar = Val(false)) + rsc = QTB.sesolve(Hc, ψ, tl; e_ops = [eop], progress_bar = Val(false)) + @test maximum(abs, rs.expect[1, :] .- rsc.expect[1, :]) < 1.0e-6 + + c_ops = [QTB.destroy(N + 1)] + rm = QTB.mesolve(H, ψ, tl, c_ops; e_ops = [eop], progress_bar = Val(false)) + rmc = QTB.mesolve(Hc, ψ, tl, c_ops; e_ops = [eop], progress_bar = Val(false)) + @test maximum(abs, rm.expect[1, :] .- rmc.expect[1, :]) < 1.0e-6 + @test abs(tr(rm.states[end]) - 1) < 1.0e-6 + end + + @testset "numeric_average: averages, scalars, vector of states" begin + h = FockSpace(:f) + @qnumbers a::Destroy(h) + α = 0.4 + 0.2im + ψt = QTB.coherent(8, α) + @test numeric_average(average(a), ψt) ≈ α atol = 1.0e-4 + @test numeric_average(average(a) + average(a'), ψt) ≈ α + conj(α) atol = 1.0e-4 + @test numeric_average(3, ψt) === ComplexF64(3) + + ψs = [QTB.coherent(8, β) for β in (0.1 + 0.2im, -0.2 + 0.1im)] + @test numeric_average(a, ψs) ≈ [QTB.expect(QTB.destroy(8), ψ) for ψ in ψs] atol = 1.0e-8 + end + + @testset "materialize on demand" begin + h = FockSpace(:f) + @qnumbers a::Destroy(h) + H = to_numeric(a' * a, h, 4; backend = QTBB, op_type = identity) + M = SO.concretize(H.data) + @test M ≈ QTB.create(5).data * QTB.destroy(5).data + end + + @testset "empty QAdd assembles to zero" begin + # A `QAdd` that cancels to zero terms still assembles to a valid zero operator on + # the full space (the empty-terms branch of `numeric_assemble`). + h = FockSpace(:f) + @qnumbers a::Destroy(h) + z = 1 * a - 1 * a + M = to_numeric(z, h, 5; backend = QTBB, op_type = identity) + @test SO.concretize(M.data) ≈ zeros(ComplexF64, 6, 6) + end + + @testset "time-dependent with constant term" begin + # A TD assembly mixing a constant term (`Δ*a'a`) with a genuinely time-dependent + # term (`f*a`) exercises both `_td_scalar` branches. + h = FockSpace(:f) + @qnumbers a::Destroy(h) + @variables Δ f + H = to_numeric( + Δ * a' * a + f * a, h, 5; + backend = QTBB, parameter = Dict(Δ => 1.0), time_parameter = Dict(f => t -> 2.0 + 0im), + ) + M = H(nothing, 0.0).data + Href = QTB.create(6).data * QTB.destroy(6).data + 2.0 * QTB.destroy(6).data + @test M ≈ Href + end + + @testset "unsupported operator role errors" begin + # A collective transition has no QuantumToolbox realisation: the closed operator + # ladder falls through to the throwing `Val` extension point. + hc = CollectiveNLevelSpace(:atom, 3) + S12 = CollectiveTransition(hc, :S, 1, 2) + @test_throws ArgumentError to_numeric(S12, 3) + end + + @testset "default backend errors when both are loaded" begin + # With QuantumOpticsBase and QuantumToolbox both loaded, the `HilbertSpace` form + # cannot pick a backend and must ask for an explicit `backend =`. + h = FockSpace(:f) + @qnumbers a::Destroy(h) + @test_throws ArgumentError to_numeric(a' * a, h, 5) + end + + @testset "Indexed path works for QuantumToolbox" begin + # Sum of excited-state projectors over two atoms is the excited-atom-count operator. + ha = NLevelSpace(:atom, 2) + i = Index(ha, :i, 2, ha) + σ(a, b, k) = IndexedOperator(Transition(ha, :σ, a, b), k) + H = Σ(σ(2, 2, i), i) + sites = Dict{Int, Vector{Int}}(1 => [1, 2]) + + Mt = to_numeric(H, [2, 2], sites) # QTB (Vector{Int} dims) + @test Mt isa QTB.QuantumObject + + P = QTB.projection(2, 1, 1) # |1><1|, 0-indexed excited level + ref = QTB.tensor(P, QTB.qeye(2)) + QTB.tensor(QTB.qeye(2), P) + @test Mt.data ≈ ref.data + + # Cross-check the eager expectation via a scalar (matrices differ by kron convention + # across backends, so compare scalars, not raw matrices). + ψt = QTB.tensor(QTB.basis(2, 1), QTB.basis(2, 1)) # both excited -> eigenvalue 2 + @test real(QTB.expect(Mt, ψt)) ≈ 2 atol = 1.0e-10 + + @test_throws ArgumentError to_numeric(H, 2, sites) + @test_throws ArgumentError to_numeric(H, [2, 2], Dict(1 => [1, 1])) + end + + @testset "Fock/PhaseSpace dims match QuantumOptics" begin + h = FockSpace(:f) + @qnumbers a::Destroy(h) + # cutoff 5 gives dimension 6 (occupations 0:5) on BOTH backends. + @test SecondQuantizedAlgebra.numeric_basis(QuantumToolboxBackend(), FockSpace(:f), 5) == 6 + @test SecondQuantizedAlgebra.numeric_basis(QuantumToolboxBackend(), PhaseSpace(:x), 5) == 6 + + Ht = to_numeric(a' * a, h, 5; backend = QuantumToolboxBackend()) + @test size(Ht.data) == (6, 6) + # highest occupation is 5 (would be 4 under the old off-by-one convention). + @test maximum(real, diag(Matrix(Ht.data))) ≈ 5 + + # Direct QTB dimensions are raw matrix dimensions; the HilbertSpace form receives + # a symbolic Fock cutoff. + @test size(to_numeric(a, 5).data) == (5, 5) + @test size(to_numeric(a, h, 5; backend = QTBB).data) == (6, 6) + + # Bare QAdd on raw QTB dimensions (positional QAdd-on-dims form). + qadd_raw = to_numeric(a' + a, 5) + @test size(qadd_raw.data) == (5, 5) + @test Matrix(qadd_raw.data) ≈ Matrix((to_numeric(a', 5) + to_numeric(a, 5)).data) + end + + @testset "explicit op_type materializes eager" begin + h = FockSpace(:f) + @qnumbers a::Destroy(h) + ref = QTB.create(6).data * QTB.destroy(6).data # cutoff 5 gives dim 6 + + sp = to_numeric(a' * a, h, 5; backend = QTBB, op_type = QTB.to_sparse) + @test sp isa QTB.QuantumObject + @test sp.data ≈ ref + + de = to_numeric(a' * a, h, 5; backend = QTBB, op_type = QTB.to_dense) + @test de isa QTB.QuantumObject + @test de.data ≈ ref + + co = to_numeric(a' * a, h, 5; backend = QTBB, op_type = SO.concretize) + @test co isa AbstractMatrix + @test co ≈ ref + + # identity keeps the lazy assembly. + lz = to_numeric(a' * a, h, 5; backend = QTBB, op_type = identity) + @test lz isa QTB.QuantumObjectEvolution + end + + @testset "vector API for QuantumToolbox" begin + h = FockSpace(:f) + @qnumbers a::Destroy(h) + + # HilbertSpace form (backend-neutral entry point). + vhs = to_numeric([a, a'], h, 5; backend = QTBB) + @test vhs isa AbstractVector + @test length(vhs) == 2 + @test all(x -> x isa QTB.QuantumObject, vhs) + + # Direct QTB-dims form. + vd = to_numeric([a, a'], 5) + @test vd isa AbstractVector + @test length(vd) == 2 + @test all(x -> x isa QTB.QuantumObject, vd) + end + + @testset "public backend hooks and product dims validation" begin + ψ = QTB.basis(3, 0) + @test numeric_backend(ψ) isa QuantumToolboxBackend + @test numeric_basis(ψ) == 3 + @test SecondQuantizedAlgebra.numeric_num_subsystems(QTBB, 3) == 1 + @test SecondQuantizedAlgebra.numeric_num_subsystems(QTBB, [2, 3]) == 2 + + h = FockSpace(:a) ⊗ FockSpace(:b) + a = Destroy(h, :a, 1) + @test_throws ArgumentError to_numeric(a, h, [2]; backend = QTBB) + @test_throws ArgumentError to_numeric(a, h, [2, 3, 99]; backend = QTBB) + # dims that are not an indexable collection are rejected up front. + @test_throws ArgumentError to_numeric(a, h, nothing; backend = QTBB) + + # A valid ProductSpace HilbertSpace conversion builds per-subspace dims (Fock + # cutoffs 2 and 3 give dimensions 3 and 4, so the composite is 12×12). + Mp = to_numeric(a, h, [2, 3]; backend = QTBB) + @test size(Mp.data) == (12, 12) + end +end diff --git a/test/numeric_test.jl b/test/numeric_test.jl index bb86c68f..c6bea351 100644 --- a/test/numeric_test.jl +++ b/test/numeric_test.jl @@ -1,12 +1,60 @@ using SecondQuantizedAlgebra -import SecondQuantizedAlgebra: QAdd, QSym, _single_qadd, _to_cnum, _to_complex, _fold_const +import SecondQuantizedAlgebra: QAdd, QSym, _single_qadd, _to_cnum, _to_complex, _fold_const, + _to_numeric_static, NumericContext using QuantumOpticsBase using Symbolics: @variables, substitute import SymbolicUtils using Test +# `to_numeric` materialises via `op_type` (default `sparse`), so the return type depends only +# on `op_type`, never on the shape of the expression: a bare operator, a product, and a sum +# all return a concrete `Operator`. `op_type = identity` opts into the natural lazy form +# (`LazyTensor`/`LazyProduct`/`LazySum`). The time-dependent form returns a native +# `TimeDependentSum`, evaluated via `H(t)` and compared through `expect`. +_dat(x) = dense(x).data + +# Tiny backend used as a conformance test for the documented third-party static interface. +struct MockNumericBackend <: SecondQuantizedAlgebra.NumericBackend end +struct MockEagerOperator + data::Matrix{ComplexF64} +end +_mock_identity(n) = ComplexF64[i == j for i in 1:n, j in 1:n] +SecondQuantizedAlgebra.numeric_basis(::MockNumericBackend, ::FockSpace, cutoff) = Int(cutoff) + 1 +SecondQuantizedAlgebra.numeric_num_subsystems(::MockNumericBackend, ::Int) = 1 +SecondQuantizedAlgebra.numeric_subbasis(::MockNumericBackend, n::Int, slot::Int) = + slot == 1 ? n : throw(ArgumentError("mock basis has no subsystem slot $slot")) +SecondQuantizedAlgebra.numeric_operator( + ::MockNumericBackend, ::SecondQuantizedAlgebra.Op, n::Int, +) = _mock_identity(n) +SecondQuantizedAlgebra.numeric_embed(::MockNumericBackend, ::Int, ::Int, leaf) = leaf +SecondQuantizedAlgebra.numeric_identity(::MockNumericBackend, n::Int) = _mock_identity(n) +function SecondQuantizedAlgebra.numeric_assemble(::MockNumericBackend, n::Int, terms) + result = zeros(ComplexF64, n, n) + for (coefficient, factors) in terms + term = isempty(factors) ? _mock_identity(n) : foldl(*, factors) + result .+= coefficient .* term + end + return result +end +SecondQuantizedAlgebra.numeric_materialize(::MockNumericBackend, assembled, ::Nothing) = + MockEagerOperator(assembled) +SecondQuantizedAlgebra.numeric_materialize( + ::MockNumericBackend, assembled, ::typeof(identity), +) = assembled + @testset "numeric conversion" begin - @testset "Single space — basic" begin + @testset "third-party static backend contract" begin + h = FockSpace(:mock) + @qnumbers a::Destroy(h) + + eager = to_numeric(2 * a, h, 3; backend = MockNumericBackend()) + lazy = to_numeric(2 * a, h, 3; backend = MockNumericBackend(), op_type = identity) + @test eager isa MockEagerOperator + @test eager.data == 2 .* _mock_identity(4) + @test lazy == eager.data + end + + @testset "Single space basic" begin h = FockSpace(:fock) @qnumbers a::Destroy(h) b = FockBasis(7) @@ -20,8 +68,10 @@ using Test @qnumbers a::Destroy(h) b = FockBasis(7) - @test to_numeric(a' * a, b) == create(b) * destroy(b) - @test to_numeric(2 * a, b) == 2 * destroy(b) + @test to_numeric(a' * a, b) isa Operator + @test to_numeric(a' * a, b; op_type = identity) isa LazySum + @test _dat(to_numeric(a' * a, b)) == _dat(create(b) * destroy(b)) + @test _dat(to_numeric(2 * a, b)) == _dat(2 * destroy(b)) end @testset "QAdd" begin @@ -30,8 +80,9 @@ using Test b = FockBasis(7) result = to_numeric(a + a', b) - expected = destroy(b) + create(b) - @test result == expected + @test result isa Operator + @test to_numeric(a + a', b; op_type = identity) isa LazySum + @test _dat(result) == _dat(destroy(b) + create(b)) end @testset "Scalar" begin @@ -44,12 +95,10 @@ using Test @qnumbers a1::Destroy(h, 1) a2::Destroy(h, 2) b = FockBasis(3) ⊗ FockBasis(3) - a1_lazy = to_numeric(a1, b; op_type = identity) - @test a1_lazy isa LazyTensor a1_num = to_numeric(a1, b) @test a1_num isa Operator - @test a1_num == sparse(a1_lazy) - @test to_numeric(a1, b; op_type = dense) == dense(a1_lazy) + # The lazy form is the vector-backed LazySum for every shape (single op included). + @test to_numeric(a1, b; op_type = identity) isa LazySum end @testset "numeric_average" begin @@ -100,8 +149,8 @@ using Test bf = FockBasis(3) bn = NLevelBasis(3) bc = bf ⊗ bn - @test to_numeric(σ12, bc; op_type = identity) isa LazyTensor @test to_numeric(σ12, bc) isa Operator + @test to_numeric(σ12, bc; op_type = identity) isa LazySum end @testset "Type stability" begin @@ -110,11 +159,20 @@ using Test b = FockBasis(7) ψ = fockstate(b, 2) - # to_numeric: leaf + QAdd + dict-substitution paths - @test @inferred(to_numeric(a, b)) isa AbstractOperator - @test @inferred(to_numeric(a' * a, b)) isa AbstractOperator - @test @inferred(to_numeric(2 * a + 3 * a' + 5, b)) isa AbstractOperator - @test @inferred(to_numeric(a, b, Dict{QSym, Union{}}())) isa AbstractOperator + # The lazy assembly is inference-stable and one concrete type for any term count: the + # property the 5-arg vector-backed `LazySum` constructor buys. `sparse` materialization + # (the default) is a top-level convenience and NOT on the `@inferred` contract, since + # `sparse` of the abstract-eltype `LazySum` widens. + ctx = NumericContext(QuantumOpticsBackend(), b, Dict{QSym, Union{}}()) + @test @inferred(_to_numeric_static(a' * a, ctx)) isa AbstractOperator + @test @inferred(_to_numeric_static(2 * a + 3 * a' + 5, ctx)) isa AbstractOperator + @test typeof(_to_numeric_static(a + a', ctx)) === typeof(_to_numeric_static(a + a' + a' * a, ctx)) + + # to_numeric: leaf + QAdd + dict-substitution paths materialise a concrete operator. + @test to_numeric(a, b) isa AbstractOperator + @test to_numeric(a' * a, b) isa Operator + @test to_numeric(2 * a + 3 * a' + 5, b) isa Operator + @test to_numeric(a, b, Dict{QSym, Union{}}()) isa AbstractOperator # numeric_average: every BasicSymbolic branch infers to ComplexF64 @test @inferred(numeric_average(a' * a, ψ)) isa ComplexF64 @@ -125,9 +183,24 @@ using Test @test @inferred(numeric_average(3, ψ)) isa ComplexF64 @test @inferred(numeric_average(3.5 + 1im, ψ)) isa ComplexF64 - # `_to_complex(::Any) → ComplexF64` is the keystone; canary against a + # `_to_complex(::Any) -> ComplexF64` is the keystone; canary against a # 7th overload tripping Julia's union-split budget. @test Base.return_types(_to_complex, (Any,))[1] === ComplexF64 + @test length(methods(_to_complex)) == 2 + end + + @testset "op_type materialization" begin + h = FockSpace(:fock) + @qnumbers a::Destroy(h) + b = FockBasis(7) + # The default materialises `sparse`; `op_type = identity` keeps the lazy assembly. + @test to_numeric(a' * a, b) isa Operator + @test to_numeric(2 * a + 3 * a', b) isa Operator + @test to_numeric(a, b) isa Operator + @test to_numeric(a' * a, b; op_type = identity) isa LazySum + @test to_numeric(2 * a + 3 * a', b; op_type = identity) isa LazySum + @test to_numeric(a, b) == sparse(to_numeric(a, b; op_type = identity)) + @test _dat(to_numeric(a' * a, b; op_type = dense)) == _dat(to_numeric(a' * a, b)) end @testset "numeric_average: Average expressions" begin @@ -137,20 +210,17 @@ using Test α = 0.1 + 0.2im ψ = coherentstate(b, α) - # numeric_average on average(a) should give ⟨a⟩ = α avg_a = average(a) @test numeric_average(avg_a, ψ) ≈ α - # Sum of averages avg_sum = average(a) + average(a') @test numeric_average(avg_sum, ψ) ≈ α + conj(α) - # Scalar times average avg_scaled = 2 * average(a) @test numeric_average(avg_scaled, ψ) ≈ 2α end - @testset "to_numeric — Dict substitution" begin + @testset "to_numeric: Dict substitution" begin h = FockSpace(:fock) @qnumbers a::Destroy(h) b = FockBasis(7) @@ -160,30 +230,27 @@ using Test # Dict substitution replaces operator @test to_numeric(a, b, d) == custom_op - # Adjoint not in dict → falls back to normal + # Adjoint not in dict, falls back to normal @test to_numeric(a', b, d) == create(b) # QAdd with Dict result_mul = to_numeric(a' * a, b, d) - expected_mul = create(b) * custom_op - @test result_mul == expected_mul + @test _dat(result_mul) == _dat(create(b) * custom_op) end - @testset "numeric_average — Dict substitution" begin + @testset "numeric_average: Dict substitution" begin h = FockSpace(:fock) @qnumbers a::Destroy(h) b = FockBasis(7) α = 0.1 + 0.2im ψ = coherentstate(b, α) - d = Dict{QSym, Any}() # empty dict — same as no dict + d = Dict{QSym, Any}() # empty dict, same as no dict @test numeric_average(a, ψ, d) ≈ α - # Average expression with dict avg_a = average(a) @test numeric_average(avg_a, ψ, d) ≈ α - # Number passthrough — coerced to ComplexF64 like every other branch @test numeric_average(3, ψ, d) === ComplexF64(3) end @@ -202,19 +269,16 @@ using Test i == j == 1 && continue op1 = a * σprod_gap(i, j) op2 = a' * σprod_gap(i, j) - @test to_numeric(op1, bprod_gap; op_type = identity) == LazyTensor( - bprod_gap, - [1, 3], + ref1 = LazyTensor( + bprod_gap, [1, 3], (destroy(bfock), QuantumOpticsBase.transition(bnlevel, i, j)), ) - @test to_numeric(op2, bprod_gap; op_type = identity) == LazyTensor( - bprod_gap, - [1, 3], + ref2 = LazyTensor( + bprod_gap, [1, 3], (create(bfock), QuantumOpticsBase.transition(bnlevel, i, j)), ) - @test to_numeric(op1, bprod_gap) == sparse( - to_numeric(op1, bprod_gap; op_type = identity), - ) + @test _dat(to_numeric(op1, bprod_gap)) == _dat(ref1) + @test _dat(to_numeric(op2, bprod_gap)) == _dat(ref2) end end @@ -223,18 +287,14 @@ using Test @qnumbers a::Destroy(hfock) bfock = FockBasis(100) - @test isequal( - 2 * create(bfock) + 2 * destroy(bfock), - to_numeric(2 * a + 2 * a', bfock), - ) - @test iszero( - (2 * create(bfock) + 2 * destroy(bfock)) - - to_numeric(2 * a + 2 * a', bfock), - ) - @test isequal(to_numeric(2 * a, bfock), 2 * to_numeric(a, bfock)) + ref = 2 * create(bfock) + 2 * destroy(bfock) + got = to_numeric(2 * a + 2 * a', bfock) + @test isequal(_dat(ref), _dat(got)) + @test iszero(_dat(ref) - _dat(got)) + @test isequal(_dat(to_numeric(2 * a, bfock)), _dat(2 * destroy(bfock))) end - @testset "numeric_average — product state" begin + @testset "numeric_average: product state" begin hfock = FockSpace(:fock) hnlevel = NLevelSpace(:nlevel, 3, 1) hprod = hfock ⊗ hnlevel @@ -256,33 +316,7 @@ using Test end end - @testset "numeric_average — LazyKet state" begin - # `expect` has no method for a sparse operator paired with a `LazyKet`. - hfock = FockSpace(:fock) - hnlevel = NLevelSpace(:nlevel, (:a, :b, :c)) - hprod = hfock ⊗ hnlevel - bfock = FockBasis(10) - bnlevel = NLevelBasis(3) - bprod = bfock ⊗ bnlevel - - @qnumbers a::Destroy(hprod, 1) - σ(i, j) = Transition(hprod, :σ, i, j, 2) - - α = 0.3 + 0.0im - ket_fock = coherentstate(bfock, α) - ket_nlevel = (nlevelstate(bnlevel, 1) + nlevelstate(bnlevel, 3)) / sqrt(2) - ψ_lazy = LazyKet(bprod, (ket_fock, ket_nlevel)) - ψ_dense = ket_fock ⊗ ket_nlevel - - for op in (a, a' * σ(:a, :c), a + a' * σ(:a, :c)) - op_num = to_numeric(op, bprod) - @test numeric_average(op, ψ_lazy) ≈ expect(op_num, ψ_dense) - end - @test numeric_average(average(a' * a), ψ_lazy) ≈ - expect(to_numeric(a' * a, bprod), ψ_dense) - end - - @testset "numeric_average — comprehensive expressions" begin + @testset "numeric_average: comprehensive expressions" begin h = FockSpace(:fock) @qnumbers a::Destroy(h) b = FockBasis(7) @@ -297,42 +331,39 @@ using Test @test numeric_average(3, ψ) ≈ 3 end - @testset "numeric_average — vector of states" begin + @testset "numeric_average: vector of states" begin h = FockSpace(:fock) @qnumbers a::Destroy(h) b = FockBasis(7) αs = (0.1 + 0.2im, -0.3 + 0.4im, 0.5 + 0.0im) ψs = [coherentstate(b, α) for α in αs] - # Vector input → vector output, same as broadcasting the scalar form @test numeric_average(a, ψs) ≈ [α for α in αs] @test numeric_average(a' * a, ψs) ≈ [abs2(α) for α in αs] @test numeric_average(average(a), ψs) ≈ [α for α in αs] - # expect alias matches @test expect(a, ψs) ≈ numeric_average(a, ψs) @test expect(a' * a, ψs[1]) ≈ numeric_average(a' * a, ψs[1]) @test numeric_average(average(a), ψs[1]) ≈ αs[1] - # Empty vector errors empty_ψs = typeof(ψs[1])[] @test_throws ArgumentError numeric_average(a, empty_ψs) end - @testset "numeric_average — vector of states with Dict" begin + @testset "numeric_average: vector of states with Dict" begin h = FockSpace(:fock) @qnumbers a::Destroy(h) b = FockBasis(7) αs = (0.1 + 0.2im, -0.3 + 0.4im) ψs = [coherentstate(b, α) for α in αs] - d = Dict{QSym, Any}() # empty dict — passthrough to base case + d = Dict{QSym, Any}() @test numeric_average(a, ψs, d) ≈ numeric_average(a, ψs) @test numeric_average(average(a' * a), ψs, d) ≈ [abs2(α) for α in αs] @test expect(a, ψs, d) ≈ numeric_average(a, ψs, d) end - @testset "numeric_average — Dict comprehensive" begin + @testset "numeric_average: Dict comprehensive" begin nQDs = 2 h_qc1 = FockSpace(:ada) h_qc2 = FockSpace(:n) @@ -349,10 +380,10 @@ using Test dd = Dict([ad, a] .=> [s(2, 2, 1), s(2, 1, 2)]) @test to_numeric(a, b_test, Dict{QSym, Any}()) == to_numeric(a, b_test) - @test to_numeric(ad * n, b_test, Dict{QSym, Any}()) == to_numeric(ad * n, b_test) - @test to_numeric(2 * ad * a, b_test, Dict{QSym, Any}()) == to_numeric(2 * ad * a, b_test) + @test _dat(to_numeric(ad * n, b_test, Dict{QSym, Any}())) == _dat(to_numeric(ad * n, b_test)) + @test _dat(to_numeric(2 * ad * a, b_test, Dict{QSym, Any}())) == _dat(to_numeric(2 * ad * a, b_test)) @test to_numeric(ad, b_all, dd) == s(2, 2, 1) - @test to_numeric(2 * ad * a, b_all, dd) == 2 * s(2, 2, 1) * s(2, 1, 2) + @test _dat(to_numeric(2 * ad * a, b_all, dd)) == _dat(2 * s(2, 2, 1) * s(2, 1, 2)) @test dense(to_numeric(3, b_all, dd)) == one(b_all) * 3 ψ0 = tensor([nlevelstate(bs, 2) for i in 1:nQDs]...) @@ -366,7 +397,7 @@ using Test expect(s(2, 2, 1) * s(2, 1, 2), ψ0) + expect(s(2, 1, 2), ψ0) end - @testset "Allocations — to_numeric" begin + @testset "Allocations: to_numeric" begin h = FockSpace(:fock) @qnumbers a::Destroy(h) b = FockBasis(7) @@ -376,7 +407,6 @@ using Test to_numeric(a', b) to_numeric(a' * a, b) - # Single operators should be bounded @test @allocations(to_numeric(a, b)) < 50 @test @allocations(to_numeric(a', b)) < 50 @test @allocations(to_numeric(a' * a, b)) < 1500 @@ -389,19 +419,18 @@ using Test @test to_numeric(a, b) == destroy(b) @test to_numeric(a', b) == create(b) - @test to_numeric(a' * a, b) == create(b) * destroy(b) + @test _dat(to_numeric(a' * a, b)) == _dat(create(b) * destroy(b)) α = 0.1 + 0.2im ψ = coherentstate(b, α) @test numeric_average(a, ψ) ≈ α @test numeric_average(a' * a, ψ) ≈ abs2(α) - # numeric_average on QAdd expr = a + a' * a @test numeric_average(expr, ψ) ≈ α + abs2(α) - # to_numeric with scalar QAdd (empty operators) - @test to_numeric(_single_qadd(_to_cnum(3), Op[]), b) == 3 * one(b) + # to_numeric with scalar QAdd (empty operators), now a lazy identity + @test dense(to_numeric(_single_qadd(_to_cnum(3), Op[]), b)) == 3 * one(b) end @testset "Number-symtype coefficient round-trip" begin @@ -409,14 +438,12 @@ using Test b = FockBasis(3) a = Destroy(FockSpace(:c), :a) - # g†g must reduce to |g|², i.e. conj(g)*g, after substituting a complex value. for v in (2 + 3im, 1 + 1im, 0.5 - 0.25im) op = substitute((g * a)' * (g * a), Dict(g => v)) - @test to_numeric(op, b) ≈ abs2(v) * to_numeric(a' * a, b) + @test _dat(to_numeric(op, b)) ≈ abs2(v) * _dat(to_numeric(a' * a, b)) end - # A bare complex coupling carries through (no conj): g·a → value·a. - @test to_numeric(substitute(g * a, Dict(g => 2 + 3im)), b) ≈ (2 + 3im) * to_numeric(a, b) + @test _dat(to_numeric(substitute(g * a, Dict(g => 2 + 3im)), b)) ≈ (2 + 3im) * _dat(to_numeric(a, b)) end @testset "Public coefficient lowering" begin @@ -436,21 +463,24 @@ using Test b = FockBasis(4) A = destroy(b) Ad = create(b) + ψ = coherentstate(b, 0.3 - 0.1im) @variables x::Real ϕ::Real E::Number - @test to_numeric(substitute(sqrt(x) * a, Dict(x => 2.0)), b) ≈ sqrt(2.0) * A - @test to_numeric(substitute(exp(im * ϕ) * a, Dict(ϕ => 0.5)), b) ≈ exp(0.5im) * A + @test _dat(to_numeric(substitute(sqrt(x) * a, Dict(x => 2.0)), b)) ≈ sqrt(2.0) * _dat(A) + @test _dat(to_numeric(substitute(exp(im * ϕ) * a, Dict(ϕ => 0.5)), b)) ≈ exp(0.5im) * _dat(A) - @test to_numeric(sqrt(x) * a, b; parameter = Dict(x => 2.0)) ≈ sqrt(2.0) * A - @test to_numeric(exp(im * ϕ) * a, b; parameter = Dict(ϕ => 0.5)) ≈ exp(0.5im) * A + @test _dat(to_numeric(sqrt(x) * a, b; parameter = Dict(x => 2.0))) ≈ sqrt(2.0) * _dat(A) + @test _dat(to_numeric(exp(im * ϕ) * a, b; parameter = Dict(ϕ => 0.5))) ≈ exp(0.5im) * _dat(A) custom = 2 * A - @test to_numeric(a', b; operators = Dict(a => custom)) == custom' - @test to_numeric(a', b; operators = Dict(a => custom), adjoint_ops = false) == Ad - @test to_numeric([a, a'], b; operators = Dict(a => custom)) == [custom, custom'] + @test _dat(to_numeric(a', b; operators = Dict(a => custom))) == _dat(custom') + @test _dat(to_numeric(a', b; operators = Dict(a => custom), adjoint_ops = false)) == _dat(Ad) + @test [_dat(x) for x in to_numeric([a, a'], b; operators = Dict(a => custom))] == [_dat(custom), _dat(custom')] + # Time-dependent: native TimeDependentSum, evaluated at a time, compared via expect. f = to_numeric(E * a + conj(E) * a', b; time_parameter = Dict(E => t -> 1 + im * t)) - @test f(0.5) ≈ (1 + 0.5im) * A + (1 - 0.5im) * Ad + @test f isa TimeDependentSum + @test expect(f(0.5), ψ) ≈ expect((1 + 0.5im) * A + (1 - 0.5im) * Ad, ψ) end @testset "keyword to_numeric on composite basis" begin @@ -461,24 +491,20 @@ using Test Ia = identityoperator(b) an = to_numeric(a, b) - # Regression: a scalar/constant term must build the full-system identity, not - # the identity of a single subsystem (the `_lazy_one(b)` fix). Before, the - # `+ 5` term produced a wrongly-sized identity on a composite basis. + # Regression: a scalar/constant term must build the full-system identity, not the + # identity of a single subsystem (the `numeric_identity(composite)` fix). H = to_numeric(Δ * a' * a + 5, b; parameter = Dict(Δ => 2.0)) @test size(dense(H).data) == (length(b), length(b)) @test dense(H).data ≈ dense(2.0 * an' * an + 5 * Ia).data - # Custom operator on one subsystem; the constant term is still full-sized. custom = 2 * an H2 = to_numeric(a + 3, b; operators = Dict(a => custom)) @test dense(H2).data ≈ dense(custom + 3 * Ia).data - # A term that cancels to zero still yields a full-sized zero operator. z = to_numeric(a - a, b) @test size(dense(z).data) == (length(b), length(b)) @test dense(z).data ≈ zeros(length(b), length(b)) - # Missing adjoint rule is auto-added for custom operators. r = to_numeric(a', b; operators = Dict(a => custom)) @test dense(r).data ≈ dense(custom').data r2 = to_numeric(a', b; operators = Dict(a => custom), adjoint_ops = false) @@ -490,46 +516,60 @@ using Test @qnumbers a::Destroy(h) b = FockBasis(4) A = destroy(b) + ψ = coherentstate(b, 0.3 - 0.1im) @variables E::Number - # A plain number value becomes a constant-in-time closure. + # A plain number value becomes a constant-in-time term. f0 = to_numeric(E * a, b; time_parameter = Dict(E => 3.0)) - @test f0(0.0) ≈ 3.0 * A - @test f0(10.0) ≈ 3.0 * A + @test expect(f0(0.0), ψ) ≈ expect(3.0 * A, ψ) + @test expect(f0(10.0), ψ) ≈ expect(3.0 * A, ψ) - # A constant-coefficient term still returns a callable when time_parameter is set. + # A constant-coefficient term still returns a native TD operator. f1 = to_numeric(2.0 * a, b; time_parameter = Dict(E => t -> 1.0 + 0im)) - @test f1(7.0) ≈ 2.0 * A + @test expect(f1(7.0), ψ) ≈ expect(2.0 * A, ψ) # `conj(v)` is accepted as a time_parameter key. f2 = to_numeric(conj(E) * a, b; time_parameter = Dict(conj(E) => t -> 2 + im * t)) - @test f2(1.0) ≈ (2 + 1im) * A + @test expect(f2(1.0), ψ) ≈ expect((2 + 1im) * A, ψ) + + @test to_numeric(E * a, b; time_parameter = Dict(E => t -> 1), op_type = identity) isa TimeDependentSum + @test_throws ArgumentError to_numeric(E * a, b; time_parameter = Dict(E => t -> 1), op_type = dense) + end + + @testset "public backend hooks and product dims validation" begin + b = FockBasis(3) + ψ = fockstate(b, 0) + @test numeric_backend(ψ) isa QuantumOpticsBackend + @test numeric_basis(ψ) == b + @test SecondQuantizedAlgebra.numeric_num_subsystems(QuantumOpticsBackend(), b) == 1 + + h = FockSpace(:a) ⊗ FockSpace(:b) + a = Destroy(h, :a, 1) + @test_throws ArgumentError to_numeric(a, h, [2]; backend = QuantumOpticsBackend()) + @test_throws ArgumentError to_numeric(a, h, [2, 3, 99]; backend = QuantumOpticsBackend()) end @testset "keyword to_numeric, scalar argument" begin b = FockBasis(4) Ib = one(b) + ψ = coherentstate(b, 0.3 - 0.1im) @variables x::Real E::Number - # Plain number scalar routes through the keyword path to a scaled identity. @test to_numeric(3, b; parameter = Dict(x => 1.0)) == 3 * Ib @test to_numeric(2.0 + 0im, b; parameter = Dict{Any, Any}()) == (2.0 + 0im) * Ib - # A symbolic scalar resolved by `parameter` becomes a constant identity. - @test to_numeric(2 * x, b; parameter = Dict(x => 1.5)) ≈ 3.0 * Ib + @test dense(to_numeric(2 * x, b; parameter = Dict(x => 1.5))) ≈ 3.0 * Ib - # A symbolic scalar with no value cannot be translated. @test_throws ArgumentError to_numeric(x, b) - # With `time_parameter`, a constant scalar yields a constant-in-time closure. + # With `time_parameter`, a constant scalar yields a constant-in-time TD operator. fconst = to_numeric(2.0, b; time_parameter = Dict(E => t -> 1.0 + 0im)) - @test fconst(0.0) == 2.0 * Ib - @test fconst(9.0) == 2.0 * Ib + @test expect(fconst(0.0), ψ) ≈ expect(2.0 * Ib, ψ) + @test expect(fconst(9.0), ψ) ≈ expect(2.0 * Ib, ψ) - # A genuinely time-dependent scalar yields a time-varying closure. ft = to_numeric(E, b; time_parameter = Dict(E => t -> 1 + im * t)) - @test ft(0.0) ≈ (1 + 0im) * Ib - @test ft(2.0) ≈ (1 + 2im) * Ib + @test expect(ft(0.0), ψ) ≈ expect((1 + 0im) * Ib, ψ) + @test expect(ft(2.0), ψ) ≈ expect((1 + 2im) * Ib, ψ) end @testset "keyword to_numeric, state argument" begin @@ -540,9 +580,8 @@ using Test ψ = coherentstate(b, α) @variables Δ::Real - # The state form derives the basis from the state and forwards keywords. op_state = to_numeric(Δ * a, ψ; parameter = Dict(Δ => 2.0)) - @test op_state == 2.0 * destroy(b) + @test _dat(op_state) == _dat(2.0 * destroy(b)) end @testset "constant symbolic coefficient reduction" begin @@ -552,10 +591,6 @@ using Test A = destroy(b) @variables z::Number - # Coefficients that stay symbolic-but-constant (built from `real`/`imag`/`conj`, - # which are not folded into the native coefficient tier) are reduced to a - # concrete number when lowered for `to_numeric`. Each shape exercises a distinct - # arithmetic node in the constant folder. c0 = 1.0 + 2.0im cases = ( ("real", real(conj(z)), real(conj(c0))), @@ -568,24 +603,19 @@ using Test ) for (label, coeff, expected) in cases op = substitute(coeff * a, Dict(z => c0)) - @test to_numeric(op, b) ≈ expected * A + @test _dat(to_numeric(op, b)) ≈ expected * _dat(A) end - # An unrecognized constant function (`sin`) is not handled by the folder and - # falls back to the compile-based reduction path. op_sin = substitute(sin(z) * a, Dict(z => 0.5)) - @test to_numeric(op_sin, b) ≈ sin(0.5) * A + @test _dat(to_numeric(op_sin, b)) ≈ sin(0.5) * _dat(A) - # The folder also handles raw subtraction nodes. Symbolics canonicalizes - # subtraction and negation to `+`/`*`, so a literal `-` node never reaches - # the folder through normal arithmetic; build it directly to cover both arms. neg_unary = SymbolicUtils.term(-, 5.0 + 0im; type = Number) neg_binary = SymbolicUtils.term(-, 7.0 + 0im, 2.0 + 0im; type = Number) @test _fold_const(neg_unary) == -(5.0 + 0im) @test _fold_const(neg_binary) == (7.0 + 0im) - (2.0 + 0im) end - @testset "keyword to_numeric: op_type, vector, complex params, errors" begin + @testset "keyword to_numeric: vector, complex params, errors" begin h = FockSpace(:c) @qnumbers a::Destroy(h) b = FockBasis(4) @@ -593,16 +623,17 @@ using Test Ad = create(b) @variables x::Real z::Number E::Number - # op_type transforms each emitted operator. + # `op_type` selects the materialization; `dense`/`sparse` on the default result work too. Hd = to_numeric(2.0 * a' * a, b; op_type = dense) @test Hd isa QuantumOpticsBase.Operator @test Hd ≈ dense(2.0 * Ad * A) + @test dense(to_numeric(2.0 * a' * a, b)) ≈ Hd # Vector form forwards keywords to each element. - @test to_numeric([a, a'], b; parameter = Dict(x => 1.0)) == [A, Ad] + @test [_dat(x) for x in to_numeric([a, a'], b; parameter = Dict(x => 1.0))] == [_dat(A), _dat(Ad)] # A complex-valued parameter key is split into real/imaginary substitutions. - @test to_numeric(z * a, b; parameter = Dict(z => 2 + 3im)) ≈ (2 + 3im) * A + @test _dat(to_numeric(z * a, b; parameter = Dict(z => 2 + 3im))) ≈ (2 + 3im) * _dat(A) # Error paths. @test_throws ArgumentError to_numeric(x * a, b) # symbolic, no value @@ -612,6 +643,24 @@ using Test ) # untimed variable end + @testset "HilbertSpace form (uniform entry)" begin + h = FockSpace(:fock) + @qnumbers a::Destroy(h) + b = FockBasis(7) + ψ = coherentstate(b, 0.2 + 0.1im) + + # Backend defaults to the single loaded backend (QuantumOpticsBase here). + @test _dat(to_numeric(a' * a, h, 7)) == _dat(to_numeric(a' * a, b)) + @test numeric_average(a' * a, ψ) ≈ expect(to_numeric(a' * a, h, 7), ψ) + + # Composite HilbertSpace with per-subspace dims. + hp = FockSpace(:c) ⊗ NLevelSpace(:atom, 3, 1) + ap = Destroy(hp, :a, 1) + σp = Transition(hp, :σ, 1, 2, 2) + bp = FockBasis(4) ⊗ NLevelBasis(3) + @test _dat(to_numeric(ap' * σp, hp, (4, 3))) == _dat(to_numeric(ap' * σp, bp)) + end + @testset "op_type is shape-independent" begin # The return type depends only on op_type, not on the expression shape. h = FockSpace(:a) ⊗ FockSpace(:b) @@ -624,7 +673,7 @@ using Test @test to_numeric(expr, b) == to_numeric(expr, b; op_type = sparse) @test to_numeric(expr, b; op_type = dense) isa Operator end - @test to_numeric(exprs[1], b; op_type = identity) isa LazyTensor + @test to_numeric(exprs[1], b; op_type = identity) isa LazySum @test to_numeric(exprs[3], b; op_type = identity) isa LazySum for expr in exprs @test to_numeric(expr, b) == sparse(to_numeric(expr, b; op_type = identity)) @@ -632,4 +681,64 @@ using Test end end + @testset "numeric_average — LazyKet state" begin + # `expect` of a sparse operator has no method for a `LazyKet`; numeric_average + # assembles the lazy form, so a `LazyKet` state works without materializing. + hfock = FockSpace(:fock) + hnlevel = NLevelSpace(:nlevel, (:a, :b, :c)) + hprod = hfock ⊗ hnlevel + bfock = FockBasis(10) + bnlevel = NLevelBasis(3) + bprod = bfock ⊗ bnlevel + + @qnumbers a::Destroy(hprod, 1) + σ(i, j) = Transition(hprod, :σ, i, j, 2) + + α = 0.3 + 0.0im + ket_fock = coherentstate(bfock, α) + ket_nlevel = (nlevelstate(bnlevel, 1) + nlevelstate(bnlevel, 3)) / sqrt(2) + ψ_lazy = LazyKet(bprod, (ket_fock, ket_nlevel)) + ψ_dense = ket_fock ⊗ ket_nlevel + + for op in (a, a' * σ(:a, :c), a + a' * σ(:a, :c)) + @test numeric_average(op, ψ_lazy) ≈ expect(to_numeric(op, bprod), ψ_dense) + end + @test numeric_average(average(a' * a), ψ_lazy) ≈ + expect(to_numeric(a' * a, bprod), ψ_dense) + end + + @testset "numeric_average — unsupported symbolic operation" begin + # A symbolic average expression whose top operation is neither `+`, `*`, nor `^` + # (here `sqrt`) is not reducible to an expectation value and must error. + h = FockSpace(:fock) + @qnumbers a::Destroy(h) + b = FockBasis(5) + ψ = fockstate(b, 2) + @test_throws ArgumentError numeric_average(sqrt(average(a' * a)), ψ) + @test_throws ArgumentError numeric_average(average(a' * a) / (1 + average(a)), ψ) + end + + @testset "complex parameter key" begin + # A `Complex` parameter key `κ` splits into real/imag substitutions, so both `κ` + # and `conj(κ)` in the expression resolve from a single complex value. + h = FockSpace(:fock) + @qnumbers a::Destroy(h) + b = FockBasis(6) + @variables κ::Complex + H = κ * a + conj(κ) * a' + M = to_numeric(H, b; parameter = Dict(κ => 1.0 + 2.0im)) + Mref = to_numeric((1.0 + 2.0im) * a + (1.0 - 2.0im) * a', b) + @test _dat(M) ≈ _dat(Mref) + end + + @testset "time_parameter — unsupported key" begin + # A single-variable `time_parameter` key that is neither a bare variable nor + # `conj(v)` (here `2*E`) is rejected. + h = FockSpace(:fock) + @qnumbers a::Destroy(h) + b = FockBasis(6) + @variables E::Number + @test_throws ArgumentError to_numeric(E * a, b; time_parameter = Dict(2 * E => t -> 1.0 + 0im)) + end + end diff --git a/test/operators/collective_transition_test.jl b/test/operators/collective_transition_test.jl index cf423897..82d2d253 100644 --- a/test/operators/collective_transition_test.jl +++ b/test/operators/collective_transition_test.jl @@ -115,6 +115,10 @@ using Test @test_throws ArgumentError to_numeric(S(1, 2), wrong) @test_throws ArgumentError to_numeric(Destroy(FockSpace(:f), :a), b) + # A collective transition has no realisation on a plain (non-ManyBody) basis: the + # closed operator ladder falls through to the throwing `Val` extension point. + @test_throws ArgumentError to_numeric(S(1, 2), b1) + # su(2) collective transitions reproduce spin-N/2 generators. h2 = CollectiveNLevelSpace(:twolevel, 2) C(i, j) = CollectiveTransition(h2, :S, i, j) diff --git a/test/operators/phase_space_test.jl b/test/operators/phase_space_test.jl index 75b91999..0459958f 100644 --- a/test/operators/phase_space_test.jl +++ b/test/operators/phase_space_test.jl @@ -142,7 +142,7 @@ import SecondQuantizedAlgebra: simplify, QAdd, QSym, HilbertSpace bf = FockBasis(3) bq = FockBasis(5) bc = bf ⊗ bq - @test to_numeric(x, bc; op_type = identity) isa LazyTensor + @test to_numeric(x, bc; op_type = identity) isa LazySum @test to_numeric(x, bc) isa Operator end