quantrs2_device/

optimization.rs

//! Circuit optimization using device calibration data
//!
//! This module provides optimization strategies that leverage device-specific
//! calibration data to improve circuit performance on real hardware.

use std::cmp::Ordering;
use std::collections::{HashMap, HashSet};

use quantrs2_circuit::builder::Circuit;
use quantrs2_core::{
    error::{QuantRS2Error, QuantRS2Result},
    gate::GateOp,
    qubit::QubitId,
};

use crate::calibration::{CalibrationManager, DeviceCalibration, PulseShape};

/// Circuit optimizer that uses device calibration data
pub struct CalibrationOptimizer {
    /// Calibration manager
    calibration_manager: CalibrationManager,
    /// Optimization configuration
    config: OptimizationConfig,
}

/// Configuration for calibration-based optimization
#[derive(Debug, Clone)]
pub struct OptimizationConfig {
    /// Optimize for gate fidelity
    pub optimize_fidelity: bool,
    /// Optimize for circuit duration
    pub optimize_duration: bool,
    /// Allow gate substitutions
    pub allow_substitutions: bool,
    /// Maximum acceptable fidelity loss for substitutions
    pub fidelity_threshold: f64,
    /// Consider crosstalk in optimization
    pub consider_crosstalk: bool,
    /// Prefer native gates
    pub prefer_native_gates: bool,
    /// Maximum circuit depth increase allowed
    pub max_depth_increase: f64,
}

impl Default for OptimizationConfig {
    fn default() -> Self {
        Self {
            optimize_fidelity: true,
            optimize_duration: true,
            allow_substitutions: true,
            fidelity_threshold: 0.99,
            consider_crosstalk: true,
            prefer_native_gates: true,
            max_depth_increase: 1.5,
        }
    }
}

/// Result of circuit optimization
#[derive(Debug, Clone)]
pub struct OptimizationResult<const N: usize> {
    /// Optimized circuit
    pub circuit: Circuit<N>,
    /// Estimated fidelity
    pub estimated_fidelity: f64,
    /// Estimated duration (ns)
    pub estimated_duration: f64,
    /// Number of gates before optimization
    pub original_gate_count: usize,
    /// Number of gates after optimization
    pub optimized_gate_count: usize,
    /// Optimization decisions made
    pub decisions: Vec<OptimizationDecision>,
}

/// Individual optimization decision
#[derive(Debug, Clone)]
pub enum OptimizationDecision {
    /// Gate was substituted
    GateSubstitution {
        original: String,
        replacement: String,
        qubits: Vec<QubitId>,
        fidelity_change: f64,
        duration_change: f64,
    },
    /// Gates were reordered
    GateReordering { gates: Vec<String>, reason: String },
    /// Gate was moved to different qubits
    QubitRemapping {
        gate: String,
        original_qubits: Vec<QubitId>,
        new_qubits: Vec<QubitId>,
        reason: String,
    },
    /// Gate decomposition was changed
    DecompositionChange {
        gate: String,
        qubits: Vec<QubitId>,
        original_depth: usize,
        new_depth: usize,
    },
}

impl CalibrationOptimizer {
    /// Create a new calibration-based optimizer
    pub const fn new(calibration_manager: CalibrationManager, config: OptimizationConfig) -> Self {
        Self {
            calibration_manager,
            config,
        }
    }

    /// Optimize a circuit for a specific device
    pub fn optimize_circuit<const N: usize>(
        &self,
        circuit: &Circuit<N>,
        device_id: &str,
    ) -> QuantRS2Result<OptimizationResult<N>> {
        // Check if calibration is available and valid
        if !self.calibration_manager.is_calibration_valid(device_id) {
            return Err(QuantRS2Error::InvalidInput(format!(
                "No valid calibration for device {device_id}"
            )));
        }

        let calibration = self
            .calibration_manager
            .get_calibration(device_id)
            .ok_or_else(|| QuantRS2Error::InvalidInput("Calibration not found".into()))?;

        // Clone the circuit for optimization
        let mut optimized_circuit = circuit.clone();
        let mut decisions = Vec::new();

        // Apply optimization strategies based on configuration
        if self.config.optimize_fidelity {
            self.optimize_for_fidelity(&mut optimized_circuit, calibration, &mut decisions)?;
        }

        if self.config.optimize_duration {
            Self::optimize_for_duration(&mut optimized_circuit, calibration, &mut decisions)?;
        }

        if self.config.allow_substitutions {
            Self::apply_gate_substitutions(&mut optimized_circuit, calibration, &mut decisions)?;
        }

        if self.config.consider_crosstalk {
            Self::mitigate_crosstalk(&mut optimized_circuit, calibration, &mut decisions)?;
        }

        // Estimate final metrics
        let estimated_fidelity = Self::estimate_circuit_fidelity(&optimized_circuit, calibration)?;
        let estimated_duration = Self::estimate_circuit_duration(&optimized_circuit, calibration)?;

        Ok(OptimizationResult {
            circuit: optimized_circuit,
            estimated_fidelity,
            estimated_duration,
            original_gate_count: circuit.gates().len(),
            optimized_gate_count: circuit.gates().len(), // This would be updated by actual optimization
            decisions,
        })
    }

    /// Optimize circuit for maximum fidelity
    fn optimize_for_fidelity<const N: usize>(
        &self,
        circuit: &mut Circuit<N>,
        calibration: &DeviceCalibration,
        decisions: &mut Vec<OptimizationDecision>,
    ) -> QuantRS2Result<()> {
        // Strategy 1: Use highest fidelity qubits for critical gates
        let qubit_qualities = Self::rank_qubits_by_quality(calibration);

        // Strategy 2: Minimize two-qubit gate count by decomposing into single-qubit gates where possible
        let mut optimized_gates: Vec<
            std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>,
        > = Vec::new();
        let original_gates = circuit.gates();

        for gate in original_gates {
            let qubits = gate.qubits();

            if qubits.len() == 2 {
                // Check if this two-qubit gate can be decomposed or substituted
                let (q1, q2) = (qubits[0], qubits[1]);

                // Get gate fidelities
                let single_q1_fidelity = calibration
                    .single_qubit_gates
                    .get(gate.name())
                    .and_then(|gate_cal| gate_cal.qubit_data.get(&q1))
                    .map_or(0.999, |data| data.fidelity);

                let single_q2_fidelity = calibration
                    .single_qubit_gates
                    .get(gate.name())
                    .and_then(|gate_cal| gate_cal.qubit_data.get(&q2))
                    .map_or(0.999, |data| data.fidelity);

                let two_qubit_fidelity = calibration
                    .two_qubit_gates
                    .get(&(q1, q2))
                    .map_or(0.99, |gate_cal| gate_cal.fidelity);

                // If decomposition into single-qubit gates would be more faithful
                let decomposition_fidelity = single_q1_fidelity * single_q2_fidelity;

                if decomposition_fidelity > two_qubit_fidelity && gate.name() != "CNOT" {
                    // Attempt decomposition for certain gates
                    if let Some(decomposed_gates) = Self::try_decompose_gate(gate.as_ref()) {
                        let new_depth = decomposed_gates.len();
                        // Convert Box<dyn GateOp> to Arc<dyn GateOp + Send + Sync>
                        for decomposed_gate in decomposed_gates {
                            // Skip conversion for now since try_decompose_gate returns None anyway
                        }

                        decisions.push(OptimizationDecision::DecompositionChange {
                            gate: gate.name().to_string(),
                            qubits: qubits.clone(),
                            original_depth: 1,
                            new_depth,
                        });
                        continue;
                    }
                }

                // Strategy 3: Remap to higher fidelity qubit pairs if available
                if let Some((better_q1, better_q2)) =
                    Self::find_better_qubit_pair(&(q1, q2), calibration)
                {
                    decisions.push(OptimizationDecision::QubitRemapping {
                        gate: gate.name().to_string(),
                        original_qubits: vec![q1, q2],
                        new_qubits: vec![better_q1, better_q2],
                        reason: format!(
                            "Higher fidelity pair: {:.4} vs {:.4}",
                            calibration
                                .two_qubit_gates
                                .get(&(better_q1, better_q2))
                                .map_or(0.99, |g| g.fidelity),
                            two_qubit_fidelity
                        ),
                    });
                }
            }

            optimized_gates.push(gate.clone());
        }

        // Strategy 4: Reorder gates to use highest quality qubits for most critical operations
        if self.config.prefer_native_gates {
            Self::prioritize_native_gates(&mut optimized_gates, calibration, decisions)?;
        }

        Ok(())
    }

    /// Optimize circuit for minimum duration
    fn optimize_for_duration<const N: usize>(
        circuit: &mut Circuit<N>,
        calibration: &DeviceCalibration,
        decisions: &mut Vec<OptimizationDecision>,
    ) -> QuantRS2Result<()> {
        // Strategy 1: Parallelize gates where possible
        let parallel_groups = Self::identify_parallelizable_gates(circuit, calibration)?;

        if parallel_groups.len() > 1 {
            decisions.push(OptimizationDecision::GateReordering {
                gates: parallel_groups
                    .iter()
                    .flat_map(|group| group.iter().map(|g| g.name().to_string()))
                    .collect(),
                reason: format!("Parallelized {} gate groups", parallel_groups.len()),
            });
        }

        // Strategy 2: Use faster gate implementations
        let original_gates = circuit.gates();
        for (i, gate) in original_gates.iter().enumerate() {
            let qubits = gate.qubits();

            // For single-qubit gates, check if there's a faster implementation
            if qubits.len() == 1 {
                let qubit = qubits[0];
                if let Some(gate_cal) = calibration.single_qubit_gates.get(gate.name()) {
                    if let Some(qubit_data) = gate_cal.qubit_data.get(&qubit) {
                        // Check for alternative implementations
                        let faster_alternatives =
                            Self::find_faster_gate_alternatives(gate.name(), qubit_data);

                        if let Some((alt_name, duration_improvement)) = faster_alternatives {
                            decisions.push(OptimizationDecision::GateSubstitution {
                                original: gate.name().to_string(),
                                replacement: alt_name,
                                qubits: qubits.clone(),
                                fidelity_change: -0.001, // Slight fidelity trade-off for speed
                                duration_change: -duration_improvement,
                            });
                        }
                    }
                }
            }

            // For two-qubit gates, optimize timing
            if qubits.len() == 2 {
                let (q1, q2) = (qubits[0], qubits[1]);
                if let Some(gate_cal) = calibration.two_qubit_gates.get(&(q1, q2)) {
                    // Check if there's a faster coupling direction
                    if let Some(reverse_cal) = calibration.two_qubit_gates.get(&(q2, q1)) {
                        if reverse_cal.duration < gate_cal.duration {
                            decisions.push(OptimizationDecision::QubitRemapping {
                                gate: gate.name().to_string(),
                                original_qubits: vec![q1, q2],
                                new_qubits: vec![q2, q1],
                                reason: format!(
                                    "Faster coupling direction: {:.1}ns vs {:.1}ns",
                                    reverse_cal.duration, gate_cal.duration
                                ),
                            });
                        }
                    }
                }
            }
        }

        // Strategy 3: Minimize circuit depth by removing redundant operations
        Self::remove_redundant_gates(circuit, decisions)?;

        // Strategy 4: Optimize gate scheduling based on hardware timing constraints
        Self::optimize_gate_scheduling(circuit, calibration, decisions)?;

        Ok(())
    }

    /// Apply gate substitutions based on calibration
    fn apply_gate_substitutions<const N: usize>(
        circuit: &mut Circuit<N>,
        calibration: &DeviceCalibration,
        decisions: &mut Vec<OptimizationDecision>,
    ) -> QuantRS2Result<()> {
        let original_gates = circuit.gates();

        for gate in original_gates {
            let qubits = gate.qubits();
            let gate_name = gate.name();

            // Strategy 1: Virtual Z gates (RZ gates can often be implemented virtually)
            if gate_name.starts_with("RZ") || gate_name.starts_with("Rz") {
                // Z rotations can often be implemented as virtual gates with zero duration
                decisions.push(OptimizationDecision::GateSubstitution {
                    original: gate_name.to_string(),
                    replacement: "Virtual_RZ".to_string(),
                    qubits: qubits.clone(),
                    fidelity_change: 0.001, // Virtual gates are typically more faithful
                    duration_change: -30.0, // Save gate duration
                });
                continue;
            }

            // Strategy 2: Native gate substitutions
            if qubits.len() == 1 {
                let qubit = qubits[0];

                // Check if this gate can be replaced with a native gate
                if let Some(native_replacement) =
                    Self::find_native_replacement(gate_name, calibration)
                {
                    if let Some(gate_cal) = calibration.single_qubit_gates.get(&native_replacement)
                    {
                        if let Some(qubit_data) = gate_cal.qubit_data.get(&qubit) {
                            // Compare fidelities
                            let original_fidelity = calibration
                                .single_qubit_gates
                                .get(gate_name)
                                .and_then(|g| g.qubit_data.get(&qubit))
                                .map_or(0.999, |d| d.fidelity);

                            if qubit_data.fidelity > original_fidelity {
                                decisions.push(OptimizationDecision::GateSubstitution {
                                    original: gate_name.to_string(),
                                    replacement: native_replacement.clone(),
                                    qubits: qubits.clone(),
                                    fidelity_change: qubit_data.fidelity - original_fidelity,
                                    duration_change: qubit_data.duration
                                        - calibration
                                            .single_qubit_gates
                                            .get(gate_name)
                                            .and_then(|g| g.qubit_data.get(&qubit))
                                            .map_or(30.0, |d| d.duration),
                                });
                            }
                        }
                    }
                }

                // Strategy 3: Composite gate decomposition
                if let Some(decomposition) = Self::find_composite_decomposition(gate_name) {
                    // Check if decomposition improves overall fidelity
                    let mut total_decomp_fidelity = 1.0;
                    let mut total_decomp_duration = 0.0;

                    for decomp_gate in &decomposition {
                        if let Some(gate_cal) = calibration.single_qubit_gates.get(decomp_gate) {
                            if let Some(qubit_data) = gate_cal.qubit_data.get(&qubit) {
                                total_decomp_fidelity *= qubit_data.fidelity;
                                total_decomp_duration += qubit_data.duration;
                            }
                        }
                    }

                    let original_fidelity = calibration
                        .single_qubit_gates
                        .get(gate_name)
                        .and_then(|g| g.qubit_data.get(&qubit))
                        .map_or(0.999, |d| d.fidelity);

                    // Only substitute if it improves fidelity or duration significantly
                    if total_decomp_fidelity > original_fidelity + 0.001
                        || (total_decomp_fidelity > original_fidelity - 0.001
                            && decomposition.len() < 3)
                    {
                        decisions.push(OptimizationDecision::DecompositionChange {
                            gate: gate_name.to_string(),
                            qubits: qubits.clone(),
                            original_depth: 1,
                            new_depth: decomposition.len(),
                        });
                    }
                }
            }

            // Strategy 4: Two-qubit gate optimizations
            if qubits.len() == 2 {
                let (q1, q2) = (qubits[0], qubits[1]);

                // Check for equivalent gates with better calibration
                if gate_name == "CNOT" || gate_name == "CX" {
                    // Check if CZ + single qubit rotations would be better
                    if let Some(cz_cal) = calibration.two_qubit_gates.get(&(q1, q2)) {
                        let cnot_fidelity = calibration
                            .two_qubit_gates
                            .get(&(q1, q2))
                            .map_or(0.99, |g| g.fidelity);

                        // CZ + H decomposition might be more faithful
                        let h_fidelity = calibration
                            .single_qubit_gates
                            .get("H")
                            .and_then(|g| g.qubit_data.get(&q2))
                            .map_or(0.999, |d| d.fidelity);

                        let decomp_fidelity = cz_cal.fidelity * h_fidelity * h_fidelity;

                        if decomp_fidelity > cnot_fidelity + 0.005 {
                            decisions.push(OptimizationDecision::GateSubstitution {
                                original: "CNOT".to_string(),
                                replacement: "CZ_H_decomposition".to_string(),
                                qubits: qubits.clone(),
                                fidelity_change: decomp_fidelity - cnot_fidelity,
                                duration_change: cz_cal.duration + 60.0
                                    - calibration
                                        .two_qubit_gates
                                        .get(&(q1, q2))
                                        .map_or(300.0, |g| g.duration),
                            });
                        }
                    }
                }
            }
        }

        Ok(())
    }

    /// Mitigate crosstalk effects
    fn mitigate_crosstalk<const N: usize>(
        circuit: &mut Circuit<N>,
        calibration: &DeviceCalibration,
        decisions: &mut Vec<OptimizationDecision>,
    ) -> QuantRS2Result<()> {
        let original_gates = circuit.gates();

        // Strategy 1: Analyze crosstalk patterns and identify problematic gate pairs
        let crosstalk_matrix = &calibration.crosstalk_matrix;
        let mut problematic_pairs = Vec::new();

        // Find gates that would execute simultaneously and have high crosstalk
        for (i, gate1) in original_gates.iter().enumerate() {
            for (j, gate2) in original_gates.iter().enumerate() {
                if i >= j {
                    continue;
                }

                // Check if these gates could execute in parallel
                let gate1_qubits = gate1.qubits();
                let gate2_qubits = gate2.qubits();

                // If gates don't share qubits, they could potentially be parallel
                let mut overlap = false;
                for &q1 in &gate1_qubits {
                    for &q2 in &gate2_qubits {
                        if q1 == q2 {
                            overlap = true;
                            break;
                        }
                    }
                    if overlap {
                        break;
                    }
                }

                if !overlap {
                    // Check crosstalk between all qubit pairs
                    let mut max_crosstalk: f32 = 0.0;
                    for &q1 in &gate1_qubits {
                        for &q2 in &gate2_qubits {
                            let q1_idx = q1.0 as usize;
                            let q2_idx = q2.0 as usize;

                            if q1_idx < crosstalk_matrix.matrix.len()
                                && q2_idx < crosstalk_matrix.matrix[q1_idx].len()
                            {
                                let crosstalk = crosstalk_matrix.matrix[q1_idx][q2_idx] as f32;
                                max_crosstalk = max_crosstalk.max(crosstalk);
                            }
                        }
                    }

                    // If crosstalk is significant, record this pair
                    if max_crosstalk > 0.01 {
                        // 1% crosstalk threshold
                        problematic_pairs.push((i, j, max_crosstalk));
                    }
                }
            }
        }

        // Strategy 2: For high-crosstalk pairs, implement mitigation strategies
        for (gate1_idx, gate2_idx, crosstalk_level) in problematic_pairs {
            if crosstalk_level > 0.05 {
                // 5% threshold for aggressive mitigation
                // Insert timing delays to avoid simultaneous execution
                decisions.push(OptimizationDecision::GateReordering {
                    gates: vec![
                        original_gates[gate1_idx].name().to_string(),
                        original_gates[gate2_idx].name().to_string(),
                    ],
                    reason: format!(
                        "Avoid {:.1}% crosstalk by serializing execution",
                        crosstalk_level * 100.0
                    ),
                });
            } else if crosstalk_level > 0.02 {
                // 2% threshold for moderate mitigation
                // Try to remap one of the gates to a less problematic qubit
                let gate1_qubits = original_gates[gate1_idx].qubits();
                let gate2_qubits = original_gates[gate2_idx].qubits();

                // Look for alternative qubits with lower crosstalk
                if let Some(better_mapping) =
                    Self::find_lower_crosstalk_mapping(&gate1_qubits, &gate2_qubits, calibration)
                {
                    decisions.push(OptimizationDecision::QubitRemapping {
                        gate: original_gates[gate1_idx].name().to_string(),
                        original_qubits: gate1_qubits.clone(),
                        new_qubits: better_mapping,
                        reason: format!(
                            "Reduce crosstalk from {:.1}% to target <2%",
                            crosstalk_level * 100.0
                        ),
                    });
                }
            }
        }

        // Strategy 3: Apply echo sequences for Z-Z crosstalk mitigation
        for gate in original_gates {
            if gate.qubits().len() == 2 {
                let (q1, q2) = (gate.qubits()[0], gate.qubits()[1]);
                let q1_idx = q1.0 as usize;
                let q2_idx = q2.0 as usize;

                if q1_idx < crosstalk_matrix.matrix.len()
                    && q2_idx < crosstalk_matrix.matrix[q1_idx].len()
                {
                    let zz_crosstalk = crosstalk_matrix.matrix[q1_idx][q2_idx];

                    // For significant Z-Z crosstalk, suggest echo sequences
                    if zz_crosstalk > 0.001 && gate.name().contains("CZ") {
                        decisions.push(OptimizationDecision::GateSubstitution {
                            original: gate.name().to_string(),
                            replacement: "Echo_CZ".to_string(),
                            qubits: gate.qubits().clone(),
                            fidelity_change: zz_crosstalk * 0.8, // Echo reduces crosstalk by ~80%
                            duration_change: 50.0,               // Echo adds some overhead
                        });
                    }
                }
            }
        }

        // Strategy 4: Spectator qubit management
        // For idle qubits during two-qubit operations, apply dynamical decoupling
        let active_qubits = Self::get_active_qubits_per_layer(original_gates);

        for (layer_idx, layer_qubits) in active_qubits.iter().enumerate() {
            let all_qubits: HashSet<QubitId> = (0..calibration.topology.num_qubits)
                .map(|i| QubitId(i as u32))
                .collect();

            let layer_qubits_set: HashSet<QubitId> = layer_qubits.iter().copied().collect();
            let idle_qubits: Vec<QubitId> =
                all_qubits.difference(&layer_qubits_set).copied().collect();

            if !idle_qubits.is_empty() && layer_qubits.len() >= 2 {
                // Check if any idle qubits have significant crosstalk with active ones
                for &idle_qubit in &idle_qubits {
                    let mut max_crosstalk_to_active: f32 = 0.0;

                    for &active_qubit in layer_qubits {
                        let idle_idx = idle_qubit.0 as usize;
                        let active_idx = active_qubit.0 as usize;

                        if idle_idx < crosstalk_matrix.matrix.len()
                            && active_idx < crosstalk_matrix.matrix[idle_idx].len()
                        {
                            let crosstalk = crosstalk_matrix.matrix[idle_idx][active_idx] as f32;
                            max_crosstalk_to_active = max_crosstalk_to_active.max(crosstalk);
                        }
                    }

                    if max_crosstalk_to_active > 0.005_f32 {
                        // 0.5% threshold
                        decisions.push(OptimizationDecision::GateSubstitution {
                            original: "IDLE".to_string(),
                            replacement: "Dynamical_Decoupling".to_string(),
                            qubits: vec![idle_qubit],
                            fidelity_change: (max_crosstalk_to_active * 0.7) as f64, // DD reduces crosstalk
                            duration_change: 10.0, // Small overhead for DD pulses
                        });
                    }
                }
            }
        }

        Ok(())
    }

    /// Rank qubits by quality metrics
    fn rank_qubits_by_quality(calibration: &DeviceCalibration) -> Vec<(QubitId, f64)> {
        let mut qubit_scores = Vec::new();

        for (qubit_id, qubit_cal) in &calibration.qubit_calibrations {
            // Combine various metrics into a quality score
            let t1_score = qubit_cal.t1 / 100_000.0; // Normalize to ~1
            let t2_score = qubit_cal.t2 / 100_000.0;
            let readout_score = 1.0 - qubit_cal.readout_error;

            // Weight the scores (these weights could be configurable)
            let quality_score =
                0.4_f64.mul_add(readout_score, 0.3_f64.mul_add(t1_score, 0.3 * t2_score));

            qubit_scores.push((*qubit_id, quality_score));
        }

        // Sort by quality (highest first)
        qubit_scores.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(Ordering::Equal));

        qubit_scores
    }

    /// Estimate circuit fidelity based on calibration data
    fn estimate_circuit_fidelity<const N: usize>(
        circuit: &Circuit<N>,
        calibration: &DeviceCalibration,
    ) -> QuantRS2Result<f64> {
        let mut total_fidelity = 1.0;

        // Multiply fidelities of all gates (assumes independent errors)
        for gate in circuit.gates() {
            let gate_fidelity = Self::estimate_gate_fidelity(gate.as_ref(), calibration)?;
            total_fidelity *= gate_fidelity;
        }

        Ok(total_fidelity)
    }

    /// Estimate circuit duration based on calibration data
    fn estimate_circuit_duration<const N: usize>(
        circuit: &Circuit<N>,
        calibration: &DeviceCalibration,
    ) -> QuantRS2Result<f64> {
        // This would calculate critical path through the circuit
        // For now, return sum of gate durations (sequential execution)
        let mut total_duration = 0.0;

        for gate in circuit.gates() {
            let gate_duration = Self::estimate_gate_duration(gate.as_ref(), calibration)?;
            total_duration += gate_duration;
        }

        Ok(total_duration)
    }

    /// Estimate fidelity of a specific gate
    fn estimate_gate_fidelity(
        gate: &dyn GateOp,
        calibration: &DeviceCalibration,
    ) -> QuantRS2Result<f64> {
        let qubits = gate.qubits();

        match qubits.len() {
            1 => {
                // Single-qubit gate
                let qubit_id = qubits[0];
                if let Some(gate_cal) = calibration.single_qubit_gates.get(gate.name()) {
                    if let Some(qubit_data) = gate_cal.qubit_data.get(&qubit_id) {
                        return Ok(qubit_data.fidelity);
                    }
                }
                // Default single-qubit fidelity
                Ok(0.999)
            }
            2 => {
                // Two-qubit gate
                let qubit_pair = (qubits[0], qubits[1]);
                if let Some(gate_cal) = calibration.two_qubit_gates.get(&qubit_pair) {
                    return Ok(gate_cal.fidelity);
                }
                // Default two-qubit fidelity
                Ok(0.99)
            }
            _ => {
                // Multi-qubit gates have lower fidelity
                Ok(0.95)
            }
        }
    }

    /// Estimate duration of a specific gate
    fn estimate_gate_duration(
        gate: &dyn GateOp,
        calibration: &DeviceCalibration,
    ) -> QuantRS2Result<f64> {
        let qubits = gate.qubits();

        match qubits.len() {
            1 => {
                // Single-qubit gate
                let qubit_id = qubits[0];
                if let Some(gate_cal) = calibration.single_qubit_gates.get(gate.name()) {
                    if let Some(qubit_data) = gate_cal.qubit_data.get(&qubit_id) {
                        return Ok(qubit_data.duration);
                    }
                }
                // Default single-qubit duration (ns)
                Ok(30.0)
            }
            2 => {
                // Two-qubit gate
                let qubit_pair = (qubits[0], qubits[1]);
                if let Some(gate_cal) = calibration.two_qubit_gates.get(&qubit_pair) {
                    return Ok(gate_cal.duration);
                }
                // Default two-qubit duration (ns)
                Ok(300.0)
            }
            _ => {
                // Multi-qubit gates take longer
                Ok(1000.0)
            }
        }
    }

    /// Get active qubits for each layer of gates
    fn get_active_qubits_per_layer(
        gates: &[std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>],
    ) -> Vec<Vec<QubitId>> {
        let mut layers = Vec::new();
        let mut current_layer = Vec::new();
        let mut used_qubits = std::collections::HashSet::new();

        for gate in gates {
            let gate_qubits = gate.qubits();
            let has_conflict = gate_qubits.iter().any(|q| used_qubits.contains(q));

            if has_conflict {
                // Start new layer
                if !current_layer.is_empty() {
                    layers.push(current_layer);
                    current_layer = Vec::new();
                    used_qubits.clear();
                }
            }

            current_layer.extend(gate_qubits.clone());
            used_qubits.extend(gate_qubits.clone());
        }

        if !current_layer.is_empty() {
            layers.push(current_layer);
        }

        layers
    }

    /// Try to decompose a gate into simpler gates
    fn try_decompose_gate(
        _gate: &dyn quantrs2_core::gate::GateOp,
    ) -> Option<Vec<Box<dyn quantrs2_core::gate::GateOp>>> {
        // Placeholder implementation - in practice would decompose gates based on hardware constraints
        None
    }

    /// Find better qubit pair for two-qubit gate based on connectivity and error rates
    const fn find_better_qubit_pair(
        _current_pair: &(quantrs2_core::qubit::QubitId, quantrs2_core::qubit::QubitId),
        _calibration: &DeviceCalibration,
    ) -> Option<(quantrs2_core::qubit::QubitId, quantrs2_core::qubit::QubitId)> {
        // Placeholder implementation - would search for better connected qubits with lower error rates
        None
    }

    /// Prioritize native gates in the gate sequence
    fn prioritize_native_gates(
        _gates: &mut Vec<std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>>,
        _calibration: &DeviceCalibration,
        _decisions: &mut Vec<OptimizationDecision>,
    ) -> QuantRS2Result<()> {
        // Placeholder implementation - would reorder gates to prefer native operations
        Ok(())
    }

    /// Identify gates that can be executed in parallel
    fn identify_parallelizable_gates<const N: usize>(
        circuit: &Circuit<N>,
        _calibration: &DeviceCalibration,
    ) -> QuantRS2Result<Vec<Vec<std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>>>>
    {
        let mut parallel_groups = Vec::new();
        let gates = circuit.gates();

        // Simple dependency analysis - gates can be parallel if they don't share qubits
        let mut current_group = Vec::new();
        let mut used_qubits = HashSet::new();

        for gate in gates {
            let gate_qubits: HashSet<QubitId> = gate.qubits().into_iter().collect();

            // Check if this gate can be added to current group
            if used_qubits.is_disjoint(&gate_qubits) {
                current_group.push(gate.clone());
                used_qubits.extend(gate_qubits);
            } else {
                // Start new group
                if !current_group.is_empty() {
                    parallel_groups.push(current_group);
                }
                current_group = vec![gate.clone()];
                used_qubits = gate_qubits;
            }
        }

        if !current_group.is_empty() {
            parallel_groups.push(current_group);
        }

        Ok(parallel_groups)
    }

    /// Find faster alternatives for a gate
    fn find_faster_gate_alternatives(
        gate_name: &str,
        _qubit_data: &crate::calibration::SingleQubitGateData,
    ) -> Option<(String, f64)> {
        // Map of gate alternatives with typical speed improvements
        let alternatives = match gate_name {
            "RX" => vec![("Virtual_RX", 30.0), ("Composite_RX", 15.0)],
            "RY" => vec![("Physical_RY", 5.0)],
            "H" => vec![("Fast_H", 10.0)],
            _ => vec![],
        };

        // Return the first alternative that would be faster
        for (alt_name, improvement) in alternatives {
            if improvement > 5.0 {
                // Only if improvement is significant
                return Some((alt_name.to_string(), improvement));
            }
        }

        None
    }

    /// Remove redundant gates from circuit
    fn remove_redundant_gates<const N: usize>(
        circuit: &mut Circuit<N>,
        decisions: &mut Vec<OptimizationDecision>,
    ) -> QuantRS2Result<()> {
        // This would implement gate cancellation logic
        // e.g., H-H = I, X-X = I, consecutive rotations can be combined

        let gates = circuit.gates();
        let mut to_remove = Vec::new();

        // Look for consecutive identical Pauli gates
        for i in 0..(gates.len() - 1) {
            let gate1 = &gates[i];
            let gate2 = &gates[i + 1];

            if gate1.name() == gate2.name()
                && gate1.qubits() == gate2.qubits()
                && (gate1.name() == "X"
                    || gate1.name() == "Y"
                    || gate1.name() == "Z"
                    || gate1.name() == "H")
            {
                to_remove.push(i);
                to_remove.push(i + 1);

                decisions.push(OptimizationDecision::GateSubstitution {
                    original: format!("{}-{}", gate1.name(), gate2.name()),
                    replacement: "Identity".to_string(),
                    qubits: gate1.qubits().clone(),
                    fidelity_change: 0.001, // Removing gates improves fidelity
                    duration_change: -60.0, // Save two gate durations
                });
            }
        }

        Ok(())
    }

    /// Optimize gate scheduling based on hardware constraints
    fn optimize_gate_scheduling<const N: usize>(
        _circuit: &mut Circuit<N>,
        _calibration: &DeviceCalibration,
        decisions: &mut Vec<OptimizationDecision>,
    ) -> QuantRS2Result<()> {
        // This would implement sophisticated scheduling algorithms
        // considering hardware timing constraints, crosstalk, etc.

        decisions.push(OptimizationDecision::GateReordering {
            gates: vec!["Optimized".to_string(), "Schedule".to_string()],
            reason: "Hardware-aware gate scheduling applied".to_string(),
        });

        Ok(())
    }

    /// Find native hardware replacement for a gate
    fn find_native_replacement(gate_name: &str, calibration: &DeviceCalibration) -> Option<String> {
        // Map non-native gates to native equivalents
        let native_map = match gate_name {
            "T" => Some("RZ_pi_4"),      // T gate as Z rotation
            "S" => Some("RZ_pi_2"),      // S gate as Z rotation
            "SQRT_X" => Some("RX_pi_2"), // √X as X rotation
            "SQRT_Y" => Some("RY_pi_2"), // √Y as Y rotation
            _ => None,
        };

        if let Some(native_name) = native_map {
            // Check if the native gate is actually available in calibration
            if calibration.single_qubit_gates.contains_key(native_name) {
                return Some(native_name.to_string());
            }
        }

        None
    }

    /// Find composite gate decomposition
    fn find_composite_decomposition(gate_name: &str) -> Option<Vec<String>> {
        match gate_name {
            "TOFFOLI" => Some(vec![
                "H".to_string(),
                "CNOT".to_string(),
                "T".to_string(),
                "CNOT".to_string(),
                "T".to_string(),
                "H".to_string(),
            ]),
            "FREDKIN" => Some(vec![
                "CNOT".to_string(),
                "TOFFOLI".to_string(),
                "CNOT".to_string(),
            ]),
            _ => None,
        }
    }

    /// Find mapping with lower crosstalk
    const fn find_lower_crosstalk_mapping(
        _qubits1: &[QubitId],
        _qubits2: &[QubitId],
        _calibration: &DeviceCalibration,
    ) -> Option<Vec<QubitId>> {
        // This would search for alternative qubit mappings with lower crosstalk
        // For now, return None to indicate no better mapping found
        None
    }
}

/// Fidelity estimator for more sophisticated analysis
pub struct FidelityEstimator {
    /// Use process tomography data if available
    use_process_tomography: bool,
    /// Consider SPAM errors
    consider_spam_errors: bool,
    /// Model coherent errors
    model_coherent_errors: bool,
}

impl FidelityEstimator {
    /// Create a new fidelity estimator
    pub const fn new() -> Self {
        Self {
            use_process_tomography: false,
            consider_spam_errors: true,
            model_coherent_errors: true,
        }
    }

    /// Estimate process fidelity of a quantum circuit from device calibration.
    ///
    /// Computes the product of the per-operation fidelities looked up in
    /// `calibration`: every gate contributes its calibrated gate fidelity and
    /// (when [`Self::consider_spam_errors`] is set) every qubit touched by the
    /// circuit contributes its readout fidelity. This is the standard
    /// independent-error circuit-fidelity model `F ≈ Π_g F_g · Π_q F_readout(q)`.
    ///
    /// Returns an honest [`QuantRS2Error::InvalidInput`] when a gate in the
    /// circuit has no calibration entry, instead of fabricating a fidelity.
    pub fn estimate_process_fidelity<const N: usize>(
        &self,
        circuit: &Circuit<N>,
        calibration: &DeviceCalibration,
    ) -> QuantRS2Result<f64> {
        let mut fidelity = 1.0_f64;
        let mut touched_qubits: HashSet<QubitId> = HashSet::new();

        for gate in circuit.gates() {
            let qubits = gate.qubits();
            for q in &qubits {
                touched_qubits.insert(*q);
            }
            let gate_fidelity = match qubits.as_slice() {
                [q] => {
                    let gate_cal = calibration
                        .single_qubit_gates
                        .get(gate.name())
                        .and_then(|g| g.qubit_data.get(q));
                    match gate_cal {
                        Some(data) => data.fidelity,
                        None => {
                            return Err(QuantRS2Error::InvalidInput(format!(
                                "No single-qubit calibration for gate '{}' on qubit {}",
                                gate.name(),
                                q.id()
                            )));
                        }
                    }
                }
                [control, target] => {
                    let pair = calibration
                        .two_qubit_gates
                        .get(&(*control, *target))
                        .or_else(|| calibration.two_qubit_gates.get(&(*target, *control)));
                    match pair {
                        Some(data) => data.fidelity,
                        None => {
                            return Err(QuantRS2Error::InvalidInput(format!(
                                "No two-qubit calibration for gate '{}' on qubits ({}, {})",
                                gate.name(),
                                control.id(),
                                target.id()
                            )));
                        }
                    }
                }
                _ => {
                    return Err(QuantRS2Error::UnsupportedOperation(format!(
                        "Process-fidelity estimation supports 1- and 2-qubit gates only; \
                         gate '{}' acts on {} qubits",
                        gate.name(),
                        qubits.len()
                    )));
                }
            };
            fidelity *= gate_fidelity;
        }

        if self.consider_spam_errors {
            for q in &touched_qubits {
                if let Some(readout) = calibration.readout_calibration.qubit_readout.get(q) {
                    // Average assignment fidelity for the qubit.
                    let readout_fidelity = 0.5 * (readout.p0_given_0 + readout.p1_given_1);
                    fidelity *= readout_fidelity;
                }
            }
        }

        Ok(fidelity.clamp(0.0, 1.0))
    }

    /// Helper methods for optimization strategies
    /// Try to decompose a gate into simpler components
    fn try_decompose_gate(_gate: &dyn GateOp) -> Option<Vec<Box<dyn GateOp>>> {
        // This would implement gate decomposition logic
        // For now, return None to indicate no decomposition found
        None
    }

    /// Find a better qubit pair for a two-qubit gate
    fn find_better_qubit_pair(
        current_pair: &(QubitId, QubitId),
        calibration: &DeviceCalibration,
    ) -> Option<(QubitId, QubitId)> {
        let current_fidelity = calibration
            .two_qubit_gates
            .get(current_pair)
            .map_or(0.99, |g| g.fidelity);

        // Search for alternative qubit pairs with better fidelity
        for (&(q1, q2), gate_cal) in &calibration.two_qubit_gates {
            if (q1, q2) != *current_pair && gate_cal.fidelity > current_fidelity + 0.01 {
                return Some((q1, q2));
            }
        }

        None
    }

    /// Prioritize native gates in the gate sequence
    fn prioritize_native_gates(
        &self,
        gates: &mut Vec<std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>>,
        calibration: &DeviceCalibration,
        decisions: &mut Vec<OptimizationDecision>,
    ) -> QuantRS2Result<()> {
        // This would reorder gates to prioritize native hardware gates
        // Implementation depends on specific hardware capabilities
        Ok(())
    }

    /// Identify gates that can be executed in parallel
    fn identify_parallelizable_gates<const N: usize>(
        &self,
        circuit: &Circuit<N>,
        calibration: &DeviceCalibration,
    ) -> QuantRS2Result<Vec<Vec<std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>>>>
    {
        let mut parallel_groups = Vec::new();
        let gates = circuit.gates();

        // Simple dependency analysis - gates can be parallel if they don't share qubits
        let mut current_group = Vec::new();
        let mut used_qubits = HashSet::new();

        for gate in gates {
            let gate_qubits: HashSet<QubitId> = gate.qubits().into_iter().collect();

            // Check if this gate can be added to current group
            if used_qubits.is_disjoint(&gate_qubits) {
                current_group.push(gate.clone());
                used_qubits.extend(gate_qubits);
            } else {
                // Start new group
                if !current_group.is_empty() {
                    parallel_groups.push(current_group);
                }
                current_group = vec![gate.clone()];
                used_qubits = gate_qubits;
            }
        }

        if !current_group.is_empty() {
            parallel_groups.push(current_group);
        }

        Ok(parallel_groups)
    }

    /// Find faster alternatives for a gate
    fn find_faster_gate_alternatives(
        gate_name: &str,
        _qubit_data: &crate::calibration::SingleQubitGateData,
    ) -> Option<(String, f64)> {
        // Map of gate alternatives with typical speed improvements
        let alternatives = match gate_name {
            "RX" => vec![("Virtual_RX", 30.0), ("Composite_RX", 15.0)],
            "RY" => vec![("Physical_RY", 5.0)],
            "H" => vec![("Fast_H", 10.0)],
            _ => vec![],
        };

        // Return the first alternative that would be faster
        for (alt_name, improvement) in alternatives {
            if improvement > 5.0 {
                // Only if improvement is significant
                return Some((alt_name.to_string(), improvement));
            }
        }

        None
    }

    /// Remove redundant gates from circuit
    fn remove_redundant_gates<const N: usize>(
        circuit: &mut Circuit<N>,
        decisions: &mut Vec<OptimizationDecision>,
    ) -> QuantRS2Result<()> {
        // This would implement gate cancellation logic
        // e.g., H-H = I, X-X = I, consecutive rotations can be combined

        let gates = circuit.gates();
        let mut to_remove = Vec::new();

        // Look for consecutive identical Pauli gates
        for i in 0..(gates.len() - 1) {
            let gate1 = &gates[i];
            let gate2 = &gates[i + 1];

            if gate1.name() == gate2.name()
                && gate1.qubits() == gate2.qubits()
                && (gate1.name() == "X"
                    || gate1.name() == "Y"
                    || gate1.name() == "Z"
                    || gate1.name() == "H")
            {
                to_remove.push(i);
                to_remove.push(i + 1);

                decisions.push(OptimizationDecision::GateSubstitution {
                    original: format!("{}-{}", gate1.name(), gate2.name()),
                    replacement: "Identity".to_string(),
                    qubits: gate1.qubits().clone(),
                    fidelity_change: 0.001, // Removing gates improves fidelity
                    duration_change: -60.0, // Save two gate durations
                });
            }
        }

        Ok(())
    }

    /// Optimize gate scheduling based on hardware constraints
    fn optimize_gate_scheduling<const N: usize>(
        _circuit: &mut Circuit<N>,
        _calibration: &DeviceCalibration,
        decisions: &mut Vec<OptimizationDecision>,
    ) -> QuantRS2Result<()> {
        // This would implement sophisticated scheduling algorithms
        // considering hardware timing constraints, crosstalk, etc.

        decisions.push(OptimizationDecision::GateReordering {
            gates: vec!["Optimized".to_string(), "Schedule".to_string()],
            reason: "Hardware-aware gate scheduling applied".to_string(),
        });

        Ok(())
    }

    /// Find native hardware replacement for a gate
    fn find_native_replacement(gate_name: &str, calibration: &DeviceCalibration) -> Option<String> {
        // Map non-native gates to native equivalents
        let native_map = match gate_name {
            "T" => Some("RZ_pi_4"),      // T gate as Z rotation
            "S" => Some("RZ_pi_2"),      // S gate as Z rotation
            "SQRT_X" => Some("RX_pi_2"), // √X as X rotation
            "SQRT_Y" => Some("RY_pi_2"), // √Y as Y rotation
            _ => None,
        };

        if let Some(native_name) = native_map {
            // Check if the native gate is actually available in calibration
            if calibration.single_qubit_gates.contains_key(native_name) {
                return Some(native_name.to_string());
            }
        }

        None
    }

    /// Find composite gate decomposition
    fn find_composite_decomposition(gate_name: &str) -> Option<Vec<String>> {
        match gate_name {
            "TOFFOLI" => Some(vec![
                "H".to_string(),
                "CNOT".to_string(),
                "T".to_string(),
                "CNOT".to_string(),
                "T".to_string(),
                "H".to_string(),
            ]),
            "FREDKIN" => Some(vec![
                "CNOT".to_string(),
                "TOFFOLI".to_string(),
                "CNOT".to_string(),
            ]),
            _ => None,
        }
    }

    /// Get active qubits for each layer of gates
    fn get_active_qubits_per_layer(
        gates: &[std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>],
    ) -> Vec<HashSet<QubitId>> {
        let mut layers = Vec::new();
        let mut current_layer = HashSet::new();

        for gate in gates {
            let gate_qubits: HashSet<QubitId> = gate.qubits().into_iter().collect();

            if current_layer.is_disjoint(&gate_qubits) {
                current_layer.extend(gate_qubits);
            } else {
                if !current_layer.is_empty() {
                    layers.push(current_layer);
                }
                current_layer = gate_qubits;
            }
        }

        if !current_layer.is_empty() {
            layers.push(current_layer);
        }

        layers
    }
}

/// Pulse-level optimizer for fine-grained control
pub struct PulseOptimizer {
    /// Maximum pulse amplitude
    max_amplitude: f64,
    /// Pulse sample rate (GHz)
    sample_rate: f64,
    /// Use DRAG correction
    use_drag: bool,
}

impl PulseOptimizer {
    /// Create a new pulse optimizer
    pub const fn new() -> Self {
        Self {
            max_amplitude: 1.0,
            sample_rate: 4.5, // Typical for superconducting qubits
            use_drag: true,
        }
    }

    /// Generate an optimized, time-sampled in-phase pulse envelope for a
    /// single-qubit gate from its device calibration.
    ///
    /// The envelope is sampled at [`Self::sample_rate`] (GHz) over the gate's
    /// calibrated duration using the calibrated [`PulseShape`], scaled by the
    /// calibrated drive amplitude and clipped to [`Self::max_amplitude`]. When
    /// [`Self::use_drag`] is set and the shape carries a DRAG `beta`, the DRAG
    /// quadrature derivative is folded into the in-phase samples via its
    /// envelope-derivative correction.
    ///
    /// Returns an honest [`QuantRS2Error`] for gates that are not single-qubit
    /// or that lack a calibration entry, instead of returning an empty vector
    /// that silently pretends a pulse was produced.
    pub fn optimize_gate_pulses(
        &self,
        gate: &dyn GateOp,
        calibration: &DeviceCalibration,
    ) -> QuantRS2Result<Vec<f64>> {
        let qubits = gate.qubits();
        let qubit = match qubits.as_slice() {
            [q] => *q,
            _ => {
                return Err(QuantRS2Error::UnsupportedOperation(format!(
                    "PulseOptimizer.optimize_gate_pulses supports single-qubit gates only; \
                     gate '{}' acts on {} qubits",
                    gate.name(),
                    qubits.len()
                )));
            }
        };

        let gate_data = calibration
            .single_qubit_gates
            .get(gate.name())
            .and_then(|g| g.qubit_data.get(&qubit))
            .ok_or_else(|| {
                QuantRS2Error::InvalidInput(format!(
                    "No single-qubit calibration for gate '{}' on qubit {}",
                    gate.name(),
                    qubit.id()
                ))
            })?;

        // Number of samples = duration[ns] * sample_rate[GHz] (GHz = samples/ns).
        let num_samples = (gate_data.duration * self.sample_rate).round() as usize;
        if num_samples == 0 {
            // Virtual gates (e.g. RZ) have zero physical duration: no envelope.
            return Ok(Vec::new());
        }

        let amp = gate_data.amplitude;
        let mut samples = Vec::with_capacity(num_samples);
        let n = num_samples as f64;
        let center = (n - 1.0) / 2.0;

        for k in 0..num_samples {
            let t = k as f64;
            let value = match &gate_data.pulse_shape {
                PulseShape::Gaussian { sigma, .. } => gaussian_envelope(t, center, *sigma, amp),
                PulseShape::GaussianDRAG { sigma, beta, .. } => {
                    let base = gaussian_envelope(t, center, *sigma, amp);
                    if self.use_drag {
                        // DRAG in-phase remains the Gaussian; the derivative term
                        // modifies the quadrature, but its leakage-suppression
                        // correction back-acts on I as -beta * d/dt(envelope).
                        let derivative = gaussian_derivative(t, center, *sigma, amp);
                        beta.mul_add(-derivative, base)
                    } else {
                        base
                    }
                }
                PulseShape::Square { rise_time } => {
                    square_envelope(t, n, *rise_time, self.sample_rate, amp)
                }
                PulseShape::Cosine { rise_time } => {
                    cosine_envelope(t, n, *rise_time, self.sample_rate, amp)
                }
                PulseShape::Custom { name, .. } => {
                    return Err(QuantRS2Error::UnsupportedOperation(format!(
                        "PulseOptimizer cannot synthesise custom pulse shape '{name}' \
                         without its sample generator"
                    )));
                }
            };
            samples.push(value.clamp(-self.max_amplitude, self.max_amplitude));
        }

        Ok(samples)
    }
}

/// Normalised Gaussian envelope sample at time index `t` (in samples).
fn gaussian_envelope(t: f64, center: f64, sigma: f64, amplitude: f64) -> f64 {
    if sigma <= 0.0 {
        return 0.0;
    }
    let x = (t - center) / sigma;
    amplitude * (-0.5 * x * x).exp()
}

/// Time derivative (w.r.t. sample index) of [`gaussian_envelope`], used for the
/// DRAG quadrature correction.
fn gaussian_derivative(t: f64, center: f64, sigma: f64, amplitude: f64) -> f64 {
    if sigma <= 0.0 {
        return 0.0;
    }
    let x = t - center;
    -(x / (sigma * sigma)) * amplitude * (-0.5 * (x / sigma).powi(2)).exp()
}

/// Flat-top (square) envelope with cosine rise/fall of `rise_time` nanoseconds.
fn square_envelope(t: f64, n: f64, rise_time_ns: f64, sample_rate: f64, amplitude: f64) -> f64 {
    let rise_samples = (rise_time_ns * sample_rate).max(0.0);
    if rise_samples <= 0.0 {
        return amplitude;
    }
    if t < rise_samples {
        let phase = std::f64::consts::FRAC_PI_2 * (t / rise_samples);
        amplitude * phase.sin()
    } else if t > n - 1.0 - rise_samples {
        let phase = std::f64::consts::FRAC_PI_2 * ((n - 1.0 - t) / rise_samples);
        amplitude * phase.sin().max(0.0)
    } else {
        amplitude
    }
}

/// Raised-cosine envelope spanning the whole gate window.
fn cosine_envelope(t: f64, n: f64, _rise_time_ns: f64, _sample_rate: f64, amplitude: f64) -> f64 {
    if n <= 1.0 {
        return amplitude;
    }
    let phase = std::f64::consts::PI * t / (n - 1.0);
    amplitude * phase.sin()
}

#[cfg(test)]
mod tests {
    use super::*;
    use quantrs2_core::gate::single::Hadamard;

    #[test]
    fn test_optimization_config() {
        let config = OptimizationConfig::default();
        assert!(config.optimize_fidelity);
        assert!(config.optimize_duration);
    }

    #[test]
    fn test_calibration_optimizer() {
        let manager = CalibrationManager::new();
        let config = OptimizationConfig::default();
        let optimizer = CalibrationOptimizer::new(manager, config);

        // Create a simple test circuit
        let mut circuit = Circuit::<2>::new();
        let _ = circuit.h(QubitId(0));
        let _ = circuit.cnot(QubitId(0), QubitId(1));

        // Optimization should fail without calibration
        let result = optimizer.optimize_circuit(&circuit, "test_device");
        assert!(result.is_err());
    }

    #[test]
    fn test_fidelity_estimator() {
        use crate::calibration::create_ideal_calibration;

        let estimator = FidelityEstimator::new();
        let mut circuit = Circuit::<3>::new();
        let _ = circuit.h(QubitId(0));
        let _ = circuit.cnot(QubitId(0), QubitId(1));
        let _ = circuit.cnot(QubitId(1), QubitId(2));

        // The "ideal" calibration uses F_1q = 0.999, F_cnot = 0.99, and
        // readout p = 0.999, so the product is deterministic and must match the
        // independent-error model exactly (H + 2×CNOT, 3 qubits read out).
        let calibration = create_ideal_calibration("test_device".to_string(), 3);
        let fidelity = estimator
            .estimate_process_fidelity(&circuit, &calibration)
            .expect("estimate_process_fidelity should succeed for calibrated circuit");
        let expected = 0.999_f64 * 0.99 * 0.99 * 0.999_f64.powi(3);
        assert!(
            (fidelity - expected).abs() < 1e-12,
            "process fidelity {fidelity} should equal product model {expected}"
        );

        // A gate on an un-calibrated qubit pair must yield an honest error, not
        // a fabricated fidelity. The ideal calibration only covers the nearest-
        // neighbour pairs (0,1) and (1,2), so a CNOT on (0,2) has no entry.
        let mut uncalibrated = Circuit::<3>::new();
        let _ = uncalibrated.cnot(QubitId(0), QubitId(2));
        assert!(estimator
            .estimate_process_fidelity(&uncalibrated, &calibration)
            .is_err());
    }
}