stabilizer_qec 2.1.0

use std::fmt;
use std::ops::{Add, AddAssign, Mul, MulAssign};
use std::str::FromStr;

use super::{BitTable, Error, Result};
use crate::batch_sampler::{BatchSampler, Counts};
use crate::circuit;
use crate::clifford::ir;
use crate::dem::Dem;
use crate::dem_program::{DemProgram, Stats};
use std::sync::OnceLock;
use crate::m2d::M2d;

/// A stabilizer circuit in Stim's circuit language.
///
/// Parse one from Stim's text (`str::parse` or [`Circuit::parse`]); [`Display`](fmt::Display)
/// prints it back as written. Every Clifford gate of Stim's is read (H, S and CX natively, the
/// rest as their exact decompositions), with resets and measurements in all three bases,
/// inverted targets, Pauli-product measurements and rotations, every noise channel,
/// measurement feedback and sweep bits, `MPAD`, detectors, observables (of records and of Pauli
/// targets), coordinates, `TICK`, `REPEAT` and instruction tags. Counts are taken through
/// `REPEAT` blocks without unrolling them.
///
/// ```
/// use stabilizer_qec::Circuit;
///
/// let c: Circuit = "R 0 1\nH 0\nCX 0 1\nM 0 1\nDETECTOR rec[-1] rec[-2]".parse()?;
/// assert_eq!((c.num_qubits(), c.num_measurements(), c.num_detectors()), (2, 2, 1));
/// assert_eq!(c.to_string().parse::<Circuit>()?, c);
/// # Ok::<(), stabilizer_qec::Error>(())
/// ```
///
/// Circuits are also built in code, as in Stim: [`append`](Circuit::append) an instruction at
/// a time, `+` one circuit after another, and `*` a circuit into a `REPEAT` block. A piece may
/// read measurements made before it (a round's detectors comparing with the last round's);
/// sampling or analysing a circuit that reads before its first measurement is the error.
///
/// ```
/// use stabilizer_qec::{Circuit, Target};
///
/// let mut round = Circuit::new();
/// round.append("CX", &[Target::Qubit(0), Target::Qubit(2), Target::Qubit(1), Target::Qubit(2)], &[])?;
/// round.append("MR", &[Target::Qubit(2)], &[0.01])?;
/// round.append("DETECTOR", &[Target::Rec(1), Target::Rec(2)], &[])?;
/// let first: Circuit = "R 0 1 2\nCX 0 2 1 2\nMR 2".parse()?;
/// let memory = &first + &(&round * 10);
/// assert_eq!((memory.num_measurements(), memory.num_detectors()), (11, 10));
/// assert!(memory.to_string().contains("REPEAT 10 {"));
/// # Ok::<(), stabilizer_qec::Error>(())
/// ```
#[derive(Clone, Debug, PartialEq, Default)]
pub struct Circuit {
    pub(crate) inner: circuit::Circuit,
    counts: Counts,
}

/// A target of an instruction, for [`Circuit::append`].
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
#[non_exhaustive]
pub enum Target {
    /// A qubit (`3`).
    Qubit(u32),
    /// A qubit whose measurement result is inverted (`!3`).
    Inverted(u32),
    /// The measurement record this many measurements back (`rec[-1]` is `Rec(1)`, the latest).
    Rec(u32),
    /// A sweep bit (`sweep[0]`).
    Sweep(u32),
    /// A Pauli on a qubit (`X3`, or `!X3` when `inverted`), for `MPP`, `SPP`, `E` and
    /// `OBSERVABLE_INCLUDE`.
    Pauli {
        /// The Pauli.
        pauli: Pauli,
        /// Its qubit.
        qubit: u32,
        /// Whether it is negated.
        inverted: bool,
    },
    /// Joins the Pauli targets either side into one product (`X0*Z1`).
    Combiner,
}

/// A single-qubit Pauli.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum Pauli {
    /// X.
    X,
    /// Y.
    Y,
    /// Z.
    Z,
}

impl fmt::Display for Target {
    /// The target as Stim writes it.
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match *self {
            Target::Qubit(q) => write!(f, "{q}"),
            Target::Inverted(q) => write!(f, "!{q}"),
            Target::Rec(k) => write!(f, "rec[-{k}]"),
            Target::Sweep(k) => write!(f, "sweep[{k}]"),
            Target::Pauli { pauli, qubit, inverted } => {
                let p = match pauli {
                    Pauli::X => 'X',
                    Pauli::Y => 'Y',
                    Pauli::Z => 'Z',
                };
                write!(f, "{}{p}{qubit}", if inverted { "!" } else { "" })
            }
            Target::Combiner => f.write_str("*"),
        }
    }
}

impl Circuit {
    /// The empty circuit.
    pub fn new() -> Circuit {
        Circuit::default()
    }

    /// Parse Stim's circuit text. As Stim's parser, neighbouring instructions that can be
    /// one are fused into one (`H 0` then `H 1` is `H 0 1`).
    ///
    /// ```
    /// use stabilizer_qec::Circuit;
    ///
    /// let c: Circuit = "H 0\nH 1\nTICK\nH 2".parse()?;
    /// assert_eq!(c.to_string(), "H 0 1\nTICK\nH 2\n");
    /// # Ok::<(), stabilizer_qec::Error>(())
    /// ```
    pub fn parse(text: &str) -> Result<Circuit> {
        let engine = circuit::Circuit::parse(text)?;
        // Errors are the engine's, at their lines as written; the text is rewritten only when
        // something fuses, which the item structure shows.
        match ir::Circuit::parse(text) {
            Ok(fused) if !same_shape(&engine.instrs, &fused.items) => Circuit::from_engine(circuit::Circuit::parse(&fused.exact_text())?),
            _ => Circuit::from_engine(engine),
        }
    }

    /// Parse the text with each line its own instruction, nothing fused: the text of a circuit
    /// written from its items, which edits such as Stim's `pop` leave unfused.
    #[cfg_attr(not(feature = "python"), allow(dead_code))]
    pub(crate) fn parse_unfused(text: &str) -> Result<Circuit> {
        Circuit::from_engine(circuit::Circuit::parse(text)?)
    }

    /// Where the faults of the circuit's error model come from (Stim's
    /// `explain_detector_error_model_errors`): for each fault class (each fault of `filter`
    /// when given), every place in the circuit such a fault arises, or with `reduce` one
    /// representative. Loops are walked in full. The text of each is Stim's.
    ///
    /// ```
    /// use stabilizer_qec::Circuit;
    ///
    /// let c: Circuit = "QUBIT_COORDS(1, 2) 0\nR 0\nX_ERROR(0.1) 0\nM 0\nDETECTOR(3) rec[-1]".parse()?;
    /// let explained = c.explain_errors(None, false)?;
    /// assert_eq!(explained[0].dem_error_terms(), vec!["D0".to_string()]);
    /// assert_eq!(explained[0].terms_with_coords()[0].to_string(), "D0[coords 3]");
    /// let location = &explained[0].circuit_error_locations()[0];
    /// assert_eq!(location.instruction_text(), "X_ERROR(0.1) 0[coords 1,2]");
    /// assert!(location.to_string().starts_with("CircuitErrorLocation {"));
    /// assert!(explained[0].to_string().contains("resolving to X_ERROR(0.1) 0[coords 1,2]"));
    /// # Ok::<(), stabilizer_qec::Error>(())
    /// ```
    pub fn explain_errors(&self, filter: Option<&DetectorErrorModel>, reduce: bool) -> Result<Vec<ExplainedError>> {
        let terms = match filter {
            Some(f) => Some(
                f.flat()?
                    .mechanisms
                    .iter()
                    .map(|m| {
                        let mut t: Vec<u64> = m.detectors.iter().map(|&d| u64::from(d)).collect();
                        t.extend((0..64).filter(|k| m.observables >> k & 1 == 1).map(|k| crate::dem_program::OBS | k));
                        t
                    })
                    .collect::<Vec<_>>(),
            ),
            None => None,
        };
        self.explained(terms.as_deref(), reduce)
    }

    fn explained(&self, terms: Option<&[Vec<u64>]>, reduce: bool) -> Result<Vec<ExplainedError>> {
        let coords = std::sync::Arc::new(crate::explain::Coords::of(&self.inner)?);
        let out = crate::explain::explain(&self.inner, terms, reduce)?;
        Ok(out.into_iter().map(|inner| ExplainedError { inner, coords: coords.clone() }).collect())
    }

    /// The circuit's graph-like distance, as Stim's `shortest_graphlike_error` finds it: the
    /// fewest graph-like pieces of its decomposed error model that together flip an observable
    /// and no detector, each explained (with `canonicalize`, by one location). Pieces of three
    /// or more detectors are skipped with `ignore_ungraphlike`, refused otherwise.
    ///
    /// ```
    /// use stabilizer_qec::{Circuit, GeneratedNoise};
    ///
    /// let noise = GeneratedNoise::new().after_clifford_depolarization(0.001);
    /// let c = Circuit::generated("surface_code:rotated_memory_z", 5, 5, &noise)?;
    /// assert_eq!(c.shortest_graphlike_error(true, false)?.len(), 5);
    /// # Ok::<(), stabilizer_qec::Error>(())
    /// ```
    pub fn shortest_graphlike_error(&self, ignore_ungraphlike: bool, canonicalize: bool) -> Result<Vec<ExplainedError>> {
        let dem = self.detector_error_model(&DemOptions::new().decompose_errors(true).approximate_disjoint_errors(Some(1.0)).ignore_decomposition_failures(true))?;
        let faults = crate::distance::shortest_graphlike(dem.flat()?, ignore_ungraphlike)?;
        self.explained(Some(&crate::distance::as_terms(&faults)), canonicalize)
    }

    /// The circuit's distance through hyperedges too, as Stim's
    /// `search_for_undetectable_logical_errors` finds it: a breadth-first search over sets of
    /// fired detectors in its error model (undecomposed), adding a fault at the set's lowest
    /// detector each step; faults of more than `max_degree` detectors are left out, sets larger
    /// than `max_symptoms` are not explored, nor (with `no_increase`) larger sets at all. Each
    /// fault found is explained (with `canonicalize`, by one location). It is the fewest faults
    /// within those limits: an upper bound on the distance.
    pub fn search_for_undetectable_logical_errors(&self, max_symptoms: usize, max_degree: usize, no_increase: bool, canonicalize: bool) -> Result<Vec<ExplainedError>> {
        let dem = self.detector_error_model(&DemOptions::new().approximate_disjoint_errors(Some(1.0)))?;
        let faults = crate::distance::search_undetectable(dem.flat()?, max_symptoms, max_degree, no_increase)?;
        self.explained(Some(&crate::distance::as_terms(&faults)), canonicalize)
    }

    /// One of Stim's generated memory experiments (`stim.Circuit.generated`), character for
    /// character as Stim writes it: `task` is `"repetition_code:memory"`,
    /// `"surface_code:rotated_memory_x"` or `_z`, `"surface_code:unrotated_memory_x"` or `_z`,
    /// or `"color_code:memory_xyz"`.
    ///
    /// ```
    /// use stabilizer_qec::{Circuit, GeneratedNoise};
    ///
    /// let noise = GeneratedNoise::new().after_clifford_depolarization(0.001);
    /// let c = Circuit::generated("surface_code:rotated_memory_z", 3, 5, &noise)?;
    /// assert_eq!((c.num_qubits(), c.num_detectors()), (26, 40));
    /// # Ok::<(), stabilizer_qec::Error>(())
    /// ```
    pub fn generated(task: &str, distance: u32, rounds: u64, noise: &GeneratedNoise) -> Result<Circuit> {
        Circuit::parse(&crate::generated::generate(task, distance, rounds, &noise.inner)?)
    }

    /// Stim's generated colour code (`"color_code:memory_xyz"`) flattened, each detector
    /// `(x, y, t)` given Chromobius's colour and basis as a 4th coordinate, `(y + t) mod 3`, as
    /// Chromobius's authors annotate it. Chromobius does not decode these circuits under
    /// circuit noise (its tests leave them out: their faults cannot all be split into basic
    /// ones), and neither does [`ColorMatching`](crate::ColorMatching); it does decode the
    /// [`CssCode`](crate::CssCode) colour memories.
    ///
    /// ```
    /// use stabilizer_qec::{Circuit, GeneratedNoise};
    ///
    /// let noise = GeneratedNoise::new().after_clifford_depolarization(0.001);
    /// let c = Circuit::generated_with_colors("color_code:memory_xyz", 3, 3, &noise)?;
    /// assert!(c.to_string().contains("DETECTOR(2, 0, 0, 0)"));
    /// # Ok::<(), stabilizer_qec::Error>(())
    /// ```
    pub fn generated_with_colors(task: &str, distance: u32, rounds: u64, noise: &GeneratedNoise) -> Result<Circuit> {
        Circuit::parse(&crate::generated::generate_annotated(task, distance, rounds, &noise.inner)?)
    }

    pub(crate) fn from_engine(inner: circuit::Circuit) -> Result<Circuit> {
        let counts = Counts::of(&inner.instrs)?;
        fits(&counts)?;
        Ok(Circuit { inner, counts })
    }

    /// Append one instruction: `name` (any of the circuit language's), its targets, and its
    /// arguments (probabilities, coordinates, an observable's index). As Stim's
    /// `Circuit.append`; an instruction the language does not allow is the error, and the
    /// circuit is left as it was.
    pub fn append(&mut self, name: &str, targets: &[Target], args: &[f64]) -> Result<()> {
        self.append_tagged(name, targets, args, "")
    }

    /// [`append`](Circuit::append), with Stim's instruction tag (`H[tag] 0`).
    pub fn append_tagged(&mut self, name: &str, targets: &[Target], args: &[f64], tag: &str) -> Result<()> {
        let targets: Vec<String> = targets.iter().map(Target::to_string).collect();
        let piece = Circuit::from_engine(circuit::Circuit::instruction(name, tag, args, &targets)?)?;
        self.append_circuit(&piece)
    }

    /// Append Stim's circuit text, which may hold several instructions and `REPEAT` blocks.
    pub fn append_text(&mut self, text: &str) -> Result<()> {
        let piece = Circuit::parse(text)?;
        self.append_circuit(&piece)
    }

    /// Append another circuit's instructions.
    /// Its first instruction is fused onto this circuit's last when they can be one, as Stim
    /// fuses them.
    pub fn append_circuit(&mut self, other: &Circuit) -> Result<()> {
        let counts = self.counts.then(other.counts)?;
        fits(&counts)?;
        let mut rest = other.inner.instrs.iter();
        if let (Some(last), Some(first)) = (self.inner.instrs.last(), other.inner.instrs.first()) {
            if let Some(both) = fused(last, first) {
                *self.inner.instrs.last_mut().expect("nonempty") = both;
                rest.next();
            }
        }
        self.inner.instrs.extend(rest.cloned());
        self.counts = counts;
        Ok(())
    }

    /// `REPEAT count { self }`, as Stim's `circuit * count`: the empty circuit for 0, the
    /// circuit itself for 1. The error is a count whose measurements or detectors overflow.
    pub fn repeated(&self, count: u64) -> Result<Circuit> {
        let counts = if count == 0 { Counts::default() } else { self.counts.times(count)? };
        fits(&counts)?;
        Ok(Circuit { inner: self.inner.repeated(count), counts })
    }

    /// One more than the largest qubit index any instruction names.
    pub fn num_qubits(&self) -> usize {
        self.counts.qubits
    }

    /// Measurement records per shot.
    pub fn num_measurements(&self) -> usize {
        self.counts.measurements as usize
    }

    /// Detectors per shot.
    pub fn num_detectors(&self) -> usize {
        self.counts.detectors as usize
    }

    /// Logical observables.
    pub fn num_observables(&self) -> usize {
        self.counts.observables
    }

    /// Sweep bits the circuit reads (one more than the largest `sweep[k]`).
    pub fn num_sweep_bits(&self) -> usize {
        self.counts.sweep_bits
    }

    /// The circuit's detector error model, built by walking it backwards as Stim's error analyzer
    /// does (see [`DemOptions`]). Its loops are folded into `repeat` blocks where they repeat,
    /// as Stim folds them, so a long memory's model costs one period of the loop.
    ///
    /// ```
    /// use stabilizer_qec::{Circuit, DemOptions};
    ///
    /// let c: Circuit = "R 0\nX_ERROR(0.125) 0\nM 0\nDETECTOR rec[-1]".parse()?;
    /// assert_eq!(c.detector_error_model(&DemOptions::new())?.to_string().trim(), "error(0.125) D0");
    /// let long: Circuit = "R 0\nREPEAT 1000000 {\n X_ERROR(0.01) 0\n MR 0\n DETECTOR rec[-1]\n}".parse()?;
    /// let dem = long.detector_error_model(&DemOptions::new())?;
    /// assert_eq!((dem.num_detectors(), dem.num_errors()), (1_000_000, 1_000_000));
    /// assert!(dem.to_string().starts_with("repeat 999999 {"));
    /// # Ok::<(), stabilizer_qec::Error>(())
    /// ```
    pub fn detector_error_model(&self, options: &DemOptions) -> Result<DetectorErrorModel> {
        let o = options;
        let build = if o.deltakit_leakage { crate::dem_build::build_deltakit } else { crate::dem_build::build_full };
        let program = build(&self.inner, o.decompose_errors, o.approximate_disjoint_errors, !o.flatten_loops, o.ignore_decomposition_failures, o.allow_gauge_detectors, o.block_decomposition_from_introducing_remnant_edges)?;
        Ok(DetectorErrorModel::from_program(program))
    }

    /// A sampler of detection events and observable flips, seeded: the same seed gives the same
    /// shots on any machine and any number of threads.
    pub fn detector_sampler(&self, seed: u64) -> Result<DetectorSampler> {
        Ok(DetectorSampler { sampler: BatchSampler::new(&self.inner)?, seed, next: 0 })
    }

    /// A sampler for circuits with coherent errors, written as tagged identities Stim reads as
    /// identities (`I_ERROR[R_Z(theta=0.01)] 0`, `II_ERROR[R_ZZ(theta=0.02)] 0 1`, any Pauli
    /// rotation as `I_ERROR[R_PAULI(theta=…, pauli=…)]`). Its shots are the Pauli-twirled
    /// circuit's, each with a weight that restores the interference the twirl leaves out:
    /// weighted, they are distributed as the coherent circuit's. Seeded like
    /// [`Circuit::detector_sampler`].
    ///
    /// ```
    /// use stabilizer_qec::{Circuit, CoherentOptions};
    ///
    /// // Two rotations either side of a gate that leaves them alone add up: sin²(0.35), not
    /// // the twirl's sin²(0.2) + sin²(0.15) (to first order).
    /// let c = Circuit::parse("RX 0\nR 1\nI_ERROR[R_Z(theta=0.2)] 0\nCZ 0 1\nI_ERROR[R_Z(theta=0.15)] 0\nMX 0\nDETECTOR rec[-1]")?;
    /// let s = c.coherent_sampler(1, CoherentOptions::new())?.sample(200_000, 0);
    /// let (mut hit, mut total) = (0.0, 0.0);
    /// for (k, w) in s.weights.iter().enumerate() {
    ///     total += w;
    ///     if s.detectors.get(k, 0) { hit += w; }
    /// }
    /// assert!((hit / total - 0.35f64.sin().powi(2)).abs() < 0.003);
    /// # Ok::<(), stabilizer_qec::Error>(())
    /// ```
    pub fn coherent_sampler(&self, seed: u64, options: CoherentOptions) -> Result<CoherentSampler> {
        let inner = crate::coherent::sampler::CoherentSampler::new(&self.inner, options.inner)?;
        Ok(CoherentSampler { inner, seed, next: 0 })
    }

    /// An exact sampler: the full state vector (up to 24 qubits in use), running every
    /// instruction, Clifford or not (the coherent rotations, `I[T]`, `I[U3(…)]`, amplitude
    /// damping, leakage), its noise channels drawn as quantum trajectories.
    pub fn exact_sampler(&self, seed: u64) -> Result<ExactSampler> {
        Ok(ExactSampler { program: crate::statevec::Program::new(&self.inner)?, seed, next: 0 })
    }

    /// The joint sampler for Clifford+T circuits with coherent errors: the non-Clifford
    /// gates exact (as [`Circuit::stabilizer_rank_sampler`] runs them), the errors drawn as
    /// their twirl and each shot reweighted by the coherent probability of its class, as
    /// [`Circuit::coherent_sampler`] reweights a Clifford circuit's.
    pub fn joint_sampler(&self, seed: u64, options: CoherentOptions) -> Result<JointSampler> {
        let inner = crate::coherent::joint::JointSampler::new(&self.inner, options.inner)?;
        Ok(JointSampler { inner, seed, next: 0 })
    }

    /// The stabilizer-rank sampler: Clifford+T circuits (`I[T]`, `I[R_Z(theta=…)]`,
    /// `I[U3(…)]`, any Pauli rotation written as a gate) at any qubit count, exactly, its cost
    /// growing with the number of non-Clifford gates rather than of qubits. Pauli noise runs
    /// as the frame sampler's; damping and leakage are refused.
    pub fn stabilizer_rank_sampler(&self, seed: u64) -> Result<StabilizerRankSampler> {
        Ok(StabilizerRankSampler { program: crate::stabrank::Program::new(&self.inner)?, seed, next: 0, peak: 0 })
    }

    /// The circuit with each coherent rotation exp(−iθP) replaced by its Pauli twirl, P with
    /// probability sin²θ: the model Stim and every decoder assume. With `merge`, rotations that
    /// are the same fault in different places become one with their angles added (the
    /// coherence-aware model); merging writes the circuit without loops.
    pub fn twirled(&self, merge: bool) -> Result<Circuit> {
        Circuit::from_engine(crate::coherent::twirled(&self.inner, merge)?)
    }

    /// A picture of the circuit, after Stim's `diagram` (see [`DiagramKind`]): its timeline as
    /// text or SVG, what its detectors compare at a moment, or its matching graph.
    ///
    /// ```
    /// use stabilizer_qec::{Circuit, DiagramKind};
    ///
    /// let c: Circuit = "R 0 1\nH 0\nTICK\nCX 0 1\nTICK\nM 0 1\nDETECTOR rec[-1] rec[-2]".parse()?;
    /// let text = c.diagram(DiagramKind::TimelineText)?;
    /// assert!(text.contains("q0: -R-H-@-M:rec[0]-DETECTOR:D0=rec[1]*rec[0]-"));
    /// assert!(text.contains("q1: -R---X-M:rec[1]-"));
    /// // Just before the measurements, the detector compares Z0 Z1 (Stim's detector slice).
    /// assert_eq!(c.diagram(DiagramKind::DetectorSliceText { tick: 2 })?, "q0: -Z:D0-\n     |\nq1: -Z:D0-");
    /// # Ok::<(), stabilizer_qec::Error>(())
    /// ```
    pub fn diagram(&self, kind: DiagramKind) -> Result<String> {
        self.diagram_with_filter(kind, &[crate::clifford::detslice::CoordFilter::default()])
    }

    /// `diagram`, with Stim's `filter_coords` choosing the detectors and observables of slice
    /// diagrams: each filter a detector or observable, or a coordinate prefix (NaN for any
    /// value) matching detectors. The default (one empty prefix) is every detector and no
    /// observable, as in Stim.
    pub fn diagram_with_filter(&self, kind: DiagramKind, filter: &[crate::clifford::detslice::CoordFilter]) -> Result<String> {
        use crate::diagram as d;
        let ir = || crate::clifford::ir::Circuit::parse(&self.inner.to_stim()).map_err(Error::new);
        let keep = || -> Result<std::collections::BTreeSet<String>> {
            let c = ir()?;
            let n = crate::clifford::transform::count_detectors(&c);
            let coords = crate::clifford::transform::detector_coordinates(&c, &(0..n).collect()).map_err(Error::new)?;
            let mut out = std::collections::BTreeSet::new();
            for (k, v) in &coords {
                if filter.iter().any(|f| f.matches(v, *k)) {
                    out.insert(format!("D{k}"));
                }
            }
            for f in filter {
                if let Some(t) = f.target.filter(|&t| crate::clifford::rev_tracker::is_observable(t)) {
                    out.insert(crate::clifford::rev_tracker::dem_target_str(t));
                }
            }
            Ok(out)
        };
        Ok(match kind {
            DiagramKind::TimelineText => d::timeline_text(&self.inner)?,
            DiagramKind::TimelineSvg => d::timeline_svg(&self.inner)?,
            DiagramKind::DetectorSliceText { tick } => crate::clifford::detslice::SliceSet::from_circuit_ticks(&ir()?, tick, 1, filter).map_err(Error::new)?.text(),
            DiagramKind::DetectorSliceSvg { tick } => d::detslice_svg_filtered(&self.inner, tick, Some(&keep()?))?,
            DiagramKind::DetectorSlicesSvg { ticks, rows } => d::slices_svg_filtered(&self.inner, ticks.0..ticks.1, rows, true, false, Some(&keep()?))?,
            DiagramKind::TimeSliceSvg { ticks, rows } => d::slices_svg(&self.inner, ticks.0..ticks.1, rows, false, true)?,
            DiagramKind::DetectorSliceWithOpsSvg { ticks, rows } => d::slices_svg_filtered(&self.inner, ticks.0..ticks.1, rows, true, true, Some(&keep()?))?,
            DiagramKind::MatchGraphSvg => {
                let dem = self.detector_error_model(&DemOptions::new().decompose_errors(true).approximate_disjoint_errors(Some(1.0)))?;
                dem.matchgraph_svg()?
            }
        })
    }

    /// A converter from raw measurement records (and sweep bits) to detection events and
    /// observable flips, as `stim m2d`. Its reference run keeps a dense tableau: at most 16,384
    /// qubits.
    pub fn measurement_converter(&self) -> Result<MeasurementConverter> {
        Ok(MeasurementConverter { m2d: M2d::new(&self.inner)? })
    }
}

/// Counts this machine can index.
fn fits(c: &Counts) -> Result<()> {
    if usize::try_from(c.measurements).is_err() || usize::try_from(c.detectors).is_err() {
        return Err(Error::new("too many measurements or detectors for this machine"));
    }
    Ok(())
}

/// One circuit after another. Panics where [`Circuit::append_circuit`] would fail: when the
/// measurements or detectors overflow a count.
impl Add<&Circuit> for &Circuit {
    type Output = Circuit;

    fn add(self, rhs: &Circuit) -> Circuit {
        let mut out = self.clone();
        out += rhs;
        out
    }
}

/// One circuit after another (see `&Circuit + &Circuit`).
impl Add<&Circuit> for Circuit {
    type Output = Circuit;

    fn add(mut self, rhs: &Circuit) -> Circuit {
        self += rhs;
        self
    }
}

/// Append a circuit (see `&Circuit + &Circuit`).
impl AddAssign<&Circuit> for Circuit {
    fn add_assign(&mut self, rhs: &Circuit) {
        if let Err(e) = self.append_circuit(rhs) {
            panic!("{e}");
        }
    }
}

/// [`Circuit::repeated`]. Panics where it would fail: when the measurements or detectors
/// overflow a count.
impl Mul<u64> for &Circuit {
    type Output = Circuit;

    fn mul(self, count: u64) -> Circuit {
        self.repeated(count).unwrap_or_else(|e| panic!("{e}"))
    }
}

/// [`Circuit::repeated`] (see `&Circuit * u64`).
impl Mul<u64> for Circuit {
    type Output = Circuit;

    fn mul(self, count: u64) -> Circuit {
        &self * count
    }
}

/// [`Circuit::repeated`], in place (see `&Circuit * u64`).
impl MulAssign<u64> for Circuit {
    fn mul_assign(&mut self, count: u64) {
        *self = &*self * count;
    }
}

impl FromStr for Circuit {
    type Err = Error;

    fn from_str(text: &str) -> Result<Circuit> {
        Circuit::parse(text)
    }
}

/// Whether the engine's instructions are the parsed items one for one, loops alike: that is,
/// whether parsing fused nothing.
fn same_shape(instrs: &[circuit::Instr], items: &[ir::Item]) -> bool {
    instrs.len() == items.len()
        && instrs.iter().zip(items).all(|(i, it)| match (i, it) {
            (circuit::Instr::Repeat { body, .. }, ir::Item::Repeat { body: b, .. }) => same_shape(body, &b.items),
            (circuit::Instr::Repeat { .. }, _) | (_, ir::Item::Repeat { .. }) => false,
            _ => true,
        })
}

/// Two instructions as one, when Stim fuses them (the same gate, arguments and tag, of a kind
/// that fuses): `H 0` and `H 1` are `H 0 1`.
fn fused(a: &circuit::Instr, b: &circuit::Instr) -> Option<circuit::Instr> {
    if matches!(a, circuit::Instr::Repeat { .. }) || matches!(b, circuit::Instr::Repeat { .. }) {
        return None;
    }
    let (la, lb) = (circuit::instr_line(a), circuit::instr_line(b));
    let head = |l: &str| l.split(|c: char| c == ' ' || c == '(' || c == '[').next().unwrap_or("").to_string();
    if head(&la) != head(&lb) {
        return None;
    }
    let one = |l: &str| match ir::Circuit::parse(l).ok()?.items.pop()? {
        ir::Item::Op(op) => Some(op),
        _ => None,
    };
    let (mut x, y) = (one(&la)?, one(&lb)?);
    if !x.fuse(&y) {
        return None;
    }
    let mut c = circuit::Circuit::parse(&x.exact_line()).ok()?;
    (c.instrs.len() == 1).then(|| c.instrs.pop().expect("one"))
}

impl fmt::Display for Circuit {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.write_str(&self.inner.to_stim())
    }
}

/// A kind of [`Circuit::diagram`], after Stim's.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
#[non_exhaustive]
pub enum DiagramKind {
    /// The timeline as text (Stim's `timeline-text`): a row per qubit, each operation in the
    /// first column its wires are free, `TICK` groups bracketed, loops drawn once.
    TimelineText,
    /// The timeline as an SVG picture (`timeline-svg`).
    TimelineSvg,
    /// What each detector compares after `tick` `TICK`s, as Stim's text diagram
    /// (`detslice-text`): each qubit's line with each detector's Pauli on it.
    DetectorSliceText {
        /// The moment: after this many `TICK`s (0 is the start).
        tick: u64,
    },
    /// The same as an SVG picture over the qubits' coordinates (`detslice-svg`).
    DetectorSliceSvg {
        /// The moment.
        tick: u64,
    },
    /// The decomposed model's matching graph as an SVG picture (`matchgraph-svg`).
    MatchGraphSvg,
    /// Detector slices after each tick of `[start, end)`, a panel each, in `rows` rows
    /// (`detslice-svg` over a range).
    DetectorSlicesSvg {
        /// The ticks, `[start, end)`.
        ticks: (u64, u64),
        /// Rows of panels (`None`: about square).
        rows: Option<u32>,
    },
    /// The operations of each tick of `[start, end)` over the qubits' coordinates
    /// (`timeslice-svg`).
    TimeSliceSvg {
        /// The ticks, `[start, end)`.
        ticks: (u64, u64),
        /// Rows of panels (`None`: about square).
        rows: Option<u32>,
    },
    /// Each tick's operations with the detector slice after them (`detslice-with-ops-svg`).
    DetectorSliceWithOpsSvg {
        /// The ticks, `[start, end)`.
        ticks: (u64, u64),
        /// Rows of panels (`None`: about square).
        rows: Option<u32>,
    },
}

/// A fault class of a circuit's error model and the places in the circuit it arises (see
/// [`Circuit::explain_errors`]). Its `Display` is Stim's text.
#[derive(Clone, Debug)]
pub struct ExplainedError {
    inner: crate::explain::Explained,
    coords: std::sync::Arc<crate::explain::Coords>,
}

/// A target of an explained error with its coordinates (empty for none). Its `Display` is
/// Stim's text: `D3[coords 1,2,0]`, or the bare target without coordinates.
#[derive(Clone, Debug, PartialEq)]
#[non_exhaustive]
pub struct TargetWithCoords {
    /// As Stim writes it: `D3`, `L0`, `X5`, `!Z2`, `7`.
    pub target: String,
    /// Its coordinates: a detector's, or a qubit's last `QUBIT_COORDS`.
    pub coords: Vec<f64>,
}

/// One place a fault arises (Stim's `CircuitErrorLocation`). Its `Display` is Stim's text.
#[derive(Clone, Debug, PartialEq)]
#[non_exhaustive]
pub struct CircuitErrorLocation {
    /// `TICK`s before it.
    pub tick_offset: u64,
    /// The Paulis it applies, in target order.
    pub flipped_pauli_product: Vec<TargetWithCoords>,
    /// The measurement it flips: its record index and the observable measured.
    pub flipped_measurement: Option<(u64, Vec<TargetWithCoords>)>,
    /// The instruction's name.
    pub gate: String,
    /// The instruction's tag (empty for none).
    pub tag: String,
    /// The instruction's arguments.
    pub args: Vec<f64>,
    /// Its targets' range in the instruction, `[start, end)`.
    pub target_range: (u32, u32),
    /// Those targets.
    pub targets_in_range: Vec<TargetWithCoords>,
    /// From the outermost block in: (instruction offset, completed iterations, repetitions).
    pub stack_frames: Vec<(u64, u64, u64)>,
    text: String,
    instruction: String,
}

impl CircuitErrorLocation {
    /// The instruction and the targets the fault covers, with coordinates, as Stim writes them
    /// (`DEPOLARIZE2(0.001) 1[coords 1,1] 9[coords 2,2]`).
    pub fn instruction_text(&self) -> &str {
        &self.instruction
    }
}

impl fmt::Display for TargetWithCoords {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        if self.coords.is_empty() {
            f.write_str(&self.target)
        } else {
            write!(f, "{}[coords {}]", self.target, crate::explain::Coords::at(&self.coords))
        }
    }
}

impl fmt::Display for CircuitErrorLocation {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.write_str(&self.text)
    }
}

impl ExplainedError {
    /// What the fault sets off: detectors (`D3`), then observables (`L0`).
    pub fn dem_error_terms(&self) -> Vec<String> {
        self.terms_with_coords().into_iter().map(|t| t.target).collect()
    }

    /// The same, with each detector's coordinates.
    pub fn terms_with_coords(&self) -> Vec<TargetWithCoords> {
        self.inner
            .terms
            .iter()
            .map(|&t| {
                if t & crate::dem_program::OBS != 0 {
                    TargetWithCoords { target: format!("L{}", t & !crate::dem_program::OBS), coords: Vec::new() }
                } else {
                    TargetWithCoords { target: format!("D{t}"), coords: self.coords.of_detector(t) }
                }
            })
            .collect()
    }

    /// Every place it arises (empty when no single fault of the circuit has these symptoms).
    pub fn circuit_error_locations(&self) -> Vec<CircuitErrorLocation> {
        let c = &self.coords;
        let paulis = |ps: &[(u32, u8)]| -> Vec<TargetWithCoords> {
            ps.iter()
                .map(|&(q, p)| TargetWithCoords {
                    target: format!("{}{}{q}", if p & 16 != 0 { "!" } else { "" }, ["I", "X", "Z", "Y"][(p & 3) as usize]),
                    coords: c.of_qubit(q),
                })
                .collect()
        };
        self.inner
            .locations
            .iter()
            .map(|l| CircuitErrorLocation {
                tick_offset: l.tick_offset,
                flipped_pauli_product: paulis(&l.pauli),
                flipped_measurement: l.measurement.as_ref().map(|(i, obs)| (*i, paulis(obs))),
                gate: l.gate.clone(),
                tag: l.tag.clone(),
                args: l.args.clone(),
                target_range: l.range,
                targets_in_range: l
                    .targets
                    .iter()
                    .map(|t| {
                        let q = t.trim_start_matches('!').trim_start_matches(['X', 'Y', 'Z']).parse::<u32>().ok();
                        TargetWithCoords { target: t.clone(), coords: q.map(|q| c.of_qubit(q)).unwrap_or_default() }
                    })
                    .collect(),
                stack_frames: l.stack.clone(),
                text: l.to_stim(c),
                instruction: l.instruction_text(c),
            })
            .collect()
    }
}

impl fmt::Display for ExplainedError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.write_str(&self.inner.to_stim(&self.coords))
    }
}

/// The noise of [`Circuit::generated`]: Stim's four strengths, each zero (left out) unless set.
#[derive(Clone, Copy, Debug, Default, PartialEq)]
pub struct GeneratedNoise {
    inner: crate::generated::Noise,
}

impl GeneratedNoise {
    /// No noise.
    pub fn new() -> GeneratedNoise {
        GeneratedNoise::default()
    }

    /// `DEPOLARIZE1`/`DEPOLARIZE2` after every Clifford gate.
    pub fn after_clifford_depolarization(mut self, p: f64) -> GeneratedNoise {
        self.inner.after_clifford_depolarization = p;
        self
    }

    /// `DEPOLARIZE1` on the data at the start of every round.
    pub fn before_round_data_depolarization(mut self, p: f64) -> GeneratedNoise {
        self.inner.before_round_data_depolarization = p;
        self
    }

    /// A flip before every measurement.
    pub fn before_measure_flip_probability(mut self, p: f64) -> GeneratedNoise {
        self.inner.before_measure_flip_probability = p;
        self
    }

    /// A flip after every reset.
    pub fn after_reset_flip_probability(mut self, p: f64) -> GeneratedNoise {
        self.inner.after_reset_flip_probability = p;
        self
    }
}

/// How [`Circuit::detector_error_model`] builds a model.
///
/// ```
/// use stabilizer_qec::DemOptions;
///
/// let options = DemOptions::new().decompose_errors(true).approximate_disjoint_errors(Some(0.1));
/// # let _ = options;
/// ```
#[derive(Clone, Copy, Debug, PartialEq, Default)]
pub struct DemOptions {
    decompose_errors: bool,
    approximate_disjoint_errors: Option<f64>,
    flatten_loops: bool,
    ignore_decomposition_failures: bool,
    allow_gauge_detectors: bool,
    block_decomposition_from_introducing_remnant_edges: bool,
    deltakit_leakage: bool,
}

impl DemOptions {
    /// The model Deltakit's Stim builds (`deltakit_stim`, Riverlane's fork): with the leakage
    /// its `HERALD_LEAKAGE_EVENT`s read, each herald firing with its qubit depolarized, and with
    /// its partners' depolarized at each `CX`, `CY`, `CZ`, `SQRT_XX` or `SQRT_XX_DAG` after the
    /// leakage came.
    pub fn deltakit_leakage(mut self, yes: bool) -> DemOptions {
        self.deltakit_leakage = yes;
        self
    }

    /// Undecomposed, no approximation.
    pub fn new() -> DemOptions {
        DemOptions::default()
    }

    /// Stim's `allow_gauge_detectors`: a detector made random by a reset or measurement (a
    /// gauge) becomes a 50% error on it, as Stim models it, instead of refusing the circuit.
    /// Observables must still be deterministic.
    pub fn allow_gauge_detectors(mut self, yes: bool) -> DemOptions {
        self.allow_gauge_detectors = yes;
        self
    }

    /// Stim's `block_decomposition_from_introducing_remnant_edges`: with `decompose_errors`,
    /// split a fault only into pieces other faults already make, never adding a new edge.
    pub fn block_decomposition_from_introducing_remnant_edges(mut self, yes: bool) -> DemOptions {
        self.block_decomposition_from_introducing_remnant_edges = yes;
        self
    }

    /// Split each fault into graph-like pieces (at most two detectors each) as Stim splits them,
    /// which matching needs.
    pub fn decompose_errors(mut self, yes: bool) -> DemOptions {
        self.decompose_errors = yes;
        self
    }

    /// Stim's `approximate_disjoint_errors`: `None` refuses channels whose cases are disjoint
    /// rather than independent (`PAULI_CHANNEL_2`, `ELSE_CORRELATED_ERROR`, heralded errors, a
    /// `PAULI_CHANNEL_1` without an independent equivalent); `Some(t)` approximates them case by
    /// case, refusing any with a probability above `t` (`Some(1.0)` is Stim's `True`).
    pub fn approximate_disjoint_errors(mut self, threshold: Option<f64>) -> DemOptions {
        self.approximate_disjoint_errors = threshold;
        self
    }

    /// Stim's `ignore_decomposition_failures`: with `decompose_errors`, keep a fault that cannot
    /// be split into graph-like pieces whole, rather than refuse the model.
    pub fn ignore_decomposition_failures(mut self, yes: bool) -> DemOptions {
        self.ignore_decomposition_failures = yes;
        self
    }

    /// Stim's `flatten_loops`: walk every pass of every loop and write the model without
    /// `repeat` blocks, rather than folding the loops that repeat (the default).
    pub fn flatten_loops(mut self, yes: bool) -> DemOptions {
        self.flatten_loops = yes;
        self
    }
}

/// A detector error model in Stim's format: independent faults, each with a probability, the
/// detectors it flips and the observables it flips. Parse one from Stim's text, or build one
/// with [`Circuit::detector_error_model`].
///
/// It is held as written, `repeat` blocks and all, and counted through its loops without
/// unrolling them; a decoder made from it unrolls it once (up to 2²⁴ faults and declarations
/// and detector 2²⁴ − 1, past which making the decoder is the error), taking its faults in the
/// order they are written, so a model and its text decode alike.
#[derive(Clone, Debug)]
pub struct DetectorErrorModel {
    program: DemProgram,
    stats: OnceLock<Stats>,
    flat: OnceLock<std::result::Result<Dem, String>>,
}

impl DetectorErrorModel {
    pub(crate) fn from_program(program: DemProgram) -> DetectorErrorModel {
        DetectorErrorModel { stats: OnceLock::new(), program, flat: OnceLock::new() }
    }

    fn stats(&self) -> &Stats {
        self.stats.get_or_init(|| self.program.stats())
    }

    /// Appends the instructions of a piece of model text (parsed alone, so a model can be
    /// built instruction by instruction without reading it again each time).
    #[cfg_attr(not(feature = "python"), allow(dead_code))]
    pub(crate) fn append_text(&mut self, text: &str) -> Result<()> {
        let piece = DemProgram::parse(text)?;
        self.program.instrs.extend(piece.instrs);
        self.stats = OnceLock::new();
        self.flat = OnceLock::new();
        Ok(())
    }

    /// The model as written.
    pub(crate) fn program(&self) -> &DemProgram {
        &self.program
    }

    /// The model unrolled, for a decoder.
    pub(crate) fn flat(&self) -> Result<&Dem> {
        self.flat.get_or_init(|| self.program.to_dem()).as_ref().map_err(|e| Error::new(e.clone()))
    }

    /// Parse Stim's detector-error-model text.
    pub fn parse(text: &str) -> Result<DetectorErrorModel> {
        Ok(DetectorErrorModel::from_program(DemProgram::parse(text)?))
    }

    /// Detectors: one more than the largest index named, counted through loops.
    pub fn num_detectors(&self) -> usize {
        usize::try_from(self.stats().num_detectors).unwrap_or(usize::MAX)
    }

    /// Logical observables.
    pub fn num_observables(&self) -> usize {
        self.stats().num_observables
    }

    /// Faults: `error` lines, counted through loops.
    pub fn num_errors(&self) -> usize {
        usize::try_from(self.stats().num_errors).unwrap_or(usize::MAX)
    }


    /// The fewest of the model's graph-like pieces (each fault's `^`-separated pieces) that
    /// together flip an observable and no detector, as a model of those faults, each with
    /// probability 1: Stim's `shortest_graphlike_error`. Its fault count is the graph-like
    /// distance. Pieces of three or more detectors are skipped with `ignore_ungraphlike`,
    /// refused otherwise.
    pub fn shortest_graphlike_error(&self, ignore_ungraphlike: bool) -> Result<DetectorErrorModel> {
        let faults = crate::distance::shortest_graphlike(self.flat()?, ignore_ungraphlike)?;
        DetectorErrorModel::parse(&crate::distance::to_dem_text(&faults))
    }

    /// Stim's `search_for_undetectable_logical_errors` on the model (see
    /// [`Circuit::search_for_undetectable_logical_errors`]): the faults found, each with
    /// probability 1.
    pub fn search_for_undetectable_logical_errors(&self, max_symptoms: usize, max_degree: usize, no_increase: bool) -> Result<DetectorErrorModel> {
        let faults = crate::distance::search_undetectable(self.flat()?, max_symptoms, max_degree, no_increase)?;
        DetectorErrorModel::parse(&crate::distance::to_dem_text(&faults))
    }

    /// A sampler of the model's faults (Stim's `compile_sampler`): each fault fires
    /// independently with its probability; the same seed gives the same shots on any number
    /// of threads.
    ///
    /// ```
    /// use stabilizer_qec::DetectorErrorModel;
    ///
    /// let dem: DetectorErrorModel = "error(0.5) D0 L0\nerror(0.1) D0 D1".parse()?;
    /// let shots = dem.sampler(7)?.sample(1000, 1);
    /// assert_eq!((shots.detectors.num_rows(), shots.errors.num_bits()), (1000, 2));
    /// # Ok::<(), stabilizer_qec::Error>(())
    /// ```
    pub fn sampler(&self, seed: u64) -> Result<DemSampler> {
        Ok(DemSampler { inner: crate::dem_sampler::DemSampler::new(self.flat()?), seed, next: 0 })
    }

    /// The matching graph as an SVG picture: each detector at its coordinates, each graph-like
    /// fault an edge, those flipping an observable heavier. The model must be decomposed, and
    /// unroll.
    pub fn matchgraph_svg(&self) -> Result<String> {
        Ok(crate::diagram::matchgraph_svg(self.flat()?)?)
    }

    /// The model without `repeat` blocks or `shift_detectors`, as Stim's `flattened`: every
    /// detector absolute, every coordinate shifted. The error is a model that unrolls past 2²⁴
    /// faults and declarations.
    ///
    /// ```
    /// use stabilizer_qec::DetectorErrorModel;
    ///
    /// let folded: DetectorErrorModel = "repeat 2 {\n error(0.125) D0\n shift_detectors 1\n}".parse()?;
    /// assert_eq!(folded.flattened()?.to_string(), "error(0.125) D0\nerror(0.125) D1\n");
    /// # Ok::<(), stabilizer_qec::Error>(())
    /// ```
    pub fn flattened(&self) -> Result<DetectorErrorModel> {
        Ok(DetectorErrorModel::from_program(self.program.flattened()?))
    }
}

impl FromStr for DetectorErrorModel {
    type Err = Error;

    fn from_str(text: &str) -> Result<DetectorErrorModel> {
        DetectorErrorModel::parse(text)
    }
}

impl fmt::Display for DetectorErrorModel {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.write_str(&self.program.to_stim())
    }
}

impl DetectorErrorModel {
    /// The model's text with every number exact (Stim's text rounds them to its precision).
    #[cfg_attr(not(feature = "python"), allow(dead_code))]
    pub(crate) fn exact_text(&self) -> String {
        self.program.to_exact_text()
    }
}

/// Shots of detection events and of observable flips, one row per shot.
#[derive(Clone, Debug, PartialEq, Eq)]
#[non_exhaustive]
pub struct Samples {
    /// One bit per detector.
    pub detectors: BitTable,
    /// One bit per observable.
    pub observables: BitTable,
}

/// Detection events and observable flips, 64 shots to a machine word. Made by
/// [`Circuit::detector_sampler`].
///
/// ```
/// use stabilizer_qec::Circuit;
///
/// let c: Circuit = "R 0\nX_ERROR(0.5) 0\nM 0\nDETECTOR rec[-1]".parse()?;
/// let a = c.detector_sampler(7)?.sample(1000, 1);
/// let b = c.detector_sampler(7)?.sample(1000, 4);
/// assert_eq!(a, b); // the same seed, the same shots, on any number of threads
/// # Ok::<(), stabilizer_qec::Error>(())
/// ```
pub struct DetectorSampler {
    sampler: BatchSampler,
    seed: u64,
    next: u64,
}

impl DetectorSampler {
    /// Detectors per shot.
    pub fn num_detectors(&self) -> usize {
        self.sampler.num_detectors
    }

    /// Observables per shot.
    pub fn num_observables(&self) -> usize {
        self.sampler.num_observables
    }

    /// `shots` more shots. `threads = 0` uses every core; the shots do not depend on it. Shots
    /// are drawn 64 at a time from streams seeded by the sampler's seed and the batch, and a
    /// call that ends partway through 64 discards the rest.
    pub fn sample(&mut self, shots: usize, threads: usize) -> Samples {
        let (d, o) = self.sampler.sample_seeded(self.seed, self.next, shots, threads);
        self.next += shots.div_ceil(64) as u64;
        Samples {
            detectors: BitTable::from_packed(shots, self.num_detectors(), d).expect("the sampler writes whole rows"),
            observables: BitTable::from_packed(shots, self.num_observables(), o).expect("the sampler writes whole rows"),
        }
    }
}

/// Options for [`Circuit::coherent_sampler`].
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct CoherentOptions {
    inner: crate::coherent::kernel::CoherentOptions,
}

impl Default for CoherentOptions {
    fn default() -> Self {
        CoherentOptions::new()
    }
}

impl CoherentOptions {
    /// The defaults: products of up to three interference generators, one step past the
    /// generators a shot's faults touch.
    pub fn new() -> CoherentOptions {
        CoherentOptions { inner: crate::coherent::kernel::CoherentOptions::default() }
    }

    /// The most interference generators a shot's class sum multiplies together.
    pub fn order(mut self, order: usize) -> Self {
        self.inner.order = order.max(1);
        self
    }

    /// How many steps of shared locations a shot's clusters reach past its faults' generators.
    pub fn hops(mut self, hops: usize) -> Self {
        self.inner.hops = hops;
        self
    }
}

/// Weighted shots of a circuit with coherent errors. Made by [`Circuit::coherent_sampler`].
pub struct CoherentSampler {
    inner: crate::coherent::sampler::CoherentSampler,
    seed: u64,
    next: u64,
}

/// Weighted shots: detection events, observable flips and each shot's weight. A rate is
/// Σ w·(event) / Σ w.
#[derive(Clone, Debug, PartialEq)]
pub struct WeightedSamples {
    /// One row per shot, one bit per detector.
    pub detectors: BitTable,
    /// One row per shot, one bit per observable.
    pub observables: BitTable,
    /// Each shot's weight.
    pub weights: Vec<f64>,
}

impl CoherentSampler {
    /// Detectors per shot.
    pub fn num_detectors(&self) -> usize {
        self.inner.program.num_detectors()
    }

    /// Observables per shot.
    pub fn num_observables(&self) -> usize {
        self.inner.program.num_observables()
    }

    /// The coherent rotations in the circuit.
    pub fn num_locations(&self) -> usize {
        self.inner.program.locations.len()
    }

    /// The generators of the circuit's interference (sets of rotations no record tells apart):
    /// the local ones and those completing them to a basis of the whole kernel.
    pub fn num_generators(&self) -> usize {
        self.inner.kernel.generators.len()
    }

    /// The generators no local search finds, such as a logical error at one time with the same
    /// logical at another.
    pub fn num_nonlocal_generators(&self) -> usize {
        self.inner.kernel.num_nonlocal
    }

    /// Whether the generators span the whole kernel: false past the completion's limit (about
    /// 50,000 rotations), where interference between distant faults is left out.
    pub fn kernel_complete(&self) -> bool {
        self.inner.kernel.complete
    }

    /// Whether every shot sums the whole kernel (complete, at most 12 generators besides the
    /// groups of equivalent rotations): the weights are then exact, not truncated.
    pub fn exact_weights(&self) -> bool {
        self.inner.kernel.is_whole()
    }

    /// `shots` more weighted shots, each from its own stream; `threads = 0` uses every core and
    /// never changes them.
    pub fn sample(&mut self, shots: usize, threads: usize) -> WeightedSamples {
        let b = self.inner.sample_seeded(self.seed, self.next, shots, threads);
        self.next += shots as u64;
        WeightedSamples {
            detectors: BitTable::from_packed(shots, self.num_detectors(), b.detectors).expect("the sampler writes whole rows"),
            observables: BitTable::from_packed(shots, self.num_observables(), b.observables).expect("the sampler writes whole rows"),
            weights: b.weights,
        }
    }
}

/// Exact shots from the state vector. Made by [`Circuit::exact_sampler`].
pub struct ExactSampler {
    program: crate::statevec::Program,
    seed: u64,
    next: u64,
}

impl ExactSampler {
    /// Detectors per shot.
    pub fn num_detectors(&self) -> usize {
        self.program.num_detectors()
    }

    /// Observables per shot.
    pub fn num_observables(&self) -> usize {
        self.program.num_observables()
    }

    /// The qubits in use, which the state vector holds.
    pub fn num_qubits(&self) -> usize {
        self.program.num_qubits()
    }

    /// `shots` more shots, each from its own stream; `threads = 0` uses every core and never
    /// changes them.
    pub fn sample(&mut self, shots: usize, threads: usize) -> Samples {
        let (d, o) = crate::statevec::sample_seeded(&self.program, self.seed, self.next, shots, threads);
        self.next += shots as u64;
        Samples {
            detectors: BitTable::from_packed(shots, self.num_detectors(), d).expect("the sampler writes whole rows"),
            observables: BitTable::from_packed(shots, self.num_observables(), o).expect("the sampler writes whole rows"),
        }
    }
}

/// Weighted shots of a Clifford+T circuit with coherent errors. Made by
/// [`Circuit::joint_sampler`].
pub struct JointSampler {
    inner: crate::coherent::joint::JointSampler,
    seed: u64,
    next: u64,
}

impl JointSampler {
    /// Detectors per shot.
    pub fn num_detectors(&self) -> usize {
        self.inner.num_detectors()
    }

    /// Observables per shot.
    pub fn num_observables(&self) -> usize {
        self.inner.num_observables()
    }

    /// The coherent errors, drawn as their twirl.
    pub fn num_locations(&self) -> usize {
        self.inner.num_locations()
    }

    /// The non-Clifford gates, run exactly.
    pub fn num_exact(&self) -> usize {
        self.inner.num_exact()
    }

    /// The generators of the interference between error sets, local and non-local.
    pub fn num_generators(&self) -> usize {
        self.inner.num_generators()
    }

    /// The generators no local search finds (a logical error at one time with one at another).
    pub fn num_nonlocal_generators(&self) -> usize {
        self.inner.num_nonlocal_generators()
    }

    /// Whether the generators span the whole projected kernel (false past the completion's
    /// limit, where interference between distant errors is left out).
    pub fn kernel_complete(&self) -> bool {
        self.inner.kernel_complete()
    }

    /// Whether every shot sums the whole projected kernel: the weights are then exact.
    pub fn exact_weights(&self) -> bool {
        self.inner.is_whole()
    }

    /// `shots` more weighted shots, each from its own stream; `threads = 0` uses every core
    /// and never changes them.
    pub fn sample(&mut self, shots: usize, threads: usize) -> WeightedSamples {
        let b = self.inner.sample_seeded(self.seed, self.next, shots, threads);
        self.next += shots as u64;
        WeightedSamples {
            detectors: BitTable::from_packed(shots, self.num_detectors(), b.detectors).expect("the sampler writes whole rows"),
            observables: BitTable::from_packed(shots, self.num_observables(), b.observables).expect("the sampler writes whole rows"),
            weights: b.weights,
        }
    }
}

/// Shots of a Clifford+T circuit from the stabilizer-rank sampler. Made by
/// [`Circuit::stabilizer_rank_sampler`].
pub struct StabilizerRankSampler {
    program: crate::stabrank::Program,
    seed: u64,
    next: u64,
    peak: usize,
}

impl StabilizerRankSampler {
    /// Detectors per shot.
    pub fn num_detectors(&self) -> usize {
        self.program.num_detectors()
    }

    /// Observables per shot.
    pub fn num_observables(&self) -> usize {
        self.program.num_observables()
    }

    /// Measurements per shot.
    pub fn num_measurements(&self) -> usize {
        self.program.num_measurements()
    }

    /// The circuit's qubits.
    pub fn num_qubits(&self) -> usize {
        self.program.num_qubits()
    }

    /// How many of the circuit's gates split the state: its non-Clifford gates.
    pub fn num_non_clifford(&self) -> usize {
        self.program.num_non_clifford()
    }

    /// The most terms any shot so far has held at once (the rank the sampler reached).
    pub fn peak_terms(&self) -> usize {
        self.peak
    }

    /// `shots` more shots' detection events and observable flips, each from its own stream;
    /// `threads = 0` uses every core and never changes them.
    pub fn sample(&mut self, shots: usize, threads: usize) -> Samples {
        let (d, o, peak) = crate::stabrank::sample_seeded(&self.program, self.seed, self.next, shots, threads);
        self.next += shots as u64;
        self.peak = self.peak.max(peak);
        Samples {
            detectors: BitTable::from_packed(shots, self.num_detectors(), d).expect("the sampler writes whole rows"),
            observables: BitTable::from_packed(shots, self.num_observables(), o).expect("the sampler writes whole rows"),
        }
    }

    /// `shots` more shots' measurement records, from the same stream of shots as `sample`.
    pub fn sample_measurements(&mut self, shots: usize, threads: usize) -> BitTable {
        let (m, peak) = crate::stabrank::measurements_seeded(&self.program, self.seed, self.next, shots, threads);
        self.next += shots as u64;
        self.peak = self.peak.max(peak);
        BitTable::from_packed(shots, self.num_measurements(), m).expect("the sampler writes whole rows")
    }

    /// The probability of a measurement record: exact for a circuit without noise channels or
    /// measurement flips; with them, an average over `samples` draws of the noise.
    pub fn probability_of(&self, record: &[bool], samples: usize) -> Result<f64> {
        Ok(self.program.probability_of(record, samples, self.seed)?)
    }
}

/// Shots drawn from a detector error model's faults. Made by [`DetectorErrorModel::sampler`].
pub struct DemSampler {
    inner: crate::dem_sampler::DemSampler,
    seed: u64,
    next: u64,
}

/// A detector error model's shots: detection events, observable flips, and the faults that
/// fired (one bit per fault, in the model's order).
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct DemSamples {
    /// One bit per detector.
    pub detectors: BitTable,
    /// One bit per observable.
    pub observables: BitTable,
    /// One bit per fault of the model.
    pub errors: BitTable,
}

impl DemSampler {
    /// `shots` more shots. `threads = 0` uses every core; the shots do not depend on it.
    pub fn sample(&mut self, shots: usize, threads: usize) -> DemSamples {
        let (d, o, e) = self.inner.sample_seeded(self.seed, self.next, shots, threads, true);
        self.next += shots.div_ceil(64) as u64;
        DemSamples {
            detectors: BitTable::from_packed(shots, self.inner.num_detectors, d).expect("the sampler writes whole rows"),
            observables: BitTable::from_packed(shots, self.inner.num_observables, o).expect("the sampler writes whole rows"),
            errors: BitTable::from_packed(shots, self.inner.num_errors(), e).expect("the sampler writes whole rows"),
        }
    }
}

/// Raw measurement records to detection events and observable flips, as `stim m2d`: each
/// detector and observable compared with a noiseless reference run. Made by
/// [`Circuit::measurement_converter`].
pub struct MeasurementConverter {
    m2d: M2d,
}

impl MeasurementConverter {
    /// Measurement records per shot.
    pub fn num_measurements(&self) -> usize {
        self.m2d.num_measurements
    }

    /// Sweep bits per shot.
    pub fn num_sweep_bits(&self) -> usize {
        self.m2d.num_sweep_bits
    }

    /// Convert shots of measurement records (and, if the circuit reads any, sweep bits).
    pub fn convert(&self, measurements: &BitTable, sweep_bits: Option<&BitTable>) -> Result<Samples> {
        let shots = measurements.num_rows();
        if measurements.num_bits() != self.num_measurements() {
            return Err(Error::new(format!("{} measurement bits per shot, not {}", measurements.num_bits(), self.num_measurements())));
        }
        let zeros;
        let sweeps = match sweep_bits {
            Some(s) => {
                if s.num_rows() != shots || s.num_bits() != self.num_sweep_bits() {
                    return Err(Error::new(format!(
                        "sweep bits are {} shots of {}, not {shots} of {}",
                        s.num_rows(),
                        s.num_bits(),
                        self.num_sweep_bits()
                    )));
                }
                s
            }
            None => {
                zeros = BitTable::zeros(shots, self.num_sweep_bits());
                &zeros
            }
        };
        let (d, o) = self.m2d.convert_b8(measurements.as_bytes(), sweeps.as_bytes(), shots)?;
        Ok(Samples {
            detectors: BitTable::from_packed(shots, self.m2d.num_detectors, d)?,
            observables: BitTable::from_packed(shots, self.m2d.num_observables, o)?,
        })
    }
}