quantrs2_sim/

fault_tolerant_synthesis.rs

//! Fault-Tolerant Gate Synthesis with Logical Operations
//!
//! This module implements fault-tolerant quantum computation by synthesizing logical gates
//! using quantum error correction codes. It provides tools for converting arbitrary logical
//! operations into fault-tolerant implementations using various error correction codes like
//! surface codes, color codes, and topological codes.
//!
//! Key features:
//! - Logical gate synthesis for various quantum error correction codes
//! - Fault-tolerant gate decomposition with minimal resource overhead
//! - Magic state distillation for non-Clifford gates
//! - Surface code compilation with optimal routing
//! - Topological quantum computation synthesis
//! - Resource estimation for fault-tolerant circuits
//! - Adaptive code distance selection
//! - Logical measurement and state preparation protocols

use scirs2_core::ndarray::{Array1, Array2};
use std::collections::HashMap;

use crate::circuit_interfaces::{InterfaceCircuit, InterfaceGate, InterfaceGateType};
use crate::error::{Result, SimulatorError};

/// Fault-tolerant synthesis configuration
#[derive(Debug, Clone)]
pub struct FaultTolerantConfig {
    /// Target logical error rate
    pub target_logical_error_rate: f64,
    /// Physical error rate of the hardware
    pub physical_error_rate: f64,
    /// Error correction code to use
    pub error_correction_code: ErrorCorrectionCode,
    /// Code distance
    pub code_distance: usize,
    /// Enable magic state distillation
    pub enable_magic_state_distillation: bool,
    /// Enable adaptive code distance
    pub enable_adaptive_distance: bool,
    /// Resource optimization level
    pub optimization_level: FTOptimizationLevel,
    /// Maximum synthesis depth
    pub max_synthesis_depth: usize,
    /// Parallelization threshold
    pub parallel_threshold: usize,
}

impl Default for FaultTolerantConfig {
    fn default() -> Self {
        Self {
            target_logical_error_rate: 1e-6,
            physical_error_rate: 1e-3,
            error_correction_code: ErrorCorrectionCode::SurfaceCode,
            code_distance: 5,
            enable_magic_state_distillation: true,
            enable_adaptive_distance: true,
            optimization_level: FTOptimizationLevel::Balanced,
            max_synthesis_depth: 1000,
            parallel_threshold: 100,
        }
    }
}

/// Error correction codes for fault-tolerant synthesis
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum ErrorCorrectionCode {
    /// Surface code (2D topological)
    SurfaceCode,
    /// Color code (2D topological)
    ColorCode,
    /// Steane code (7,1,3)
    SteaneCode,
    /// Shor code (9,1,3)
    ShorCode,
    /// Reed-Muller code
    ReedMullerCode,
    /// Bacon-Shor code
    BaconShorCode,
    /// Subsystem surface code
    SubsystemSurfaceCode,
}

/// Fault-tolerant optimization levels
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FTOptimizationLevel {
    /// Minimize resource usage
    Space,
    /// Minimize computation time
    Time,
    /// Balance space and time
    Balanced,
    /// Minimize logical error rate
    ErrorRate,
    /// Custom optimization
    Custom,
}

/// Logical gate types for fault-tolerant synthesis
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum LogicalGateType {
    /// Logical Pauli-X
    LogicalX,
    /// Logical Pauli-Y
    LogicalY,
    /// Logical Pauli-Z
    LogicalZ,
    /// Logical Hadamard
    LogicalH,
    /// Logical S gate
    LogicalS,
    /// Logical T gate (requires magic states)
    LogicalT,
    /// Logical CNOT
    LogicalCNOT,
    /// Logical CZ
    LogicalCZ,
    /// Logical Toffoli (requires magic states)
    LogicalToffoli,
    /// Logical rotation (parametric)
    LogicalRotation(f64),
    /// Logical measurement
    LogicalMeasurement,
    /// Logical state preparation
    LogicalPreparation,
}

impl Eq for LogicalGateType {}

impl std::hash::Hash for LogicalGateType {
    fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
        std::mem::discriminant(self).hash(state);
        if let Self::LogicalRotation(angle) = self {
            // Convert float to bits for consistent hashing
            angle.to_bits().hash(state);
        }
    }
}

/// Logical gate implementation
#[derive(Debug, Clone)]
pub struct LogicalGate {
    /// Gate type
    pub gate_type: LogicalGateType,
    /// Target logical qubits
    pub logical_qubits: Vec<usize>,
    /// Physical implementation
    pub physical_implementation: InterfaceCircuit,
    /// Resource requirements
    pub resources: ResourceRequirements,
    /// Error rate estimate
    pub error_rate: f64,
}

/// Resource requirements for fault-tolerant operations
#[derive(Debug, Clone, Default)]
pub struct ResourceRequirements {
    /// Number of physical qubits
    pub physical_qubits: usize,
    /// Number of physical gates
    pub physical_gates: usize,
    /// Number of measurement rounds
    pub measurement_rounds: usize,
    /// Magic states required
    pub magic_states: usize,
    /// Computation time (in time steps)
    pub time_steps: usize,
    /// Memory requirements (ancilla qubits)
    pub ancilla_qubits: usize,
}

/// Fault-tolerant synthesis result
#[derive(Debug, Clone)]
pub struct FaultTolerantSynthesisResult {
    /// Synthesized fault-tolerant circuit
    pub fault_tolerant_circuit: InterfaceCircuit,
    /// Logical error rate achieved
    pub logical_error_rate: f64,
    /// Resource usage
    pub resources: ResourceRequirements,
    /// Synthesis statistics
    pub synthesis_stats: SynthesisStatistics,
    /// Code distance used
    pub code_distance: usize,
    /// Error correction overhead
    pub overhead_factor: f64,
}

/// Synthesis statistics
#[derive(Debug, Clone, Default)]
pub struct SynthesisStatistics {
    /// Number of logical gates synthesized
    pub logical_gates_synthesized: usize,
    /// Average gate synthesis time
    pub avg_synthesis_time_ms: f64,
    /// Total synthesis time
    pub total_synthesis_time_ms: f64,
    /// Magic state consumption
    pub magic_states_consumed: usize,
    /// Code distance adaptations
    pub distance_adaptations: usize,
    /// Optimization passes
    pub optimization_passes: usize,
}

/// Magic state types for non-Clifford gates
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum MagicStateType {
    /// T-state for T gate implementation
    TState,
    /// Y-state for Y-rotation
    YState,
    /// CCZ-state for Toffoli implementation
    CCZState,
    /// Custom magic state
    Custom(usize),
}

/// Magic state distillation protocol
#[derive(Debug, Clone)]
pub struct MagicStateProtocol {
    /// Input magic state type
    pub input_state: MagicStateType,
    /// Output magic state type
    pub output_state: MagicStateType,
    /// Distillation circuit
    pub distillation_circuit: InterfaceCircuit,
    /// Error reduction factor
    pub error_reduction: f64,
    /// Resource overhead
    pub overhead: usize,
}

/// Surface code implementation for fault-tolerant synthesis
#[derive(Debug, Clone)]
pub struct SurfaceCodeSynthesizer {
    /// Code distance
    pub distance: usize,
    /// Surface code layout
    pub layout: SurfaceCodeLayout,
    /// Stabilizer generators
    pub stabilizers: Vec<Array1<i8>>,
    /// Logical operators
    pub logical_operators: HashMap<LogicalGateType, Array2<i8>>,
    /// Error correction schedule
    pub error_correction_schedule: Vec<ErrorCorrectionRound>,
}

/// Surface code layout
#[derive(Debug, Clone)]
pub struct SurfaceCodeLayout {
    /// Data qubit positions
    pub data_qubits: Array2<usize>,
    /// X-stabilizer positions
    pub x_stabilizers: Array2<usize>,
    /// Z-stabilizer positions
    pub z_stabilizers: Array2<usize>,
    /// Boundary conditions
    pub boundaries: BoundaryConditions,
}

/// Boundary conditions for surface codes
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum BoundaryConditions {
    /// Open boundaries
    Open,
    /// Periodic boundaries
    Periodic,
    /// Twisted boundaries
    Twisted,
    /// Rough-smooth boundaries
    RoughSmooth,
}

/// Error correction round
#[derive(Debug, Clone)]
pub struct ErrorCorrectionRound {
    /// Stabilizer measurements
    pub stabilizer_measurements: Vec<StabilizerMeasurement>,
    /// Syndrome extraction
    pub syndrome_extraction: InterfaceCircuit,
    /// Error correction
    pub error_correction: InterfaceCircuit,
    /// Round duration
    pub duration: usize,
}

/// Stabilizer measurement
#[derive(Debug, Clone)]
pub struct StabilizerMeasurement {
    /// Stabilizer index
    pub stabilizer_index: usize,
    /// Measurement circuit
    pub measurement_circuit: InterfaceCircuit,
    /// Syndrome qubit
    pub syndrome_qubit: usize,
    /// Data qubits involved
    pub data_qubits: Vec<usize>,
}

/// Main fault-tolerant gate synthesizer
pub struct FaultTolerantSynthesizer {
    /// Configuration
    config: FaultTolerantConfig,
    /// Surface code synthesizer
    surface_code: Option<SurfaceCodeSynthesizer>,
    /// Magic state protocols
    magic_state_protocols: HashMap<LogicalGateType, MagicStateProtocol>,
    /// Logical gate library
    gate_library: HashMap<LogicalGateType, LogicalGate>,
    /// Resource estimator
    resource_estimator: ResourceEstimator,
    /// Synthesis cache
    synthesis_cache: HashMap<String, FaultTolerantSynthesisResult>,
}

/// Resource estimator for fault-tolerant circuits
#[derive(Debug, Clone, Default)]
pub struct ResourceEstimator {
    /// Physical error model
    pub error_model: PhysicalErrorModel,
    /// Code parameters
    pub code_parameters: HashMap<ErrorCorrectionCode, CodeParameters>,
    /// Magic state costs
    pub magic_state_costs: HashMap<MagicStateType, usize>,
}

/// Physical error model
#[derive(Debug, Clone, Default)]
pub struct PhysicalErrorModel {
    /// Gate error rates
    pub gate_errors: HashMap<String, f64>,
    /// Measurement error rate
    pub measurement_error: f64,
    /// Memory error rate (per time step)
    pub memory_error: f64,
    /// Correlated error probability
    pub correlated_error: f64,
}

/// Code parameters for different error correction codes
#[derive(Debug, Clone, Default)]
pub struct CodeParameters {
    /// Encoding rate (k/n)
    pub encoding_rate: f64,
    /// Threshold error rate
    pub threshold: f64,
    /// Resource scaling
    pub resource_scaling: f64,
    /// Logical error suppression
    pub error_suppression: f64,
}

impl FaultTolerantSynthesizer {
    /// Create new fault-tolerant synthesizer
    pub fn new(config: FaultTolerantConfig) -> Result<Self> {
        let mut synthesizer = Self {
            config: config.clone(),
            surface_code: None,
            magic_state_protocols: HashMap::new(),
            gate_library: HashMap::new(),
            resource_estimator: ResourceEstimator::default(),
            synthesis_cache: HashMap::new(),
        };

        // Initialize based on error correction code
        if config.error_correction_code == ErrorCorrectionCode::SurfaceCode {
            synthesizer.surface_code = Some(synthesizer.create_surface_code()?);
        } else {
            // Initialize other codes as needed
        }

        // Initialize magic state protocols
        synthesizer.initialize_magic_state_protocols()?;

        // Initialize gate library
        synthesizer.initialize_gate_library()?;

        // Initialize resource estimator
        synthesizer.initialize_resource_estimator()?;

        Ok(synthesizer)
    }

    /// Synthesize fault-tolerant implementation of a logical circuit
    pub fn synthesize_logical_circuit(
        &mut self,
        logical_circuit: &InterfaceCircuit,
    ) -> Result<FaultTolerantSynthesisResult> {
        let start_time = std::time::Instant::now();

        // Check cache first
        let cache_key = self.generate_cache_key(logical_circuit);
        if let Some(cached_result) = self.synthesis_cache.get(&cache_key) {
            return Ok(cached_result.clone());
        }

        // Adapt code distance if enabled
        let optimal_distance = if self.config.enable_adaptive_distance {
            self.calculate_optimal_distance(logical_circuit)?
        } else {
            self.config.code_distance
        };

        // Initialize synthesis result
        let mut result = FaultTolerantSynthesisResult {
            fault_tolerant_circuit: InterfaceCircuit::new(0, 0),
            logical_error_rate: 0.0,
            resources: ResourceRequirements::default(),
            synthesis_stats: SynthesisStatistics::default(),
            code_distance: optimal_distance,
            overhead_factor: 0.0,
        };

        // Synthesize each logical gate
        for gate in &logical_circuit.gates {
            let logical_gate_type = self.map_interface_gate_to_logical(gate)?;
            let synthesized_gate = self.synthesize_logical_gate(logical_gate_type, &gate.qubits)?;

            // Add to fault-tolerant circuit
            self.append_synthesized_gate(&mut result.fault_tolerant_circuit, &synthesized_gate)?;

            // Update resource requirements
            self.update_resources(&mut result.resources, &synthesized_gate.resources);

            result.synthesis_stats.logical_gates_synthesized += 1;
        }

        // Add error correction rounds
        self.add_error_correction_rounds(&mut result.fault_tolerant_circuit, optimal_distance)?;

        // Calculate logical error rate
        result.logical_error_rate = self.calculate_logical_error_rate(&result)?;

        // Calculate overhead factor
        result.overhead_factor =
            result.resources.physical_qubits as f64 / logical_circuit.num_qubits as f64;

        // Update synthesis statistics
        result.synthesis_stats.total_synthesis_time_ms = start_time.elapsed().as_millis() as f64;
        result.synthesis_stats.avg_synthesis_time_ms =
            result.synthesis_stats.total_synthesis_time_ms
                / result.synthesis_stats.logical_gates_synthesized as f64;

        // Cache result
        self.synthesis_cache.insert(cache_key, result.clone());

        Ok(result)
    }

    /// Synthesize a single logical gate
    pub fn synthesize_logical_gate(
        &mut self,
        gate_type: LogicalGateType,
        logical_qubits: &[usize],
    ) -> Result<LogicalGate> {
        // Check if gate is in library
        if let Some(template) = self.gate_library.get(&gate_type) {
            let mut synthesized = template.clone();
            synthesized.logical_qubits = logical_qubits.to_vec();
            return Ok(synthesized);
        }

        // Synthesize gate based on type
        match gate_type {
            LogicalGateType::LogicalX | LogicalGateType::LogicalY | LogicalGateType::LogicalZ => {
                self.synthesize_logical_pauli(gate_type, logical_qubits)
            }
            LogicalGateType::LogicalH => self.synthesize_logical_hadamard(logical_qubits),
            LogicalGateType::LogicalS => self.synthesize_logical_s(logical_qubits),
            LogicalGateType::LogicalT => {
                self.synthesize_logical_t_with_magic_states(logical_qubits)
            }
            LogicalGateType::LogicalCNOT => self.synthesize_logical_cnot(logical_qubits),
            LogicalGateType::LogicalToffoli => {
                self.synthesize_logical_toffoli_with_magic_states(logical_qubits)
            }
            LogicalGateType::LogicalRotation(angle) => {
                self.synthesize_logical_rotation(logical_qubits, angle)
            }
            _ => Err(SimulatorError::InvalidConfiguration(format!(
                "Unsupported logical gate type: {gate_type:?}"
            ))),
        }
    }

    /// Create surface code synthesizer
    fn create_surface_code(&self) -> Result<SurfaceCodeSynthesizer> {
        let distance = self.config.code_distance;

        // Create surface code layout
        let layout = self.create_surface_code_layout(distance)?;

        // Generate stabilizer generators
        let stabilizers = self.generate_surface_code_stabilizers(distance)?;

        // Create logical operators
        let logical_operators = self.create_logical_operators(distance)?;

        // Create temporary surface code to generate error correction schedule
        let temp_surface_code = SurfaceCodeSynthesizer {
            distance,
            layout: layout.clone(),
            stabilizers: stabilizers.clone(),
            logical_operators: logical_operators.clone(),
            error_correction_schedule: Vec::new(), // Will be filled below
        };

        // Create error correction schedule using the temporary surface code
        let error_correction_schedule =
            self.create_error_correction_schedule_with_surface_code(distance, &temp_surface_code)?;

        Ok(SurfaceCodeSynthesizer {
            distance,
            layout,
            stabilizers,
            logical_operators,
            error_correction_schedule,
        })
    }

    /// Create surface code layout
    pub fn create_surface_code_layout(&self, distance: usize) -> Result<SurfaceCodeLayout> {
        let size = 2 * distance - 1;

        // Initialize qubit arrays
        let mut data_qubits = Array2::zeros((size, size));
        let mut x_stabilizers = Array2::zeros((distance - 1, distance));
        let mut z_stabilizers = Array2::zeros((distance, distance - 1));

        // Assign data qubit indices
        let mut qubit_index = 0;
        for i in 0..size {
            for j in 0..size {
                if (i + j) % 2 == 0 {
                    data_qubits[[i, j]] = qubit_index;
                    qubit_index += 1;
                }
            }
        }

        // Assign stabilizer indices
        for i in 0..distance - 1 {
            for j in 0..distance {
                x_stabilizers[[i, j]] = qubit_index;
                qubit_index += 1;
            }
        }

        for i in 0..distance {
            for j in 0..distance - 1 {
                z_stabilizers[[i, j]] = qubit_index;
                qubit_index += 1;
            }
        }

        Ok(SurfaceCodeLayout {
            data_qubits,
            x_stabilizers,
            z_stabilizers,
            boundaries: BoundaryConditions::Open,
        })
    }

    /// Generate stabilizer generators for surface code
    pub fn generate_surface_code_stabilizers(&self, distance: usize) -> Result<Vec<Array1<i8>>> {
        let mut stabilizers = Vec::new();
        let total_qubits = distance * distance;

        // X-type stabilizers
        for i in 0..distance - 1 {
            for j in 0..distance {
                let mut stabilizer = Array1::zeros(2 * total_qubits); // X and Z parts

                // Add X operations on neighboring data qubits
                let neighbors = self.get_x_stabilizer_neighbors(i, j, distance);
                for &qubit in &neighbors {
                    stabilizer[qubit] = 1; // X operation
                }

                stabilizers.push(stabilizer);
            }
        }

        // Z-type stabilizers
        for i in 0..distance {
            for j in 0..distance - 1 {
                let mut stabilizer = Array1::zeros(2 * total_qubits);

                // Add Z operations on neighboring data qubits
                let neighbors = self.get_z_stabilizer_neighbors(i, j, distance);
                for &qubit in &neighbors {
                    stabilizer[total_qubits + qubit] = 1; // Z operation
                }

                stabilizers.push(stabilizer);
            }
        }

        Ok(stabilizers)
    }

    /// Data-qubit flat index for grid position `(row, col)` on the
    /// `distance × distance` data lattice, or `None` if the position falls off
    /// the grid (used for boundary truncation).
    fn data_qubit_on_grid(row: isize, col: isize, distance: usize) -> Option<usize> {
        if row < 0 || col < 0 {
            return None;
        }
        let (r, c) = (row as usize, col as usize);
        if r >= distance || c >= distance {
            return None;
        }
        Some(r * distance + c)
    }

    /// Data qubits acted on by the X-type stabilizer at lattice site `(i, j)`.
    ///
    /// X stabilizers are modelled as **star (vertex) operators**: the operator
    /// centred on data-lattice vertex `(i, j)` acts on the four orthogonally
    /// adjacent data qubits — the up/down/left/right cross `{(i-1,j), (i+1,j),
    /// (i,j-1), (i,j+1)}` — the classic surface-code "plus" pattern. Neighbours
    /// that fall outside the `distance × distance` grid are truncated, giving
    /// weight-4 operators in the bulk and reduced-weight operators on the
    /// boundary. This is genuine lattice geometry (adjacent data qubits on the
    /// grid), replacing the former modular-arithmetic placeholder that wrapped
    /// around `distance²` with no geometric meaning.
    fn get_x_stabilizer_neighbors(&self, i: usize, j: usize, distance: usize) -> Vec<usize> {
        let (i, j) = (i as isize, j as isize);
        [(i - 1, j), (i + 1, j), (i, j - 1), (i, j + 1)]
            .into_iter()
            .filter_map(|(r, c)| Self::data_qubit_on_grid(r, c, distance))
            .collect()
    }

    /// Data qubits acted on by the Z-type stabilizer at lattice site `(i, j)`.
    ///
    /// Z stabilizers are modelled as **plaquette (face) operators**: the
    /// operator on the face whose top-left corner is data qubit `(i, j)` acts
    /// on the four corners of that unit cell — `{(i,j), (i,j+1), (i+1,j),
    /// (i+1,j+1)}` — the dual of the X star. Corners outside the grid are
    /// truncated for boundary faces, again replacing the meaningless modular
    /// placeholder with real grid adjacency.
    fn get_z_stabilizer_neighbors(&self, i: usize, j: usize, distance: usize) -> Vec<usize> {
        let (i, j) = (i as isize, j as isize);
        [(i, j), (i, j + 1), (i + 1, j), (i + 1, j + 1)]
            .into_iter()
            .filter_map(|(r, c)| Self::data_qubit_on_grid(r, c, distance))
            .collect()
    }

    /// Create logical operators
    fn create_logical_operators(
        &self,
        distance: usize,
    ) -> Result<HashMap<LogicalGateType, Array2<i8>>> {
        let mut logical_operators = HashMap::new();
        let total_qubits = distance * distance;

        // Logical X operator
        let mut logical_x = Array2::zeros((1, 2 * total_qubits));
        for i in 0..distance {
            logical_x[[0, i]] = 1; // X operation on first row
        }
        logical_operators.insert(LogicalGateType::LogicalX, logical_x);

        // Logical Z operator
        let mut logical_z = Array2::zeros((1, 2 * total_qubits));
        for i in 0..distance {
            logical_z[[0, total_qubits + i * distance]] = 1; // Z operation on first column
        }
        logical_operators.insert(LogicalGateType::LogicalZ, logical_z);

        Ok(logical_operators)
    }

    /// Create error correction schedule with provided surface code
    fn create_error_correction_schedule_with_surface_code(
        &self,
        distance: usize,
        surface_code: &SurfaceCodeSynthesizer,
    ) -> Result<Vec<ErrorCorrectionRound>> {
        let mut schedule = Vec::new();

        // Create syndrome extraction round
        let mut round = ErrorCorrectionRound {
            stabilizer_measurements: Vec::new(),
            syndrome_extraction: InterfaceCircuit::new(distance * distance + 100, 0), // Extra for ancillas
            error_correction: InterfaceCircuit::new(distance * distance, 0),
            duration: 1,
        };

        // Add stabilizer measurements
        for (i, stabilizer) in surface_code.stabilizers.iter().enumerate() {
            let measurement = StabilizerMeasurement {
                stabilizer_index: i,
                measurement_circuit: self.create_stabilizer_measurement_circuit(stabilizer)?,
                syndrome_qubit: distance * distance + i,
                data_qubits: self.get_stabilizer_data_qubits(stabilizer),
            };
            round.stabilizer_measurements.push(measurement);
        }

        schedule.push(round);
        Ok(schedule)
    }

    /// Create error correction schedule
    fn create_error_correction_schedule(
        &self,
        distance: usize,
    ) -> Result<Vec<ErrorCorrectionRound>> {
        let mut schedule = Vec::new();

        // Create syndrome extraction round
        let mut round = ErrorCorrectionRound {
            stabilizer_measurements: Vec::new(),
            syndrome_extraction: InterfaceCircuit::new(distance * distance + 100, 0), // Extra for ancillas
            error_correction: InterfaceCircuit::new(distance * distance, 0),
            duration: 1,
        };

        // Add stabilizer measurements
        let surface_code = self.surface_code.as_ref().ok_or_else(|| {
            crate::error::SimulatorError::InvalidConfiguration(
                "Surface code not initialized".to_string(),
            )
        })?;

        for (i, stabilizer) in surface_code.stabilizers.iter().enumerate() {
            let measurement = StabilizerMeasurement {
                stabilizer_index: i,
                measurement_circuit: self.create_stabilizer_measurement_circuit(stabilizer)?,
                syndrome_qubit: distance * distance + i,
                data_qubits: self.get_stabilizer_data_qubits(stabilizer),
            };
            round.stabilizer_measurements.push(measurement);
        }

        schedule.push(round);
        Ok(schedule)
    }

    /// Create stabilizer measurement circuit
    fn create_stabilizer_measurement_circuit(
        &self,
        stabilizer: &Array1<i8>,
    ) -> Result<InterfaceCircuit> {
        let mut circuit = InterfaceCircuit::new(stabilizer.len() + 1, 0); // +1 for ancilla
        let ancilla_qubit = stabilizer.len();

        // Initialize ancilla in |+⟩ state
        circuit.add_gate(InterfaceGate::new(
            InterfaceGateType::Hadamard,
            vec![ancilla_qubit],
        ));

        // Apply controlled operations
        for (i, &op) in stabilizer.iter().enumerate() {
            if op == 1 {
                if i < stabilizer.len() / 2 {
                    // X operation
                    circuit.add_gate(InterfaceGate::new(
                        InterfaceGateType::CNOT,
                        vec![ancilla_qubit, i],
                    ));
                } else {
                    // Z operation
                    let data_qubit = i - stabilizer.len() / 2;
                    circuit.add_gate(InterfaceGate::new(
                        InterfaceGateType::CZ,
                        vec![ancilla_qubit, data_qubit],
                    ));
                }
            }
        }

        // Measure ancilla
        circuit.add_gate(InterfaceGate::new(
            InterfaceGateType::Hadamard,
            vec![ancilla_qubit],
        ));

        Ok(circuit)
    }

    /// Get data qubits involved in stabilizer
    fn get_stabilizer_data_qubits(&self, stabilizer: &Array1<i8>) -> Vec<usize> {
        let mut data_qubits = Vec::new();
        let half_len = stabilizer.len() / 2;

        for i in 0..half_len {
            if stabilizer[i] == 1 || stabilizer[i + half_len] == 1 {
                data_qubits.push(i);
            }
        }

        data_qubits
    }

    /// Initialize magic state protocols
    fn initialize_magic_state_protocols(&mut self) -> Result<()> {
        // T-state protocol for T gates
        let t_protocol = MagicStateProtocol {
            input_state: MagicStateType::TState,
            output_state: MagicStateType::TState,
            distillation_circuit: self.create_t_state_distillation_circuit()?,
            error_reduction: 0.1, // 10x error reduction
            overhead: 15,         // 15 T-states input for 1 T-state output
        };
        self.magic_state_protocols
            .insert(LogicalGateType::LogicalT, t_protocol);

        // CCZ-state protocol for Toffoli gates
        let ccz_protocol = MagicStateProtocol {
            input_state: MagicStateType::CCZState,
            output_state: MagicStateType::CCZState,
            distillation_circuit: self.create_ccz_state_distillation_circuit()?,
            error_reduction: 0.05, // 20x error reduction
            overhead: 25,          // 25 CCZ-states input for 1 CCZ-state output
        };
        self.magic_state_protocols
            .insert(LogicalGateType::LogicalToffoli, ccz_protocol);

        Ok(())
    }

    /// Create T-state distillation circuit
    fn create_t_state_distillation_circuit(&self) -> Result<InterfaceCircuit> {
        let mut circuit = InterfaceCircuit::new(15, 0); // 15-to-1 distillation

        // Simplified 15-to-1 T-state distillation
        // In practice, this would be a complex multi-level protocol

        // First level: 7-to-1 distillation (using Steane code)
        for i in 0..7 {
            circuit.add_gate(InterfaceGate::new(InterfaceGateType::Hadamard, vec![i]));
        }

        // Add stabilizer measurements for error detection
        for i in 0..3 {
            circuit.add_gate(InterfaceGate::new(InterfaceGateType::CNOT, vec![i, i + 7]));
            circuit.add_gate(InterfaceGate::new(
                InterfaceGateType::CNOT,
                vec![i + 3, i + 7],
            ));
        }

        // Second level: 15-to-1 using two 7-to-1 outputs
        circuit.add_gate(InterfaceGate::new(InterfaceGateType::CNOT, vec![13, 14]));

        Ok(circuit)
    }

    /// Create CCZ-state distillation circuit
    fn create_ccz_state_distillation_circuit(&self) -> Result<InterfaceCircuit> {
        let mut circuit = InterfaceCircuit::new(25, 0); // 25-to-1 distillation

        // Simplified CCZ-state distillation
        // This would typically involve multiple rounds of error detection and correction

        for i in 0..5 {
            for j in 0..5 {
                let qubit = i * 5 + j;
                circuit.add_gate(InterfaceGate::new(InterfaceGateType::Hadamard, vec![qubit]));
            }
        }

        // Add CCZ gates for entanglement
        for i in 0..20 {
            circuit.add_gate(InterfaceGate::new(
                InterfaceGateType::Toffoli,
                vec![i, i + 1, i + 2],
            ));
        }

        Ok(circuit)
    }

    /// Initialize gate library with common logical gates
    fn initialize_gate_library(&mut self) -> Result<()> {
        // Logical Pauli-X
        let logical_x = LogicalGate {
            gate_type: LogicalGateType::LogicalX,
            logical_qubits: vec![0],
            physical_implementation: self.create_logical_pauli_x_circuit()?,
            resources: ResourceRequirements {
                physical_qubits: self.config.code_distance * self.config.code_distance,
                physical_gates: 1,
                measurement_rounds: 0,
                magic_states: 0,
                time_steps: 1,
                ancilla_qubits: 0,
            },
            error_rate: self.calculate_logical_gate_error_rate(LogicalGateType::LogicalX)?,
        };
        self.gate_library
            .insert(LogicalGateType::LogicalX, logical_x);

        // Similar for other Clifford gates...

        Ok(())
    }

    /// Create logical Pauli-X circuit
    fn create_logical_pauli_x_circuit(&self) -> Result<InterfaceCircuit> {
        let distance = self.config.code_distance;
        let mut circuit = InterfaceCircuit::new(distance * distance, 0);

        // Apply X gates to logical X string
        for i in 0..distance {
            circuit.add_gate(InterfaceGate::new(InterfaceGateType::PauliX, vec![i]));
        }

        Ok(circuit)
    }

    /// Calculate logical gate error rate
    fn calculate_logical_gate_error_rate(&self, gate_type: LogicalGateType) -> Result<f64> {
        let p_phys = self.config.physical_error_rate;
        let d = self.config.code_distance;

        // Simplified error rate calculation
        // Real calculation would depend on specific error correction protocol
        match gate_type {
            LogicalGateType::LogicalX | LogicalGateType::LogicalY | LogicalGateType::LogicalZ => {
                // Pauli gates: error suppression ~ (p_phys)^((d+1)/2)
                Ok(p_phys.powf((d + 1) as f64 / 2.0))
            }
            LogicalGateType::LogicalH | LogicalGateType::LogicalS => {
                // Clifford gates: similar suppression but with overhead
                Ok(2.0 * p_phys.powf((d + 1) as f64 / 2.0))
            }
            LogicalGateType::LogicalT => {
                // T gate requires magic states: higher error rate
                Ok(10.0 * p_phys.powf((d + 1) as f64 / 2.0))
            }
            _ => Ok(p_phys), // Conservative estimate
        }
    }

    /// Synthesize logical Pauli gates
    pub fn synthesize_logical_pauli(
        &self,
        gate_type: LogicalGateType,
        logical_qubits: &[usize],
    ) -> Result<LogicalGate> {
        let distance = self.config.code_distance;
        let mut circuit = InterfaceCircuit::new(distance * distance, 0);

        let physical_gate = match gate_type {
            LogicalGateType::LogicalX => InterfaceGateType::PauliX,
            LogicalGateType::LogicalY => InterfaceGateType::PauliY,
            LogicalGateType::LogicalZ => InterfaceGateType::PauliZ,
            _ => {
                return Err(SimulatorError::InvalidConfiguration(
                    "Invalid Pauli gate".to_string(),
                ))
            }
        };

        // Apply gate to logical string
        for i in 0..distance {
            circuit.add_gate(InterfaceGate::new(physical_gate.clone(), vec![i]));
        }

        Ok(LogicalGate {
            gate_type,
            logical_qubits: logical_qubits.to_vec(),
            physical_implementation: circuit,
            resources: ResourceRequirements {
                physical_qubits: distance * distance,
                physical_gates: distance,
                measurement_rounds: 0,
                magic_states: 0,
                time_steps: 1,
                ancilla_qubits: 0,
            },
            error_rate: self.calculate_logical_gate_error_rate(gate_type)?,
        })
    }

    /// Synthesize logical Hadamard gate
    pub fn synthesize_logical_hadamard(&self, logical_qubits: &[usize]) -> Result<LogicalGate> {
        let distance = self.config.code_distance;
        let mut circuit = InterfaceCircuit::new(distance * distance, 0);

        // Logical Hadamard: transversal for many codes
        for i in 0..distance * distance {
            circuit.add_gate(InterfaceGate::new(InterfaceGateType::Hadamard, vec![i]));
        }

        Ok(LogicalGate {
            gate_type: LogicalGateType::LogicalH,
            logical_qubits: logical_qubits.to_vec(),
            physical_implementation: circuit,
            resources: ResourceRequirements {
                physical_qubits: distance * distance,
                physical_gates: distance * distance,
                measurement_rounds: 0,
                magic_states: 0,
                time_steps: 1,
                ancilla_qubits: 0,
            },
            error_rate: self.calculate_logical_gate_error_rate(LogicalGateType::LogicalH)?,
        })
    }

    /// Synthesize logical S gate
    fn synthesize_logical_s(&self, logical_qubits: &[usize]) -> Result<LogicalGate> {
        let distance = self.config.code_distance;
        let mut circuit = InterfaceCircuit::new(distance * distance, 0);

        // Logical S: transversal for CSS codes
        for i in 0..distance * distance {
            circuit.add_gate(InterfaceGate::new(InterfaceGateType::S, vec![i]));
        }

        Ok(LogicalGate {
            gate_type: LogicalGateType::LogicalS,
            logical_qubits: logical_qubits.to_vec(),
            physical_implementation: circuit,
            resources: ResourceRequirements {
                physical_qubits: distance * distance,
                physical_gates: distance * distance,
                measurement_rounds: 0,
                magic_states: 0,
                time_steps: 1,
                ancilla_qubits: 0,
            },
            error_rate: self.calculate_logical_gate_error_rate(LogicalGateType::LogicalS)?,
        })
    }

    /// Synthesize logical T gate using magic states
    fn synthesize_logical_t_with_magic_states(
        &self,
        logical_qubits: &[usize],
    ) -> Result<LogicalGate> {
        let distance = self.config.code_distance;
        let mut circuit = InterfaceCircuit::new(distance * distance + 10, 0); // Extra qubits for magic state

        // Magic state injection protocol
        // 1. Prepare magic state |T⟩ = T|+⟩
        circuit.add_gate(InterfaceGate::new(
            InterfaceGateType::Hadamard,
            vec![distance * distance],
        ));
        circuit.add_gate(InterfaceGate::new(
            InterfaceGateType::T,
            vec![distance * distance],
        ));

        // 2. Teleport T gate using magic state
        for i in 0..distance {
            circuit.add_gate(InterfaceGate::new(
                InterfaceGateType::CNOT,
                vec![i, distance * distance + 1],
            ));
        }

        // 3. Measure magic state and apply corrections
        circuit.add_gate(InterfaceGate::new(
            InterfaceGateType::Hadamard,
            vec![distance * distance],
        ));

        Ok(LogicalGate {
            gate_type: LogicalGateType::LogicalT,
            logical_qubits: logical_qubits.to_vec(),
            physical_implementation: circuit,
            resources: ResourceRequirements {
                physical_qubits: distance * distance + 10,
                physical_gates: distance + 3,
                measurement_rounds: 1,
                magic_states: 1,
                time_steps: 5,
                ancilla_qubits: 10,
            },
            error_rate: self.calculate_logical_gate_error_rate(LogicalGateType::LogicalT)?,
        })
    }

    /// Synthesize logical CNOT gate
    pub fn synthesize_logical_cnot(&self, logical_qubits: &[usize]) -> Result<LogicalGate> {
        if logical_qubits.len() != 2 {
            return Err(SimulatorError::InvalidConfiguration(
                "CNOT requires exactly 2 qubits".to_string(),
            ));
        }

        let distance = self.config.code_distance;
        let mut circuit = InterfaceCircuit::new(2 * distance * distance, 0);

        // Transversal CNOT for CSS codes
        for i in 0..distance * distance {
            circuit.add_gate(InterfaceGate::new(
                InterfaceGateType::CNOT,
                vec![i, i + distance * distance],
            ));
        }

        Ok(LogicalGate {
            gate_type: LogicalGateType::LogicalCNOT,
            logical_qubits: logical_qubits.to_vec(),
            physical_implementation: circuit,
            resources: ResourceRequirements {
                physical_qubits: 2 * distance * distance,
                physical_gates: distance * distance,
                measurement_rounds: 0,
                magic_states: 0,
                time_steps: 1,
                ancilla_qubits: 0,
            },
            error_rate: self.calculate_logical_gate_error_rate(LogicalGateType::LogicalCNOT)?,
        })
    }

    /// Synthesize logical Toffoli gate using magic states
    fn synthesize_logical_toffoli_with_magic_states(
        &self,
        logical_qubits: &[usize],
    ) -> Result<LogicalGate> {
        if logical_qubits.len() != 3 {
            return Err(SimulatorError::InvalidConfiguration(
                "Toffoli requires exactly 3 qubits".to_string(),
            ));
        }

        let distance = self.config.code_distance;
        let mut circuit = InterfaceCircuit::new(3 * distance * distance + 20, 0);

        // CCZ magic state injection protocol
        // Complex protocol involving multiple magic states and measurements

        // Prepare CCZ magic state
        for i in 0..3 {
            circuit.add_gate(InterfaceGate::new(
                InterfaceGateType::Hadamard,
                vec![3 * distance * distance + i],
            ));
        }

        // Apply CCZ to magic state
        circuit.add_gate(InterfaceGate::new(
            InterfaceGateType::Toffoli,
            vec![
                3 * distance * distance,
                3 * distance * distance + 1,
                3 * distance * distance + 2,
            ],
        ));

        // Teleportation protocol (simplified)
        for i in 0..distance * distance {
            circuit.add_gate(InterfaceGate::new(
                InterfaceGateType::CNOT,
                vec![i, 3 * distance * distance + 3],
            ));
            circuit.add_gate(InterfaceGate::new(
                InterfaceGateType::CNOT,
                vec![i + distance * distance, 3 * distance * distance + 4],
            ));
            circuit.add_gate(InterfaceGate::new(
                InterfaceGateType::CNOT,
                vec![i + 2 * distance * distance, 3 * distance * distance + 5],
            ));
        }

        Ok(LogicalGate {
            gate_type: LogicalGateType::LogicalToffoli,
            logical_qubits: logical_qubits.to_vec(),
            physical_implementation: circuit,
            resources: ResourceRequirements {
                physical_qubits: 3 * distance * distance + 20,
                physical_gates: 4 + 3 * distance * distance,
                measurement_rounds: 3,
                magic_states: 1, // CCZ magic state
                time_steps: 10,
                ancilla_qubits: 20,
            },
            error_rate: self.calculate_logical_gate_error_rate(LogicalGateType::LogicalToffoli)?,
        })
    }

    /// Helper methods and remaining implementation details...
    fn map_interface_gate_to_logical(&self, gate: &InterfaceGate) -> Result<LogicalGateType> {
        match gate.gate_type {
            InterfaceGateType::PauliX => Ok(LogicalGateType::LogicalX),
            InterfaceGateType::PauliY => Ok(LogicalGateType::LogicalY),
            InterfaceGateType::PauliZ => Ok(LogicalGateType::LogicalZ),
            InterfaceGateType::Hadamard => Ok(LogicalGateType::LogicalH),
            InterfaceGateType::S => Ok(LogicalGateType::LogicalS),
            InterfaceGateType::T => Ok(LogicalGateType::LogicalT),
            InterfaceGateType::CNOT => Ok(LogicalGateType::LogicalCNOT),
            InterfaceGateType::Toffoli => Ok(LogicalGateType::LogicalToffoli),
            InterfaceGateType::RY(angle) => Ok(LogicalGateType::LogicalRotation(angle)),
            InterfaceGateType::RX(angle) => Ok(LogicalGateType::LogicalRotation(angle)),
            InterfaceGateType::RZ(angle) => Ok(LogicalGateType::LogicalRotation(angle)),
            _ => Err(SimulatorError::InvalidConfiguration(format!(
                "Unsupported gate type for logical synthesis: {:?}",
                gate.gate_type
            ))),
        }
    }

    fn append_synthesized_gate(
        &self,
        circuit: &mut InterfaceCircuit,
        gate: &LogicalGate,
    ) -> Result<()> {
        // Append the physical implementation to the fault-tolerant circuit
        for physical_gate in &gate.physical_implementation.gates {
            circuit.add_gate(physical_gate.clone());
        }
        Ok(())
    }

    fn update_resources(&self, total: &mut ResourceRequirements, gate: &ResourceRequirements) {
        total.physical_qubits = total.physical_qubits.max(gate.physical_qubits);
        total.physical_gates += gate.physical_gates;
        total.measurement_rounds += gate.measurement_rounds;
        total.magic_states += gate.magic_states;
        total.time_steps += gate.time_steps;
        total.ancilla_qubits = total.ancilla_qubits.max(gate.ancilla_qubits);
    }

    fn add_error_correction_rounds(
        &self,
        circuit: &mut InterfaceCircuit,
        distance: usize,
    ) -> Result<()> {
        // Add periodic error correction rounds
        let rounds_needed = circuit.gates.len() / 10; // Every 10 gates

        for _ in 0..rounds_needed {
            // Add syndrome extraction
            for i in 0..distance * distance {
                circuit.add_gate(InterfaceGate::new(
                    InterfaceGateType::Hadamard,
                    vec![circuit.num_qubits + i % 10],
                ));
            }
        }

        Ok(())
    }

    fn calculate_logical_error_rate(&self, result: &FaultTolerantSynthesisResult) -> Result<f64> {
        let p_phys = self.config.physical_error_rate;
        let d = result.code_distance;
        let gate_count = result.synthesis_stats.logical_gates_synthesized;

        // Simplified calculation
        let base_error_rate = p_phys.powf((d + 1) as f64 / 2.0);
        let total_error_rate = gate_count as f64 * base_error_rate;

        Ok(total_error_rate.min(1.0))
    }

    fn calculate_optimal_distance(&self, circuit: &InterfaceCircuit) -> Result<usize> {
        let gate_count = circuit.gates.len();
        let target_error = self.config.target_logical_error_rate;
        let p_phys = self.config.physical_error_rate;

        // Find minimum distance that achieves target error rate
        for d in (3..20).step_by(2) {
            let logical_error = gate_count as f64 * p_phys.powf((d + 1) as f64 / 2.0);
            if logical_error < target_error {
                return Ok(d);
            }
        }

        Ok(19) // Maximum reasonable distance
    }

    fn generate_cache_key(&self, circuit: &InterfaceCircuit) -> String {
        // Simple cache key based on circuit structure
        format!(
            "{}_{}_{}_{}",
            circuit.num_qubits,
            circuit.gates.len(),
            self.config.code_distance,
            format!("{:?}", self.config.error_correction_code)
        )
    }

    fn initialize_resource_estimator(&mut self) -> Result<()> {
        // Initialize error model
        let mut gate_errors = HashMap::new();
        gate_errors.insert("CNOT".to_string(), 1e-3);
        gate_errors.insert("H".to_string(), 5e-4);
        gate_errors.insert("T".to_string(), 1e-3);

        self.resource_estimator.error_model = PhysicalErrorModel {
            gate_errors,
            measurement_error: 1e-3,
            memory_error: 1e-5,
            correlated_error: 1e-4,
        };

        // Initialize code parameters
        let mut code_params = HashMap::new();
        code_params.insert(
            ErrorCorrectionCode::SurfaceCode,
            CodeParameters {
                encoding_rate: 1.0 / (self.config.code_distance.pow(2) as f64),
                threshold: 1e-2,
                resource_scaling: 2.0,
                error_suppression: (self.config.code_distance + 1) as f64 / 2.0,
            },
        );

        self.resource_estimator.code_parameters = code_params;

        // Initialize magic state costs
        let mut magic_costs = HashMap::new();
        magic_costs.insert(MagicStateType::TState, 15);
        magic_costs.insert(MagicStateType::CCZState, 25);

        self.resource_estimator.magic_state_costs = magic_costs;

        Ok(())
    }

    /// Synthesize logical T with magic states (public version)
    pub fn synthesize_logical_t_with_magic_states_public(
        &self,
        logical_qubits: &[usize],
    ) -> Result<LogicalGate> {
        self.synthesize_logical_t_with_magic_states(logical_qubits)
    }

    /// Create T state distillation circuit (public version)
    pub fn create_t_state_distillation_circuit_public(&self) -> Result<InterfaceCircuit> {
        self.create_t_state_distillation_circuit()
    }

    /// Create CCZ state distillation circuit (public version)
    pub fn create_ccz_state_distillation_circuit_public(&self) -> Result<InterfaceCircuit> {
        self.create_ccz_state_distillation_circuit()
    }

    /// Synthesize logical rotation gate
    fn synthesize_logical_rotation(
        &self,
        logical_qubits: &[usize],
        angle: f64,
    ) -> Result<LogicalGate> {
        let distance = self.config.code_distance;
        let mut circuit = InterfaceCircuit::new(distance * distance + 10, 0);

        // Decompose rotation into Clifford + T gates (Solovay-Kitaev decomposition)
        // For simplicity, we'll use a basic decomposition into a few T gates
        // In practice, this would be a more sophisticated decomposition

        // Apply logical Z rotations using T gates
        // R_z(θ) ≈ sequence of T gates and Clifford operations
        let num_t_gates = ((angle.abs() / (std::f64::consts::PI / 4.0)).ceil() as usize).max(1);

        for i in 0..distance * distance {
            // Apply the rotation as a sequence of elementary operations
            if angle.abs() > 1e-10 {
                // Apply Hadamard to convert between X and Z rotations if needed
                if angle.abs() > std::f64::consts::PI / 8.0 {
                    circuit.add_gate(InterfaceGate::new(InterfaceGateType::Hadamard, vec![i]));
                }

                // Apply T gates to approximate the rotation
                for _ in 0..num_t_gates {
                    circuit.add_gate(InterfaceGate::new(InterfaceGateType::T, vec![i]));
                }

                // Apply inverse Hadamard if needed
                if angle.abs() > std::f64::consts::PI / 8.0 {
                    circuit.add_gate(InterfaceGate::new(InterfaceGateType::Hadamard, vec![i]));
                }
            }
        }

        Ok(LogicalGate {
            gate_type: LogicalGateType::LogicalRotation(angle),
            logical_qubits: logical_qubits.to_vec(),
            physical_implementation: circuit,
            resources: ResourceRequirements {
                physical_qubits: distance * distance,
                physical_gates: distance * distance * num_t_gates * 2,
                measurement_rounds: distance,
                magic_states: num_t_gates * distance * distance,
                time_steps: num_t_gates * 2,
                ancilla_qubits: 10,
            },
            error_rate: 0.001 * (num_t_gates as f64).mul_add(0.001, 1.0),
        })
    }

    /// Update resources (public version)
    pub fn update_resources_public(
        &self,
        total: &mut ResourceRequirements,
        gate: &ResourceRequirements,
    ) {
        self.update_resources(total, gate);
    }

    /// Calculate optimal distance (public version)
    pub fn calculate_optimal_distance_public(&self, circuit: &InterfaceCircuit) -> Result<usize> {
        self.calculate_optimal_distance(circuit)
    }

    /// Calculate logical gate error rate (public version)
    pub fn calculate_logical_gate_error_rate_public(
        &self,
        gate_type: LogicalGateType,
    ) -> Result<f64> {
        self.calculate_logical_gate_error_rate(gate_type)
    }
}

/// Benchmark function for fault-tolerant synthesis
pub fn benchmark_fault_tolerant_synthesis() -> Result<()> {
    println!("Benchmarking Fault-Tolerant Gate Synthesis...");

    let config = FaultTolerantConfig::default();
    let mut synthesizer = FaultTolerantSynthesizer::new(config)?;

    // Create test logical circuit
    let mut logical_circuit = InterfaceCircuit::new(2, 0);
    logical_circuit.add_gate(InterfaceGate::new(InterfaceGateType::Hadamard, vec![0]));
    logical_circuit.add_gate(InterfaceGate::new(InterfaceGateType::CNOT, vec![0, 1]));
    logical_circuit.add_gate(InterfaceGate::new(InterfaceGateType::T, vec![1]));

    let start_time = std::time::Instant::now();

    // Synthesize fault-tolerant implementation
    let result = synthesizer.synthesize_logical_circuit(&logical_circuit)?;

    let duration = start_time.elapsed();

    println!("✅ Fault-Tolerant Synthesis Results:");
    println!(
        "   Logical Gates Synthesized: {}",
        result.synthesis_stats.logical_gates_synthesized
    );
    println!(
        "   Physical Qubits Required: {}",
        result.resources.physical_qubits
    );
    println!(
        "   Physical Gates Required: {}",
        result.resources.physical_gates
    );
    println!(
        "   Magic States Consumed: {}",
        result.resources.magic_states
    );
    println!("   Code Distance: {}", result.code_distance);
    println!("   Logical Error Rate: {:.2e}", result.logical_error_rate);
    println!("   Overhead Factor: {:.1}x", result.overhead_factor);
    println!("   Synthesis Time: {:.2}ms", duration.as_millis());

    Ok(())
}

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

    #[test]
    fn test_fault_tolerant_synthesizer_creation() {
        let config = FaultTolerantConfig::default();
        let synthesizer = FaultTolerantSynthesizer::new(config);
        assert!(synthesizer.is_ok());
    }

    #[test]
    fn test_surface_code_layout_creation() {
        let config = FaultTolerantConfig::default();
        let synthesizer =
            FaultTolerantSynthesizer::new(config).expect("Failed to create synthesizer");
        let layout = synthesizer.create_surface_code_layout(3);
        assert!(layout.is_ok());
    }

    #[test]
    fn test_logical_pauli_synthesis() {
        let config = FaultTolerantConfig::default();
        let synthesizer =
            FaultTolerantSynthesizer::new(config).expect("Failed to create synthesizer");
        let result = synthesizer.synthesize_logical_pauli(LogicalGateType::LogicalX, &[0]);
        assert!(result.is_ok());
    }

    #[test]
    fn test_logical_hadamard_synthesis() {
        let config = FaultTolerantConfig::default();
        let synthesizer =
            FaultTolerantSynthesizer::new(config).expect("Failed to create synthesizer");
        let result = synthesizer.synthesize_logical_hadamard(&[0]);
        assert!(result.is_ok());
    }

    #[test]
    fn test_logical_cnot_synthesis() {
        let config = FaultTolerantConfig::default();
        let synthesizer =
            FaultTolerantSynthesizer::new(config).expect("Failed to create synthesizer");
        let result = synthesizer.synthesize_logical_cnot(&[0, 1]);
        assert!(result.is_ok());
    }

    #[test]
    fn test_resource_requirements_update() {
        let config = FaultTolerantConfig::default();
        let synthesizer =
            FaultTolerantSynthesizer::new(config).expect("Failed to create synthesizer");

        let mut total = ResourceRequirements::default();
        let gate_resources = ResourceRequirements {
            physical_qubits: 10,
            physical_gates: 5,
            measurement_rounds: 1,
            magic_states: 2,
            time_steps: 3,
            ancilla_qubits: 4,
        };

        synthesizer.update_resources(&mut total, &gate_resources);

        assert_eq!(total.physical_qubits, 10);
        assert_eq!(total.physical_gates, 5);
        assert_eq!(total.magic_states, 2);
    }

    #[test]
    fn test_optimal_distance_calculation() {
        let config = FaultTolerantConfig {
            target_logical_error_rate: 1e-10,
            physical_error_rate: 1e-3,
            ..FaultTolerantConfig::default()
        };
        let synthesizer =
            FaultTolerantSynthesizer::new(config).expect("Failed to create synthesizer");

        let circuit = InterfaceCircuit::new(2, 0);
        let distance = synthesizer.calculate_optimal_distance(&circuit);
        assert!(distance.is_ok());
        let distance_value = distance.expect("Failed to calculate optimal distance");
        assert!(distance_value >= 3);
    }

    /// Regression: the stabilizer neighbour lists used a modular-arithmetic
    /// placeholder `(i*d + j + offset) % (d*d)` with no relation to lattice
    /// geometry. They now return genuine star / plaquette adjacency on the
    /// data-qubit grid, with boundary truncation.
    #[test]
    fn test_surface_code_stabilizer_neighbors_are_geometric() {
        let synthesizer = FaultTolerantSynthesizer::new(FaultTolerantConfig::default())
            .expect("Failed to create synthesizer");
        let distance = 3usize;

        // Bulk X-star at vertex (1,1): the plus of orthogonal neighbours.
        let mut x_bulk = synthesizer.get_x_stabilizer_neighbors(1, 1, distance);
        x_bulk.sort_unstable();
        assert_eq!(x_bulk, vec![1, 3, 5, 7]); // (0,1),(1,0),(1,2),(2,1)

        // Corner X-star (0,0) truncates to two in-grid neighbours (no wrap).
        let mut x_corner = synthesizer.get_x_stabilizer_neighbors(0, 0, distance);
        x_corner.sort_unstable();
        assert_eq!(x_corner, vec![1, 3]); // (0,1),(1,0)

        // Bulk Z-plaquette at face (0,0): its four cell corners.
        let mut z_bulk = synthesizer.get_z_stabilizer_neighbors(0, 0, distance);
        z_bulk.sort_unstable();
        assert_eq!(z_bulk, vec![0, 1, 3, 4]); // (0,0),(0,1),(1,0),(1,1)

        // Boundary Z-plaquette at the bottom edge truncates to weight 2.
        let mut z_edge = synthesizer.get_z_stabilizer_neighbors(2, 0, distance);
        z_edge.sort_unstable();
        assert_eq!(z_edge, vec![6, 7]); // (2,0),(2,1); (3,*) off-grid

        // The old modular placeholder produced neighbours like {4,5,6,7} for
        // X-stab (1,1); assert we no longer see that meaningless wraparound.
        assert_ne!(x_bulk, vec![4, 5, 6, 7]);

        // Every neighbour is a valid data-qubit index and geometrically
        // adjacent to the stabilizer cell (Chebyshev distance <= 1), and the
        // lists contain no duplicates.
        let total = distance * distance;
        for i in 0..distance - 1 {
            for j in 0..distance {
                let n = synthesizer.get_x_stabilizer_neighbors(i, j, distance);
                let mut seen = std::collections::HashSet::new();
                for &q in &n {
                    assert!(q < total, "neighbour {q} out of range");
                    let (r, c) = (q / distance, q % distance);
                    let dr = (r as isize - i as isize).unsigned_abs();
                    let dc = (c as isize - j as isize).unsigned_abs();
                    assert!(dr <= 1 && dc <= 1, "neighbour not adjacent");
                    assert!(seen.insert(q), "duplicate neighbour {q}");
                }
            }
        }
    }
}