prism_q/sim/

mod.rs

//! Simulation orchestration: routes a circuit to a backend and answers the
//! terminal queries. Entry points: [`simulate`], [`run_qasm`], [`run_on`].

pub mod braket;
mod branching;
pub mod calibration;
mod clifford_register;
pub mod compiled;
mod decomposed;
mod dispatch;
mod dynamic;
pub mod gradient;
pub mod homological;
mod metadata;
pub mod noise;
mod observable;
mod probability;
pub(crate) mod shots;
pub mod stabilizer_rank;
mod terminal_sampling;
mod trajectory;
pub mod unified_pauli;

pub use braket::ResultValue;
pub(crate) use decomposed::merge_probabilities;
use decomposed::{
    MIN_DECOMPOSITION_QUBITS, PreparedBlocks, run_decomposed, run_decomposed_prefused,
    should_decompose,
};
pub use dispatch::BackendKind;
use dispatch::{
    AUTO_MPS_BOND_DIM, BackendPlan, ExecutionPlan, Family, MAX_AUTO_T_COUNT_SHOTS,
    MAX_STABILIZER_RANK_QUBITS, MIN_BLOCK_FOR_FACTORED_STAB, MIN_FACTORED_STABILIZER_QUBITS,
    MIN_QUBITS_FOR_SPD_AUTO, accel_for, approximate_route_name, auto_selects_cpu_statevector,
    auto_spd_work_budget, build_statevector, has_temporal_clifford_opportunity, initial_state_plan,
    mps_apply_cost, plan_for_family, plan_temporal_clifford, resolve, resolve_backend,
    run_temporal_clifford, stabilizer_rank_budget, validate_explicit_backend,
};
pub use dynamic::{SimulateProgram, simulate_program};
pub use metadata::{
    BondReport, Engine, Exactness, ExpectationResult, Placement, ResolvedBackend, RunMetadata,
};
#[cfg(feature = "distributed")]
pub(crate) use observable::pauli_sandwich;
pub use observable::{ObservableExpectation, PauliObservable};
pub(crate) use observable::{
    finish_expectations, i_pow, pauli_expectation_from_masks, pauli_expectations_from_masks,
    pauli_masks, pauli_sandwiches_from_masks, validate_observable,
};
pub use probability::{FactoredBlock, Probabilities, ProbabilitiesIter};
pub use shots::{ShotsResult, bitstring};

use std::collections::HashMap;

use num_complex::Complex64;

use crate::backend::sparse::MAX_SPARSE_INDEX_QUBITS;
use crate::backend::statevector::StatevectorBackend;
use crate::backend::{Backend, max_statevector_qubits};
use crate::circuit::{Circuit, Instruction, SaveRecord, SaveSpec, SavedValue};
use crate::error::{PrismError, Result};
use crate::sim::noise::NoiseModel;
use shots::{
    counts_from_basis_samples, packed_shots_to_classical_bits, record_counts_to_classical_bits,
    sample_shots, shots_from_basis_samples,
};
use terminal_sampling::{
    sample_counts_from_probs, sample_counts_from_state, sample_shots_from_probs,
    sample_shots_from_state,
};
use unified_pauli::{PauliTerm, SpdTruncation};

type TerminalStatevector = (StatevectorBackend, Vec<(usize, usize)>);

#[derive(Debug, Clone, Copy)]
pub(crate) struct SimOptions {
    pub(crate) probabilities: bool,
}

impl Default for SimOptions {
    fn default() -> Self {
        Self {
            probabilities: true,
        }
    }
}

impl SimOptions {
    pub(crate) fn classical_only() -> Self {
        Self {
            probabilities: false,
        }
    }
}

/// Result of a generic simulation run.
#[derive(Debug, Clone)]
pub struct RunOutcome {
    /// Classical measurement outcomes, indexed by classical bit number.
    /// `true` = measured |1⟩.
    pub classical_bits: Vec<bool>,
    /// Probability of each computational basis state (length 2^n).
    ///
    /// `None` means the selected backend cannot expose a dense probability
    /// distribution for this circuit. Other probability extraction failures
    /// are returned as errors by the query that produced this result.
    pub probabilities: Option<Probabilities>,
    /// Which engine ran, whether the answer is exact, and where the state lived.
    pub metadata: RunMetadata,
    /// What each [`Instruction::Save`] point recorded, in the order the points
    /// were reached. Empty unless the circuit carries save points.
    pub saves: Vec<SaveRecord>,
}

/// Frequency histogram returned by query-aware count sampling.
#[derive(Debug, Clone)]
pub struct CountsResult {
    /// Histogram keyed by packed classical bits; same key layout as
    /// [`ShotsResult::counts`], formattable with [`bitstring`].
    pub counts: HashMap<Vec<u64>, u64>,
    pub num_classical_bits: usize,
    pub metadata: RunMetadata,
}

impl CountsResult {
    pub fn into_counts(self) -> HashMap<Vec<u64>, u64> {
        self.counts
    }
}

/// Per-qubit marginal probabilities returned by query-aware marginal sampling.
#[derive(Debug, Clone)]
pub struct MarginalsResult {
    /// `(P(0), P(1))` per qubit, indexed by qubit number.
    pub marginals: Vec<(f64, f64)>,
    pub metadata: RunMetadata,
}

impl MarginalsResult {
    pub fn into_vec(self) -> Vec<(f64, f64)> {
        self.marginals
    }
}

/// Reduced density matrix of a qubit subset, returned by
/// [`Simulate::reduced_density_matrix`].
#[derive(Debug, Clone)]
pub struct ReducedDensityMatrix {
    /// The subsystem as it was requested, which fixes the index order below.
    pub qubits: Vec<usize>,
    /// Row major with side `2^k` over `k = qubits.len()`: `data[t * 2^k + t']`
    /// is `<t|rho|t'>`, where bit `i` of `t` is the state of `qubits[i]`, so
    /// `qubits[0]` is the lowest bit as `q[0]` is in a basis index. Trace one,
    /// Hermitian to rounding.
    pub data: Vec<Complex64>,
    pub metadata: RunMetadata,
}

impl ReducedDensityMatrix {
    /// `Tr(rho^2)`: 1 for a pure marginal, `2^-k` for the maximally mixed one.
    pub fn purity(&self) -> f64 {
        self.data.iter().map(|entry| entry.norm_sqr()).sum()
    }
}

/// Operator variance of a weighted Pauli observable, returned by
/// [`Simulate::observable_variance`].
#[derive(Debug, Clone)]
pub struct ObservableVariance {
    /// `<H^2> - <H>^2` on the output state.
    pub variance: f64,
    /// `<H>` on the same state, evaluated on the way to the variance.
    pub mean: f64,
    pub metadata: RunMetadata,
}

/// Entanglement entropy of a subsystem, returned by
/// [`Simulate::entanglement_entropy`].
#[derive(Debug, Clone)]
pub struct EntropyResult {
    /// The subsystem as it was requested, the side of the cut the entropy is
    /// read on.
    pub subsystem: Vec<usize>,
    /// Von Neumann entropy of the subsystem in nats: a Bell pair reads `ln 2`.
    pub entropy: f64,
    /// Schmidt values across the cut, descending, with squares summing to 1.
    /// `None` where the backend holds the entropy without the spectrum that
    /// stands behind it, as a stabilizer cut past the export cap does: its
    /// `2^r` equal values do not fit.
    pub schmidt_values: Option<Vec<f64>>,
    pub metadata: RunMetadata,
}

/// Overlap between the output states of two runs, returned by
/// [`Simulate::overlap`].
#[derive(Debug, Clone)]
pub struct OverlapResult {
    /// `|<a|b>|^2` over the two normalized states: 1 for the same state up to
    /// phase, 0 for orthogonal ones. The amplitude itself is not reported,
    /// since a tableau keeps no global phase and every MPS truncation moves
    /// one.
    pub fidelity: f64,
    /// Provenance of the run the terminal was called on, the left of the
    /// inner product.
    pub left: RunMetadata,
    /// Provenance of the run passed as the argument.
    pub right: RunMetadata,
}

/// Typestate marker: [`Simulate`] builder with no seed chosen yet.
#[derive(Debug, Clone, Copy)]
pub struct Unseeded;

/// Typestate marker: [`Simulate`] builder with its RNG seed fixed.
#[derive(Debug, Clone, Copy)]
pub struct Seeded {
    seed: u64,
}

/// A caller-supplied start: a `2^n` amplitude vector, or a `4^n` mixture in
/// the density-matrix layout.
#[derive(Debug, Clone, Copy)]
pub(crate) enum StartState<'c> {
    Amplitudes(&'c [Complex64]),
    DensityMatrix(&'c [Complex64]),
}

impl StartState<'_> {
    /// Load the start into `backend`, whose register width the buffer sets.
    fn load(self, backend: &mut dyn Backend, num_classical_bits: usize) -> Result<()> {
        match self {
            Self::Amplitudes(amplitudes) => {
                backend.init_from_amplitudes(amplitudes.to_vec(), num_classical_bits)
            }
            Self::DensityMatrix(rho) => {
                backend.init_from_density_matrix(rho.to_vec(), num_classical_bits)
            }
        }
    }
}

/// Builder for query-aware simulation requests.
pub struct Simulate<'c, SeedState> {
    circuit: &'c Circuit,
    kind: BackendKind,
    seed: SeedState,
    noise_model: Option<&'c noise::NoiseModel>,
    initial_state: Option<StartState<'c>>,
    require_exact: bool,
}

impl<'c, SeedState> Simulate<'c, SeedState> {
    /// Select an explicit backend kind instead of [`BackendKind::Auto`] routing.
    #[inline]
    pub fn backend(mut self, kind: BackendKind) -> Self {
        self.kind = kind;
        self
    }

    /// Reject a route that could return an approximate answer, rather than
    /// taking it and saying so in the result.
    ///
    /// [`BackendKind::Auto`] sends a circuit past the statevector cap to an MPS
    /// at a bounded bond dimension, which is the only route those circuits have;
    /// the result reports [`Exactness::Approximate`] either way. Call this when
    /// an approximate answer is worse than no answer, and the run returns
    /// `IncompatibleBackend` naming the engine it would have used.
    ///
    /// A circuit whose gates all preserve computational basis states is the
    /// exception: it stays on one basis state, so that route holds it at bond 1
    /// and answers exactly.
    ///
    /// Routes that can be decided from the circuit are rejected before any state
    /// is allocated; sparse Pauli dynamics only learns that it truncated while
    /// propagating, so that one is caught on the finished result instead.
    #[inline]
    pub fn require_exact(mut self) -> Self {
        self.require_exact = true;
        self
    }

    /// Attach a noise model.
    ///
    /// [`Simulate::shots`] and [`Simulate::sample_counts`] accept one on any
    /// backend with a per-shot pure state, averaging trajectories.
    /// [`Simulate::run`], [`Simulate::marginals`],
    /// [`Simulate::expectation_values`] and
    /// [`Simulate::reduced_density_matrix`] answer from the exact mixture
    /// instead, which only [`BackendKind::DensityMatrix`] and its device
    /// sibling hold, so they require one of those, as does
    /// [`Simulate::expectation_gradient_shift`] and, only to decline on it,
    /// [`Simulate::entanglement_entropy`]: a mixture has no Schmidt
    /// decomposition. [`Simulate::expectation_gradient`] rejects a noise model
    /// on every backend.
    #[inline]
    pub fn noise(mut self, model: &'c noise::NoiseModel) -> Self {
        self.noise_model = Some(model);
        self
    }

    /// Start from `amplitudes` instead of |0...0⟩.
    ///
    /// Indexed with qubit 0 in the least significant bit, length `2^n` for the
    /// circuit's `n` qubits, and normalized. A vector failing any of those is
    /// rejected with `InvalidParameter` before the run.
    ///
    /// A start state also constrains the route, because shape-based dispatch
    /// reads the circuit alone and its shortcuts hold only from |0...0⟩:
    /// [`BackendKind::Auto`] resolves to the statevector, and every backend
    /// other than the statevector (dense, device, or distributed) and
    /// [`BackendKind::DensityMatrix`] reports `IncompatibleBackend`.
    /// [`Simulate::expectation_gradient`] declines a start
    /// state, as do [`Simulate::shots`] and [`Simulate::sample_counts`] with a
    /// noise model attached, since trajectory replay has no start-state path.
    #[inline]
    pub fn initial_state(mut self, amplitudes: &'c [Complex64]) -> Self {
        self.initial_state = Some(StartState::Amplitudes(amplitudes));
        self
    }

    /// Start from the mixture `rho` instead of |0...0⟩.
    ///
    /// Row-major `2^n x 2^n` in the layout
    /// [`DensityMatrixBackend::density_matrix`] exports, qubit 0 in the least
    /// significant bit of both indices, so entry `(r << n) | c` is `<r|rho|c>`.
    /// The buffer must be Hermitian, each entry pair agreeing to 1e-12, and
    /// have unit trace to 1e-9; one failing either check, or of a length other
    /// than `4^n`, is rejected with `InvalidParameter` naming the check.
    /// Positive semidefiniteness is not checked, because the backend carries no
    /// eigendecomposition to check it with.
    ///
    /// Only [`BackendKind::DensityMatrix`] and its device sibling hold a
    /// mixture, so every other kind, [`BackendKind::Auto`] included, reports
    /// `IncompatibleBackend`. The terminals that decline a start state decline
    /// this one for the same reasons. Replaces an earlier
    /// [`Simulate::initial_state`], as that call replaces this one.
    ///
    /// [`DensityMatrixBackend::density_matrix`]: crate::backend::density_matrix::DensityMatrixBackend::density_matrix
    #[inline]
    pub fn initial_density_matrix(mut self, rho: &'c [Complex64]) -> Self {
        self.initial_state = Some(StartState::DensityMatrix(rho));
        self
    }

    /// Shortcut for [`Simulate::backend`] with [`BackendKind::StatevectorGpu`].
    #[cfg(feature = "gpu")]
    #[inline]
    pub fn gpu(self, context: std::sync::Arc<crate::gpu::GpuContext>) -> Self {
        self.backend(BackendKind::StatevectorGpu { context })
    }

    /// Automatic backend selection with GPU acceleration opted in via `context`.
    ///
    /// Routes like [`BackendKind::Auto`], but a selected statevector or
    /// stabilizer block that clears the qubit crossover with VRAM to spare runs
    /// on the device. Unsupported cases fall back to the identical CPU path.
    #[cfg(feature = "gpu")]
    #[inline]
    pub fn gpu_auto(self, context: std::sync::Arc<crate::gpu::GpuContext>) -> Self {
        self.backend(BackendKind::AutoGpu { context })
    }

    /// Distribute the exact state vector across the ranks of `context`.
    ///
    /// With a single rank this behaves like [`Simulate::backend`] with
    /// [`BackendKind::Statevector`].
    #[cfg(feature = "distributed")]
    pub fn distributed(
        self,
        context: std::sync::Arc<crate::distributed::DistributedContext>,
    ) -> Self {
        self.backend(BackendKind::StatevectorDistributed { context })
    }
}

impl<'c> Simulate<'c, Unseeded> {
    /// Query methods exist only on the seeded builder.
    #[inline]
    pub fn seed(self, seed: u64) -> Simulate<'c, Seeded> {
        Simulate {
            circuit: self.circuit,
            kind: self.kind,
            seed: Seeded { seed },
            noise_model: self.noise_model,
            initial_state: self.initial_state,
            require_exact: self.require_exact,
        }
    }
}

impl<'c> Simulate<'c, Seeded> {
    #[inline]
    fn seed_value(&self) -> u64 {
        self.seed.seed
    }

    /// Trajectory replay reinitializes a pure state per shot and the compiled
    /// noisy samplers are tableau based, so neither carries a start state.
    fn require_no_initial_state_under_noise(&self, terminal: &str) -> Result<()> {
        if self.initial_state.is_some() {
            return Err(reject_initial_state(
                &self.kind,
                terminal,
                "noisy trajectory replay starts every shot from |0...0>; read the exact mixture \
                 with `run`, `marginals`, or `expectation_values` on the density-matrix backend",
            ));
        }
        Ok(())
    }

    /// Execute the circuit once.
    ///
    /// With a noise model attached the probabilities are the exact noisy
    /// distribution rather than one trajectory, so the run needs the
    /// density-matrix backend; the classical bits are one draw, matching
    /// `shots(1)`. Readout error reaches the draw and not the state, so a
    /// model carrying it is rejected rather than answered with two fields from
    /// different distributions; `shots` and `sample_counts` apply it.
    #[inline]
    pub fn run(self) -> Result<RunOutcome> {
        let seed = self.seed_value();
        if self.require_exact {
            reject_approximate_route(&self.kind, self.circuit)?;
        }
        if let Some(noise_model) = self.noise_model {
            require_exact_mixture(&self.kind, "a single run")?;
            reject_readout_at(self.circuit, noise_model, "a single run")?;
            let probabilities = exact_noisy_probabilities(
                &self.kind,
                self.circuit,
                noise_model,
                self.initial_state,
                seed,
            )?;
            let classical_bits =
                sample_exact_noisy_shots(&probabilities, self.circuit, noise_model, 1, seed)
                    .swap_remove(0);
            let outcome = RunOutcome {
                classical_bits,
                probabilities: Some(probabilities),
                metadata: exact_mixture_metadata(&self.kind),
                saves: Vec::new(),
            };
            ensure_saves_recorded(self.circuit, &outcome)?;
            return Ok(outcome);
        }
        if let Some(state) = self.initial_state {
            return run_from_initial_state(
                &self.kind,
                self.circuit,
                state,
                seed,
                &SimOptions::default(),
            );
        }
        let outcome = run_with_internal(self.kind, self.circuit, seed, SimOptions::default())?;
        ensure_exact_result(self.require_exact, &outcome.metadata)?;
        Ok(outcome)
    }

    /// Refuse a circuit carrying save points on a terminal that cannot return
    /// them. Only [`Simulate::run`] has a place to put the records.
    fn reject_saves(&self, terminal: &str) -> Result<()> {
        if self.circuit.save_count() == 0 {
            return Ok(());
        }
        Err(PrismError::IncompatibleBackend {
            backend: format!("{:?}", self.kind),
            reason: format!(
                "the circuit carries {} save point(s), which only `run` returns, so \
                 `{terminal}` declines it",
                self.circuit.save_count()
            ),
        })
    }

    /// [`reject_approximate_route`] for the shot terminals, which without a noise
    /// model or a start state can sample through the stabilizer-rank expansion.
    fn reject_approximate_shot_route(&self) -> Result<()> {
        let sampling = self.noise_model.is_none() && self.initial_state.is_none();
        reject_predicted_route(approximate_route_name(&self.kind, self.circuit, sampling))
    }

    /// Execute `num_shots` times, collecting per-shot classical bits. Accepts
    /// an attached noise model.
    #[inline]
    pub fn shots(self, num_shots: usize) -> Result<ShotsResult> {
        self.reject_saves("shots")?;
        let seed = self.seed_value();
        let require_exact = self.require_exact;
        if require_exact {
            self.reject_approximate_shot_route()?;
        }
        let result = if let Some(noise_model) = self.noise_model {
            self.require_no_initial_state_under_noise("shot sampling")?;
            run_shots_with_noise(self.kind, self.circuit, noise_model, num_shots, seed)?
        } else if let Some(state) = self.initial_state {
            shots_from_initial_state(&self.kind, self.circuit, state, num_shots, seed)?
        } else {
            run_shots_with(self.kind, self.circuit, num_shots, seed)?
        };
        ensure_exact_result(require_exact, &result.metadata)?;
        Ok(result)
    }

    /// Sample a frequency histogram over `num_shots` executions. Accepts an
    /// attached noise model.
    ///
    /// Counts may be sampled directly from the output distribution, so seeded
    /// counts can differ from [`Simulate::shots`] plus [`ShotsResult::counts`]
    /// while drawing from the identical distribution.
    #[inline]
    pub fn sample_counts(self, num_shots: usize) -> Result<CountsResult> {
        self.reject_saves("sample_counts")?;
        let seed = self.seed_value();
        if self.require_exact {
            self.reject_approximate_shot_route()?;
        }
        let (counts, metadata) = if let Some(noise_model) = self.noise_model {
            self.require_no_initial_state_under_noise("count sampling")?;
            let shots =
                run_shots_with_noise(self.kind, self.circuit, noise_model, num_shots, seed)?;
            (shots.counts(), shots.metadata)
        } else if let Some(state) = self.initial_state {
            let shots = shots_from_initial_state(&self.kind, self.circuit, state, num_shots, seed)?;
            (shots.counts(), shots.metadata)
        } else {
            run_counts_with(self.kind, self.circuit, num_shots, seed)?
        };
        ensure_exact_result(self.require_exact, &metadata)?;
        Ok(CountsResult {
            counts,
            num_classical_bits: self.circuit.num_classical_bits,
            metadata,
        })
    }

    /// Per-qubit marginal probabilities as `(P(0), P(1))` pairs. Rejects
    /// backends without probability output, and with a noise model attached
    /// answers exactly from the mixture, which needs the density-matrix
    /// backend, and rejects a model carrying readout error, since
    /// `sample_counts` is the terminal that applies it.
    #[inline]
    pub fn marginals(self) -> Result<MarginalsResult> {
        self.reject_saves("marginals")?;
        let seed = self.seed_value();
        if self.require_exact {
            reject_approximate_route(&self.kind, self.circuit)?;
        }
        if let Some(noise_model) = self.noise_model {
            require_exact_mixture(&self.kind, "marginals")?;
            reject_readout_at(self.circuit, noise_model, "marginals")?;
            let probs = exact_noisy_probabilities(
                &self.kind,
                self.circuit,
                noise_model,
                self.initial_state,
                seed,
            )?;
            return Ok(MarginalsResult {
                marginals: probs.marginals(),
                metadata: exact_mixture_metadata(&self.kind),
            });
        }
        let result = if let Some(state) = self.initial_state {
            marginals_from_initial_state(&self.kind, self.circuit, state, seed)?
        } else {
            run_marginals_result_with(self.kind, self.circuit, seed)?
        };
        ensure_exact_result(self.require_exact, &result.metadata)?;
        Ok(result)
    }

    /// Compute `⟨ψ|P|ψ⟩` for each joint Pauli observable on the circuit's
    /// output state, honoring the selected backend.
    ///
    /// Each observable is a product of single-qubit Paulis (identity factors
    /// omitted). The circuit must be unitary. Clifford circuits propagate each
    /// observable exactly. Non-Clifford circuits use the state vector while they
    /// fit it; above that cap the selected backend evaluates the observable on
    /// its own representation, and a backend without one reports
    /// `BackendUnsupported` naming itself.
    ///
    /// With a noise model attached the value is the exact `Tr(rho P)` on the
    /// evolved mixture, which needs the density-matrix backend. A model
    /// carrying readout error is rejected: readout acts on the measurement
    /// record, which no observable sees, and `shots` on the same model would
    /// disagree by the readout rate.
    #[inline]
    pub fn expectation_values(self, observables: &[Vec<PauliTerm>]) -> Result<Vec<f64>> {
        self.reject_saves("expectation_values")?;
        self.expectation_values_reported(observables)
            .map(ExpectationResult::into_values)
    }

    /// [`Simulate::expectation_values`] with the provenance of the run and, for
    /// a route that estimates rather than evaluates, a standard error per value.
    pub fn expectation_values_reported(
        self,
        observables: &[Vec<PauliTerm>],
    ) -> Result<ExpectationResult> {
        let seed = self.seed_value();
        if self.require_exact {
            reject_approximate_route(&self.kind, self.circuit)?;
        }
        if let Some(noise_model) = self.noise_model {
            reject_readout_at(self.circuit, noise_model, "expectation values")?;
        }
        if let BackendKind::PauliPath { epsilon, max_terms } = self.kind {
            reject_pauli_path_initial_state(self.initial_state)?;
            return pauli_path_expectations(
                self.circuit,
                self.noise_model,
                observables,
                epsilon,
                max_terms,
            );
        }
        if let Some(noise_model) = self.noise_model {
            require_exact_mixture(&self.kind, "expectation values")?;
            require_unitary_circuit(&self.kind, self.circuit, "expectation values require")?;
            let values = noise::dm_expectation_values(
                &self.kind,
                self.circuit,
                observables,
                Some(noise_model),
                self.initial_state,
                seed,
            )?;
            return Ok(analytic_expectations(
                values,
                exact_mixture_metadata(&self.kind),
            ));
        }
        if let Some(state) = self.initial_state {
            require_unitary_circuit(&self.kind, self.circuit, "expectation values require")?;
            return expectation_values_from_initial_state(
                &self.kind,
                self.circuit,
                state,
                observables,
                seed,
            );
        }
        let result = run_expectation_values_reported(self.kind, self.circuit, observables, seed)?;
        ensure_exact_result(self.require_exact, &result.metadata)?;
        Ok(result)
    }

    /// Compute `⟨H⟩` and its grouped-measurement variance for a weighted
    /// Pauli observable on the circuit's output state.
    ///
    /// The statevector family evaluates one traversal per qubit-wise-commuting
    /// group and reports the variance; see [`ObservableExpectation::variance`]
    /// for what the number means. Every other route, including runs with a
    /// noise model or start state attached, evaluates term by term through
    /// [`Simulate::expectation_values`] semantics and reports the weighted
    /// mean with no variance. A noise model carrying readout error is rejected
    /// for the same reason as there.
    pub fn observable_expectation(
        self,
        observable: &PauliObservable,
    ) -> Result<ObservableExpectation> {
        self.reject_saves("observable_expectation")?;
        self.observable_expectation_ref(observable)
    }

    /// [`Simulate::observable_expectation`] without consuming the builder, so
    /// [`Simulate::observable_variance`] can evaluate `H` and `H^2` on one
    /// request.
    fn observable_expectation_ref(
        &self,
        observable: &PauliObservable,
    ) -> Result<ObservableExpectation> {
        let seed = self.seed_value();
        if self.require_exact {
            reject_approximate_route(&self.kind, self.circuit)?;
        }
        if let Some(noise_model) = self.noise_model {
            reject_readout_at(self.circuit, noise_model, "observable expectation")?;
        }
        if let BackendKind::PauliPath { epsilon, max_terms } = self.kind {
            reject_pauli_path_initial_state(self.initial_state)?;
            let result = pauli_path_expectations(
                self.circuit,
                self.noise_model,
                &observable_vecs(observable),
                epsilon,
                max_terms,
            )?;
            let metadata = result.metadata;
            return Ok(weighted_observable_result(
                observable,
                &result.values,
                None,
                metadata,
            ));
        }
        if let Some(noise_model) = self.noise_model {
            require_exact_mixture(&self.kind, "expectation values")?;
            require_unitary_circuit(&self.kind, self.circuit, "expectation values require")?;
            let values = noise::dm_expectation_values(
                &self.kind,
                self.circuit,
                &observable_vecs(observable),
                Some(noise_model),
                self.initial_state,
                seed,
            )?;
            return Ok(weighted_observable_result(
                observable,
                &values,
                None,
                exact_mixture_metadata(&self.kind),
            ));
        }
        if let Some(state) = self.initial_state {
            require_unitary_circuit(&self.kind, self.circuit, "expectation values require")?;
            let result = expectation_values_from_initial_state(
                &self.kind,
                self.circuit,
                state,
                &observable_vecs(observable),
                seed,
            )?;
            return Ok(weighted_observable_result(
                observable,
                &result.values,
                result.std_errors.as_deref(),
                result.metadata,
            ));
        }
        let result =
            run_observable_expectation_reported(self.kind.clone(), self.circuit, observable, seed)?;
        ensure_exact_result(self.require_exact, &result.metadata)?;
        Ok(result)
    }

    /// `Var(H) = <H^2> - <H>^2` for a weighted Pauli observable on the
    /// circuit's output state.
    ///
    /// This is the spread of the operator itself, the number a shot-based
    /// estimate of `<H>` converges on dividing by the shot count. It is not
    /// [`ObservableExpectation::variance`], which sums per-group variances and
    /// so drops the covariance between measurement groups.
    ///
    /// Evaluates `H` and the square of its traceless part through
    /// [`Simulate::observable_expectation`], so backend routing, noise, and
    /// start states behave as they do there. The constant term is held out of
    /// the square rather than cancelled inside it; see
    /// [`PauliObservable::split_identity`]. The square carries up to `T^2`
    /// terms over `H`'s `T`; see [`PauliObservable::square`].
    pub fn observable_variance(self, observable: &PauliObservable) -> Result<ObservableVariance> {
        self.reject_saves("observable_variance")?;
        let (offset, traceless) = observable.split_identity();
        let mean = self.observable_expectation_ref(observable)?;
        let second = self.observable_expectation_ref(&traceless.square())?;
        let centered = mean.mean - offset;
        Ok(ObservableVariance {
            variance: second.mean - centered * centered,
            mean: mean.mean,
            metadata: mean.metadata,
        })
    }

    /// Joint probability distribution over `qubits`, `2^k` entries with
    /// `qubits[0]` in the lowest bit.
    ///
    /// The subset generalizes [`Simulate::marginals`], which reports each
    /// qubit on its own and so cannot show correlation: a Bell pair reads
    /// `(0.5, 0.5)` twice there and `[0.5, 0, 0, 0.5]` here. Routing follows
    /// [`Simulate::run`], including the exact mixture a noise model asks for
    /// and its rejection of readout error, which acts on the measurement
    /// record rather than on the state.
    ///
    /// A backend that exposes no distribution for the circuit reports
    /// `BackendUnsupported` naming itself.
    pub fn probabilities_of(self, qubits: &[usize]) -> Result<Vec<f64>> {
        self.reject_saves("probabilities_of")?;
        crate::backend::schmidt::validate_qubit_set(qubits, self.circuit.num_qubits)?;
        let kind = format!("{:?}", self.kind);
        let outcome = self.run()?;
        let probabilities = outcome
            .probabilities
            .ok_or(PrismError::BackendUnsupported {
                backend: kind,
                operation: "a probability distribution to marginalize".into(),
            })?;
        Ok(probabilities.subset_marginal(qubits))
    }

    /// Full amplitude vector of the circuit's output state, honoring the
    /// selected backend.
    ///
    /// Indexed with qubit 0 in the least significant bit, so `x q[0]` puts the
    /// amplitude at index 1. The circuit must be unitary, for the reason
    /// [`Simulate::reduced_density_matrix`] gives.
    ///
    /// A noise model declines: a mixture has no single amplitude vector, and
    /// [`Simulate::reduced_density_matrix`] over the whole register is the
    /// terminal that answers there. The density-matrix backend declines for
    /// the same reason whether or not noise is attached.
    ///
    /// The vector holds `2^n` amplitudes, so a register past the dense export
    /// cap reports `IncompatibleBackend` before allocating rather than after.
    pub fn state_vector(self) -> Result<Vec<Complex64>> {
        self.reject_saves("state_vector")?;
        let seed = self.seed_value();
        let diagnostic = Diagnostic::StateVector;
        require_unitary_circuit(&self.kind, self.circuit, "a statevector requires")?;
        if self.require_exact {
            reject_approximate_route(&self.kind, self.circuit)?;
        }
        if self.noise_model.is_some() {
            return Err(PrismError::IncompatibleBackend {
                backend: format!("{:?}", self.kind),
                reason: format!(
                    "{} is a pure state; a noise model evolves a mixture, which \
                     `reduced_density_matrix` over the whole register reports",
                    diagnostic.terminal()
                ),
            });
        }
        let backend = diagnostic_backend(
            &self.kind,
            self.circuit,
            self.initial_state,
            seed,
            diagnostic,
            self.circuit.num_qubits,
        )?;
        ensure_exact_result(self.require_exact, &backend_metadata(&*backend))?;
        backend.export_statevector()
    }

    /// Reduced density matrix of `qubits` on the circuit's output state,
    /// honoring the selected backend.
    ///
    /// Row major with side `2^k`; [`ReducedDensityMatrix::data`] states the
    /// index order. The subsystem is named once and may be the whole register.
    /// The circuit must be unitary: the answer is read off one state, and a
    /// measurement, reset or conditional leaves one seeded branch of several.
    ///
    /// An explicitly selected backend that holds no partial trace reports
    /// `BackendUnsupported` naming itself. Under [`BackendKind::Auto`] a route
    /// that cannot answer falls back to the statevector while the circuit fits
    /// its cap, so the diagnostic is served rather than declined by a choice
    /// the caller did not make. With a noise model attached the answer is the
    /// marginal of the exact mixture, which needs the density-matrix backend.
    pub fn reduced_density_matrix(self, qubits: &[usize]) -> Result<ReducedDensityMatrix> {
        self.reject_saves("reduced_density_matrix")?;
        let seed = self.seed_value();
        let diagnostic = Diagnostic::ReducedDensityMatrix;
        let terminal = diagnostic.terminal();
        crate::backend::schmidt::validate_qubit_set(qubits, self.circuit.num_qubits)?;
        require_unitary_circuit(
            &self.kind,
            self.circuit,
            "a reduced density matrix requires",
        )?;
        if self.require_exact {
            reject_approximate_route(&self.kind, self.circuit)?;
        }
        if let Some(noise_model) = self.noise_model {
            reject_readout_at(self.circuit, noise_model, terminal)?;
            require_exact_mixture(&self.kind, terminal)?;
            let mut mixture = noise::evolve_density_matrix(
                &self.kind,
                self.circuit,
                Some(noise_model),
                self.initial_state,
                seed,
            )?;
            return Ok(ReducedDensityMatrix {
                qubits: qubits.to_vec(),
                data: mixture.reduced_density_matrix(qubits)?,
                metadata: exact_mixture_metadata(&self.kind),
            });
        }
        let mut backend = diagnostic_backend(
            &self.kind,
            self.circuit,
            self.initial_state,
            seed,
            diagnostic,
            qubits.len(),
        )?;
        let metadata = backend_metadata(&*backend);
        ensure_exact_result(self.require_exact, &metadata)?;
        Ok(ReducedDensityMatrix {
            qubits: qubits.to_vec(),
            data: backend.reduced_density_matrix(qubits)?,
            metadata,
        })
    }

    /// Entanglement entropy of `subsystem` across its cut with the rest of the
    /// register, in nats, honoring the selected backend.
    ///
    /// `subsystem` must leave both sides of the cut non-empty, and the circuit
    /// must be unitary, for the reason [`Simulate::reduced_density_matrix`]
    /// gives. The Schmidt values come back with the entropy, descending and
    /// normalized.
    ///
    /// A backend that holds the entropy without the spectrum behind it, a
    /// stabilizer cut past the export cap, answers with
    /// [`EntropyResult::schmidt_values`] set to `None`.
    /// An explicitly selected backend that holds neither reports
    /// `BackendUnsupported` naming itself; under [`BackendKind::Auto`] such a
    /// route falls back to the statevector while the circuit fits its cap. A
    /// noise model declines outright: it sends the run to the density matrix,
    /// whose mixed state has no Schmidt decomposition.
    pub fn entanglement_entropy(self, subsystem: &[usize]) -> Result<EntropyResult> {
        self.reject_saves("entanglement_entropy")?;
        let seed = self.seed_value();
        let diagnostic = Diagnostic::Entropy;
        let terminal = diagnostic.terminal();
        crate::backend::schmidt::validate_subsystem(subsystem, self.circuit.num_qubits)?;
        require_unitary_circuit(&self.kind, self.circuit, "entanglement entropy requires")?;
        if self.require_exact {
            reject_approximate_route(&self.kind, self.circuit)?;
        }
        if let Some(noise_model) = self.noise_model {
            reject_readout_at(self.circuit, noise_model, terminal)?;
            require_exact_mixture(&self.kind, terminal)?;
            let mut mixture = noise::evolve_density_matrix(
                &self.kind,
                self.circuit,
                Some(noise_model),
                self.initial_state,
                seed,
            )?;
            return Ok(EntropyResult {
                subsystem: subsystem.to_vec(),
                entropy: mixture.entanglement_entropy(subsystem)?,
                schmidt_values: None,
                metadata: exact_mixture_metadata(&self.kind),
            });
        }
        let mut backend = diagnostic_backend(
            &self.kind,
            self.circuit,
            self.initial_state,
            seed,
            diagnostic,
            subsystem.len(),
        )?;
        let metadata = backend_metadata(&*backend);
        ensure_exact_result(self.require_exact, &metadata)?;
        // A tableau past the export cap holds the entropy as a rank while its
        // `2^r` equal values do not fit, so the spectrum declining is not the
        // terminal declining. A backend holding neither reports the spectrum's
        // own error, which is what its entropy raises too.
        let (entropy, schmidt_values) = match backend.schmidt_values(subsystem) {
            Ok(values) => (
                crate::backend::schmidt::entropy_of_schmidt_values(&values),
                Some(values),
            ),
            Err(declined) => match backend.entanglement_entropy(subsystem) {
                Ok(entropy) => (entropy, None),
                Err(_) => return Err(declined),
            },
        };
        Ok(EntropyResult {
            subsystem: subsystem.to_vec(),
            entropy,
            schmidt_values,
            metadata,
        })
    }

    /// `|<a|b>|^2` between this circuit's output state and `other`'s, honoring
    /// the backend each side selected.
    ///
    /// The two circuits must declare the same width, and both must be unitary
    /// for the reason [`Simulate::reduced_density_matrix`] gives. Each side
    /// carries its own backend, seed and start state, and each resolves to a
    /// single backend rather than the decomposed route, since two circuits
    /// need not split into the same independent blocks.
    ///
    /// A pair of unlike representations is served by a dense export of both
    /// states, so it reaches as far as the export cap does. A pair that shares
    /// one answers natively at any width: two chains in the same site order,
    /// two tableaux, two product states, or two sparse maps. A noise model on
    /// either side is rejected, since the fidelity of two mixtures is not an
    /// inner product.
    pub fn overlap(self, other: Simulate<'_, Seeded>) -> Result<OverlapResult> {
        self.reject_saves("overlap")?;
        let diagnostic = Diagnostic::Overlap;
        if self.circuit.num_qubits != other.circuit.num_qubits {
            return Err(PrismError::InvalidParameter {
                message: format!(
                    "{} needs two circuits of the same width; got {} and {} qubits",
                    diagnostic.terminal(),
                    self.circuit.num_qubits,
                    other.circuit.num_qubits
                ),
            });
        }
        let (left_backend, left) = self.overlap_side(diagnostic)?;
        let (right_backend, right) = other.overlap_side(diagnostic)?;
        Ok(OverlapResult {
            fidelity: left_backend.overlap_sq(&*right_backend)?,
            left,
            right,
        })
    }

    /// One side of [`Simulate::overlap`]: the guards both sides answer to,
    /// then the backend this side's kind resolves to, run and handed back with
    /// its provenance.
    fn overlap_side(self, diagnostic: Diagnostic) -> Result<(Box<dyn Backend>, RunMetadata)> {
        let seed = self.seed_value();
        let terminal = diagnostic.terminal();
        require_unitary_circuit(&self.kind, self.circuit, "a state overlap requires")?;
        if self.require_exact {
            reject_approximate_route(&self.kind, self.circuit)?;
        }
        if self.noise_model.is_some() {
            return Err(PrismError::IncompatibleBackend {
                backend: format!("{:?}", self.kind),
                reason: format!(
                    "{terminal} is an inner product of two pure states, and a noise model \
                     evolves a mixture, whose fidelity is a different computation; drop \
                     the model, or compare the mixtures through `reduced_density_matrix`"
                ),
            });
        }
        let backend = diagnostic_backend(
            &self.kind,
            self.circuit,
            self.initial_state,
            seed,
            diagnostic,
            0,
        )?;
        let metadata = backend_metadata(&*backend);
        ensure_exact_result(self.require_exact, &metadata)?;
        Ok((backend, metadata))
    }

    /// Compute `⟨H⟩` and its exact gradient with respect to the bound
    /// parameters using the adjoint method.
    ///
    /// `hamiltonian` is a weighted Pauli sum `Σ c_k P_k` with real
    /// coefficients. `params` declares which gate instructions carry parameters.
    /// Runs on the statevector backend; the selected backend must be `Auto` or
    /// `Statevector`. The circuit must be unitary. See
    /// [`gradient::run_expectation_gradient`].
    #[inline]
    pub fn expectation_gradient(
        self,
        hamiltonian: &[(f64, Vec<PauliTerm>)],
        params: &crate::circuit::Parameters,
    ) -> Result<gradient::ExpectationGradient> {
        self.reject_saves("expectation_gradient")?;
        let seed = self.seed_value();
        if self.require_exact {
            reject_approximate_route(&self.kind, self.circuit)?;
        }
        if self.noise_model.is_some() {
            return Err(PrismError::IncompatibleBackend {
                backend: format!("{:?}", self.kind),
                reason: "the adjoint method backpropagates through a pure state, so no backend \
                         has a noisy adjoint path; drop the noise model, or take \
                         `expectation_gradient_shift` on the density-matrix backend"
                    .into(),
            });
        }
        if self.initial_state.is_some() {
            return Err(reject_initial_state(
                &self.kind,
                "the adjoint gradient",
                "the backward pass reconstructs the input register by inverting the circuit from \
                 |0...0>, so a start state would have to be inverted with it",
            ));
        }
        if !(self.kind.is_auto() || matches!(self.kind, BackendKind::Statevector)) {
            return Err(PrismError::IncompatibleBackend {
                backend: format!("{:?}", self.kind),
                reason:
                    "adjoint gradients run on the statevector backend; select Auto or Statevector"
                        .into(),
            });
        }
        gradient::run_expectation_gradient(self.circuit, hamiltonian, params, seed)
    }

    /// Compute `⟨H⟩` and its gradient by the parameter-shift rule on the
    /// selected backend.
    ///
    /// Serves the cases [`Simulate::expectation_gradient`] declines: any
    /// backend with a native observable path, circuits containing `QftBlock`,
    /// widths past the statevector cap, and a noise model. It differentiates the
    /// same gate set (`Rx`, `Ry`, `Rz`, `Rzz`, `P`, `PauliRot`) at `1 + 2 * links`
    /// circuit evaluations against the adjoint's one, so the adjoint stays the
    /// better choice where it applies. A backend with no native observable path
    /// reports `BackendUnsupported` naming itself. Under a noise model every
    /// evaluation reads the exact mixture, so the backend must be
    /// [`BackendKind::DensityMatrix`] or its device sibling; the shift stays
    /// exact because the channels do not depend on the shifted angle. See
    /// [`gradient::run_expectation_gradient_shift`].
    #[inline]
    pub fn expectation_gradient_shift(
        self,
        hamiltonian: &[(f64, Vec<PauliTerm>)],
        params: &crate::circuit::Parameters,
    ) -> Result<gradient::ExpectationGradient> {
        self.reject_saves("expectation_gradient_shift")?;
        let seed = self.seed_value();
        if self.require_exact {
            reject_approximate_route(&self.kind, self.circuit)?;
        }
        if self.noise_model.is_some()
            && !self.kind.is_density_matrix()
            && !matches!(self.kind, BackendKind::PauliPath { .. })
        {
            return Err(PrismError::IncompatibleBackend {
                backend: format!("{:?}", self.kind),
                reason: "a parameter-shift gradient under a noise model evaluates the exact \
                         mixed state, which only the density-matrix backend holds; select \
                         `BackendKind::DensityMatrix` or its device sibling, or drop the noise \
                         model"
                    .into(),
            });
        }
        gradient::shift_gradient(
            &self.kind,
            self.circuit,
            hamiltonian,
            params,
            self.noise_model,
            self.initial_state,
            seed,
        )
    }
}

/// Start a query-aware simulation request for `circuit`.
///
/// The returned builder defaults to automatic backend selection; chain
/// [`Simulate::seed`] to unlock the query methods.
///
/// # Examples
///
/// ```
/// use prism_q::{Circuit, Gate, simulate};
///
/// let mut circuit = Circuit::new(2, 0);
/// circuit.add_gate(Gate::H, &[0]);
/// circuit.add_gate(Gate::Cx, &[0, 1]);
///
/// let result = simulate(&circuit).seed(42).run()?;
/// let probs = result.probabilities.expect("no probabilities").to_vec();
/// // Bell state: ~50% |00>, ~50% |11>
/// assert!((probs[0] - 0.5).abs() < 1e-10);
/// assert!((probs[3] - 0.5).abs() < 1e-10);
/// # Ok::<(), prism_q::PrismError>(())
/// ```
#[inline]
pub fn simulate(circuit: &Circuit) -> Simulate<'_, Unseeded> {
    Simulate {
        circuit,
        kind: BackendKind::Auto,
        seed: Unseeded,
        noise_model: None,
        initial_state: None,
        require_exact: false,
    }
}

/// Rejection naming the two terminals the Pauli path engine serves.
fn reject_pauli_path(terminal: &str) -> PrismError {
    PrismError::IncompatibleBackend {
        backend: "PauliPath".into(),
        reason: format!(
            "{terminal} is not served by Pauli path propagation, which answers \
             `expectation_values` and `observable_expectation` only"
        ),
    }
}

fn reject_pauli_path_initial_state(state: Option<StartState<'_>>) -> Result<()> {
    match state {
        Some(_) => Err(reject_pauli_path("a start state")),
        None => Ok(()),
    }
}

/// One Pauli path evaluation per observable, under `noise` when one is
/// attached and on the bare circuit otherwise.
///
/// The route is exact only when nothing was truncated, so exactness is decided
/// by the discarded mass the run reports rather than by the parameters it was
/// given: a budget that never binds still returns an exact value.
fn pauli_path_expectations(
    circuit: &Circuit,
    noise: Option<&NoiseModel>,
    observables: &[Vec<PauliTerm>],
    epsilon: f64,
    max_terms: usize,
) -> Result<ExpectationResult> {
    let empty;
    let noise = match noise {
        Some(model) => model,
        None => {
            empty = NoiseModel {
                after_gate: vec![Vec::new(); circuit.instructions.len()],
                readout: vec![None; circuit.num_classical_bits],
            };
            &empty
        }
    };
    let mut values = Vec::with_capacity(observables.len());
    let mut discarded = 0.0f64;
    for obs in observables {
        let result =
            unified_pauli::run_pauli_path_observable(circuit, noise, obs, epsilon, max_terms)?;
        discarded = discarded.max(result.total_discarded);
        values.push(result.mean);
    }
    let metadata = if discarded > 0.0 {
        RunMetadata::approximate(ResolvedBackend::PauliPath)
    } else {
        RunMetadata::exact(ResolvedBackend::PauliPath)
    };
    Ok(analytic_expectations(values, metadata))
}

fn require_exact_mixture(kind: &BackendKind, terminal: &str) -> Result<()> {
    if kind.is_density_matrix() {
        return Ok(());
    }
    Err(PrismError::IncompatibleBackend {
        backend: format!("{kind:?}"),
        reason: format!(
            "{terminal} under a noise model reads the exact mixed state, which only the \
             density-matrix backend holds; select it, or average trajectories through `shots` \
             or `sample_counts`"
        ),
    })
}

/// Gate for the terminals that answer from the state, mixed or pure, which
/// readout error is not part of.
///
/// The condition is the one under which a draw would actually differ. Entries
/// past the circuit's bit count never apply, matching what
/// `sample_exact_noisy_shots` consumes, so a model built for a wider circuit is
/// not rejected on bits this one does not have, and a zero-rate entry flips
/// nothing, so a sweep that starts at zero is not rejected at its first point.
fn reject_readout_at(
    circuit: &Circuit,
    noise_model: &noise::NoiseModel,
    terminal: &str,
) -> Result<()> {
    let inert = noise_model
        .readout
        .iter()
        .take(circuit.num_classical_bits)
        .all(|entry| entry.as_ref().is_none_or(|readout| readout.is_inert()));
    if inert {
        return Ok(());
    }
    Err(PrismError::InvalidParameter {
        message: format!(
            "{terminal} answers from the mixed state, which readout error is not part of: it \
             acts on the measurement record and is indexed by classical bit, not qubit. Drop \
             it from the model, or use `shots` or `sample_counts`, which apply it"
        ),
    })
}

/// Gate for the terminals that cannot carry a start state through their own
/// machinery.
fn reject_initial_state(kind: &BackendKind, terminal: &str, instead: &str) -> PrismError {
    PrismError::IncompatibleBackend {
        backend: format!("{kind:?}"),
        reason: format!("{terminal} does not accept a start state; {instead}"),
    }
}

/// Reject a start state whose width disagrees with the circuit's.
///
/// The circuit's declared register wins: the buffer sets the backend's width,
/// so a shorter or longer one would silently simulate a different register
/// than the one the instructions index. A register too wide to index with a
/// `usize` has no dense start state at all.
fn check_initial_state_len(state: StartState<'_>, num_qubits: usize) -> Result<()> {
    let (buffer, bits, entries) = match state {
        StartState::Amplitudes(amplitudes) => (amplitudes, num_qubits, "amplitudes"),
        StartState::DensityMatrix(rho) => (rho, 2 * num_qubits, "density matrix entries"),
    };
    let want = (bits < usize::BITS as usize).then(|| 1usize << bits);
    if want == Some(buffer.len()) {
        return Ok(());
    }
    let needs = match want {
        Some(count) => count.to_string(),
        None => format!("2^{bits}"),
    };
    Err(PrismError::InvalidParameter {
        message: format!(
            "start state has {} {entries}, but a {num_qubits}-qubit circuit needs {needs}",
            buffer.len()
        ),
    })
}

/// Build the constrained backend for a start state and load it.
fn backend_from_initial_state(
    kind: &BackendKind,
    circuit: &Circuit,
    state: StartState<'_>,
    seed: u64,
) -> Result<Box<dyn Backend>> {
    if !kind.is_auto() {
        validate_explicit_backend(kind, circuit)?;
    }
    check_initial_state_len(state, circuit.num_qubits)?;
    let mut backend = initial_state_plan(kind, circuit.num_qubits, state)?.build(seed);
    state.load(&mut *backend, circuit.num_classical_bits)?;
    Ok(backend)
}

/// Expand the gate forms `backend` has no native kernel for, `QftBlock` and
/// `PauliRot`, or on a backend with the `QftBlock` kernel fold expanded QFTs into
/// blocks. The stream stays borrowed when nothing changes.
fn expand_for_backend<'c>(
    backend: &dyn Backend,
    circuit: &'c Circuit,
) -> std::borrow::Cow<'c, Circuit> {
    use std::borrow::Cow;
    let expanded = if backend.supports_qft_block() {
        crate::circuit::recognize_qft_blocks(circuit)
    } else {
        crate::circuit::expand_qft_blocks(circuit)
    };
    if backend.supports_pauli_rotation() {
        return expanded;
    }
    match expanded {
        Cow::Borrowed(borrowed) => crate::circuit::expand_pauli_rotations(borrowed),
        Cow::Owned(owned) => {
            let rotations = crate::circuit::expand_pauli_rotations(&owned);
            if let Cow::Owned(expanded_rotations) = rotations {
                return Cow::Owned(expanded_rotations);
            }
            Cow::Owned(owned)
        }
    }
}

/// Fuse `circuit` against what `backend` accepts and how wide a buffer it
/// sweeps, the two backend facts the pass pipeline is gated on.
fn fuse_for_backend<'a>(
    backend: &dyn Backend,
    circuit: &'a Circuit,
) -> std::borrow::Cow<'a, Circuit> {
    crate::circuit::fusion::fuse_circuit_for_width(
        circuit,
        backend.supports_fused_gates(),
        backend.fusion_state_qubits(circuit.num_qubits),
    )
}

/// [`apply_fused_circuit`] for a terminal that has nowhere to put a save.
///
/// Only [`Simulate::run`] returns save records, so every other terminal
/// declines a circuit carrying a save point rather than running it and
/// discarding what the points recorded.
fn apply_fused_without_saves(
    backend: &mut dyn Backend,
    circuit: &Circuit,
    terminal: &str,
) -> Result<()> {
    if circuit.save_count() > 0 {
        return Err(PrismError::IncompatibleBackend {
            backend: backend.name().to_string(),
            reason: format!(
                "save points are returned by `run`, and `{terminal}` has nowhere to put them"
            ),
        });
    }
    apply_fused_circuit(backend, circuit).map(|_| ())
}

/// Fuse `circuit` for `backend` and apply it, leaving initialization to the
/// caller. The start-state analogue of [`execute`], which owns the |0...0⟩ init.
fn apply_fused_circuit(backend: &mut dyn Backend, circuit: &Circuit) -> Result<Vec<SaveRecord>> {
    let expanded = expand_for_backend(&*backend, circuit);
    let fused = fuse_for_backend(&*backend, &expanded);
    apply_recording_saves(backend, &fused.instructions)
}

fn run_from_initial_state(
    kind: &BackendKind,
    circuit: &Circuit,
    state: StartState<'_>,
    seed: u64,
    opts: &SimOptions,
) -> Result<RunOutcome> {
    let mut backend = backend_from_initial_state(kind, circuit, state, seed)?;
    let saves = apply_fused_circuit(&mut *backend, circuit)?;

    let probabilities = if opts.probabilities {
        try_backend_probabilities(&*backend)?
    } else {
        None
    };
    Ok(RunOutcome {
        classical_bits: backend.classical_results().to_vec(),
        probabilities,
        metadata: backend_metadata(&*backend),
        saves,
    })
}

/// Shots for a start state. Terminal measurements sample one evolved
/// distribution; anything else replays the start state per shot, which is what
/// mid-circuit collapse and classical feedback need.
fn shots_from_initial_state(
    kind: &BackendKind,
    circuit: &Circuit,
    state: StartState<'_>,
    num_shots: usize,
    seed: u64,
) -> Result<ShotsResult> {
    let bits = circuit.num_classical_bits;
    if circuit.has_terminal_measurements_only() && !circuit.has_resets() {
        let stripped = circuit.without_measurements();
        let mut backend = backend_from_initial_state(kind, &stripped, state, seed)?;
        apply_fused_circuit(&mut *backend, &stripped)?;
        let meas_map = circuit.measurement_map();
        let shots = if backend.supports_native_sampling() {
            let samples = backend.sample_basis_states(num_shots, seed)?;
            Some(shots_from_basis_samples(&samples, &meas_map, bits))
        } else {
            try_backend_probabilities(&*backend)?
                .map(|probs| sample_shots(&probs, &meas_map, bits, num_shots, seed))
        };
        if let Some(shots) = shots {
            return Ok(
                ShotsResult::from_shots(shots, bits).with_metadata(backend_metadata(&*backend))
            );
        }
    }

    // Plan, expansion, and fusion are shot independent, so they are hoisted out
    // of the replay loop, matching `run_shots_per_shot`.
    if !kind.is_auto() {
        validate_explicit_backend(kind, circuit)?;
    }
    check_initial_state_len(state, circuit.num_qubits)?;
    let plan = initial_state_plan(kind, circuit.num_qubits, state)?;
    let probe = plan.build(seed);
    let expanded = expand_for_backend(&*probe, circuit);
    let fused = fuse_for_backend(&*probe, &expanded);

    let route = plan.resolved();
    let states = [(route, circuit.num_qubits)];
    collect_shots(
        circuit,
        num_shots,
        seed,
        route,
        kind,
        &states,
        |shot_seed| {
            let mut backend = plan.build(shot_seed);
            state.load(&mut *backend, circuit.num_classical_bits)?;
            backend.apply_instructions(&fused.instructions)?;
            Ok((
                backend.classical_results().to_vec(),
                backend_metadata(&*backend),
            ))
        },
    )
}

fn marginals_from_initial_state(
    kind: &BackendKind,
    circuit: &Circuit,
    state: StartState<'_>,
    seed: u64,
) -> Result<MarginalsResult> {
    let mut backend = backend_from_initial_state(kind, circuit, state, seed)?;
    apply_fused_without_saves(&mut *backend, circuit, "marginals")?;
    if backend.supports_pauli_expectation() {
        return marginals_from_pauli_expectations(&*backend, circuit.num_qubits);
    }
    Ok(MarginalsResult {
        marginals: Probabilities::Dense(backend.probabilities()?).marginals(),
        metadata: backend_metadata(&*backend),
    })
}

fn expectation_values_from_initial_state(
    kind: &BackendKind,
    circuit: &Circuit,
    state: StartState<'_>,
    observables: &[Vec<PauliTerm>],
    seed: u64,
) -> Result<ExpectationResult> {
    // Before the run, matching the dense statevector path.
    let masks = observables
        .iter()
        .map(|obs| pauli_masks(obs, circuit.num_qubits))
        .collect::<Result<Vec<_>>>()?;

    let mut backend = backend_from_initial_state(kind, circuit, state, seed)?;
    apply_fused_without_saves(&mut *backend, circuit, "expectation_values")?;
    let metadata = backend_metadata(&*backend);
    if backend.supports_pauli_expectation() {
        let values = backend.pauli_expectations(observables)?;
        return Ok(analytic_expectations(values, metadata));
    }

    let evolved = backend.export_statevector()?;
    let (values, _) = pauli_expectations_from_masks(&evolved, &masks);
    Ok(analytic_expectations(values, metadata))
}

/// Reject a result that turned out approximate, for a caller that opted out.
///
/// The backstop behind [`reject_approximate_route`], which decides from the
/// circuit alone. Sparse Pauli dynamics truncates on coefficient magnitudes it
/// only learns while propagating, so whether it stayed exact is not knowable
/// before the run.
fn ensure_exact_result(require_exact: bool, metadata: &RunMetadata) -> Result<()> {
    if require_exact && !metadata.is_exact() {
        return Err(PrismError::IncompatibleBackend {
            backend: format!("{:?}", metadata.backend),
            reason: "require_exact rejects an approximate result; this engine discarded state                      weight while running, which the route could not predict"
                .into(),
        });
    }
    Ok(())
}

/// Reject an approximate route for a caller that opted out of one.
fn reject_approximate_route(kind: &BackendKind, circuit: &Circuit) -> Result<()> {
    reject_predicted_route(approximate_route_name(kind, circuit, false))
}

/// The engine is named in the error, since a caller who asked for exactness
/// wants to know which one would have answered.
fn reject_predicted_route(engine: Option<&'static str>) -> Result<()> {
    match engine {
        Some(engine) => Err(PrismError::IncompatibleBackend {
            backend: engine.into(),
            reason: "require_exact rejects a route that can discard state weight; drop the \
                     requirement to accept the approximation, which the result reports, or \
                     select a backend that represents this circuit exactly"
                .into(),
        }),
        None => Ok(()),
    }
}

/// `subject` names what needs the unitary circuit with its verb, so the
/// sentence reads for a plural terminal ("expectation values require") and for
/// a singular one ("a reduced density matrix requires").
fn require_unitary_circuit(kind: &BackendKind, circuit: &Circuit, subject: &str) -> Result<()> {
    if has_nonunitary_or_classical_ops(circuit) {
        return Err(PrismError::IncompatibleBackend {
            backend: format!("{kind:?}"),
            reason: format!(
                "{subject} a unitary circuit without measurements, resets, or conditionals"
            ),
        });
    }
    Ok(())
}

/// Provenance of a result read off the exact mixture: the density matrix,
/// placed wherever `kind` holds it.
fn exact_mixture_metadata(kind: &BackendKind) -> RunMetadata {
    let metadata = RunMetadata::exact(ResolvedBackend::DensityMatrix);
    #[cfg(feature = "gpu")]
    if matches!(kind, BackendKind::DensityMatrixGpu { .. }) {
        let mut on_device = metadata;
        on_device.placement = Placement::Device;
        return on_device;
    }
    #[cfg(not(feature = "gpu"))]
    let _ = kind;
    metadata
}

/// Exact output distribution of `circuit` under `noise_model`, evolved once on
/// the density matrix. The mixture carries every measurement branch at once,
/// so it answers for a whole shot only when the measurements are terminal.
fn exact_noisy_probabilities(
    kind: &BackendKind,
    circuit: &Circuit,
    noise_model: &noise::NoiseModel,
    initial_state: Option<StartState<'_>>,
    seed: u64,
) -> Result<Probabilities> {
    Ok(Probabilities::Dense(noise::density_matrix_probabilities(
        kind,
        circuit,
        noise_model,
        initial_state,
        seed,
    )?))
}

/// Draw classical bits from the exact noisy distribution, then apply readout
/// error. Readout acts on the outcome rather than on the state, so it is not
/// carried by the distribution and has to be applied per draw, on a stream of
/// its own so it does not track the state draws.
fn sample_exact_noisy_shots(
    probs: &Probabilities,
    circuit: &Circuit,
    noise_model: &noise::NoiseModel,
    num_shots: usize,
    seed: u64,
) -> Vec<Vec<bool>> {
    let bits = circuit.num_classical_bits;
    let mut shots = sample_shots(probs, &circuit.measurement_map(), bits, num_shots, seed);
    if noise_model.readout.iter().any(Option::is_some) {
        let readout = trajectory::written_readout(circuit, &noise_model.readout);
        let mut rng = trajectory::noise_rng(seed);
        for shot in &mut shots {
            trajectory::apply_readout_errors(shot, &readout, &mut rng);
        }
    }
    shots
}

#[inline]
fn probs_only_result(probs: Vec<f64>, metadata: RunMetadata) -> RunOutcome {
    RunOutcome {
        probabilities: Some(Probabilities::Dense(probs)),
        classical_bits: vec![],
        metadata,
        saves: Vec::new(),
    }
}

/// Fail rather than return a run that quietly skipped a save point.
///
/// Routes that cannot read a state mid-circuit decline by name where they meet
/// the instruction. This is the backstop for a route that instead drops the
/// instruction while rewriting the circuit, which is what the decomposed and
/// exact-mixture paths do.
fn ensure_saves_recorded(circuit: &Circuit, outcome: &RunOutcome) -> Result<()> {
    let wanted = circuit.save_count();
    if wanted == outcome.saves.len() {
        return Ok(());
    }
    Err(PrismError::IncompatibleBackend {
        backend: format!("{:?}", outcome.metadata.backend),
        reason: format!(
            "circuit carries {wanted} save point(s) and the route that ran recorded {}",
            outcome.saves.len()
        ),
    })
}

fn try_backend_probabilities(backend: &dyn Backend) -> Result<Option<Probabilities>> {
    if let Some(factored) = backend.block_probabilities() {
        return Ok(Some(factored));
    }
    match backend.probabilities() {
        Ok(probs) => Ok(Some(Probabilities::Dense(probs))),
        Err(PrismError::BackendUnsupported { .. }) => Ok(None),
        Err(err) => Err(err),
    }
}

/// Core execution: fuse, init, apply, extract.
fn execute(backend: &mut dyn Backend, circuit: &Circuit, opts: &SimOptions) -> Result<RunOutcome> {
    let expanded = expand_for_backend(&*backend, circuit);
    let fused = fuse_for_backend(&*backend, &expanded);
    execute_circuit(backend, &fused, opts)
}

/// The execution route a template settles on, for a caller that holds one
/// circuit across many bindings and fuses it itself.
///
/// Resolving the route from the template rather than from a bound circuit
/// matters twice over: it is the work being amortized, and a fused stream
/// misreports the circuit to dispatch (a fused Clifford circuit no longer
/// looks Clifford).
pub(crate) struct PreparedRoute {
    plan: BackendPlan,
    supports_fused: bool,
    fusion_state_qubits: usize,
    /// Backend held across points. `init` reuses its state buffer when the
    /// width matches, so a sweep pays one `2^n` allocation rather than one per
    /// point. Rebuilt when the seed changes, since the seed feeds its RNG.
    held: Option<(u64, Box<dyn Backend + Send>)>,
    /// Set when the template reads the RNG, so every point gets a fresh backend
    /// and measures what a solo run with its seed measures.
    draws_randomness: bool,
    values: PreparedValues,
    observable: PreparedObservable,
}

/// Who answers [`Simulate::expectation_values`] for a prepared template.
#[derive(Clone, Copy)]
enum PreparedValues {
    Held,
    /// `Auto` from [`TENSOR_ROUTE_MIN_QUBITS`]: the bound circuit is offered to
    /// the scalar tensor route first, as `Simulate` offers it.
    TensorFirst,
    /// `Simulate` would take a route the held backend does not reproduce.
    Simulate,
}

/// Who answers [`Simulate::observable_expectation`] for a prepared template.
#[derive(Clone, Copy)]
enum PreparedObservable {
    /// Term by term on the held backend, with no variance.
    Held,
    /// The grouped evaluator on the held statevector, with its variance.
    Grouped,
    /// `Simulate` would take a route the held backend does not reproduce.
    Simulate,
}

impl PreparedRoute {
    /// True when the chosen backend accepts fused gates, so the caller should
    /// hand `run` a fused stream rather than the bound template.
    pub(crate) fn supports_fused(&self) -> bool {
        self.supports_fused
    }

    /// The width `simulate` fuses at for the held backend, so a stream fused
    /// here matches the one a direct run applies.
    pub(crate) fn fusion_state_qubits(&self) -> usize {
        self.fusion_state_qubits
    }

    /// Apply `circuit` verbatim, with no further fusion.
    pub(crate) fn run(&mut self, circuit: &Circuit, seed: u64) -> Result<RunOutcome> {
        if self.draws_randomness {
            let mut backend = self.plan.build(seed);
            return execute_circuit(&mut *backend, circuit, &SimOptions::default());
        }
        execute_circuit(self.held_backend(seed), circuit, &SimOptions::default())
    }

    /// [`Simulate::expectation_values`] on `bound`, evaluated by applying
    /// `applied` verbatim, or `None` when the caller should ask `Simulate`.
    pub(crate) fn expectation_values(
        &mut self,
        bound: &Circuit,
        applied: &Circuit,
        observables: &[Vec<PauliTerm>],
        seed: u64,
    ) -> Option<Result<Vec<f64>>> {
        match self.values {
            PreparedValues::Simulate => return None,
            PreparedValues::TensorFirst => {
                if let Some(result) = tensor_route_expectations(bound, observables) {
                    return Some(result.map(ExpectationResult::into_values));
                }
            }
            PreparedValues::Held => {}
        }
        Some(
            self.held_expectation_values(applied, observables, seed)
                .map(|(values, _)| values),
        )
    }

    /// [`Simulate::observable_expectation`] on the circuit `applied` was bound
    /// from, or `None` when the caller should ask `Simulate`.
    pub(crate) fn observable_expectation(
        &mut self,
        applied: &Circuit,
        observable: &PauliObservable,
        seed: u64,
    ) -> Option<Result<ObservableExpectation>> {
        match self.observable {
            PreparedObservable::Simulate => None,
            PreparedObservable::Held => Some(
                self.held_expectation_values(applied, &observable_vecs(observable), seed)
                    .map(|(values, metadata)| {
                        weighted_observable_result(observable, &values, None, metadata)
                    }),
            ),
            PreparedObservable::Grouped => Some(self.held_grouped(applied, observable, seed)),
        }
    }

    fn held_backend(&mut self, seed: u64) -> &mut dyn Backend {
        if !matches!(&self.held, Some((s, _)) if *s == seed) {
            self.held = Some((seed, self.plan.build(seed)));
        }
        let (_, backend) = self.held.as_mut().expect("just built");
        &mut **backend
    }

    fn held_expectation_values(
        &mut self,
        circuit: &Circuit,
        observables: &[Vec<PauliTerm>],
        seed: u64,
    ) -> Result<(Vec<f64>, RunMetadata)> {
        for observable in observables {
            validate_observable(observable, circuit.num_qubits)?;
        }
        let backend = self.held_backend(seed);
        execute_circuit(&mut *backend, circuit, &SimOptions::classical_only())?;
        Ok((
            backend.pauli_expectations(observables)?,
            backend_metadata(backend),
        ))
    }

    fn held_grouped(
        &mut self,
        circuit: &Circuit,
        observable: &PauliObservable,
        seed: u64,
    ) -> Result<ObservableExpectation> {
        for (_, factors) in observable.terms() {
            validate_observable(factors, circuit.num_qubits)?;
        }
        let backend = self.held_backend(seed);
        execute_circuit(&mut *backend, circuit, &SimOptions::classical_only())?;
        let statevector = backend
            .as_any()
            .and_then(|any| any.downcast_ref::<StatevectorBackend>())
            .expect("a grouped answer is settled only on a statevector plan");
        grouped_expectation_on_state(statevector, observable)
    }
}

/// Settle the route for `template`, or `None` when it takes one that reshapes
/// execution (decomposition, stabilizer rank, temporal Clifford) or expands
/// `QftBlock`, where a caller-supplied stream has nowhere to go.
pub(crate) fn prepared_route(kind: &BackendKind, template: &Circuit) -> Option<PreparedRoute> {
    if !kind.is_auto() && validate_explicit_backend(kind, template).is_err() {
        return None;
    }
    let ProbabilityRoute::Direct {
        has_partial_independence,
    } = plan_probability_route(kind, template)
    else {
        return None;
    };
    let ExecutionPlan::Backend(plan) = resolve(kind, template, has_partial_independence) else {
        return None;
    };
    let probe = plan.build(0);
    let has_qft_block = crate::circuit::any_gate(&template.instructions, &mut |gate| {
        matches!(gate, crate::gates::Gate::QftBlock { .. })
    });
    if has_qft_block && !probe.supports_qft_block() {
        return None;
    }
    let has_pauli_rot = crate::circuit::any_gate(&template.instructions, &mut |gate| {
        matches!(gate, crate::gates::Gate::PauliRot(_))
    });
    if has_pauli_rot && !probe.supports_pauli_rotation() {
        return None;
    }
    let (values, observable) = prepared_expectations(kind, template, &*probe);
    Some(PreparedRoute {
        supports_fused: probe.supports_fused_gates(),
        fusion_state_qubits: probe.fusion_state_qubits(template.num_qubits),
        plan,
        held: None,
        draws_randomness: draws_randomness(template),
        values,
        observable,
    })
}

/// Settle which expectation terminals the held backend `probe` answers the
/// way `Simulate` answers them.
///
/// Both decline a template the terminals reject, whose error `Simulate`
/// reports, and a Clifford template under `Auto` or a stabilizer kind, which
/// `Simulate` propagates through the Pauli route rather than a backend. The
/// weighted observable takes the grouped evaluator wherever `Simulate` would,
/// so its variance survives, and declines when the held backend is not the
/// statevector that evaluator reads.
fn prepared_expectations(
    kind: &BackendKind,
    template: &Circuit,
    probe: &dyn Backend,
) -> (PreparedValues, PreparedObservable) {
    let clifford_route =
        (kind.is_auto() || kind.is_stabilizer_family()) && template.is_clifford_only();
    if draws_randomness(template)
        || template.save_count() > 0
        || clifford_route
        || !probe.supports_pauli_expectation()
    {
        return (PreparedValues::Simulate, PreparedObservable::Simulate);
    }
    let fits_dense = template.num_qubits <= max_statevector_qubits();
    let values = if kind.is_auto() && fits_dense && template.num_qubits >= TENSOR_ROUTE_MIN_QUBITS {
        PreparedValues::TensorFirst
    } else {
        PreparedValues::Held
    };
    let grouped = match kind {
        BackendKind::Statevector => true,
        #[cfg(feature = "gpu")]
        BackendKind::StatevectorGpu { .. } => true,
        _ => kind.is_auto() && fits_dense,
    };
    let holds_statevector = probe
        .as_any()
        .is_some_and(|any| any.is::<StatevectorBackend>());
    let observable = match (grouped, holds_statevector) {
        (false, _) => PreparedObservable::Held,
        (true, true) => PreparedObservable::Grouped,
        (true, false) => PreparedObservable::Simulate,
    };
    (values, observable)
}

/// Shared init → apply → extract logic.
fn execute_circuit(
    backend: &mut dyn Backend,
    circuit: &Circuit,
    opts: &SimOptions,
) -> Result<RunOutcome> {
    backend.init(circuit.num_qubits, circuit.num_classical_bits)?;
    let saves = apply_recording_saves(backend, &circuit.instructions)?;

    let probabilities = if opts.probabilities {
        try_backend_probabilities(backend)?
    } else {
        None
    };

    Ok(RunOutcome {
        classical_bits: backend.classical_results().to_vec(),
        probabilities,
        metadata: backend_metadata(backend),
        saves,
    })
}

/// Apply `instructions`, reading the state at each save point.
///
/// The stream is handed to the backend in the segments between save points, so
/// no backend implements a save and none can batch across one. A circuit with
/// no save points takes one call.
fn apply_recording_saves(
    backend: &mut dyn Backend,
    instructions: &[Instruction],
) -> Result<Vec<SaveRecord>> {
    let first = instructions
        .iter()
        .position(|i| matches!(i, Instruction::Save { .. }));
    let Some(first) = first else {
        backend.apply_instructions(instructions)?;
        return Ok(Vec::new());
    };

    let mut saves = Vec::new();
    let mut start = 0;
    for (i, inst) in instructions.iter().enumerate().skip(first) {
        let Instruction::Save { spec, label, .. } = inst else {
            continue;
        };
        backend.apply_instructions(&instructions[start..i])?;
        saves.push(SaveRecord {
            label: label.clone(),
            value: read_save(backend, *spec, label)?,
        });
        start = i + 1;
    }
    backend.apply_instructions(&instructions[start..])?;
    Ok(saves)
}

/// Read one save off the backend that is mid-circuit.
///
/// Every spec maps to an export the backend already offers, so a backend that
/// cannot produce the requested form fails here with its own message rather
/// than needing a save implementation of its own.
fn read_save(backend: &mut dyn Backend, spec: SaveSpec, label: &str) -> Result<SavedValue> {
    let named = |err: crate::error::PrismError| match err {
        crate::error::PrismError::IncompatibleBackend { backend, reason } => {
            crate::error::PrismError::IncompatibleBackend {
                backend,
                reason: format!("save point `{label}`: {reason}"),
            }
        }
        other => other,
    };
    match spec {
        SaveSpec::StateVector => backend
            .export_statevector()
            .map(SavedValue::StateVector)
            .map_err(named),
        SaveSpec::Probabilities => backend
            .probabilities()
            .map(SavedValue::Probabilities)
            .map_err(named),
        SaveSpec::DensityMatrix => {
            let whole: Vec<usize> = (0..backend.num_qubits()).collect();
            backend
                .reduced_density_matrix(&whole)
                .map(SavedValue::DensityMatrix)
                .map_err(named)
        }
    }
}

/// Provenance read off the engine that ran, after it ran. Exactness and
/// placement are reports rather than predictions: the MPS bound reflects what
/// this run discarded, and the placement reflects where the amplitudes ended up
/// after any device fallback.
pub(crate) fn backend_metadata(backend: &dyn Backend) -> RunMetadata {
    let mut metadata =
        RunMetadata::new(backend.resolved(), backend.exactness(), backend.placement());
    metadata.bond = backend.bond_report();
    metadata
}

#[cfg(test)]
fn run(circuit: &Circuit, seed: u64) -> Result<RunOutcome> {
    run_with(BackendKind::Auto, circuit, seed)
}

/// Build the backend `kind` selects and run the circuit on it. For a
/// pre-constructed backend, use [`run_on`].
pub(crate) fn run_with(kind: BackendKind, circuit: &Circuit, seed: u64) -> Result<RunOutcome> {
    run_with_internal(kind, circuit, seed, SimOptions::default())
}

fn run_with_internal(
    kind: BackendKind,
    circuit: &Circuit,
    seed: u64,
    opts: SimOptions,
) -> Result<RunOutcome> {
    if !kind.is_auto() {
        validate_explicit_backend(&kind, circuit)?;
    }
    // The distributed backend runs the whole circuit across ranks in lockstep.
    // Subsystem decomposition, Clifford+T, and temporal-Clifford shortcuts all
    // reshape execution per sub-block, which would desynchronize the collective
    // calls every rank must issue in the same order. Dispatch directly.
    #[cfg(feature = "distributed")]
    if matches!(kind, BackendKind::StatevectorDistributed { .. }) {
        let mut backend = resolve_backend(&kind, circuit, false).build(seed);
        if opts.probabilities {
            // The gather cap follows from the register width alone, so a
            // register past it is rejected here rather than after the run has
            // been paid for. Same check `probabilities()` applies, so the
            // message does not depend on where it surfaced.
            crate::backend::dense_probability_len(backend.name(), circuit.num_qubits)?;
        }
        return execute(&mut *backend, circuit, &opts);
    }
    let route = plan_probability_route(&kind, circuit);
    let outcome = run_route(&kind, circuit, seed, opts, &route)?;
    ensure_saves_recorded(circuit, &outcome)?;
    Ok(outcome)
}

/// Execute one seed of a planned probability route. Shot loops plan once and
/// call this per seed.
fn run_route(
    kind: &BackendKind,
    circuit: &Circuit,
    seed: u64,
    opts: SimOptions,
    route: &ProbabilityRoute,
) -> Result<RunOutcome> {
    match route {
        ProbabilityRoute::FactoredStabilizer => {
            let mut backend =
                crate::backend::factored_stabilizer::FactoredStabilizerBackend::new(seed);
            let fs_opts = if circuit.num_qubits > 64 {
                SimOptions {
                    probabilities: false,
                }
            } else {
                opts
            };
            execute(&mut backend, circuit, &fs_opts)
        }
        ProbabilityRoute::Decomposed(components) => {
            run_decomposed(kind, components, circuit, seed, &opts)
        }
        ProbabilityRoute::TemporalClifford {
            has_partial_independence,
        } => match plan_temporal_clifford(kind, circuit) {
            Some(tc) => run_temporal_clifford(&tc, seed, opts.probabilities),
            None => run_direct(kind, circuit, seed, opts, *has_partial_independence),
        },
        ProbabilityRoute::Direct {
            has_partial_independence,
        } => run_direct(kind, circuit, seed, opts, *has_partial_independence),
    }
}

fn run_direct(
    kind: &BackendKind,
    circuit: &Circuit,
    seed: u64,
    opts: SimOptions,
    has_partial_independence: bool,
) -> Result<RunOutcome> {
    match resolve(kind, circuit, has_partial_independence) {
            ExecutionPlan::Backend(plan) => {
                let mut backend = plan.build(seed);
                execute(&mut *backend, circuit, &opts)
            }
            ExecutionPlan::StabilizerRank => {
                let sr = stabilizer_rank::run_stabilizer_rank(circuit, seed)?;
                Ok(probs_only_result(
                    sr.probabilities,
                    RunMetadata::exact(ResolvedBackend::StabilizerRank),
                ))
            }
            ExecutionPlan::StochasticPauli { num_samples } => {
                Err(crate::error::PrismError::IncompatibleBackend {
                    backend: format!(
                        "{:?}",
                        BackendKind::StochasticPauli { num_samples }
                    ),
                    reason: "StochasticPauli produces marginal estimates only; use `simulate(...).marginals()`".into(),
                })
            }
            ExecutionPlan::DeterministicPauli { truncation } => {
                Err(crate::error::PrismError::IncompatibleBackend {
                    backend: format!(
                        "{:?}",
                        BackendKind::DeterministicPauli { truncation }
                    ),
                    reason: "DeterministicPauli produces marginals only; use `simulate(...).marginals()`".into(),
                })
            }
            ExecutionPlan::PauliPath => Err(reject_pauli_path("a single run")),
    }
}

/// Execute a circuit on a pre-constructed backend. For automatic dispatch,
/// use [`simulate`].
pub fn run_on(backend: &mut dyn Backend, circuit: &Circuit) -> Result<RunOutcome> {
    execute(backend, circuit, &SimOptions::default())
}

/// Execute a circuit on a pre-constructed backend from a start state other than
/// |0...0⟩, the [`run_on`] sibling for a caller holding its own amplitudes.
///
/// The backend must accept one; see [`Backend::init_from_amplitudes`] for the
/// validation applied to `initial_state` and for which backends decline it.
pub fn run_on_state(
    backend: &mut dyn Backend,
    circuit: &Circuit,
    initial_state: &[Complex64],
) -> Result<RunOutcome> {
    check_initial_state_len(StartState::Amplitudes(initial_state), circuit.num_qubits)?;
    backend.init_from_amplitudes(initial_state.to_vec(), circuit.num_classical_bits)?;
    let saves = apply_fused_circuit(backend, circuit)?;
    Ok(RunOutcome {
        classical_bits: backend.classical_results().to_vec(),
        probabilities: try_backend_probabilities(backend)?,
        metadata: backend_metadata(backend),
        saves,
    })
}

/// Run several circuits, holding one backend across those that can share it.
///
/// Under the `parallel` feature, a batch whose circuits all sit at 16 qubits or
/// fewer splits across Rayon workers, each holding its own backend: a 200-point
/// sweep of a two-layer hardware-efficient ansatz ran 3.2x to 4.8x faster than a
/// loop over [`simulate`] at 4 to 12 qubits on a four-core i7-6700K, and one circuit
/// per core beat the kernels' own threads at 14 and 16. A wider circuit, a density
/// matrix from 9 qubits, or a GPU or distributed kind keeps the whole batch on one
/// thread.
///
/// On one thread the saving is the backend construction and its `2^n`
/// allocation, which `init` reuses when the next circuit has the same width;
/// route analysis is planned per circuit either way. On the same sweep that read
/// within about 10% of the loop, in both directions.
///
/// Results are identical to running each circuit on its own with the same seed.
/// A circuit that draws randomness (a measurement, a reset, or a classical
/// condition) gets a backend of its own, because a held backend would carry its
/// RNG forward into the next circuit and change what the next one measures.
///
/// # Errors
/// Returns the error of the first failing circuit in list order. A caller that
/// wants every circuit's result should call this per circuit.
pub fn run_batch(circuits: &[Circuit], kind: BackendKind, seed: u64) -> Result<Vec<RunOutcome>> {
    #[cfg(feature = "parallel")]
    if circuits.len() > 1
        && circuits
            .iter()
            .all(|c| runs_split_across_workers(&kind, c.num_qubits))
    {
        use rayon::prelude::*;
        let outcomes: Vec<Result<RunOutcome>> = circuits
            .par_iter()
            .map_init(
                || None,
                |held, circuit| run_batch_entry(&kind, circuit, seed, held),
            )
            .collect();
        return outcomes.into_iter().collect();
    }

    let mut held = None;
    circuits
        .iter()
        .map(|circuit| run_batch_entry(&kind, circuit, seed, &mut held))
        .collect()
}

type HeldBackend = Option<(dispatch::BackendPlan, usize, Box<dyn Backend + Send>)>;

fn run_batch_entry(
    kind: &BackendKind,
    circuit: &Circuit,
    seed: u64,
    held: &mut HeldBackend,
) -> Result<RunOutcome> {
    let Some(plan) = batch_plan(kind, circuit) else {
        *held = None;
        return run_with_internal(kind.clone(), circuit, seed, SimOptions::default());
    };
    if !kind.is_auto() {
        validate_explicit_backend(kind, circuit)?;
    }

    let reusable = matches!(
        &*held,
        Some((p, width, _)) if *width == circuit.num_qubits && dispatch::same_plan(p, &plan)
    );
    if !reusable {
        *held = Some((plan.clone(), circuit.num_qubits, plan.build(seed)));
    }
    let (_, _, backend) = held.as_mut().expect("just built");
    execute(&mut **backend, circuit, &SimOptions::default())
}

/// Width below which separate runs split across Rayon workers, above the kernels'
/// own parallel floor: on a 4-core host a 200-circuit two-layer ansatz batch split
/// one circuit per core took 18% of the unsplit time at 14 qubits and 53% at 16, and
/// 10% more at 18. Under miri it follows the kernels' floor down.
#[cfg(all(feature = "parallel", not(miri)))]
const RUN_SPLIT_QUBITS: usize = 17;
#[cfg(all(feature = "parallel", miri))]
const RUN_SPLIT_QUBITS: usize = crate::backend::PARALLEL_THRESHOLD_QUBITS;

/// Whether separate runs of `num_qubits`-wide circuits on `kind` should split across
/// Rayon workers: only below `RUN_SPLIT_QUBITS`, counting a density matrix at
/// twice its width since each worker holds a state, and never on a kind bound to one
/// device or rank context.
#[cfg(feature = "parallel")]
pub(crate) fn runs_split_across_workers(kind: &BackendKind, num_qubits: usize) -> bool {
    if bound_to_one_context(kind) {
        return false;
    }
    let width = if kind.is_density_matrix() {
        2 * num_qubits
    } else {
        num_qubits
    };
    width < RUN_SPLIT_QUBITS
}

/// Whether per-shot runs on `kind` should split across Rayon workers, given the
/// resolved family and width of each state one shot holds. See
/// [`state_splits_across_workers`] for the per-state rule.
#[cfg(feature = "parallel")]
fn shots_split_across_workers(kind: &BackendKind, states: &[(ResolvedBackend, usize)]) -> bool {
    !bound_to_one_context(kind)
        && states
            .iter()
            .all(|&(family, width)| state_splits_across_workers(family, width))
}

/// Whether separate runs each holding one `width`-qubit state of `family` should
/// split across Rayon workers: only below the width where that family's own
/// kernels go parallel, since each worker holds a state and a shot task nested
/// inside kernel joins can pile stolen frames onto one worker stack.
///
/// A tableau parallelizes its row loops only from `MIN_QUBITS_FOR_PAR_GATES`
/// qubits, and a product state parallelizes no gate, measurement or reset. A
/// density matrix counts at twice its width.
#[cfg(feature = "parallel")]
pub(crate) fn state_splits_across_workers(family: ResolvedBackend, width: usize) -> bool {
    use crate::backend::{MIN_QUBITS_FOR_PAR_GATES, PARALLEL_THRESHOLD_QUBITS};
    match family {
        ResolvedBackend::Stabilizer | ResolvedBackend::FactoredStabilizer => {
            width < MIN_QUBITS_FOR_PAR_GATES
        }
        ResolvedBackend::ProductState => true,
        ResolvedBackend::DensityMatrix => 2 * width < PARALLEL_THRESHOLD_QUBITS,
        _ => width < PARALLEL_THRESHOLD_QUBITS,
    }
}

#[cfg(feature = "parallel")]
fn bound_to_one_context(kind: &BackendKind) -> bool {
    match kind {
        #[cfg(feature = "gpu")]
        BackendKind::AutoGpu { .. }
        | BackendKind::StatevectorGpu { .. }
        | BackendKind::DensityMatrixGpu { .. }
        | BackendKind::StabilizerGpu { .. } => true,
        #[cfg(feature = "distributed")]
        BackendKind::StatevectorDistributed { .. } => true,
        _ => false,
    }
}

/// The plan a circuit can share, or `None` when it must run on its own.
///
/// A circuit that draws randomness needs a fresh RNG to match its solo run, and
/// anything that is not a direct backend route has no single plan to hold.
fn batch_plan(kind: &BackendKind, circuit: &Circuit) -> Option<dispatch::BackendPlan> {
    if draws_randomness(circuit) {
        return None;
    }
    let ProbabilityRoute::Direct {
        has_partial_independence,
    } = plan_probability_route(kind, circuit)
    else {
        return None;
    };
    match resolve(kind, circuit, has_partial_independence) {
        ExecutionPlan::Backend(plan) => Some(plan),
        _ => None,
    }
}

/// True when a run of `circuit` reads the backend RNG, so a backend carried over
/// from an earlier run would change its outcome.
fn draws_randomness(circuit: &Circuit) -> bool {
    circuit.instructions.iter().any(|i| {
        matches!(
            i,
            Instruction::Measure { .. }
                | Instruction::Reset { .. }
                | Instruction::Conditional { .. }
                | Instruction::Region(_)
        )
    })
}

/// Parse an OpenQASM string and execute with automatic backend selection.
pub fn run_qasm(qasm: &str, seed: u64) -> Result<RunOutcome> {
    let circuit = crate::circuit::openqasm::parse(qasm)?;
    simulate(&circuit).seed(seed).run()
}

#[cfg(test)]
fn run_shots(circuit: &Circuit, num_shots: usize, seed: u64) -> Result<ShotsResult> {
    run_shots_with(BackendKind::Auto, circuit, num_shots, seed)
}

pub(crate) fn supports_compiled_measurement_sampling(circuit: &Circuit) -> bool {
    circuit.is_clifford_only()
        && !circuit.has_resets()
        && circuit.has_terminal_measurements_only()
        && circuit
            .instructions
            .iter()
            .any(|inst| matches!(inst, Instruction::Measure { .. }))
}

fn supports_deferred_measurement_sampling(circuit: &Circuit) -> bool {
    circuit.is_clifford_only()
        && (circuit.has_resets() || !circuit.has_terminal_measurements_only())
        && circuit
            .instructions
            .iter()
            .any(|inst| matches!(inst, Instruction::Measure { .. }))
        && !circuit.instructions.iter().any(|inst| {
            matches!(
                inst,
                Instruction::Conditional { .. } | Instruction::Region(_)
            )
        })
}

fn is_clifford_sampler_kind(kind: &BackendKind) -> bool {
    if kind.is_auto() {
        return true;
    }
    match kind {
        BackendKind::Stabilizer | BackendKind::FactoredStabilizer => true,
        #[cfg(feature = "gpu")]
        BackendKind::StabilizerGpu { .. } => true,
        _ => false,
    }
}

fn should_use_compiled_clifford_sampling(
    kind: &BackendKind,
    circuit: &Circuit,
    num_shots: usize,
) -> bool {
    num_shots >= 1
        && supports_compiled_measurement_sampling(circuit)
        && is_clifford_sampler_kind(kind)
}

fn should_use_deferred_clifford_sampling(
    kind: &BackendKind,
    circuit: &Circuit,
    num_shots: usize,
) -> bool {
    num_shots >= 2
        && supports_deferred_measurement_sampling(circuit)
        && is_clifford_sampler_kind(kind)
}

fn compile_measurements_for_kind(
    kind: &BackendKind,
    circuit: &Circuit,
    seed: u64,
) -> Result<compiled::CompiledSampler> {
    #[cfg(not(feature = "gpu"))]
    let _ = kind;

    let sampler = compiled::compile_measurements(circuit, seed)?;

    #[cfg(feature = "gpu")]
    if let BackendKind::StabilizerGpu { context } = kind {
        return Ok(sampler.with_gpu(context.clone()));
    }

    Ok(sampler)
}

/// Components to decompose with, and whether every gate is Clifford.
type Decomposition = (Vec<Vec<usize>>, bool);

/// Independence analysis shared by the routing prelude in
/// `run_with_internal`, the shots slow path, and the terminal fast-path
/// candidacy. Returns the decomposition when full decomposition should fire,
/// plus the partial-independence flag otherwise.
fn analyze_independence(circuit: &Circuit) -> (Option<Decomposition>, bool) {
    if circuit.num_qubits >= MIN_DECOMPOSITION_QUBITS {
        let (components, clifford_only) = circuit.subsystems_and_clifford();
        if components.len() > 1 {
            if should_decompose(&components, circuit.num_qubits) {
                return (Some((components, clifford_only)), false);
            }
            return (None, true);
        }
    }
    (None, false)
}

/// `(t_count, stabilizer_rank_budget)` when the auto Clifford+T family gate
/// passes. Callers apply their own per-entry-point T-count ceilings.
fn auto_clifford_t_budget(circuit: &Circuit) -> Option<(usize, usize)> {
    (circuit.is_clifford_plus_t() && circuit.has_t_gates()).then(|| {
        (
            circuit.t_count(),
            stabilizer_rank_budget(circuit.num_qubits),
        )
    })
}

/// Auto-dispatch gate for the Clifford+T stabilizer-rank shortcut: the T
/// count must fit both the caller's ceiling and the size-derived
/// stabilizer-rank budget. Returns the T count when the shortcut applies.
pub(super) fn auto_stabilizer_rank_t_count(circuit: &Circuit, max_t: usize) -> Option<usize> {
    let (t, sr_budget) = auto_clifford_t_budget(circuit)?;
    (t <= max_t && t <= sr_budget).then_some(t)
}

/// Routing precedence for the probability path: decomposition (with the
/// large sparse-Clifford factored-stabilizer override), then temporal Clifford,
/// then direct family resolution. `run_with_internal` executes this plan and
/// `auto_terminal_statevector_candidate` consults it, so the two cannot
/// drift apart.
enum ProbabilityRoute {
    FactoredStabilizer,
    Decomposed(Vec<Vec<usize>>),
    /// The temporal-Clifford predicate holds; `run_route` builds the plan and
    /// falls back to direct resolution should the split come back empty.
    TemporalClifford {
        has_partial_independence: bool,
    },
    Direct {
        has_partial_independence: bool,
    },
}

fn plan_probability_route(kind: &BackendKind, circuit: &Circuit) -> ProbabilityRoute {
    let (decompose, has_partial_independence) = analyze_independence(circuit);
    if let Some((components, clifford_only)) = decompose {
        let max_block = components.iter().map(|c| c.len()).max().unwrap_or(0);
        if kind.is_auto()
            && clifford_only
            && circuit.num_qubits >= MIN_FACTORED_STABILIZER_QUBITS
            && max_block >= MIN_BLOCK_FOR_FACTORED_STAB
        {
            return ProbabilityRoute::FactoredStabilizer;
        }
        return ProbabilityRoute::Decomposed(components);
    }
    if has_temporal_clifford_opportunity(kind, circuit) {
        return ProbabilityRoute::TemporalClifford {
            has_partial_independence,
        };
    }
    ProbabilityRoute::Direct {
        has_partial_independence,
    }
}

/// True when the auto probability route falls through to direct family
/// resolution and that resolver picks the CPU statevector.
fn auto_terminal_statevector_candidate(circuit: &Circuit) -> bool {
    match plan_probability_route(&BackendKind::Auto, circuit) {
        ProbabilityRoute::Direct {
            has_partial_independence,
        } => auto_selects_cpu_statevector(circuit, has_partial_independence),
        _ => false,
    }
}

fn terminal_statevector_candidate(kind: &BackendKind, circuit: &Circuit) -> bool {
    if kind.is_auto() {
        return auto_terminal_statevector_candidate(circuit);
    }
    match kind {
        BackendKind::Statevector => true,
        #[cfg(feature = "gpu")]
        BackendKind::StatevectorGpu { .. } => true,
        _ => false,
    }
}

fn try_terminal_statevector_backend(
    kind: &BackendKind,
    circuit: &Circuit,
    seed: u64,
) -> Result<Option<TerminalStatevector>> {
    if !circuit.has_terminal_measurements_only() {
        return Ok(None);
    }

    let meas_map = circuit.measurement_map();
    if meas_map.is_empty() {
        return Ok(None);
    }

    let stripped = circuit.without_measurements();
    if !terminal_statevector_candidate(kind, &stripped) {
        return Ok(None);
    }

    let accel = accel_for(kind, Family::Statevector, stripped.num_qubits);
    let mut backend = build_statevector(&accel, seed);
    let expanded = expand_for_backend(&backend, &stripped);
    let fused = fuse_for_backend(&backend, &expanded);
    backend.init(fused.num_qubits, fused.num_classical_bits)?;
    backend.apply_instructions(&fused.instructions)?;

    Ok(Some((backend, meas_map)))
}

/// Build and run the backend for a terminal-measurement circuit when routing
/// lands on a single backend that samples from its own representation.
///
/// Returns `None` when the route is not a direct single backend, or when that
/// backend has no native sampler, leaving the dense probability path untouched.
/// The capability is probed before `init`, so a backend without one costs an
/// allocation and nothing else.
///
/// The product state is the one route taken past subsystem decomposition: it
/// already stores one factor per qubit, so splitting the circuit into
/// independent blocks pays a backend, a partition, and a merge per block to
/// rebuild what one native draw reads straight off the state, and above 64
/// qubits the merged block distribution does not exist at all. Every other
/// backend keeps the block split.
fn try_native_terminal_backend(
    kind: &BackendKind,
    stripped: &Circuit,
    seed: u64,
) -> Result<Option<Box<dyn Backend>>> {
    if !kind.is_auto() {
        validate_explicit_backend(kind, stripped)?;
    }
    let (decomposed, has_partial_independence) = match plan_probability_route(kind, stripped) {
        ProbabilityRoute::Direct {
            has_partial_independence,
        } => (false, has_partial_independence),
        ProbabilityRoute::Decomposed(_) => (true, false),
        _ => return Ok(None),
    };
    let ExecutionPlan::Backend(plan) = resolve(kind, stripped, has_partial_independence) else {
        return Ok(None);
    };
    if decomposed && !matches!(plan, BackendPlan::ProductState) {
        return Ok(None);
    }
    let mut backend = plan.build(seed);
    if !backend.supports_native_sampling() {
        return Ok(None);
    }
    execute(&mut *backend, stripped, &SimOptions::classical_only())?;
    Ok(Some(backend))
}

/// Build and run the backend for a marginal query when routing lands on a
/// single backend that evaluates observables on its own representation.
///
/// Returns `None` when the route is not a direct single backend, or when that
/// backend has no native observable path, leaving the dense probability route
/// untouched. The capability is probed before `init`, so a backend without one
/// costs an allocation and nothing else.
///
/// The product state is carried past subsystem decomposition for the reason
/// [`try_native_terminal_backend`] carries it: it already holds one factor per
/// qubit, and the decomposed route would build the `2^n` merged distribution to
/// read marginals a per-qubit expectation answers directly. Every other backend
/// keeps the block split, which has no single backend holding the joint state.
fn try_native_marginal_backend(
    kind: &BackendKind,
    circuit: &Circuit,
    seed: u64,
) -> Result<Option<Box<dyn Backend>>> {
    if !kind.is_auto() {
        validate_explicit_backend(kind, circuit)?;
    }
    let (decomposed, has_partial_independence) = match plan_probability_route(kind, circuit) {
        ProbabilityRoute::Direct {
            has_partial_independence,
        } => (false, has_partial_independence),
        ProbabilityRoute::Decomposed(_) => (true, false),
        _ => return Ok(None),
    };
    let ExecutionPlan::Backend(plan) = resolve(kind, circuit, has_partial_independence) else {
        return Ok(None);
    };
    if decomposed && !matches!(plan, BackendPlan::ProductState) {
        return Ok(None);
    }
    let mut backend = plan.build(seed);
    if !backend.supports_pauli_expectation() {
        return Ok(None);
    }
    execute(&mut *backend, circuit, &SimOptions::classical_only())?;
    Ok(Some(backend))
}

/// A source that answers every measurement in a circuit from a single state,
/// in the precedence [`prepare_shot_source`] applies. Shots and counts both
/// consume it, so the two cannot disagree about which shortcut applies.
enum ShotSource {
    /// `deferred` marks a sampler built from the measure/reset-deferred
    /// rewrite, whose renumbered measurements only read through `meas_map`.
    Compiled {
        sampler: Box<compiled::CompiledSampler>,
        meas_map: Vec<(usize, usize)>,
        deferred: bool,
    },
    TerminalStatevector {
        backend: Box<StatevectorBackend>,
        meas_map: Vec<(usize, usize)>,
    },
    /// A backend that draws basis states from its own representation.
    Native {
        backend: Box<dyn Backend>,
        meas_map: Vec<(usize, usize)>,
    },
    /// Dense output distribution of the measurement-stripped circuit, with the
    /// provenance of the run that produced it.
    TerminalProbabilities {
        probs: Probabilities,
        meas_map: Vec<(usize, usize)>,
        metadata: RunMetadata,
    },
    StabilizerRank,
    PerShot,
}

impl ShotSource {
    /// Which engine will answer, decided by `prepare_shot_source` and read here
    /// so the shot and count entry points do not re-derive it. `None` for the
    /// per-shot route, which builds one backend per shot and stamps its own.
    fn metadata(&self) -> Option<RunMetadata> {
        match self {
            ShotSource::Compiled { .. } => Some(
                RunMetadata::exact(ResolvedBackend::CompiledStabilizer)
                    .with_engine(Engine::CompiledSampler),
            ),
            ShotSource::TerminalStatevector { backend, .. } => Some(backend_metadata(&**backend)),
            ShotSource::Native { backend, .. } => Some(backend_metadata(&**backend)),
            ShotSource::TerminalProbabilities { metadata, .. } => Some(metadata.clone()),
            ShotSource::StabilizerRank => Some(RunMetadata::exact(ResolvedBackend::StabilizerRank)),
            ShotSource::PerShot => None,
        }
    }
}

/// Narrowest register `Auto` samples through the Clifford register rather than the
/// statevector. At 12 qubits the statevector still wins, 0.8 against 1.3 ms for 1024
/// shots; at 16 the register takes 1.5 ms against 3.3.
const MIN_AUTO_CLIFFORD_REGISTER_QUBITS: usize = 16;

/// Whether terminal shots of `circuit` go to the Clifford register ahead of the
/// statevector: Clifford+T, at least [`MIN_AUTO_CLIFFORD_REGISTER_QUBITS`] wide, and
/// fewer T gates than qubits, within the register's width cap.
fn auto_prefers_clifford_register(circuit: &Circuit) -> bool {
    circuit.num_qubits >= MIN_AUTO_CLIFFORD_REGISTER_QUBITS
        && circuit.has_terminal_measurements_only()
        && !circuit.has_resets()
        && auto_clifford_t_budget(circuit).is_some_and(|(t, _)| {
            t < circuit.num_qubits && t <= clifford_register::MAX_REGISTER_QUBITS
        })
}

/// Whether terminal shots of a circuit wider than [`MAX_STABILIZER_RANK_QUBITS`]
/// go to the stabilizer-rank expansion once the statevector has declined them.
fn auto_samples_wide_by_stabilizer_rank(circuit: &Circuit) -> bool {
    circuit.has_terminal_measurements_only()
        && circuit.num_qubits > MAX_STABILIZER_RANK_QUBITS
        && auto_stabilizer_rank_t_count(circuit, MAX_AUTO_T_COUNT_SHOTS).is_some()
}

/// Whether `Auto` shots of `circuit`, without noise or a start state, run on the
/// exact stabilizer-rank expansion: the register and wide-circuit tests in
/// [`prepare_shot_source`], or the per-shot test in [`run_shots_per_shot`].
fn auto_samples_by_stabilizer_rank(circuit: &Circuit) -> bool {
    let folded = circuit.fold_static_guards();
    let circuit = folded.as_ref();
    if circuit.has_terminal_measurements_only() {
        auto_prefers_clifford_register(circuit) || auto_samples_wide_by_stabilizer_rank(circuit)
    } else {
        auto_stabilizer_rank_t_count(circuit, MAX_AUTO_T_COUNT_SHOTS).is_some()
    }
}

/// Select and prepare the sampling source for `circuit`.
///
/// Preparation is real work: compiling a sampler, building and running a
/// backend, or executing the stripped circuit once. Call once per entry point
/// and match on the result rather than re-deriving the choice.
fn prepare_shot_source(
    kind: &BackendKind,
    circuit: &Circuit,
    num_shots: usize,
    seed: u64,
) -> Result<ShotSource> {
    if should_use_compiled_clifford_sampling(kind, circuit, num_shots) {
        return Ok(ShotSource::Compiled {
            sampler: Box::new(compile_measurements_for_kind(kind, circuit, seed)?),
            meas_map: circuit.measurement_map(),
            deferred: false,
        });
    }

    if should_use_deferred_clifford_sampling(kind, circuit, num_shots)
        && let Ok(deferred) = compiled::defer_measure_reset_circuit(circuit)
    {
        return Ok(ShotSource::Compiled {
            sampler: Box::new(compile_measurements_for_kind(kind, &deferred, seed)?),
            meas_map: deferred.measurement_map(),
            deferred: true,
        });
    }

    if kind.is_auto() && auto_prefers_clifford_register(circuit) {
        return Ok(ShotSource::StabilizerRank);
    }

    if let Some((backend, meas_map)) = try_terminal_statevector_backend(kind, circuit, seed)? {
        return Ok(ShotSource::TerminalStatevector {
            backend: Box::new(backend),
            meas_map,
        });
    }

    if matches!(kind, BackendKind::StabilizerRank) && circuit.has_t_gates() {
        return Ok(ShotSource::StabilizerRank);
    }
    if kind.is_auto() && auto_samples_wide_by_stabilizer_rank(circuit) {
        return Ok(ShotSource::StabilizerRank);
    }

    if circuit.has_terminal_measurements_only()
        || (kind.is_density_matrix() && circuit.has_terminal_measurements_under_reset_channel())
    {
        let stripped = circuit.without_measurements();
        if let Some(backend) = try_native_terminal_backend(kind, &stripped, seed)? {
            return Ok(ShotSource::Native {
                backend,
                meas_map: circuit.measurement_map(),
            });
        }
        let result = run_with_internal(kind.clone(), &stripped, seed, SimOptions::default())?;
        if let Some(probs) = result.probabilities {
            return Ok(ShotSource::TerminalProbabilities {
                probs,
                meas_map: circuit.measurement_map(),
                metadata: result.metadata,
            });
        }
    }

    Ok(ShotSource::PerShot)
}

#[cfg(test)]
fn run_counts(circuit: &Circuit, num_shots: usize, seed: u64) -> Result<HashMap<Vec<u64>, u64>> {
    run_counts_with(BackendKind::Auto, circuit, num_shots, seed).map(|(counts, _)| counts)
}

/// Execute a circuit multiple times with explicit backend selection and return counts.
///
/// For Clifford circuits with terminal measurements and no resets, Auto,
/// Stabilizer, FactoredStabilizer, and explicit `StabilizerGpu` route through
/// the compiled sampler's optimized counting path. Explicit `StabilizerGpu`
/// carries its GPU context into the compiled sampler so large shot runs avoid
/// the raw tableau measurement round-trips. Other circuits fall back to
/// per-shot simulation with counting.
///
/// Optimized terminal statevector paths sample counts directly from the output
/// distribution. The distribution is equivalent to materializing shots first,
/// but finite seeded counts may differ from `run_shots_with(...).counts()`.
pub(crate) fn run_counts_with(
    kind: BackendKind,
    circuit: &Circuit,
    num_shots: usize,
    seed: u64,
) -> Result<(HashMap<Vec<u64>, u64>, RunMetadata)> {
    #[cfg(feature = "distributed")]
    if matches!(kind, BackendKind::StatevectorDistributed { .. }) {
        let shots = run_shots_with(kind, circuit, num_shots, seed)?;
        return Ok((shots.counts(), shots.metadata));
    }

    let folded = circuit.fold_static_guards();
    let circuit = folded.as_ref();

    let bits = circuit.num_classical_bits;
    let source = prepare_shot_source(&kind, circuit, num_shots, seed)?;
    let Some(metadata) = source.metadata() else {
        let shots = run_shots_per_shot(kind, circuit, num_shots, seed, ShotRows::CountedOnly)?;
        return Ok((shots.counts(), shots.metadata));
    };
    let counts = match source {
        ShotSource::Compiled {
            mut sampler,
            meas_map,
            deferred,
        } => {
            if deferred {
                let packed = sampler.try_sample_bulk_packed(num_shots)?;
                counts_of(
                    packed_shots_to_classical_bits(&packed, &meas_map, bits),
                    bits,
                )
            } else {
                record_counts_to_classical_bits(
                    sampler.try_sample_counts(num_shots)?,
                    &meas_map,
                    bits,
                )
            }
        }
        ShotSource::TerminalStatevector { backend, meas_map } => {
            if backend.is_gpu_resident() {
                let probs = backend.probabilities()?;
                sample_counts_from_probs(&probs, &meas_map, bits, num_shots, seed)
            } else {
                sample_counts_from_state(
                    backend.state_vector(),
                    backend.probability_scale(),
                    &meas_map,
                    bits,
                    num_shots,
                    seed,
                )
            }
        }
        ShotSource::Native {
            mut backend,
            meas_map,
        } => {
            let samples = backend.sample_basis_states(num_shots, seed)?;
            counts_from_basis_samples(&samples, &meas_map, bits)
        }
        ShotSource::TerminalProbabilities {
            probs, meas_map, ..
        } => counts_of(sample_shots(&probs, &meas_map, bits, num_shots, seed), bits),
        ShotSource::StabilizerRank => {
            stabilizer_rank::run_stabilizer_rank_shots(circuit, num_shots, seed)?.counts()
        }
        ShotSource::PerShot => unreachable!("handled above"),
    };
    Ok((counts, metadata.with_shots(num_shots)))
}

fn counts_of(shots: Vec<Vec<bool>>, num_classical_bits: usize) -> HashMap<Vec<u64>, u64> {
    ShotsResult::from_shots(shots, num_classical_bits).counts()
}

#[cfg(test)]
fn run_marginals(circuit: &Circuit, seed: u64) -> Result<Vec<(f64, f64)>> {
    run_marginals_result_with(BackendKind::Auto, circuit, seed).map(MarginalsResult::into_vec)
}

#[cfg(test)]
fn run_marginals_with(kind: BackendKind, circuit: &Circuit, seed: u64) -> Result<Vec<(f64, f64)>> {
    run_marginals_result_with(kind, circuit, seed).map(MarginalsResult::into_vec)
}

/// Per-qubit marginals from native single-qubit Z expectations, for a backend
/// whose own representation answers them without a dense probability vector.
fn marginals_from_pauli_expectations(
    backend: &dyn Backend,
    num_qubits: usize,
) -> Result<MarginalsResult> {
    let observables: Vec<Vec<PauliTerm>> = (0..num_qubits).map(|q| vec![PauliTerm::z(q)]).collect();
    let expectations = backend.pauli_expectations(&observables)?;
    Ok(MarginalsResult {
        marginals: expectations_to_marginals(&expectations),
        metadata: backend_metadata(backend),
    })
}

/// Route-level exactness per the [`Exactness`] convention: a policy that can
/// discard marks the route approximate even on a run that discarded nothing.
/// Only a zero-threshold policy stays exact, since it overflows into an error
/// instead of an approximation; a budget always drops once the sum passes it.
/// The realized bound is a coefficient magnitude rather than a state overlap,
/// so it stays `total_discarded` on the engine result and never a fidelity
/// bound.
fn spd_metadata(truncation: &SpdTruncation) -> RunMetadata {
    if truncation.can_discard() {
        RunMetadata::approximate(ResolvedBackend::DeterministicPauli)
    } else {
        RunMetadata::exact(ResolvedBackend::DeterministicPauli)
    }
}

pub(crate) fn expectations_to_marginals(expectations: &[f64]) -> Vec<(f64, f64)> {
    expectations
        .iter()
        .map(|ez| {
            let p0 = ((1.0 + ez) / 2.0).clamp(0.0, 1.0);
            (p0, 1.0 - p0)
        })
        .collect()
}

pub(super) fn has_nonunitary_or_classical_ops(circuit: &Circuit) -> bool {
    circuit.instructions.iter().any(|inst| {
        matches!(
            inst,
            Instruction::Measure { .. }
                | Instruction::Reset { .. }
                | Instruction::Conditional { .. }
                | Instruction::Region(_)
        )
    })
}

fn supports_pauli_marginal_backend(circuit: &Circuit) -> bool {
    circuit.is_clifford_plus_t() && !has_nonunitary_or_classical_ops(circuit)
}

/// Gate-set support is left to the engines, which accept Clifford gates and
/// Pauli rotations and report the offending gate by name.
fn validate_pauli_marginal_backend(kind: &BackendKind, circuit: &Circuit) -> Result<()> {
    if has_nonunitary_or_classical_ops(circuit) {
        return Err(PrismError::IncompatibleBackend {
            backend: format!("{kind:?}"),
            reason: "Pauli marginal backends require a unitary circuit without measurements, resets, or conditionals".into(),
        });
    }
    Ok(())
}

fn run_marginals_result_with(
    kind: BackendKind,
    circuit: &Circuit,
    seed: u64,
) -> Result<MarginalsResult> {
    let n = circuit.num_qubits;

    match &kind {
        BackendKind::StochasticPauli { num_samples } => {
            validate_pauli_marginal_backend(&kind, circuit)?;
            let spp = unified_pauli::run_spp(circuit, *num_samples, seed)?;
            return Ok(MarginalsResult {
                marginals: expectations_to_marginals(&spp.expectations),
                metadata: RunMetadata::approximate(ResolvedBackend::StochasticPauli)
                    .with_shots(*num_samples),
            });
        }
        BackendKind::DeterministicPauli { truncation } => {
            validate_pauli_marginal_backend(&kind, circuit)?;
            let spd = unified_pauli::run_spd_with(circuit, truncation)?;
            return Ok(MarginalsResult {
                marginals: expectations_to_marginals(&spd.expectations),
                metadata: spd_metadata(truncation),
            });
        }
        _ => {}
    }

    if kind.is_auto()
        && supports_pauli_marginal_backend(circuit)
        && circuit.has_t_gates()
        && n >= MIN_QUBITS_FOR_SPD_AUTO
    {
        let exact = SpdTruncation::Threshold {
            epsilon: 0.0,
            max_terms: 0,
        };
        let spd = if n <= max_statevector_qubits() {
            unified_pauli::run_spd_exact_within(circuit, auto_spd_work_budget(circuit))?
        } else {
            Some(unified_pauli::run_spd_with(circuit, &exact)?)
        };
        if let Some(spd) = spd {
            return Ok(MarginalsResult {
                marginals: expectations_to_marginals(&spd.expectations),
                metadata: spd_metadata(&exact),
            });
        }
    }

    // The distributed backend answers a marginal from rank-local sums plus one
    // `Allreduce`, so it never needs the 2^n gather the fallback below takes.
    #[cfg(feature = "distributed")]
    if let BackendKind::StatevectorDistributed { context } = &kind {
        let mut backend =
            crate::backend::distributed_statevector::DistributedStatevectorBackend::new(
                context.clone(),
                seed,
            );
        execute(&mut backend, circuit, &SimOptions::classical_only())?;
        return marginals_from_pauli_expectations(&backend, n);
    }

    // A marginal is a `Z_q` expectation, so a unitary circuit gets the same
    // tensor routes the expectation terminal takes.
    if kind.is_auto() && !has_nonunitary_or_classical_ops(circuit) {
        let observables: Vec<Vec<PauliTerm>> = (0..n).map(|q| vec![PauliTerm::z(q)]).collect();
        let tensor_route = if n <= max_statevector_qubits() {
            tensor_route_expectations(circuit, &observables)
        } else {
            tensor_route_over_mps(circuit, &observables)
        };
        if let Some(result) = tensor_route {
            let result = result?;
            return Ok(MarginalsResult {
                marginals: expectations_to_marginals(&result.values),
                metadata: result.metadata,
            });
        }
    }

    if let Some(backend) = try_native_marginal_backend(&kind, circuit, seed)? {
        return marginals_from_pauli_expectations(&*backend, n);
    }

    let result = run_with(kind, circuit, seed)?;
    if let Some(probs) = &result.probabilities {
        Ok(MarginalsResult {
            marginals: probs.marginals(),
            metadata: result.metadata.clone(),
        })
    } else {
        Err(PrismError::BackendUnsupported {
            backend: "simulate".into(),
            operation: format!(
                "marginals for {} qubits without backend probability output",
                circuit.num_qubits
            ),
        })
    }
}

/// Compute `⟨ψ|P|ψ⟩` for each joint Pauli observable on a unitary circuit's
/// output state, using automatic backend selection. See
/// [`Simulate::expectation_values`] for explicit backend control.
///
/// # Examples
///
/// ```
/// use prism_q::{Circuit, Gate, PauliTerm, run_expectation_values};
///
/// let mut bell = Circuit::new(2, 0);
/// bell.add_gate(Gate::H, &[0]);
/// bell.add_gate(Gate::Cx, &[0, 1]);
///
/// let observables = vec![
///     vec![PauliTerm::z(0)],
///     vec![PauliTerm::z(0), PauliTerm::z(1)],
/// ];
/// let values = run_expectation_values(&bell, &observables, 42)?;
/// assert!(values[0].abs() < 1e-10); // <Z0> = 0
/// assert!((values[1] - 1.0).abs() < 1e-10); // <Z0 Z1> = 1
/// # Ok::<(), prism_q::PrismError>(())
/// ```
pub fn run_expectation_values(
    circuit: &Circuit,
    observables: &[Vec<PauliTerm>],
    seed: u64,
) -> Result<Vec<f64>> {
    run_expectation_values_with(BackendKind::Auto, circuit, observables, seed)
}

fn run_expectation_values_with(
    kind: BackendKind,
    circuit: &Circuit,
    observables: &[Vec<PauliTerm>],
    seed: u64,
) -> Result<Vec<f64>> {
    run_expectation_values_reported(kind, circuit, observables, seed)
        .map(ExpectationResult::into_values)
}

fn run_expectation_values_reported(
    kind: BackendKind,
    circuit: &Circuit,
    observables: &[Vec<PauliTerm>],
    seed: u64,
) -> Result<ExpectationResult> {
    require_unitary_circuit(&kind, circuit, "expectation values require")?;

    match &kind {
        BackendKind::StochasticPauli { num_samples } => {
            let mut values = Vec::with_capacity(observables.len());
            let mut std_errors = Vec::with_capacity(observables.len());
            for (i, obs) in observables.iter().enumerate() {
                let r = unified_pauli::run_spp_observable(
                    circuit,
                    obs,
                    *num_samples,
                    mix_seed(seed, i),
                )?;
                values.push(r.mean);
                std_errors.push(r.std_error);
            }
            Ok(ExpectationResult {
                values,
                std_errors: Some(std_errors),
                metadata: RunMetadata::approximate(ResolvedBackend::StochasticPauli)
                    .with_shots(*num_samples),
            })
        }
        BackendKind::DeterministicPauli { truncation } => {
            let mut values = Vec::with_capacity(observables.len());
            for obs in observables {
                let r = unified_pauli::run_spd_observable_with(circuit, obs, truncation)?;
                values.push(r.mean);
            }
            Ok(analytic_expectations(values, spd_metadata(truncation)))
        }
        _ if kind.is_auto() || kind.is_stabilizer_family() => {
            if circuit.is_clifford_only() {
                let mut values = Vec::with_capacity(observables.len());
                for obs in observables {
                    let r = unified_pauli::run_spd_observable(circuit, obs, 0.0, 0)?;
                    values.push(r.mean);
                }
                Ok(analytic_expectations(
                    values,
                    spd_metadata(&SpdTruncation::Threshold {
                        epsilon: 0.0,
                        max_terms: 0,
                    }),
                ))
            } else if kind.is_auto() {
                if circuit.num_qubits > max_statevector_qubits() {
                    if let Some(result) = tensor_route_over_mps(circuit, observables) {
                        return result;
                    }
                    return expectation_values_native(&kind, circuit, observables, seed);
                }
                if let Some(result) = tensor_route_expectations(circuit, observables) {
                    return result;
                }
                expectation_values_statevector(&kind, circuit, observables, seed)
            } else {
                Err(PrismError::IncompatibleBackend {
                    backend: format!("{kind:?}"),
                    reason: "stabilizer backends require a Clifford-only circuit".into(),
                })
            }
        }
        BackendKind::Statevector => {
            expectation_values_statevector(&kind, circuit, observables, seed)
        }
        #[cfg(feature = "gpu")]
        BackendKind::StatevectorGpu { .. } => {
            expectation_values_statevector(&kind, circuit, observables, seed)
        }
        // The scalar network prunes each observable to its light cone, which
        // the backend's stored network cannot: it no longer knows its gates.
        BackendKind::TensorNetwork => {
            validate_explicit_backend(&kind, circuit)?;
            for observable in observables {
                validate_observable(observable, circuit.num_qubits)?;
            }
            let values =
                crate::backend::tensornetwork::expectations_zero_state(circuit, observables)?;
            Ok(analytic_expectations(
                values,
                RunMetadata::exact(ResolvedBackend::TensorNetwork),
            ))
        }
        other => expectation_values_native(other, circuit, observables, seed),
    }
}

/// Compute `⟨H⟩` and its grouped-measurement variance for a weighted Pauli
/// observable, using automatic backend selection. See
/// [`Simulate::observable_expectation`] for explicit backend control and
/// [`ObservableExpectation::variance`] for the variance contract.
///
/// # Examples
///
/// ```
/// use prism_q::{Circuit, Gate, PauliObservable, PauliTerm, run_observable_expectation};
///
/// let mut circuit = Circuit::new(2, 0);
/// circuit.add_gate(Gate::H, &[0]);
/// circuit.add_gate(Gate::Cx, &[0, 1]);
/// circuit.add_gate(Gate::T, &[0]);
///
/// let hamiltonian = PauliObservable::from_terms([
///     (1.0, vec![PauliTerm::z(0)]),
///     (1.0, vec![PauliTerm::z(0), PauliTerm::z(1)]),
/// ])?;
///
/// // The T phase moves nothing here: both terms are Z-only, one commuting
/// // group covers them, and the variance is exactly Var(H) with outcomes
/// // 2 and 0 at probability 1/2 each.
/// let result = run_observable_expectation(&circuit, &hamiltonian, 42)?;
/// assert!((result.mean - 1.0).abs() < 1e-10);
/// assert!((result.variance.unwrap() - 1.0).abs() < 1e-10);
/// # Ok::<(), prism_q::PrismError>(())
/// ```
pub fn run_observable_expectation(
    circuit: &Circuit,
    observable: &PauliObservable,
    seed: u64,
) -> Result<ObservableExpectation> {
    run_observable_expectation_reported(BackendKind::Auto, circuit, observable, seed)
}

fn run_observable_expectation_reported(
    kind: BackendKind,
    circuit: &Circuit,
    observable: &PauliObservable,
    seed: u64,
) -> Result<ObservableExpectation> {
    require_unitary_circuit(&kind, circuit, "expectation values require")?;

    let grouped_statevector = match &kind {
        BackendKind::Statevector => true,
        #[cfg(feature = "gpu")]
        BackendKind::StatevectorGpu { .. } => true,
        _ => {
            kind.is_auto()
                && !circuit.is_clifford_only()
                && circuit.num_qubits <= max_statevector_qubits()
        }
    };
    if grouped_statevector {
        return grouped_expectation_statevector(&kind, circuit, observable, seed);
    }

    let result =
        run_expectation_values_reported(kind, circuit, &observable_vecs(observable), seed)?;
    Ok(weighted_observable_result(
        observable,
        &result.values,
        result.std_errors.as_deref(),
        result.metadata,
    ))
}

fn observable_vecs(observable: &PauliObservable) -> Vec<Vec<PauliTerm>> {
    observable
        .terms()
        .iter()
        .map(|(_, factors)| factors.clone())
        .collect()
}

/// Fold per-term values into the weighted mean; no grouped traversal ran, so
/// there is no variance to report. Independent per-term estimates combine
/// into the weighted-sum standard error.
fn weighted_observable_result(
    observable: &PauliObservable,
    values: &[f64],
    std_errors: Option<&[f64]>,
    metadata: RunMetadata,
) -> ObservableExpectation {
    let coefficients = observable.terms().iter().map(|(c, _)| *c);
    let mean = coefficients.clone().zip(values).map(|(c, v)| c * v).sum();
    let std_error = std_errors.map(|errors| {
        coefficients
            .zip(errors)
            .map(|(c, e)| (c * e).powi(2))
            .sum::<f64>()
            .sqrt()
    });
    ObservableExpectation {
        mean,
        variance: None,
        group_variances: None,
        std_error,
        metadata,
    }
}

/// Evaluate a weighted observable on the statevector: small groups from one
/// shared batched traversal, each large group's mean and variance from a
/// dedicated moments pass.
///
/// `Var(H_g) = <H_g^2> - <H_g>^2` per commuting group. For a small group the
/// square expands into pairwise product strings appended to the same
/// traversal that serves its term means; a group past the pair budget takes a
/// single-pass moment accumulation instead, on the state as run when the
/// group is Z-only and on a basis-rotated copy otherwise.
fn grouped_expectation_statevector(
    kind: &BackendKind,
    circuit: &Circuit,
    observable: &PauliObservable,
    seed: u64,
) -> Result<ObservableExpectation> {
    let terms = observable.terms();
    // Validate before the 2^n simulation so bad observables fail cheaply. The
    // mask reduction shifts by qubit index, so it follows the `init` that
    // rejects a width no mask can address.
    for (_, factors) in terms {
        validate_observable(factors, circuit.num_qubits)?;
    }

    let accel = accel_for(kind, Family::Statevector, circuit.num_qubits);
    let mut backend = build_statevector(&accel, seed);
    let expanded = expand_for_backend(&backend, circuit);
    let fused = fuse_for_backend(&backend, &expanded);
    backend.init(fused.num_qubits, fused.num_classical_bits)?;
    backend.apply_instructions(&fused.instructions)?;
    grouped_expectation_on_state(&backend, observable)
}

/// The reduction half of [`grouped_expectation_statevector`], on a backend
/// that has already run the circuit.
fn grouped_expectation_on_state(
    backend: &StatevectorBackend,
    observable: &PauliObservable,
) -> Result<ObservableExpectation> {
    let num_qubits = backend.num_qubits();
    let terms = observable.terms();
    let masks = terms
        .iter()
        .map(|(_, factors)| pauli_masks(factors, num_qubits))
        .collect::<Result<Vec<_>>>()?;
    let metadata = backend_metadata(backend);

    let grouping = observable.grouping();

    // Small groups get `<H_g^2>` from pairwise product strings appended to the
    // shared traversal: qubit-wise-commuting strings multiply phase-free
    // (shared qubits carry equal axes and cancel to identity), so each pair is
    // one more mask. Large groups fall back to a dedicated moments pass, which
    // costs a fixed number of state sweeps where the pair expansion grows
    // quadratically. Each group's mean comes out of whichever pass serves its
    // variance, so a large group's terms stay out of the shared traversal.
    let mut combined = Vec::with_capacity(masks.len());
    let mut pair_blocks: Vec<(usize, usize, usize, Vec<f64>)> = Vec::new();
    let mut deferred: Vec<usize> = Vec::new();
    for (gi, group) in grouping.groups.iter().enumerate() {
        let members = &group.term_indices;
        if members.len() * (members.len() - 1) / 2 > MAX_PAIR_MASKS_PER_GROUP {
            deferred.push(gi);
            continue;
        }
        let first_term = combined.len();
        combined.extend(members.iter().map(|&i| masks[i]));
        let first_pair = combined.len();
        let mut pair_coefficients = Vec::with_capacity(members.len() * (members.len() - 1) / 2);
        for (pos, &i) in members.iter().enumerate() {
            for &j in &members[pos + 1..] {
                let product_x = masks[i].0 ^ masks[j].0;
                let product_z = masks[i].1 ^ masks[j].1;
                combined.push((product_x, product_z, (product_x & product_z).count_ones()));
                pair_coefficients.push(2.0 * terms[i].0 * terms[j].0);
            }
        }
        pair_blocks.push((gi, first_term, first_pair, pair_coefficients));
    }

    // A device-resident state reduces every mask on the card; the host keeps
    // its state and norm for the moments pass below.
    let (values, host) = match backend.pauli_expectations_on_device(&combined) {
        Some(values) => (values?, None),
        None => {
            let state = backend.state_vector();
            let (values, norm) = pauli_expectations_from_masks(state, &combined);
            (values, Some((state, norm)))
        }
    };

    let mut mean: f64 = terms
        .iter()
        .filter(|(_, factors)| factors.is_empty())
        .map(|(c, _)| c)
        .sum();
    let mut group_variances = vec![0.0; grouping.groups.len()];

    for (gi, first_term, first_pair, pair_coefficients) in &pair_blocks {
        let group = &grouping.groups[*gi];
        let m1: f64 = group
            .term_indices
            .iter()
            .zip(&values[*first_term..])
            .map(|(&i, v)| terms[i].0 * v)
            .sum();
        let square_diag: f64 = group
            .term_indices
            .iter()
            .map(|&i| terms[i].0 * terms[i].0)
            .sum();
        let square_cross: f64 = pair_coefficients
            .iter()
            .zip(&values[*first_pair..])
            .map(|(c, v)| c * v)
            .sum();
        mean += m1;
        group_variances[*gi] = (square_diag + square_cross - m1 * m1).max(0.0);
    }

    // A Z-only group reads the state as run, any other group a copy rotated so
    // its terms read as Z strings (`x | z` per term). A device-resident state
    // takes both on the card and is exported only when the rotated copy does
    // not fit there.
    if !deferred.is_empty() {
        let requests: Vec<(Vec<u64>, Vec<f64>, u64, u64)> = deferred
            .iter()
            .map(|&gi| {
                let group = &grouping.groups[gi];
                let (x_bits, y_bits) = if group.is_z_only() {
                    (0, 0)
                } else {
                    group.rotation_masks()
                };
                (
                    group
                        .term_indices
                        .iter()
                        .map(|&i| (masks[i].0 | masks[i].1) as u64)
                        .collect(),
                    group.term_indices.iter().map(|&i| terms[i].0).collect(),
                    x_bits as u64,
                    y_bits as u64,
                )
            })
            .collect();
        let on_device = match host {
            Some(_) => None,
            None => backend.group_moments_on_device(&requests).transpose()?,
        };
        let moments = match on_device {
            Some(moments) => moments,
            None => {
                let exported;
                let (state, norm): (&[Complex64], f64) = match host {
                    Some(host) => host,
                    None => {
                        exported = backend.export_statevector()?;
                        (&exported, crate::backend::state_norm_sqr(&exported))
                    }
                };
                host_group_moments(state, norm, &requests)
            }
        };
        for (&gi, (m1, m2)) in deferred.iter().zip(moments) {
            mean += m1;
            group_variances[gi] = (m2 - m1 * m1).max(0.0);
        }
    }

    let variance = group_variances.iter().sum();
    Ok(ObservableExpectation {
        mean,
        variance: Some(variance),
        group_variances: Some(group_variances),
        std_error: None,
        metadata,
    })
}

/// The host moments pass over each `(zmasks, coefficients, x_bits, y_bits)`
/// request of [`grouped_expectation_on_state`].
pub(crate) fn host_group_moments(
    state: &[Complex64],
    norm: f64,
    requests: &[(Vec<u64>, Vec<f64>, u64, u64)],
) -> Vec<(f64, f64)> {
    let mut rotated = Vec::new();
    requests
        .iter()
        .map(|(zmasks, coefficients, x_bits, y_bits)| {
            let zmasks: Vec<usize> = zmasks.iter().map(|&z| z as usize).collect();
            let rotation = x_bits | y_bits;
            if rotation == 0 {
                return observable::weighted_group_moments(state, &zmasks, coefficients, norm);
            }
            observable::rotate_to_z_basis(state, &mut rotated, *x_bits as usize, *y_bits as usize);
            let growth = 2f64.powi(rotation.count_ones() as i32);
            observable::weighted_group_moments(&rotated, &zmasks, coefficients, norm * growth)
        })
        .collect()
}

/// Pair-expansion budget per commuting group. The dedicated pass pays a scratch
/// copy and a basis rotation per group that the pair expansion avoids. On the
/// 2000-string Jordan-Wigner fixture, sending every multi-term group through it
/// cost about 3x this budget's reduction at n=16 and n=20, and about 0.6x at
/// n=12.
const MAX_PAIR_MASKS_PER_GROUP: usize = 20;

/// Values from a route that evaluates rather than samples, so there is no
/// interval to report.
fn analytic_expectations(values: Vec<f64>, metadata: RunMetadata) -> ExpectationResult {
    ExpectationResult {
        values,
        std_errors: None,
        metadata,
    }
}

/// Which state diagnostic a terminal is asking for, carried into dispatch so
/// the route can be judged before the circuit runs.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum Diagnostic {
    StateVector,
    ReducedDensityMatrix,
    Entropy,
    Overlap,
}

impl Diagnostic {
    fn terminal(self) -> &'static str {
        match self {
            Diagnostic::StateVector => "a statevector",
            Diagnostic::ReducedDensityMatrix => "a reduced density matrix",
            Diagnostic::Entropy => "entanglement entropy",
            Diagnostic::Overlap => "a state overlap",
        }
    }
}

/// Whether the backend `plan` builds implements `diagnostic` on its own
/// representation.
///
/// Read from the plan rather than from a built state so [`BackendKind::Auto`]
/// can re-route before paying for the run. A plan added without the kernel
/// reads as declining here, which is what its trait default does.
fn plan_answers(plan: &BackendPlan, diagnostic: Diagnostic) -> bool {
    match diagnostic {
        Diagnostic::ReducedDensityMatrix => matches!(
            plan,
            BackendPlan::Statevector { .. }
                | BackendPlan::Sparse
                | BackendPlan::Factored
                | BackendPlan::ProductState
                | BackendPlan::DensityMatrix { .. }
                | BackendPlan::Stabilizer { .. }
                | BackendPlan::FactoredStabilizer
        ),
        Diagnostic::Entropy => matches!(
            plan,
            BackendPlan::Statevector { .. }
                | BackendPlan::Mps { .. }
                | BackendPlan::ProductState
                | BackendPlan::Stabilizer { .. }
                | BackendPlan::FactoredStabilizer
        ),
        Diagnostic::StateVector | Diagnostic::Overlap => {
            !matches!(plan, BackendPlan::DensityMatrix { .. })
        }
    }
}

/// The error for a route that propagates an observable and holds no state,
/// naming the route so the caller knows which engine `kind` reached.
fn stateless_route(kind: &BackendKind, diagnostic: Diagnostic, route: &str) -> PrismError {
    PrismError::IncompatibleBackend {
        backend: format!("{kind:?}"),
        reason: format!(
            "{} needs a backend that holds a state; the {route} route returns \
             probabilities only",
            diagnostic.terminal()
        ),
    }
}

/// Price a diagnostic's output against the dense export cap before the circuit
/// runs: a reduced density matrix holds `4^k` entries, which depends on `k`
/// and the backend's name alone, both known once the plan is built.
fn check_diagnostic_width(backend: &dyn Backend, diagnostic: Diagnostic, k: usize) -> Result<()> {
    match diagnostic {
        Diagnostic::ReducedDensityMatrix => {
            crate::backend::reduced_density::reduced_density_side(backend.name(), k)?;
            Ok(())
        }
        Diagnostic::StateVector => {
            if k > crate::backend::schmidt::export_cap() {
                return Err(crate::backend::schmidt::export_cap_exceeded(
                    backend.name(),
                    "a dense statevector".to_string(),
                    k,
                ));
            }
            Ok(())
        }
        Diagnostic::Entropy | Diagnostic::Overlap => Ok(()),
    }
}

/// Build the backend `kind` resolves to, run `circuit` on it, and hand it back
/// for one terminal read of its state.
///
/// Resolution goes straight to a single backend, as the native expectation
/// path does: a diagnostic is read off one state, and the decomposed route
/// holds one per independent block. Under [`BackendKind::Auto`] a resolved
/// plan that cannot answer `diagnostic` is replaced by the statevector while
/// the circuit fits its cap, since the route was the dispatcher's choice and
/// not the caller's; an explicit kind keeps its backend and declines.
fn diagnostic_backend(
    kind: &BackendKind,
    circuit: &Circuit,
    initial_state: Option<StartState<'_>>,
    seed: u64,
    diagnostic: Diagnostic,
    subsystem_len: usize,
) -> Result<Box<dyn Backend>> {
    if let Some(state) = initial_state {
        let mut backend = backend_from_initial_state(kind, circuit, state, seed)?;
        check_diagnostic_width(&*backend, diagnostic, subsystem_len)?;
        apply_fused_without_saves(&mut *backend, circuit, "state diagnostics")?;
        return Ok(backend);
    }
    if !kind.is_auto() {
        validate_explicit_backend(kind, circuit)?;
    }
    let (_, has_partial_independence) = analyze_independence(circuit);
    let mut plan = match resolve(kind, circuit, has_partial_independence) {
        ExecutionPlan::Backend(plan) => plan,
        ExecutionPlan::StabilizerRank => {
            return Err(stateless_route(kind, diagnostic, "stabilizer-rank"));
        }
        ExecutionPlan::StochasticPauli { .. } => {
            return Err(stateless_route(kind, diagnostic, "stochastic Pauli"));
        }
        ExecutionPlan::DeterministicPauli { .. } => {
            return Err(stateless_route(kind, diagnostic, "deterministic Pauli"));
        }
        ExecutionPlan::PauliPath => {
            return Err(stateless_route(kind, diagnostic, "Pauli path"));
        }
    };
    if kind.is_auto()
        && !plan_answers(&plan, diagnostic)
        && circuit.num_qubits <= max_statevector_qubits()
    {
        plan = plan_for_family(kind, Family::Statevector, circuit.num_qubits);
    }
    let mut backend: Box<dyn Backend> = plan.build(seed);
    check_diagnostic_width(&*backend, diagnostic, subsystem_len)?;
    execute(&mut *backend, circuit, &SimOptions::classical_only())?;
    Ok(backend)
}

/// Width from which an `Auto` expectation terminal tries the scalar tensor
/// path first. Below it the statevector run takes a few milliseconds and the
/// bounded dry run that decides the route would cost as much as the run.
const TENSOR_ROUTE_MIN_QUBITS: usize = 18;

/// The scalar tensor route for an `Auto` expectation terminal under the
/// statevector cap: `None` below [`TENSOR_ROUTE_MIN_QUBITS`], on a circuit the
/// scalar network cannot hold, or when any observable's bounded greedy plan
/// crosses the peak bound, so the caller takes its dense route. Observables
/// are validated first, so a bad one errors before any dry run.
fn tensor_route_expectations(
    circuit: &Circuit,
    observables: &[Vec<PauliTerm>],
) -> Option<Result<ExpectationResult>> {
    if circuit.num_qubits < TENSOR_ROUTE_MIN_QUBITS || observables.is_empty() {
        return None;
    }
    for observable in observables {
        if let Err(e) = validate_observable(observable, circuit.num_qubits) {
            return Some(Err(e));
        }
    }
    let planned = match crate::backend::tensornetwork::plan_bounded_expectations(
        circuit,
        observables,
        crate::backend::tensornetwork::AUTO_EXPECTATION_PEAK_BOUND,
    )? {
        Ok(planned) => planned,
        Err(e) => return Some(Err(e)),
    };
    Some(planned.contract().map(|values| {
        analytic_expectations(values, RunMetadata::exact(ResolvedBackend::TensorNetwork))
    }))
}

/// MPS work per observable, in [`mps_apply_cost`] units, from which an `Auto`
/// expectation above the statevector cap plans the scalar tensor route before
/// running MPS. About 2 ms of MPS time, the most one observable's dry run took
/// across the calibration circuits (0.1 to 1.9 ms), so a declined or losing
/// plan costs at most about the MPS run it tried to replace.
const TENSOR_ROUTE_MPS_COST_PER_OBSERVABLE: f64 = 1.5e5;

/// [`mps_apply_cost`] units one planned tensor-network element costs: about
/// 40 ns per element contracted against 15 ns per unit.
const TENSOR_ELEMENT_MPS_COST: f64 = 3.0;

/// The scalar tensor route for an `Auto` expectation above the statevector cap,
/// where the alternative is the MPS `Auto` resolves to rather than a dense run.
/// `None` when `Auto` resolves to another backend, when the MPS estimate is
/// under [`TENSOR_ROUTE_MPS_COST_PER_OBSERVABLE`] per observable, when a plan
/// crosses the peak bound, or when the planned contraction costs more than the
/// MPS estimate.
fn tensor_route_over_mps(
    circuit: &Circuit,
    observables: &[Vec<PauliTerm>],
) -> Option<Result<ExpectationResult>> {
    if observables.is_empty() {
        return None;
    }
    let mps_cost = mps_apply_cost(circuit, AUTO_MPS_BOND_DIM);
    if mps_cost < TENSOR_ROUTE_MPS_COST_PER_OBSERVABLE * observables.len() as f64 {
        return None;
    }
    let (_, has_partial_independence) = analyze_independence(circuit);
    let ExecutionPlan::Backend(plan) =
        resolve(&BackendKind::Auto, circuit, has_partial_independence)
    else {
        return None;
    };
    if !matches!(plan, BackendPlan::Mps { .. }) {
        return None;
    }
    for observable in observables {
        if let Err(e) = validate_observable(observable, circuit.num_qubits) {
            return Some(Err(e));
        }
    }
    let planned = match crate::backend::tensornetwork::plan_bounded_expectations(
        circuit,
        observables,
        crate::backend::tensornetwork::AUTO_EXPECTATION_PEAK_BOUND,
    )? {
        Ok(planned) => planned,
        Err(e) => return Some(Err(e)),
    };
    if planned.total_elements() as f64 * TENSOR_ELEMENT_MPS_COST >= mps_cost {
        return None;
    }
    Some(planned.contract().map(|values| {
        analytic_expectations(values, RunMetadata::exact(ResolvedBackend::TensorNetwork))
    }))
}

/// Evaluate `observables` on the backend `kind` resolves to, using that
/// backend's own representation.
///
/// Backends without a native Pauli path report `BackendUnsupported` naming
/// themselves, so a request that cannot be served says which engine could not
/// serve it rather than blaming the route that picked it.
fn expectation_values_native(
    kind: &BackendKind,
    circuit: &Circuit,
    observables: &[Vec<PauliTerm>],
    seed: u64,
) -> Result<ExpectationResult> {
    if !kind.is_auto() {
        validate_explicit_backend(kind, circuit)?;
    }
    // Before the run, matching the statevector path, so a typo in an observable
    // does not cost a 40-qubit simulation first.
    for observable in observables {
        validate_observable(observable, circuit.num_qubits)?;
    }

    let (_, has_partial_independence) = analyze_independence(circuit);
    let ExecutionPlan::Backend(plan) = resolve(kind, circuit, has_partial_independence) else {
        return Err(PrismError::IncompatibleBackend {
            backend: format!("{kind:?}"),
            reason: "expectation values need a backend that holds a state; the stabilizer-rank \
                     route returns probabilities only"
                .into(),
        });
    };

    let mut backend = plan.build(seed);
    if !backend.supports_pauli_expectation() {
        return Err(PrismError::BackendUnsupported {
            backend: backend.name().to_string(),
            operation: "Pauli expectation values".to_string(),
        });
    }
    execute(&mut *backend, circuit, &SimOptions::classical_only())?;
    let values = backend.pauli_expectations(observables)?;
    Ok(analytic_expectations(values, backend_metadata(&*backend)))
}

fn expectation_values_statevector(
    kind: &BackendKind,
    circuit: &Circuit,
    observables: &[Vec<PauliTerm>],
    seed: u64,
) -> Result<ExpectationResult> {
    // Validate before the 2^n simulation so bad observables fail cheaply. The
    // mask reduction shifts by qubit index, so it follows the `init` that
    // rejects a width no mask can address.
    for obs in observables {
        validate_observable(obs, circuit.num_qubits)?;
    }

    let accel = accel_for(kind, Family::Statevector, circuit.num_qubits);
    let mut backend = build_statevector(&accel, seed);
    let expanded = expand_for_backend(&backend, circuit);
    let fused = fuse_for_backend(&backend, &expanded);
    backend.init(fused.num_qubits, fused.num_classical_bits)?;
    let masks = observables
        .iter()
        .map(|obs| pauli_masks(obs, circuit.num_qubits))
        .collect::<Result<Vec<_>>>()?;
    backend.apply_instructions(&fused.instructions)?;

    let values = match backend.pauli_expectations_on_device(&masks) {
        Some(values) => values?,
        None => pauli_expectations_from_masks(backend.state_vector(), &masks).0,
    };
    let metadata = backend_metadata(&backend);
    Ok(analytic_expectations(values, metadata))
}

/// Multi-shot execution for the distributed statevector backend.
///
/// Every rank runs this function in lockstep. Circuits with only terminal
/// measurements run once and sample basis indices without gathering the dense
/// state on any rank. Circuits with mid-circuit measurements run once per shot,
/// prefused, with per-shot seeds matching the generic slow path.
#[cfg(feature = "distributed")]
fn run_shots_distributed(
    context: std::sync::Arc<crate::distributed::DistributedContext>,
    circuit: &Circuit,
    num_shots: usize,
    seed: u64,
) -> Result<ShotsResult> {
    use crate::backend::distributed_statevector::DistributedStatevectorBackend;

    let meas_map = circuit.measurement_map();
    if meas_map.is_empty() {
        // No measurements means every shot is all false, but init must still
        // run so invalid rank counts and local qubit floor violations surface as
        // errors instead of fabricated output.
        let mut backend = DistributedStatevectorBackend::new(context, seed);
        backend.init(circuit.num_qubits, circuit.num_classical_bits)?;
        return Ok(ShotsResult::from_shots(
            vec![vec![false; circuit.num_classical_bits]; num_shots],
            circuit.num_classical_bits,
        )
        .with_metadata(backend_metadata(&backend)));
    }

    if circuit.has_terminal_measurements_only() {
        let stripped = circuit.without_measurements();
        let mut backend = DistributedStatevectorBackend::new(context, seed);
        execute(&mut backend, &stripped, &SimOptions::classical_only())?;
        let samples = backend.sample_basis_states(num_shots, seed)?;
        return Ok(ShotsResult::from_shots(
            shots_from_basis_samples(&samples, &meas_map, circuit.num_classical_bits),
            circuit.num_classical_bits,
        )
        .with_metadata(backend_metadata(&backend)));
    }

    let probe = DistributedStatevectorBackend::new(context.clone(), seed);
    let expanded = expand_for_backend(&probe, circuit);
    let fused = fuse_for_backend(&probe, &expanded);
    let opts = SimOptions::classical_only();
    let mut shots = Vec::with_capacity(num_shots);
    let mut metadata = RunMetadata::exact(ResolvedBackend::Distributed);
    for i in 0..num_shots {
        let shot_seed = mix_seed(seed, i);
        let mut backend = DistributedStatevectorBackend::new(context.clone(), shot_seed);
        let result = execute_circuit(&mut backend, &fused, &opts)?;
        metadata.weaken_with(&result.metadata);
        shots.push(result.classical_bits);
    }
    Ok(ShotsResult::from_shots(shots, circuit.num_classical_bits).with_metadata(metadata))
}

/// Execute a circuit multiple times with explicit backend selection.
pub(crate) fn run_shots_with(
    kind: BackendKind,
    circuit: &Circuit,
    num_shots: usize,
    seed: u64,
) -> Result<ShotsResult> {
    // The distributed backend runs every rank in lockstep, so shot execution
    // must not route through shortcuts that reshape the collective call
    // sequence. Dispatch directly.
    #[cfg(feature = "distributed")]
    if let BackendKind::StatevectorDistributed { context } = &kind {
        return run_shots_distributed(context.clone(), circuit, num_shots, seed);
    }

    // Once per shots call, not per shot: a guard that cannot depend on a
    // measurement is resolved here so the sampling predicates below see the
    // circuit that actually runs. A circuit with no guard borrows through.
    let folded = circuit.fold_static_guards();
    let circuit = folded.as_ref();

    let bits = circuit.num_classical_bits;
    let source = prepare_shot_source(&kind, circuit, num_shots, seed)?;
    let Some(metadata) = source.metadata() else {
        return run_shots_per_shot(kind, circuit, num_shots, seed, ShotRows::Exchangeable);
    };
    let result = match source {
        ShotSource::Compiled {
            mut sampler,
            meas_map,
            ..
        } => {
            let packed = sampler.try_sample_bulk_packed(num_shots)?;
            ShotsResult::from_shots(
                packed_shots_to_classical_bits(&packed, &meas_map, bits),
                bits,
            )
        }
        ShotSource::TerminalStatevector { backend, meas_map } => {
            let shots = if backend.is_gpu_resident() {
                let probs = backend.probabilities()?;
                sample_shots_from_probs(&probs, &meas_map, bits, num_shots, seed)
            } else {
                sample_shots_from_state(
                    backend.state_vector(),
                    backend.probability_scale(),
                    &meas_map,
                    bits,
                    num_shots,
                    seed,
                )
            };
            ShotsResult::from_shots(shots, bits)
        }
        ShotSource::Native {
            mut backend,
            meas_map,
        } => {
            let samples = backend.sample_basis_states(num_shots, seed)?;
            ShotsResult::from_shots(shots_from_basis_samples(&samples, &meas_map, bits), bits)
        }
        ShotSource::TerminalProbabilities {
            probs, meas_map, ..
        } => ShotsResult::from_shots(sample_shots(&probs, &meas_map, bits, num_shots, seed), bits),
        ShotSource::StabilizerRank => {
            stabilizer_rank::run_stabilizer_rank_shots(circuit, num_shots, seed)?
        }
        ShotSource::PerShot => unreachable!("handled above"),
    };
    Ok(result.with_metadata(metadata))
}

/// Whether a caller of [`run_shots_per_shot`] reads the rows themselves or only
/// counts them, which spares the shuffle that keeps branched rows exchangeable.
#[derive(Clone, Copy, PartialEq, Eq)]
enum ShotRows {
    Exchangeable,
    CountedOnly,
}

/// Run `circuit` once per shot, which mid-circuit measurements force, or on
/// the host statevector once per outcome history
/// ([`branching::branch_shots`]).
fn run_shots_per_shot(
    kind: BackendKind,
    circuit: &Circuit,
    num_shots: usize,
    seed: u64,
    rows: ShotRows,
) -> Result<ShotsResult> {
    // Pre-compute seed-independent analysis to avoid redundant work.
    if !kind.is_auto() {
        validate_explicit_backend(&kind, circuit)?;
    }

    let (decompose, has_partial_independence) = analyze_independence(circuit);

    if matches!(kind, BackendKind::StabilizerRank) {
        return stabilizer_rank::run_stabilizer_rank_shots(circuit, num_shots, seed);
    }
    if matches!(
        kind,
        BackendKind::StochasticPauli { .. } | BackendKind::DeterministicPauli { .. }
    ) {
        return Err(crate::error::PrismError::IncompatibleBackend {
            backend: format!("{kind:?}"),
            reason: "Pauli propagation backends do not support mid-circuit measurements".into(),
        });
    }
    if kind.is_auto() && auto_stabilizer_rank_t_count(circuit, MAX_AUTO_T_COUNT_SHOTS).is_some() {
        return stabilizer_rank::run_stabilizer_rank_shots(circuit, num_shots, seed);
    }

    if has_temporal_clifford_opportunity(&kind, circuit) {
        if decompose.is_none()
            && let Some(tc) = plan_temporal_clifford(&kind, circuit)
        {
            let route = ResolvedBackend::Statevector;
            let states = [(route, circuit.num_qubits)];
            return collect_shots(
                circuit,
                num_shots,
                seed,
                route,
                &kind,
                &states,
                |shot_seed| {
                    let outcome = run_temporal_clifford(&tc, shot_seed, false)?;
                    Ok((outcome.classical_bits, outcome.metadata))
                },
            );
        }
        // Decomposable circuits with a temporal prefix keep the per-shot
        // full-pipeline route; the prefix spans blocks that decomposition
        // would otherwise split. Its blocks are partitioned and planned once.
        let opts = SimOptions::classical_only();
        let route = resolve_backend(&kind, circuit, has_partial_independence).resolved();
        let plan = plan_probability_route(&kind, circuit);
        let blocks = match &plan {
            ProbabilityRoute::Decomposed(components) => {
                Some(PreparedBlocks::new(&kind, components.clone(), circuit))
            }
            _ => None,
        };
        let states = [(route, circuit.num_qubits)];
        return collect_shots(
            circuit,
            num_shots,
            seed,
            route,
            &kind,
            &states,
            |shot_seed| {
                let outcome = match &blocks {
                    Some(blocks) => blocks.run(&kind, circuit, shot_seed, &opts)?,
                    None => run_route(&kind, circuit, shot_seed, opts, &plan)?,
                };
                Ok((outcome.classical_bits, outcome.metadata))
            },
        );
    }

    let opts = SimOptions::classical_only();

    if let Some((ref comps, _)) = decompose {
        let partitions = circuit.partition_subcircuits(comps);
        let block_plans: Vec<BackendPlan> = partitions
            .iter()
            .map(|(sub, _, _)| {
                if !kind.is_auto() {
                    validate_explicit_backend(&kind, sub)?;
                }
                Ok(resolve_backend(&kind, sub, false))
            })
            .collect::<Result<_>>()?;
        let fused_blocks: Vec<std::borrow::Cow<'_, Circuit>> = partitions
            .iter()
            .zip(&block_plans)
            .map(|((sub, _, _), plan)| {
                let probe = plan.build(seed);
                let expanded = expand_for_backend(&*probe, sub);
                std::borrow::Cow::Owned(fuse_for_backend(&*probe, &expanded).into_owned())
            })
            .collect();

        let block_states: Vec<(ResolvedBackend, usize)> = block_plans
            .iter()
            .zip(comps)
            .map(|(plan, comp)| (plan.resolved(), comp.len()))
            .collect();
        collect_shots(
            circuit,
            num_shots,
            seed,
            ResolvedBackend::Decomposed,
            &kind,
            &block_states,
            |shot_seed| {
                let result = run_decomposed_prefused(
                    &block_plans,
                    comps,
                    &partitions,
                    &fused_blocks,
                    shot_seed,
                    &opts,
                    circuit,
                )?;
                Ok((result.classical_bits, result.metadata))
            },
        )
    } else {
        let plan = resolve_backend(&kind, circuit, has_partial_independence);
        let probe = plan.build(seed);
        let expanded = expand_for_backend(&*probe, circuit);
        let fused = fuse_for_backend(&*probe, &expanded);

        let route = plan.resolved();
        let states = [(route, circuit.num_qubits)];
        if plan.is_host_statevector() && circuit.num_qubits < max_statevector_qubits() {
            #[cfg(feature = "parallel")]
            let parallel = shots_split_across_workers(&kind, &states)
                && circuit.num_qubits >= branching::MIN_QUBITS_FOR_PARALLEL_SPLITS;
            #[cfg(not(feature = "parallel"))]
            let parallel = false;
            let limits = branching::BranchLimits::for_width(circuit.num_qubits, parallel);
            return branching::branch_shots(
                &fused,
                num_shots,
                seed,
                &limits,
                rows == ShotRows::Exchangeable,
            );
        }
        collect_shots(
            circuit,
            num_shots,
            seed,
            route,
            &kind,
            &states,
            |shot_seed| {
                let mut backend = plan.build(shot_seed);
                let outcome = execute_circuit(&mut *backend, &fused, &opts)?;
                Ok((outcome.classical_bits, outcome.metadata))
            },
        )
    }
}

pub(crate) const fn splitmix64(mut z: u64) -> u64 {
    z = z.wrapping_add(0x9E37_79B9_7F4A_7C15);
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    z ^ (z >> 31)
}

/// Seed for item `index` (a shot, trajectory, block, qubit or observable) of a run
/// seeded `seed`. Both inputs pass through SplitMix64, so runs on nearby seeds draw
/// unrelated streams rather than shifted copies of one another.
pub(crate) const fn mix_seed(seed: u64, index: usize) -> u64 {
    splitmix64(splitmix64(seed) ^ index as u64)
}

/// Each shot evolves its own state, so `shot` returns the provenance of its own
/// run and the ensemble keeps the weakest claim across them. `route` names the
/// engine for a zero-shot request, which runs nothing to read provenance off.
///
/// Under the `parallel` feature the shots split across Rayon workers when
/// [`shots_split_across_workers`] accepts `kind` for `states`, the family and width
/// of each state one shot holds. Shot `i` always runs on `mix_seed(seed, i)` and
/// the fold runs in shot order, so the result does not depend on the thread
/// count.
fn collect_shots(
    circuit: &Circuit,
    num_shots: usize,
    seed: u64,
    route: ResolvedBackend,
    kind: &BackendKind,
    states: &[(ResolvedBackend, usize)],
    shot: impl Fn(u64) -> Result<(Vec<bool>, RunMetadata)> + Sync,
) -> Result<ShotsResult> {
    #[cfg(feature = "parallel")]
    if num_shots > 1 && shots_split_across_workers(kind, states) {
        use rayon::prelude::*;
        let runs: Vec<Result<(Vec<bool>, RunMetadata)>> = (0..num_shots)
            .into_par_iter()
            .map(|i| shot(mix_seed(seed, i)))
            .collect();
        return fold_shots(circuit, route, runs);
    }
    #[cfg(not(feature = "parallel"))]
    let _ = (kind, states);
    fold_shots(
        circuit,
        route,
        (0..num_shots).map(|i| shot(mix_seed(seed, i))),
    )
}

fn fold_shots(
    circuit: &Circuit,
    route: ResolvedBackend,
    runs: impl IntoIterator<Item = Result<(Vec<bool>, RunMetadata)>>,
) -> Result<ShotsResult> {
    let runs = runs.into_iter();
    let mut shots = Vec::with_capacity(runs.size_hint().0);
    let mut metadata = RunMetadata::exact(route);
    for (i, run) in runs.enumerate() {
        let (bits, shot_metadata) = run?;
        if i == 0 {
            metadata = shot_metadata;
        } else {
            metadata.weaken_with(&shot_metadata);
        }
        shots.push(bits);
    }
    Ok(ShotsResult::from_shots(shots, circuit.num_classical_bits).with_metadata(metadata))
}

/// Family choice for auto-routed non-Pauli noise trajectories. Restricted to
/// families whose trajectory operations (1q Kraus, qubit probability, reduced
/// density matrix, reset) are supported; the statevector leaf carries the
/// kind's acceleration. A two-qubit Kraus event entangles its pair, so it keeps
/// the circuit off the product state as an entangling gate would.
fn general_noise_plan(
    kind: &BackendKind,
    circuit: &Circuit,
    noise_model: &NoiseModel,
) -> BackendPlan {
    let family = if !circuit.has_entangling_gates() && !noise_model.has_two_qubit_kraus() {
        Family::ProductState
    } else if circuit.num_qubits > max_statevector_qubits() {
        if circuit.is_sparse_friendly() && circuit.num_qubits <= MAX_SPARSE_INDEX_QUBITS {
            Family::Sparse
        } else {
            Family::Mps
        }
    } else {
        Family::Statevector
    };
    plan_for_family(kind, family, circuit.num_qubits)
}

/// Execute a noisy circuit for multiple shots with explicit backend selection.
///
/// For Clifford circuits with Auto/Stabilizer/FactoredStabilizer backends,
/// uses the compiled noisy sampler (fast O(n²·m) compile + O(events·m/64) per shot).
/// For all other cases, falls back to per-shot simulation with noise injection.
/// The compiled noisy path is limited to terminal measurements with no resets
/// or classical conditionals. Pauli noise on the host statevector, or on a
/// tensor network whose probabilities fit the dense cap, under the same limits
/// draws every shot's errors first and simulates each distinct pattern once;
/// see [`trajectory::PauliGroups`]. On the host statevector, a circuit with a
/// mid-circuit measurement, reset or condition shares the gates before the first
/// of them the same way, and each shot finishes its own trajectory from there.
pub(crate) fn run_shots_with_noise(
    kind: BackendKind,
    circuit: &Circuit,
    noise_model: &noise::NoiseModel,
    num_shots: usize,
    seed: u64,
) -> Result<ShotsResult> {
    noise_model.validate_for(circuit)?;

    // Trajectory execution runs shots on Rayon worker threads, whose
    // scheduling order differs per rank. Per-shot distributed backends would
    // issue collectives out of lockstep and deadlock or corrupt exchanges.
    // Reject until a lockstep noisy path exists.
    #[cfg(feature = "distributed")]
    if matches!(kind, BackendKind::StatevectorDistributed { .. }) {
        return Err(crate::error::PrismError::IncompatibleBackend {
            backend: format!("{kind:?}"),
            reason: "noisy shot sampling is not supported on the distributed backend; \
                     trajectory execution cannot keep rank collectives in lockstep"
                .into(),
        });
    }

    if kind.is_density_matrix() {
        let probs = exact_noisy_probabilities(&kind, circuit, noise_model, None, seed)?;
        return Ok(ShotsResult::from_shots(
            sample_exact_noisy_shots(&probs, circuit, noise_model, num_shots, seed),
            circuit.num_classical_bits,
        )
        .with_metadata(exact_mixture_metadata(&kind)));
    }

    if !kind.supports_noisy_per_shot() {
        return Err(crate::error::PrismError::IncompatibleBackend {
            backend: format!("{kind:?}"),
            reason: "this backend holds no per-shot pure state to inject noise into; select \
                     DensityMatrix for the exact mixed state, or a backend that evolves one \
                     state per trajectory"
                .into(),
        });
    }

    let is_stabilizer_kind = kind.is_stabilizer_family();

    if is_stabilizer_kind && !noise_model.has_only_pauli_channels() {
        return Err(crate::error::PrismError::IncompatibleBackend {
            backend: format!("{kind:?}"),
            reason: format!(
                "stabilizer backends only support Pauli/depolarizing noise; use {} for amplitude damping, phase damping, thermal relaxation, or custom Kraus",
                BackendKind::general_noise_backend_names()
            ),
        });
    }

    if !noise_model.has_only_pauli_channels() && !kind.supports_general_noise() {
        return Err(crate::error::PrismError::IncompatibleBackend {
            backend: format!("{kind:?}"),
            reason: format!(
                "non-Pauli noise requires {}",
                BackendKind::general_noise_backend_names()
            ),
        });
    }

    if is_stabilizer_kind && !circuit.is_clifford_only() {
        return Err(crate::error::PrismError::IncompatibleBackend {
            backend: format!("{kind:?}"),
            reason: "circuit contains non-Clifford gates".into(),
        });
    }

    if !kind.is_auto() {
        validate_explicit_backend(&kind, circuit)?;
    }

    if noise_model.has_only_pauli_channels() {
        let use_compiled = (kind.is_auto()
            || matches!(
                kind,
                BackendKind::Stabilizer | BackendKind::FactoredStabilizer
            ))
            && supports_compiled_measurement_sampling(circuit)
            || {
                #[cfg(feature = "gpu")]
                {
                    matches!(kind, BackendKind::StabilizerGpu { .. })
                        && supports_compiled_measurement_sampling(circuit)
                }
                #[cfg(not(feature = "gpu"))]
                {
                    false
                }
            };

        if use_compiled {
            #[cfg(feature = "gpu")]
            if let BackendKind::StabilizerGpu { context } = &kind {
                return noise::run_shots_noisy_with_gpu(
                    circuit,
                    noise_model,
                    num_shots,
                    seed,
                    context.clone(),
                );
            }
            return noise::run_shots_noisy(circuit, noise_model, num_shots, seed);
        }
    }

    let plan = if kind.is_auto() && !noise_model.has_only_pauli_channels() {
        general_noise_plan(&kind, circuit, noise_model)
    } else {
        resolve_backend(&kind, circuit, false)
    };
    // Noise events are indexed per instruction, so trajectories apply the
    // stream raw and the Pauli-rotation lowering pass cannot run. The
    // statevector applies the gate natively (its device path lowers inline);
    // any other backend without the kernel is rejected before a shot starts.
    let has_pauli_rot = crate::circuit::any_gate(&circuit.instructions, &mut |gate| {
        matches!(gate, crate::gates::Gate::PauliRot(_))
    });
    if has_pauli_rot
        && !matches!(plan, BackendPlan::Statevector { .. })
        && !plan.build(seed).supports_pauli_rotation()
    {
        return Err(crate::error::PrismError::IncompatibleBackend {
            backend: format!("{:?}", plan.resolved()),
            reason: "noisy trajectories apply the instruction stream raw so noise events \
                     stay aligned to it, which leaves no room for the Pauli-rotation \
                     lowering this backend needs; run on the statevector, or expand the \
                     rotations with circuit::expand_pauli_rotations and attach the noise \
                     model to the expanded circuit"
                .into(),
        });
    }
    // A correlated two-qubit Kraus branch reads a two-qubit reduced density
    // matrix, which only the host statevector answers. Without this the model
    // clears every routing gate and fails part way through the first shot, on a
    // backend an Auto route may have picked for the caller.
    if noise_model.has_two_qubit_kraus() && !plan.build(seed).supports_two_qubit_kraus() {
        return Err(crate::error::PrismError::IncompatibleBackend {
            backend: format!("{:?}", plan.resolved()),
            reason: "a two-qubit Kraus channel needs the pair's reduced density matrix \
                     its branch probabilities are drawn from, and a kernel for the branch \
                     operator, which is not unitary. A tableau and a product state have no \
                     such kernel, and a tensor network traces out no arbitrary pair; run on \
                     a statevector, sparse, factored or MPS route, or evaluate the channel \
                     exactly on BackendKind::DensityMatrix"
                .into(),
        });
    }
    let host_statevector = matches!(plan, BackendPlan::Statevector { .. }) && !plan.is_gpu();
    let dense_tensor_network = matches!(plan, BackendPlan::TensorNetwork { .. })
        && crate::backend::tensor_probability_len("tensor network", circuit.num_qubits).is_ok();
    if (host_statevector || dense_tensor_network)
        && let Some(groups) = trajectory::PauliGroups::sample(circuit, noise_model, num_shots, seed)
    {
        let build = |s| plan.build(s);
        return trajectory::run_pauli_groups(
            &groups,
            circuit,
            noise_model,
            seed,
            dense_tensor_network.then_some(&build as trajectory::GroupBackendFactory<'_>),
            plan.resolved(),
        );
    }
    // A checkpoint, a group state and a shot state at once.
    if host_statevector
        && circuit.num_qubits + 1 < max_statevector_qubits()
        && let Some(groups) =
            trajectory::PauliGroups::sample_replayed(circuit, noise_model, num_shots, seed)
    {
        return trajectory::run_replayed_groups(
            &groups,
            circuit,
            noise_model,
            seed,
            plan.resolved(),
        );
    }
    let route = plan.resolved();
    trajectory::run_trajectories(
        |s| plan.build(s),
        circuit,
        noise_model,
        num_shots,
        seed,
        plan.is_gpu(),
        route,
    )
}

#[cfg(test)]
mod tests;

#[cfg(all(test, feature = "gpu"))]
mod gpu_stub_tests;

#[cfg(test)]
mod terminal_candidate_matrix_tests;

#[cfg(test)]
mod diagnostic_terminal_tests;