quantrs2-anneal 0.2.1

//! Error Syndrome Detection and Correction
//!
//! This module implements quantum error syndrome detection and correction algorithms
//! specifically designed for quantum annealing systems. It provides functionality for:
//! - Real-time syndrome detection during annealing
//! - Error classification and pattern recognition
//! - Correction protocol execution
//! - Integration with annealing schedules

use scirs2_core::ndarray::{Array1, Array2};
use scirs2_core::Complex64;
use std::collections::HashMap;
use std::time::{Duration, Instant};

use super::codes::{CodeParameters, ErrorCorrectionCode};
use super::config::{QECResult, QuantumErrorCorrectionError};
use super::logical_operations::LogicalOperation;
use crate::ising::IsingModel;
use crate::qaoa::QuantumState;
use crate::simulator::AnnealingResult;

/// Syndrome detection and correction engine
#[derive(Debug, Clone)]
pub struct SyndromeDetector {
    /// Error correction code being used
    pub code: ErrorCorrectionCode,
    /// Code parameters
    pub parameters: CodeParameters,
    /// Stabilizer generators for syndrome detection
    pub stabilizer_generators: Array2<u8>,
    /// Parity check matrix
    pub parity_check_matrix: Array2<u8>,
    /// Correction lookup table
    pub correction_table: HashMap<Vec<u8>, CorrectionOperation>,
    /// Detection statistics
    pub detection_stats: DetectionStatistics,
    /// Configuration parameters
    pub config: SyndromeDetectorConfig,
}

/// Configuration for syndrome detection
#[derive(Debug, Clone)]
pub struct SyndromeDetectorConfig {
    /// Syndrome measurement frequency (Hz)
    pub measurement_frequency: f64,
    /// Error probability threshold for correction
    pub correction_threshold: f64,
    /// Maximum syndrome history to keep
    pub max_history_length: usize,
    /// Enable adaptive correction timing
    pub adaptive_timing: bool,
    /// Noise model parameters
    pub noise_model: NoiseModel,
    /// Decoding algorithm
    pub decoder: DecodingAlgorithm,
}

/// Noise model for syndrome detection
#[derive(Debug, Clone)]
pub struct NoiseModel {
    /// Single-qubit depolarizing error rate
    pub single_qubit_error_rate: f64,
    /// Two-qubit gate error rate
    pub two_qubit_error_rate: f64,
    /// Measurement error rate
    pub measurement_error_rate: f64,
    /// Coherence time (microseconds)
    pub coherence_time: f64,
    /// Temperature effects
    pub thermal_noise_rate: f64,
}

/// Decoding algorithms for error correction
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum DecodingAlgorithm {
    /// Minimum weight perfect matching
    MinimumWeight,
    /// Maximum likelihood decoding
    MaximumLikelihood,
    /// Neural network decoder
    NeuralNetwork,
    /// Belief propagation
    BeliefPropagation,
    /// Lookup table decoder
    LookupTable,
    /// Union-Find decoder
    UnionFind,
}

/// Correction operation to apply
#[derive(Debug, Clone)]
pub struct CorrectionOperation {
    /// Pauli corrections to apply
    pub pauli_corrections: Vec<PauliCorrection>,
    /// Confidence level of correction
    pub confidence: f64,
    /// Estimated success probability
    pub success_probability: f64,
    /// Required resources
    pub required_resources: CorrectionResources,
}

/// Individual Pauli correction
#[derive(Debug, Clone)]
pub struct PauliCorrection {
    /// Qubit index
    pub qubit: usize,
    /// Pauli operation (I, X, Y, Z)
    pub operation: PauliType,
    /// Correction weight
    pub weight: f64,
}

/// Pauli operation types
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum PauliType {
    I, // Identity
    X, // Pauli-X
    Y, // Pauli-Y
    Z, // Pauli-Z
}

/// Resources required for correction
#[derive(Debug, Clone)]
pub struct CorrectionResources {
    /// Number of syndrome measurements
    pub num_measurements: usize,
    /// Number of correction gates
    pub num_correction_gates: usize,
    /// Estimated time overhead (microseconds)
    pub time_overhead: f64,
    /// Ancilla qubits required
    pub ancilla_qubits: usize,
}

/// Detection and correction statistics
#[derive(Debug, Clone)]
pub struct DetectionStatistics {
    /// Total number of syndrome measurements
    pub total_measurements: usize,
    /// Number of non-trivial syndromes detected
    pub errors_detected: usize,
    /// Number of successful corrections
    pub successful_corrections: usize,
    /// Number of failed corrections
    pub failed_corrections: usize,
    /// Average syndrome detection time
    pub avg_detection_time: Duration,
    /// Average correction time
    pub avg_correction_time: Duration,
    /// Syndrome history
    pub syndrome_history: Vec<SyndromeRecord>,
}

/// Individual syndrome measurement record
#[derive(Debug, Clone)]
pub struct SyndromeRecord {
    /// Timestamp of measurement
    pub timestamp: Instant,
    /// Measured syndrome
    pub syndrome: Vec<u8>,
    /// Error locations (if decoded)
    pub error_locations: Option<Vec<usize>>,
    /// Correction applied
    pub correction_applied: Option<CorrectionOperation>,
    /// Measurement confidence
    pub confidence: f64,
}

/// Result of syndrome detection and correction
#[derive(Debug, Clone)]
pub struct SyndromeResult {
    /// Detected syndrome
    pub syndrome: Vec<u8>,
    /// Recommended correction
    pub correction: Option<CorrectionOperation>,
    /// Detection confidence
    pub confidence: f64,
    /// Error locations
    pub error_locations: Vec<usize>,
    /// Detection metadata
    pub metadata: SyndromeMetadata,
}

/// Metadata for syndrome detection
#[derive(Debug, Clone)]
pub struct SyndromeMetadata {
    /// Time taken for detection
    pub detection_time: Duration,
    /// Number of measurements used
    pub num_measurements: usize,
    /// Decoder used
    pub decoder_used: DecodingAlgorithm,
    /// Noise level estimate
    pub estimated_noise_level: f64,
}

impl SyndromeDetector {
    /// Create new syndrome detector
    pub fn new(
        code: ErrorCorrectionCode,
        parameters: CodeParameters,
        config: SyndromeDetectorConfig,
    ) -> QECResult<Self> {
        let stabilizer_generators = Self::generate_stabilizers(&code, &parameters)?;
        let parity_check_matrix = Self::generate_parity_check_matrix(&code, &parameters)?;
        let correction_table = Self::build_correction_table(&code, &parameters, &config.decoder)?;

        Ok(Self {
            code,
            parameters,
            stabilizer_generators,
            parity_check_matrix,
            correction_table,
            detection_stats: DetectionStatistics::new(),
            config,
        })
    }

    /// Detect syndrome in quantum state
    pub fn detect_syndrome(&mut self, state: &QuantumState) -> QECResult<SyndromeResult> {
        let start_time = Instant::now();

        // Perform syndrome measurements
        let syndrome = self.measure_syndrome(state)?;

        // Decode errors if syndrome is non-trivial
        let (correction, error_locations) = if self.is_trivial_syndrome(&syndrome) {
            (None, Vec::new())
        } else {
            let correction = self.decode_syndrome(&syndrome)?;
            let error_locations = self.extract_error_locations(&correction);
            (Some(correction), error_locations)
        };

        // Calculate confidence based on syndrome stability and noise model
        let confidence = self.calculate_detection_confidence(&syndrome, &error_locations)?;

        // Update statistics
        self.update_detection_stats(
            &syndrome,
            correction.as_ref(),
            confidence,
            start_time.elapsed(),
        );

        let metadata = SyndromeMetadata {
            detection_time: start_time.elapsed(),
            num_measurements: self.stabilizer_generators.nrows(),
            decoder_used: self.config.decoder.clone(),
            estimated_noise_level: self.estimate_noise_level(&syndrome),
        };

        Ok(SyndromeResult {
            syndrome,
            correction,
            confidence,
            error_locations,
            metadata,
        })
    }

    /// Apply correction to quantum state
    pub fn apply_correction(
        &mut self,
        state: &mut QuantumState,
        correction: &CorrectionOperation,
    ) -> QECResult<()> {
        let start_time = Instant::now();

        for pauli_correction in &correction.pauli_corrections {
            self.apply_pauli_correction(state, pauli_correction)?;
        }

        // Update correction statistics
        self.detection_stats.successful_corrections += 1;
        self.detection_stats.avg_correction_time = Self::update_average_time(
            self.detection_stats.avg_correction_time,
            start_time.elapsed(),
            self.detection_stats.successful_corrections,
        );

        Ok(())
    }

    /// Perform full detection and correction cycle
    pub fn detect_and_correct(&mut self, state: &mut QuantumState) -> QECResult<SyndromeResult> {
        let syndrome_result = self.detect_syndrome(state)?;

        if let Some(ref correction) = syndrome_result.correction {
            if syndrome_result.confidence >= self.config.correction_threshold {
                self.apply_correction(state, correction)?;
            }
        }

        Ok(syndrome_result)
    }

    /// Generate stabilizer generators for the error correction code
    fn generate_stabilizers(
        code: &ErrorCorrectionCode,
        parameters: &CodeParameters,
    ) -> QECResult<Array2<u8>> {
        match code {
            ErrorCorrectionCode::SurfaceCode => Self::generate_surface_code_stabilizers(parameters),
            ErrorCorrectionCode::RepetitionCode => {
                Self::generate_repetition_code_stabilizers(parameters)
            }
            ErrorCorrectionCode::SteaneCode => Self::generate_steane_code_stabilizers(parameters),
            ErrorCorrectionCode::ShorCode => Self::generate_shor_code_stabilizers(parameters),
            _ => Err(QuantumErrorCorrectionError::CodeError(format!(
                "Stabilizer generation not implemented for code: {code:?}"
            ))),
        }
    }

    /// Generate surface code stabilizers
    fn generate_surface_code_stabilizers(parameters: &CodeParameters) -> QECResult<Array2<u8>> {
        let d = parameters.distance;
        let num_qubits = parameters.num_physical_qubits;
        let num_stabilizers = num_qubits - parameters.num_logical_qubits;

        let mut stabilizers = Array2::zeros((num_stabilizers, 2 * num_qubits));

        // Generate X-type and Z-type stabilizers for surface code
        let mut stabilizer_idx = 0;

        // X-type stabilizers (plaquette checks)
        for row in 0..(d - 1) {
            for col in 0..(d - 1) {
                if (row + col) % 2 == 0 {
                    // X-type plaquettes
                    let qubits = Self::get_plaquette_qubits(row, col, d);
                    for &qubit in &qubits {
                        if qubit < num_qubits {
                            stabilizers[[stabilizer_idx, qubit]] = 1; // X part
                        }
                    }
                    stabilizer_idx += 1;
                }
            }
        }

        // Z-type stabilizers (vertex checks)
        for row in 0..d {
            for col in 0..d {
                if (row + col) % 2 == 1 {
                    // Z-type vertices
                    let qubits = Self::get_vertex_qubits(row, col, d);
                    for &qubit in &qubits {
                        if qubit < num_qubits {
                            stabilizers[[stabilizer_idx, num_qubits + qubit]] = 1;
                            // Z part
                        }
                    }
                    stabilizer_idx += 1;
                    if stabilizer_idx >= num_stabilizers {
                        break;
                    }
                }
            }
            if stabilizer_idx >= num_stabilizers {
                break;
            }
        }

        Ok(stabilizers)
    }

    /// Generate repetition code stabilizers
    fn generate_repetition_code_stabilizers(parameters: &CodeParameters) -> QECResult<Array2<u8>> {
        let n = parameters.num_physical_qubits;
        let num_stabilizers = n - 1;

        let mut stabilizers = Array2::zeros((num_stabilizers, 2 * n));

        // Generate ZZ stabilizers for repetition code
        for i in 0..num_stabilizers {
            stabilizers[[i, n + i]] = 1; // Z on qubit i
            stabilizers[[i, n + i + 1]] = 1; // Z on qubit i+1
        }

        Ok(stabilizers)
    }

    /// Generate Steane code stabilizers
    fn generate_steane_code_stabilizers(parameters: &CodeParameters) -> QECResult<Array2<u8>> {
        let n = 7; // Steane code is always 7 qubits
        let mut stabilizers = Array2::zeros((6, 14)); // 6 stabilizers for [[7,1,3]] code

        // X-type stabilizers
        stabilizers[[0, 0]] = 1;
        stabilizers[[0, 2]] = 1;
        stabilizers[[0, 4]] = 1;
        stabilizers[[0, 6]] = 1;
        stabilizers[[1, 1]] = 1;
        stabilizers[[1, 2]] = 1;
        stabilizers[[1, 5]] = 1;
        stabilizers[[1, 6]] = 1;
        stabilizers[[2, 3]] = 1;
        stabilizers[[2, 4]] = 1;
        stabilizers[[2, 5]] = 1;
        stabilizers[[2, 6]] = 1;

        // Z-type stabilizers
        stabilizers[[3, 7]] = 1;
        stabilizers[[3, 9]] = 1;
        stabilizers[[3, 11]] = 1;
        stabilizers[[3, 13]] = 1;
        stabilizers[[4, 8]] = 1;
        stabilizers[[4, 9]] = 1;
        stabilizers[[4, 12]] = 1;
        stabilizers[[4, 13]] = 1;
        stabilizers[[5, 10]] = 1;
        stabilizers[[5, 11]] = 1;
        stabilizers[[5, 12]] = 1;
        stabilizers[[5, 13]] = 1;

        Ok(stabilizers)
    }

    /// Generate Shor code stabilizers
    fn generate_shor_code_stabilizers(parameters: &CodeParameters) -> QECResult<Array2<u8>> {
        let n = 9; // Shor code is always 9 qubits
        let mut stabilizers = Array2::zeros((8, 18)); // 8 stabilizers for [[9,1,3]] code

        // X-type stabilizers (phase error detection)
        stabilizers[[0, 0]] = 1;
        stabilizers[[0, 1]] = 1;
        stabilizers[[0, 2]] = 1;
        stabilizers[[1, 3]] = 1;
        stabilizers[[1, 4]] = 1;
        stabilizers[[1, 5]] = 1;
        stabilizers[[2, 6]] = 1;
        stabilizers[[2, 7]] = 1;
        stabilizers[[2, 8]] = 1;

        // Z-type stabilizers (bit error detection)
        stabilizers[[3, 9]] = 1;
        stabilizers[[3, 10]] = 1;
        stabilizers[[4, 10]] = 1;
        stabilizers[[4, 11]] = 1;
        stabilizers[[5, 12]] = 1;
        stabilizers[[5, 13]] = 1;
        stabilizers[[6, 13]] = 1;
        stabilizers[[6, 14]] = 1;
        stabilizers[[7, 15]] = 1;
        stabilizers[[7, 16]] = 1;

        Ok(stabilizers)
    }

    /// Generate parity check matrix
    fn generate_parity_check_matrix(
        code: &ErrorCorrectionCode,
        parameters: &CodeParameters,
    ) -> QECResult<Array2<u8>> {
        // For most stabilizer codes, parity check matrix is derived from stabilizers
        let stabilizers = Self::generate_stabilizers(code, parameters)?;
        Ok(stabilizers) // Simplified - in practice would extract proper parity check
    }

    /// Build correction lookup table
    fn build_correction_table(
        code: &ErrorCorrectionCode,
        parameters: &CodeParameters,
        decoder: &DecodingAlgorithm,
    ) -> QECResult<HashMap<Vec<u8>, CorrectionOperation>> {
        let mut table = HashMap::new();

        match decoder {
            DecodingAlgorithm::LookupTable => {
                Self::build_lookup_table_decoder(code, parameters, &mut table)?;
            }
            DecodingAlgorithm::MinimumWeight => {
                // For minimum weight, we'll compute corrections on-demand
                // Just add trivial syndrome entry
                table.insert(
                    vec![0; parameters.num_physical_qubits - parameters.num_logical_qubits],
                    CorrectionOperation::identity(),
                );
            }
            _ => {
                // Other decoders compute corrections algorithmically
                table.insert(
                    vec![0; parameters.num_physical_qubits - parameters.num_logical_qubits],
                    CorrectionOperation::identity(),
                );
            }
        }

        Ok(table)
    }

    /// Build lookup table for small codes
    fn build_lookup_table_decoder(
        code: &ErrorCorrectionCode,
        parameters: &CodeParameters,
        table: &mut HashMap<Vec<u8>, CorrectionOperation>,
    ) -> QECResult<()> {
        let num_syndrome_bits = parameters.num_physical_qubits - parameters.num_logical_qubits;

        // For small codes, enumerate all possible syndromes
        if num_syndrome_bits <= 10 {
            // Limit to reasonable table size
            for syndrome_int in 0..(1 << num_syndrome_bits) {
                let syndrome = Self::int_to_syndrome(syndrome_int, num_syndrome_bits);
                let correction = Self::compute_minimum_weight_correction(&syndrome, parameters)?;
                table.insert(syndrome, correction);
            }
        }

        Ok(())
    }

    /// Convert integer to syndrome vector
    fn int_to_syndrome(syndrome_int: usize, num_bits: usize) -> Vec<u8> {
        (0..num_bits)
            .map(|i| ((syndrome_int >> i) & 1) as u8)
            .collect()
    }

    /// Compute minimum weight correction for syndrome
    fn compute_minimum_weight_correction(
        syndrome: &[u8],
        parameters: &CodeParameters,
    ) -> QECResult<CorrectionOperation> {
        // Simplified minimum weight decoding
        let mut corrections = Vec::new();

        // Find non-zero syndrome bits and create corresponding corrections
        for (i, &bit) in syndrome.iter().enumerate() {
            if bit == 1 {
                corrections.push(PauliCorrection {
                    qubit: i % parameters.num_physical_qubits,
                    operation: if i < parameters.num_physical_qubits {
                        PauliType::X
                    } else {
                        PauliType::Z
                    },
                    weight: 1.0,
                });
            }
        }

        // Derive decode-quality metrics from the correction weight relative to
        // the code's correctable radius t = (d-1)/2 rather than reporting flat
        // literals. A correction within the radius is fully trusted; trust
        // decays as the weight approaches and exceeds t (where the simplified
        // minimum-weight decoder is no longer guaranteed to succeed).
        let correction_weight = corrections.len();
        let correctable_radius = parameters.distance.saturating_sub(1) / 2;
        let (confidence, success_probability) = if correction_weight == 0 {
            (1.0, 1.0)
        } else if correctable_radius == 0 {
            // Distance-1 "code" cannot correct any error.
            (0.0, 0.0)
        } else {
            let ratio = correction_weight as f64 / correctable_radius as f64;
            // Within the radius (ratio <= 1) confidence stays high; beyond it the
            // metric decays smoothly towards zero.
            let confidence = (1.0 / (1.0 + (ratio - 1.0).max(0.0))).clamp(0.0, 1.0);
            let success_probability = (1.0 - 0.5 * ratio).clamp(0.0, 1.0);
            (confidence, success_probability)
        };

        Ok(CorrectionOperation {
            pauli_corrections: corrections.clone(),
            confidence,
            success_probability,
            required_resources: CorrectionResources {
                num_measurements: syndrome.len(),
                num_correction_gates: corrections.len(),
                time_overhead: 10.0 * corrections.len() as f64,
                ancilla_qubits: 0,
            },
        })
    }

    /// Measure the stabilizer syndrome of a quantum state.
    ///
    /// Each stabilizer generator is stored in binary symplectic form: row `i`
    /// has an `X` part in columns `[0, n)` and a `Z` part in columns
    /// `[n, 2n)`, so the per-qubit Pauli is `I`, `X`, `Z`, or `Y = XZ` depending
    /// on which bits are set. The syndrome bit is the genuine measurement
    /// outcome of that multi-qubit Pauli `P` on the state vector: it is `0` when
    /// the state lies in the `+1` eigenspace (`⟨ψ|P|ψ⟩ > 0`) and `1` when it
    /// lies in the `−1` eigenspace, i.e. an error anticommuting with the
    /// stabilizer has occurred. This reads the actual `state` rather than
    /// fabricating coin-flips from the noise rate.
    fn measure_syndrome(&self, state: &QuantumState) -> QECResult<Vec<u8>> {
        let num_stabilizers = self.stabilizer_generators.nrows();
        let num_columns = self.stabilizer_generators.ncols();
        let num_qubits = num_columns / 2;
        let mut syndrome = vec![0u8; num_stabilizers];

        for i in 0..num_stabilizers {
            let expectation = self.stabilizer_expectation(state, i, num_qubits)?;
            // +1 eigenvalue -> trivial (0); -1 eigenvalue -> error flagged (1).
            // Threshold at zero so partially-mixed/noisy states still classify.
            syndrome[i] = u8::from(expectation < 0.0);
        }

        Ok(syndrome)
    }

    /// Compute `Re⟨ψ|P_i|ψ⟩` for stabilizer generator `i` (a Hermitian
    /// multi-qubit Pauli) on the state vector. Applies the Pauli string to a
    /// working copy of the amplitudes and contracts with the original state.
    fn stabilizer_expectation(
        &self,
        state: &QuantumState,
        stabilizer_index: usize,
        num_qubits: usize,
    ) -> QECResult<f64> {
        if num_qubits != state.num_qubits {
            return Err(QuantumErrorCorrectionError::SyndromeError(format!(
                "Stabilizer acts on {num_qubits} qubits but state has {}",
                state.num_qubits
            )));
        }

        // Apply the multi-qubit Pauli to a copy of the state, qubit by qubit.
        let mut transformed = state.amplitudes.clone();
        for qubit in 0..num_qubits {
            let x_bit = self.stabilizer_generators[[stabilizer_index, qubit]] != 0;
            let z_bit = self.stabilizer_generators[[stabilizer_index, num_qubits + qubit]] != 0;
            let pauli = match (x_bit, z_bit) {
                (false, false) => PauliType::I,
                (true, false) => PauliType::X,
                (false, true) => PauliType::Z,
                (true, true) => PauliType::Y,
            };
            apply_pauli_to_amplitudes(&mut transformed, qubit, pauli);
        }

        // ⟨ψ|Pψ⟩ = Σ conj(ψ_k) · (Pψ)_k. For a Hermitian Pauli the result is
        // real up to floating-point error, so take the real part.
        let expectation: f64 = state
            .amplitudes
            .iter()
            .zip(transformed.iter())
            .map(|(psi, p_psi)| (psi.conj() * p_psi).re)
            .sum();

        Ok(expectation)
    }

    /// Check if syndrome is trivial (all zeros)
    fn is_trivial_syndrome(&self, syndrome: &[u8]) -> bool {
        syndrome.iter().all(|&bit| bit == 0)
    }

    /// Decode syndrome to find correction
    fn decode_syndrome(&self, syndrome: &[u8]) -> QECResult<CorrectionOperation> {
        match self.config.decoder {
            DecodingAlgorithm::LookupTable => {
                if let Some(correction) = self.correction_table.get(syndrome) {
                    Ok(correction.clone())
                } else {
                    Self::compute_minimum_weight_correction(syndrome, &self.parameters)
                }
            }
            DecodingAlgorithm::MinimumWeight => {
                Self::compute_minimum_weight_correction(syndrome, &self.parameters)
            }
            _ => {
                // Fallback to minimum weight for other decoders
                Self::compute_minimum_weight_correction(syndrome, &self.parameters)
            }
        }
    }

    /// Extract error locations from correction
    fn extract_error_locations(&self, correction: &CorrectionOperation) -> Vec<usize> {
        correction
            .pauli_corrections
            .iter()
            .map(|pc| pc.qubit)
            .collect()
    }

    /// Calculate detection confidence
    fn calculate_detection_confidence(
        &self,
        syndrome: &[u8],
        error_locations: &[usize],
    ) -> QECResult<f64> {
        // Base confidence on syndrome weight and consistency
        let syndrome_weight = syndrome.iter().map(|&b| b as usize).sum::<usize>();
        let max_weight = syndrome.len();

        // Lower confidence for higher syndrome weights (more errors)
        let weight_factor = 1.0 - (syndrome_weight as f64 / max_weight as f64);

        // Factor in noise model
        let noise_factor = 1.0 - self.config.noise_model.measurement_error_rate;

        Ok(weight_factor * noise_factor * 0.9) // Cap at 0.9
    }

    /// Estimate noise level from syndrome
    fn estimate_noise_level(&self, syndrome: &[u8]) -> f64 {
        let syndrome_weight = syndrome.iter().map(|&b| f64::from(b)).sum::<f64>();
        syndrome_weight / syndrome.len() as f64
    }

    /// Apply a single-qubit Pauli correction to the state vector in place.
    ///
    /// Acts on the amplitude vector (indexed so that qubit `q` is bit `q`, the
    /// same convention as [`QuantumState::expectation_z`]):
    ///
    /// * `X` swaps the amplitudes of basis states that differ in qubit `q`;
    /// * `Z` negates the amplitude of every basis state with qubit `q` set;
    /// * `Y = iXZ`: `Y|0⟩ = i|1⟩`, `Y|1⟩ = −i|0⟩`;
    /// * `I` is a no-op.
    ///
    /// Returns [`QuantumErrorCorrectionError::SyndromeError`] if the target
    /// qubit is out of range for the state.
    fn apply_pauli_correction(
        &self,
        state: &mut QuantumState,
        correction: &PauliCorrection,
    ) -> QECResult<()> {
        let qubit = correction.qubit;
        if qubit >= state.num_qubits {
            return Err(QuantumErrorCorrectionError::SyndromeError(format!(
                "Correction targets qubit {qubit} but state has {} qubits",
                state.num_qubits
            )));
        }

        let bit = 1usize << qubit;
        let dim = state.amplitudes.len();

        match correction.operation {
            PauliType::I => {}
            PauliType::X => {
                // Swap |..0..⟩ <-> |..1..⟩ for the target bit.
                for index in 0..dim {
                    if index & bit == 0 {
                        let partner = index | bit;
                        state.amplitudes.swap(index, partner);
                    }
                }
            }
            PauliType::Z => {
                // Phase flip on the |1⟩ component of the target qubit.
                for (index, amplitude) in state.amplitudes.iter_mut().enumerate() {
                    if index & bit != 0 {
                        *amplitude = -*amplitude;
                    }
                }
            }
            PauliType::Y => {
                // Y|0> = i|1>, Y|1> = -i|0>.
                let imag = scirs2_core::Complex64::new(0.0, 1.0);
                for index in 0..dim {
                    if index & bit == 0 {
                        let partner = index | bit;
                        let a0 = state.amplitudes[index];
                        let a1 = state.amplitudes[partner];
                        state.amplitudes[index] = -imag * a1;
                        state.amplitudes[partner] = imag * a0;
                    }
                }
            }
        }

        Ok(())
    }

    /// Update detection statistics
    fn update_detection_stats(
        &mut self,
        syndrome: &[u8],
        correction: Option<&CorrectionOperation>,
        confidence: f64,
        detection_time: Duration,
    ) {
        self.detection_stats.total_measurements += 1;

        if !self.is_trivial_syndrome(syndrome) {
            self.detection_stats.errors_detected += 1;
        }

        self.detection_stats.avg_detection_time = Self::update_average_time(
            self.detection_stats.avg_detection_time,
            detection_time,
            self.detection_stats.total_measurements,
        );

        // Record syndrome
        let record = SyndromeRecord {
            timestamp: Instant::now(),
            syndrome: syndrome.to_vec(),
            error_locations: correction.map(|c| self.extract_error_locations(c)),
            correction_applied: correction.cloned(),
            // The genuine per-detection confidence computed in `detect_syndrome`
            // from the syndrome weight and noise model.
            confidence,
        };

        self.detection_stats.syndrome_history.push(record);

        // Limit history size
        if self.detection_stats.syndrome_history.len() > self.config.max_history_length {
            self.detection_stats.syndrome_history.remove(0);
        }
    }

    /// Update average time calculation
    fn update_average_time(current_avg: Duration, new_time: Duration, count: usize) -> Duration {
        let current_total = current_avg * (count - 1) as u32;
        let new_total = current_total + new_time;
        new_total / count as u32
    }

    /// Get plaquette qubits for surface code
    fn get_plaquette_qubits(row: usize, col: usize, d: usize) -> Vec<usize> {
        let mut qubits = Vec::new();

        // Add qubits around the plaquette
        if row > 0 {
            qubits.push(row * d + col);
        }
        if col > 0 {
            qubits.push(row * d + col - 1);
        }
        if row < d - 1 {
            qubits.push((row + 1) * d + col);
        }
        if col < d - 1 {
            qubits.push(row * d + col + 1);
        }

        qubits
    }

    /// Get vertex qubits for surface code
    fn get_vertex_qubits(row: usize, col: usize, d: usize) -> Vec<usize> {
        // Similar to plaquette but for vertex checks
        Self::get_plaquette_qubits(row, col, d)
    }
}

impl CorrectionOperation {
    /// Create identity correction (no operation)
    #[must_use]
    pub const fn identity() -> Self {
        Self {
            pauli_corrections: Vec::new(),
            confidence: 1.0,
            success_probability: 1.0,
            required_resources: CorrectionResources {
                num_measurements: 0,
                num_correction_gates: 0,
                time_overhead: 0.0,
                ancilla_qubits: 0,
            },
        }
    }
}

impl DetectionStatistics {
    /// Create new detection statistics
    #[must_use]
    pub const fn new() -> Self {
        Self {
            total_measurements: 0,
            errors_detected: 0,
            successful_corrections: 0,
            failed_corrections: 0,
            avg_detection_time: Duration::from_nanos(0),
            avg_correction_time: Duration::from_nanos(0),
            syndrome_history: Vec::new(),
        }
    }

    /// Get error detection rate
    #[must_use]
    pub fn error_detection_rate(&self) -> f64 {
        if self.total_measurements == 0 {
            0.0
        } else {
            self.errors_detected as f64 / self.total_measurements as f64
        }
    }

    /// Get correction success rate
    #[must_use]
    pub fn correction_success_rate(&self) -> f64 {
        let total_corrections = self.successful_corrections + self.failed_corrections;
        if total_corrections == 0 {
            0.0
        } else {
            self.successful_corrections as f64 / total_corrections as f64
        }
    }
}

impl Default for SyndromeDetectorConfig {
    fn default() -> Self {
        Self {
            measurement_frequency: 1000.0, // 1 kHz
            correction_threshold: 0.7,
            max_history_length: 1000,
            adaptive_timing: true,
            noise_model: NoiseModel::default(),
            decoder: DecodingAlgorithm::MinimumWeight,
        }
    }
}

impl Default for NoiseModel {
    fn default() -> Self {
        Self {
            single_qubit_error_rate: 0.001,
            two_qubit_error_rate: 0.01,
            measurement_error_rate: 0.02,
            coherence_time: 100.0, // microseconds
            thermal_noise_rate: 0.0001,
        }
    }
}

/// Apply a single-qubit Pauli to a raw amplitude vector in place.
///
/// Uses the same bit convention as
/// [`SyndromeDetector::apply_pauli_correction`] (qubit `q` is bit `q`):
/// * `X` swaps amplitudes whose target bit differs;
/// * `Z` negates amplitudes whose target bit is set;
/// * `Y = iXZ` with `Y|0⟩ = i|1⟩`, `Y|1⟩ = −i|0⟩`;
/// * `I` is a no-op.
fn apply_pauli_to_amplitudes(amplitudes: &mut [Complex64], qubit: usize, pauli: PauliType) {
    let bit = 1usize << qubit;
    let dim = amplitudes.len();

    match pauli {
        PauliType::I => {}
        PauliType::X => {
            for index in 0..dim {
                if index & bit == 0 {
                    amplitudes.swap(index, index | bit);
                }
            }
        }
        PauliType::Z => {
            for (index, amplitude) in amplitudes.iter_mut().enumerate() {
                if index & bit != 0 {
                    *amplitude = -*amplitude;
                }
            }
        }
        PauliType::Y => {
            let imag = Complex64::new(0.0, 1.0);
            for index in 0..dim {
                if index & bit == 0 {
                    let partner = index | bit;
                    let a0 = amplitudes[index];
                    let a1 = amplitudes[partner];
                    amplitudes[index] = -imag * a1;
                    amplitudes[partner] = imag * a0;
                }
            }
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::quantum_error_correction::codes::*;

    #[test]
    fn test_syndrome_detector_creation() {
        let code = ErrorCorrectionCode::RepetitionCode;
        let parameters = CodeParameters {
            distance: 3,
            num_logical_qubits: 1,
            num_physical_qubits: 3,
            num_ancilla_qubits: 2,
            code_rate: 1.0 / 3.0,
            threshold_probability: 0.1,
        };
        let config = SyndromeDetectorConfig::default();

        let detector = SyndromeDetector::new(code, parameters, config);
        assert!(detector.is_ok());
    }

    #[test]
    fn test_repetition_code_stabilizers() {
        let parameters = CodeParameters {
            distance: 3,
            num_logical_qubits: 1,
            num_physical_qubits: 3,
            num_ancilla_qubits: 2,
            code_rate: 1.0 / 3.0,
            threshold_probability: 0.1,
        };

        let stabilizers = SyndromeDetector::generate_repetition_code_stabilizers(&parameters)
            .expect("should generate repetition code stabilizers");
        assert_eq!(stabilizers.nrows(), 2); // n-1 stabilizers
        assert_eq!(stabilizers.ncols(), 6); // 2n columns (X and Z parts)
    }

    #[test]
    fn test_trivial_syndrome_detection() {
        let detector = create_test_detector();
        let syndrome = vec![0, 0, 0];
        assert!(detector.is_trivial_syndrome(&syndrome));

        let non_trivial = vec![1, 0, 1];
        assert!(!detector.is_trivial_syndrome(&non_trivial));
    }

    #[test]
    fn test_minimum_weight_correction() {
        let parameters = CodeParameters {
            distance: 3,
            num_logical_qubits: 1,
            num_physical_qubits: 3,
            num_ancilla_qubits: 2,
            code_rate: 1.0 / 3.0,
            threshold_probability: 0.1,
        };

        let syndrome = vec![1, 0];
        let correction =
            SyndromeDetector::compute_minimum_weight_correction(&syndrome, &parameters)
                .expect("should compute minimum weight correction");
        assert_eq!(correction.pauli_corrections.len(), 1);
        assert_eq!(correction.pauli_corrections[0].qubit, 0);
    }

    fn create_test_detector() -> SyndromeDetector {
        let code = ErrorCorrectionCode::RepetitionCode;
        let parameters = CodeParameters {
            distance: 3,
            num_logical_qubits: 1,
            num_physical_qubits: 3,
            num_ancilla_qubits: 2,
            code_rate: 1.0 / 3.0,
            threshold_probability: 0.1,
        };
        let config = SyndromeDetectorConfig::default();

        SyndromeDetector::new(code, parameters, config)
            .expect("should create syndrome detector for testing")
    }
}