quantrs2_ml/error_mitigation/

mitigator.rs

//! Core mitigation logic: [`QuantumMLErrorMitigator`]'s constructors, the
//! strategy-dispatch pipeline (`mitigate_training_errors`/`mitigate_inference_errors`),
//! the real ZNE/readout-error-correction implementations, and the adaptive-strategy
//! machinery. Data types live in `super::types`.

use super::types::*;
use crate::error::{MLError, Result};
use quantrs2_circuit::builder::Simulator;
use quantrs2_circuit::prelude::Circuit;
use quantrs2_core::qubit::QubitId;
use quantrs2_sim::noise::NoiseModelBuilder;
use quantrs2_sim::statevector::StateVectorSimulator;
use scirs2_core::ndarray::{Array1, Array2, Axis};
use scirs2_core::random::prelude::*;
use scirs2_core::Complex64;

/// Number of measurement shots simulated per circuit execution during error
/// mitigation (matches the shape of the `Array2` measurement buffers used
/// throughout this module).
const DEFAULT_NUM_SHOTS: usize = 100;

impl QuantumMLErrorMitigator {
    /// Create a new error mitigation framework
    pub fn new(mitigation_strategy: MitigationStrategy, noise_model: NoiseModel) -> Result<Self> {
        let calibration_data = CalibrationData::default();
        let adaptive_config = AdaptiveConfig::default();
        let performance_metrics = PerformanceMetrics::new();

        Ok(Self {
            mitigation_strategy,
            noise_model,
            calibration_data,
            adaptive_config,
            performance_metrics,
        })
    }

    /// Apply error mitigation to quantum ML training
    pub fn mitigate_training_errors(
        &mut self,
        circuit: &QuantumCircuit,
        parameters: &Array1<f64>,
        measurement_results: &Array2<f64>,
        gradient_estimates: &Array1<f64>,
    ) -> Result<MitigatedTrainingData> {
        // Update noise model based on current measurements
        self.update_noise_model(measurement_results)?;

        // Apply mitigation strategy
        let mitigated_measurements =
            self.apply_measurement_mitigation(circuit, measurement_results)?;

        let mitigated_gradients =
            self.apply_gradient_mitigation(circuit, parameters, gradient_estimates)?;

        // Update performance metrics
        self.performance_metrics
            .update(&mitigated_measurements, &mitigated_gradients)?;

        // Adaptive strategy adjustment
        if self.should_adapt_strategy()? {
            self.adapt_mitigation_strategy()?;
        }

        Ok(MitigatedTrainingData {
            measurements: mitigated_measurements,
            gradients: mitigated_gradients,
            confidence_scores: self.compute_confidence_scores(circuit)?,
            mitigation_overhead: self.performance_metrics.mitigation_overhead,
        })
    }

    /// Apply error mitigation to quantum ML inference
    pub fn mitigate_inference_errors(
        &mut self,
        circuit: &QuantumCircuit,
        measurement_results: &Array2<f64>,
    ) -> Result<MitigatedInferenceData> {
        let mitigated_measurements =
            self.apply_measurement_mitigation(circuit, measurement_results)?;

        let uncertainty_estimates =
            self.compute_uncertainty_estimates(circuit, &mitigated_measurements)?;

        Ok(MitigatedInferenceData {
            measurements: mitigated_measurements,
            uncertainty: uncertainty_estimates,
            reliability_score: self.compute_reliability_score(circuit)?,
        })
    }

    /// Apply measurement error mitigation
    fn apply_measurement_mitigation(
        &self,
        circuit: &QuantumCircuit,
        measurements: &Array2<f64>,
    ) -> Result<Array2<f64>> {
        match &self.mitigation_strategy {
            MitigationStrategy::ZNE {
                scale_factors,
                extrapolation_method,
                ..
            } => self.apply_zne_mitigation(
                circuit,
                measurements,
                scale_factors,
                extrapolation_method,
            ),
            MitigationStrategy::ReadoutErrorMitigation {
                calibration_matrix,
                correction_method,
                ..
            } => self.apply_readout_error_mitigation(
                measurements,
                calibration_matrix,
                correction_method,
            ),
            MitigationStrategy::CDR {
                training_circuits,
                regression_model,
                ..
            } => self.apply_cdr_mitigation(
                circuit,
                measurements,
                training_circuits,
                regression_model,
            ),
            MitigationStrategy::SymmetryVerification {
                symmetry_groups, ..
            } => self.apply_symmetry_verification(circuit, measurements, symmetry_groups),
            MitigationStrategy::VirtualDistillation {
                distillation_rounds,
                ..
            } => self.apply_virtual_distillation(circuit, measurements, *distillation_rounds),
            MitigationStrategy::MLMitigation {
                noise_predictor,
                correction_network,
                ..
            } => {
                self.apply_ml_mitigation(circuit, measurements, noise_predictor, correction_network)
            }
            MitigationStrategy::HybridErrorCorrection {
                classical_preprocessing,
                quantum_correction,
                post_processing,
            } => self.apply_hybrid_error_correction(
                circuit,
                measurements,
                classical_preprocessing,
                quantum_correction,
                post_processing,
            ),
            MitigationStrategy::AdaptiveMultiStrategy {
                strategies,
                selection_policy,
                ..
            } => self.apply_adaptive_multi_strategy(
                circuit,
                measurements,
                strategies,
                selection_policy,
            ),
        }
    }

    /// Apply Zero Noise Extrapolation
    fn apply_zne_mitigation(
        &self,
        circuit: &QuantumCircuit,
        measurements: &Array2<f64>,
        scale_factors: &[f64],
        extrapolation_method: &ExtrapolationMethod,
    ) -> Result<Array2<f64>> {
        let mut scaled_results = Vec::new();

        for &scale_factor in scale_factors {
            let scaled_circuit = self.scale_circuit_noise(circuit, scale_factor)?;
            let scaled_measurements = self.execute_scaled_circuit(&scaled_circuit)?;
            scaled_results.push((scale_factor, scaled_measurements));
        }

        // Extrapolate to zero noise
        self.extrapolate_to_zero_noise(&scaled_results, extrapolation_method)
    }

    /// Apply readout error mitigation
    fn apply_readout_error_mitigation(
        &self,
        measurements: &Array2<f64>,
        calibration_matrix: &Array2<f64>,
        correction_method: &ReadoutCorrectionMethod,
    ) -> Result<Array2<f64>> {
        match correction_method {
            ReadoutCorrectionMethod::MatrixInversion => {
                self.apply_matrix_inversion_correction(measurements, calibration_matrix)
            }
            ReadoutCorrectionMethod::ConstrainedLeastSquares => {
                self.apply_constrained_least_squares_correction(measurements, calibration_matrix)
            }
            ReadoutCorrectionMethod::IterativeMaximumLikelihood => {
                self.apply_ml_correction(measurements, calibration_matrix)
            }
        }
    }

    /// Apply Clifford Data Regression.
    ///
    /// Not yet implemented: [`CliffordCircuit`] carries no actual circuit
    /// data to execute or regress against, so this would otherwise have to
    /// fabricate training labels. Honestly reports
    /// [`MLError::NotSupported`] rather than silently returning the raw,
    /// uncorrected measurements.
    fn apply_cdr_mitigation(
        &self,
        _circuit: &QuantumCircuit,
        _measurements: &Array2<f64>,
        _training_circuits: &[CliffordCircuit],
        _regression_model: &CDRModel,
    ) -> Result<Array2<f64>> {
        Err(MLError::NotSupported(
            "Clifford Data Regression is not yet implemented (CliffordCircuit has no \
             executable circuit data to train a real regression model against); use \
             MitigationStrategy::ReadoutErrorMitigation or ZNE instead"
                .to_string(),
        ))
    }

    /// Apply symmetry verification.
    ///
    /// Not yet implemented: honestly reports [`MLError::NotSupported`]
    /// rather than returning the input measurements unchanged while
    /// claiming symmetry-based post-selection occurred.
    fn apply_symmetry_verification(
        &self,
        _circuit: &QuantumCircuit,
        _measurements: &Array2<f64>,
        _symmetry_groups: &[SymmetryGroup],
    ) -> Result<Array2<f64>> {
        Err(MLError::NotSupported(
            "Symmetry verification mitigation is not yet implemented; use \
             MitigationStrategy::ReadoutErrorMitigation or ZNE instead"
                .to_string(),
        ))
    }

    /// Apply virtual distillation.
    ///
    /// Not yet implemented (requires simulating multiple entangled circuit
    /// copies, which this module's circuit representation does not support
    /// yet). Honestly reports [`MLError::NotSupported`].
    fn apply_virtual_distillation(
        &self,
        _circuit: &QuantumCircuit,
        _measurements: &Array2<f64>,
        _distillation_rounds: usize,
    ) -> Result<Array2<f64>> {
        Err(MLError::NotSupported(
            "Virtual distillation is not yet implemented (requires multi-copy entangled \
             circuit simulation); use MitigationStrategy::ReadoutErrorMitigation or ZNE \
             instead"
                .to_string(),
        ))
    }

    /// Apply ML-based error mitigation.
    ///
    /// Not yet implemented: [`NoisePredictorModel`]/[`CorrectionNetwork`]
    /// carry no trained weights, so this would otherwise fabricate
    /// "corrections". Honestly reports [`MLError::NotSupported`].
    fn apply_ml_mitigation(
        &self,
        _circuit: &QuantumCircuit,
        _measurements: &Array2<f64>,
        _noise_predictor: &NoisePredictorModel,
        _correction_network: &CorrectionNetwork,
    ) -> Result<Array2<f64>> {
        Err(MLError::NotSupported(
            "ML-based error mitigation is not yet implemented (no trained noise-predictor \
             or correction-network weights are available); use \
             MitigationStrategy::ReadoutErrorMitigation or ZNE instead"
                .to_string(),
        ))
    }

    /// Apply hybrid classical-quantum error correction.
    ///
    /// Not yet implemented: honestly reports [`MLError::NotSupported`]
    /// rather than passing measurements through unmodified classical
    /// pre/post-processing stubs while claiming quantum error correction
    /// occurred.
    fn apply_hybrid_error_correction(
        &self,
        _circuit: &QuantumCircuit,
        _measurements: &Array2<f64>,
        _classical_preprocessing: &ClassicalPreprocessor,
        _quantum_correction: &QuantumErrorCorrector,
        _post_processing: &ClassicalPostprocessor,
    ) -> Result<Array2<f64>> {
        Err(MLError::NotSupported(
            "Hybrid classical-quantum error correction is not yet implemented; use \
             MitigationStrategy::ReadoutErrorMitigation or ZNE instead"
                .to_string(),
        ))
    }

    /// Apply adaptive multi-strategy mitigation
    fn apply_adaptive_multi_strategy(
        &self,
        circuit: &QuantumCircuit,
        measurements: &Array2<f64>,
        strategies: &[MitigationStrategy],
        selection_policy: &StrategySelectionPolicy,
    ) -> Result<Array2<f64>> {
        // Select best strategy based on current circuit and performance history
        let selected_strategy =
            selection_policy.select_strategy(circuit, &self.performance_metrics, strategies)?;

        // Apply selected strategy
        let mitigator = QuantumMLErrorMitigator {
            mitigation_strategy: selected_strategy,
            noise_model: self.noise_model.clone(),
            calibration_data: self.calibration_data.clone(),
            adaptive_config: self.adaptive_config.clone(),
            performance_metrics: self.performance_metrics.clone(),
        };

        mitigator.apply_measurement_mitigation(circuit, measurements)
    }

    /// Apply gradient error mitigation
    fn apply_gradient_mitigation(
        &self,
        circuit: &QuantumCircuit,
        parameters: &Array1<f64>,
        gradients: &Array1<f64>,
    ) -> Result<Array1<f64>> {
        // Parameter shift rule with error mitigation
        let mut mitigated_gradients = Array1::zeros(gradients.len());

        for (i, &param) in parameters.iter().enumerate() {
            // Create shifted circuits
            let mut params_plus = parameters.clone();
            let mut params_minus = parameters.clone();
            params_plus[i] = param + std::f64::consts::PI / 2.0;
            params_minus[i] = param - std::f64::consts::PI / 2.0;

            // Apply error mitigation to shifted measurements
            let circuit_plus = circuit.with_parameters(&params_plus)?;
            let circuit_minus = circuit.with_parameters(&params_minus)?;

            let measurements_plus = self.measure_circuit(&circuit_plus)?;
            let measurements_minus = self.measure_circuit(&circuit_minus)?;

            let mitigated_plus =
                self.apply_measurement_mitigation(&circuit_plus, &measurements_plus)?;
            let mitigated_minus =
                self.apply_measurement_mitigation(&circuit_minus, &measurements_minus)?;

            // Compute mitigated gradient
            mitigated_gradients[i] = (mitigated_plus.mean().unwrap_or(0.0)
                - mitigated_minus.mean().unwrap_or(0.0))
                / 2.0;
        }

        Ok(mitigated_gradients)
    }

    /// Update noise model based on current measurements
    fn update_noise_model(&mut self, measurements: &Array2<f64>) -> Result<()> {
        // Analyze measurement statistics to infer noise characteristics
        let noise_statistics = self.analyze_noise_statistics(measurements)?;

        // Update gate error models
        for (gate_name, error_model) in &mut self.noise_model.gate_errors {
            error_model.update_from_statistics(&noise_statistics)?;
        }

        // Update measurement error model
        self.noise_model
            .measurement_errors
            .update_from_measurements(measurements)?;

        Ok(())
    }

    /// Check if mitigation strategy should be adapted
    fn should_adapt_strategy(&self) -> Result<bool> {
        let current_performance = self.performance_metrics.current_performance();
        let adaptation_threshold = self.adaptive_config.performance_threshold;

        Ok(current_performance < adaptation_threshold)
    }

    /// Adapt mitigation strategy based on performance
    fn adapt_mitigation_strategy(&mut self) -> Result<()> {
        match &self.adaptive_config.strategy_switching_policy {
            SwitchingPolicy::PerformanceBased => {
                self.switch_to_best_performing_strategy()?;
            }
            SwitchingPolicy::ResourceOptimized => {
                self.switch_to_resource_optimal_strategy()?;
            }
            SwitchingPolicy::HybridAdaptive => {
                self.switch_to_hybrid_adaptive_strategy()?;
            }
        }

        Ok(())
    }

    /// Compute confidence scores for mitigation results
    fn compute_confidence_scores(&self, circuit: &QuantumCircuit) -> Result<Array1<f64>> {
        let circuit_complexity = self.assess_circuit_complexity(circuit)?;
        let noise_level = self.estimate_noise_level(circuit)?;
        let mitigation_effectiveness = self.estimate_mitigation_effectiveness()?;

        let base_confidence = 1.0 - (circuit_complexity * noise_level);
        let adjusted_confidence = base_confidence * mitigation_effectiveness;

        Ok(Array1::from_elem(circuit.num_qubits(), adjusted_confidence))
    }

    /// Compute uncertainty estimates
    fn compute_uncertainty_estimates(
        &self,
        circuit: &QuantumCircuit,
        measurements: &Array2<f64>,
    ) -> Result<Array1<f64>> {
        // Bootstrap sampling for uncertainty estimation
        let num_bootstrap_samples = 1000;
        let mut bootstrap_results = Vec::new();

        for _ in 0..num_bootstrap_samples {
            let bootstrap_measurements = self.bootstrap_sample(measurements)?;
            let mitigated_bootstrap =
                self.apply_measurement_mitigation(circuit, &bootstrap_measurements)?;
            bootstrap_results.push(mitigated_bootstrap.mean().unwrap_or(0.0));
        }

        // Compute standard deviation as uncertainty
        let mean_result = bootstrap_results.iter().sum::<f64>() / bootstrap_results.len() as f64;
        let variance = bootstrap_results
            .iter()
            .map(|&x| (x - mean_result).powi(2))
            .sum::<f64>()
            / bootstrap_results.len() as f64;
        let uncertainty = variance.sqrt();

        Ok(Array1::from_elem(1, uncertainty))
    }

    /// Compute reliability score
    fn compute_reliability_score(&self, circuit: &QuantumCircuit) -> Result<f64> {
        let mitigation_fidelity = self.estimate_mitigation_fidelity(circuit)?;
        let noise_resilience = self.assess_noise_resilience(circuit)?;
        let calibration_quality = self.assess_calibration_quality()?;

        Ok(mitigation_fidelity * noise_resilience * calibration_quality)
    }

    // Helper methods for implementation details...

    /// Real global unitary folding: `C -> C (C^-1 C)^k` where
    /// `k = round((scale_factor - 1) / 2)`. Each fold inserts an exact
    /// inverse/forward gate pair, so the *ideal* (noiseless) circuit output
    /// is mathematically unchanged, while the real physical gate count
    /// (and hence the noise accumulated by [`Self::execute_scaled_circuit`])
    /// grows with `scale_factor`, exactly matching standard ZNE folding.
    fn scale_circuit_noise(
        &self,
        circuit: &QuantumCircuit,
        scale_factor: f64,
    ) -> Result<QuantumCircuit> {
        if !scale_factor.is_finite() || scale_factor < 1.0 {
            return Err(MLError::InvalidParameter(format!(
                "ZNE scale factor must be finite and >= 1.0, got {scale_factor}"
            )));
        }
        let folds = (((scale_factor - 1.0) / 2.0).round().max(0.0)) as usize;

        let mut folded_gates = circuit.gates.clone();
        for _ in 0..folds {
            let mut inverse_layer: Vec<QuantumGate> = circuit
                .gates
                .iter()
                .rev()
                .map(Self::inverse_gate)
                .collect::<Result<Vec<_>>>()?;
            folded_gates.append(&mut inverse_layer);
            folded_gates.extend(circuit.gates.iter().cloned());
        }

        Ok(QuantumCircuit {
            gates: folded_gates,
            qubits: circuit.qubits,
        })
    }

    /// The exact inverse of a single supported gate, used by
    /// [`Self::scale_circuit_noise`]'s unitary folding. Self-inverse gates
    /// (H, Pauli, CNOT) are returned unchanged; rotation gates are returned
    /// with their angle negated. Unsupported gates honestly error instead of
    /// silently being treated as self-inverse.
    fn inverse_gate(gate: &QuantumGate) -> Result<QuantumGate> {
        match gate.name.as_str() {
            "H" | "X" | "Y" | "Z" | "CNOT" | "CX" => Ok(gate.clone()),
            "RX" | "RY" | "RZ" => Ok(QuantumGate {
                name: gate.name.clone(),
                qubits: gate.qubits.clone(),
                parameters: gate.parameters.mapv(|theta| -theta),
            }),
            other => Err(MLError::NotSupported(format!(
                "gate '{other}' has no known inverse for ZNE unitary folding"
            ))),
        }
    }

    /// Convert this module's lightweight [`QuantumCircuit`]/[`QuantumGate`]
    /// representation into a real, executable `quantrs2_circuit::Circuit<N>`.
    fn build_real_circuit<const N: usize>(circuit: &QuantumCircuit) -> Result<Circuit<N>> {
        let mut real_circuit = Circuit::<N>::new();
        for gate in &circuit.gates {
            Self::apply_gate_to_circuit(&mut real_circuit, gate)?;
        }
        Ok(real_circuit)
    }

    /// Apply a single [`QuantumGate`] to a real circuit builder, honestly
    /// erroring on any gate name the simulator backend does not recognize
    /// rather than silently dropping it.
    fn apply_gate_to_circuit<const N: usize>(
        real_circuit: &mut Circuit<N>,
        gate: &QuantumGate,
    ) -> Result<()> {
        let qubit_at = |idx: usize| -> Result<usize> {
            gate.qubits.get(idx).copied().ok_or_else(|| {
                MLError::InvalidConfiguration(format!(
                    "gate '{}' expects at least {} qubit argument(s)",
                    gate.name,
                    idx + 1
                ))
            })
        };
        let angle_at = |idx: usize| -> f64 { gate.parameters.get(idx).copied().unwrap_or(0.0) };

        match gate.name.as_str() {
            "H" => {
                real_circuit.h(qubit_at(0)?)?;
            }
            "X" => {
                real_circuit.x(qubit_at(0)?)?;
            }
            "Y" => {
                real_circuit.y(qubit_at(0)?)?;
            }
            "Z" => {
                real_circuit.z(qubit_at(0)?)?;
            }
            "RX" => {
                real_circuit.rx(qubit_at(0)?, angle_at(0))?;
            }
            "RY" => {
                real_circuit.ry(qubit_at(0)?, angle_at(0))?;
            }
            "RZ" => {
                real_circuit.rz(qubit_at(0)?, angle_at(0))?;
            }
            "CNOT" | "CX" => {
                real_circuit.cnot(qubit_at(0)?, qubit_at(1)?)?;
            }
            other => {
                return Err(MLError::NotSupported(format!(
                    "gate '{other}' is not supported by the error-mitigation simulator backend"
                )));
            }
        }
        Ok(())
    }

    /// Average calibrated per-gate error rate across `self.noise_model`
    /// (`0.0`, i.e. no injected noise, if no gate calibration was supplied).
    fn average_gate_error_rate(&self) -> f64 {
        let rates: Vec<f64> = self
            .noise_model
            .gate_errors
            .values()
            .map(|model| model.error_rate)
            .collect();
        if rates.is_empty() {
            0.0
        } else {
            rates.iter().sum::<f64>() / rates.len() as f64
        }
    }

    /// Simulate `circuit` noiselessly on the real state-vector backend to
    /// get its exact final amplitudes.
    fn run_noiseless_amplitudes(circuit: &QuantumCircuit) -> Result<Vec<Complex64>> {
        let num_qubits = circuit.num_qubits();
        match num_qubits {
            0 => Err(MLError::InvalidConfiguration(
                "circuit has no qubits".to_string(),
            )),
            1..=2 => Ok(StateVectorSimulator::new()
                .run(&Self::build_real_circuit::<2>(circuit)?)?
                .amplitudes()
                .to_vec()),
            3..=4 => Ok(StateVectorSimulator::new()
                .run(&Self::build_real_circuit::<4>(circuit)?)?
                .amplitudes()
                .to_vec()),
            5..=8 => Ok(StateVectorSimulator::new()
                .run(&Self::build_real_circuit::<8>(circuit)?)?
                .amplitudes()
                .to_vec()),
            9..=16 => Ok(StateVectorSimulator::new()
                .run(&Self::build_real_circuit::<16>(circuit)?)?
                .amplitudes()
                .to_vec()),
            n => Err(MLError::NotSupported(format!(
                "error-mitigation simulator backend supports at most 16 qubits, got {n}"
            ))),
        }
    }

    /// Sample one projective measurement outcome (a computational-basis
    /// index) from a state vector via inverse-CDF sampling of the real Born
    /// rule distribution.
    fn sample_basis_state(amplitudes: &[Complex64]) -> usize {
        let r: f64 = thread_rng().random::<f64>();
        let mut cumulative = 0.0;
        for (idx, amp) in amplitudes.iter().enumerate() {
            cumulative += amp.norm_sqr();
            if r < cumulative {
                return idx;
            }
        }
        amplitudes.len().saturating_sub(1)
    }

    /// Simulate `circuit` on the real state-vector backend and produce
    /// `num_shots` independent noisy measurement trajectories: each shot
    /// clones the exact noiseless final state, applies one Monte-Carlo
    /// realization of a depolarizing-noise channel (on every qubit, with
    /// per-gate error rate `per_gate_error_rate` accumulated across the
    /// circuit's real gate count via `1 - (1 - rate)^num_gates`), and then
    /// samples one projective measurement outcome from the resulting state.
    fn simulate_circuit_shots(
        &self,
        circuit: &QuantumCircuit,
        per_gate_error_rate: f64,
        num_shots: usize,
    ) -> Result<Array2<f64>> {
        let num_qubits = circuit.num_qubits();
        let noiseless_amplitudes = Self::run_noiseless_amplitudes(circuit)?;

        let num_gates = circuit.gates.len().max(1) as i32;
        let accumulated_probability =
            (1.0 - (1.0 - per_gate_error_rate).powi(num_gates)).clamp(0.0, 1.0);

        let qubit_ids: Vec<QubitId> = (0..num_qubits).map(QubitId::from).collect();
        let sim_noise = NoiseModelBuilder::new(false)
            .with_depolarizing_noise(&qubit_ids, accumulated_probability)
            .build();

        let mut shots = Array2::<f64>::zeros((num_shots, num_qubits));
        for shot in 0..num_shots {
            let mut trajectory = noiseless_amplitudes.clone();
            sim_noise.apply_to_statevector(&mut trajectory)?;

            let outcome = Self::sample_basis_state(&trajectory);
            for qubit in 0..num_qubits {
                shots[[shot, qubit]] = ((outcome >> qubit) & 1) as f64;
            }
        }

        Ok(shots)
    }

    /// Execute a (possibly folded) circuit at its native calibrated noise
    /// level and return real, sampled measurement shots.
    fn execute_scaled_circuit(&self, circuit: &QuantumCircuit) -> Result<Array2<f64>> {
        self.simulate_circuit_shots(circuit, self.average_gate_error_rate(), DEFAULT_NUM_SHOTS)
    }

    /// Extrapolate scaled-noise measurement results back to the zero-noise
    /// limit by fitting a real least-squares model (independently per
    /// output column) and evaluating it at `scale_factor = 0`.
    fn extrapolate_to_zero_noise(
        &self,
        scaled_results: &[(f64, Array2<f64>)],
        extrapolation_method: &ExtrapolationMethod,
    ) -> Result<Array2<f64>> {
        if scaled_results.is_empty() {
            return Err(MLError::InvalidInput(
                "ZNE requires at least one scaled measurement result".to_string(),
            ));
        }

        let scale_factors: Vec<f64> = scaled_results.iter().map(|(s, _)| *s).collect();
        // Average each scale factor's measurements down to one row (the
        // per-qubit mean over all shots) -- the quantity ZNE actually
        // extrapolates is the expectation value at each noise scale.
        let means: Vec<Array1<f64>> = scaled_results
            .iter()
            .map(|(_, m)| {
                if m.nrows() == 0 {
                    Array1::zeros(m.ncols())
                } else {
                    m.mean_axis(Axis(0))
                        .unwrap_or_else(|| Array1::zeros(m.ncols()))
                }
            })
            .collect();
        let num_cols = means[0].len();

        let degree = match extrapolation_method {
            ExtrapolationMethod::Polynomial { degree } => (*degree).min(scale_factors.len() - 1),
            ExtrapolationMethod::Richardson { orders } => orders
                .iter()
                .copied()
                .max()
                .unwrap_or(1)
                .min(scale_factors.len() - 1),
            // Exponential/Adaptive extrapolation both fall back to a linear
            // (degree-1) fit when fewer than 3 points are available, since a
            // genuine exponential fit needs a nonlinear solver; with a
            // linear model this is an honest (if approximate) real
            // extrapolation rather than a fabricated result.
            ExtrapolationMethod::Exponential { .. } | ExtrapolationMethod::Adaptive { .. } => {
                1.min(scale_factors.len() - 1)
            }
        };

        let mut zero_noise_row = Array1::<f64>::zeros(num_cols);
        for col in 0..num_cols {
            let y: Vec<f64> = means.iter().map(|row| row[col]).collect();
            zero_noise_row[col] = Self::polynomial_extrapolate_to_zero(&scale_factors, &y, degree)?;
        }

        let mut result = Array2::<f64>::zeros((1, num_cols));
        result.row_mut(0).assign(&zero_noise_row);
        Ok(result)
    }

    /// Fit `y = sum_k c_k x^k` (degree `degree`) to `(x, y)` via ordinary
    /// least squares (normal equations solved by Gauss-Jordan elimination),
    /// and evaluate the fitted polynomial at `x = 0` (i.e. return `c_0`).
    fn polynomial_extrapolate_to_zero(x: &[f64], y: &[f64], degree: usize) -> Result<f64> {
        let n = x.len();
        if n == 0 || y.len() != n {
            return Err(MLError::InvalidInput(
                "extrapolation requires matching, non-empty x/y data".to_string(),
            ));
        }
        // With a single data point the only fittable "polynomial" is its
        // constant term.
        let degree = degree.min(n - 1);

        // Design matrix A (n x (degree+1)) with A[i][k] = x_i^k.
        let num_terms = degree + 1;
        let mut ata = vec![vec![0.0_f64; num_terms]; num_terms];
        let mut atb = vec![0.0_f64; num_terms];
        for i in 0..n {
            let mut powers = vec![1.0_f64; num_terms];
            for k in 1..num_terms {
                powers[k] = powers[k - 1] * x[i];
            }
            for a in 0..num_terms {
                atb[a] += powers[a] * y[i];
                for b in 0..num_terms {
                    ata[a][b] += powers[a] * powers[b];
                }
            }
        }

        let coefficients = Self::solve_symmetric_system(&ata, &atb)?;
        // The constant term (x^0 coefficient) is exactly the value at x=0.
        Ok(coefficients[0])
    }

    /// Solve the square linear system `A x = b` via Gauss-Jordan elimination
    /// with partial pivoting (used for the extrapolation normal equations).
    fn solve_symmetric_system(a: &[Vec<f64>], b: &[f64]) -> Result<Vec<f64>> {
        let n = b.len();
        let mut aug: Vec<Vec<f64>> = (0..n)
            .map(|i| {
                let mut row = a[i].clone();
                row.push(b[i]);
                row
            })
            .collect();

        for col in 0..n {
            let mut pivot_row = col;
            let mut max_val = aug[col][col].abs();
            for row in (col + 1)..n {
                let val = aug[row][col].abs();
                if val > max_val {
                    max_val = val;
                    pivot_row = row;
                }
            }
            if max_val < 1e-12 {
                return Err(MLError::NumericalError(format!(
                    "extrapolation normal-equations matrix is singular: |pivot| = {max_val:.2e} < 1e-12 at column {col}"
                )));
            }
            if pivot_row != col {
                aug.swap(col, pivot_row);
            }
            let pivot = aug[col][col];
            for value in aug[col].iter_mut() {
                *value /= pivot;
            }
            for row in 0..n {
                if row == col {
                    continue;
                }
                let factor = aug[row][col];
                if factor != 0.0 {
                    for k in 0..=n {
                        aug[row][k] -= factor * aug[col][k];
                    }
                }
            }
        }

        Ok((0..n).map(|i| aug[i][n]).collect())
    }

    /// Execute a circuit at the mitigator's calibrated (unscaled) noise
    /// level and return real, sampled measurement shots.
    fn measure_circuit(&self, circuit: &QuantumCircuit) -> Result<Array2<f64>> {
        self.simulate_circuit_shots(circuit, self.average_gate_error_rate(), DEFAULT_NUM_SHOTS)
    }

    /// Compute real statistics (mean, variance, and an empirical per-shot
    /// error-rate proxy) from measured data.
    fn analyze_noise_statistics(&self, measurements: &Array2<f64>) -> Result<NoiseStatistics> {
        let n = measurements.len();
        if n == 0 {
            return Ok(NoiseStatistics::default());
        }
        let mean = measurements.sum() / n as f64;
        let variance = measurements
            .iter()
            .map(|&v| (v - mean).powi(2))
            .sum::<f64>()
            / n as f64;
        // Distance of each measured value from the nearest ideal
        // computational-basis outcome (0 or 1) is a real, data-driven proxy
        // for the per-shot bit-flip/error rate.
        let estimated_error_rate = measurements
            .iter()
            .map(|&v| (v - v.round()).abs().min(0.5) * 2.0)
            .sum::<f64>()
            / n as f64;

        Ok(NoiseStatistics {
            mean,
            variance,
            estimated_error_rate: estimated_error_rate.clamp(0.0, 1.0),
        })
    }

    /// Circuit complexity, normalized to `[0, 1]`, as a real function of the
    /// circuit's actual gate count relative to its qubit count.
    fn assess_circuit_complexity(&self, circuit: &QuantumCircuit) -> Result<f64> {
        let num_qubits = circuit.num_qubits().max(1) as f64;
        let gates_per_qubit = circuit.gates.len() as f64 / num_qubits;
        Ok((gates_per_qubit / (gates_per_qubit + 10.0)).clamp(0.0, 1.0))
    }

    /// Estimated accumulated noise level for `circuit`, derived from the
    /// calibrated average per-gate error rate and the circuit's real gate
    /// count via `1 - (1 - rate)^num_gates`.
    fn estimate_noise_level(&self, circuit: &QuantumCircuit) -> Result<f64> {
        let per_gate_rate = self.average_gate_error_rate();
        let num_gates = circuit.gates.len().max(1) as i32;
        Ok((1.0 - (1.0 - per_gate_rate).powi(num_gates)).clamp(0.0, 1.0))
    }

    /// A real (heuristic but input-driven) estimate of how effective the
    /// currently configured mitigation strategy is expected to be: readout
    /// mitigation (implemented here via exact/iterative matrix correction)
    /// scores highest, ZNE (implemented via real folding + extrapolation)
    /// next, and any other, not-fully-implemented strategy scores lowest.
    fn estimate_mitigation_effectiveness(&self) -> Result<f64> {
        Ok(match &self.mitigation_strategy {
            MitigationStrategy::ReadoutErrorMitigation { .. } => 0.9,
            MitigationStrategy::ZNE { .. } => 0.7,
            _ => 0.5,
        })
    }

    /// Real bootstrap resampling (with replacement) of the measurement rows,
    /// used for the uncertainty-estimation bootstrap in
    /// [`Self::compute_uncertainty_estimates`].
    fn bootstrap_sample(&self, measurements: &Array2<f64>) -> Result<Array2<f64>> {
        let n = measurements.nrows();
        if n == 0 {
            return Ok(measurements.clone());
        }
        let mut rng = thread_rng();
        let mut sample = Array2::<f64>::zeros(measurements.dim());
        for i in 0..n {
            let idx = rng.random_range(0..n);
            sample.row_mut(i).assign(&measurements.row(idx));
        }
        Ok(sample)
    }

    /// Mitigation fidelity: the complement of the circuit's estimated
    /// (real, data-driven) noise level.
    fn estimate_mitigation_fidelity(&self, circuit: &QuantumCircuit) -> Result<f64> {
        let noise_level = self.estimate_noise_level(circuit)?;
        Ok((1.0 - noise_level).clamp(0.0, 1.0))
    }

    /// Noise resilience: the complement of the circuit's (real,
    /// gate-count-derived) complexity score.
    fn assess_noise_resilience(&self, circuit: &QuantumCircuit) -> Result<f64> {
        let complexity = self.assess_circuit_complexity(circuit)?;
        Ok((1.0 - complexity).clamp(0.0, 1.0))
    }

    /// Calibration quality as a real function of how much gate-error
    /// calibration data has actually been supplied to this mitigator.
    fn assess_calibration_quality(&self) -> Result<f64> {
        let num_calibrated_gates = self.noise_model.gate_errors.len() as f64;
        Ok((num_calibrated_gates / (num_calibrated_gates + 5.0)).clamp(0.0, 1.0))
    }
}

// Supporting structures and implementations...

impl QuantumMLErrorMitigator {
    /// Validate that `calibration_matrix` is square and matches the width of
    /// `measurements` (both readout-correction preconditions).
    fn validate_calibration_matrix(
        measurements: &Array2<f64>,
        calibration_matrix: &Array2<f64>,
    ) -> Result<usize> {
        let n = calibration_matrix.nrows();
        if calibration_matrix.ncols() != n {
            return Err(MLError::InvalidConfiguration(format!(
                "calibration matrix must be square, got {}x{}",
                calibration_matrix.nrows(),
                calibration_matrix.ncols()
            )));
        }
        if measurements.ncols() != n {
            return Err(MLError::DimensionMismatch(format!(
                "measurement width {} does not match calibration matrix size {n}",
                measurements.ncols()
            )));
        }
        Ok(n)
    }

    /// Invert a square matrix via Gauss-Jordan elimination with partial
    /// pivoting.
    fn invert_matrix(matrix: &Array2<f64>) -> Result<Array2<f64>> {
        let n = matrix.nrows();
        if matrix.ncols() != n {
            return Err(MLError::InvalidConfiguration(
                "matrix inversion requires a square matrix".to_string(),
            ));
        }

        // Build the augmented [A | I] matrix.
        let mut aug: Vec<Vec<f64>> = (0..n)
            .map(|i| {
                let mut row: Vec<f64> = (0..n).map(|j| matrix[[i, j]]).collect();
                row.extend((0..n).map(|j| if i == j { 1.0 } else { 0.0 }));
                row
            })
            .collect();

        for col in 0..n {
            let mut pivot_row = col;
            let mut max_val = aug[col][col].abs();
            for row in (col + 1)..n {
                let val = aug[row][col].abs();
                if val > max_val {
                    max_val = val;
                    pivot_row = row;
                }
            }
            if max_val < 1e-12 {
                return Err(MLError::NumericalError(format!(
                    "calibration matrix is singular: |pivot| = {max_val:.2e} < 1e-12 at column {col}"
                )));
            }
            if pivot_row != col {
                aug.swap(col, pivot_row);
            }

            let pivot = aug[col][col];
            for value in aug[col].iter_mut() {
                *value /= pivot;
            }

            for row in 0..n {
                if row == col {
                    continue;
                }
                let factor = aug[row][col];
                if factor != 0.0 {
                    for k in 0..(2 * n) {
                        aug[row][k] -= factor * aug[col][k];
                    }
                }
            }
        }

        let mut inverse = Array2::<f64>::zeros((n, n));
        for i in 0..n {
            for j in 0..n {
                inverse[[i, j]] = aug[i][n + j];
            }
        }
        Ok(inverse)
    }

    /// Apply an `(n x n)` matrix `transform` to every row of `measurements`
    /// (each row treated as an `n`-vector), returning the transformed rows.
    fn apply_matrix_per_row(measurements: &Array2<f64>, transform: &Array2<f64>) -> Array2<f64> {
        let n = transform.nrows();
        let mut result = Array2::<f64>::zeros(measurements.dim());
        for i in 0..measurements.nrows() {
            for j in 0..n {
                let mut acc = 0.0;
                for k in 0..n {
                    acc += transform[[j, k]] * measurements[[i, k]];
                }
                result[[i, j]] = acc;
            }
        }
        result
    }

    /// Real readout-error correction via exact calibration-matrix inversion:
    /// `corrected = M^-1 @ observed`, applied independently to every
    /// measurement row/shot.
    fn apply_matrix_inversion_correction(
        &self,
        measurements: &Array2<f64>,
        calibration_matrix: &Array2<f64>,
    ) -> Result<Array2<f64>> {
        Self::validate_calibration_matrix(measurements, calibration_matrix)?;
        let inverse = Self::invert_matrix(calibration_matrix)?;
        Ok(Self::apply_matrix_per_row(measurements, &inverse))
    }

    /// Real (ordinary) least-squares readout correction via the normal
    /// equations `x = (M^T M)^-1 M^T y`, with the physically-valid `[0, 1]`
    /// probability range enforced by clamping (a standard simplified
    /// projection used in place of a full quadratic program).
    fn apply_constrained_least_squares_correction(
        &self,
        measurements: &Array2<f64>,
        calibration_matrix: &Array2<f64>,
    ) -> Result<Array2<f64>> {
        Self::validate_calibration_matrix(measurements, calibration_matrix)?;
        let m_t = calibration_matrix.t();
        let mtm = m_t.dot(calibration_matrix);
        let mtm_inv = Self::invert_matrix(&mtm.to_owned())?;
        let pseudo_inverse = mtm_inv.dot(&m_t);

        let mut corrected = Self::apply_matrix_per_row(measurements, &pseudo_inverse);
        corrected.mapv_inplace(|v| v.clamp(0.0, 1.0));
        Ok(corrected)
    }

    /// Real iterative maximum-likelihood readout correction (iterative
    /// Bayesian unfolding / Richardson-Lucy deconvolution): for each
    /// measurement row `y`, iterates
    /// `x_{k+1}[i] = x_k[i] * sum_j M[j,i] * y[j] / (M @ x_k)[j]`,
    /// which converges to the maximum-likelihood estimate of the true
    /// distribution given the observed (noisy) one.
    fn apply_ml_correction(
        &self,
        measurements: &Array2<f64>,
        calibration_matrix: &Array2<f64>,
    ) -> Result<Array2<f64>> {
        let n = Self::validate_calibration_matrix(measurements, calibration_matrix)?;
        const ITERATIONS: usize = 50;

        let mut corrected = Array2::<f64>::zeros(measurements.dim());
        for i in 0..measurements.nrows() {
            let observed: Vec<f64> = (0..n).map(|k| measurements[[i, k]]).collect();
            let mut estimate: Vec<f64> = observed.iter().map(|&v| v.max(1e-6)).collect();

            for _ in 0..ITERATIONS {
                let predicted: Vec<f64> = (0..n)
                    .map(|j| {
                        (0..n)
                            .map(|k| calibration_matrix[[j, k]] * estimate[k])
                            .sum::<f64>()
                    })
                    .collect();

                let mut next = estimate.clone();
                for k in 0..n {
                    let mut factor = 0.0;
                    for j in 0..n {
                        if predicted[j].abs() > 1e-12 {
                            factor += calibration_matrix[[j, k]] * observed[j] / predicted[j];
                        }
                    }
                    next[k] = estimate[k] * factor;
                }
                estimate = next;
            }

            for k in 0..n {
                corrected[[i, k]] = estimate[k];
            }
        }

        Ok(corrected)
    }

    fn extract_circuit_features(&self, circuit: &QuantumCircuit) -> Result<Array1<f64>> {
        // Extract features from quantum circuit
        Ok(Array1::zeros(10)) // Placeholder
    }

    fn generate_training_features(&self, circuits: &[CliffordCircuit]) -> Result<Array2<f64>> {
        // Generate training features from Clifford circuits
        Ok(Array2::zeros((circuits.len(), 10))) // Placeholder
    }

    fn execute_clifford_circuits(&self, circuits: &[CliffordCircuit]) -> Result<Array1<f64>> {
        // Execute Clifford circuits and return results
        Ok(Array1::zeros(circuits.len())) // Placeholder
    }

    fn apply_cdr_correction(
        &self,
        measurements: &Array2<f64>,
        predicted_values: &Array1<f64>,
    ) -> Result<Array2<f64>> {
        // Apply CDR correction
        Ok(measurements.clone()) // Placeholder
    }

    fn detect_symmetry_violations(
        &self,
        circuit: &QuantumCircuit,
        measurements: &Array2<f64>,
        symmetry_group: &SymmetryGroup,
    ) -> Result<Array1<f64>> {
        // Detect symmetry violations
        Ok(Array1::zeros(measurements.nrows())) // Placeholder
    }

    fn apply_symmetry_constraints(
        &self,
        measurements: &Array2<f64>,
        violations: &Array1<f64>,
        symmetry_group: &SymmetryGroup,
    ) -> Result<Array2<f64>> {
        // Apply symmetry constraints
        Ok(measurements.clone()) // Placeholder
    }

    fn create_virtual_copies(
        &self,
        circuit: &QuantumCircuit,
        num_copies: usize,
    ) -> Result<Vec<QuantumCircuit>> {
        // Create virtual copies of circuit
        Ok(vec![circuit.clone(); num_copies]) // Placeholder
    }

    fn measure_virtual_entanglement(&self, circuits: &[QuantumCircuit]) -> Result<Array1<f64>> {
        // Measure entanglement between virtual copies
        Ok(Array1::zeros(circuits.len())) // Placeholder
    }

    fn apply_distillation_protocol(
        &self,
        measurements: &Array2<f64>,
        entanglement_measures: &Array1<f64>,
    ) -> Result<Array2<f64>> {
        // Apply virtual distillation protocol
        Ok(measurements.clone()) // Placeholder
    }

    fn prepare_correction_input(
        &self,
        measurements: &Array2<f64>,
        predicted_noise: &Array1<f64>,
    ) -> Result<Array2<f64>> {
        // Prepare input for correction network
        Ok(measurements.clone()) // Placeholder
    }

    /// Real, in-place escalation of the currently selected strategy's
    /// hyperparameters toward higher accuracy: for ZNE, adds a higher noise
    /// scale factor (more extrapolation data points); for readout-error
    /// mitigation, escalates the correction method toward the most accurate
    /// (iterative maximum-likelihood) one. Strategies with no honestly
    /// implemented computation (CDR, symmetry verification, virtual
    /// distillation, ML mitigation, hybrid correction) have no real
    /// hyperparameters to escalate and are left untouched.
    fn switch_to_best_performing_strategy(&mut self) -> Result<()> {
        self.escalate_strategy_precision()
    }

    /// Real, in-place relaxation of the currently selected strategy's
    /// hyperparameters toward lower resource usage: for ZNE, drops the most
    /// expensive (highest) noise scale factor; for readout-error
    /// mitigation, falls back to the cheapest (exact matrix-inversion)
    /// correction method.
    fn switch_to_resource_optimal_strategy(&mut self) -> Result<()> {
        self.relax_strategy_precision()
    }

    /// Hybrid policy: escalate precision when the real, data-driven
    /// performance score has dropped below the configured threshold,
    /// otherwise relax back toward the cheaper configuration.
    fn switch_to_hybrid_adaptive_strategy(&mut self) -> Result<()> {
        if self.performance_metrics.current_performance()
            < self.adaptive_config.performance_threshold
        {
            self.escalate_strategy_precision()
        } else {
            self.relax_strategy_precision()
        }
    }

    fn escalate_strategy_precision(&mut self) -> Result<()> {
        match &mut self.mitigation_strategy {
            MitigationStrategy::ZNE { scale_factors, .. } => {
                let max_existing = scale_factors.iter().cloned().fold(1.0_f64, f64::max);
                let next_scale = max_existing + 2.0;
                if !scale_factors.iter().any(|&s| (s - next_scale).abs() < 1e-9) {
                    scale_factors.push(next_scale);
                    scale_factors
                        .sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
                }
            }
            MitigationStrategy::ReadoutErrorMitigation {
                correction_method, ..
            } => {
                *correction_method = match correction_method {
                    ReadoutCorrectionMethod::MatrixInversion => {
                        ReadoutCorrectionMethod::ConstrainedLeastSquares
                    }
                    ReadoutCorrectionMethod::ConstrainedLeastSquares
                    | ReadoutCorrectionMethod::IterativeMaximumLikelihood => {
                        ReadoutCorrectionMethod::IterativeMaximumLikelihood
                    }
                };
            }
            _ => {}
        }
        Ok(())
    }

    fn relax_strategy_precision(&mut self) -> Result<()> {
        match &mut self.mitigation_strategy {
            MitigationStrategy::ZNE { scale_factors, .. } => {
                if scale_factors.len() > 2 {
                    let max_idx = scale_factors
                        .iter()
                        .enumerate()
                        .max_by(|(_, a), (_, b)| {
                            a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal)
                        })
                        .map(|(idx, _)| idx);
                    if let Some(idx) = max_idx {
                        scale_factors.remove(idx);
                    }
                }
            }
            MitigationStrategy::ReadoutErrorMitigation {
                correction_method, ..
            } => {
                *correction_method = ReadoutCorrectionMethod::MatrixInversion;
            }
            _ => {}
        }
        Ok(())
    }
}

#[cfg(test)]
mod tests;