Private / Internal API

Largest spin for which "S=k/2" site types are pre-defined (S = 3/2 to 10 in half-integer steps; "S=1/2" and "S=1" are native to ITensors).

AbstractBackend

Abstract base type for all simulation backends. A backend is the mutable struct that owns the numerical representation of the quantum state; the user-facing SimulationState is parametrically typed on it (SimulationState{B<:AbstractBackend}), so backend selection is pure multiple dispatch.

Concrete backends: MPSBackend (ITensor MPS), StateVectorBackend (dense Vector{ComplexF64}), CliffordBackend (QuantumClifford.jl stabilizer tableau).

A new backend must implement, for SimulationState{MyBackend}:

• initialize!(state, init) for each supported AbstractInitialState
• _apply_single!(state, gate, phy_sites) — the gate-application primitive (either one generic method that resolves gates to matrices/operators, or per-gate methods plus an informative-ArgumentError fallback)
• born_probability(state, site, outcome) — required by the default measurement path and the BornProbability observable
• _measure_single_site!(state, site) — only if the default (Born-sample + Projection) implementation in src/Core/apply.jl cannot work for the representation
• a method (or an explicit ArgumentError rejection) for every exported observable callable — unhandled observables fall through to MPS-assumed generics and crash with a FieldError

The full contract — required fields, method tables per backend, RNG-stream expectations, and indexing conventions — is documented in docs/src/devdocs/backend_interface.md.

AbstractFeedback

Supertype for measurement-feedback specifications attached to a Measure gate via its feedback= keyword.

Concrete subtypes:

• OnOutcome: declarative outcome → gate map
• CallbackFeedback: wraps a raw function (state, sites::Vector{Int}, outcome) -> Any (raw functions passed to Measure(...; feedback=f) are auto-wrapped)

Feedback runs INSIDE with_guarded_stream(registry, :gates_spacetime): any attempt to draw from :gates_spacetime during feedback throws (fixed-draw contract — see RNGRegistry). Feedback randomness must come from :gates_realization (e.g. HaarRandom) or :born_measurement (nested measurements). Feedback gates do NOT advance engine counters (gate_idx) and emit no GateApplied events; a nested Measure still emits its MeasurementOutcome.

Abstract type for initial state specifications.

AbstractStaircase <: AbstractGeometry

Base type for staircase geometries with internal pointer.

CallbackFeedback(f)

Closure escape hatch for measurement feedback. f is called as f(state, sites::Vector{Int}, outcome) after the Born sample; its return value is ignored. sites is always a Vector{Int} (single-site Measure passes [site]).

Raw functions passed to Measure(...; feedback=f) are wrapped in this type automatically — you rarely construct it yourself.

Recursion (feedback applying another Measure) is allowed and is the user's responsibility to terminate.

CircuitBuilder

Internal mutable structure for recording circuit operations during do-block construction.

Users never see this type directly - they interact via:

circuit = Circuit(L=4, bc=:periodic) do c
    apply!(c, Reset(), SingleSite(1))
    apply_with_prob!(c; outcomes=[
        (probability=0.5, gate=PauliX(), geometry=SingleSite(1))
    ])
end

The builder records operations as NamedTuples which are then passed to the Circuit constructor.

Fields

• L::Int: Number of physical sites
• bc::Symbol: Boundary conditions (:periodic or :open)
• operations::Vector{NamedTuple}: Accumulated operation records
• params::Dict{Symbol,Any}: User-defined parameters passed from outer constructor

CircuitEvent

Abstract supertype for all typed circuit events recorded in a SimulationState's opt-in event log (see SimulationState(...; log_events=true)).

Concrete subtypes:

• GateApplied: a gate was executed on specific sites
• MeasurementOutcome: a Born-sampled measurement produced an outcome

Access recorded events via events and measurements.

Clifford (stabilizer) backend: holds a QuantumClifford.jl stabilizer tableau.

Uses MixedDestabilizer (tracks both stabilizer and destabilizer generators), which enables efficient O(n²) measurement via project!.

GateApplied(step, op_idx, element_idx, gate_label, sites) <: CircuitEvent

Recorded when a gate is applied during circuit execution.

Fields:

• step::Int: circuit step index (1-based; 0 if unavailable at emission site)
• op_idx::Int: index of the operation within circuit.operations
• element_idx::Int: index of the geometry element within the operation (1-based)
• gate_label::String: human-readable gate label (see gate_label)
• sites::Vector{Int}: physical sites the gate acted on

MPS-based backend: holds the underlying matrix product state, site indices, and truncation parameters (SVD cutoff, maximum bond dimension).

MeasurementOutcome(step, op_idx, sites, outcome) <: CircuitEvent

Recorded when a Born-sampled projective measurement produces an outcome.

Fields:

• step::Int: circuit step index (0 when emitted from the low-level measurement primitive, which has no step context — Task 9 engine rewire threads real indices)
• op_idx::Int: operation index (0 when unavailable, see above)
• sites::Vector{Int}: measured physical sites
• outcome::Int: measurement outcome (0 or 1 for qubits)

SentinelRNG <: Random.AbstractRNG

Poison-pill RNG installed by with_guarded_stream in place of a guarded stream. Any attempt to draw from it throws an ErrorException explaining that the stream is forbidden in the current context (e.g. :gates_spacetime during measurement feedback, whose draw count must stay data-independent — see the fixed-draw contract in RNGRegistry).

State-vector backend: holds the full statevector as a dense complex vector.

engine selects the gate-application algorithm used by _apply_single!:

• :builtin (default): Tier 1, reshape + permutedims + matmul (see apply_gate_sv! in src/StateVector/StateVector.jl). Ground truth.
• :optimized: Tier 2, hand-written stride-loop kernel (see apply_gate_sv_optimized! in src/StateVector/optimized.jl). Faster, numerically verified to match Tier 1 bitwise/to <1e-13.

Base.copy(circuit::Circuit) -> Circuit

Deep copy of a circuit, including its geometries. This is THE per-trajectory pattern for thread safety: staircase/Pointer geometries are mutable (their positions advance during simulate!), so concurrent trajectories must each own a private copy of the circuit.

copy is implemented as deepcopy, which preserves INTRA-circuit geometry aliasing: if the same staircase object appears in several operations of the original, the copies also share one (new) object — required for CIPT-style shared random-walk positions.

Threads.@threads for seed in seeds
    c = copy(circuit)                       # private geometry state
    st = SimulationState(...; rng=RNGRegistry(gates_spacetime=seed, ...))
    initialize!(st, ProductState(binary_int=0))
    simulate!(c, st; n_steps=n)
end

Base.getproperty(state::SimulationState{MPSBackend}, name::Symbol)

Property forwarding for the MPS backend — SUPPORTED API (not a deprecation shim). For state::SimulationState{MPSBackend}, the four property names

• state.mps → state.backend.mps
• state.sites → state.backend.sites
• state.cutoff → state.backend.cutoff
• state.maxdim → state.backend.maxdim

forward transparently to the MPSBackend fields, for both reads and writes (Base.setproperty! forwards the same four names). All other property names resolve to SimulationState's own fields.

This convenience is MPS-only by design: SimulationState{StateVectorBackend} and SimulationState{CliffordBackend} have no such forwarding (their payloads are reached explicitly via state.backend.ψ / state.backend.tableau), so accessing state.mps on them raises a FieldError — a loud signal that MPS-specific code received a non-MPS state.

rand(registry::RNGRegistry, stream::Symbol) -> Float64

Draw a random Float64 from the specified stream. Mid-level API.

randn(registry::RNGRegistry, stream::Symbol, dims...) -> Array{Float64}

Draw random normal array from the specified stream. For Haar random matrices.

randn(registry::RNGRegistry, stream::Symbol) -> Float64

Draw a random normal Float64 from the specified stream.

Define exp(iπ Sz) operator for S=1 sites.

For spin-1: Sz = diag(+1, 0, -1) in basis |+1⟩, |0⟩, |-1⟩ exp(iπ Sz) = diag(exp(iπ), exp(0), exp(-iπ)) = diag(-1, 1, -1)

This is a diagonal operator in the Sz basis.

_apply_single!(state::SimulationState, gate::AbstractGate, phy_sites::Vector{Int})

Apply gate to specific physical sites. Internal workhorse.

Steps:

1. Validate support matches site count
2. Convert physical sites to RAM indices
3. Build operator with physical site indices
4. Apply operator to MPS
5. Normalize + truncate iff needs_normalization(gate) (trait, Contract 3.5)

_apply_single!(state::SimulationState, gate::BondParity, phy_sites::Vector{Int})

Generic (MPS/state-vector-backend) rejection fallback: BondParity has no gate_matrix/build_operator representation and is only implemented on backend=:gaussian. A later task adds the Gaussian-specific _apply_single!/execute! override (more specific on state), which Julia's dispatch prefers whenever the state actually uses GaussianBackend. The CliffordBackend's own AbstractGate fallback independently rejects BondParity on the Clifford backend with its own message.

_apply_single!(state::SimulationState, gate::GaussianHaar, phy_sites::Vector{Int})

Generic (MPS/state-vector-backend) rejection fallback: GaussianHaar has no gate_matrix/build_operator representation and is only implemented on backend=:gaussian. A later task adds the Gaussian-specific _apply_single!(state::SimulationState{GaussianBackend}, ::GaussianHaar, ...) override; Julia's dispatch prefers that method (more specific on state) over this fallback whenever the state actually uses GaussianBackend. The CliffordBackend's own AbstractGate fallback (src/Clifford/Clifford.jl) independently rejects GaussianHaar on the Clifford backend with its own message (more specific on state there too).

_apply_single!(state::SimulationState{GaussianBackend}, gate::AbstractGate, phy_sites::Vector{Int})

Fallback for any gate NOT handled by one of the specific _apply_single! methods above. The Gaussian (free-fermion covariance-matrix) backend can only represent fermionic Gaussian operations; generic qubit gates (e.g. Hadamard, CNOT, Haar-random qubit unitaries) are not Gaussian and have no covariance-matrix representation. Throws an informative ArgumentError naming the offending gate type and suggesting the dense-backend alternatives.

_apply_single!(state::SimulationState{GaussianBackend}, gate::BondParity, phy_sites::Vector{Int})

Dispatch disambiguator (always throws). BondParity is a projective measurement: on the Gaussian backend it is executed through the execute! measurement protocol (Gaussian execute! override, added by the measurement task), never through the _apply_single! gate path. This method exists so the AbstractGate catch-all below (specializing on state's type parameter) is not ambiguous against the un-parameterized _apply_single!(state::SimulationState, gate::BondParity, ...) rejection fallback in src/Gates/bond_parity.jl (specializing on gate's type) — the same ambiguity bug class previously fixed for the StateVector/Clifford backends (see the disambiguating overrides in that file).

_apply_single!(state::SimulationState{GaussianBackend}, gate::GaussianHaar, phy_sites::Vector{Int})

Apply a Haar-random SO(n) Majorana rotation to the two sites in phy_sites, where n is the total number of Majorana indices carried by the two sites (resolved granularity-aware via site_majoranas):

• fermionic-mode granularity (default, majoranas_per_site == 2): the two sites carry 4 Majoranas ix = [2a-1, 2a, 2b-1, 2b] → Haar-SO(4).
• Majorana-chain granularity (site_type="Majorana", majoranas_per_site == 1): the two sites ARE two Majoranas ix = [a, b] → Haar-SO(2). Haar on SO(2) is EXACTLY the uniform-angle rotation exp(θ γ_a γ_b) with θ ~ U[0, 2π) (SO(2) ≅ U(1), Haar = uniform angle) — the class-DIII unitary K_U of Pan, Shapourian, Jian, arXiv:2411.04191 (Eq. S-III.1; Python reference GTN.measure_all_tri_op's Υ = kraus((0, cos φ, sin φ)) branch). The φ ↔ rotation convention, DERIVED EMPIRICALLY from the Python golden cross-check (test/gaussian/testmajoranachain.jl): contracting kraus((0, cos φ, sin φ)) on the Majorana pair (a, b) equals direct conjugation Γ ← R Γ Rᵀ with R = [[cos φ, −sin φ], [sin φ, cos φ]] on rows/columns (a, b) (exact to machine precision). Since φ ~ U[0, 2π) makes R uniform over SO(2), the two parameterizations define the SAME ensemble.

The orthogonal matrix O = haar_orthogonal(rng, n) is drawn from the :gates_realization RNG stream (one draw per application, mirroring RandomClifford on the Clifford backend) and conjugated DIRECTLY onto the covariance matrix at the Majorana rows/columns ix:

Γ[ix, :] = O · Γ[ix, :]
Γ[:, ix] = Γ[:, ix] · Oᵀ

i.e. Γ ← R Γ Rᵀ with R = O ⊕ I — exact for Gaussian unitaries (no contraction kernel involved). The result is re-antisymmetrized and, if the purity diagnostic max |diag(Γ²) + 1| exceeds state.backend.purify_tol, re-purified via purify!.

_apply_single!(state::SimulationState{GaussianBackend}, gate::PauliX, phy_sites::Vector{Int})

Fermionic OCCUPATION FLIP on the single site in phy_sites — NOT the Jordan-Wigner qubit Pauli-X. Implemented as the reflection of one on-site Majorana: the sign of row and column 2i of Γ is flipped (i = RAM index of the site), which negates the on-site element Γ[2i-1, 2i] = ⟨iγ_{2i-1}γ_{2i}⟩ and therefore flips the occupation ⟨c†c⟩ = (1 - Γ[2i-1,2i])/2 between 0 and 1. A reflection Γ ← RΓR with R = diag(1,…,-1,…,1) is exactly orthogonal, so the pure-state invariant Γ² = -I is preserved to machine precision (no re-purification needed).

This is the operation that lets the generic Reset gate (measure, then flip if occupied — src/Core/apply.jl) work on the Gaussian backend: vacuum + PauliX(site i) ⇒ site i measures occupied with probability 1.

On the Majorana chain (site_type="Majorana", majoranas_per_site == 1) this throws an informative ArgumentError: a site is a single Majorana mode and no single-Majorana occupation flip exists.

_apply_single!(state::SimulationState{CliffordBackend}, gate::AbstractGate, phy_sites::Vector{Int})

Fallback for any gate NOT handled by one of the specific _apply_single! methods above. The Clifford (stabilizer-tableau) backend can only represent Clifford-group operations; non-Clifford gates (e.g. arbitrary rotations, Haar-random unitaries, or non-Clifford projections) have no native tableau representation. Throws an informative ArgumentError naming the offending gate type and suggesting the dense-backend alternatives.

_apply_single!(state::SimulationState{CliffordBackend}, gate::CNOT, phy_sites::Vector{Int})

Apply a CNOT gate to the stabilizer tableau via QuantumClifford's sCNOT.

_apply_single!(state::SimulationState{CliffordBackend}, gate::CZ, phy_sites::Vector{Int})

Apply a controlled-Z gate to the stabilizer tableau. NOTE: QuantumClifford.jl does NOT have a function named sCZ — the correct symbolic gate is sCPHASE, verified (via stabilizer conjugation of a +X_ state) to match the standard CZ = diag(1,1,1,-1) convention.

_apply_single!(state::SimulationState{CliffordBackend}, gate::Hadamard, phy_sites::Vector{Int})

Apply a Hadamard gate to the stabilizer tableau via QuantumClifford's sHadamard.

_apply_single!(state::SimulationState{CliffordBackend}, gate::PauliX, phy_sites::Vector{Int})

Apply a Pauli-X gate to the stabilizer tableau via QuantumClifford's sX.

_apply_single!(state::SimulationState{CliffordBackend}, gate::PauliY, phy_sites::Vector{Int})

Apply a Pauli-Y gate to the stabilizer tableau via QuantumClifford's sY.

_apply_single!(state::SimulationState{CliffordBackend}, gate::PauliZ, phy_sites::Vector{Int})

Apply a Pauli-Z gate to the stabilizer tableau via QuantumClifford's sZ.

_apply_single!(state::SimulationState{CliffordBackend}, gate::PhaseGate, phy_sites::Vector{Int})

Apply an S (phase) gate to the stabilizer tableau via QuantumClifford's sPhase.

_apply_single!(state::SimulationState{CliffordBackend}, gate::RandomClifford, phy_sites::Vector{Int})

Sample a random n-qubit Clifford operator from the :gates_realization RNG stream via QuantumClifford.random_clifford(rng, gate.n), then apply it NATIVELY to the stabilizer tableau (via QuantumClifford.apply!(tableau, qubit_indices, op)) — no dense-matrix conversion (no QuantumOpticsBase involvement), unlike the MPS/state-vector backends' gate_matrix path.

_apply_single!(state::SimulationState{CliffordBackend}, gate::SWAP, phy_sites::Vector{Int})

Apply a SWAP gate to the stabilizer tableau via QuantumClifford's sSWAP.

_apply_single!(state::SimulationState{StateVectorBackend}, gate::AbstractGate, phy_sites::Vector{Int})

Apply gate to specific physical sites on a dense state-vector backend. Internal workhorse (Tier 1 / vanilla engine — see applygatesv!).

Steps:

1. Validate support matches site count
2. Convert physical sites to RAM indices (identity for SV, kept for code-path consistency with the MPS backend)
3. Resolve the gate's dense matrix (HaarRandom needs rng + local_dim)
4. Apply via applygatesv! (Tier 1) or applygatesv_optimized! (Tier 2, selected via state.backend.engine == :optimized — see src/StateVector/optimized.jl), reassigning state.backend.ψ
5. Normalize iff needs_normalization(gate) (trait, Contract 3.5) — there is NO truncate! equivalent for state vectors (no bond dimension)

_bit_pattern_string(init::ProductState, L::Int, local_dim::Int) -> String

Derive the physical-site bit-pattern string (site 1 = MSB, site L = LSB) from a ProductState's binary_int/binary_decimal/bitstring specification, padded or truncated to L characters. Shared verbatim by all three backends' initialize!(::SimulationState, ::ProductState) methods (MPS: src/State/initialization.jl, state-vector: src/StateVector/initialization.jl, Clifford: src/Clifford/initialization.jl) — previously copy-pasted identically in all three.

Does NOT handle init.spin_state — callers must branch on spin_state before calling this function, since each backend treats it differently (MPS/state-vector build a uniform product state from it directly; Clifford rejects it since it is qubit-only). Assumes exactly one of binary_int/binary_decimal/bitstring is set (enforced by the ProductState inner constructor); the final else branch is an unreachable-in-practice defensive fallback preserved verbatim from the pre-refactor code.

local_dim is currently UNUSED by the derivation itself (only L affects padding/truncation) but is accepted as part of the signature for forward compatibility with an arbitrary spin-S generalization of ProductState (planned future work) — kept so that extension does not require touching call sites again.

_build_gate_operator(state, gate, phy_sites, ram_sites) -> ITensor

Build the operator tensor for the gate.

_evaluate_recording(record_when, ctx, step_idx, n_steps) -> (set_flag::Bool, record_now::Bool)

Evaluate recording criteria and return whether to set the recording flag and/or record immediately.

For compound geometries (called inside element loops), :every_gate triggers immediate recording. For simple geometries, the caller handles the immediate recording after setting the flag.

Returns a tuple:

• set_flag: Whether to set should_record_this_step = true
• record_now: Whether to call record!(state) immediately

_gaussian_region_majoranas(state, region) -> Vector{Int}

Map a (sorted, duplicate-free) collection of PHYSICAL sites to the sorted vector of its Majorana indices, granularity-aware via site_majoranas: fermionic-mode granularity — each site s contributes (2r−1, 2r) with r = state.phy_ram[s] (2·|region| indices); Majorana-chain granularity (site_type="Majorana") — each site contributes its own single Majorana index (|region| indices).

_haar_unitary(N::Int, rng::AbstractRNG) -> Matrix{ComplexF64}

Haar random N×N unitary via QR of a complex Ginibre matrix. EXACT reproduction of the CT.jl U() algorithm (lines 585-592): two separate randn(rng, N, N) calls (real then imaginary parts), QR, then phase correction by diag(R)/|diag(R)|. Do NOT change the call shape/order — the :gates_realization stream consumption is an API contract (goldens).

_infer_matrix_gate_dims(N::Int) -> (d, n)

Infer (local dimension, site count) from a MatrixGate size N per the documented rule: 2^n (n ≥ 1) → qubits; 3^n (n ≥ 2) → spin-1; otherwise ArgumentError.

_is_static_geometry(geo::AbstractGeometry) -> Bool

Trait: true for geometries whose elements(geo, L, bc) enumeration is provably STEP-INVARIANT (immutable structs with no advancing/movable position): Bricklayer, AllSites, EachSite, Sites, SingleSite, AdjacentPair.

false for mutable geometries — StaircaseLeft/StaircaseRight (positions advance during simulate!) and Pointer (moved explicitly via move!) — and, CONSERVATIVELY, for any geometry type not explicitly listed here (the fallback is false: correctness over performance). Any future compound geometry containing a member must only return true if EVERY member is static.

Used by simulate! to decide whether an op's elements() result may be cached across steps within one simulate! call. Never used to cache mutable geometries — no invalidation scheme exists by design.

_kraus(n::NTuple{3,<:Real}) -> Matrix{Float64}

4×4 Gaussian gate tensor Υ parameterized by n = (n₁,n₂,n₃) — port of GTN.kraus (GTN.py:105-118) with c = [1,1,1]. Unitary 2-Majorana rotation for ‖n‖ = 1; parity projection for n = (±1,0,0).

_local_basis_vector(site_type::String, local_dim::Int, name::String) -> Vector{Float64}

Resolve a single-site ITensor state name (e.g. "0", "1", "Up", "Z0", "Dn") into a dense local basis vector, WITHOUT constructing a full MPS. Used to build a state-vector-backend product state site-by-site.

Verified outputs: Qubit "0"->[1,0], "1"->[0,1]; S=1 "Up"->[1,0,0], "Z0"->[0,1,0], "Dn"->[0,0,1].

_matrix_gate_site_count(N::Int, d::Int) -> n

Exact site count n with d^n == N (n ≥ 1), or ArgumentError if N is not a positive power of d. Integer arithmetic only (no float log rounding).

_measure_single_site!(state::SimulationState, site::Int) -> Int

Perform Born-sampled projective measurement on a single site. Returns the measurement outcome (a level index 0 .. local_dim-1; 0 or 1 for qubits).

This is the FUNDAMENTAL measurement operation:

1. Compute Born probabilities P(k|ψ) for levels k = 0 .. local_dim-2
2. Sample ONE categorical outcome using a single scalar draw from the :bornmeasurement RNG stream (at localdim=2 this reduces EXACTLY to the historical binary draw rand < P(0) ? 0 : 1 — same draw count, same float comparison, bitwise-identical qubit trajectories)
3. Apply the per-level Projection operator
4. Return outcome (for conditional logic in Reset / feedback)

_measure_single_site!(state::SimulationState{GaussianBackend}, site::Int) -> Int

Override the default (Projection-based) _measure_single_site! for the Gaussian backend. Born-samples the occupation of site and collapses the covariance matrix onto the observed parity sector via gaussian_contraction! with the projector parity_projection_upsilon, mutating state.backend.corr in place.

Draw contract (REDUNDANT-DRAW, cross-backend lockstep): exactly ONE scalar :born_measurement draw is consumed per measured site — always, unconditionally, BEFORE any probability computation — matching the generic MPS/SV _measure_single_site! (Core/apply.jl) and the Clifford override (src/Clifford/measurement.jl) draw-for-draw. The outcome uses the shared threshold convention outcome = r < p₀ ? 0 : 1; probabilities are continuous here, so no deterministic branch is needed — the draw is always genuinely consumed by the comparison (for a deterministic site p₀ ∈ {0,1} and the comparison is vacuous, i.e. the draw is redundant, exactly as on the Clifford backend).

Outcome → parity sign mapping (VERIFIED empirically, T2/T8): the kernel projector parity_projection_upsilon(s) leaves the post-measurement state with Γ[2i−1, 2i] = −s. To end with Γ[2i−1,2i] = +1 (outcome 0, unoccupied) we contract with s = −1; for Γ[2i−1,2i] = −1 (outcome 1, occupied), s = +1. I.e. s = 2·outcome − 1.

Throws ArgumentError if the sampled outcome has probability ≤ 1e-15 (probability-zero branch — cannot be projected onto); gaussian_contraction! additionally raises its own singular-matrix ArgumentError as a second line of defense.

Logs a MeasurementOutcome event exactly like the default implementation (Core/apply.jl) does, and returns the outcome (0 or 1).

_measure_single_site!(state::SimulationState{CliffordBackend}, site::Int) -> Int

Override the default (Projection-based) _measure_single_site! for the Clifford backend. Uses QuantumClifford.projectZ! to measure site site in the Z basis, mutating state.backend.tableau in place.

Draw contract (REDUNDANT-DRAW, cross-backend lockstep): exactly ONE scalar :born_measurement draw is consumed per measured site — always, unconditionally — matching the generic MPS/SV _measure_single_site! (Core/apply.jl) draw-for-draw:

• Deterministic case (anticom_index == 0): the outcome is read directly off the returned phase byte. The draw is REDUNDANT here (the stabilizer formalism already fixes the outcome) and its value is DISCARDED — it is consumed purely to keep the :born_measurement stream position identical to the MPS/SV backends, where the (equally ignored) draw is structurally unavoidable.
• Undetermined case (anticom_index > 0): the outcome is genuinely random (50/50) and uses the same draw via the shared rand(rng) < p₀ ? 0 : 1 threshold convention (p₀ = 0.5 exactly for a stabilizer state); the phase of the anticommuting stabilizer row is set accordingly.

This deliberate redundant draw is what guarantees ABSOLUTE cross-backend reproducibility: same seeds ⇒ same measurement record on MPS, state-vector, and Clifford alike (audit T7 + T11, test/audit/cross_backend.jl (b), test/audit/born_measurement.jl (e)).

Logs a MeasurementOutcome event exactly like the default implementation (Core/apply.jl) does, and returns the outcome (0 or 1).

_mi_contiguous_region(r, which::String) -> UnitRange{Int}

Normalize a region spec (Int, range, or vector of integers) into a UnitRange{Int}, throwing an informative ArgumentError unless it is a non-empty, positive, strictly contiguous ascending set of sites. Used where contiguity is a CONSTRUCTION-time requirement (e.g. TripartiteMutualInformation's region layout).

_mi_entropy_from_probs(p, n, base, threshold) -> Float64

Entropy from a Schmidt/RDM probability spectrum p (need not be normalized; tiny negative eigenvalues from numerical RDMs are clamped at threshold^2). Mirrors _von_neumann_entropy's formula cases exactly: n=1 gives von Neumann −Σ p·logb(p); n≥2 gives Rényi logb(Σ pⁿ)/(1−n).

_mi_is_contiguous(v::Vector{Int}) -> Bool

true iff the (sorted) site vector is a plain ascending unit-step range.

_mi_region_vector(r, which::String) -> Vector{Int}

Normalize a region spec (Int, range, or vector of integers) into a sorted, duplicate-free Vector{Int} of positive sites, throwing an informative ArgumentError if empty, repeated, or non-positive. Contiguity is NOT required here — it is enforced per backend at evaluation time (see _validate_mutual_information); the Gaussian backend accepts arbitrary site subsets.

_mi_validate_bounds(mi::MutualInformation, state) -> nothing

Evaluation-time bounds check shared by ALL backends (including Gaussian): every region site must lie within 1:L. Regions are stored sorted, so checking last suffices.

_mps_subset_rdm_probs(state, phys_sites) -> Vector{Float64}

Eigenvalues of the reduced density matrix of the given PHYSICAL sites of an MPS state. Sites are mapped to RAM positions via phy_ram (handling the folded-PBC layout transparently); the RDM is built by sweeping the RAM span lo:hi, keeping open (primed) site indices on the requested sites and tracing all others. Left/right environments collapse to identity because the MPS is orthogonalized to lo (sites < lo left-orthonormal, sites > lo right-orthonormal).

_parse_spin_site_type(site_type::AbstractString) -> Union{Nothing, Rational{Int}}

Parse a spin site-type string of the form "S=<n>" or "S=<k>/2" and return the spin S as a Rational{Int} (e.g. "S=1" → 1//1, "S=3/2" → 3//2). Returns nothing for any string that is not a spin site type (e.g. "Qubit", "Qudit"), so callers can use it as a pattern test.

_product_region(pg::ProductGate, L::Int, bc::Symbol) -> Vector{Int}

Canonical region of a ProductGate: the sorted union of elements(pg.region_geometry, L, bc). This is the ONLY region a ProductGate may be applied to (API contract; see the ProductGate docstring).

_product_state_vector(site_type, local_dim, state_names_physical) -> Vector{ComplexF64}

Build the full dense product-state vector via Kronecker product of the L per-site local basis vectors, with site 1 FIRST in the kron chain (i.e. kron(v1, v2, ..., vL)). This makes physical site 1 the slowest-varying/MSB tensor index, matching the MPS backend's documented MSB convention (bit_pattern_str[1] is site 1's bit) and the SAME Kronecker convention documented in src/Gates/matrix_gate.jl (used by gate_matrix).

_projection_matrix(outcome::Int, d::Int) -> Matrix{ComplexF64}

Dense d×d projector |outcome⟩⟨outcome| (0-based level index), validated against the local dimension d.

_random_clifford_unitary(n::Int, rng::AbstractRNG; local_dim::Int=2) -> Matrix{ComplexF64}

Sample a Haar-random n-qubit Clifford operator via QuantumClifford.random_clifford(rng, n) and convert it to a dense local_dim^n × local_dim^n unitary matrix.

Conversion path: QuantumClifford.CliffordOperator → QuantumOpticsBase.Operator (via the QuantumCliffordQOpticsExt package extension, triggered by having both QuantumClifford and QuantumOpticsBase loaded) → dense .data matrix. This is the officially documented conversion route in QuantumClifford.jl (no direct Matrix(::CliffordOperator) method exists in QuantumClifford.jl itself).

_reset_circuit_geometries!(circuit::Circuit)

Reset all staircase geometry positions to their initial values. Called at the start of simulate! to ensure deterministic behavior.

_s_dot_s(s::Rational) -> Matrix{Float64}

The S₁·S₂ operator for two spin-s particles as a (2s+1)² × (2s+1)² matrix in the descending-m tensor product basis (generic version of s1_dot_s2, built from spin_operators(s)).

_spin_m_label(m::Rational) -> String

Format a magnetic quantum number m for use in "Z<m>" state names: integers render without denominator ("1", "-2"), half-integers as "<2m>/2" fractions ("3/2", "-1/2").

_sv_digit(n0, stride, d) -> Int

Extract the base-d digit of the 0-indexed basis integer n0 at the physical site whose precomputed stride is stride = d^(L - s) (site 1 = MSB convention shared by all state-vector observables). Callers hoist the loop-invariant d^(L - s) power out of their inner loops and pass it in.

_sv_spin_magnetization(state, s::Rational) -> Float64

Spin-site state-vector Magnetization: Mz = (1/L) Σᵢ ⟨Szᵢ⟩ with ⟨Sz⟩ = Σₖ (S-k)·P(digit=k), summed over all 2S+1 levels (level k, 0-based, has eigenvalue m = S-k; site 1 = MSB digit). One O(L·dᴸ) pass — spin SV states are small, so the simple double loop is fine.

_sv_subset_probs(state, phys_sites) -> Vector{Float64}

Schmidt probabilities (RDM eigenvalues) of an arbitrary subset of physical sites of a dense state vector: permute the kept sites' tensor dimensions to the front, reshape to (d^m × d^(L-m)), and square the singular values. The ordering of kept dimensions does not affect the spectrum.

_validate_mutual_information(mi::MutualInformation, state) -> nothing

Shared evaluation-time validation for the MPS, state-vector, and Clifford backends: region sites within 1:L AND each region CONTIGUOUS (these backends keep the historical contiguous-region contract). The Gaussian override calls _mi_validate_bounds only and accepts arbitrary site subsets (see src/Gaussian/mutual_information.jl).

_validate_pauli_string(obs::PauliString, state) -> nothing

Shared evaluation-time validation for all backends: qubit-only (v0.4.0) and site indices within 1:L.

_von_neumann_entropy(mps::MPS, i::Int; n::Int=1, threshold::Float64=1e-16, base::Float64=2.0) -> Float64

Compute entanglement entropy at bond i of an MPS.

Arguments:

• mps: The MPS state
• i: The bond index (site index) where entropy is computed
• n: Rényi index (1=von Neumann, n=Rényi-n)
• threshold: Minimum threshold for singular values to avoid log(0)
• base: Base of logarithm for entropy computation (default: 2.0 for bits)

Returns:

• Entanglement entropy value

The function:

1. Orthogonalizes the MPS to site i
2. Performs SVD on the tensor to extract Schmidt values
3. Computes probabilities from Schmidt values (squared)
4. Returns entropy based on Rényi index:
   • n=1: von Neumann entropy S₁ = -Σ p log_b(p)
   • n=0: Hartley entropy S₀ = log_b(rank)
   • n≠1: Rényi entropy Sₙ = log_b(Σ pⁿ) / (1-n)

_warn_bricklayer_odd_pbc(geo::Bricklayer, L::Int, bc::Symbol)

Emit a one-time warning (once per (parity, L) combination per session, via maxlog=1 with a per-combination log _id) when a single brickwork layer (:odd or :even) is used with odd L under periodic boundary conditions.

An odd ring has no valid brickwork tiling: the layer leaves one site unpaired (:odd → site L, :even → site 1 — no wrap pair is added, see elements(::Bricklayer, ...)), and the wrap bond (L, 1) is gated by neither single layer, so an alternating :odd/:even circuit is effectively open across that bond. Enumeration behavior is NOT changed by this helper — it only warns.

Called from the circuit-builder recording path (apply!(builder, ...), apply_with_prob!(builder; ...) — fires at circuit-definition time) and the immediate-mode dispatch path (_apply_dispatch!(state, gate, ::Bricklayer)). Deliberately NOT called inside elements() itself, which sits in performance-critical loops (benchmarks, per-element expansion).

:nn (ALL NN bonds — a bond enumeration, not a single layer) and the NNN sublayers are exempt; NNN odd-L coverage policy is a separate open question.

_xlogx(x::Real) -> Float64

x * log(x) with the exact limit 0 · log(0) = 0 (returns 0.0 for any x <= 0). Natural log. Used by subsystem_entropy so that exactly (un)occupied modes (λ ∈ {0, 1}) contribute zero entropy with no NaN/Inf and no imaginary-epsilon regularization.

advance!(geo::StaircaseLeft, L::Int, bc::Symbol)

Advance staircase left by one position. Internal use by apply!.

• StaircaseLeft: pos -= 1, wraps 1 → L (PBC) or 1 → L-1 (OBC)

advance!(geo::StaircaseRight, L::Int, bc::Symbol)

Advance staircase right by one position. Internal use by apply!.

• StaircaseRight: pos += 1, wraps L → 1 (PBC) or L-1 → 1 (OBC)

apply_feedback!(fb::AbstractFeedback, state, sites::Vector{Int}, outcome)

Internal: dispatch a feedback specification after a measurement. OnOutcome looks up the gate registered for outcome (no-op when absent) and executes it on sites; CallbackFeedback calls the wrapped function.

apply_gate_sv!(ψ, U, target_sites, L, d) -> Vector{ComplexF64}

Apply a d^n × d^n gate matrix U to target_sites (physical site numbers, 1-indexed, site 1 = MSB) of a length-d^L state vector ψ. Returns the new state vector (caller reassigns state.backend.ψ = result, does NOT mutate ψ's contents in-place — matches the existing state.backend.mps = ... reassignment pattern used elsewhere in the codebase).

Algorithm (verified numerically across multiple test cases, including an asymmetric CNOT gate applied in both site-argument orders):

1. Reshape ψ into an L-dimensional tensor of shape (d,d,...,d). Due to Julia's column-major array layout, tensor dimension k corresponds to PHYSICAL SITE L - k + 1 (site 1/MSB is the LAST dimension, site L/LSB is the FIRST dimension).
2. Convert each target physical site s to its Julia array dimension via L - s + 1, then REVERSE the resulting list (required so the merged multi-site index matches gate_matrix's "first site in the gate's argument list = slowest/most-significant digit" convention).
3. Permute dimensions so the (reversed) target dims come first, followed by all other dims (in ascending order).
4. Reshape the permuted tensor to 2D: (d^n, d^(L-n)) where n = length(target_sites).
5. Left-multiply by U: newmat = U * Amat.
6. Reshape back to L dims (d,d,...,d), then permutedims with the INVERSE of the original permutation.
7. vec(...) the result back to a flat Vector{ComplexF64}.

Handles 1-site AND n-site gates (adjacent or non-adjacent) UNIFORMLY — no special-casing needed for adjacent vs non-adjacent target sites.

apply_gate_sv_optimized!(ψ, U, target_sites, L, d) -> ψ

Tier 2 entry point: dispatches to the zero-allocation stride-loop kernel for 1-site gates, or the reshape/permutedims fallback for n-site (n>=2) gates. Mutates ψ in place and returns it (caller reassigns state.backend.ψ = apply_gate_sv_optimized!(...), matching Tier 1's reassignment call-site pattern even though Tier 2 itself mutates in place).

apply_gate_sv_optimized_1site!(ψ, U, site, L, d) -> ψ

In-place, zero-inner-loop-allocation 1-site gate application via the bit/digit-stride pattern (Yao.jl's u1apply!, generalized from d=2 to arbitrary local dimension d).

step = d^(L - site) is the stride between the d amplitudes that get mixed together by U at each combined "outer" index — this matches Tier 1's convention that physical site s corresponds to Julia tensor dimension L - s + 1 (site 1 = MSB = slowest-varying digit).

Mutates ψ in place and also returns it (for chaining / consistency with the rest of the Tier 2 API).

apply_gate_sv_optimized_nsite!(ψ, U, target_sites, L, d) -> ψ

n-site (n >= 2) fallback for the optimized engine: adjacent AND non-adjacent target sites, using the SAME reshape + permutedims + matmul algorithm as Tier 1's apply_gate_sv! (guaranteeing bitwise-identical results by construction), but mutating ψ IN PLACE via ψ .= vec(new_A) rather than allocating and returning a fresh vector. A hand-written stride-loop generalization to arbitrary n-site arity did not prove measurably beneficial over this approach (per plan Task 14 guidance: "the 1-site stride-loop alone... can deliver the bulk of the speedup" for the dominant gate types — Pauli, Hadamard, Rz, single-qubit Haar — so a hand-optimized kernel for every n-site arity is not required).

apply_op_internal!(mps::MPS, op::ITensor, sites::Vector{Index}, cutoff::Float64, maxdim::Int)

Apply operator to MPS following CT.jl algorithm (lines 147-172).

Contract 3.6: Index matching via Index comparison, NOT tag parsing.

build_operator(gate, site_or_sites, local_dim; kwargs...) -> ITensor

Build the operator tensor for a gate acting on given site(s).

build_operator(gate::CNOT, sites::Vector{Index}, local_dim::Int) -> ITensor

Build CNOT gate operator. Control = sites[1], target = sites[2]. Qubit-only (local_dim == 2).

Vectorized construction (T21): reshape the fixed gate_matrix(::CNOT) into a 4-index tensor, same convention as CZ/MatrixGate above. Verified bitwise identical to the prior loop-based construction (element-by-element ==).

build_operator(gate::CZ, sites::Vector{Index}, local_dim::Int) -> ITensor

Build CZ gate operator. Qubit-only (local_dim == 2).

Vectorized construction (T21): reshape the fixed gate_matrix(::CZ) into a 4-index tensor and hand it directly to ITensor(...), following the same output-primed-first, input-unprimed-second, reverse-site-order convention as MatrixGate (see matrix_gate.jl:98-110) — eliminates the 16 scalar setindex! calls of the previous per-application loop. Verified bitwise identical to the prior loop-based construction (element-by-element ==).

build_operator(gate::HaarRandom, site::Index, local_dim::Int; rng) -> ITensor

Single-site (n = 1) Haar random unitary. Same Ginibre-QR algorithm on a d × d matrix, consuming from :gates_realization.

build_operator(gate::HaarRandom, sites::Vector{Index}, local_dim::Int; rng::RNGRegistry) -> ITensor

Build an n-site Haar random unitary operator from the :gates_realization stream. For gate.n == 2 this is bit-identical to the historical implementation (same randn call shape/order, same index ordering).

build_operator(gate::MatrixGate, site::Index, local_dim::Int) -> ITensor

Single-site MatrixGate: ITensor(U, site', site) with primed = output.

build_operator(gate::MatrixGate, sites::Vector{Index}, local_dim::Int) -> ITensor

Reshape the stored d^n × d^n matrix into a 2n-index ITensor. Column-major reshape puts the LAST site of the region on the fastest matrix digit, so the tensor dims map to (sₙ', …, s₁', sₙ, …, s₁) with primed = output (row) and unprimed = input (column).

build_operator(gate::PauliX, site::Index, local_dim::Int) -> ITensor

Build Pauli X operator tensor.

build_operator(gate::PauliY, site::Index, local_dim::Int) -> ITensor

Build Pauli Y operator tensor.

build_operator(gate::PauliZ, site::Index, local_dim::Int) -> ITensor

Build Pauli Z operator tensor.

build_operator(gate::PhaseGate, site::Index, local_dim::Int) -> ITensor

Build Phase (S) operator tensor.

build_operator(gate::Projection, site::Index, local_dim::Int) -> ITensor

Build projection operator |outcome⟩⟨outcome| via the site type's per-level "Proj<k>" op string (defined for Qubit/"S=1/2" by ITensors and for all spin site types by src/Core/spin_sites.jl).

build_operator(gate::RandomClifford, site::Index, local_dim::Int; rng) -> ITensor

Single-site (n = 1) random Clifford unitary. Same conversion path on a d × d matrix, consuming from :gates_realization.

build_operator(gate::RandomClifford, sites::Vector{Index}, local_dim::Int; rng) -> ITensor

Build an n-site random Clifford unitary operator from the :gates_realization stream. Follows the same index-ordering convention as HaarRandom's build_operator (see two_qubit.jl): output-primed-first, input-unprimed-second, reverse-site-order.

build_operator(gate::SWAP, sites::Vector{Index}, local_dim::Int) -> ITensor

Build SWAP gate operator. Qubit-only (local_dim == 2).

Vectorized construction (T21): reshape the fixed gate_matrix(::SWAP) into a 4-index tensor, same convention as CZ/CNOT/MatrixGate above. Verified bitwise identical to the prior loop-based construction (element-by-element ==).

build_operator(gate::SpinSectorMeasurement, sites::Vector{Index}, local_dim::Int; rng, mps, ram_sites) -> ITensor

Build measurement operator for two spin-1 sites with proper Born sampling.

Implements the Born rule measurement:

1. Compute Born probabilities p(S) = ⟨ψ|P_S|ψ⟩ for each allowed sector S
2. Normalize probabilities over allowed sectors
3. Sample outcome S with probability p(S) using the provided RNG
4. Return the projector P_S onto the sampled sector

Arguments

• gate: SpinSectorMeasurement specifying allowed sectors
• sites: ITensor site indices for the two sites
• local_dim: Local Hilbert space dimension (d = 2s+1; 3 for spin-1)
• rng: RNG registry for reproducible sampling (uses :born_measurement stream)
• mps: Current MPS state for computing Born probabilities
• ram_sites: RAM indices of the two sites

Returns

• ITensor projector onto the randomly sampled spin sector

build_operator(gate::SpinSectorProjection, sites::Vector{Index}, local_dim::Int; kwargs...) -> ITensor

Build projector operator for two spin-1 sites. Returns ITensor representation of the projector matrix.

build_operator(gate::Union{Rx,Ry,Rz,Hadamard}, site::Index, local_dim::Int) -> ITensor

Build the exact 2×2 matrix operator: ITensor(M, site', site) with primed = output. Qubit-only (local_dim == 2).

build_template_groups_ascii(circuit; n_steps::Int=1) -> Vector{Vector{Vector{ExpandedOp}}}

Build operation groups for ASCII rendering by iterating circuit.operations directly, showing ALL outcomes of stochastic operations unconditionally (circuit template).

For deterministic ops: one group per apply! call (same as expandcircuitgrouped). For stochastic ops: one group per outcome (ALL outcomes shown, no random selection).

Returns the same steps → groups → ops structure as expandcircuitgrouped.

Compute physical-to-RAM and RAM-to-physical site mappings.

For OBC: identity mapping (1:L -> 1:L) For PBC: folded mapping - interleaves sites zig-zagging outward from pbc_fold_start

Parameters:

• L: system size
• bc: boundary condition (:open or :periodic)
• pbcfoldstart: physical site the PBC zig-zag fold starts from (default L÷4+1, the middle-aligned choice giving a contiguous half-cut). Ignored for bc == :open. Must satisfy 1 <= pbc_fold_start <= L for bc == :periodic.

Returns: (phyram, ramphy) where:

• phyram[physicalsite] = ram_index
• ramphy[ramindex] = physical_site

compute_pair_staircase(pos::Int, L::Int, bc::Symbol) -> Vector{Int}

Convert a single staircase position to an adjacent pair [pos, pos+1].

Arguments

• pos: Position (1-indexed)
• L: System size (number of sites)
• bc: Boundary condition (:periodic or :open)

Returns

• Two-element vector [pos, second] where second wraps at L→1 for PBC

Validation

• For OBC, pos must be in range 1:(L-1) (cannot form pair at position L)
• For PBC, pos can be 1:L (position L pairs with 1)

Example

compute_pair_staircase(3, 5, :periodic)  # Returns [3, 4]
compute_pair_staircase(5, 5, :periodic)  # Returns [5, 1] (wrap)
compute_pair_staircase(4, 5, :open)      # Returns [4, 5]
compute_pair_staircase(5, 5, :open)      # Throws error (no site 6)

compute_projector_product_expectation(state, sites_zero::Vector{Int}, site_one::Int) -> Float64

Compute ⟨ψ| (∏k P0k) P1 |ψ⟩ where P0 projects to |0⟩ and P1 projects to |1⟩. siteszero: physical sites that should be "0" siteone: physical site that should be "1"

This uses MPO construction for the projector product.

compute_site_staircase_left(start::Int, step::Int, L::Int, bc::Symbol) -> Int

Compute the position of a StaircaseLeft geometry after step iterations.

Arguments

• start: Starting position (1-indexed)
• step: Step number (1 = initial position, 2 = after first advance, etc.)
• L: System size (number of sites)
• bc: Boundary condition (:periodic or :open)

Returns

• Position after (step-1) advances

Wrapping behavior

• PBC: Position cycles L → ... → 2 → 1 → L
• OBC: Position cycles (L-1) → ... → 2 → 1 → (L-1)

Example

# Start at position 3, system size 5, periodic BC
compute_site_staircase_left(3, 1, 5, :periodic)  # Returns 3 (initial)
compute_site_staircase_left(3, 2, 5, :periodic)  # Returns 2 (after 1 advance)
compute_site_staircase_left(3, 3, 5, :periodic)  # Returns 1 (after 2 advances)
compute_site_staircase_left(3, 4, 5, :periodic)  # Returns 5 (after 3 advances, wrapped)

compute_site_staircase_right(start::Int, step::Int, L::Int, bc::Symbol) -> Int

Compute the position of a StaircaseRight geometry after step iterations.

Arguments

• start: Starting position (1-indexed)
• step: Step number (1 = initial position, 2 = after first advance, etc.)
• L: System size (number of sites)
• bc: Boundary condition (:periodic or :open)

Returns

• Position after (step-1) advances

Wrapping behavior

• PBC: Position cycles 1 → 2 → ... → L → 1
• OBC: Position cycles 1 → 2 → ... → (L-1) → 1

Example

# Start at position 3, system size 5, periodic BC
compute_site_staircase_right(3, 1, 5, :periodic)  # Returns 3 (initial)
compute_site_staircase_right(3, 2, 5, :periodic)  # Returns 4 (after 1 advance)
compute_site_staircase_right(3, 3, 5, :periodic)  # Returns 5 (after 2 advances)
compute_site_staircase_right(3, 4, 5, :periodic)  # Returns 1 (after 3 advances, wrapped)

compute_sites(geo::AdjacentPair, step::Int, L::Int, bc::Symbol) -> Vector{Int}

Compute sites for AdjacentPair geometry.

Arguments

• geo: AdjacentPair geometry
• step: Step number (unused for static geometry)
• L: System size
• bc: Boundary condition (:periodic or :open)

Returns

• Two-element vector [geo.first, second] where second wraps at L→1 for PBC

compute_sites(geo::SingleSite, step::Int, L::Int, bc::Symbol) -> Vector{Int}

Compute sites for SingleSite geometry (always returns [geo.site]).

Arguments

• geo: SingleSite geometry
• step: Step number (unused for static geometry)
• L: System size (unused for static geometry)
• bc: Boundary condition (unused for static geometry)

Returns

• Single-element vector [geo.site]

compute_sites(geo::StaircaseLeft, step::Int, L::Int, bc::Symbol, gate::AbstractGate) -> Vector{Int}

Compute sites for StaircaseLeft geometry based on gate support.

Arguments

• geo: StaircaseLeft geometry
• step: Step number
• L: System size
• bc: Boundary condition (:periodic or :open)
• gate: Gate to determine support (1-site or 2-site)

Returns

• For single-site gates (support == 1): [pos]
• For two-site gates (support == 2): [pos, pos+1] with PBC wrapping

compute_sites(geo::StaircaseRight, step::Int, L::Int, bc::Symbol, gate::AbstractGate) -> Vector{Int}

Compute sites for StaircaseRight geometry based on gate support.

Arguments

• geo: StaircaseRight geometry
• step: Step number
• L: System size
• bc: Boundary condition (:periodic or :open)
• gate: Gate to determine support (1-site or 2-site)

Returns

• For single-site gates (support == 1): [pos]
• For two-site gates (support == 2): [pos, pos+1] with PBC wrapping

compute_sites_dispatch(geo::AbstractGeometry, gate::AbstractGate, step::Int, L::Int, bc::Symbol) -> Vector{Int}

Dispatch compute_sites with appropriate arguments based on geometry type.

For StaircaseRight/StaircaseLeft: requires gate parameter to determine support. For SingleSite/AdjacentPair/Sites: gate parameter not needed.

compute_two_site_born_probability(mps::MPS, projector::Matrix{Float64}, ram_sites::Vector{Int}, local_dim::Int) -> Float64

Compute the Born probability ⟨ψ|P|ψ⟩ for a two-site projector.

Computes ⟨ψ|P|ψ⟩ by:

1. Orthogonalizing MPS to the target region
2. Contracting the two-site block into a local tensor
3. Computing ⟨T|P|T⟩ using proper index contraction

Arguments

• mps: Current MPS state
• projector: 9×9 projector matrix (d²×d² where d=local_dim)
• ram_sites: RAM indices of the two sites [i, j]
• local_dim: Local Hilbert space dimension (3 for spin-1)

Returns

• Probability p = ⟨ψ|P|ψ⟩ (real, non-negative)

current_position(p::Pointer) -> Int

Get the current position (read-only).

current_position(geo::AbstractStaircase) -> Int

Get the current position of the staircase (READ-ONLY accessor).

domain_wall(state::SimulationState, i1::Int, order::Int) -> Float64

Compute domain wall observable at sampling site i1 with given order. Ports CT.jl's dw_FM algorithm for xj=Set([0]) case.

The domain wall measures where the first "1" appears in the bit string when scanning cyclically from position i1.

domain_wall(state::SimulationState{StateVectorBackend}, i1::Int, order::Int) -> Float64

State-vector implementation of the domain_wall function (mirrors the formula of the existing MPS domain_wall in src/Observables/domain_wall.jl exactly: cyclic scan from i1, weight (L-j+1)^order, term j is ⟨ψ| (∏{k<j} P0k) P1_j |ψ⟩ over the cyclic site list).

Implementation: a single O(d^L) pass instead of L independent O(d^L) scans (one per term j). The L projector-product conditions are mutually exclusive in the computational basis: basis state n satisfies the term-j condition iff j is the FIRST position in the cyclic site list whose digit is nonzero AND that digit equals 1. So one sweep over the amplitudes finds, for each basis state, its unique contributing term (if any) by scanning for the first nonzero digit (amortized O(1) per basis state — a depth-k scan occurs for a d^-k fraction of basis indices). Each per-term accumulator receives exactly the same addends in the same ascending basis-index order as the former per-term scans, and the final weighted sum runs in the same j order, so the result is bitwise identical to the O(L·d^L) implementation.

draw(registry::RNGRegistry, stream::Symbol) -> Float64
draw(state::SimulationState, stream::Symbol) -> Float64

Draw ONE scalar uniform Float64 from the named RNG stream. This is the canonical way to consume a stream (replaces the removed, type-pirating rand(state, stream) extension; behavior is identical).

All engine coin draws are scalar — see the SCALAR-DRAW CONTRACT note in the RNGRegistry docstring.

Example: if draw(state, :gatesspacetime) < pctrl # control branch end

draw(state::SimulationState, stream::Symbol) -> Float64

Draw ONE scalar uniform Float64 from the named stream of the state's RNG registry (see draw(::RNGRegistry, ::Symbol)). Errors if no registry is attached.

execute!(state::SimulationState, gate::AbstractGate, region::Vector{Int})

Execute gate on the physical sites region (already resolved from a geometry). This is the v0.1 gate-execution protocol:

• Default implementation: build_operator → apply_op_internal!, then normalize + truncate iff needs_normalization(gate) (see _apply_single!).
• Gate-specific overrides: gates with non-operator semantics (Born sampling + classical logic) override this method — see Measure and Reset below. User-defined gates may override it the same way:

function QuantumCircuitsMPS.execute!(state::SimulationState, g::MyGate, region::Vector{Int})
    # custom behavior; region is a vector of physical sites
end

Related traits: needs_normalization(gate) (post-apply renormalization) and is_measurement(gate) (gate Born-samples via :born_measurement).

execute!(state::SimulationState, gate::Measure, region::Vector{Int})

Measure (v0.1, feedback-capable): Born-sample the single site in region (via _measure_single_site!, which emits the MeasurementOutcome event with the engine's real step/op indices), then dispatch the gate's feedback — if any — with the observed outcome.

Feedback runs INSIDE with_guarded_stream(registry, :gates_spacetime): it can never consume spacetime coins (fixed-draw contract), while :gates_realization and :born_measurement remain freely usable. Feedback gates execute through this same execute! protocol but are NOT engine op slots: they advance no counters and emit no GateApplied events. Recursion (feedback applying another Measure) is allowed — user responsibility.

execute!(state::SimulationState, pg::ProductGate, region::Vector{Int})

Execute a ProductGate: validate that region covers exactly the union of elements(pg.region_geometry, state.L, state.bc), then apply pg.inner to each element in canonical order via the uniform execute! protocol (one fresh :gates_realization draw per element for random inner gates). Logs one GateApplied per element (inner gate's label) when the event log is enabled.

execute!(state::SimulationState, gate::Reset, region::Vector{Int})

Reset (DERIVED - measurement + conditional X): Born-sample the single site in region; if the outcome is 1, flip it back to |0⟩ with PauliX.

execute!(state::SimulationState{GaussianBackend}, gate::BondParity, phy_sites::Vector{Int})

Projective measurement of the BOND fermion parity i γ_a γ_b on the two INNER Majoranas of an adjacent-site bond, for the fermionic Gaussian backend.

Inner-Majorana convention: site i carries Majoranas (2i−1, 2i) (RAM-mapped). For the bond between sites (i, i+1) the measured pair is ix = [2i, 2i+1] — site i's SECOND Majorana and site i+1's FIRST Majorana. For the periodic wrap bond (L, 1) (only valid when state.bc == :periodic) the pair is ix = [2L, 1] — site L's second Majorana and site 1's first. Sites may be passed in either order ([i, i+1] or [i+1, i]); non-adjacent pairs throw ArgumentError, as does the wrap bond under bc = :open.

Majorana-chain granularity (site_type="Majorana", majoranas_per_site == 1): each site IS one Majorana, so the bond (i, i+1) measures the pair ix = [i, i+1] directly — this is the i γ_i γ_{i+1} parity measurement of the class-DIII monitored Majorana chain (wrap bond (L, 1) → ix = [L, 1] under PBC). Same Born rule, same kernel, same draw contract — only the site→Majorana index mapping (via site_majoranas) differs.

Born rule (VERIFIED empirically vs the T5 ED/Pfaffian oracle, T9): the covariance element g = Γ[ix₁, ix₂] satisfies ⟨i γ̂_{ix₁} γ̂_{ix₂}⟩ = −g, so with the outcome encoding outcome ∈ (0, 1) ↔ parity eigenvalue s = 2·outcome − 1 ∈ (−1, +1):

P(outcome = 0) = (1 + g)/2        P(outcome = 1) = (1 − g)/2

— the SAME affine structure as the on-site occupation measurement (_measure_single_site!), because on-site occupation is itself the bond parity of the intra-mode pair (n = (1 + iγ₁γ₂)/2). The kernel s IS the measured parity eigenvalue: contracting parity_projection_upsilon(s) leaves Γ[ix₁, ix₂] = −s, i.e. ⟨i γ̂ γ̂⟩ = s post-measurement.

Draw contract (REDUNDANT-DRAW, cross-backend lockstep): exactly ONE scalar :born_measurement draw is consumed per BondParity application — always, unconditionally, BEFORE any probability computation — matching _measure_single_site! draw-for-draw.

Throws ArgumentError if the sampled outcome has probability ≤ 1e-15. Logs a MeasurementOutcome event exactly like Measure's execute! (via _measure_single_site!) does. Returns nothing.

gate_label(gate::AbstractGate) -> String

Return a short visualization label for a gate type.

Labels

• Reset → "Rst"
• HaarRandom → "Haar"
• RandomClifford → "Cl"
• Projection → "Prj"
• PauliX → "X"
• PauliY → "Y"
• PauliZ → "Z"
• CZ → "CZ"
• Other → Type name as string

Examples

gate_label(Reset())       # Returns "Rst"
gate_label(HaarRandom())  # Returns "Haar"
gate_label(CZ())          # Returns "CZ"

gate_matrix(gate::AbstractGate) -> Matrix{ComplexF64}

Return the explicit matrix of a gate whose action is defined by a fixed matrix (e.g. MatrixGate, Rx, Ry, Rz, Hadamard). Matrix convention: U[out, in] = ⟨out|U|in⟩ with Kronecker (row-major site) ordering — the FIRST site of the region is the slowest basis digit, so kron(A, B) acts with A on the first site.

gate_matrix(g::HaarRandom, rng::AbstractRNG; local_dim::Int=2) -> Matrix{ComplexF64}

State-vector-path equivalent of build_operator(gate::HaarRandom, ...): draws a fresh d^n × d^n Haar random unitary (d = local_dim, n = g.n) by reusing the same _haar_unitary core used by the MPS build_operator path. Consumes from whichever RNG stream rng is (caller is responsible for passing the appropriate stream, e.g. :gates_realization).

gate_matrix(g::RandomClifford, rng::AbstractRNG; local_dim::Int=2) -> Matrix{ComplexF64}

State-vector-path equivalent of build_operator(gate::RandomClifford, ...): draws a fresh d^n × d^n Haar-random Clifford unitary (d = local_dim, n = g.n) by reusing the same _random_clifford_unitary core used by the MPS build_operator path. Consumes from whichever RNG stream rng is (caller is responsible for passing the appropriate stream, e.g. :gates_realization).

gate_matrix(gate::SpinSectorProjection) -> Matrix{ComplexF64}

State-vector-path equivalent of build_operator(gate::SpinSectorProjection, ...): the 9×9 projector as a dense complex matrix (audit fix, T17 — this method was missing, so apply! on backend=:statevector threw MethodError and the README AKLT protocol could not run on the SV backend).

Basis ordering matches the state-vector engine convention ("first site in the gate's argument list = slowest/most-significant digit", see apply_gate_sv!): projector row/col = (m1-1)*3 + m2 with m1 = first site (slow), m2 = second site (fast). Renormalization after application is provided by the needs_normalization trait above (honored by the SV _apply_single! path).

gaussian_contraction!(Γ, Υ, ix; scratch=nothing, purify_tol=1e-10) -> Γ

Contract the Gaussian gate tensor Υ (2m×2m for m = length(ix) Majorana legs) into the covariance matrix Γ at Majorana indices ix (1-based), in place — line-by-line port of P_contraction_2 (GTN.py:1063-1102):

C = (Γ_RR·Υ_LL + I)⁻¹
Γ[ix̄,ix̄] += Γ_LR·(Υ_LL·C)·Γ_LRᵀ
Γ[ix,ix̄]  = Υ_RL·C·Γ_LRᵀ
Γ[ix,ix]  = Υ_RR + Υ_RL·(Γ_RR·Cᵀ)·Υ_RLᵀ
Γ[ix̄,ix]  = −Γ[ix,ix̄]ᵀ

Throws ArgumentError when Γ_RR·Υ_LL + I is singular (a probability-zero measurement outcome) — no least-squares fallback. If the result drifts off the pure-state manifold beyond purify_tol (max |diag(Γ²) + 1|), calls purify! and re-antisymmetrizes. scratch (same size as Γ) is an optional preallocated buffer for the ix̄-block update.

get_all_sites(geo::AllSites, state) -> Vector{Int}

Get all physical sites (1:L).

get_compound_elements(geo::AbstractGeometry, L::Int, bc::Symbol) -> Vector{Vector{Int}}

Get elements for compound geometry iteration.

Legacy name — delegates to the canonical elements; enumeration order is identical (bit-for-bit API contract).

Examples

# Bricklayer with odd parity on L=4 system
geo = Bricklayer(:odd)
get_compound_elements(geo, 4, :open)  # [[1, 2], [3, 4]]

# AllSites on L=3 system
geo = AllSites()
get_compound_elements(geo, 3, :open)  # [[1], [2], [3]]

get_op_ram_sites(op::ITensor, sites::Vector{Index}) -> Vector{Int}

Get RAM site indices from operator indices using Index comparison (Contract 3.6). Does NOT parse tags.

get_pairs(geo::Bricklayer, state) -> Vector{Tuple{Int,Int}}

Get all pairs for bricklayer pattern. Returns pairs of physical sites.

Legacy name — delegates to the canonical elements(geo, L, bc) (see Geometry/elements.jl); enumeration order is identical (bit-for-bit API contract).

get_sites(geo::AbstractGeometry, state) -> Vector{Int}

Get the physical sites for this geometry. Returns physical site indices. For iterating geometries (Bricklayer, AllSites), returns the first/next set. For staircases, returns current position pair.

get_sites(p::Pointer, state) -> Vector{Int}

Get the current pair of sites [pos, pos+1] for two-qubit gates. Handles PBC wrap when pos == L.

get_sites(geo::AbstractStaircase, state) -> Vector{Int}

Get current pair of physical sites for the staircase. Returns [pos, pos+range] with proper boundary condition handling.

haar_orthogonal(rng::AbstractRNG, n::Int) -> Matrix{Float64}

Exact Haar-random special-orthogonal matrix Q ∈ SO(n): QR-decompose a Ginibre matrix, fix the R-diagonal sign ambiguity (Q ← Q·diag(sign(diag(R))), required for Haar on O(n)), then flip the first column if det(Q) < 0 to land in SO(n).

is_aliased(registry::RNGRegistry) -> Bool

true when the :gates_spacetime and :gates_realization streams are the SAME RNG object (===), as in RNGRegistry(Val(:ct_compat); ...). Aliased registries are exempt from the fixed-draw invariant and from the with_guarded_stream sentinel guard (guarding one alias would also block the other; interleaved consumption is CT.jl-faithful by design).

is_compound_geometry(geo::AbstractGeometry) -> Bool

Check if geometry requires element-by-element iteration.

Compound geometries (Bricklayer, AllSites, EachSite) need to be expanded into multiple individual gate applications, one for each element/site pair.

Equivalent to the v0.1 trait is_broadcast; kept as a legacy alias until engine call sites are rewired.

is_measurement(gate::AbstractGate) -> Bool

Trait: does this gate perform a Born-rule measurement (consuming the :born_measurement RNG stream and collapsing the state)? Default false.

true for Measurement and SpinSectorMeasurement. Note Reset is a derived operation (measure + conditional flip) implemented via its own execute! method; it reports false here because the gate as a whole is a deterministic-outcome channel, not an outcome-reporting measurement.

is_record_mark(op::ExpandedOp) -> Bool

True when op is a record! marker pseudo-op (gate-less annotation produced from a (type=:record_mark, ...) circuit operation). Marker ops carry no gate (gate === nothing) and no sites; visualization renders them as a marker row (▽ in Unicode mode, [R] in ASCII mode).

log_event!(state::SimulationState, ev::CircuitEvent)

Append a CircuitEvent to the state's event log. Silent no-op when the log is disabled (the default), so emission sites can call this unconditionally. Returns nothing.

majorana_indices(site::Int) -> (Int, Int)

Majorana indices (2site-1, 2site) for 1-based fermionic mode site.

needs_normalization(gate::AbstractGate) -> Bool

Trait: does applying this gate require renormalizing (and truncating) the MPS afterwards? Default false (unitaries preserve the norm).

true for projective/collapsing gates (Projection, SpinSectorProjection, SpinSectorMeasurement). User-defined non-unitary gates opt in with:

QuantumCircuitsMPS.needs_normalization(::MyProjectiveGate) = true

occupation_covariance(bits::AbstractVector{Bool}) -> Matrix{Float64}

Product-state covariance matrix for occupations bits: mode i gets block [[0,1],[−1,0]] when bits[i] == false (⟨c†c⟩ = 0) and the sign-flipped block [[0,−1],[1,0]] when bits[i] == true (⟨c†c⟩ = 1). Sign ↔ occupation mapping verified empirically against Python get_C_f (GTN.py:1518-1523).

pack_ops_into_layers(ops::Vector{ExpandedOp}) -> Vector{Vector{ExpandedOp}}

Greedily pack a flat list of expanded operations into non-overlapping layers.

Each layer contains operations that act on disjoint sets of qubits. Operations are placed using a first-fit strategy: each op is assigned to the first existing layer where it does not conflict with any already-placed op, or a new layer is opened if no such layer exists.

Conflict detection uses set intersection on op.sites, so wrapped pairs like [8, 1] are handled correctly.

Arguments

• ops: Flat list of ExpandedOp values to pack (e.g., from expand_circuit).

Returns

A Vector{Vector{ExpandedOp}} where each inner vector is one non-overlapping layer. Returns an empty vector when ops is empty.

Examples

# Bricklayer(:nn) on L=8 periodic → 2 layers of 4 ops each
ops = [ExpandedOp(1, HaarRandom(), [1,2], "Haar"),
       ExpandedOp(1, HaarRandom(), [3,4], "Haar"),
       ExpandedOp(1, HaarRandom(), [5,6], "Haar"),
       ExpandedOp(1, HaarRandom(), [7,8], "Haar"),
       ExpandedOp(1, HaarRandom(), [2,3], "Haar"),
       ExpandedOp(1, HaarRandom(), [4,5], "Haar"),
       ExpandedOp(1, HaarRandom(), [6,7], "Haar"),
       ExpandedOp(1, HaarRandom(), [8,1], "Haar")]
layers = pack_ops_into_layers(ops)
length(layers)       # 2
length(layers[1])    # 4
length(layers[2])    # 4

parity_projection_upsilon(s::Int) -> Matrix{Float64}

4×4 Υ projecting a Majorana pair (i,j) onto parity outcome s ∈ {+1,−1}, i.e. _kraus((s,0,0)). Post-measurement state has Γ[i,j] = −s (verified against Python: vacuum Γ[2i−1,2i] = +1 is the s = −1 outcome; contracting the vacuum pair with s = +1 is the probability-zero outcome and throws).

purify!(Γ::Matrix{Float64}) -> Matrix{Float64}

Project a (numerically drifted) covariance matrix back onto the manifold of pure Gaussian states (Γ² = −I) — port of purify (GTN.py:1131-1136, see App. B2 of PhysRevB.106.134206). Eigendecomposes the Hermitian matrix Γ/i, clamps eigenvalues to ±1, reconstructs −Im(V·diag(vals)·V†), then antisymmetrizes. In-place; returns Γ.

record_value(obs, state; i1=nothing) -> Number | AbstractVector

Observable-level recording hook: compute the value record! should append for obs. The default is simply obs(state) — both for AbstractObservable instances and for generic callables (the untyped fallback), so ANY callable registered via track! is recorded without further ceremony.

Observables that need extra recording-time context override this method instead of being special-cased inside record! — e.g. DomainWall resolves its sampling site from i1_fn (registration-time closure) or the explicit i1 keyword. User-defined observables may override it the same way:

QuantumCircuitsMPS.record_value(obs::MyObs, state; i1=nothing) = ...

This is a MECHANISM hook only — it must return the value the observable's call protocol defines: a Float64 for the built-in scalar observables, or more generally any Number or AbstractVector (see track! for the storage contract).

reset!(::AbstractGeometry)

No-op base method for reset! on general geometries.

reset!(geo::AbstractStaircase)

Reset staircase position to its initial value.

s1_dot_s2()

Compute the S₁·S₂ operator for two spin-1 particles. Returns a 9×9 matrix in the tensor product basis.

S₁·S₂ = Sz₁⊗Sz₂ + (1/2)(S+₁⊗S-₂ + S-₁⊗S+₂)

select_outcome_index(rng::AbstractRNG, probs::Vector{Float64}) -> Int

SINGLE SOURCE OF TRUTH for the v0.1 unified stochastic rule's categorical selection. Draws exactly ONE scalar coin from rng and returns the 1-based index of the selected outcome, or 0 for the identity remainder (r >= Σp).

Semantics (bit-identical to test/reference_rule.jl's reference_select for a single element):

• one scalar rand(rng) per call — never vectorized draws
• cumulative walk over probs with strict <
• cumsum snapping: when abs(sum(probs) - 1) <= 1e-10, the LAST cumulative boundary is snapped to exactly 1.0, so float dust in Σp cannot leak spurious identity selections

expand.jl's visualization path (Task 15) also delegates to this function — there is exactly one categorical-selection implementation in the package.

set_event_context!(state::SimulationState, step::Int, op_idx::Int, element_idx::Int)

Record the engine's current execution position on the state. simulate! calls this before every execute! invocation so that events emitted from deep call sites (e.g. MeasurementOutcome from the Born-measurement primitive) carry real (step, op_idx) indices. Feedback gates (Task 10) run inside their measuring gate's execute! and inherit its context — they never advance counters or context themselves.

Costs three integer stores; never allocates.

site_majoranas(state, site::Int) -> NTuple{N,Int}

Majorana indices carried by physical site of a Gaussian-backend SimulationState, RAM-mapped through state.phy_ram (identity on the Gaussian backend, kept for protocol uniformity). THE single source of truth for the site → Majorana index mapping — every Gaussian gate / measurement / observable resolves sites through this helper:

• fermionic-mode granularity (state.backend.majoranas_per_site == 2, default): returns (2r−1, 2r) — the pair of one fermionic mode.
• Majorana-chain granularity (majoranas_per_site == 1, site_type="Majorana"): returns (r,) — the site IS one Majorana mode.

Duck-typed on state (needs only phy_ram and backend.majoranas_per_site) so this pure-kernel file keeps no SimulationState dependency.

spin1_operators()

Return the spin-1 operators Sz, S+, S- as 3×3 matrices. Basis ordering: |+1⟩, |0⟩, |-1⟩ (descending m).

spin_operators(s) -> (Sz, Sp, Sm)

Return the spin-s operators Sz, S+, S- as dense (2s+1)×(2s+1) matrices in the descending-m basis |s,s⟩, |s,s-1⟩, ..., |s,-s⟩ (matching ITensors' spin site conventions and spin1_operators).

Standard ladder formulas: Sz|s,m⟩ = m|s,m⟩ and S±|s,m⟩ = √(s(s+1) − m(m±1)) |s,m±1⟩.

s may be any positive integer or half-integer (3//2, 1.5, 2, ...).

subsystem_entropy(Γ::AbstractMatrix{<:Real}, majorana_idx::AbstractVector{Int}) -> Float64

Von Neumann entropy IN NATS (natural log) of the reduced state of a fermionic Gaussian state on the subsystem spanned by the Majorana indices majorana_idx, computed from the full 2L×2L covariance matrix Γ.

Port of von_Neumann_entropy_m (GTN.py:753-759):

1. Γ_A = Γ[majorana_idx, majorana_idx] — covariance submatrix (plain fancy indexing; majorana_idx need not be contiguous, which is what makes this helper reusable for MutualInformation on disjoint regions).
2. ξ = eigvals(Hermitian(im .* Γ_A)) — real spectrum in [−1, 1] (i·Γ_A is exactly Hermitian since Γ is exactly antisymmetric).
3. λ = clamp.((1 .- ξ) ./ 2, 0.0, 1.0) — occupation eigenvalues in [0, 1] (clamped: float noise can push ξ marginally outside [−1, 1]).
4. S = -Σ [λ log λ + (1−λ) log(1−λ)] / 2 — the division by 2 compensates the double-counting of the 2|A| Majorana eigenvalues (they come in ±ξ pairs; Python divides the full sum by 2 the same way).

Callers wanting a different log base divide the result by log(base).

support(gate::AbstractGate) -> Int

Return the number of sites this gate acts on (1 for single-qubit, 2 for two-qubit).

sync_staircase_positions!(outcomes, selected_geo::AbstractStaircase)

After advancing the selected staircase in a stochastic group, sync all other AbstractStaircase geometries in the group to the selected geometry's new position. This implements shared random-walk position for CIPT circuits.

vacuum_covariance(L::Int) -> Matrix{Float64}

2L×2L Majorana covariance matrix Γ₀ = ⊕ᵢ [[0,1],[−1,0]] of the all-modes- unoccupied vacuum (port of correlation_matrix, GTN.py:29-55, non-op branch). Verified via Python get_C_f: every mode has ⟨c†c⟩ = 0.

validate_geometry(geo::AbstractGeometry)

Validate that a geometry type is supported for circuit expansion.

Supported:

• StaircaseRight, StaircaseLeft
• SingleSite, AdjacentPair, Sites
• Bricklayer, AllSites, EachSite

Throws

ArgumentError if geometry is not supported (e.g. Pointer, whose position depends on runtime measurement outcomes and cannot be expanded statically).

validate_support(gate::AbstractGate, geo::Sites, L::Int, bc::Symbol)

Validate that support(gate) equals the region size of a Sites geometry; throws ArgumentError naming both numbers on mismatch. For all other geometries this is a no-op — their regions are resolved from the gate support itself (staircases) or are fixed-size by construction.

with_guarded_stream(f, registry::RNGRegistry, stream::Symbol)

Run f() with registry's stream temporarily replaced by a SentinelRNG: any draw from the guarded stream inside f throws an ErrorException ("... forbidden ..."), while all other streams remain fully usable. The original stream object is ALWAYS restored (try/finally), even if f throws. Returns f()'s value.

Used by the feedback system to guarantee that measurement feedback cannot consume :gates_spacetime coins (fixed-draw contract).