prism-q 0.35.0

//! Error chain complex of a noisy Clifford circuit and the exact analytic noisy
//! marginals read off it.

use crate::circuit::{Circuit, Instruction};
use crate::error::Result;
use crate::sim::compiled::batch_propagate_backward;
use crate::sim::noise::NoiseModel;

/// Dense binary matrix over GF(2) stored as packed u64 words per row.
/// Row-major: row i is stored in words[i * row_words .. (i+1) * row_words].
struct F2DenseMatrix {
    #[cfg(test)]
    num_rows: usize,
    #[cfg(test)]
    num_cols: usize,
    row_words: usize,
    data: Vec<u64>,
}

impl F2DenseMatrix {
    fn new(num_rows: usize, num_cols: usize) -> Self {
        let row_words = num_cols.div_ceil(64);
        Self {
            #[cfg(test)]
            num_rows,
            #[cfg(test)]
            num_cols,
            row_words,
            data: vec![0u64; num_rows * row_words],
        }
    }

    #[inline(always)]
    fn set(&mut self, row: usize, col: usize) {
        self.data[row * self.row_words + col / 64] |= 1u64 << (col % 64);
    }

    #[cfg(test)]
    fn get(&self, row: usize, col: usize) -> bool {
        (self.data[row * self.row_words + col / 64] >> (col % 64)) & 1 != 0
    }

    #[cfg(test)]
    fn row(&self, row: usize) -> &[u64] {
        let start = row * self.row_words;
        &self.data[start..start + self.row_words]
    }

    #[cfg(test)]
    fn xor_row(&mut self, dst_row: usize, src_row: usize) {
        let rw = self.row_words;
        let (dst_start, src_start) = (dst_row * rw, src_row * rw);
        if dst_start < src_start {
            let (left, right) = self.data.split_at_mut(src_start);
            for w in 0..rw {
                left[dst_start + w] ^= right[w];
            }
        } else {
            let (left, right) = self.data.split_at_mut(dst_start);
            for w in 0..rw {
                right[w] ^= left[src_start + w];
            }
        }
    }

    #[cfg(test)]
    fn swap_rows(&mut self, a: usize, b: usize) {
        if a == b {
            return;
        }
        let rw = self.row_words;
        let (a_start, b_start) = (a * rw, b * rw);
        for w in 0..rw {
            self.data.swap(a_start + w, b_start + w);
        }
    }
}

/// Basis of ker(M) over GF(2) for M ∈ F₂^{m×n}, each vector packed into `u64`
/// words of length n, read off a row reduction of `M^T` augmented with `I_n`.
#[cfg(test)]
fn gf2_kernel(matrix: &F2DenseMatrix) -> Vec<Vec<u64>> {
    let m = matrix.num_rows;
    let n = matrix.num_cols;
    let n_words = n.div_ceil(64);

    let aug_cols = m + n;
    let mut aug = F2DenseMatrix::new(n, aug_cols);

    for r in 0..m {
        for c in 0..n {
            if matrix.get(r, c) {
                aug.set(c, r);
            }
        }
    }
    for i in 0..n {
        aug.set(i, m + i);
    }

    let mut pivot_row = 0;
    for col in 0..m {
        let mut found = None;
        for r in pivot_row..n {
            if aug.get(r, col) {
                found = Some(r);
                break;
            }
        }
        let Some(pr) = found else { continue };

        aug.swap_rows(pivot_row, pr);

        for r in 0..n {
            if r != pivot_row && aug.get(r, col) {
                aug.xor_row(r, pivot_row);
            }
        }
        pivot_row += 1;
    }

    let mut kernel = Vec::new();
    let m_words = m.div_ceil(64);
    for r in 0..n {
        let row = aug.row(r);
        let mt_zero = row[..m_words].iter().enumerate().all(|(w, &val)| {
            if w == m_words - 1 && !m.is_multiple_of(64) {
                val & ((1u64 << (m % 64)) - 1) == 0
            } else {
                val == 0
            }
        });
        if mt_zero {
            let mut kv = vec![0u64; n_words];
            for c in 0..n {
                if aug.get(r, m + c) {
                    kv[c / 64] |= 1u64 << (c % 64);
                }
            }
            kernel.push(kv);
        }
    }

    kernel
}

/// Error chain complex for a Clifford circuit with noise.
///
/// Represents the chain complex C₂ →∂₂→ C₁ →∂₁→ C₀ where:
/// - C₀ = F₂^m (measurement/detector space)
/// - C₁ = F₂^p (error location space)
/// - C₂ = F₂^s (stabilizer space)
/// - ∂₁ = E (the m×p error propagation matrix)
pub struct ErrorChainComplex {
    /// E-matrix: m × p binary matrix. E[d][e] = 1 if error e flips measurement d.
    e_matrix: F2DenseMatrix,
    /// Error probabilities: p_total[e] = px + py + pz for error location e.
    error_probs: Vec<f64>,
    num_measurements: usize,
    num_errors: usize,
    /// dim(im(∂₂) ∩ ker(∂₁)): stabilizer generators undetectable by measurements
    boundary_dim: usize,
    /// dim(H₁) = dim(ker(∂₁)/im(∂₂)): independent logical error classes
    homology_dim: usize,
}

impl ErrorChainComplex {
    /// Build the error chain complex from a Clifford circuit and noise model.
    ///
    /// Uses backward Pauli propagation (same as the compiled noisy sampler)
    /// to determine which measurements are sensitive to each error location.
    pub fn build(circuit: &Circuit, noise: &NoiseModel, _seed: u64) -> Result<Self> {
        let m = circuit
            .instructions
            .iter()
            .filter(|i| matches!(i, Instruction::Measure { .. }))
            .count();
        if m == 0 {
            return Ok(Self {
                e_matrix: F2DenseMatrix::new(0, 0),
                error_probs: Vec::new(),
                num_measurements: 0,
                num_errors: 0,
                boundary_dim: circuit.num_qubits,
                homology_dim: 0,
            });
        }

        let m_words = m.div_ceil(64);
        let n = circuit.num_qubits;

        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![0u64; m_words];

        let mut meas_idx = m;
        for instr in circuit.instructions.iter().rev() {
            if let Instruction::Measure { qubit, .. } = instr {
                meas_idx -= 1;
                let word = meas_idx / 64;
                let bit = meas_idx % 64;
                z_packed[*qubit][word] |= 1u64 << bit;
            }
        }

        let mut error_probs = Vec::new();
        let mut e_cols: Vec<Vec<u64>> = Vec::new();

        for (instr_idx, instr) in circuit.instructions.iter().enumerate().rev() {
            match instr {
                Instruction::Gate { gate, targets } => {
                    let noise_events = &noise.after_gate[instr_idx];
                    for event in noise_events {
                        // An inert channel contributes no column, whatever its
                        // variant.
                        if event.channel.is_inert() {
                            continue;
                        }
                        let (px, py, pz) = event.pauli_probs();
                        let q = event.qubit();
                        let p_total = px + py + pz;
                        if p_total < 1e-15 {
                            continue;
                        }

                        let x_sens = &z_packed[q];
                        let z_sens = &x_packed[q];

                        if px > 1e-15 && x_sens.iter().any(|&w| w != 0) {
                            error_probs.push(px);
                            e_cols.push(x_sens.clone());
                        }

                        if pz > 1e-15 && z_sens.iter().any(|&w| w != 0) {
                            error_probs.push(pz);
                            e_cols.push(z_sens.clone());
                        }

                        if py > 1e-15 {
                            let mut y_sens = vec![0u64; m_words];
                            for w in 0..m_words {
                                y_sens[w] = x_sens[w] ^ z_sens[w];
                            }
                            if y_sens.iter().any(|&w| w != 0) {
                                error_probs.push(py);
                                e_cols.push(y_sens);
                            }
                        }
                    }

                    batch_propagate_backward(
                        &mut x_packed,
                        &mut z_packed,
                        &mut sign_packed,
                        gate,
                        targets.as_slice(),
                        m_words,
                    );
                }
                Instruction::Measure { .. }
                | Instruction::Reset { .. }
                | Instruction::Barrier { .. }
                | Instruction::Save { .. } => {}
                Instruction::Conditional { gate, targets, .. } => {
                    batch_propagate_backward(
                        &mut x_packed,
                        &mut z_packed,
                        &mut sign_packed,
                        gate,
                        targets.as_slice(),
                        m_words,
                    );
                }
                Instruction::Region(_) => {
                    return Err(crate::error::PrismError::IncompatibleBackend {
                        backend: "HomologicalDetectorModel".to_string(),
                        reason: "error propagation does not support guarded regions".to_string(),
                    });
                }
            }
        }

        let p = error_probs.len();
        let mut e_matrix = F2DenseMatrix::new(m, p);

        for (col, col_data) in e_cols.iter().enumerate() {
            for (w, &word) in col_data.iter().enumerate() {
                if word == 0 {
                    continue;
                }
                let base = w * 64;
                let mut bits = word;
                while bits != 0 {
                    let bit = bits.trailing_zeros() as usize;
                    let row = base + bit;
                    if row < m {
                        e_matrix.set(row, col);
                    }
                    bits &= bits - 1;
                }
            }
        }

        let (boundary_dim, homology_dim) = Self::compute_boundary_space(circuit, n);

        Ok(Self {
            e_matrix,
            error_probs,
            num_measurements: m,
            num_errors: p,
            boundary_dim,
            homology_dim,
        })
    }

    /// Forward-propagate stabilizer generators and compute ∂₂ boundary space.
    ///
    /// Returns (boundary_dim, homology_dim) where:
    /// - boundary_dim = dim(im(∂₂) ∩ ker(∂₁)) = stabilizers undetectable by measurements
    /// - homology_dim = dim(H₁) = independent logical error classes
    ///
    /// Algorithm: forward-propagate Z_0,...,Z_{n-1} through the circuit to get
    /// output stabilizer generators. Build X-projection onto measured qubits.
    /// rank(X_proj) counts stabilizers with detectable X components;
    /// H₁ = ker(σ) / (S ∩ ker(σ)) has dim = n - num_measured + rank(X_proj).
    fn compute_boundary_space(circuit: &Circuit, n: usize) -> (usize, usize) {
        if n == 0 {
            return (0, 0);
        }

        let n_words = n.div_ceil(64);
        let mut stab_x: Vec<Vec<u64>> = vec![vec![0u64; n_words]; n];
        let mut stab_z: Vec<Vec<u64>> = vec![vec![0u64; n_words]; n];
        let mut stab_sign = vec![0u64; n_words];

        for i in 0..n {
            stab_z[i][i / 64] |= 1u64 << (i % 64);
        }

        for instr in circuit.instructions.iter() {
            match instr {
                Instruction::Gate { gate, targets } => {
                    batch_propagate_backward(
                        &mut stab_x,
                        &mut stab_z,
                        &mut stab_sign,
                        gate,
                        targets.as_slice(),
                        n_words,
                    );
                }
                Instruction::Conditional { gate, targets, .. } => {
                    batch_propagate_backward(
                        &mut stab_x,
                        &mut stab_z,
                        &mut stab_sign,
                        gate,
                        targets.as_slice(),
                        n_words,
                    );
                }
                _ => {}
            }
        }

        let mut measured = vec![false; n];
        for instr in &circuit.instructions {
            if let Instruction::Measure { qubit, .. } = instr {
                measured[*qubit] = true;
            }
        }
        let num_measured = measured.iter().filter(|&&b| b).count();
        let measured_indices: Vec<usize> = (0..n).filter(|&q| measured[q]).collect();

        if num_measured == 0 {
            return (n, 0);
        }

        let proj_words = num_measured.div_ceil(64);
        let mut proj = vec![0u64; n * proj_words];

        for stab_idx in 0..n {
            for (proj_col, &q) in measured_indices.iter().enumerate() {
                let x_bit = (stab_x[q][stab_idx / 64] >> (stab_idx % 64)) & 1;
                if x_bit != 0 {
                    proj[stab_idx * proj_words + proj_col / 64] |= 1u64 << (proj_col % 64);
                }
            }
        }

        let mut rank = 0;
        let mut pivot_row = 0;
        for col in 0..num_measured {
            let mut found = None;
            for r in pivot_row..n {
                if (proj[r * proj_words + col / 64] >> (col % 64)) & 1 != 0 {
                    found = Some(r);
                    break;
                }
            }
            let Some(pr) = found else { continue };

            if pr != pivot_row {
                for w in 0..proj_words {
                    proj.swap(pivot_row * proj_words + w, pr * proj_words + w);
                }
            }

            for r in 0..n {
                if r != pivot_row && (proj[r * proj_words + col / 64] >> (col % 64)) & 1 != 0 {
                    for w in 0..proj_words {
                        proj[r * proj_words + w] ^= proj[pivot_row * proj_words + w];
                    }
                }
            }

            pivot_row += 1;
            rank += 1;
        }

        let boundary_dim = n - rank;
        let homology_dim = n - num_measured + rank;
        (boundary_dim, homology_dim)
    }

    /// dim(im(∂₂) ∩ ker(∂₁)): stabilizer generators undetectable by the
    /// measurements.
    pub fn boundary_dim(&self) -> usize {
        self.boundary_dim
    }

    /// dim(H₁): independent logical error classes.
    pub fn homology_dim(&self) -> usize {
        self.homology_dim
    }

    /// Compute exact noisy marginals analytically. No sampling, no rank limit.
    ///
    /// For each measurement j, the noisy probability is:
    ///   p_j^noisy = p_j + (1 - 2·p_j) · (1 - f_j) / 2
    /// where f_j = Π_{e: E(j,e)=1} (1 - 2·p_e) is the flip attenuation factor
    /// and p_j is the noiseless marginal (0, 0.5, or 1).
    ///
    /// Cost: O(nnz(E)). Works for any qubit count.
    pub fn noisy_marginals(&self, noiseless_marginals: &[f64]) -> Vec<f64> {
        let m = self.num_measurements;
        let p = self.num_errors;
        if m == 0 || p == 0 {
            return noiseless_marginals.to_vec();
        }

        let mut flip_factor = vec![1.0f64; m];
        let rw = self.e_matrix.row_words;

        for e in 0..p {
            let factor = 1.0 - 2.0 * self.error_probs[e];
            if (factor - 1.0).abs() < 1e-15 {
                continue;
            }

            let col_word = e / 64;
            let col_bit = 1u64 << (e % 64);

            for (j, ff) in flip_factor.iter_mut().enumerate() {
                if self.e_matrix.data[j * rw + col_word] & col_bit != 0 {
                    *ff *= factor;
                }
            }
        }

        let mut result = Vec::with_capacity(m);
        for j in 0..m {
            let p_j = noiseless_marginals[j];
            let p_flip = (1.0 - flip_factor[j]) / 2.0;
            result.push(p_j + (1.0 - 2.0 * p_j) * p_flip);
        }
        result
    }
}

/// Compute exact noisy marginals analytically. No sampling, no rank limit.
///
/// Builds the error chain complex and compiled sampler, then computes
/// exact per-measurement noisy probabilities in O(nnz(E)) time.
/// Works for any qubit count, not limited by syndrome rank.
///
/// Readout error needs no syndrome class on a marginal: a bit reported 1 was
/// either measured 1 and not flipped or measured 0 and flipped, so the reported
/// probability is `p(1 - p10) + (1 - p) p01`. Per-shot sampling is still the
/// only route for correlations between bits.
pub fn noisy_marginals_analytical(
    circuit: &Circuit,
    noise: &NoiseModel,
    seed: u64,
) -> Result<Vec<f64>> {
    noise.ensure_chain_complex_channels()?;
    let ecc = ErrorChainComplex::build(circuit, noise, seed)?;
    let compiled = crate::sim::compiled::compile_measurements(circuit, seed)?;
    let noiseless = compiled.marginal_probabilities();
    let noisy = ecc.noisy_marginals(&noiseless);

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

    let mut result = vec![0.5f64; num_classical];
    for (mi, &cbit) in classical_bit_order.iter().enumerate() {
        if cbit < num_classical && mi < noisy.len() {
            let p = noisy[mi];
            result[cbit] = match noise.readout.get(cbit).and_then(Option::as_ref) {
                Some(err) => p * (1.0 - err.p10) + (1.0 - p) * err.p01,
                None => p,
            };
        }
    }
    Ok(result)
}

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