prism-q 0.35.0

//! Noise channels, noise models, and noisy samplers. Pauli-only models on
//! Clifford circuits compile into a [`NoisyCompiledSampler`]; general
//! channels run through the trajectory engine.

use num_complex::Complex64;
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;
use smallvec::smallvec;

use crate::backend::Backend;
use crate::backend::density_matrix::{
    DeferredSuperoperators, DensityMatrixBackend, MIN_DEFERRED_QUBITS, PairChannel,
};
use crate::backend::stabilizer::StabilizerBackend;
use crate::circuit::{Circuit, Instruction, SmallVec};
use crate::error::Result;
use crate::gates::Gate;
use crate::sim::compiled::{
    PackedShots, batch_propagate_backward, compile_measurements, xor_words,
};
use crate::sim::unified_pauli::PauliAxis;
use crate::sim::{BackendKind, ShotsResult};

/// A noise channel on one qubit, or on two for `TwoQubitDepolarizing` and
/// `Kraus2q`.
///
/// Fires after each instruction whose `NoiseModel::after_gate` slot holds a
/// `NoiseEvent` carrying it. On a Clifford circuit, channels that sample as a
/// Pauli frame draw take the compiled stabilizer sampler; see
/// [`NoiseChannel::is_pauli`] and
/// [`NoiseChannel::TwoQubitDepolarizing`].
#[derive(Debug, Clone, PartialEq)]
#[non_exhaustive]
pub enum NoiseChannel {
    /// Independent Pauli X/Y/Z error with given per-branch probabilities.
    Pauli { px: f64, py: f64, pz: f64 },
    /// Symmetric depolarizing: `I` with probability `1-p`, each of `X/Y/Z`
    /// with probability `p/3`.
    Depolarizing { p: f64 },
    /// T1 relaxation: `|1⟩ → |0⟩` with amplitude `sqrt(gamma)`.
    AmplitudeDamping { gamma: f64 },
    /// Pure dephasing: phase randomisation with amplitude `sqrt(gamma)` on the
    /// `|1⟩` state. Populations are preserved.
    PhaseDamping { gamma: f64 },
    /// Combined T1 + T2 relaxation over a gate of duration `gate_time`,
    /// requiring `t2 <= 2·t1`. Populations decay as `exp(-gate_time/t1)` and
    /// coherences as `exp(-gate_time/t2)`; both the exact and the trajectory
    /// route realize that as generalized amplitude damping composed with pure
    /// dephasing.
    ///
    /// `excited_population` is the steady state the qubit relaxes toward: 0 for
    /// a cold qubit that settles in the ground state, 0.5 for one whose steady
    /// state is maximally mixed. A device's value is
    /// `1 / (1 + exp(h f / k T))` at frequency `f` and temperature `T`.
    ThermalRelaxation {
        t1: f64,
        t2: f64,
        gate_time: f64,
        excited_population: f64,
    },
    /// Symmetric two-qubit depolarizing: each of the 15 non-identity
    /// two-qubit Paulis occurs with probability `p/15`.
    TwoQubitDepolarizing { p: f64 },
    /// General single-qubit channel described by a set of Kraus operators.
    ///
    /// Dense Kraus operators are supported by trajectory backends that expose
    /// `Backend::reduced_density_matrix_1q`. Branch probabilities use the full
    /// one-qubit reduced density matrix, so channels may depend on coherence
    /// as well as populations. Call [`NoiseModel::validate`] in setup code to
    /// catch empty or non-finite custom channels before sampling.
    Custom { kraus: Vec<[[Complex64; 2]; 2]> },
    /// General two-qubit channel described by a set of 4x4 Kraus operators.
    ///
    /// Each operator is indexed `K[t][t']` with `t = 2 * bit(qubits[0]) +
    /// bit(qubits[1])`, the packing [`crate::gates::Gate::matrix_4x4`] uses, so
    /// the first target is the high bit. Runs exactly on the density matrix and
    /// on trajectory backends reporting [`Backend::supports_two_qubit_kraus`].
    ///
    /// Correlated `ZZ` dephasing of strength `p`, as the diagonal pair
    /// `sqrt(1-p) I` and `sqrt(p) (Z (x) Z)`:
    ///
    /// ```
    /// use num_complex::Complex64;
    /// use prism_q::NoiseChannel;
    ///
    /// let p = 0.02_f64;
    /// let diag = |scale: f64, signs: [f64; 4]| {
    ///     let mut m = [[Complex64::new(0.0, 0.0); 4]; 4];
    ///     for (t, sign) in signs.iter().enumerate() {
    ///         m[t][t] = Complex64::new(scale * sign, 0.0);
    ///     }
    ///     m
    /// };
    /// let zz = NoiseChannel::Kraus2q {
    ///     kraus: vec![
    ///         diag((1.0 - p).sqrt(), [1.0, 1.0, 1.0, 1.0]),
    ///         diag(p.sqrt(), [1.0, -1.0, -1.0, 1.0]),
    ///     ],
    /// };
    /// zz.validate()?;
    /// # Ok::<(), prism_q::PrismError>(())
    /// ```
    Kraus2q { kraus: Vec<[[Complex64; 4]; 4]> },
    /// With probability `p`, an unleaked qubit leaves the computational subspace.
    ///
    /// A leaked qubit is a classical flag in the trajectory: gates on it do not
    /// act, a gate pairing it with unleaked qubits applies a uniformly random
    /// Pauli (`I`, `X`, `Y`, `Z` at a quarter each) to each of them instead,
    /// measuring it reports 1, and a reset or a [`NoiseChannel::Seepage`] clears
    /// the flag. Every other channel skips an event naming a leaked qubit.
    /// [`ShotsResult::leaked`](crate::ShotsResult::leaked) records which qubits
    /// leaked in each shot. Runs on the per-shot trajectory engines only.
    Leakage { p: f64 },
    /// With probability `p`, a leaked qubit returns to the computational
    /// subspace in a uniformly random basis state, `|0>` or `|1>` at one half
    /// each and uncorrelated with the rest of the register.
    Seepage { p: f64 },
    /// With probability `p`, leakage spreads across a pair: when exactly one of
    /// the two qubits is leaked as the event fires, the other leaks too.
    ///
    /// The flags are read when the event fires, so declare it ahead of a
    /// [`NoiseChannel::Leakage`] event in the same slot to spread only the
    /// leakage the gate itself saw.
    LeakageTransport { p: f64 },
    /// Coherent rotation `exp(-i angle P / 2)` about `axis`, with an angle
    /// drawn once per shot as `sum_k weight_k * z[source_k]` over the pairs in
    /// `weights`.
    ///
    /// `z` holds independent standard normal draws, one per source index, made
    /// at the start of each shot from its noise stream and shared by every event
    /// of the shot, so two events naming one source rotate by correlated angles.
    /// That makes each shot a different unitary, which only the per-shot
    /// trajectory engines can sample; [`NoiseBuilder`] lowers detuning and
    /// amplitude drift onto it.
    QuasiStatic {
        axis: PauliAxis,
        weights: Vec<(usize, f64)>,
    },
}

impl NoiseChannel {
    /// Number of qubits an event carrying this channel must name.
    pub fn num_qubits(&self) -> usize {
        match self {
            NoiseChannel::TwoQubitDepolarizing { .. }
            | NoiseChannel::Kraus2q { .. }
            | NoiseChannel::LeakageTransport { .. } => 2,
            _ => 1,
        }
    }

    /// Per-Pauli `(px, py, pz)` probabilities for a Pauli or Depolarizing
    /// channel, `None` for every other channel.
    pub fn as_pauli(&self) -> Option<(f64, f64, f64)> {
        match self {
            NoiseChannel::Pauli { px, py, pz } => Some((*px, *py, *pz)),
            NoiseChannel::Depolarizing { p } => {
                let pp = p / 3.0;
                Some((pp, pp, pp))
            }
            _ => None,
        }
    }

    /// Total error probability of a two-qubit Pauli channel, `None` for every
    /// other channel. The 15 non-identity products share it equally.
    pub fn pauli_pair_rate(&self) -> Option<f64> {
        match self {
            NoiseChannel::TwoQubitDepolarizing { p } => Some(*p),
            _ => None,
        }
    }

    /// True when no branch of the channel can fire, so it moves no state and
    /// decides no route, the reading [`ReadoutError::is_inert`] gives a
    /// zero-rate readout entry. Only the Pauli family is considered: a damping
    /// or Kraus channel stays outside the Pauli frame whatever its parameters.
    pub(crate) fn is_inert(&self) -> bool {
        match self {
            NoiseChannel::Pauli { px, py, pz } => *px == 0.0 && *py == 0.0 && *pz == 0.0,
            NoiseChannel::Depolarizing { p } | NoiseChannel::TwoQubitDepolarizing { p } => {
                *p == 0.0
            }
            _ => false,
        }
    }

    /// True for the Pauli and Depolarizing channels.
    pub fn is_pauli(&self) -> bool {
        matches!(
            self,
            NoiseChannel::Pauli { .. } | NoiseChannel::Depolarizing { .. }
        )
    }

    /// True for the channels a Pauli frame can carry, which is
    /// [`NoiseChannel::is_pauli`] plus
    /// [`NoiseChannel::TwoQubitDepolarizing`]. The pair channel needs a joint
    /// draw over its 15 branches rather than one per qubit, so the samplers
    /// hold it apart from the single-qubit rows.
    pub(crate) fn is_pauli_frame(&self) -> bool {
        self.is_pauli() || self.pauli_pair_rate().is_some()
    }

    /// True when `validate` accepts the channel.
    pub fn is_exactly_samplable(&self) -> bool {
        self.validate().is_ok()
    }

    /// Validate that this channel is usable by the trajectory engine.
    pub fn validate(&self) -> Result<()> {
        match self {
            NoiseChannel::Pauli { px, py, pz } => {
                validate_probability("px", *px)?;
                validate_probability("py", *py)?;
                validate_probability("pz", *pz)?;
                if *px + *py + *pz > 1.0 + 1e-12 {
                    return Err(crate::error::PrismError::InvalidParameter {
                        message: "Pauli channel probabilities must sum to <= 1".into(),
                    });
                }
            }
            NoiseChannel::Depolarizing { p } | NoiseChannel::TwoQubitDepolarizing { p } => {
                validate_probability("p", *p)?;
            }
            NoiseChannel::AmplitudeDamping { gamma } | NoiseChannel::PhaseDamping { gamma } => {
                validate_probability("gamma", *gamma)?;
            }
            NoiseChannel::ThermalRelaxation {
                t1,
                t2,
                gate_time,
                excited_population,
            } => {
                if !t1.is_finite() || *t1 <= 0.0 {
                    return Err(crate::error::PrismError::InvalidParameter {
                        message: "thermal relaxation t1 must be finite and positive".into(),
                    });
                }
                if !t2.is_finite() || *t2 <= 0.0 {
                    return Err(crate::error::PrismError::InvalidParameter {
                        message: "thermal relaxation t2 must be finite and positive".into(),
                    });
                }
                if !gate_time.is_finite() || *gate_time < 0.0 {
                    return Err(crate::error::PrismError::InvalidParameter {
                        message: "thermal relaxation gate_time must be finite and non-negative"
                            .into(),
                    });
                }
                if *t2 > 2.0 * *t1 * (1.0 + 1e-9) {
                    return Err(crate::error::PrismError::InvalidParameter {
                        message: "thermal relaxation requires t2 <= 2*t1".into(),
                    });
                }
                validate_probability("excited_population", *excited_population)?;
            }
            NoiseChannel::Custom { kraus } => validate_kraus_set("Custom", kraus)?,
            NoiseChannel::Kraus2q { kraus } => validate_kraus_set("Kraus2q", kraus)?,
            NoiseChannel::Leakage { p }
            | NoiseChannel::Seepage { p }
            | NoiseChannel::LeakageTransport { p } => validate_probability("p", *p)?,
            NoiseChannel::QuasiStatic { weights, .. } => {
                if weights.iter().any(|(_, weight)| !weight.is_finite()) {
                    return Err(crate::error::PrismError::InvalidParameter {
                        message: "quasi-static weights must be finite".into(),
                    });
                }
            }
        }
        Ok(())
    }
}

/// Deviation from `sum_k Kdagger_k K_k = I` tolerated in a declared Custom
/// channel. Loose enough for a set assembled from square roots, tight enough
/// that a channel scaled by a constant is rejected.
const CPTP_TOLERANCE: f64 = 1e-9;

/// Reject an empty operator set, a non-finite entry, and any deviation from
/// `sum_k Kdagger_k K_k = I` past [`CPTP_TOLERANCE`]. Shared by the one-qubit
/// and two-qubit dense channels, which differ only in `D`.
fn validate_kraus_set<const D: usize>(name: &str, kraus: &[[[Complex64; D]; D]]) -> Result<()> {
    if kraus.is_empty() {
        return Err(crate::error::PrismError::InvalidParameter {
            message: format!("{name} Kraus channel must contain at least one operator"),
        });
    }
    for (op_idx, op) in kraus.iter().enumerate() {
        for (row_idx, row) in op.iter().enumerate() {
            for (col_idx, val) in row.iter().enumerate() {
                if !val.re.is_finite() || !val.im.is_finite() {
                    return Err(crate::error::PrismError::InvalidParameter {
                        message: format!(
                            "{name} Kraus operator {op_idx} entry ({row_idx},{col_idx}) must be finite"
                        ),
                    });
                }
            }
        }
    }

    let mut sum = [[Complex64::new(0.0, 0.0); D]; D];
    for op in kraus {
        for (row, sum_row) in sum.iter_mut().enumerate() {
            for (col, entry) in sum_row.iter_mut().enumerate() {
                for op_row in op.iter() {
                    *entry += op_row[row].conj() * op_row[col];
                }
            }
        }
    }
    for (row, sum_row) in sum.iter().enumerate() {
        for (col, entry) in sum_row.iter().enumerate() {
            let expected = if row == col { 1.0 } else { 0.0 };
            if (entry.re - expected).abs() > CPTP_TOLERANCE || entry.im.abs() > CPTP_TOLERANCE {
                return Err(crate::error::PrismError::InvalidParameter {
                    message: format!(
                        "{name} Kraus operators must be trace preserving: \
                         sum of Kdagger K at ({row},{col}) is {entry}, expected {expected}"
                    ),
                });
            }
        }
    }
    Ok(())
}

fn validate_probability(name: &str, value: f64) -> Result<()> {
    if !value.is_finite() || !(0.0..=1.0).contains(&value) {
        return Err(crate::error::PrismError::InvalidParameter {
            message: format!("{name} must be finite and in [0, 1]"),
        });
    }
    Ok(())
}

/// A noise channel bound to its target qubits; fires after the instruction
/// its [`NoiseModel::after_gate`] slot is attached to.
#[derive(Debug, Clone, PartialEq)]
pub struct NoiseEvent {
    pub channel: NoiseChannel,
    /// As many entries as [`NoiseChannel::num_qubits`] requires, in the order
    /// the channel indexes them.
    pub qubits: SmallVec<[usize; 2]>,
}

impl NoiseEvent {
    pub fn pauli(qubit: usize, px: f64, py: f64, pz: f64) -> Self {
        Self {
            channel: NoiseChannel::Pauli { px, py, pz },
            qubits: smallvec![qubit],
        }
    }

    /// First target qubit (the only one for single-qubit channels).
    pub fn qubit(&self) -> usize {
        self.qubits[0]
    }

    /// Return the per-Pauli probabilities for a Pauli/Depolarizing channel.
    ///
    /// # Panics
    ///
    /// Panics if the channel is not Pauli or Depolarizing. Callers on the
    /// compiled stabilizer sampler path must guard with
    /// [`NoiseModel::is_pauli_only`] before invoking this method, and read a
    /// two-qubit channel through [`NoiseEvent::pauli_pair`] instead.
    pub fn pauli_probs(&self) -> (f64, f64, f64) {
        self.channel
            .as_pauli()
            .expect("pauli_probs called on non-Pauli channel; caller must check is_pauli_only")
    }

    /// Both targets and the total rate of a two-qubit Pauli channel, `None`
    /// when the channel names one qubit. An event answers this or
    /// [`NoiseEvent::pauli_probs`], never both.
    pub fn pauli_pair(&self) -> Option<(usize, usize, f64)> {
        self.channel
            .pauli_pair_rate()
            .map(|p| (self.qubits[0], self.qubits[1], p))
    }
}

/// Classical bit-flip error applied to a measurement outcome.
#[derive(Debug, Clone, PartialEq)]
pub struct ReadoutError {
    /// Probability a measured 0 reads out as 1.
    pub p01: f64,
    /// Probability a measured 1 reads out as 0.
    pub p10: f64,
}

impl ReadoutError {
    /// True when neither rate can flip a bit, so the entry changes no draw and
    /// a sweep that starts at zero is not a readout model yet.
    pub(crate) fn is_inert(&self) -> bool {
        self.p01 == 0.0 && self.p10 == 0.0
    }
}

/// Per-instruction noise attached to a specific circuit.
///
/// # Examples
///
/// ```
/// use prism_q::{Circuit, Gate, NoiseModel, run_shots_noisy};
///
/// let mut circuit = Circuit::new(2, 2);
/// circuit.add_gate(Gate::H, &[0]);
/// circuit.add_gate(Gate::Cx, &[0, 1]);
/// circuit.add_measure(0, 0);
/// circuit.add_measure(1, 1);
///
/// let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
/// let result = run_shots_noisy(&circuit, &noise, 100, 42)?;
/// assert_eq!(result.shots.len(), 100);
/// # Ok::<(), prism_q::PrismError>(())
/// ```
#[derive(Debug, Clone)]
pub struct NoiseModel {
    /// Events fired after each instruction, indexed by instruction position;
    /// length must equal the circuit's instruction count.
    pub after_gate: Vec<Vec<NoiseEvent>>,
    /// Per-classical-bit readout error; `None` means ideal readout. Either empty,
    /// for no readout error at all, or one entry per classical bit. A bit no
    /// measurement writes is never flipped.
    pub readout: Vec<Option<ReadoutError>>,
}

impl NoiseModel {
    /// Depolarizing noise with total probability `p` (not per branch) on
    /// every target of every gate, ideal readout.
    pub fn uniform_depolarizing(circuit: &Circuit, p: f64) -> Self {
        let px = p / 3.0;
        let py = p / 3.0;
        let pz = p / 3.0;

        let mut after_gate = Vec::with_capacity(circuit.instructions.len());
        for instr in &circuit.instructions {
            match instr {
                Instruction::Gate { targets, .. } => {
                    let ops: Vec<NoiseEvent> = targets
                        .iter()
                        .map(|&q| NoiseEvent::pauli(q, px, py, pz))
                        .collect();
                    after_gate.push(ops);
                }
                _ => {
                    after_gate.push(Vec::new());
                }
            }
        }

        Self {
            after_gate,
            readout: vec![None; circuit.num_classical_bits],
        }
    }

    /// Amplitude damping with rate `gamma` on every target of every gate,
    /// ideal readout.
    pub fn with_amplitude_damping(circuit: &Circuit, gamma: f64) -> Self {
        let mut after_gate = Vec::with_capacity(circuit.instructions.len());
        for instr in &circuit.instructions {
            match instr {
                Instruction::Gate { targets, .. } => {
                    let ops: Vec<NoiseEvent> = targets
                        .iter()
                        .map(|&q| NoiseEvent {
                            channel: NoiseChannel::AmplitudeDamping { gamma },
                            qubits: smallvec![q],
                        })
                        .collect();
                    after_gate.push(ops);
                }
                _ => {
                    after_gate.push(Vec::new());
                }
            }
        }

        Self {
            after_gate,
            readout: vec![None; circuit.num_classical_bits],
        }
    }

    /// Set the same readout error on every classical bit.
    pub fn with_readout_error(&mut self, p01: f64, p10: f64) -> &mut Self {
        self.readout.fill(Some(ReadoutError { p01, p10 }));
        self
    }

    /// Set the readout error of one classical bit, leaving the rest alone.
    ///
    /// Indexed by classical bit rather than by qubit: a qubit measured into
    /// several bits carries an independent rate per bit.
    ///
    /// # Panics
    ///
    /// Panics if `bit` is outside the classical register the model was built
    /// for, matching the index-bounds policy of the circuit builders.
    pub fn set_bit_readout_error(&mut self, bit: usize, p01: f64, p10: f64) -> &mut Self {
        assert!(
            bit < self.readout.len(),
            "classical bit {bit} is outside the {}-bit register",
            self.readout.len()
        );
        self.readout[bit] = Some(ReadoutError { p01, p10 });
        self
    }

    /// True when every channel is a single-qubit Pauli or Depolarizing channel
    /// and no readout entry can flip a bit. A live two-qubit channel answers
    /// `false` here and still runs on the Clifford samplers; an inert one
    /// blocks nothing.
    pub fn is_pauli_only(&self) -> bool {
        self.has_chain_complex_channels()
            && self.readout.iter().flatten().all(ReadoutError::is_inert)
    }

    /// True when every channel fits one column of the error chain complex: a
    /// single-qubit Pauli draw, or inert. Readout is not a channel and is
    /// judged separately, since only some of the routes over the complex need
    /// it folded into a syndrome class.
    pub(crate) fn has_chain_complex_channels(&self) -> bool {
        self.after_gate
            .iter()
            .flat_map(|events| events.iter())
            .all(|e| e.channel.is_pauli() || e.channel.is_inert())
    }

    /// True when every channel samples as a Pauli frame draw, whatever the
    /// readout entries hold. Readout acts on the measurement record after
    /// sampling, so the Clifford samplers apply it themselves and only the
    /// channels decide whether a model can take that route.
    pub(crate) fn has_only_pauli_channels(&self) -> bool {
        self.after_gate
            .iter()
            .flat_map(|events| events.iter())
            .all(|e| e.channel.is_pauli_frame())
    }

    /// True when any event carries a [`NoiseChannel::Kraus2q`], the one channel
    /// needing a two-qubit reduction from the backend running the trajectory.
    pub(crate) fn has_two_qubit_kraus(&self) -> bool {
        self.after_gate
            .iter()
            .flat_map(|events| events.iter())
            .any(|event| matches!(event.channel, NoiseChannel::Kraus2q { .. }))
    }

    /// True when any event carries a leakage channel, which gives each
    /// trajectory a leak flag per qubit and each shot a leaked-qubit record.
    pub(crate) fn has_leakage(&self) -> bool {
        self.after_gate.iter().flatten().any(|event| {
            matches!(
                event.channel,
                NoiseChannel::Leakage { .. }
                    | NoiseChannel::Seepage { .. }
                    | NoiseChannel::LeakageTransport { .. }
            )
        })
    }

    /// Standard normal draws each shot makes for its
    /// [`NoiseChannel::QuasiStatic`] events: one past the highest source index.
    pub(crate) fn quasi_static_sources(&self) -> usize {
        self.after_gate
            .iter()
            .flatten()
            .filter_map(|event| match &event.channel {
                NoiseChannel::QuasiStatic { weights, .. } => {
                    weights.iter().map(|&(source, _)| source + 1).max()
                }
                _ => None,
            })
            .max()
            .unwrap_or(0)
    }

    /// Return an error naming `engine` if any channel needs per-shot memory:
    /// a leak flag or a quasi-static draw. A mixture averages each event on its
    /// own, which neither survives.
    pub(crate) fn ensure_memoryless(&self, engine: &str) -> Result<()> {
        if self.has_leakage() || self.quasi_static_sources() > 0 {
            return Err(crate::error::PrismError::IncompatibleBackend {
                backend: engine.into(),
                reason: "leakage and quasi-static channels carry state from one event to the \
                         next within a shot, which a mixture evolved event by event does not \
                         hold; sample them through per-shot trajectories"
                    .into(),
            });
        }
        Ok(())
    }

    pub fn has_noise(&self) -> bool {
        self.after_gate.iter().any(|events| !events.is_empty())
            || self.readout.iter().any(|r| r.is_some())
    }

    /// Return an error if any channel falls outside the error chain complex,
    /// which holds one column per single-qubit Pauli error.
    pub(crate) fn ensure_chain_complex_channels(&self) -> Result<()> {
        if !self.has_chain_complex_channels() {
            return Err(crate::error::PrismError::IncompatibleBackend {
                backend: "error chain complex".into(),
                reason: "non-Pauli noise channels (amplitude damping, phase damping, \
                         thermal relaxation, custom Kraus) carry no Pauli frame, and a \
                         two-qubit channel spans two error columns where the complex \
                         holds one; sample through `run_shots_noisy`"
                    .into(),
            });
        }
        Ok(())
    }

    /// Return an error if any channel is non-Pauli, the precondition of the
    /// Clifford samplers. Readout error is allowed: they apply it to the
    /// measurement record after the reference bits are folded in.
    pub(crate) fn ensure_only_pauli_channels(&self) -> Result<()> {
        if !self.has_only_pauli_channels() {
            return Err(crate::error::PrismError::IncompatibleBackend {
                backend: "compiled stabilizer sampler".into(),
                reason: "non-Pauli noise channels (amplitude damping, phase damping, \
                         thermal relaxation, custom Kraus) require the trajectory \
                         engine; use a non-stabilizer backend or a Pauli-only noise \
                         model"
                    .into(),
            });
        }
        Ok(())
    }

    /// Validate as [`NoiseModel::validate`] does, plus the checks that need the
    /// circuit: one event slot per instruction, a readout table that is empty
    /// or holds one entry per classical bit, and every target qubit inside the
    /// register.
    ///
    /// Every noisy entry point calls this before allocating state. Qubit bounds
    /// cannot be checked without the circuit, and an out-of-range target reaches
    /// kernels that index amplitudes without one.
    pub fn validate_for(&self, circuit: &Circuit) -> Result<()> {
        if self.after_gate.len() != circuit.instructions.len() {
            return Err(crate::error::PrismError::InvalidParameter {
                message: format!(
                    "noise model carries {} event slots for {} instructions",
                    self.after_gate.len(),
                    circuit.instructions.len()
                ),
            });
        }
        if !self.readout.is_empty() && self.readout.len() != circuit.num_classical_bits {
            return Err(crate::error::PrismError::InvalidParameter {
                message: format!(
                    "noise model carries {} readout entries for {} classical bits",
                    self.readout.len(),
                    circuit.num_classical_bits
                ),
            });
        }
        self.validate()?;

        // One event slot per top-level instruction, so a region's body has no
        // slots to carry noise and would run noiselessly without this. A model
        // carrying no quantum events has nothing to lose there, so readout-only
        // noise stays legal on a circuit with regions.
        if self.after_gate.iter().any(|events| !events.is_empty())
            && let Some(index) = circuit
                .instructions
                .iter()
                .position(|inst| matches!(inst, Instruction::Region(_)))
        {
            return Err(crate::error::PrismError::IncompatibleBackend {
                backend: "NoiseModel".to_string(),
                reason: format!(
                    "noise events are indexed per instruction, which leaves the body of the \
                         guarded region at instruction {index} without noise slots"
                ),
            });
        }

        for (instr_idx, events) in self.after_gate.iter().enumerate() {
            for (event_idx, event) in events.iter().enumerate() {
                for &qubit in &event.qubits {
                    if qubit >= circuit.num_qubits {
                        return Err(crate::error::PrismError::InvalidParameter {
                            message: format!(
                                "noise event {event_idx} after instruction {instr_idx}: qubit {qubit} is outside the {}-qubit register",
                                circuit.num_qubits
                            ),
                        });
                    }
                }
            }
        }
        Ok(())
    }

    /// Validate all channels and readout errors in this noise model.
    pub fn validate(&self) -> Result<()> {
        for (instr_idx, events) in self.after_gate.iter().enumerate() {
            for (event_idx, event) in events.iter().enumerate() {
                event.channel.validate().map_err(|err| {
                    crate::error::PrismError::InvalidParameter {
                        message: format!(
                            "noise event {event_idx} after instruction {instr_idx}: {err}"
                        ),
                    }
                })?;

                let expected_qubits = event.channel.num_qubits();
                if event.qubits.len() != expected_qubits {
                    return Err(crate::error::PrismError::InvalidParameter {
                        message: format!(
                            "noise event {event_idx} after instruction {instr_idx}: expected {expected_qubits} qubit(s), got {}",
                            event.qubits.len()
                        ),
                    });
                }
                if expected_qubits == 2 && event.qubits[0] == event.qubits[1] {
                    return Err(crate::error::PrismError::InvalidParameter {
                        message: format!(
                            "noise event {event_idx} after instruction {instr_idx}: two-qubit channel needs distinct targets, got qubit {} twice",
                            event.qubits[0]
                        ),
                    });
                }
            }
        }

        for (bit_idx, readout) in self.readout.iter().enumerate() {
            if let Some(readout) = readout {
                validate_probability("readout p01", readout.p01).map_err(|err| {
                    crate::error::PrismError::InvalidParameter {
                        message: format!("readout error {bit_idx}: {err}"),
                    }
                })?;
                validate_probability("readout p10", readout.p10).map_err(|err| {
                    crate::error::PrismError::InvalidParameter {
                        message: format!("readout error {bit_idx}: {err}"),
                    }
                })?;
            }
        }

        Ok(())
    }
}

/// Rows a two-qubit event stores: the propagated X and Z components of each
/// of its two targets, in the order `x0, z0, x1, z1`. The 15 branches are
/// drawn jointly at sample time and read these four rows, rather than being
/// precomputed as 15 rows of their own.
const PAIR_ROWS: usize = 4;

struct FlatNoiseSensitivity {
    x_data: Vec<u64>,
    z_data: Vec<u64>,
    probs: Vec<[f64; 3]>,
    /// [`PAIR_ROWS`] rows of `m_words` per two-qubit event. Kept apart from
    /// the single-qubit rows, which the LUT build, the device CSR transpose
    /// and the threshold table all index positionally.
    pair_data: Vec<u64>,
    pair_probs: Vec<f64>,
    m_words: usize,
}

impl FlatNoiseSensitivity {
    fn new(m_words: usize, capacity: usize) -> Self {
        Self {
            x_data: Vec::with_capacity(capacity * m_words),
            z_data: Vec::with_capacity(capacity * m_words),
            probs: Vec::with_capacity(capacity),
            pair_data: Vec::new(),
            pair_probs: Vec::new(),
            m_words,
        }
    }

    fn push(&mut self, x_flip: &[u64], z_flip: &[u64], px: f64, py: f64, pz: f64) {
        self.x_data.extend_from_slice(x_flip);
        self.z_data.extend_from_slice(z_flip);
        self.probs.push([px, py, pz]);
    }

    fn push_pair(&mut self, x0: &[u64], z0: &[u64], x1: &[u64], z1: &[u64], p: f64) {
        for row in [x0, z0, x1, z1] {
            self.pair_data.extend_from_slice(row);
        }
        self.pair_probs.push(p);
    }

    /// The device noise kernel derives each thread's bit from
    /// `(seed, event, batch)` for the one measurement row it owns, so a branch
    /// drawn there would decorrelate the two qubits. A non-empty pair table
    /// therefore keeps noise application on the host.
    #[inline(always)]
    fn has_pairs(&self) -> bool {
        !self.pair_probs.is_empty()
    }

    /// One row of pair event `idx`, indexed as [`PAIR_ROWS`] describes.
    #[inline(always)]
    fn pair_row(&self, idx: usize, row: usize) -> &[u64] {
        let off = (idx * PAIR_ROWS + row) * self.m_words;
        &self.pair_data[off..off + self.m_words]
    }

    #[inline(always)]
    fn len(&self) -> usize {
        self.probs.len()
    }

    #[inline(always)]
    fn is_empty(&self) -> bool {
        self.probs.is_empty()
    }

    #[inline(always)]
    fn x_flip(&self, idx: usize) -> &[u64] {
        let off = idx * self.m_words;
        &self.x_data[off..off + self.m_words]
    }

    #[inline(always)]
    fn z_flip(&self, idx: usize) -> &[u64] {
        let off = idx * self.m_words;
        &self.z_data[off..off + self.m_words]
    }
}

/// XOR the measurements one Pauli letter flips into `shot`. An X error
/// anticommutes with the propagated Z component, a Z error with the X
/// component, and a Y error with both.
#[inline(always)]
fn xor_pair_letter(shot: &mut [u64], x_row: &[u64], z_row: &[u64], letter: usize) {
    match letter {
        1 => xor_words(shot, z_row),
        2 => {
            xor_words(shot, x_row);
            xor_words(shot, z_row);
        }
        3 => xor_words(shot, x_row),
        _ => {}
    }
}

/// `u` is uniform on `[0, 1)` conditioned on the event having fired.
#[inline(always)]
fn xor_pair_branch(events: &FlatNoiseSensitivity, idx: usize, u: f64, shot: &mut [u64]) {
    let (l0, l1) = pair_branch(u);
    xor_pair_letter(shot, events.pair_row(idx, 0), events.pair_row(idx, 1), l0);
    xor_pair_letter(shot, events.pair_row(idx, 2), events.pair_row(idx, 3), l1);
}

/// The Pauli letter on each target of one uniformly drawn two-qubit branch,
/// as 0 for identity through 3 for Z, from a uniform draw on `[0, 1)`.
///
/// The identity product is not a branch: a firing shot picks one of the 15
/// others, so the index runs 1 through 15 and every firing shot moves.
#[inline(always)]
fn pair_branch(u: f64) -> (usize, usize) {
    let idx = 1 + ((u * 15.0) as usize).min(14);
    (idx >> 2, idx & 3)
}

#[inline(always)]
fn geometric_sample(rng: &mut ChaCha8Rng, ln_1mp: f64) -> usize {
    let u: f64 = 1.0 - rand::RngExt::random::<f64>(rng);
    (u.ln() / ln_1mp) as usize
}

/// [`geometric_sample`] on the sampler's own generator. `ln_1mp` is `(1 - p).ln()`,
/// precomputed by the caller.
#[inline(always)]
pub(crate) fn geometric_sample_xoshiro(
    rng: &mut crate::sim::compiled::rng::Xoshiro256PlusPlus,
    ln_1mp: f64,
) -> usize {
    // `next_f64` returns [0, 1); shift to (0, 1] so `ln` stays finite.
    let u: f64 = 1.0 - rng.next_f64();
    (u.ln() / ln_1mp) as usize
}

/// RNG source for [`fill_one_event`]. Bundles uniform and geometric draws
/// onto a single trait so the shared body can use one `&mut` borrow instead
/// of two concurrent closures over the underlying generator.
#[cfg(feature = "gpu")]
trait EventRng {
    fn uniform(&mut self) -> f64;
    fn geom(&mut self, ln_1mp: f64) -> usize;
}

#[cfg(feature = "gpu")]
impl EventRng for ChaCha8Rng {
    #[inline]
    fn uniform(&mut self) -> f64 {
        rand::RngExt::random::<f64>(self)
    }
    #[inline]
    fn geom(&mut self, ln_1mp: f64) -> usize {
        geometric_sample(self, ln_1mp)
    }
}

#[cfg(feature = "gpu")]
impl EventRng for crate::sim::compiled::rng::Xoshiro256PlusPlus {
    #[inline]
    fn uniform(&mut self) -> f64 {
        self.next_f64()
    }
    #[inline]
    fn geom(&mut self, ln_1mp: f64) -> usize {
        geometric_sample_xoshiro(self, ln_1mp)
    }
}

/// Fill one event's portion of the noise masks. Shared by the serial master
/// ChaCha path and the parallel per-event Xoshiro path, so the bernoulli
/// and geometric branches live in exactly one place.
#[cfg(feature = "gpu")]
#[inline]
fn fill_one_event(
    x_slot: &mut [u64],
    z_slot: &mut [u64],
    probs: [f64; 3],
    chunk_shots: usize,
    rng: &mut dyn EventRng,
) {
    let [px, py, pz] = probs;
    let p_event = px + py + pz;
    if p_event == 0.0 {
        return;
    }
    let pxy = px + py;

    // High-p or tiny chunk: a per-shot bernoulli sweep dominates geometric
    // skipping overhead because most shots flip anyway.
    if p_event >= 0.5 || chunk_shots < 32 {
        for s in 0..chunk_shots {
            let r = rng.uniform();
            let slot_idx = s / 64;
            let bit = 1u64 << (s % 64);
            if r < px {
                z_slot[slot_idx] |= bit;
            } else if r < pxy {
                z_slot[slot_idx] |= bit;
                x_slot[slot_idx] |= bit;
            } else if r < p_event {
                x_slot[slot_idx] |= bit;
            }
        }
        return;
    }

    // Low-p: skip inactive shots via geometric sampling. `px_frac` / `pxy_frac`
    // rescale to the conditional distribution given an event fired.
    let ln_1mp = (1.0 - p_event).ln();
    let px_frac = px / p_event;
    let pxy_frac = pxy / p_event;
    let mut pos = rng.geom(ln_1mp);
    while pos < chunk_shots {
        let r = rng.uniform();
        let slot_idx = pos / 64;
        let bit = 1u64 << (pos % 64);
        if r < px_frac {
            z_slot[slot_idx] |= bit;
        } else if r < pxy_frac {
            z_slot[slot_idx] |= bit;
            x_slot[slot_idx] |= bit;
        } else {
            x_slot[slot_idx] |= bit;
        }
        pos += 1 + rng.geom(ln_1mp);
    }
}

const NOISE_LUT_K: usize = 8;
const NOISE_LUT_MIN_EVENTS: usize = 16;
const NOISE_LUT_TILE: usize = 4096;

fn unpack_and_remap(
    accum: &[u64],
    m_words: usize,
    num_shots: usize,
    classical_bit_order: &[usize],
    num_classical: usize,
) -> Vec<Vec<bool>> {
    fn unpack_one(src: &[u64], classical_bit_order: &[usize], num_classical: usize) -> Vec<bool> {
        let mut out = vec![false; num_classical];
        for (mi, &cbit) in classical_bit_order.iter().enumerate() {
            if cbit < num_classical {
                out[cbit] = (src[mi / 64] >> (mi % 64)) & 1 != 0;
            }
        }
        out
    }

    #[cfg(feature = "parallel")]
    if num_shots >= 256 {
        use rayon::prelude::*;
        return accum
            .par_chunks(m_words)
            .map(|src| unpack_one(src, classical_bit_order, num_classical))
            .collect();
    }

    #[cfg(not(feature = "parallel"))]
    let _ = num_shots;

    accum
        .chunks(m_words)
        .map(|src| unpack_one(src, classical_bit_order, num_classical))
        .collect()
}

fn unpack_and_remap_packed(
    packed: &PackedShots,
    num_shots: usize,
    classical_bit_order: &[usize],
    num_classical: usize,
) -> Vec<Vec<bool>> {
    fn unpack_one_packed(
        packed: &PackedShots,
        shot: usize,
        classical_bit_order: &[usize],
        num_classical: usize,
    ) -> Vec<bool> {
        let mut out = vec![false; num_classical];
        for (mi, &cbit) in classical_bit_order.iter().enumerate() {
            if cbit < num_classical {
                out[cbit] = packed.get_bit(shot, mi);
            }
        }
        out
    }

    #[cfg(feature = "parallel")]
    if num_shots >= 256 {
        use rayon::prelude::*;
        return (0..num_shots)
            .into_par_iter()
            .map(|s| unpack_one_packed(packed, s, classical_bit_order, num_classical))
            .collect();
    }

    #[cfg(not(feature = "parallel"))]
    let _ = num_shots;

    (0..packed.num_shots())
        .map(|s| unpack_one_packed(packed, s, classical_bit_order, num_classical))
        .collect()
}

/// Readout error on one measurement record, the index the packed shot buffers
/// carry. Only records carrying an entry that can flip a bit appear.
struct RecordReadout {
    record: usize,
    p01: f64,
    p10: f64,
}

fn record_readout(circuit: &Circuit, noise: &NoiseModel) -> Vec<RecordReadout> {
    circuit
        .classical_bit_order()
        .into_iter()
        .enumerate()
        .filter_map(|(record, cbit)| {
            let err = noise.readout.get(cbit)?.as_ref()?;
            (!err.is_inert()).then_some(RecordReadout {
                record,
                p01: err.p01,
                p10: err.p10,
            })
        })
        .collect()
}

/// Flip readout errors into a shot-major packed record, in place.
///
/// Runs after the reference bits are folded in, so a set bit is a measured one
/// and takes `p10`. Asymmetric rates rule out a single flip mask: candidates
/// are thinned at `max(p01, p10)`, then accepted at the rate the live bit
/// selects. A record with no readout entry draws nothing.
fn apply_readout_packed(
    accum: &mut [u64],
    num_shots: usize,
    m_words: usize,
    readout: &[RecordReadout],
    rng: &mut ChaCha8Rng,
) {
    for entry in readout {
        let p_max = entry.p01.max(entry.p10);
        let word = entry.record / 64;
        let bit = 1u64 << (entry.record % 64);

        if p_max >= 0.5 || num_shots < 32 {
            for s in 0..num_shots {
                let slot = &mut accum[s * m_words + word];
                let p = if *slot & bit != 0 {
                    entry.p10
                } else {
                    entry.p01
                };
                let r: f64 = rand::RngExt::random(rng);
                if r < p {
                    *slot ^= bit;
                }
            }
            continue;
        }

        let ln_1mp = (1.0 - p_max).ln();
        let accept_set = entry.p10 / p_max;
        let accept_clear = entry.p01 / p_max;
        let mut pos = geometric_sample(rng, ln_1mp);
        while pos < num_shots {
            let slot = &mut accum[pos * m_words + word];
            let accept = if *slot & bit != 0 {
                accept_set
            } else {
                accept_clear
            };
            let r: f64 = rand::RngExt::random(rng);
            if r < accept {
                *slot ^= bit;
            }
            pos += 1 + geometric_sample(rng, ln_1mp);
        }
    }
}

struct NoiseFlipLut {
    data: Vec<u64>,
    m_words: usize,
    num_full_groups: usize,
    remainder_size: usize,
}

impl NoiseFlipLut {
    fn build_from_flat(flat_data: &[u64], num_rows: usize, m_words: usize) -> Self {
        let num_full_groups = num_rows / NOISE_LUT_K;
        let remainder_size = num_rows % NOISE_LUT_K;
        let total_groups = num_full_groups + usize::from(remainder_size > 0);
        let entries = 1 << NOISE_LUT_K;

        let mut data = vec![0u64; total_groups * entries * m_words];

        for g in 0..total_groups {
            let group_start = g * NOISE_LUT_K;
            let k = if g < num_full_groups {
                NOISE_LUT_K
            } else {
                remainder_size
            };
            let lut_off = g * entries * m_words;

            for byte in 1..(1usize << k) {
                let lowest = byte & byte.wrapping_neg();
                let row_idx = group_start + lowest.trailing_zeros() as usize;
                let prev = byte ^ lowest;

                let dst = lut_off + byte * m_words;
                let src = lut_off + prev * m_words;
                let row_off = row_idx * m_words;

                for w in 0..m_words {
                    data[dst + w] = data[src + w] ^ flat_data[row_off + w];
                }
            }
        }

        Self {
            data,
            m_words,
            num_full_groups,
            remainder_size,
        }
    }

    #[inline(always)]
    fn total_groups(&self) -> usize {
        self.num_full_groups + usize::from(self.remainder_size > 0)
    }

    #[inline(always)]
    fn apply_masked(&self, accum: &mut [u64], mask: &[u64]) {
        for g in 0..self.total_groups() {
            let byte = ((mask[g / 8] >> ((g % 8) * 8)) & 0xFF) as usize;
            if byte != 0 {
                let offset = (g * (1 << NOISE_LUT_K) + byte) * self.m_words;
                xor_words(accum, &self.data[offset..offset + self.m_words]);
            }
        }
    }
}

fn build_noise_luts(events: &FlatNoiseSensitivity) -> (Option<NoiseFlipLut>, Option<NoiseFlipLut>) {
    let ne = events.len();
    if ne < NOISE_LUT_MIN_EVENTS {
        return (None, None);
    }

    let avg_p: f64 = events
        .probs
        .iter()
        .map(|[px, py, pz]| px + py + pz)
        .sum::<f64>()
        / ne as f64;
    if avg_p < 0.05 {
        return (None, None);
    }

    let mw = events.m_words;
    let z_lut = NoiseFlipLut::build_from_flat(&events.z_data, ne, mw);
    let x_lut = NoiseFlipLut::build_from_flat(&events.x_data, ne, mw);
    (Some(z_lut), Some(x_lut))
}

/// Pack `(px, py, pz)` per event into three `u64` thresholds for the fused
/// noise kernel. Each threshold is `p * 2^64` rounded down, so a uniform
/// `u64` random compares below it with probability `p`. Returns the flat
/// buffer along with `false` when any event's probabilities are out of range
/// (negative, NaN, or sum greater than 1). The caller then falls back to the
/// host mask path.
#[cfg(feature = "gpu")]
fn flatten_event_thresholds_u64(probs: &[[f64; 3]]) -> (Vec<u64>, bool) {
    // Nearest f64 below 2^64. Scaling by 2^64 exactly rounds to 2^64 which
    // doesn't fit in a u64; this preserves `u64_rand < threshold` at the
    // intended probability.
    const SCALE: f64 = 18446744073709549568.0_f64;
    const PROB_SUM_SLACK: f64 = 1e-12;

    let mut out = Vec::with_capacity(probs.len() * 3);
    let mut valid = true;
    for &[px, py, pz] in probs {
        let finite = px.is_finite() && py.is_finite() && pz.is_finite();
        let sum = px + py + pz;
        if !finite || px < 0.0 || py < 0.0 || pz < 0.0 || sum > 1.0 + PROB_SUM_SLACK {
            valid = false;
        }
        let to_u64 = |p: f64| -> u64 {
            if !p.is_finite() || p <= 0.0 {
                0
            } else if p >= 1.0 {
                u64::MAX
            } else {
                (p * SCALE) as u64
            }
        };
        out.push(to_u64(px));
        out.push(to_u64(px + py));
        out.push(to_u64(sum));
    }
    (out, valid)
}

/// Transpose the event-major `(x_data, z_data)` bitmaps in a
/// `FlatNoiseSensitivity` into a row-major CSR. Each row's entry list names
/// the events that touch it, packed as `event << 2 | flag`: flag bit 0 for
/// an X contribution, bit 1 for Z, both set for a Y contribution. Lets the
/// per-(row, batch) kernel accumulate XORs in a register without atomics.
/// Paired with `bts_generate_and_apply_noise_meas_major_by_row`.
#[cfg(feature = "gpu")]
fn transpose_events_to_row_major(
    events: &FlatNoiseSensitivity,
    num_meas: usize,
) -> (Vec<u32>, Vec<u32>) {
    let num_events = events.len();
    let m_words = events.m_words;

    // First pass: count per-row entries.
    let mut row_counts = vec![0u32; num_meas];
    for e in 0..num_events {
        let x = &events.x_data[e * m_words..(e + 1) * m_words];
        let z = &events.z_data[e * m_words..(e + 1) * m_words];
        for w in 0..m_words {
            let union = x[w] | z[w];
            let mut bits = union;
            while bits != 0 {
                let row = w * 64 + bits.trailing_zeros() as usize;
                if row < num_meas {
                    row_counts[row] += 1;
                }
                bits &= bits - 1;
            }
        }
    }

    let mut offsets = Vec::with_capacity(num_meas + 1);
    offsets.push(0u32);
    let mut running = 0u32;
    for c in &row_counts {
        running += c;
        offsets.push(running);
    }
    let total = running as usize;

    // Second pass: write entries through a cursor starting at `offsets[row]`.
    let mut entries = vec![0u32; total];
    let mut cursor: Vec<u32> = offsets[..num_meas].to_vec();
    for e in 0..num_events {
        let x = &events.x_data[e * m_words..(e + 1) * m_words];
        let z = &events.z_data[e * m_words..(e + 1) * m_words];
        for w in 0..m_words {
            let union = x[w] | z[w];
            let mut bits = union;
            while bits != 0 {
                let b = bits.trailing_zeros() as usize;
                bits &= bits - 1;
                let row = w * 64 + b;
                if row >= num_meas {
                    continue;
                }
                let bit_mask = 1u64 << b;
                let flag: u32 = (if x[w] & bit_mask != 0 { 1 } else { 0 })
                    | (if z[w] & bit_mask != 0 { 2 } else { 0 });
                let slot = cursor[row] as usize;
                entries[slot] = ((e as u32) << 2) | flag;
                cursor[row] += 1;
            }
        }
    }

    (offsets, entries)
}

#[cfg(feature = "gpu")]
fn flatten_event_rows_to_csr(
    flat_data: &[u64],
    num_rows: usize,
    m_words: usize,
) -> (Vec<u32>, Vec<u32>) {
    let mut row_offsets = Vec::with_capacity(num_rows + 1);
    let mut row_indices = Vec::new();

    for row in 0..num_rows {
        row_offsets.push(row_indices.len() as u32);
        let base = row * m_words;
        for word in 0..m_words {
            let mut bits = flat_data[base + word];
            while bits != 0 {
                let bit = bits.trailing_zeros() as usize;
                row_indices.push((word * 64 + bit) as u32);
                bits &= bits - 1;
            }
        }
    }
    row_offsets.push(row_indices.len() as u32);

    (row_offsets, row_indices)
}

#[cfg(feature = "gpu")]
struct GpuNoiseCache {
    context: std::sync::Arc<crate::gpu::GpuContext>,
    x_row_offsets_dev: crate::gpu::GpuBuffer<u32>,
    x_row_indices_dev: crate::gpu::GpuBuffer<u32>,
    z_row_offsets_dev: crate::gpu::GpuBuffer<u32>,
    z_row_indices_dev: crate::gpu::GpuBuffer<u32>,
    x_masks_dev: crate::gpu::GpuBuffer<u64>,
    z_masks_dev: crate::gpu::GpuBuffer<u64>,
    x_masks_host: Vec<u64>,
    z_masks_host: Vec<u64>,
    /// Precomputed per-event thresholds for the fused generator. Three u64s
    /// per event hold `[px, px+py, px+py+pz]` scaled to `2^64`, so the
    /// kernel's inner loop replaces a floating-point multiply and branch per
    /// shot with a `u64 <` compare. Stays resident across sampling calls.
    event_thresholds_dev: crate::gpu::GpuBuffer<u64>,
    /// True when every event's probabilities can be represented as a `u64`
    /// threshold. A NaN or a sum greater than 1 flips this to false and
    /// routes the caller to the host mask path.
    event_thresholds_valid: bool,
    /// Row-major event CSR for the per-(row, batch) fused kernel. Each entry
    /// packs `event << 2 | flag`, where bit 0 means the event contributes X
    /// to this row and bit 1 means Z. The row-major layout lets one thread
    /// accumulate every XOR for its output word in a register instead of
    /// broadcasting `(event, batch)` writes across rows with atomics.
    row_event_offsets_dev: crate::gpu::GpuBuffer<u32>,
    row_event_entries_dev: crate::gpu::GpuBuffer<u32>,
}

#[cfg(feature = "gpu")]
impl GpuNoiseCache {
    fn new(
        context: std::sync::Arc<crate::gpu::GpuContext>,
        events: &FlatNoiseSensitivity,
    ) -> Result<Self> {
        let device = context.device();
        let num_events = events.len();
        let (x_row_offsets, x_row_indices) =
            flatten_event_rows_to_csr(&events.x_data, num_events, events.m_words);
        let (z_row_offsets, z_row_indices) =
            flatten_event_rows_to_csr(&events.z_data, num_events, events.m_words);

        let mut x_row_offsets_dev =
            crate::gpu::GpuBuffer::<u32>::alloc_zeros(device, x_row_offsets.len().max(1))?;
        let mut x_row_indices_dev =
            crate::gpu::GpuBuffer::<u32>::alloc_zeros(device, x_row_indices.len().max(1))?;
        let mut z_row_offsets_dev =
            crate::gpu::GpuBuffer::<u32>::alloc_zeros(device, z_row_offsets.len().max(1))?;
        let mut z_row_indices_dev =
            crate::gpu::GpuBuffer::<u32>::alloc_zeros(device, z_row_indices.len().max(1))?;

        if !x_row_offsets.is_empty() {
            x_row_offsets_dev.copy_from_host(device, &x_row_offsets)?;
        }
        if !x_row_indices.is_empty() {
            x_row_indices_dev.copy_from_host(device, &x_row_indices)?;
        }
        if !z_row_offsets.is_empty() {
            z_row_offsets_dev.copy_from_host(device, &z_row_offsets)?;
        }
        if !z_row_indices.is_empty() {
            z_row_indices_dev.copy_from_host(device, &z_row_indices)?;
        }

        let (thresholds_host, thresholds_valid) = flatten_event_thresholds_u64(&events.probs);
        let mut event_thresholds_dev =
            crate::gpu::GpuBuffer::<u64>::alloc_zeros(device, thresholds_host.len().max(1))?;
        if !thresholds_host.is_empty() {
            event_thresholds_dev.copy_from_host(device, &thresholds_host)?;
        }

        // Row-major transpose for the per-(row, batch) fused kernel. Cheap
        // bit-scan over events.x_data | events.z_data at cache init time.
        let (row_event_offsets_host, row_event_entries_host) =
            transpose_events_to_row_major(events, events.m_words * 64);
        let mut row_event_offsets_dev =
            crate::gpu::GpuBuffer::<u32>::alloc_zeros(device, row_event_offsets_host.len().max(1))?;
        let mut row_event_entries_dev =
            crate::gpu::GpuBuffer::<u32>::alloc_zeros(device, row_event_entries_host.len().max(1))?;
        if !row_event_offsets_host.is_empty() {
            row_event_offsets_dev.copy_from_host(device, &row_event_offsets_host)?;
        }
        if !row_event_entries_host.is_empty() {
            row_event_entries_dev.copy_from_host(device, &row_event_entries_host)?;
        }

        Ok(Self {
            context: context.clone(),
            x_row_offsets_dev,
            x_row_indices_dev,
            z_row_offsets_dev,
            z_row_indices_dev,
            x_masks_dev: crate::gpu::GpuBuffer::<u64>::alloc_zeros(device, 1)?,
            z_masks_dev: crate::gpu::GpuBuffer::<u64>::alloc_zeros(device, 1)?,
            x_masks_host: Vec::new(),
            z_masks_host: Vec::new(),
            event_thresholds_dev,
            event_thresholds_valid: thresholds_valid,
            row_event_offsets_dev,
            row_event_entries_dev,
        })
    }

    fn ensure_mask_capacity(&mut self, len: usize) -> Result<()> {
        let needed = len.max(1);
        let device = self.context.device();
        if self.x_masks_dev.len() < needed {
            self.x_masks_dev = crate::gpu::GpuBuffer::<u64>::alloc_zeros(device, needed)?;
        }
        if self.z_masks_dev.len() < needed {
            self.z_masks_dev = crate::gpu::GpuBuffer::<u64>::alloc_zeros(device, needed)?;
        }
        if self.x_masks_host.len() < needed {
            self.x_masks_host.resize(needed, 0);
        }
        if self.z_masks_host.len() < needed {
            self.z_masks_host.resize(needed, 0);
        }
        Ok(())
    }
}

/// Compiled sampler for Clifford circuits with Pauli noise.
///
/// Reuses the noiseless circuit's compiled parity map and applies the Pauli
/// noise per batch of shots. With a GPU context attached (`with_gpu`, `gpu`
/// feature), parity sampling may take the GPU BTS path. Materialized shots return
/// to the host in shot-major form; counts and marginals can keep the packed
/// measurement-major buffer on the device through noise and reduction. `try_`
/// methods propagate GPU errors; their plain twins fall back to the CPU path.
pub struct NoisyCompiledSampler {
    noiseless: crate::sim::compiled::CompiledSampler,
    events: FlatNoiseSensitivity,
    /// Readout error per measurement record, applied on the host once the
    /// reference bits are folded in. A non-empty table keeps sampling off the
    /// device reductions, which read the record before readout exists.
    readout: Vec<RecordReadout>,
    num_measurements: usize,
    rng: ChaCha8Rng,
    z_lut: Option<NoiseFlipLut>,
    x_lut: Option<NoiseFlipLut>,
    #[cfg(feature = "gpu")]
    gpu_noise_cache: Option<GpuNoiseCache>,
}

impl NoisyCompiledSampler {
    /// Opt in to GPU-accelerated sampling for the underlying compiled
    /// noiseless sampler.
    ///
    /// Small shot batches may continue to use the CPU path when the compiled
    /// sampler stays below the GPU BTS threshold. On the device path, noise
    /// masks are generated by a device-seeded RNG stream when the compiled
    /// event probabilities permit, otherwise on the host and uploaded; counts
    /// and marginals then apply them and reduce on the device. Device draws
    /// use a different RNG stream than the CPU path, so GPU output matches
    /// CPU output statistically, not shot for shot.
    #[cfg(feature = "gpu")]
    pub fn with_gpu(mut self, context: std::sync::Arc<crate::gpu::GpuContext>) -> Self {
        self.noiseless = self.noiseless.with_gpu(context);
        self.gpu_noise_cache = None;
        self
    }

    fn sample_bulk_words_shot_major_cpu(&mut self, num_shots: usize) -> (Vec<u64>, usize) {
        let (mut accum, m_words) = self.noiseless.sample_bulk_words_shot_major(num_shots);
        self.apply_noise_bulk(&mut accum, num_shots, m_words);

        let ref_bits_packed = self.noiseless.ref_bits_packed();

        #[cfg(feature = "parallel")]
        if num_shots >= 256 {
            use rayon::prelude::*;
            accum
                .par_chunks_mut(m_words)
                .for_each(|shot| xor_words(shot, ref_bits_packed));
        } else {
            for s in 0..num_shots {
                let shot_base = s * m_words;
                xor_words(&mut accum[shot_base..shot_base + m_words], ref_bits_packed);
            }
        }

        #[cfg(not(feature = "parallel"))]
        for s in 0..num_shots {
            let shot_base = s * m_words;
            xor_words(&mut accum[shot_base..shot_base + m_words], ref_bits_packed);
        }

        apply_readout_packed(&mut accum, num_shots, m_words, &self.readout, &mut self.rng);

        (accum, m_words)
    }

    #[cfg(feature = "gpu")]
    fn try_sample_bulk_packed_with_gpu_context(&mut self, num_shots: usize) -> Result<PackedShots> {
        let m_words = self.num_measurements.div_ceil(64);
        if num_shots == 0 || self.num_measurements == 0 {
            return Ok(PackedShots::from_shot_major(
                vec![0u64; num_shots * m_words],
                num_shots,
                self.num_measurements,
            ));
        }

        // Prefer the GPU shot-major path: BTS plus on-device bit-transpose
        // into shot-major layout eliminates the host `into_shot_major_data()`
        // transpose that otherwise dominates noisy sampling at large shot
        // counts. Falls back to the default sampler (which may itself dispatch
        // to the GPU meas-major BTS) when the shot-major path does not apply.
        let mut accum = match self.noiseless.try_sample_bulk_shot_major_gpu(num_shots) {
            Some(Ok(data)) => data,
            Some(Err(e)) => return Err(e),
            None => {
                let (accum, _m_words) = self.sample_bulk_words_shot_major_cpu(num_shots);
                return Ok(PackedShots::from_shot_major(
                    accum,
                    num_shots,
                    self.num_measurements,
                ));
            }
        };
        self.apply_noise_bulk(&mut accum, num_shots, m_words);
        Ok(PackedShots::from_shot_major(
            accum,
            num_shots,
            self.num_measurements,
        ))
    }

    #[cfg(feature = "gpu")]
    fn take_gpu_noise_cache(&mut self) -> Result<GpuNoiseCache> {
        if let Some(cache) = self.gpu_noise_cache.take() {
            return Ok(cache);
        }
        let context = self
            .noiseless
            .gpu_context()
            .expect("gpu noise cache requested without gpu_context");
        GpuNoiseCache::new(context, &self.events)
    }

    #[cfg(feature = "gpu")]
    fn fill_noise_masks_event_major(
        events: &FlatNoiseSensitivity,
        rng: &mut ChaCha8Rng,
        x_masks: &mut [u64],
        z_masks: &mut [u64],
        chunk_shots: usize,
        chunk_s_words: usize,
    ) {
        x_masks.fill(0);
        z_masks.fill(0);

        // Events touch disjoint per-event slots of `x_masks` / `z_masks`, so
        // they parallelise cleanly once each worker has its own RNG stream.
        // Seed one Xoshiro per event from the master ChaCha so the draws are
        // independent of Rayon's thread count while matching the serial
        // reference statistically.
        #[cfg(feature = "parallel")]
        {
            use rand::Rng;
            use rayon::prelude::*;
            // Rayon split + per-event seed draw has to be dominated by real
            // work; below this threshold the serial master-rng loop wins.
            // Tuned against the noisy-sampling bench at 1M shots where
            // num_events is small and the tile is also small.
            const MIN_PAR_EVENTS: usize = 64;
            const MIN_PAR_SLOTS: usize = 32_768;
            let work_slots = events.len().saturating_mul(chunk_s_words);
            if events.len() >= MIN_PAR_EVENTS
                && work_slots >= MIN_PAR_SLOTS
                && x_masks.len() == events.len() * chunk_s_words
                && z_masks.len() == events.len() * chunk_s_words
            {
                let thread_seeds: Vec<[u64; 4]> = (0..events.len())
                    .map(|_| {
                        [
                            rng.next_u64(),
                            rng.next_u64(),
                            rng.next_u64(),
                            rng.next_u64(),
                        ]
                    })
                    .collect();

                x_masks
                    .par_chunks_mut(chunk_s_words)
                    .zip(z_masks.par_chunks_mut(chunk_s_words))
                    .zip(thread_seeds.par_iter())
                    .enumerate()
                    .for_each(|(i, ((x_slot, z_slot), seed))| {
                        let mut trng =
                            crate::sim::compiled::rng::Xoshiro256PlusPlus::from_seeds(*seed);
                        fill_one_event(x_slot, z_slot, events.probs[i], chunk_shots, &mut trng);
                    });
                return;
            }
        }

        for (i, &probs) in events.probs.iter().enumerate() {
            let base = i * chunk_s_words;
            let end = base + chunk_s_words;
            fill_one_event(
                &mut x_masks[base..end],
                &mut z_masks[base..end],
                probs,
                chunk_shots,
                rng,
            );
        }
    }

    #[cfg(feature = "gpu")]
    fn apply_noise_device_meas_major(
        &mut self,
        packed: &mut crate::sim::compiled::DevicePackedShots,
    ) -> Result<()> {
        if self.events.is_empty() || packed.num_shots() == 0 || self.num_measurements == 0 {
            return Ok(());
        }

        let num_events = self.events.len();
        let mut cache = self.take_gpu_noise_cache()?;
        let total_s_words = packed.s_words();
        let context = packed.context().clone();
        let num_shots = packed.num_shots();

        // Fused path: the kernel seeds a per-(event, batch) xoshiro stream
        // from a master seed drawn on the host, then generates and XORs the
        // noise masks in one launch per tile. Drops the per-tile mask fill
        // and ~25-50 MB PCIe upload that the fallback branch pays at large
        // shot counts. The fallback runs only when compiled event
        // probabilities fail the `u64` threshold check.
        let use_fused = cache.event_thresholds_valid;
        let master_seed: u64 = if use_fused {
            rand::Rng::next_u64(&mut self.rng)
        } else {
            0
        };

        let mut shot_start = 0;
        while shot_start < num_shots {
            let chunk_shots = (shot_start + NOISE_LUT_TILE).min(num_shots) - shot_start;
            let chunk_s_words = chunk_shots.div_ceil(64);
            let word_offset = shot_start / 64;

            if use_fused {
                crate::gpu::kernels::bts::generate_and_apply_noise_masks_meas_major_by_row(
                    &context,
                    crate::gpu::kernels::bts::NoiseDeviceGenApplyByRow {
                        meas_major: packed.data_mut(),
                        num_meas: self.num_measurements,
                        s_words: total_s_words,
                        word_offset,
                        chunk_s_words,
                        row_event_offsets: &cache.row_event_offsets_dev,
                        row_event_entries: &cache.row_event_entries_dev,
                        event_thresholds: &cache.event_thresholds_dev,
                        master_seed,
                        batch_offset: word_offset as u64,
                    },
                )?;
            } else {
                let mask_len = num_events * chunk_s_words;
                cache.ensure_mask_capacity(mask_len)?;
                Self::fill_noise_masks_event_major(
                    &self.events,
                    &mut self.rng,
                    &mut cache.x_masks_host[..mask_len],
                    &mut cache.z_masks_host[..mask_len],
                    chunk_shots,
                    chunk_s_words,
                );
                if mask_len > 0 {
                    cache
                        .x_masks_dev
                        .copy_from_host(cache.context.device(), &cache.x_masks_host[..mask_len])?;
                    cache
                        .z_masks_dev
                        .copy_from_host(cache.context.device(), &cache.z_masks_host[..mask_len])?;
                    let base = crate::gpu::kernels::bts::NoiseApplyBase {
                        meas_major: packed.data_mut(),
                        num_meas: self.num_measurements,
                        s_words: total_s_words,
                        word_offset,
                        chunk_s_words,
                        num_events,
                        x_row_offsets: &cache.x_row_offsets_dev,
                        x_row_indices: &cache.x_row_indices_dev,
                        z_row_offsets: &cache.z_row_offsets_dev,
                        z_row_indices: &cache.z_row_indices_dev,
                    };
                    crate::gpu::kernels::bts::apply_noise_masks_meas_major(
                        &context,
                        crate::gpu::kernels::bts::NoiseMaskApply {
                            base,
                            x_masks: &cache.x_masks_dev,
                            z_masks: &cache.z_masks_dev,
                        },
                    )?;
                }
            }

            shot_start += chunk_shots;
        }

        self.gpu_noise_cache = Some(cache);
        Ok(())
    }

    #[cfg(feature = "gpu")]
    fn try_sample_marginals_gpu(&mut self, total_shots: usize) -> Option<Result<Vec<f64>>> {
        if !self.readout.is_empty()
            || self.events.has_pairs()
            || !self.noiseless.should_use_gpu_bts(total_shots)
        {
            return None;
        }

        let mut packed = match self.noiseless.sample_bulk_packed_device(total_shots) {
            Ok(packed) => packed,
            Err(e) => return Some(Err(e)),
        };
        if let Err(e) = self.apply_noise_device_meas_major(&mut packed) {
            return Some(Err(e));
        }
        Some(packed.marginals())
    }

    #[cfg(feature = "gpu")]
    fn try_sample_counts_gpu(
        &mut self,
        total_shots: usize,
    ) -> Option<Result<std::collections::HashMap<Vec<u64>, u64>>> {
        if !self.readout.is_empty()
            || self.events.has_pairs()
            || !self.noiseless.should_use_gpu_bts(total_shots)
        {
            return None;
        }

        let mut packed = match self.noiseless.sample_bulk_packed_device(total_shots) {
            Ok(packed) => packed,
            Err(e) => return Some(Err(e)),
        };
        if let Err(e) = self.apply_noise_device_meas_major(&mut packed) {
            return Some(Err(e));
        }
        Some(packed.counts_with_rank_hint(self.num_measurements))
    }

    /// Sample `num_shots` noisy outcomes, packed shot-major.
    ///
    /// With a GPU context attached, parity sampling may run on the device; the
    /// packed output is materialized on the host either way.
    pub fn sample_bulk_packed(&mut self, num_shots: usize) -> PackedShots {
        match self.try_sample_bulk_packed(num_shots) {
            Ok(packed) => packed,
            Err(_) => self.sample_bulk_packed_cpu(num_shots),
        }
    }

    pub fn try_sample_bulk_packed(&mut self, num_shots: usize) -> Result<PackedShots> {
        #[cfg(feature = "gpu")]
        if self.noiseless.has_gpu_context() && self.readout.is_empty() {
            return self.try_sample_bulk_packed_with_gpu_context(num_shots);
        }

        Ok(self.sample_bulk_packed_cpu(num_shots))
    }

    fn sample_bulk_packed_cpu(&mut self, num_shots: usize) -> PackedShots {
        let m_words = self.num_measurements.div_ceil(64);
        if num_shots == 0 || self.num_measurements == 0 {
            return PackedShots::from_shot_major(
                vec![0u64; num_shots * m_words],
                num_shots,
                self.num_measurements,
            );
        }

        let (accum, _m_words) = self.sample_bulk_words_shot_major_cpu(num_shots);
        PackedShots::from_shot_major(accum, num_shots, self.num_measurements)
    }

    /// Stream noisy shots into `acc` at the default chunk size, without
    /// materializing the full shot matrix.
    pub fn sample_chunked<A: crate::sim::compiled::ShotAccumulator>(
        &mut self,
        total_shots: usize,
        acc: &mut A,
    ) {
        let chunk_size = crate::sim::compiled::default_chunk_size(self.num_measurements);
        self.sample_chunked_with_size(total_shots, chunk_size, acc);
    }

    /// Stream noisy shots into `acc` with an explicit chunk size.
    ///
    /// With a GPU context attached, each chunk may use GPU BTS sampling for the
    /// noiseless parity stage.
    pub fn sample_chunked_with_size<A: crate::sim::compiled::ShotAccumulator>(
        &mut self,
        total_shots: usize,
        chunk_size: usize,
        acc: &mut A,
    ) {
        #[cfg(feature = "gpu")]
        if self.noiseless.has_gpu_context()
            && self.readout.is_empty()
            && self
                .noiseless
                .should_use_gpu_bts(total_shots.min(chunk_size.max(1)))
        {
            crate::sim::compiled::for_each_chunk(total_shots, chunk_size, |this_batch| {
                let packed = match self.try_sample_bulk_packed_with_gpu_context(this_batch) {
                    Ok(packed) => packed,
                    Err(_) => self.sample_bulk_packed_cpu(this_batch),
                };
                acc.accumulate(&packed);
            });
            return;
        }

        let m_words = self.num_measurements.div_ceil(64);
        let mut accum_buf: Vec<u64> = Vec::new();
        let mut rand_buf: Vec<u8> = Vec::new();
        let ref_bits_packed = self.noiseless.ref_bits_packed().to_vec();
        crate::sim::compiled::for_each_chunk(total_shots, chunk_size, |this_batch| {
            self.noiseless.sample_bulk_words_shot_major_reuse(
                &mut accum_buf,
                &mut rand_buf,
                this_batch,
            );

            self.apply_noise_bulk(&mut accum_buf, this_batch, m_words);

            crate::sim::compiled::xor_mask_shot_major(
                &mut accum_buf,
                this_batch,
                m_words,
                &ref_bits_packed,
            );

            apply_readout_packed(
                &mut accum_buf,
                this_batch,
                m_words,
                &self.readout,
                &mut self.rng,
            );

            let data = std::mem::take(&mut accum_buf);
            let packed = PackedShots::from_shot_major(data, this_batch, self.num_measurements);
            acc.accumulate(&packed);
            accum_buf = packed.into_data();
        });
    }

    /// Sample noisy outcomes and return exact packed counts.
    ///
    /// A GPU failure falls back to the CPU sampler; [`Self::try_sample_counts`]
    /// surfaces it instead.
    pub fn sample_counts(
        &mut self,
        total_shots: usize,
    ) -> std::collections::HashMap<Vec<u64>, u64> {
        match self.try_sample_counts(total_shots) {
            Ok(counts) => counts,
            Err(_) => self.sample_counts_cpu(total_shots),
        }
    }

    pub fn try_sample_counts(
        &mut self,
        total_shots: usize,
    ) -> Result<std::collections::HashMap<Vec<u64>, u64>> {
        #[cfg(feature = "gpu")]
        if let Some(counts) = self.try_sample_counts_gpu(total_shots) {
            return counts;
        }

        Ok(self.sample_counts_cpu(total_shots))
    }

    fn sample_counts_cpu(
        &mut self,
        total_shots: usize,
    ) -> std::collections::HashMap<Vec<u64>, u64> {
        crate::sim::compiled::counts_from_chunks(|acc| self.sample_chunked(total_shots, acc))
    }

    /// Sample noisy outcomes and return per-measurement `P(bit = 1)`.
    pub fn sample_marginals(&mut self, total_shots: usize) -> Vec<f64> {
        match self.try_sample_marginals(total_shots) {
            Ok(marginals) => marginals,
            Err(_) => self.sample_marginals_cpu(total_shots),
        }
    }

    pub fn try_sample_marginals(&mut self, total_shots: usize) -> Result<Vec<f64>> {
        #[cfg(feature = "gpu")]
        if let Some(marginals) = self.try_sample_marginals_gpu(total_shots) {
            return marginals;
        }

        Ok(self.sample_marginals_cpu(total_shots))
    }

    fn sample_marginals_cpu(&mut self, total_shots: usize) -> Vec<f64> {
        crate::sim::compiled::marginals_from_chunks(self.num_measurements, |acc| {
            self.sample_chunked(total_shots, acc)
        })
    }

    #[cfg(test)]
    fn sample_bulk(&mut self, num_shots: usize) -> Vec<Vec<bool>> {
        self.sample_bulk_packed(num_shots).to_shots()
    }

    fn apply_noise_bulk(&mut self, accum: &mut [u64], num_shots: usize, m_words: usize) {
        if !self.events.is_empty() {
            if self.z_lut.is_some() {
                self.apply_noise_bulk_grouped(accum, num_shots, m_words);
            } else {
                self.apply_noise_bulk_scalar(accum, num_shots, m_words);
            }
        }

        // Two-qubit events run as a second pass over their own table, so the
        // single-qubit pass above keeps the row order the LUT and the device
        // buffers were built against.
        if self.events.has_pairs() {
            self.apply_pair_noise_bulk(accum, num_shots, m_words);
        }
    }

    fn apply_noise_bulk_scalar(&mut self, accum: &mut [u64], num_shots: usize, m_words: usize) {
        #[cfg(feature = "parallel")]
        if num_shots >= 4096 {
            self.apply_range_par(accum, num_shots, m_words, Self::apply_noise_range);
            return;
        }

        self.apply_noise_bulk_scalar_seq(accum, num_shots, m_words);
    }

    fn apply_noise_bulk_scalar_seq(&mut self, accum: &mut [u64], num_shots: usize, m_words: usize) {
        Self::apply_noise_range(&self.events, accum, 0, num_shots, m_words, &mut self.rng);
    }

    fn apply_pair_noise_bulk(&mut self, accum: &mut [u64], num_shots: usize, m_words: usize) {
        #[cfg(feature = "parallel")]
        if num_shots >= 4096 {
            self.apply_range_par(accum, num_shots, m_words, Self::apply_pair_noise_range);
            return;
        }

        Self::apply_pair_noise_range(&self.events, accum, 0, num_shots, m_words, &mut self.rng);
    }

    #[cfg(feature = "parallel")]
    fn apply_range_par(
        &mut self,
        accum: &mut [u64],
        num_shots: usize,
        m_words: usize,
        apply: fn(&FlatNoiseSensitivity, &mut [u64], usize, usize, usize, &mut ChaCha8Rng),
    ) {
        use rayon::prelude::*;

        let num_threads = rayon::current_num_threads().max(1);
        let shots_per_thread = (num_shots.div_ceil(num_threads) + 63) & !63;

        let seeds: Vec<u64> = (0..num_threads)
            .map(|_| rand::RngExt::random(&mut self.rng))
            .collect();

        let events = &self.events;

        accum
            .par_chunks_mut(shots_per_thread * m_words)
            .enumerate()
            .for_each(|(tid, chunk)| {
                let chunk_shots = chunk.len() / m_words;
                if chunk_shots == 0 {
                    return;
                }
                let mut rng = ChaCha8Rng::seed_from_u64(seeds[tid]);
                apply(events, chunk, 0, chunk_shots, m_words, &mut rng);
            });
    }

    fn apply_noise_range(
        events: &FlatNoiseSensitivity,
        accum: &mut [u64],
        start: usize,
        end: usize,
        m_words: usize,
        rng: &mut ChaCha8Rng,
    ) {
        let ne = events.len();
        let num_shots = end - start;

        for i in 0..ne {
            let [px, py, pz] = events.probs[i];
            let p_event = px + py + pz;
            if p_event == 0.0 {
                continue;
            }

            if p_event >= 0.5 || num_shots < 32 {
                for s in start..end {
                    let r: f64 = rand::RngExt::random(rng);
                    if r < px {
                        let b = s * m_words;
                        xor_words(&mut accum[b..b + m_words], events.z_flip(i));
                    } else if r < px + py {
                        let b = s * m_words;
                        xor_words(&mut accum[b..b + m_words], events.x_flip(i));
                        xor_words(&mut accum[b..b + m_words], events.z_flip(i));
                    } else if r < p_event {
                        let b = s * m_words;
                        xor_words(&mut accum[b..b + m_words], events.x_flip(i));
                    }
                }
            } else {
                let ln_1mp = (1.0 - p_event).ln();
                let px_frac = px / p_event;
                let pxy_frac = (px + py) / p_event;

                let mut pos = start + geometric_sample(rng, ln_1mp);
                while pos < end {
                    let r: f64 = rand::RngExt::random(rng);
                    let b = pos * m_words;
                    if r < px_frac {
                        xor_words(&mut accum[b..b + m_words], events.z_flip(i));
                    } else if r < pxy_frac {
                        xor_words(&mut accum[b..b + m_words], events.x_flip(i));
                        xor_words(&mut accum[b..b + m_words], events.z_flip(i));
                    } else {
                        xor_words(&mut accum[b..b + m_words], events.x_flip(i));
                    }
                    pos += 1 + geometric_sample(rng, ln_1mp);
                }
            }
        }
    }

    /// Draw the two-qubit events over `[start, end)`. One branch of 15 per
    /// firing shot, thinned exactly as the single-qubit pass thins its own
    /// events.
    fn apply_pair_noise_range(
        events: &FlatNoiseSensitivity,
        accum: &mut [u64],
        start: usize,
        end: usize,
        m_words: usize,
        rng: &mut ChaCha8Rng,
    ) {
        let num_shots = end - start;

        for i in 0..events.pair_probs.len() {
            let p = events.pair_probs[i];
            if p == 0.0 {
                continue;
            }

            if p >= 0.5 || num_shots < 32 {
                for s in start..end {
                    let r: f64 = rand::RngExt::random(rng);
                    if r < p {
                        let b = s * m_words;
                        xor_pair_branch(events, i, r / p, &mut accum[b..b + m_words]);
                    }
                }
            } else {
                let ln_1mp = (1.0 - p).ln();
                let mut pos = start + geometric_sample(rng, ln_1mp);
                while pos < end {
                    let r: f64 = rand::RngExt::random(rng);
                    let b = pos * m_words;
                    xor_pair_branch(events, i, r, &mut accum[b..b + m_words]);
                    pos += 1 + geometric_sample(rng, ln_1mp);
                }
            }
        }
    }

    fn apply_noise_bulk_grouped(&mut self, accum: &mut [u64], num_shots: usize, m_words: usize) {
        let num_events = self.events.len();
        let e_words = num_events.div_ceil(64);

        for tile_start in (0..num_shots).step_by(NOISE_LUT_TILE) {
            let tile_end = (tile_start + NOISE_LUT_TILE).min(num_shots);
            let tile_n = tile_end - tile_start;

            let mut z_mask = vec![0u64; tile_n * e_words];
            let mut x_mask = vec![0u64; tile_n * e_words];

            for i in 0..num_events {
                let [px, py, pz] = self.events.probs[i];
                let p_event = px + py + pz;
                if p_event == 0.0 {
                    continue;
                }

                let ew = i / 64;
                let eb = 1u64 << (i % 64);

                if p_event >= 0.5 || tile_n < 32 {
                    for s in 0..tile_n {
                        let r: f64 = rand::RngExt::random(&mut self.rng);
                        if r < px {
                            z_mask[s * e_words + ew] |= eb;
                        } else if r < px + py {
                            z_mask[s * e_words + ew] |= eb;
                            x_mask[s * e_words + ew] |= eb;
                        } else if r < p_event {
                            x_mask[s * e_words + ew] |= eb;
                        }
                    }
                } else {
                    let ln_1mp = (1.0 - p_event).ln();
                    let px_frac = px / p_event;
                    let pxy_frac = (px + py) / p_event;

                    let mut pos = geometric_sample(&mut self.rng, ln_1mp);
                    while pos < tile_n {
                        let r: f64 = rand::RngExt::random(&mut self.rng);
                        if r < px_frac {
                            z_mask[pos * e_words + ew] |= eb;
                        } else if r < pxy_frac {
                            z_mask[pos * e_words + ew] |= eb;
                            x_mask[pos * e_words + ew] |= eb;
                        } else {
                            x_mask[pos * e_words + ew] |= eb;
                        }
                        pos += 1 + geometric_sample(&mut self.rng, ln_1mp);
                    }
                }
            }

            let z_lut = self.z_lut.as_ref().unwrap();
            let x_lut = self.x_lut.as_ref().unwrap();

            for s in 0..tile_n {
                let shot_base = (tile_start + s) * m_words;
                let mask_base = s * e_words;
                let shot_accum = &mut accum[shot_base..shot_base + m_words];

                z_lut.apply_masked(shot_accum, &z_mask[mask_base..mask_base + e_words]);
                x_lut.apply_masked(shot_accum, &x_mask[mask_base..mask_base + e_words]);
            }
        }
    }
}

/// Compile a Clifford circuit with Pauli noise into a reusable noisy sampler.
///
/// Requires the same circuit subset as [`crate::compile_measurements`],
/// namely terminal measurements with no resets or classical conditionals.
pub fn compile_noisy(
    circuit: &Circuit,
    noise: &NoiseModel,
    seed: u64,
) -> Result<NoisyCompiledSampler> {
    noise.validate_for(circuit)?;
    noise.ensure_only_pauli_channels()?;
    if circuit.num_qubits >= 4 {
        let blocks = circuit.independent_subsystems();
        if blocks.len() > 1 {
            let max_block = blocks.iter().map(|b| b.len()).max().unwrap_or(0);
            if max_block < circuit.num_qubits && pairs_stay_in_one_block(circuit, noise, &blocks) {
                return compile_noisy_filtered(circuit, noise, &blocks, seed);
            }
        }
    }

    compile_noisy_monolithic(circuit, noise, seed)
}

/// True when any measurement is still sensitive to a Pauli on `q` at this
/// point in the backward sweep.
fn sensitive(x_packed: &[Vec<u64>], z_packed: &[Vec<u64>], q: usize) -> bool {
    x_packed[q].iter().any(|&w| w != 0) || z_packed[q].iter().any(|&w| w != 0)
}

/// Scatter a block-local measurement mask into the global measurement order,
/// which is what the sensitivity table is indexed by.
fn scatter_to_global(local: &[u64], bm_list: &[usize], out: &mut [u64]) {
    out.fill(0);
    for (local_mi, &global_mi) in bm_list.iter().enumerate() {
        if (local[local_mi / 64] >> (local_mi % 64)) & 1 != 0 {
            out[global_mi / 64] |= 1u64 << (global_mi % 64);
        }
    }
}

/// True when every two-qubit event lands inside the subsystem block holding
/// the gate it follows.
///
/// The filtered compile sweeps one block at a time and holds only that
/// block's propagated masks, so a pair reaching across blocks has nowhere to
/// read its second qubit from and would silently run at the rate of its first
/// qubit alone. Such models take the monolithic compile, which propagates the
/// whole register at once. A zero-rate pair stores no row and constrains
/// nothing.
fn pairs_stay_in_one_block(circuit: &Circuit, noise: &NoiseModel, blocks: &[Vec<usize>]) -> bool {
    let mut qubit_to_block = vec![usize::MAX; circuit.num_qubits];
    for (bi, block) in blocks.iter().enumerate() {
        for &q in block {
            qubit_to_block[q] = bi;
        }
    }

    for (inst, events) in circuit.instructions.iter().zip(&noise.after_gate) {
        let mut pairs = events
            .iter()
            .filter_map(NoiseEvent::pauli_pair)
            .filter(|&(_, _, p)| p > 0.0)
            .peekable();
        if pairs.peek().is_none() {
            continue;
        }

        let targets = match inst {
            Instruction::Gate { targets, .. } | Instruction::Conditional { targets, .. } => targets,
            _ => return false,
        };
        let Some(&first) = targets.first() else {
            return false;
        };
        let block = qubit_to_block[first];
        if !pairs.all(|(q0, q1, _)| qubit_to_block[q0] == block && qubit_to_block[q1] == block) {
            return false;
        }
    }
    true
}

fn compile_noisy_filtered(
    circuit: &Circuit,
    noise: &NoiseModel,
    blocks: &[Vec<usize>],
    seed: u64,
) -> Result<NoisyCompiledSampler> {
    let noiseless = compile_measurements(circuit, seed)?;

    let measurement_qubits: Vec<usize> = circuit
        .instructions
        .iter()
        .filter_map(|inst| match inst {
            Instruction::Measure { qubit, .. } => Some(*qubit),
            _ => None,
        })
        .collect();
    let num_measurements = measurement_qubits.len();
    let m_words = num_measurements.div_ceil(64);

    if num_measurements == 0 {
        return Ok(NoisyCompiledSampler {
            noiseless,
            events: FlatNoiseSensitivity::new(1, 0),
            readout: Vec::new(),
            num_measurements: 0,
            rng: ChaCha8Rng::seed_from_u64(seed.wrapping_add(0xA01CE)),
            z_lut: None,
            x_lut: None,
            #[cfg(feature = "gpu")]
            gpu_noise_cache: None,
        });
    }

    let mut qubit_to_block = vec![0usize; circuit.num_qubits];
    let mut qubit_to_local = vec![0usize; circuit.num_qubits];
    for (bi, block) in blocks.iter().enumerate() {
        for (li, &q) in block.iter().enumerate() {
            qubit_to_block[q] = bi;
            qubit_to_local[q] = li;
        }
    }

    let mut block_meas: Vec<Vec<usize>> = vec![Vec::new(); blocks.len()];
    for (mi, &q) in measurement_qubits.iter().enumerate() {
        block_meas[qubit_to_block[q]].push(mi);
    }

    let total_noise_events: usize = noise.after_gate.iter().map(|ops| ops.len()).sum();
    let mut events = FlatNoiseSensitivity::new(m_words, total_noise_events);
    let mut global_x_buf = vec![0u64; m_words];
    let mut global_z_buf = vec![0u64; m_words];
    let mut pair_x_buf = vec![0u64; m_words];
    let mut pair_z_buf = vec![0u64; m_words];

    for (bi, block) in blocks.iter().enumerate() {
        let bm_list = &block_meas[bi];
        if bm_list.is_empty() {
            continue;
        }

        let bn = block.len();
        let bm = bm_list.len();
        let bm_words = bm.div_ceil(64);

        let mut x_packed: Vec<Vec<u64>> = vec![vec![0u64; bm_words]; bn];
        let mut z_packed: Vec<Vec<u64>> = vec![vec![0u64; bm_words]; bn];
        let mut sign_packed: Vec<u64> = vec![0u64; bm_words];

        for (local_mi, &global_mi) in bm_list.iter().enumerate() {
            let q = measurement_qubits[global_mi];
            let local_q = qubit_to_local[q];
            z_packed[local_q][local_mi / 64] |= 1u64 << (local_mi % 64);
        }

        let block_gates: Vec<(usize, &Gate, SmallVec<[usize; 4]>)> = circuit
            .instructions
            .iter()
            .enumerate()
            .filter_map(|(idx, inst)| match inst {
                Instruction::Gate { gate, targets } => {
                    if targets.iter().all(|&t| qubit_to_block[t] == bi) {
                        let local_targets: SmallVec<[usize; 4]> =
                            targets.iter().map(|&t| qubit_to_local[t]).collect();
                        Some((idx, gate, local_targets))
                    } else {
                        None
                    }
                }
                Instruction::Conditional { gate, targets, .. } => {
                    if targets.iter().all(|&t| qubit_to_block[t] == bi) {
                        let local_targets: SmallVec<[usize; 4]> =
                            targets.iter().map(|&t| qubit_to_local[t]).collect();
                        Some((idx, gate, local_targets))
                    } else {
                        None
                    }
                }
                _ => None,
            })
            .collect();

        for (instr_idx, gate, local_targets) in block_gates.iter().rev() {
            let noise_ops = &noise.after_gate[*instr_idx];
            if !noise_ops.is_empty() {
                for event in noise_ops {
                    if let Some((q0, q1, p)) = event.pauli_pair() {
                        if p == 0.0 {
                            continue;
                        }
                        let (l0, l1) = (qubit_to_local[q0], qubit_to_local[q1]);
                        if !sensitive(&x_packed, &z_packed, l0)
                            && !sensitive(&x_packed, &z_packed, l1)
                        {
                            continue;
                        }
                        scatter_to_global(&x_packed[l0], bm_list, &mut global_x_buf);
                        scatter_to_global(&z_packed[l0], bm_list, &mut global_z_buf);
                        scatter_to_global(&x_packed[l1], bm_list, &mut pair_x_buf);
                        scatter_to_global(&z_packed[l1], bm_list, &mut pair_z_buf);
                        events.push_pair(&global_x_buf, &global_z_buf, &pair_x_buf, &pair_z_buf, p);
                        continue;
                    }

                    let (px, py, pz) = event.pauli_probs();
                    let local_q = qubit_to_local[event.qubit()];

                    if sensitive(&x_packed, &z_packed, local_q) {
                        global_x_buf.fill(0);
                        global_z_buf.fill(0);
                        for (local_mi, &global_mi) in bm_list.iter().enumerate() {
                            if (x_packed[local_q][local_mi / 64] >> (local_mi % 64)) & 1 != 0 {
                                global_x_buf[global_mi / 64] |= 1u64 << (global_mi % 64);
                            }
                            if (z_packed[local_q][local_mi / 64] >> (local_mi % 64)) & 1 != 0 {
                                global_z_buf[global_mi / 64] |= 1u64 << (global_mi % 64);
                            }
                        }
                        events.push(&global_x_buf, &global_z_buf, px, py, pz);
                    }
                }
            }

            batch_propagate_backward(
                &mut x_packed,
                &mut z_packed,
                &mut sign_packed,
                gate,
                local_targets.as_slice(),
                bm_words,
            );
        }
    }

    let (z_lut, x_lut) = build_noise_luts(&events);

    Ok(NoisyCompiledSampler {
        noiseless,
        events,
        readout: record_readout(circuit, noise),
        num_measurements,
        rng: ChaCha8Rng::seed_from_u64(seed.wrapping_add(0xCAFE_BABE)),
        z_lut,
        x_lut,
        #[cfg(feature = "gpu")]
        gpu_noise_cache: None,
    })
}

fn compile_noisy_monolithic(
    circuit: &Circuit,
    noise: &NoiseModel,
    seed: u64,
) -> Result<NoisyCompiledSampler> {
    let noiseless = compile_measurements(circuit, seed)?;

    let n = circuit.num_qubits;

    let gate_indices: Vec<(usize, &Gate, &[usize])> = circuit
        .instructions
        .iter()
        .enumerate()
        .filter_map(|(idx, inst)| match inst {
            Instruction::Gate { gate, targets } => Some((idx, gate, targets.as_slice())),
            Instruction::Conditional { gate, targets, .. } => Some((idx, gate, targets.as_slice())),
            _ => None,
        })
        .collect();

    let measurement_qubits: Vec<usize> = circuit
        .instructions
        .iter()
        .filter_map(|inst| match inst {
            Instruction::Measure { qubit, .. } => Some(*qubit),
            _ => None,
        })
        .collect();
    let num_measurements = measurement_qubits.len();
    let m_words = num_measurements.div_ceil(64);

    if num_measurements == 0 {
        return Ok(NoisyCompiledSampler {
            noiseless,
            events: FlatNoiseSensitivity::new(m_words, 0),
            readout: Vec::new(),
            num_measurements: 0,
            rng: ChaCha8Rng::seed_from_u64(seed.wrapping_add(0xA01CE)),
            z_lut: None,
            x_lut: None,
            #[cfg(feature = "gpu")]
            gpu_noise_cache: None,
        });
    }

    let mut x_packed: Vec<Vec<u64>> = vec![vec![0u64; m_words]; n];
    let mut z_packed: Vec<Vec<u64>> = vec![vec![0u64; m_words]; n];
    let mut sign_packed: Vec<u64> = vec![0u64; m_words];

    for (mi, &q) in measurement_qubits.iter().enumerate() {
        z_packed[q][mi / 64] |= 1u64 << (mi % 64);
    }

    let total_noise_events: usize = noise.after_gate.iter().map(|ops| ops.len()).sum();
    let mut events = FlatNoiseSensitivity::new(m_words, total_noise_events);

    for &(instr_idx, gate, targets) in gate_indices.iter().rev() {
        let noise_ops = &noise.after_gate[instr_idx];
        if !noise_ops.is_empty() {
            for event in noise_ops {
                if let Some((q0, q1, p)) = event.pauli_pair() {
                    if p == 0.0 {
                        continue;
                    }
                    if sensitive(&x_packed, &z_packed, q0) || sensitive(&x_packed, &z_packed, q1) {
                        events.push_pair(
                            &x_packed[q0],
                            &z_packed[q0],
                            &x_packed[q1],
                            &z_packed[q1],
                            p,
                        );
                    }
                    continue;
                }

                let (px, py, pz) = event.pauli_probs();
                let q = event.qubit();

                if sensitive(&x_packed, &z_packed, q) {
                    events.push(&x_packed[q], &z_packed[q], px, py, pz);
                }
            }
        }

        batch_propagate_backward(
            &mut x_packed,
            &mut z_packed,
            &mut sign_packed,
            gate,
            targets,
            m_words,
        );
    }

    let (z_lut, x_lut) = build_noise_luts(&events);

    Ok(NoisyCompiledSampler {
        noiseless,
        events,
        readout: record_readout(circuit, noise),
        num_measurements,
        rng: ChaCha8Rng::seed_from_u64(seed.wrapping_add(0xCAFE_BABE)),
        z_lut,
        x_lut,
        #[cfg(feature = "gpu")]
        gpu_noise_cache: None,
    })
}

const FRAME_BATCH_SIZE: usize = 256;

/// Gate count per qubit below which the Pauli frame sampler beats the compiled
/// sensitivity map, whose build cost amortizes only over deeper circuits. Wide
/// circuits carry a higher cutoff: the frame walks one word per qubit per gate,
/// so its per-gate cost grows with the register while the map's does not.
const FRAME_DEPTH_RATIO: f64 = 3.0;
const WIDE_FRAME_QUBITS: usize = 200;
const WIDE_FRAME_DEPTH_RATIO: f64 = 5.0;

/// True when [`run_shots_noisy`] takes the Pauli frame sampler rather than the
/// compiled one. Both answer from the same distribution and stamp the same
/// resolved backend, so a result carries nothing that separates them.
fn use_frame_sampler(circuit: &Circuit) -> bool {
    let n = circuit.num_qubits;
    let gate_count = circuit
        .instructions
        .iter()
        .filter(|i| {
            matches!(
                i,
                Instruction::Gate { .. } | Instruction::Conditional { .. }
            )
        })
        .count();

    let depth_ratio = gate_count as f64 / n.max(1) as f64;
    depth_ratio < FRAME_DEPTH_RATIO
        || (n >= WIDE_FRAME_QUBITS && depth_ratio < WIDE_FRAME_DEPTH_RATIO)
}

struct ReferenceInfo {
    outcomes: Vec<bool>,
    is_random: Vec<bool>,
    random_x_support: Vec<Vec<usize>>,
}

fn reference_simulation(circuit: &Circuit, seed: u64) -> Result<ReferenceInfo> {
    let mut stab = StabilizerBackend::new(seed);
    stab.init(circuit.num_qubits, circuit.num_classical_bits)?;

    let meas_info: Vec<(usize, usize, usize)> = circuit
        .instructions
        .iter()
        .enumerate()
        .filter_map(|(i, inst)| match inst {
            Instruction::Measure {
                qubit,
                classical_bit,
            } => Some((i, *qubit, *classical_bit)),
            _ => None,
        })
        .collect();

    let num_meas = meas_info.len();
    if num_meas == 0 {
        return Ok(ReferenceInfo {
            outcomes: Vec::new(),
            is_random: Vec::new(),
            random_x_support: Vec::new(),
        });
    }

    let first_meas_idx = meas_info[0].0;
    let all_at_end = meas_info
        .iter()
        .enumerate()
        .all(|(i, &(inst_idx, _, _))| inst_idx == first_meas_idx + i);
    let mut seen_qubits = vec![false; circuit.num_qubits];
    let mut seen_bits = vec![false; circuit.num_classical_bits];
    let distinct = meas_info.iter().all(|&(_, qubit, classical_bit)| {
        !std::mem::replace(&mut seen_qubits[qubit], true)
            && !std::mem::replace(&mut seen_bits[classical_bit], true)
    });

    if all_at_end && distinct {
        stab.apply_gates_only(&circuit.instructions[..first_meas_idx])?;

        let measurements: Vec<(usize, usize)> = meas_info
            .iter()
            .map(|&(_, qubit, classical_bit)| (qubit, classical_bit))
            .collect();
        let (is_random, random_x_support, outcomes) = stab.batch_measure_ref_info(&measurements);

        return Ok(ReferenceInfo {
            outcomes,
            is_random,
            random_x_support,
        });
    }

    let mut outcomes = Vec::with_capacity(num_meas);
    let mut is_random = Vec::with_capacity(num_meas);
    let mut random_x_support: Vec<Vec<usize>> = Vec::with_capacity(num_meas);

    let mut seg_start = 0usize;
    for &(meas_inst_idx, qubit, classical_bit) in &meas_info {
        if seg_start < meas_inst_idx {
            stab.apply_gates_only(&circuit.instructions[seg_start..meas_inst_idx])?;
        }

        let (meas_random, support) = stab.apply_measure_with_info(qubit, classical_bit);
        outcomes.push(stab.classical_results()[classical_bit]);
        is_random.push(meas_random);
        random_x_support.push(support);
        seg_start = meas_inst_idx + 1;
    }

    Ok(ReferenceInfo {
        outcomes,
        is_random,
        random_x_support,
    })
}

#[inline(always)]
fn apply_gate_to_frame(
    gate: &Gate,
    targets: &[usize],
    x_frame: &mut [Vec<u64>],
    z_frame: &mut [Vec<u64>],
    bw: usize,
) {
    match gate {
        Gate::H => {
            let q = targets[0];
            std::mem::swap(&mut x_frame[q], &mut z_frame[q]);
        }
        Gate::S | Gate::Sdg => {
            let q = targets[0];
            for w in 0..bw {
                z_frame[q][w] ^= x_frame[q][w];
            }
        }
        Gate::SX | Gate::SXdg => {
            let q = targets[0];
            for w in 0..bw {
                x_frame[q][w] ^= z_frame[q][w];
            }
        }
        Gate::X | Gate::Y | Gate::Z | Gate::Id => {}
        Gate::Cx => {
            let ctrl = targets[0];
            let tgt = targets[1];
            for w in 0..bw {
                x_frame[tgt][w] ^= x_frame[ctrl][w];
                z_frame[ctrl][w] ^= z_frame[tgt][w];
            }
        }
        Gate::Cz => {
            let q0 = targets[0];
            let q1 = targets[1];
            for w in 0..bw {
                z_frame[q0][w] ^= x_frame[q1][w];
                z_frame[q1][w] ^= x_frame[q0][w];
            }
        }
        Gate::Swap => {
            let q0 = targets[0];
            let q1 = targets[1];
            x_frame.swap(q0, q1);
            z_frame.swap(q0, q1);
        }
        _ => {
            debug_assert!(
                false,
                "apply_gate_to_frame: unhandled Clifford gate {:?}",
                gate
            );
        }
    }
}

#[inline(always)]
fn flip_frame_letter(
    x_frame: &mut [Vec<u64>],
    z_frame: &mut [Vec<u64>],
    q: usize,
    pos: usize,
    letter: usize,
) {
    let w = pos / 64;
    let bit = 1u64 << (pos % 64);
    match letter {
        1 => x_frame[q][w] ^= bit,
        2 => {
            x_frame[q][w] ^= bit;
            z_frame[q][w] ^= bit;
        }
        3 => z_frame[q][w] ^= bit,
        _ => {}
    }
}

fn run_shots_noisy_frame(
    circuit: &Circuit,
    noise: &NoiseModel,
    num_shots: usize,
    seed: u64,
) -> Result<ShotsResult> {
    let n = circuit.num_qubits;
    let num_classical = circuit.num_classical_bits;

    let ref_info = reference_simulation(circuit, seed)?;

    let classical_bit_order = circuit.classical_bit_order();
    let num_measurements = classical_bit_order.len();
    let m_words = num_measurements.div_ceil(64);

    let mut all_packed = vec![0u64; num_shots * m_words];
    let readout = record_readout(circuit, noise);
    let mut rng = ChaCha8Rng::seed_from_u64(seed.wrapping_add(0xFAAB_E001));

    let max_bw = FRAME_BATCH_SIZE.div_ceil(64);
    let mut x_frame: Vec<Vec<u64>> = vec![vec![0u64; max_bw]; n];
    let mut z_frame: Vec<Vec<u64>> = vec![vec![0u64; max_bw]; n];

    for batch_start in (0..num_shots).step_by(FRAME_BATCH_SIZE) {
        let batch_end = (batch_start + FRAME_BATCH_SIZE).min(num_shots);
        let batch_n = batch_end - batch_start;
        let bw = batch_n.div_ceil(64);

        for row in x_frame.iter_mut().chain(z_frame.iter_mut()) {
            row[..bw].fill(0);
        }

        let mut meas_idx = 0usize;

        for (idx, inst) in circuit.instructions.iter().enumerate() {
            match inst {
                Instruction::Gate { gate, targets }
                | Instruction::Conditional { gate, targets, .. } => {
                    apply_gate_to_frame(gate, targets.as_slice(), &mut x_frame, &mut z_frame, bw);

                    for event in &noise.after_gate[idx] {
                        if let Some((q0, q1, p)) = event.pauli_pair() {
                            if p == 0.0 {
                                continue;
                            }
                            if p < 0.5 && batch_n >= 32 {
                                let ln_1mp = (1.0 - p).ln();
                                let mut pos = geometric_sample(&mut rng, ln_1mp);
                                while pos < batch_n {
                                    let r: f64 = rand::RngExt::random(&mut rng);
                                    let (l0, l1) = pair_branch(r);
                                    flip_frame_letter(&mut x_frame, &mut z_frame, q0, pos, l0);
                                    flip_frame_letter(&mut x_frame, &mut z_frame, q1, pos, l1);
                                    pos += 1 + geometric_sample(&mut rng, ln_1mp);
                                }
                            } else {
                                for s in 0..batch_n {
                                    let r: f64 = rand::RngExt::random(&mut rng);
                                    if r < p {
                                        let (l0, l1) = pair_branch(r / p);
                                        flip_frame_letter(&mut x_frame, &mut z_frame, q0, s, l0);
                                        flip_frame_letter(&mut x_frame, &mut z_frame, q1, s, l1);
                                    }
                                }
                            }
                            continue;
                        }

                        let (px, py, pz) = event.pauli_probs();
                        let q = event.qubit();
                        let p_event = px + py + pz;
                        if p_event == 0.0 {
                            continue;
                        }

                        let px_frac = px / p_event;
                        let pxy_frac = (px + py) / p_event;

                        if p_event < 0.5 && batch_n >= 32 {
                            let ln_1mp = (1.0 - p_event).ln();
                            let mut pos = geometric_sample(&mut rng, ln_1mp);
                            while pos < batch_n {
                                let r: f64 = rand::RngExt::random(&mut rng);
                                let bit = 1u64 << (pos % 64);
                                let w = pos / 64;
                                if r < px_frac {
                                    x_frame[q][w] ^= bit;
                                } else if r < pxy_frac {
                                    x_frame[q][w] ^= bit;
                                    z_frame[q][w] ^= bit;
                                } else {
                                    z_frame[q][w] ^= bit;
                                }
                                pos += 1 + geometric_sample(&mut rng, ln_1mp);
                            }
                        } else {
                            for s in 0..batch_n {
                                let r: f64 = rand::RngExt::random(&mut rng);
                                if r < px {
                                    x_frame[q][s / 64] ^= 1u64 << (s % 64);
                                } else if r < px + py {
                                    x_frame[q][s / 64] ^= 1u64 << (s % 64);
                                    z_frame[q][s / 64] ^= 1u64 << (s % 64);
                                } else if r < p_event {
                                    z_frame[q][s / 64] ^= 1u64 << (s % 64);
                                }
                            }
                        }
                    }
                }
                Instruction::Measure {
                    qubit,
                    classical_bit: _,
                } => {
                    if ref_info.is_random[meas_idx] {
                        let support = &ref_info.random_x_support[meas_idx];
                        #[allow(clippy::needless_range_loop)]
                        for w in 0..bw {
                            let random_word: u64 = rand::RngExt::random(&mut rng);
                            let mask = if w == bw - 1 && !batch_n.is_multiple_of(64) {
                                random_word & ((1u64 << (batch_n % 64)) - 1)
                            } else {
                                random_word
                            };
                            if mask != 0 {
                                for &q in support {
                                    x_frame[q][w] ^= mask;
                                }
                            }
                        }
                    }

                    let ref_bit = ref_info.outcomes[meas_idx];
                    let mi_word = meas_idx / 64;
                    let mi_bit = meas_idx % 64;
                    #[allow(clippy::needless_range_loop)]
                    for w in 0..bw {
                        let frame_word = x_frame[*qubit][w];
                        let effective = if ref_bit { !frame_word } else { frame_word };
                        let num_bits = if w == bw - 1 && !batch_n.is_multiple_of(64) {
                            batch_n % 64
                        } else {
                            64
                        };
                        let mask = if num_bits == 64 {
                            effective
                        } else {
                            effective & ((1u64 << num_bits) - 1)
                        };
                        let mut bits = mask;
                        while bits != 0 {
                            let s = bits.trailing_zeros() as usize;
                            let gs = batch_start + w * 64 + s;
                            all_packed[gs * m_words + mi_word] |= 1u64 << mi_bit;
                            bits &= bits - 1;
                        }
                    }

                    meas_idx += 1;
                }
                _ => {}
            }
        }
    }

    apply_readout_packed(&mut all_packed, num_shots, m_words, &readout, &mut rng);

    let shots = unpack_and_remap(
        &all_packed,
        m_words,
        num_shots,
        &classical_bit_order,
        num_classical,
    );

    Ok(
        ShotsResult::from_shots(shots, circuit.num_classical_bits).with_metadata(
            crate::sim::RunMetadata::exact(crate::sim::ResolvedBackend::CompiledStabilizer)
                .with_engine(crate::sim::Engine::FrameSampler),
        ),
    )
}

/// Sample a Pauli-noisy Clifford circuit through the compiled noisy sampler or
/// the frame sampler.
///
/// Circuits with resets, classical conditionals, or mid-circuit measurements
/// fall back to per-shot stabilizer replay. Non-Clifford circuits take the
/// same fallback and error on the first non-Clifford gate; the public entry
/// point routes them to the trajectory engine before reaching here.
pub fn run_shots_noisy(
    circuit: &Circuit,
    noise: &NoiseModel,
    num_shots: usize,
    seed: u64,
) -> Result<ShotsResult> {
    noise.validate_for(circuit)?;
    noise.ensure_only_pauli_channels()?;
    if !circuit.is_clifford_only() {
        return run_shots_noisy_brute_with(
            |s| Box::new(StabilizerBackend::new(s)),
            circuit,
            noise,
            num_shots,
            seed,
        );
    }

    if !super::supports_compiled_measurement_sampling(circuit) {
        return run_shots_noisy_brute_with(
            |s| Box::new(StabilizerBackend::new(s)),
            circuit,
            noise,
            num_shots,
            seed,
        );
    }

    if use_frame_sampler(circuit) {
        run_shots_noisy_frame(circuit, noise, num_shots, seed)
    } else {
        run_shots_noisy_compiled(circuit, noise, num_shots, seed)
    }
}

fn finish_noisy_compiled_run(
    mut sampler: NoisyCompiledSampler,
    circuit: &Circuit,
    num_shots: usize,
) -> Result<ShotsResult> {
    let classical_bit_order = circuit.classical_bit_order();
    let num_classical = circuit.num_classical_bits;

    let packed = sampler.try_sample_bulk_packed(num_shots)?;
    let shots = unpack_and_remap_packed(&packed, num_shots, &classical_bit_order, num_classical);

    Ok(
        ShotsResult::from_shots(shots, circuit.num_classical_bits).with_metadata(
            crate::sim::RunMetadata::exact(crate::sim::ResolvedBackend::CompiledStabilizer)
                .with_engine(crate::sim::Engine::NoisyCompiledSampler),
        ),
    )
}

fn run_shots_noisy_compiled(
    circuit: &Circuit,
    noise: &NoiseModel,
    num_shots: usize,
    seed: u64,
) -> Result<ShotsResult> {
    let sampler = compile_noisy(circuit, noise, seed)?;
    finish_noisy_compiled_run(sampler, circuit, num_shots)
}

#[cfg(feature = "gpu")]
pub(crate) fn run_shots_noisy_with_gpu(
    circuit: &Circuit,
    noise: &NoiseModel,
    num_shots: usize,
    seed: u64,
    context: std::sync::Arc<crate::gpu::GpuContext>,
) -> Result<ShotsResult> {
    noise.ensure_only_pauli_channels()?;
    if !circuit.is_clifford_only() {
        return run_shots_noisy_brute_with(
            |s| Box::new(StabilizerBackend::new(s)),
            circuit,
            noise,
            num_shots,
            seed,
        );
    }

    if !super::supports_compiled_measurement_sampling(circuit) {
        return run_shots_noisy_brute_with(
            |s| Box::new(StabilizerBackend::new(s)),
            circuit,
            noise,
            num_shots,
            seed,
        );
    }

    if use_frame_sampler(circuit) {
        return run_shots_noisy_frame(circuit, noise, num_shots, seed);
    }

    let sampler = compile_noisy(circuit, noise, seed)?.with_gpu(context);
    finish_noisy_compiled_run(sampler, circuit, num_shots)
}

pub(crate) fn run_shots_noisy_brute_with(
    backend_factory: impl Fn(u64) -> Box<dyn Backend>,
    circuit: &Circuit,
    noise: &NoiseModel,
    num_shots: usize,
    seed: u64,
) -> Result<ShotsResult> {
    let mut shots = Vec::with_capacity(num_shots);
    let mut metadata = crate::sim::RunMetadata::exact(crate::sim::ResolvedBackend::Stabilizer);
    let readout = crate::sim::trajectory::written_readout(circuit, &noise.readout);

    for i in 0..num_shots {
        let shot_seed = crate::sim::mix_seed(seed, i);
        let mut rng = crate::sim::trajectory::noise_rng(shot_seed);
        let mut backend = backend_factory(shot_seed);
        backend.init(circuit.num_qubits, circuit.num_classical_bits)?;

        for (idx, instr) in circuit.instructions.iter().enumerate() {
            backend.apply(instr)?;

            let noise_events = &noise.after_gate[idx];
            for event in noise_events {
                if let Some((q0, q1, p)) = event.pauli_pair() {
                    let r: f64 = rand::RngExt::random(&mut rng);
                    if r < p {
                        let (l0, l1) = pair_branch(r / p);
                        apply_pauli_letter(backend.as_mut(), q0, l0)?;
                        apply_pauli_letter(backend.as_mut(), q1, l1)?;
                    }
                    continue;
                }
                let (px, py, pz) = event.pauli_probs();
                let q = event.qubit();
                let r: f64 = rand::RngExt::random(&mut rng);
                if r < px {
                    backend.apply(&Instruction::Gate {
                        gate: Gate::X,
                        targets: SmallVec::from_elem(q, 1),
                    })?;
                } else if r < px + py {
                    backend.apply(&Instruction::Gate {
                        gate: Gate::Y,
                        targets: SmallVec::from_elem(q, 1),
                    })?;
                } else if r < px + py + pz {
                    backend.apply(&Instruction::Gate {
                        gate: Gate::Z,
                        targets: SmallVec::from_elem(q, 1),
                    })?;
                }
            }
        }

        let shot_metadata = crate::sim::backend_metadata(backend.as_ref());
        if i == 0 {
            metadata = shot_metadata;
        } else {
            metadata.weaken_with(&shot_metadata);
        }
        let mut results = backend.classical_results().to_vec();
        crate::sim::trajectory::apply_readout_errors(&mut results, &readout, &mut rng);
        shots.push(results);
    }

    Ok(ShotsResult::from_shots(shots, circuit.num_classical_bits).with_metadata(metadata))
}

/// Apply one letter of a two-qubit branch, indexed as [`pair_branch`] returns
/// it. Identity applies nothing.
fn apply_pauli_letter(backend: &mut dyn Backend, qubit: usize, letter: usize) -> Result<()> {
    let gate = match letter {
        1 => Gate::X,
        2 => Gate::Y,
        3 => Gate::Z,
        _ => return Ok(()),
    };
    backend.apply(&Instruction::Gate {
        gate,
        targets: SmallVec::from_elem(qubit, 1),
    })
}

fn amplitude_damping_kraus(gamma: f64) -> [[[Complex64; 2]; 2]; 2] {
    let s = (1.0 - gamma).max(0.0).sqrt();
    let g = gamma.max(0.0).sqrt();
    [
        [
            [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
            [Complex64::new(0.0, 0.0), Complex64::new(s, 0.0)],
        ],
        [
            [Complex64::new(0.0, 0.0), Complex64::new(g, 0.0)],
            [Complex64::new(0.0, 0.0), Complex64::new(0.0, 0.0)],
        ],
    ]
}

/// Kraus set of generalized amplitude damping: decay toward `|0>` at rate
/// `gamma` weighted by `1 - excited`, excitation toward `|1>` at the same rate
/// weighted by `excited`. At `excited == 0` the last two operators vanish and
/// the set is the amplitude-damping pair.
fn generalized_amplitude_damping_kraus(gamma: f64, excited: f64) -> [[[Complex64; 2]; 2]; 4] {
    let s = (1.0 - gamma).max(0.0).sqrt();
    let g = gamma.max(0.0).sqrt();
    let cold = (1.0 - excited).max(0.0).sqrt();
    let hot = excited.max(0.0).sqrt();
    let c = |re: f64| Complex64::new(re, 0.0);
    [
        [[c(cold), c(0.0)], [c(0.0), c(cold * s)]],
        [[c(0.0), c(cold * g)], [c(0.0), c(0.0)]],
        [[c(hot * s), c(0.0)], [c(0.0), c(hot)]],
        [[c(0.0), c(0.0)], [c(hot * g), c(0.0)]],
    ]
}

fn phase_damping_kraus(gamma: f64) -> [[[Complex64; 2]; 2]; 2] {
    let s = (1.0 - gamma).max(0.0).sqrt();
    let g = gamma.max(0.0).sqrt();
    [
        [
            [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
            [Complex64::new(0.0, 0.0), Complex64::new(s, 0.0)],
        ],
        [
            [Complex64::new(0.0, 0.0), Complex64::new(0.0, 0.0)],
            [Complex64::new(0.0, 0.0), Complex64::new(g, 0.0)],
        ],
    ]
}

/// Amplitude-damping and phase-damping rates whose composition decays
/// populations as `exp(-t/t1)` and coherences as `exp(-t/t2)`. Shared by the
/// exact Kraus lowering and the trajectory unraveling so the two routes cannot
/// disagree on what a declared channel means. Requires `t2 <= 2*t1`, which
/// [`NoiseChannel::validate`] enforces; the clamp only absorbs rounding at the
/// boundary.
pub(crate) fn thermal_relaxation_rates(t1: f64, t2: f64, gate_time: f64) -> (f64, f64) {
    let gad = 1.0 - (-gate_time / t1).exp();
    let gpd = (1.0 - (gate_time / t1 - 2.0 * gate_time / t2).exp()).clamp(0.0, 1.0);
    (gad, gpd)
}

/// Lower a single-qubit noise channel to its Kraus operators for exact
/// density-matrix evolution. Thermal relaxation (`t2 <= 2*t1`, zero
/// temperature) is the amplitude-damping channel composed with a pure-dephasing
/// channel chosen so populations decay as `exp(-t/t1)` and coherences as
/// `exp(-t/t2)`.
pub(crate) fn kraus_1q(channel: &NoiseChannel) -> Vec<[[Complex64; 2]; 2]> {
    let id = [
        [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
        [Complex64::new(0.0, 0.0), Complex64::new(1.0, 0.0)],
    ];
    let x = [
        [Complex64::new(0.0, 0.0), Complex64::new(1.0, 0.0)],
        [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
    ];
    let y = [
        [Complex64::new(0.0, 0.0), Complex64::new(0.0, -1.0)],
        [Complex64::new(0.0, 1.0), Complex64::new(0.0, 0.0)],
    ];
    let z = [
        [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
        [Complex64::new(0.0, 0.0), Complex64::new(-1.0, 0.0)],
    ];
    let scale = |m: [[Complex64; 2]; 2], s: f64| {
        let s = Complex64::new(s, 0.0);
        [[m[0][0] * s, m[0][1] * s], [m[1][0] * s, m[1][1] * s]]
    };
    let pauli = |px: f64, py: f64, pz: f64| {
        let pi = (1.0 - px - py - pz).max(0.0);
        vec![
            scale(id, pi.sqrt()),
            scale(x, px.sqrt()),
            scale(y, py.sqrt()),
            scale(z, pz.sqrt()),
        ]
    };
    match channel {
        NoiseChannel::Pauli { px, py, pz } => pauli(*px, *py, *pz),
        NoiseChannel::Depolarizing { p } => pauli(p / 3.0, p / 3.0, p / 3.0),
        NoiseChannel::AmplitudeDamping { gamma } => amplitude_damping_kraus(*gamma).to_vec(),
        NoiseChannel::PhaseDamping { gamma } => phase_damping_kraus(*gamma).to_vec(),
        NoiseChannel::ThermalRelaxation {
            t1,
            t2,
            gate_time,
            excited_population,
        } => {
            let (gad, gpd) = thermal_relaxation_rates(*t1, *t2, *gate_time);
            let ad = generalized_amplitude_damping_kraus(gad, *excited_population);
            let pd = phase_damping_kraus(gpd);
            let mut out = Vec::with_capacity(8);
            for a in &ad {
                for p in &pd {
                    let k = crate::gates::mat_mul_2x2(a, p);
                    // A zero operator is a branch that never fires, and at zero
                    // temperature half the set is zero.
                    if k.iter().flatten().any(|e| e.norm_sqr() > 0.0) {
                        out.push(k);
                    }
                }
            }
            out
        }
        NoiseChannel::Custom { kraus } => kraus.clone(),
        NoiseChannel::TwoQubitDepolarizing { .. }
        | NoiseChannel::Kraus2q { .. }
        | NoiseChannel::LeakageTransport { .. } => {
            unreachable!("two-qubit channels have no single-qubit Kraus lowering")
        }
        NoiseChannel::Leakage { .. }
        | NoiseChannel::Seepage { .. }
        | NoiseChannel::QuasiStatic { .. } => {
            unreachable!("evolve_density_matrix rejects channels with per-shot memory")
        }
    }
}

/// Apply one noise event to a density-matrix backend exactly.
pub(crate) fn apply_noise_event_dm(dm: &mut DensityMatrixBackend, event: &NoiseEvent) {
    match &event.channel {
        NoiseChannel::TwoQubitDepolarizing { p } => {
            dm.apply_2q_depolarizing(event.qubits[0], event.qubits[1], *p);
        }
        NoiseChannel::Kraus2q { kraus } => {
            dm.apply_2q_kraus(event.qubits[0], event.qubits[1], kraus);
        }
        other => {
            let kraus = kraus_1q(other);
            dm.apply_1q_kraus(event.qubits[0], &kraus);
        }
    }
}

/// Evolve a density-matrix backend through `circuit` and optional `noise`,
/// starting from `initial_state` when one is given and |0...0> otherwise.
///
/// Measurements are not collapsed; observables and marginals are read off the
/// final mixed state, so the circuit must carry terminal measurements only and
/// no classical conditionals, which
/// [`Circuit::has_terminal_measurements_under_reset_channel`] decides and this
/// function rejects. The circuit is checked against the
/// density-matrix qubit cap and the noise model is validated before any
/// allocation.
pub(crate) fn evolve_density_matrix(
    kind: &BackendKind,
    circuit: &Circuit,
    noise: Option<&NoiseModel>,
    initial_state: Option<super::StartState<'_>>,
    seed: u64,
) -> Result<DensityMatrixBackend> {
    use super::dispatch::{Accel, Family, accel_for, build_density_matrix};

    let accel = accel_for(kind, Family::DensityMatrix, circuit.num_qubits);
    if matches!(accel, Accel::Cpu) {
        crate::backend::check_state_allocation(
            "density_matrix",
            circuit.num_qubits,
            crate::backend::max_density_matrix_qubits(),
            crate::backend::DM_QUBIT_CAP_ENV,
        )?;
    }
    if !circuit.has_terminal_measurements_under_reset_channel() {
        return Err(crate::error::PrismError::IncompatibleBackend {
            backend: "density_matrix".into(),
            reason: "the mixture holds every measurement branch at once, so mid-circuit \
                     measurement and classical conditioning cannot be replayed from it; move the \
                     measurements to the end of the circuit, or run trajectories on a backend \
                     with a per-shot pure state"
                .into(),
        });
    }
    if let Some(noise) = noise {
        noise.validate_for(circuit)?;
        noise.ensure_memoryless("density_matrix")?;
    }

    let mut dm = build_density_matrix(&accel, seed);
    match initial_state {
        Some(state) => {
            crate::sim::check_initial_state_len(state, circuit.num_qubits)?;
            state.load(&mut dm, circuit.num_classical_bits)?;
        }
        None => dm.init(circuit.num_qubits, circuit.num_classical_bits)?,
    }
    let Some(noise) = noise else {
        let expanded = super::expand_for_backend(&dm, circuit);
        let fused = super::fuse_for_backend(&dm, &expanded);
        for inst in &fused.instructions {
            if !matches!(inst, Instruction::Measure { .. }) {
                dm.apply(inst)?;
            }
        }
        return Ok(dm);
    };

    let defer = circuit.num_qubits >= MIN_DEFERRED_QUBITS;
    let mut deferred = DeferredSuperoperators::new(circuit.num_qubits);
    for (inst, events) in circuit.instructions.iter().zip(&noise.after_gate) {
        let mut trailing: &[NoiseEvent] = events;
        let mut folded_pair = None;
        match inst {
            Instruction::Gate { gate, targets } => match *targets.as_slice() {
                [q] if (defer || events.iter().any(|e| e.qubits[..] == [q]))
                    && deferred.defer_gate(gate, q) => {}
                [q0, q1] if defer && fold_2q(&mut dm, &mut deferred, gate, q0, q1, events) => {
                    folded_pair = Some((q0, q1));
                }
                _ => {
                    deferred.flush_qubits(&mut dm, targets);
                    dm.apply(inst)?;
                }
            },
            Instruction::Measure { .. } | Instruction::Barrier { .. } => {}
            _ => {
                deferred.flush(&mut dm);
                dm.apply(inst)?;
            }
        }
        if let Some((q0, q1)) = folded_pair {
            let split = events
                .iter()
                .position(|e| straddles(e, q0, q1))
                .unwrap_or(events.len());
            for event in events[..split]
                .iter()
                .filter(|e| e.qubits.iter().all(|&q| q != q0 && q != q1))
            {
                defer_event(&mut dm, &mut deferred, event);
            }
            trailing = &events[split..];
        }
        for event in trailing {
            defer_event(&mut dm, &mut deferred, event);
        }
        if !defer {
            deferred.flush(&mut dm);
        }
    }
    deferred.flush(&mut dm);
    Ok(dm)
}

/// Whether `event` touches exactly one qubit of `(q0, q1)` along with one
/// outside it, which a sweep over the pair can neither absorb nor commute past.
fn straddles(event: &NoiseEvent, q0: usize, q1: usize) -> bool {
    let inside = event.qubits.iter().filter(|&&q| q == q0 || q == q1).count();
    inside > 0 && inside < event.qubits.len()
}

/// Fold the two-qubit `gate` on `(q0, q1)` and the channels in `events` that
/// act on the pair, up to the first that straddles it, into one sweep. False,
/// with nothing applied, when [`DeferredSuperoperators::apply_2q`] declines.
fn fold_2q(
    dm: &mut DensityMatrixBackend,
    deferred: &mut DeferredSuperoperators,
    gate: &Gate,
    q0: usize,
    q1: usize,
    events: &[NoiseEvent],
) -> bool {
    let channels: Vec<PairChannel> = events
        .iter()
        .take_while(|e| !straddles(e, q0, q1))
        .filter(|e| e.qubits.iter().all(|&q| q == q0 || q == q1))
        .map(|e| match &e.channel {
            NoiseChannel::TwoQubitDepolarizing { p } => PairChannel::depolarizing_2q(*p),
            NoiseChannel::Kraus2q { kraus } => PairChannel::kraus_2q(kraus, e.qubits[0] != q0),
            other => PairChannel::block(e.qubits[0], &kraus_1q(other)),
        })
        .collect();
    deferred.apply_2q(dm, gate, q0, q1, &channels)
}

/// Hold a one-qubit channel on its qubit's pending map, or sweep a two-qubit
/// one after flushing both its qubits.
fn defer_event(
    dm: &mut DensityMatrixBackend,
    deferred: &mut DeferredSuperoperators,
    event: &NoiseEvent,
) {
    if event.channel.num_qubits() == 1 {
        deferred.defer_channel(event.qubits[0], &kraus_1q(&event.channel));
    } else {
        deferred.flush_qubits(dm, &event.qubits);
        apply_noise_event_dm(dm, event);
    }
}

/// Exact output distribution of `circuit` under `noise`, read off the diagonal
/// of the evolved mixture. Measurements are not collapsed during evolution, so
/// the distribution answers for a whole shot only when every measurement is
/// terminal; callers gate on that shape. Events past the last measurement of
/// every qubit they touch are dropped, see [`settle_measured_qubits`].
pub(crate) fn density_matrix_probabilities(
    kind: &BackendKind,
    circuit: &Circuit,
    noise: &NoiseModel,
    initial_state: Option<super::StartState<'_>>,
    seed: u64,
) -> Result<Vec<f64>> {
    let noise = if circuit.has_terminal_measurements_under_reset_channel() {
        settle_measured_qubits(circuit, noise)?
    } else {
        std::borrow::Cow::Borrowed(noise)
    };
    evolve_density_matrix(kind, circuit, Some(&noise), initial_state, seed)?.probabilities()
}

/// `noise` without the events that fire after the last measurement of every
/// qubit they act on. A shot records the bit before such an event, but the
/// mixture is read once at the end, so applying it would rewrite the outcome.
///
/// # Errors
/// `IncompatibleBackend` for an event on a qubit between two of its
/// measurements, or one coupling a measured qubit to an unmeasured one: the
/// record and the state part ways there, and the mixture holds only the state.
fn settle_measured_qubits<'a>(
    circuit: &Circuit,
    noise: &'a NoiseModel,
) -> Result<std::borrow::Cow<'a, NoiseModel>> {
    let mut first = vec![usize::MAX; circuit.num_qubits];
    let mut last = vec![usize::MAX; circuit.num_qubits];
    for (idx, inst) in circuit.instructions.iter().enumerate() {
        if let Instruction::Measure { qubit, .. } = *inst {
            first[qubit] = first[qubit].min(idx);
            last[qubit] = idx;
        }
    }
    let measured = |q: usize, idx: usize| first.get(q).is_some_and(|&m| m <= idx);
    let settled = |event: &NoiseEvent, idx: usize| {
        event
            .qubits
            .iter()
            .all(|&q| last.get(q).is_some_and(|&m| m <= idx))
    };

    let mut any_settled = false;
    for (idx, events) in noise.after_gate.iter().enumerate() {
        for event in events {
            if settled(event, idx) {
                any_settled = true;
            } else if event.qubits.iter().any(|&q| measured(q, idx)) {
                return Err(crate::error::PrismError::IncompatibleBackend {
                    backend: "density_matrix".into(),
                    reason: format!(
                        "the noise event after instruction {idx} acts on a measured qubit that \
                         is measured again or coupled to an unmeasured one, so the recorded bit \
                         and the state diverge and one mixture cannot hold both; run \
                         trajectories on a backend with a per-shot pure state"
                    ),
                });
            }
        }
    }
    if !any_settled {
        return Ok(std::borrow::Cow::Borrowed(noise));
    }
    let after_gate = noise
        .after_gate
        .iter()
        .enumerate()
        .map(|(idx, events)| {
            events
                .iter()
                .filter(|event| !settled(event, idx))
                .cloned()
                .collect()
        })
        .collect();
    Ok(std::borrow::Cow::Owned(NoiseModel {
        after_gate,
        readout: noise.readout.clone(),
    }))
}

/// Exact per-classical-bit measurement marginals under a noise model, evolved
/// on the density-matrix backend. Measurements are read as `Tr(P_1 rho)` on the
/// final state without collapse, so this assumes terminal (non-feedback)
/// measurement, matching the domain of `noisy_marginals_analytical`.
#[cfg(test)]
pub(crate) fn dm_noisy_marginals(
    circuit: &Circuit,
    noise: &NoiseModel,
    seed: u64,
) -> Result<Vec<f64>> {
    let dm = evolve_density_matrix(
        &BackendKind::DensityMatrix,
        circuit,
        Some(noise),
        None,
        seed,
    )?;
    let mut result = vec![0.5f64; circuit.num_classical_bits];
    for inst in &circuit.instructions {
        if let Instruction::Measure {
            qubit,
            classical_bit,
        } = inst
            && *classical_bit < result.len()
        {
            result[*classical_bit] = dm.qubit_probability(*qubit)?;
        }
    }
    Ok(result)
}

/// Exact `Tr(rho P_k)` for each joint Pauli `P_k`, evolving the density-matrix
/// backend through `circuit` and (optionally) `noise`: the exact (zero-variance)
/// analogue of trajectory-averaged expectation values. The circuit must fit the
/// density-matrix qubit cap and contain no measurement, which returns
/// `IncompatibleBackend`; resets are applied as channels. Observables reuse the same
/// Pauli-mask reduction as [`run_expectation_values`](crate::sim::run_expectation_values).
pub fn density_matrix_expectation_values(
    circuit: &Circuit,
    observables: &[Vec<crate::PauliTerm>],
    noise: Option<&NoiseModel>,
    seed: u64,
) -> Result<Vec<f64>> {
    dm_expectation_values(
        &BackendKind::DensityMatrix,
        circuit,
        observables,
        noise,
        None,
        seed,
    )
}

/// [`density_matrix_expectation_values`] with a caller-supplied start state,
/// which `Simulate::expectation_values` carries and the public entry point
/// leaves at |0...0>.
pub(crate) fn dm_expectation_values(
    kind: &BackendKind,
    circuit: &Circuit,
    observables: &[Vec<crate::PauliTerm>],
    noise: Option<&NoiseModel>,
    initial_state: Option<super::StartState<'_>>,
    seed: u64,
) -> Result<Vec<f64>> {
    if circuit
        .instructions
        .iter()
        .any(|inst| matches!(inst, Instruction::Measure { .. }))
    {
        return Err(crate::error::PrismError::IncompatibleBackend {
            backend: "density_matrix".into(),
            reason: "density-matrix expectation values require a circuit without \
                     measurements, since the evolution skips them and a measured qubit would \
                     report its pre-measurement value; remove the measurements"
                .into(),
        });
    }
    let dm = evolve_density_matrix(kind, circuit, noise, initial_state, seed)?;
    let masks = observables
        .iter()
        .map(|obs| crate::sim::pauli_masks(obs, circuit.num_qubits))
        .collect::<Result<Vec<_>>>()?;
    Ok(dm.expectations_pauli(&masks))
}

#[path = "noise_builder.rs"]
pub mod builder;

pub use builder::{DriftDistribution, GateFilter, NoiseBuilder};

#[cfg(test)]
#[path = "noise_tests.rs"]
mod tests;