QEC annotations, fault introspection, ML decoder, multi-decoder sweep, noise model unification

Unified PauliAnnotation system, fault introspection API, lookup table decoder,
per-gate noise model, DEM-level decoder comparison infrastructure, and p_init->p_prep rename.

Key additions:
- PauliAnnotation (Detector/Observable/Operator) replacing three separate types
- Fault introspection: gate_fault_locations(), events(), for_each_fault_combo(w)
- LookupDecoder: ML lookup table from weighted fault enumeration
- Consistent per-gate noise model (p/3 for 1q, p/15 for 2q) across all DEM paths
- DemCheckMatrix in pecos-decoder-core for DEM-to-check-matrix conversion
- DemAwareDecoder: wraps BP+OSD/LSD/UnionFind/RelayBP for DEM-level decoding
- Multi-decoder sweep: --decoder pymatching tesseract bp_osd etc. with side-by-side reports
- Parallel Tesseract decode_batch via rayon thread-local instances (~5x speedup)
- Full (non-decomposed) DEMs for Tesseract/check-matrix decoders
- Surface code circuit builder migrated to typed PauliAnnotation
- TickCircuit idle gate insertion, annotation transfer to/from DagCircuit
- T1/T2 idle noise, PauliWeights with pattern-based matching
- Rename p_init to p_prep throughout (preserve external selene_sim API)
- Fix Y anticommutation bug, prep-gate propagation stop, probability sum validation

Author: ciaranra

Cargo.lock

@@ -20,7 +20,7 @@
 
 use criterion::{BenchmarkId, Criterion, Throughput, measurement::Measurement};
 use pecos_gpu_sims::{GpuInfluenceMapData, GpuInfluenceSampler};
-use pecos_qec::fault_tolerance::noisy_sampler::{NoisySampler, UniformNoiseModel};
+use pecos_qec::fault_tolerance::dem_builder::DemSampler;
 use pecos_qec::fault_tolerance::{DagFaultInfluenceMap, InfluenceBuilder};
 use pecos_quantum::DagCircuit;
 use std::hint::black_box;
@@ -106,7 +106,7 @@ fn build_influence_maps(
     num_data: usize,
 ) -> (DagFaultInfluenceMap, GpuInfluenceMapData) {
     let logical_qubits: Vec<usize> = (0..num_data).collect();
-    let builder = InfluenceBuilder::new(circuit).with_logical_z(logical_qubits);
+    let builder = InfluenceBuilder::new(circuit).with_z(logical_qubits);
     let influence_map = builder.build();
 
     let (
@@ -155,11 +155,11 @@ fn bench_cpu_vs_gpu_surface_codes<M: Measurement>(c: &mut Criterion<M>) {
         group.throughput(Throughput::Elements(u64::from(num_shots)));
 
         // CPU benchmark
-        let noise = UniformNoiseModel::depolarizing(p_error);
-        let mut cpu_sampler = NoisySampler::new(&cpu_map, noise, seed);
+        let probs = vec![p_error; cpu_map.locations.len()];
+        let cpu_sampler = DemSampler::from_influence_map(&cpu_map, &probs);
 
         group.bench_with_input(BenchmarkId::new("CPU", &label), &(), |b, ()| {
-            b.iter(|| black_box(cpu_sampler.sample(num_shots as usize)));
+            b.iter(|| black_box(cpu_sampler.sample_statistics(num_shots as usize, seed)));
         });
 
         // GPU benchmark
@@ -196,11 +196,11 @@ fn bench_gpu_sampler_shot_scaling<M: Measurement>(c: &mut Criterion<M>) {
         group.throughput(Throughput::Elements(u64::from(num_shots)));
 
         // CPU benchmark
-        let noise = UniformNoiseModel::depolarizing(p_error);
-        let mut cpu_sampler = NoisySampler::new(&cpu_map, noise, seed);
+        let probs = vec![p_error; cpu_map.locations.len()];
+        let cpu_sampler = DemSampler::from_influence_map(&cpu_map, &probs);
 
         group.bench_with_input(BenchmarkId::new("CPU", &label), &num_shots, |b, &shots| {
-            b.iter(|| black_box(cpu_sampler.sample(shots as usize)));
+            b.iter(|| black_box(cpu_sampler.sample_statistics(shots as usize, seed)));
         });
 
         // GPU benchmark

crates/benchmarks/benches/modules/gpu_influence_sampler.rs

@@ -108,6 +108,16 @@ pub enum GateType {
     /// Free/deallocate a qubit
     QFree = 136,
     Idle = 200,
+    /// Meta-gate: Pauli operator annotation for fault tracking.
+    ///
+    /// This gate carries a Pauli string but has no effect on quantum state.
+    /// Its position in the circuit determines which faults can flip the operator
+    /// (only faults before this node are relevant). The propagator uses it as a
+    /// backward propagation start point.
+    ///
+    /// The Pauli string is encoded in `params`: each param encodes
+    /// `qubit * 4 + pauli_type` where pauli_type is 1=X, 2=Y, 3=Z.
+    PauliOperatorMeta = 210,
     MeasCrosstalkGlobalPayload = 218,
     MeasCrosstalkLocalPayload = 219,
     /// Custom/unrecognized gate type, with actual name stored in metadata
@@ -164,13 +174,23 @@ impl From<u8> for GateType {
             200 => GateType::Idle,
             218 => GateType::MeasCrosstalkGlobalPayload,
             219 => GateType::MeasCrosstalkLocalPayload,
+            210 => GateType::PauliOperatorMeta,
             255 => GateType::Custom,
             _ => panic!("Invalid gate type ID: {value}"),
         }
     }
 }
 
 impl GateType {
+    /// Returns true if this gate type is a meta-gate (annotation, not physical).
+    ///
+    /// Meta-gates have a position in the DAG but do not affect quantum state
+    /// and should not create fault locations or receive noise.
+    #[must_use]
+    pub const fn is_meta(self) -> bool {
+        matches!(self, GateType::PauliOperatorMeta)
+    }
+
     /// Returns the number of angle parameters this gate type requires
     ///
     /// # Returns
@@ -215,7 +235,8 @@ impl GateType {
             | GateType::PZ
             | GateType::QAlloc
             | GateType::QFree
-            | GateType::Custom => 0,
+            | GateType::Custom
+            | GateType::PauliOperatorMeta => 0,
 
             // Gates with one parameter
             GateType::RX
@@ -277,7 +298,11 @@ impl GateType {
             | GateType::Idle
             | GateType::MeasCrosstalkGlobalPayload
             | GateType::MeasCrosstalkLocalPayload
-            | GateType::Custom => 1,
+            | GateType::Custom
+            // PauliOperatorMeta is variable-arity but returns 1 here because
+            // gate validation checks `is_multiple_of(quantum_arity())` and any
+            // count is a multiple of 1. The actual qubit count is in the gate.
+            | GateType::PauliOperatorMeta => 1,
 
             // Two-qubit gates
             GateType::CX
@@ -405,6 +430,7 @@ impl fmt::Display for GateType {
             GateType::MeasCrosstalkGlobalPayload => write!(f, "MeasCrosstalkGlobalPayload"),
             GateType::MeasCrosstalkLocalPayload => write!(f, "MeasCrosstalkLocalPayload"),
             GateType::Custom => write!(f, "Custom"),
+            GateType::PauliOperatorMeta => write!(f, "PauliOperator"),
         }
     }
 }

crates/pecos-core/src/gate_type.rs

@@ -26,8 +26,6 @@ env_logger.workspace = true
 
 [features]
 default = []
-# Generate bindings at build time (requires cuQuantum headers)
-runtime-bindgen = []
 
 [package.metadata.docs.rs]
 # Don't try to build docs on docs.rs (no CUDA/cuQuantum available)

crates/pecos-cuquantum-sys/Cargo.toml

@@ -200,6 +200,178 @@ pub mod utils {
     }
 }
 
+/// Check matrix representation extracted from a Detector Error Model.
+///
+/// Converts a DEM string into the matrices needed by check-matrix-based
+/// decoders (BP+OSD, UnionFind, RelayBP, etc.):
+///
+/// - **check_matrix** `H[d][m]`: 1 if error mechanism `m` flips detector `d`
+/// - **observable_matrix** `L[o][m]`: 1 if mechanism `m` flips observable `o`
+/// - **error_priors** `p[m]`: probability of mechanism `m`
+///
+/// # Example
+///
+/// ```
+/// use pecos_decoder_core::dem::DemCheckMatrix;
+///
+/// let dem = "error(0.01) D0 D1 L0\nerror(0.02) D1 D2";
+/// let dcm = DemCheckMatrix::from_dem_str(dem).unwrap();
+/// assert_eq!(dcm.num_detectors, 3);
+/// assert_eq!(dcm.num_observables, 1);
+/// assert_eq!(dcm.num_mechanisms, 2);
+/// assert_eq!(dcm.error_priors, vec![0.01, 0.02]);
+/// ```
+#[derive(Debug, Clone)]
+pub struct DemCheckMatrix {
+    /// Check matrix: rows = detectors, columns = error mechanisms.
+    pub check_matrix: ndarray::Array2<u8>,
+    /// Observable matrix: rows = observables, columns = error mechanisms.
+    pub observable_matrix: ndarray::Array2<u8>,
+    /// Error probability per mechanism.
+    pub error_priors: Vec<f64>,
+    /// Number of detectors (rows of check_matrix).
+    pub num_detectors: usize,
+    /// Number of observables (rows of observable_matrix).
+    pub num_observables: usize,
+    /// Number of error mechanisms (columns of both matrices).
+    pub num_mechanisms: usize,
+}
+
+impl DemCheckMatrix {
+    /// Parse a DEM string into check matrix form.
+    ///
+    /// Each `error(p) D_i D_j ... L_k ...` line becomes one column in the
+    /// check matrix (for the D entries) and one column in the observable
+    /// matrix (for the L entries). Decomposed mechanisms (`D0 ^ D1`) are
+    /// combined by XOR.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`DecoderError`] if the DEM string is malformed.
+    pub fn from_dem_str(dem: &str) -> Result<Self, DecoderError> {
+        // First pass: collect mechanisms and find dimensions.
+        let mut mechanisms: Vec<(f64, Vec<u32>, Vec<u32>)> = Vec::new();
+        let mut max_detector: Option<u32> = None;
+        let mut max_observable: Option<u32> = None;
+
+        for line in dem.lines() {
+            let line = line.trim();
+            if line.is_empty() || line.starts_with('#') {
+                continue;
+            }
+            if !line.starts_with("error(") {
+                // Skip non-error lines (detector, logical_observable, etc.)
+                continue;
+            }
+
+            // Parse "error(p) D0 D1 ... L0 ..." or "error(p) D0 ^ D1 ..."
+            let close_paren = line.find(')').ok_or_else(|| {
+                DecoderError::InvalidConfiguration("Missing closing parenthesis in error line".into())
+            })?;
+            let prob_str = &line[6..close_paren];
+            let probability: f64 = prob_str.parse().map_err(|_| {
+                DecoderError::InvalidConfiguration(format!("Invalid probability: {prob_str}"))
+            })?;
+
+            let tokens_str = &line[close_paren + 1..];
+
+            // Handle decomposed mechanisms (with ^) by XOR-ing components.
+            let mut det_set = std::collections::BTreeSet::new();
+            let mut obs_set = std::collections::BTreeSet::new();
+
+            for component in tokens_str.split('^') {
+                for token in component.split_whitespace() {
+                    if let Some(d_str) = token.strip_prefix('D') {
+                        let d: u32 = d_str.parse().map_err(|_| {
+                            DecoderError::InvalidConfiguration(format!("Invalid detector: {token}"))
+                        })?;
+                        // XOR: toggle membership
+                        if !det_set.remove(&d) {
+                            det_set.insert(d);
+                        }
+                        max_detector = Some(max_detector.map_or(d, |m| m.max(d)));
+                    } else if let Some(l_str) = token.strip_prefix('L') {
+                        let l: u32 = l_str.parse().map_err(|_| {
+                            DecoderError::InvalidConfiguration(format!("Invalid observable: {token}"))
+                        })?;
+                        if !obs_set.remove(&l) {
+                            obs_set.insert(l);
+                        }
+                        max_observable = Some(max_observable.map_or(l, |m| m.max(l)));
+                    }
+                }
+            }
+
+            let detectors: Vec<u32> = det_set.into_iter().collect();
+            let observables: Vec<u32> = obs_set.into_iter().collect();
+            mechanisms.push((probability, detectors, observables));
+        }
+
+        let num_detectors = max_detector.map_or(0, |m| m as usize + 1);
+        let num_observables = max_observable.map_or(0, |m| m as usize + 1);
+        let num_mechanisms = mechanisms.len();
+
+        // Build matrices.
+        let mut check_matrix = ndarray::Array2::<u8>::zeros((num_detectors, num_mechanisms));
+        let mut observable_matrix = ndarray::Array2::<u8>::zeros((num_observables, num_mechanisms));
+        let mut error_priors = Vec::with_capacity(num_mechanisms);
+
+        for (col, (prob, detectors, observables)) in mechanisms.iter().enumerate() {
+            error_priors.push(*prob);
+            for &d in detectors {
+                check_matrix[[d as usize, col]] = 1;
+            }
+            for &o in observables {
+                observable_matrix[[o as usize, col]] = 1;
+            }
+        }
+
+        Ok(Self {
+            check_matrix,
+            observable_matrix,
+            error_priors,
+            num_detectors,
+            num_observables,
+            num_mechanisms,
+        })
+    }
+
+    /// Compute the observable prediction from a correction vector.
+    ///
+    /// Given a binary correction vector (one entry per mechanism, from a
+    /// check-matrix decoder), returns the observable mask as
+    /// `observable_matrix @ correction (mod 2)`.
+    #[must_use]
+    pub fn observables_from_correction(&self, correction: &[u8]) -> Vec<u8> {
+        let mut obs = vec![0u8; self.num_observables];
+        for (o, row) in self.observable_matrix.rows().into_iter().enumerate() {
+            let mut sum = 0u8;
+            for (m, &val) in row.iter().enumerate() {
+                if val != 0 && m < correction.len() && correction[m] != 0 {
+                    sum ^= 1;
+                }
+            }
+            obs[o] = sum;
+        }
+        obs
+    }
+
+    /// Pack observable predictions into a bitmask (u64).
+    ///
+    /// Bit `i` is set if observable `i` is predicted to flip.
+    #[must_use]
+    pub fn observables_mask_from_correction(&self, correction: &[u8]) -> u64 {
+        let obs = self.observables_from_correction(correction);
+        let mut mask = 0u64;
+        for (i, &v) in obs.iter().enumerate() {
+            if v != 0 {
+                mask |= 1 << i;
+            }
+        }
+        mask
+    }
+}
+
 /// Information about a detector error model
 #[derive(Debug, Clone, PartialEq)]
 pub struct DemInfo {
@@ -301,6 +473,63 @@ mod tests {
         assert_eq!(observables, 2); // L0 and L1
     }
 
+    #[test]
+    fn test_dem_check_matrix_basic() {
+        let dem = "error(0.01) D0 D1 L0\nerror(0.02) D1 D2\nerror(0.03) D0 D2 L0";
+        let dcm = DemCheckMatrix::from_dem_str(dem).unwrap();
+
+        assert_eq!(dcm.num_detectors, 3);
+        assert_eq!(dcm.num_observables, 1);
+        assert_eq!(dcm.num_mechanisms, 3);
+        assert_eq!(dcm.error_priors, vec![0.01, 0.02, 0.03]);
+
+        // Check matrix: mechanism 0 -> D0,D1; mechanism 1 -> D1,D2; mechanism 2 -> D0,D2
+        assert_eq!(dcm.check_matrix[[0, 0]], 1); // D0, mech 0
+        assert_eq!(dcm.check_matrix[[1, 0]], 1); // D1, mech 0
+        assert_eq!(dcm.check_matrix[[2, 0]], 0); // D2, mech 0
+        assert_eq!(dcm.check_matrix[[0, 1]], 0); // D0, mech 1
+        assert_eq!(dcm.check_matrix[[1, 1]], 1); // D1, mech 1
+        assert_eq!(dcm.check_matrix[[2, 1]], 1); // D2, mech 1
+
+        // Observable matrix: mechanism 0 -> L0; mechanism 1 -> none; mechanism 2 -> L0
+        assert_eq!(dcm.observable_matrix[[0, 0]], 1);
+        assert_eq!(dcm.observable_matrix[[0, 1]], 0);
+        assert_eq!(dcm.observable_matrix[[0, 2]], 1);
+    }
+
+    #[test]
+    fn test_dem_check_matrix_observables_from_correction() {
+        let dem = "error(0.01) D0 L0\nerror(0.01) D1 L1\nerror(0.01) D0 D1 L0 L1";
+        let dcm = DemCheckMatrix::from_dem_str(dem).unwrap();
+
+        // Correction activates mechanism 0 -> L0 flips
+        assert_eq!(dcm.observables_mask_from_correction(&[1, 0, 0]), 0b01);
+        // Correction activates mechanism 2 -> L0 and L1 flip
+        assert_eq!(dcm.observables_mask_from_correction(&[0, 0, 1]), 0b11);
+        // Correction activates mechanisms 0 and 2 -> L0 xor L0 = 0, L1 flips
+        assert_eq!(dcm.observables_mask_from_correction(&[1, 0, 1]), 0b10);
+    }
+
+    #[test]
+    fn test_dem_check_matrix_decomposed() {
+        // Decomposed mechanism: D0 ^ D1 means XOR
+        let dem = "error(0.01) D0 D1 ^ D1 D2";
+        let dcm = DemCheckMatrix::from_dem_str(dem).unwrap();
+
+        // D1 appears in both components -> XOR cancels it
+        assert_eq!(dcm.check_matrix[[0, 0]], 1); // D0
+        assert_eq!(dcm.check_matrix[[1, 0]], 0); // D1 cancels
+        assert_eq!(dcm.check_matrix[[2, 0]], 1); // D2
+    }
+
+    #[test]
+    fn test_dem_check_matrix_empty() {
+        let dem = "";
+        let dcm = DemCheckMatrix::from_dem_str(dem).unwrap();
+        assert_eq!(dcm.num_mechanisms, 0);
+        assert_eq!(dcm.num_detectors, 0);
+    }
+
     #[test]
     fn test_dem_config_builder() {
         let config = DemConfigBuilder::new()

crates/pecos-decoder-core/src/dem.rs

@@ -28,7 +28,7 @@ pub use advanced::{
 pub use config::{
     BatchConfig, ConfigBuilder, DecoderConfig, DecodingMethod, PerformanceConfig, SolverType,
 };
-pub use dem::{DemConfig, DemConfigBuilder, DemDecoder, DemInfo};
+pub use dem::{DemCheckMatrix, DemConfig, DemConfigBuilder, DemDecoder, DemInfo};
 pub use errors::{ConfigError, DecoderError, ErrorConvert, GraphError, MatrixError};
 pub use matrix::{CheckMatrixConfig, CheckMatrixDecoder, SparseCheckMatrix};
 pub use results::{