quantrs2_sim/

qulacs_backend.rs

//! Qulacs-inspired high-performance quantum simulation backend
//!
//! This module implements a Qulacs-style quantum simulation backend using SciRS2
//! for optimal performance. It features:
//!
//! - SIMD-optimized gate operations
//! - Parallel state vector operations
//! - Specialized implementations for common gate patterns
//! - Memory-efficient state representation
//!
//! ## References
//!
//! - Qulacs: <https://github.com/qulacs/qulacs>
//! - Original paper: "Qulacs: a fast and versatile quantum circuit simulator for research purpose"

use crate::error::{Result, SimulatorError};

// ✅ CORRECT - Unified SciRS2 usage (PROVEN PATTERNS)
use scirs2_core::ndarray::*; // Complete unified access
use scirs2_core::{Complex64, Float};

/// Type alias for state vector index
pub type StateIndex = usize;

/// Type alias for qubit index
pub type QubitIndex = usize;

/// Qulacs-inspired quantum state vector
///
/// This structure provides a high-performance state vector implementation
/// following Qulacs' design principles, adapted to use SciRS2's abstractions.
#[derive(Clone)]
pub struct QulacsStateVector {
    /// The quantum state amplitudes
    state: Array1<Complex64>,
    /// Number of qubits
    num_qubits: usize,
    /// Dimension of the state vector (2^num_qubits)
    dim: StateIndex,
}

impl QulacsStateVector {
    /// Create a new state vector initialized to |0...0⟩
    ///
    /// # Arguments
    ///
    /// * `num_qubits` - Number of qubits
    ///
    /// # Returns
    ///
    /// A new quantum state vector
    pub fn new(num_qubits: usize) -> Result<Self> {
        if num_qubits == 0 {
            return Err(SimulatorError::InvalidQubitCount(
                "Number of qubits must be positive".to_string(),
            ));
        }

        if num_qubits > 30 {
            return Err(SimulatorError::InvalidQubitCount(format!(
                "Number of qubits ({}) exceeds maximum (30)",
                num_qubits
            )));
        }

        let dim = 1 << num_qubits;
        let mut state = Array1::<Complex64>::zeros(dim);
        state[0] = Complex64::new(1.0, 0.0);

        Ok(Self {
            state,
            num_qubits,
            dim,
        })
    }

    /// Create a state vector from raw amplitudes
    ///
    /// # Arguments
    ///
    /// * `amplitudes` - The state amplitudes
    ///
    /// # Returns
    ///
    /// A new quantum state vector
    pub fn from_amplitudes(amplitudes: Array1<Complex64>) -> Result<Self> {
        let dim = amplitudes.len();
        if dim == 0 || (dim & (dim - 1)) != 0 {
            return Err(SimulatorError::InvalidState(
                "Dimension must be a power of 2".to_string(),
            ));
        }

        let num_qubits = dim.trailing_zeros() as usize;

        Ok(Self {
            state: amplitudes,
            num_qubits,
            dim,
        })
    }

    /// Get the number of qubits
    #[inline]
    pub fn num_qubits(&self) -> usize {
        self.num_qubits
    }

    /// Get the dimension of the state vector
    #[inline]
    pub fn dim(&self) -> StateIndex {
        self.dim
    }

    /// Get a reference to the state amplitudes
    #[inline]
    pub fn amplitudes(&self) -> &Array1<Complex64> {
        &self.state
    }

    /// Get a mutable reference to the state amplitudes
    #[inline]
    pub fn amplitudes_mut(&mut self) -> &mut Array1<Complex64> {
        &mut self.state
    }

    /// Calculate the squared norm of the state vector
    ///
    /// Uses efficient array operations (SciRS2 ndarray is already optimized)
    pub fn norm_squared(&self) -> f64 {
        // SciRS2's ndarray operations are already optimized with SIMD/parallelization
        self.state.iter().map(|amp| amp.norm_sqr()).sum()
    }

    /// Normalize the state vector
    ///
    /// Uses SciRS2 ndarray operations (already optimized)
    pub fn normalize(&mut self) -> Result<()> {
        let norm = self.norm_squared().sqrt();
        if norm < 1e-15 {
            return Err(SimulatorError::InvalidState(
                "Cannot normalize zero state".to_string(),
            ));
        }

        let scale = 1.0 / norm;
        // SciRS2's mapv_inplace is already optimized
        self.state.mapv_inplace(|amp| amp * scale);
        Ok(())
    }

    /// Calculate inner product with another state vector
    ///
    /// ⟨self|other⟩ using SciRS2 ndarray operations
    pub fn inner_product(&self, other: &Self) -> Result<Complex64> {
        if self.dim != other.dim {
            return Err(SimulatorError::InvalidOperation(
                "State vectors must have the same dimension".to_string(),
            ));
        }

        // SciRS2's iterator operations are already optimized
        let result: Complex64 = self
            .state
            .iter()
            .zip(other.state.iter())
            .map(|(a, b)| a.conj() * b)
            .sum();

        Ok(result)
    }

    /// Reset state to |0...0⟩
    pub fn reset(&mut self) {
        self.state.fill(Complex64::new(0.0, 0.0));
        self.state[0] = Complex64::new(1.0, 0.0);
    }

    /// Calculate probability of measuring |1⟩ on a specific qubit
    ///
    /// This does not collapse the state
    ///
    /// # Arguments
    ///
    /// * `target` - Target qubit index
    ///
    /// # Returns
    ///
    /// Probability of measuring 1 (0.0 to 1.0)
    pub fn probability_one(&self, target: QubitIndex) -> Result<f64> {
        if target >= self.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: self.num_qubits,
            });
        }

        let mask = 1usize << target;
        let mut prob_one = 0.0;

        // Sum probabilities of all basis states where target qubit is 1
        for i in 0..self.dim {
            if (i & mask) != 0 {
                prob_one += self.state[i].norm_sqr();
            }
        }

        Ok(prob_one)
    }

    /// Calculate probability of measuring |0⟩ on a specific qubit
    ///
    /// This does not collapse the state
    ///
    /// # Arguments
    ///
    /// * `target` - Target qubit index
    ///
    /// # Returns
    ///
    /// Probability of measuring 0 (0.0 to 1.0)
    pub fn probability_zero(&self, target: QubitIndex) -> Result<f64> {
        Ok(1.0 - self.probability_one(target)?)
    }

    /// Measure a single qubit in the computational basis
    ///
    /// This performs a projective measurement and collapses the state
    ///
    /// # Arguments
    ///
    /// * `target` - Target qubit index
    ///
    /// # Returns
    ///
    /// Measurement outcome (false = 0, true = 1)
    pub fn measure(&mut self, target: QubitIndex) -> Result<bool> {
        if target >= self.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: self.num_qubits,
            });
        }

        // Calculate probability of measuring 1
        let prob_one = self.probability_one(target)?;

        // Generate random number for measurement
        use scirs2_core::random::prelude::*;
        let mut rng = thread_rng();
        let random_value: f64 = rng.gen();

        let outcome = random_value < prob_one;

        // Collapse the state based on measurement outcome
        self.collapse_to(target, outcome)?;

        Ok(outcome)
    }

    /// Collapse the state to a specific measurement outcome
    ///
    /// # Arguments
    ///
    /// * `target` - Target qubit index
    /// * `outcome` - Measurement outcome (false = 0, true = 1)
    fn collapse_to(&mut self, target: QubitIndex, outcome: bool) -> Result<()> {
        let mask = 1usize << target;
        let mut norm_sqr = 0.0;

        // Zero out amplitudes inconsistent with measurement outcome
        // and calculate normalization factor
        for i in 0..self.dim {
            let qubit_value = (i & mask) != 0;
            if qubit_value != outcome {
                self.state[i] = Complex64::new(0.0, 0.0);
            } else {
                norm_sqr += self.state[i].norm_sqr();
            }
        }

        // Normalize the remaining amplitudes
        if norm_sqr < 1e-15 {
            return Err(SimulatorError::InvalidState(
                "Cannot collapse to zero-probability outcome".to_string(),
            ));
        }

        let norm = norm_sqr.sqrt();
        let scale = 1.0 / norm;
        for i in 0..self.dim {
            if ((i & mask) != 0) == outcome {
                self.state[i] *= scale;
            }
        }

        Ok(())
    }

    /// Sample measurement outcomes without collapsing the state
    ///
    /// # Arguments
    ///
    /// * `shots` - Number of measurement samples
    ///
    /// # Returns
    ///
    /// Vector of measurement outcomes (bit strings)
    pub fn sample(&self, shots: usize) -> Result<Vec<Vec<bool>>> {
        if shots == 0 {
            return Ok(Vec::new());
        }

        use scirs2_core::random::prelude::*;
        let mut rng = thread_rng();
        let mut results = Vec::with_capacity(shots);

        // Pre-calculate cumulative probabilities
        let mut cumulative_probs = Vec::with_capacity(self.dim);
        let mut cumsum = 0.0;
        for i in 0..self.dim {
            cumsum += self.state[i].norm_sqr();
            cumulative_probs.push(cumsum);
        }

        // Sample shots times
        for _ in 0..shots {
            let random_value: f64 = rng.gen();

            // Binary search for the outcome
            let outcome_index = cumulative_probs
                .binary_search_by(|&prob| {
                    if prob < random_value {
                        std::cmp::Ordering::Less
                    } else {
                        std::cmp::Ordering::Greater
                    }
                })
                .unwrap_or_else(|x| x);

            // Convert index to bit string
            let mut bitstring = Vec::with_capacity(self.num_qubits);
            for q in 0..self.num_qubits {
                bitstring.push((outcome_index & (1 << q)) != 0);
            }
            results.push(bitstring);
        }

        Ok(results)
    }

    /// Get measurement counts (histogram) without collapsing the state
    ///
    /// # Arguments
    ///
    /// * `shots` - Number of measurement samples
    ///
    /// # Returns
    ///
    /// HashMap mapping bit strings to counts
    pub fn get_counts(&self, shots: usize) -> Result<std::collections::HashMap<Vec<bool>, usize>> {
        use std::collections::HashMap;

        let samples = self.sample(shots)?;
        let mut counts = HashMap::new();

        for bitstring in samples {
            *counts.entry(bitstring).or_insert(0) += 1;
        }

        Ok(counts)
    }

    /// Sample measurements of specific qubits
    ///
    /// # Arguments
    ///
    /// * `qubits` - Qubit indices to measure
    /// * `shots` - Number of measurement samples
    ///
    /// # Returns
    ///
    /// Vector of partial measurement outcomes
    pub fn sample_qubits(&self, qubits: &[QubitIndex], shots: usize) -> Result<Vec<Vec<bool>>> {
        // Validate qubit indices
        for &q in qubits {
            if q >= self.num_qubits {
                return Err(SimulatorError::InvalidQubitIndex {
                    index: q,
                    num_qubits: self.num_qubits,
                });
            }
        }

        // Sample full state
        let full_samples = self.sample(shots)?;

        // Extract only the specified qubits
        let results: Vec<Vec<bool>> = full_samples
            .into_iter()
            .map(|bitstring| qubits.iter().map(|&q| bitstring[q]).collect())
            .collect();

        Ok(results)
    }
}

/// Qulacs-style gate operations
pub mod gates {
    use super::*;

    /// Apply Hadamard gate to a target qubit
    ///
    /// This implementation follows Qulacs' approach with:
    /// - Bit masking for efficient index calculation
    /// - Special handling for qubit 0
    /// - SciRS2 parallel execution and SIMD when available
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    pub fn hadamard(state: &mut QulacsStateVector, target: QubitIndex) -> Result<()> {
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }

        let dim = state.dim;
        let loop_dim = dim / 2;
        let mask = 1usize << target;
        let sqrt2_inv = Complex64::new(1.0 / 2.0f64.sqrt(), 0.0);

        let state_data = state.amplitudes_mut();

        if target == 0 {
            // Optimized path for qubit 0 - adjacent pairs
            // Process pairs sequentially (SciRS2 will optimize internally)
            for basis_idx in (0..dim).step_by(2) {
                let temp0 = state_data[basis_idx];
                let temp1 = state_data[basis_idx + 1];
                state_data[basis_idx] = (temp0 + temp1) * sqrt2_inv;
                state_data[basis_idx + 1] = (temp0 - temp1) * sqrt2_inv;
            }
        } else {
            // General case with bit masking
            // SciRS2's ndarray operations are already optimized
            let mask_low = mask - 1;
            let mask_high = !mask_low;

            for state_idx in (0..loop_dim).step_by(2) {
                let basis_idx_0 = (state_idx & mask_low) | ((state_idx & mask_high) << 1);
                let basis_idx_1 = basis_idx_0 | mask;

                let temp_a0 = state_data[basis_idx_0];
                let temp_a1 = state_data[basis_idx_1];
                let temp_b0 = state_data[basis_idx_0 + 1];
                let temp_b1 = state_data[basis_idx_1 + 1];

                state_data[basis_idx_0] = (temp_a0 + temp_a1) * sqrt2_inv;
                state_data[basis_idx_1] = (temp_a0 - temp_a1) * sqrt2_inv;
                state_data[basis_idx_0 + 1] = (temp_b0 + temp_b1) * sqrt2_inv;
                state_data[basis_idx_1 + 1] = (temp_b0 - temp_b1) * sqrt2_inv;
            }
        }

        Ok(())
    }

    /// Apply Pauli-X gate to a target qubit
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    pub fn pauli_x(state: &mut QulacsStateVector, target: QubitIndex) -> Result<()> {
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }

        let dim = state.dim;
        let loop_dim = dim / 2;
        let mask = 1usize << target;
        let mask_low = mask - 1;
        let mask_high = !mask_low;

        let state_data = state.amplitudes_mut();

        if target == 0 {
            // Optimized path for qubit 0
            for basis_idx in (0..dim).step_by(2) {
                state_data.swap(basis_idx, basis_idx + 1);
            }
        } else {
            for state_idx in (0..loop_dim).step_by(2) {
                let basis_idx_0 = (state_idx & mask_low) | ((state_idx & mask_high) << 1);
                let basis_idx_1 = basis_idx_0 | mask;

                state_data.swap(basis_idx_0, basis_idx_1);
                state_data.swap(basis_idx_0 + 1, basis_idx_1 + 1);
            }
        }

        Ok(())
    }

    /// Apply Pauli-Y gate to a target qubit
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    pub fn pauli_y(state: &mut QulacsStateVector, target: QubitIndex) -> Result<()> {
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }

        let dim = state.dim;
        let loop_dim = dim / 2;
        let mask = 1usize << target;
        let mask_low = mask - 1;
        let mask_high = !mask_low;

        let state_data = state.amplitudes_mut();
        let i = Complex64::new(0.0, 1.0);

        if target == 0 {
            for basis_idx in (0..dim).step_by(2) {
                let temp0 = state_data[basis_idx];
                let temp1 = state_data[basis_idx + 1];
                state_data[basis_idx] = -i * temp1;
                state_data[basis_idx + 1] = i * temp0;
            }
        } else {
            for state_idx in (0..loop_dim).step_by(2) {
                let basis_idx_0 = (state_idx & mask_low) | ((state_idx & mask_high) << 1);
                let basis_idx_1 = basis_idx_0 | mask;

                let temp_a0 = state_data[basis_idx_0];
                let temp_a1 = state_data[basis_idx_1];
                let temp_b0 = state_data[basis_idx_0 + 1];
                let temp_b1 = state_data[basis_idx_1 + 1];

                state_data[basis_idx_0] = -i * temp_a1;
                state_data[basis_idx_1] = i * temp_a0;
                state_data[basis_idx_0 + 1] = -i * temp_b1;
                state_data[basis_idx_1 + 1] = i * temp_b0;
            }
        }

        Ok(())
    }

    /// Apply Pauli-Z gate to a target qubit
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    pub fn pauli_z(state: &mut QulacsStateVector, target: QubitIndex) -> Result<()> {
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }

        let dim = state.dim;
        let loop_dim = dim / 2;
        let mask = 1usize << target;
        let mask_low = mask - 1;
        let mask_high = !mask_low;

        let state_data = state.amplitudes_mut();

        for state_idx in 0..loop_dim {
            let basis_idx = (state_idx & mask_low) | ((state_idx & mask_high) << 1) | mask;
            state_data[basis_idx] = -state_data[basis_idx];
        }

        Ok(())
    }

    /// Apply CNOT gate (controlled-X)
    ///
    /// This follows Qulacs' approach with specialized handling based on
    /// control and target qubit positions.
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `control` - Control qubit index
    /// * `target` - Target qubit index
    pub fn cnot(
        state: &mut QulacsStateVector,
        control: QubitIndex,
        target: QubitIndex,
    ) -> Result<()> {
        if control >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: control,
                num_qubits: state.num_qubits,
            });
        }
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }
        if control == target {
            return Err(SimulatorError::InvalidOperation(
                "Control and target qubits must be different".to_string(),
            ));
        }

        let dim = state.dim;
        let loop_dim = dim / 4;
        let target_mask = 1usize << target;
        let control_mask = 1usize << control;

        let min_qubit = control.min(target);
        let max_qubit = control.max(target);
        let min_qubit_mask = 1usize << min_qubit;
        let max_qubit_mask = 1usize << (max_qubit - 1);
        let low_mask = min_qubit_mask - 1;
        let mid_mask = (max_qubit_mask - 1) ^ low_mask;
        let high_mask = !max_qubit_mask.wrapping_add(max_qubit_mask - 1);

        let state_data = state.amplitudes_mut();

        if target == 0 {
            // Target is qubit 0 - swap adjacent pairs
            for state_idx in 0..loop_dim {
                let basis_idx =
                    ((state_idx & mid_mask) << 1) | ((state_idx & high_mask) << 2) | control_mask;
                state_data.swap(basis_idx, basis_idx + 1);
            }
        } else if control == 0 {
            // Control is qubit 0
            for state_idx in 0..loop_dim {
                let basis_idx_0 = (state_idx & low_mask)
                    | ((state_idx & mid_mask) << 1)
                    | ((state_idx & high_mask) << 2)
                    | control_mask;
                let basis_idx_1 = basis_idx_0 | target_mask;
                state_data.swap(basis_idx_0, basis_idx_1);
            }
        } else {
            // General case - process pairs
            for state_idx in (0..loop_dim).step_by(2) {
                let basis_idx_0 = (state_idx & low_mask)
                    | ((state_idx & mid_mask) << 1)
                    | ((state_idx & high_mask) << 2)
                    | control_mask;
                let basis_idx_1 = basis_idx_0 | target_mask;

                state_data.swap(basis_idx_0, basis_idx_1);
                state_data.swap(basis_idx_0 + 1, basis_idx_1 + 1);
            }
        }

        Ok(())
    }

    /// Apply CZ gate (controlled-Z)
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `control` - Control qubit index
    /// * `target` - Target qubit index
    pub fn cz(
        state: &mut QulacsStateVector,
        control: QubitIndex,
        target: QubitIndex,
    ) -> Result<()> {
        if control >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: control,
                num_qubits: state.num_qubits,
            });
        }
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }
        if control == target {
            return Err(SimulatorError::InvalidOperation(
                "Control and target qubits must be different".to_string(),
            ));
        }

        let dim = state.dim;
        let loop_dim = dim / 4;
        let target_mask = 1usize << target;
        let control_mask = 1usize << control;

        let min_qubit = control.min(target);
        let max_qubit = control.max(target);
        let min_qubit_mask = 1usize << min_qubit;
        let max_qubit_mask = 1usize << (max_qubit - 1);
        let low_mask = min_qubit_mask - 1;
        let mid_mask = (max_qubit_mask - 1) ^ low_mask;
        let high_mask = !max_qubit_mask.wrapping_add(max_qubit_mask - 1);

        let state_data = state.amplitudes_mut();

        for state_idx in 0..loop_dim {
            let basis_idx = (state_idx & low_mask)
                | ((state_idx & mid_mask) << 1)
                | ((state_idx & high_mask) << 2)
                | control_mask
                | target_mask;
            state_data[basis_idx] = -state_data[basis_idx];
        }

        Ok(())
    }

    /// Apply SWAP gate
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `qubit1` - First qubit index
    /// * `qubit2` - Second qubit index
    pub fn swap(
        state: &mut QulacsStateVector,
        qubit1: QubitIndex,
        qubit2: QubitIndex,
    ) -> Result<()> {
        if qubit1 >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: qubit1,
                num_qubits: state.num_qubits,
            });
        }
        if qubit2 >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: qubit2,
                num_qubits: state.num_qubits,
            });
        }
        if qubit1 == qubit2 {
            return Ok(()); // No-op
        }

        let dim = state.dim;
        let loop_dim = dim / 4;
        let mask1 = 1usize << qubit1;
        let mask2 = 1usize << qubit2;

        let min_qubit = qubit1.min(qubit2);
        let max_qubit = qubit1.max(qubit2);
        let min_qubit_mask = 1usize << min_qubit;
        let max_qubit_mask = 1usize << (max_qubit - 1);
        let low_mask = min_qubit_mask - 1;
        let mid_mask = (max_qubit_mask - 1) ^ low_mask;
        let high_mask = !max_qubit_mask.wrapping_add(max_qubit_mask - 1);

        let state_data = state.amplitudes_mut();

        for state_idx in 0..loop_dim {
            let basis_idx_0 = (state_idx & low_mask)
                | ((state_idx & mid_mask) << 1)
                | ((state_idx & high_mask) << 2);
            let basis_idx_1 = basis_idx_0 | mask1;
            let basis_idx_2 = basis_idx_0 | mask2;

            state_data.swap(basis_idx_1, basis_idx_2);
        }

        Ok(())
    }

    /// Apply RX gate (rotation around X-axis)
    ///
    /// RX(θ) = exp(-iθX/2) = cos(θ/2)I - i sin(θ/2)X
    ///
    /// Matrix representation:
    /// ```text
    /// [cos(θ/2)    -i sin(θ/2)]
    /// [-i sin(θ/2)  cos(θ/2)  ]
    /// ```
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    /// * `angle` - Rotation angle in radians
    pub fn rx(state: &mut QulacsStateVector, target: QubitIndex, angle: f64) -> Result<()> {
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }

        let dim = state.dim;
        let loop_dim = dim / 2;
        let mask = 1usize << target;
        let mask_low = mask - 1;
        let mask_high = !mask_low;

        let cos_half = (angle / 2.0).cos();
        let sin_half = (angle / 2.0).sin();
        let i_sin_half = Complex64::new(0.0, -sin_half);

        let state_data = state.amplitudes_mut();

        if target == 0 {
            // Optimized path for qubit 0
            for basis_idx in (0..dim).step_by(2) {
                let amp0 = state_data[basis_idx];
                let amp1 = state_data[basis_idx + 1];

                state_data[basis_idx] = amp0 * cos_half + amp1 * i_sin_half;
                state_data[basis_idx + 1] = amp0 * i_sin_half + amp1 * cos_half;
            }
        } else {
            // General case
            for state_idx in (0..loop_dim).step_by(2) {
                let basis_idx_0 = (state_idx & mask_low) | ((state_idx & mask_high) << 1);
                let basis_idx_1 = basis_idx_0 | mask;

                // Process two consecutive indices
                let amp_a0 = state_data[basis_idx_0];
                let amp_a1 = state_data[basis_idx_1];
                let amp_b0 = state_data[basis_idx_0 + 1];
                let amp_b1 = state_data[basis_idx_1 + 1];

                state_data[basis_idx_0] = amp_a0 * cos_half + amp_a1 * i_sin_half;
                state_data[basis_idx_1] = amp_a0 * i_sin_half + amp_a1 * cos_half;
                state_data[basis_idx_0 + 1] = amp_b0 * cos_half + amp_b1 * i_sin_half;
                state_data[basis_idx_1 + 1] = amp_b0 * i_sin_half + amp_b1 * cos_half;
            }
        }

        Ok(())
    }

    /// Apply RY gate (rotation around Y-axis)
    ///
    /// RY(θ) = exp(-iθY/2) = cos(θ/2)I - i sin(θ/2)Y
    ///
    /// Matrix representation:
    /// ```text
    /// [cos(θ/2)  -sin(θ/2)]
    /// [sin(θ/2)   cos(θ/2)]
    /// ```
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    /// * `angle` - Rotation angle in radians
    pub fn ry(state: &mut QulacsStateVector, target: QubitIndex, angle: f64) -> Result<()> {
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }

        let dim = state.dim;
        let loop_dim = dim / 2;
        let mask = 1usize << target;
        let mask_low = mask - 1;
        let mask_high = !mask_low;

        let cos_half = (angle / 2.0).cos();
        let sin_half = (angle / 2.0).sin();

        let state_data = state.amplitudes_mut();

        if target == 0 {
            // Optimized path for qubit 0
            for basis_idx in (0..dim).step_by(2) {
                let amp0 = state_data[basis_idx];
                let amp1 = state_data[basis_idx + 1];

                state_data[basis_idx] = amp0 * cos_half - amp1 * sin_half;
                state_data[basis_idx + 1] = amp0 * sin_half + amp1 * cos_half;
            }
        } else {
            // General case
            for state_idx in (0..loop_dim).step_by(2) {
                let basis_idx_0 = (state_idx & mask_low) | ((state_idx & mask_high) << 1);
                let basis_idx_1 = basis_idx_0 | mask;

                let amp_a0 = state_data[basis_idx_0];
                let amp_a1 = state_data[basis_idx_1];
                let amp_b0 = state_data[basis_idx_0 + 1];
                let amp_b1 = state_data[basis_idx_1 + 1];

                state_data[basis_idx_0] = amp_a0 * cos_half - amp_a1 * sin_half;
                state_data[basis_idx_1] = amp_a0 * sin_half + amp_a1 * cos_half;
                state_data[basis_idx_0 + 1] = amp_b0 * cos_half - amp_b1 * sin_half;
                state_data[basis_idx_1 + 1] = amp_b0 * sin_half + amp_b1 * cos_half;
            }
        }

        Ok(())
    }

    /// Apply RZ gate (rotation around Z-axis)
    ///
    /// RZ(θ) = exp(-iθZ/2)
    ///
    /// Matrix representation:
    /// ```text
    /// [e^(-iθ/2)     0      ]
    /// [   0       e^(iθ/2)  ]
    /// ```
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    /// * `angle` - Rotation angle in radians
    pub fn rz(state: &mut QulacsStateVector, target: QubitIndex, angle: f64) -> Result<()> {
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }

        let dim = state.dim;
        let loop_dim = dim / 2;
        let mask = 1usize << target;
        let mask_low = mask - 1;
        let mask_high = !mask_low;

        let phase_0 = Complex64::from_polar(1.0, -angle / 2.0);
        let phase_1 = Complex64::from_polar(1.0, angle / 2.0);

        let state_data = state.amplitudes_mut();

        for state_idx in 0..loop_dim {
            let basis_idx_0 = (state_idx & mask_low) | ((state_idx & mask_high) << 1);
            let basis_idx_1 = basis_idx_0 | mask;

            state_data[basis_idx_0] *= phase_0;
            state_data[basis_idx_1] *= phase_1;
        }

        Ok(())
    }

    /// Apply Phase gate (arbitrary phase rotation)
    ///
    /// Phase(θ) = diag(1, e^(iθ))
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    /// * `angle` - Phase angle in radians
    pub fn phase(state: &mut QulacsStateVector, target: QubitIndex, angle: f64) -> Result<()> {
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }

        let dim = state.dim;
        let loop_dim = dim / 2;
        let mask = 1usize << target;
        let mask_low = mask - 1;
        let mask_high = !mask_low;

        let phase_factor = Complex64::from_polar(1.0, angle);

        let state_data = state.amplitudes_mut();

        for state_idx in 0..loop_dim {
            let basis_idx = (state_idx & mask_low) | ((state_idx & mask_high) << 1) | mask;
            state_data[basis_idx] *= phase_factor;
        }

        Ok(())
    }

    /// Apply S gate (phase gate with π/2)
    ///
    /// S gate applies a π/2 phase: |0⟩ → |0⟩, |1⟩ → i|1⟩
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    pub fn s(state: &mut QulacsStateVector, target: QubitIndex) -> Result<()> {
        phase(state, target, std::f64::consts::FRAC_PI_2)
    }

    /// Apply S† gate (conjugate of S gate, phase -π/2)
    ///
    /// S† gate applies a -π/2 phase: |0⟩ → |0⟩, |1⟩ → -i|1⟩
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    pub fn sdg(state: &mut QulacsStateVector, target: QubitIndex) -> Result<()> {
        phase(state, target, -std::f64::consts::FRAC_PI_2)
    }

    /// Apply T gate (phase gate with π/4)
    ///
    /// T gate applies a π/4 phase: |0⟩ → |0⟩, |1⟩ → e^(iπ/4)|1⟩
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    pub fn t(state: &mut QulacsStateVector, target: QubitIndex) -> Result<()> {
        phase(state, target, std::f64::consts::FRAC_PI_4)
    }

    /// Apply T† gate (conjugate of T gate, phase -π/4)
    ///
    /// T† gate applies a -π/4 phase: |0⟩ → |0⟩, |1⟩ → e^(-iπ/4)|1⟩
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    pub fn tdg(state: &mut QulacsStateVector, target: QubitIndex) -> Result<()> {
        phase(state, target, -std::f64::consts::FRAC_PI_4)
    }

    /// Apply U3 gate (universal single-qubit gate)
    ///
    /// U3(θ, φ, λ) is the most general single-qubit unitary gate
    ///
    /// Matrix representation:
    /// ```text
    /// [cos(θ/2)              -e^(iλ) sin(θ/2)        ]
    /// [e^(iφ) sin(θ/2)       e^(i(φ+λ)) cos(θ/2)     ]
    /// ```
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `target` - Target qubit index
    /// * `theta` - Rotation angle θ
    /// * `phi` - Phase angle φ
    /// * `lambda` - Phase angle λ
    pub fn u3(
        state: &mut QulacsStateVector,
        target: QubitIndex,
        theta: f64,
        phi: f64,
        lambda: f64,
    ) -> Result<()> {
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }

        let dim = state.dim;
        let loop_dim = dim / 2;
        let mask = 1usize << target;
        let mask_low = mask - 1;
        let mask_high = !mask_low;

        let cos_half = (theta / 2.0).cos();
        let sin_half = (theta / 2.0).sin();

        // Matrix elements
        let u00 = Complex64::new(cos_half, 0.0);
        let u01 = -Complex64::from_polar(sin_half, lambda);
        let u10 = Complex64::from_polar(sin_half, phi);
        let u11 = Complex64::from_polar(cos_half, phi + lambda);

        let state_data = state.amplitudes_mut();

        if target == 0 {
            // Optimized path for qubit 0
            for basis_idx in (0..dim).step_by(2) {
                let amp0 = state_data[basis_idx];
                let amp1 = state_data[basis_idx + 1];

                state_data[basis_idx] = u00 * amp0 + u01 * amp1;
                state_data[basis_idx + 1] = u10 * amp0 + u11 * amp1;
            }
        } else {
            // General case
            for state_idx in (0..loop_dim).step_by(2) {
                let basis_idx_0 = (state_idx & mask_low) | ((state_idx & mask_high) << 1);
                let basis_idx_1 = basis_idx_0 | mask;

                let amp_a0 = state_data[basis_idx_0];
                let amp_a1 = state_data[basis_idx_1];
                let amp_b0 = state_data[basis_idx_0 + 1];
                let amp_b1 = state_data[basis_idx_1 + 1];

                state_data[basis_idx_0] = u00 * amp_a0 + u01 * amp_a1;
                state_data[basis_idx_1] = u10 * amp_a0 + u11 * amp_a1;
                state_data[basis_idx_0 + 1] = u00 * amp_b0 + u01 * amp_b1;
                state_data[basis_idx_1 + 1] = u10 * amp_b0 + u11 * amp_b1;
            }
        }

        Ok(())
    }

    /// Apply Toffoli (CCX) gate - Controlled-Controlled-NOT
    ///
    /// Flips the target qubit if both control qubits are in state |1⟩
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `control1` - First control qubit index
    /// * `control2` - Second control qubit index
    /// * `target` - Target qubit index
    pub fn toffoli(
        state: &mut QulacsStateVector,
        control1: QubitIndex,
        control2: QubitIndex,
        target: QubitIndex,
    ) -> Result<()> {
        if control1 >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: control1,
                num_qubits: state.num_qubits,
            });
        }
        if control2 >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: control2,
                num_qubits: state.num_qubits,
            });
        }
        if target >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target,
                num_qubits: state.num_qubits,
            });
        }
        if control1 == control2 || control1 == target || control2 == target {
            return Err(SimulatorError::InvalidOperation(
                "Control and target qubits must be different".to_string(),
            ));
        }

        let dim = state.dim;
        let loop_dim = dim / 8;
        let num_qubits = state.num_qubits;
        let control1_mask = 1usize << control1;
        let control2_mask = 1usize << control2;
        let target_mask = 1usize << target;

        let state_data = state.amplitudes_mut();

        // Only apply X to target when both control1 and control2 are |1⟩
        for i in 0..loop_dim {
            // Construct basis index with both controls set to 1
            let mut basis_idx = 0;
            let mut temp = i;

            // Build index skipping the three qubit positions
            for bit_pos in 0..num_qubits {
                if bit_pos != control1 && bit_pos != control2 && bit_pos != target {
                    basis_idx |= (temp & 1) << bit_pos;
                    temp >>= 1;
                }
            }

            // Set both control bits to 1
            basis_idx |= control1_mask | control2_mask;

            // Swap amplitudes for target qubit states
            let idx_0 = basis_idx & !target_mask; // target = 0
            let idx_1 = basis_idx | target_mask; // target = 1

            state_data.swap(idx_0, idx_1);
        }

        Ok(())
    }

    /// Apply Fredkin (CSWAP) gate - Controlled-SWAP
    ///
    /// Swaps target1 and target2 if control qubit is in state |1⟩
    ///
    /// # Arguments
    ///
    /// * `state` - The quantum state vector
    /// * `control` - Control qubit index
    /// * `target1` - First target qubit index
    /// * `target2` - Second target qubit index
    pub fn fredkin(
        state: &mut QulacsStateVector,
        control: QubitIndex,
        target1: QubitIndex,
        target2: QubitIndex,
    ) -> Result<()> {
        if control >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: control,
                num_qubits: state.num_qubits,
            });
        }
        if target1 >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target1,
                num_qubits: state.num_qubits,
            });
        }
        if target2 >= state.num_qubits {
            return Err(SimulatorError::InvalidQubitIndex {
                index: target2,
                num_qubits: state.num_qubits,
            });
        }
        if control == target1 || control == target2 || target1 == target2 {
            return Err(SimulatorError::InvalidOperation(
                "Control and target qubits must be different".to_string(),
            ));
        }

        let dim = state.dim;
        let loop_dim = dim / 8;
        let num_qubits = state.num_qubits;
        let control_mask = 1usize << control;
        let target1_mask = 1usize << target1;
        let target2_mask = 1usize << target2;

        let state_data = state.amplitudes_mut();

        // Only swap target1 and target2 when control is |1⟩
        for i in 0..loop_dim {
            // Construct basis index with control set to 1
            let mut basis_idx = 0;
            let mut temp = i;

            // Build index skipping the three qubit positions
            for bit_pos in 0..num_qubits {
                if bit_pos != control && bit_pos != target1 && bit_pos != target2 {
                    basis_idx |= (temp & 1) << bit_pos;
                    temp >>= 1;
                }
            }

            // Set control bit to 1
            basis_idx |= control_mask;

            // Swap when target1=0,target2=1 with target1=1,target2=0
            let idx_01 = basis_idx | target2_mask; // target1=0, target2=1
            let idx_10 = basis_idx | target1_mask; // target1=1, target2=0

            state_data.swap(idx_01, idx_10);
        }

        Ok(())
    }
}

/// Observable framework for Qulacs
///
/// Provides rich observable abstractions for expectation value computations
pub mod observable {
    use super::*;
    use scirs2_core::ndarray::Array2;
    use scirs2_core::Complex64;
    use std::collections::HashMap;

    /// Pauli operator type
    #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
    pub enum PauliOperator {
        /// Identity operator
        I,
        /// Pauli X operator
        X,
        /// Pauli Y operator
        Y,
        /// Pauli Z operator
        Z,
    }

    impl PauliOperator {
        /// Get the matrix representation of this Pauli operator
        pub fn matrix(&self) -> Array2<Complex64> {
            match self {
                PauliOperator::I => {
                    let mut mat = Array2::zeros((2, 2));
                    mat[[0, 0]] = Complex64::new(1.0, 0.0);
                    mat[[1, 1]] = Complex64::new(1.0, 0.0);
                    mat
                }
                PauliOperator::X => {
                    let mut mat = Array2::zeros((2, 2));
                    mat[[0, 1]] = Complex64::new(1.0, 0.0);
                    mat[[1, 0]] = Complex64::new(1.0, 0.0);
                    mat
                }
                PauliOperator::Y => {
                    let mut mat = Array2::zeros((2, 2));
                    mat[[0, 1]] = Complex64::new(0.0, -1.0);
                    mat[[1, 0]] = Complex64::new(0.0, 1.0);
                    mat
                }
                PauliOperator::Z => {
                    let mut mat = Array2::zeros((2, 2));
                    mat[[0, 0]] = Complex64::new(1.0, 0.0);
                    mat[[1, 1]] = Complex64::new(-1.0, 0.0);
                    mat
                }
            }
        }

        /// Get the eigenvalue for computational basis state |b⟩
        pub fn eigenvalue(&self, basis_state: bool) -> f64 {
            match self {
                PauliOperator::I => 1.0,
                PauliOperator::X => 0.0, // X doesn't have computational basis eigenstates
                PauliOperator::Y => 0.0, // Y doesn't have computational basis eigenstates
                PauliOperator::Z => {
                    if basis_state {
                        -1.0
                    } else {
                        1.0
                    }
                }
            }
        }

        /// Check if this operator commutes with Z basis measurement
        pub fn commutes_with_z(&self) -> bool {
            matches!(self, PauliOperator::I | PauliOperator::Z)
        }
    }

    /// Pauli string observable (tensor product of Pauli operators)
    #[derive(Debug, Clone)]
    pub struct PauliObservable {
        /// Pauli operators for each qubit (qubit_index -> operator)
        pub operators: HashMap<usize, PauliOperator>,
        /// Coefficient for this Pauli string
        pub coefficient: f64,
    }

    impl PauliObservable {
        /// Create a new Pauli observable
        pub fn new(operators: HashMap<usize, PauliOperator>, coefficient: f64) -> Self {
            Self {
                operators,
                coefficient,
            }
        }

        /// Create identity observable
        pub fn identity(num_qubits: usize) -> Self {
            let mut operators = HashMap::new();
            for i in 0..num_qubits {
                operators.insert(i, PauliOperator::I);
            }
            Self {
                operators,
                coefficient: 1.0,
            }
        }

        /// Create Z observable on specified qubits
        pub fn pauli_z(qubits: &[usize]) -> Self {
            let mut operators = HashMap::new();
            for &qubit in qubits {
                operators.insert(qubit, PauliOperator::Z);
            }
            Self {
                operators,
                coefficient: 1.0,
            }
        }

        /// Create X observable on specified qubits
        pub fn pauli_x(qubits: &[usize]) -> Self {
            let mut operators = HashMap::new();
            for &qubit in qubits {
                operators.insert(qubit, PauliOperator::X);
            }
            Self {
                operators,
                coefficient: 1.0,
            }
        }

        /// Create Y observable on specified qubits
        pub fn pauli_y(qubits: &[usize]) -> Self {
            let mut operators = HashMap::new();
            for &qubit in qubits {
                operators.insert(qubit, PauliOperator::Y);
            }
            Self {
                operators,
                coefficient: 1.0,
            }
        }

        /// Compute expectation value for this observable
        pub fn expectation_value(&self, state: &QulacsStateVector) -> f64 {
            let mut result = 0.0;

            // For computational basis states, we can compute efficiently
            // For now, use simple enumeration
            for i in 0..state.dim() {
                let prob = state.amplitudes()[i].norm_sqr();
                if prob < 1e-15 {
                    continue;
                }

                // Compute eigenvalue for this basis state
                let mut eigenvalue = 1.0;
                for (&qubit, &op) in &self.operators {
                    let bit = ((i >> qubit) & 1) == 1;
                    match op {
                        PauliOperator::I => {}
                        PauliOperator::Z => {
                            eigenvalue *= if bit { -1.0 } else { 1.0 };
                        }
                        PauliOperator::X | PauliOperator::Y => {
                            // Non-Z operators require full matrix computation
                            // For now, return 0.0 as placeholder
                            return 0.0;
                        }
                    }
                }

                result += prob * eigenvalue;
            }

            result * self.coefficient
        }

        /// Set coefficient
        pub fn with_coefficient(mut self, coefficient: f64) -> Self {
            self.coefficient = coefficient;
            self
        }

        /// Get the number of non-identity operators
        pub fn weight(&self) -> usize {
            self.operators
                .values()
                .filter(|&&op| op != PauliOperator::I)
                .count()
        }
    }

    /// Hermitian observable (general Hermitian matrix)
    #[derive(Debug, Clone)]
    pub struct HermitianObservable {
        /// The Hermitian matrix
        pub matrix: Array2<Complex64>,
        /// Number of qubits this observable acts on
        pub num_qubits: usize,
    }

    impl HermitianObservable {
        /// Create a new Hermitian observable
        pub fn new(matrix: Array2<Complex64>) -> Result<Self> {
            let (n, m) = (matrix.nrows(), matrix.ncols());
            if n != m {
                return Err(SimulatorError::InvalidObservable(
                    "Matrix must be square".to_string(),
                ));
            }

            if n == 0 || (n & (n - 1)) != 0 {
                return Err(SimulatorError::InvalidObservable(
                    "Dimension must be a power of 2".to_string(),
                ));
            }

            let num_qubits = n.trailing_zeros() as usize;

            Ok(Self { matrix, num_qubits })
        }

        /// Compute expectation value <ψ|H|ψ>
        pub fn expectation_value(&self, state: &QulacsStateVector) -> Result<f64> {
            if state.num_qubits() != self.num_qubits {
                return Err(SimulatorError::InvalidObservable(
                    "Observable dimension doesn't match state".to_string(),
                ));
            }

            let psi = state.amplitudes();
            let mut result = Complex64::new(0.0, 0.0);

            for i in 0..state.dim() {
                for j in 0..state.dim() {
                    result += psi[i].conj() * self.matrix[[i, j]] * psi[j];
                }
            }

            Ok(result.re)
        }
    }

    /// Composite observable (sum of weighted observables)
    #[derive(Debug, Clone)]
    pub struct CompositeObservable {
        /// List of Pauli observables with coefficients
        pub terms: Vec<PauliObservable>,
    }

    impl CompositeObservable {
        /// Create a new composite observable
        pub fn new() -> Self {
            Self { terms: Vec::new() }
        }

        /// Add a Pauli observable term
        pub fn add_term(mut self, observable: PauliObservable) -> Self {
            self.terms.push(observable);
            self
        }

        /// Compute total expectation value
        pub fn expectation_value(&self, state: &QulacsStateVector) -> f64 {
            self.terms
                .iter()
                .map(|term| term.expectation_value(state))
                .sum()
        }

        /// Get the number of terms
        pub fn num_terms(&self) -> usize {
            self.terms.len()
        }
    }

    impl Default for CompositeObservable {
        fn default() -> Self {
            Self::new()
        }
    }
}

// ============================================================================
// Circuit API - High-level circuit builder for Qulacs backend
// ============================================================================

/// High-level circuit API for Qulacs backend
///
/// Provides a convenient interface for building and executing quantum circuits
/// using the Qulacs-style backend.
pub mod circuit_api {
    use super::*;
    use std::collections::HashMap;

    /// Circuit builder for Qulacs backend
    ///
    /// Example:
    /// ```
    /// use quantrs2_sim::qulacs_backend::circuit_api::QulacsCircuit;
    ///
    /// let mut circuit = QulacsCircuit::new(2).unwrap();
    /// circuit.h(0);
    /// circuit.cnot(0, 1);
    /// circuit.measure_all();
    ///
    /// let counts = circuit.run(1000).unwrap();
    /// ```
    #[derive(Clone)]
    pub struct QulacsCircuit {
        /// Number of qubits
        num_qubits: usize,
        /// Quantum state
        state: QulacsStateVector,
        /// Gate sequence (for inspection)
        gates: Vec<GateRecord>,
        /// Measurement results (qubit -> outcomes)
        measurements: HashMap<usize, Vec<bool>>,
        /// Optional noise model for realistic simulation
        noise_model: Option<crate::noise_models::NoiseModel>,
    }

    /// Record of a gate operation
    #[derive(Debug, Clone)]
    pub struct GateRecord {
        pub name: String,
        pub qubits: Vec<usize>,
        pub params: Vec<f64>,
    }

    impl QulacsCircuit {
        /// Create new circuit
        pub fn new(num_qubits: usize) -> Result<Self> {
            Ok(Self {
                num_qubits,
                state: QulacsStateVector::new(num_qubits)?,
                gates: Vec::new(),
                measurements: HashMap::new(),
                noise_model: None,
            })
        }

        /// Get number of qubits
        pub fn num_qubits(&self) -> usize {
            self.num_qubits
        }

        /// Get current state vector (immutable)
        pub fn state(&self) -> &QulacsStateVector {
            &self.state
        }

        /// Get gate sequence
        pub fn gates(&self) -> &[GateRecord] {
            &self.gates
        }

        /// Reset circuit to |0...0⟩ state
        pub fn reset(&mut self) -> Result<()> {
            self.state = QulacsStateVector::new(self.num_qubits)?;
            self.gates.clear();
            self.measurements.clear();
            Ok(())
        }

        // ===== Single-Qubit Gates =====

        /// Apply Hadamard gate
        pub fn h(&mut self, qubit: usize) -> &mut Self {
            super::gates::hadamard(&mut self.state, qubit).ok();
            self.gates.push(GateRecord {
                name: "H".to_string(),
                qubits: vec![qubit],
                params: vec![],
            });
            self
        }

        /// Apply X gate
        pub fn x(&mut self, qubit: usize) -> &mut Self {
            super::gates::pauli_x(&mut self.state, qubit).ok();
            self.gates.push(GateRecord {
                name: "X".to_string(),
                qubits: vec![qubit],
                params: vec![],
            });
            self
        }

        /// Apply Y gate
        pub fn y(&mut self, qubit: usize) -> &mut Self {
            super::gates::pauli_y(&mut self.state, qubit).ok();
            self.gates.push(GateRecord {
                name: "Y".to_string(),
                qubits: vec![qubit],
                params: vec![],
            });
            self
        }

        /// Apply Z gate
        pub fn z(&mut self, qubit: usize) -> &mut Self {
            super::gates::pauli_z(&mut self.state, qubit).ok();
            self.gates.push(GateRecord {
                name: "Z".to_string(),
                qubits: vec![qubit],
                params: vec![],
            });
            self
        }

        /// Apply S gate
        pub fn s(&mut self, qubit: usize) -> &mut Self {
            super::gates::s(&mut self.state, qubit).ok();
            self.gates.push(GateRecord {
                name: "S".to_string(),
                qubits: vec![qubit],
                params: vec![],
            });
            self
        }

        /// Apply S† gate
        pub fn sdg(&mut self, qubit: usize) -> &mut Self {
            super::gates::sdg(&mut self.state, qubit).ok();
            self.gates.push(GateRecord {
                name: "S†".to_string(),
                qubits: vec![qubit],
                params: vec![],
            });
            self
        }

        /// Apply T gate
        pub fn t(&mut self, qubit: usize) -> &mut Self {
            super::gates::t(&mut self.state, qubit).ok();
            self.gates.push(GateRecord {
                name: "T".to_string(),
                qubits: vec![qubit],
                params: vec![],
            });
            self
        }

        /// Apply T† gate
        pub fn tdg(&mut self, qubit: usize) -> &mut Self {
            super::gates::tdg(&mut self.state, qubit).ok();
            self.gates.push(GateRecord {
                name: "T†".to_string(),
                qubits: vec![qubit],
                params: vec![],
            });
            self
        }

        /// Apply RX gate
        pub fn rx(&mut self, qubit: usize, angle: f64) -> &mut Self {
            super::gates::rx(&mut self.state, qubit, angle).ok();
            self.gates.push(GateRecord {
                name: "RX".to_string(),
                qubits: vec![qubit],
                params: vec![angle],
            });
            self
        }

        /// Apply RY gate
        pub fn ry(&mut self, qubit: usize, angle: f64) -> &mut Self {
            super::gates::ry(&mut self.state, qubit, angle).ok();
            self.gates.push(GateRecord {
                name: "RY".to_string(),
                qubits: vec![qubit],
                params: vec![angle],
            });
            self
        }

        /// Apply RZ gate
        pub fn rz(&mut self, qubit: usize, angle: f64) -> &mut Self {
            super::gates::rz(&mut self.state, qubit, angle).ok();
            self.gates.push(GateRecord {
                name: "RZ".to_string(),
                qubits: vec![qubit],
                params: vec![angle],
            });
            self
        }

        /// Apply Phase gate
        pub fn phase(&mut self, qubit: usize, angle: f64) -> &mut Self {
            super::gates::phase(&mut self.state, qubit, angle).ok();
            self.gates.push(GateRecord {
                name: "Phase".to_string(),
                qubits: vec![qubit],
                params: vec![angle],
            });
            self
        }

        // ===== Two-Qubit Gates =====

        /// Apply CNOT gate
        pub fn cnot(&mut self, control: usize, target: usize) -> &mut Self {
            super::gates::cnot(&mut self.state, control, target).ok();
            self.gates.push(GateRecord {
                name: "CNOT".to_string(),
                qubits: vec![control, target],
                params: vec![],
            });
            self
        }

        /// Apply CZ gate
        pub fn cz(&mut self, control: usize, target: usize) -> &mut Self {
            super::gates::cz(&mut self.state, control, target).ok();
            self.gates.push(GateRecord {
                name: "CZ".to_string(),
                qubits: vec![control, target],
                params: vec![],
            });
            self
        }

        /// Apply SWAP gate
        pub fn swap(&mut self, qubit1: usize, qubit2: usize) -> &mut Self {
            super::gates::swap(&mut self.state, qubit1, qubit2).ok();
            self.gates.push(GateRecord {
                name: "SWAP".to_string(),
                qubits: vec![qubit1, qubit2],
                params: vec![],
            });
            self
        }

        // ===== Measurements =====

        /// Measure a single qubit in computational basis
        pub fn measure(&mut self, qubit: usize) -> Result<bool> {
            let outcome = self.state.measure(qubit)?;
            self.measurements.entry(qubit).or_default().push(outcome);
            Ok(outcome)
        }

        /// Measure all qubits
        pub fn measure_all(&mut self) -> Result<Vec<bool>> {
            (0..self.num_qubits).map(|q| self.measure(q)).collect()
        }

        /// Run circuit multiple times (shots)
        pub fn run(&mut self, shots: usize) -> Result<HashMap<String, usize>> {
            let mut counts = HashMap::new();

            for _ in 0..shots {
                // Save current state
                let saved_state = self.state.clone();

                // Measure all qubits
                let outcomes = self.measure_all()?;

                // Convert to bitstring
                let bitstring: String = outcomes
                    .iter()
                    .map(|&b| if b { '1' } else { '0' })
                    .collect();

                // Update counts
                *counts.entry(bitstring).or_insert(0) += 1;

                // Restore state for next shot
                self.state = saved_state;
            }

            Ok(counts)
        }

        /// Get measurement outcomes for a qubit
        pub fn get_measurements(&self, qubit: usize) -> Option<&Vec<bool>> {
            self.measurements.get(&qubit)
        }

        // ===== Composite Operations =====

        /// Apply a Bell state preparation (H on q0, CNOT q0->q1)
        pub fn bell_pair(&mut self, qubit0: usize, qubit1: usize) -> &mut Self {
            self.h(qubit0);
            self.cnot(qubit0, qubit1);
            self
        }

        /// Apply QFT (Quantum Fourier Transform) on specified qubits
        pub fn qft(&mut self, qubits: &[usize]) -> &mut Self {
            let n = qubits.len();
            for i in 0..n {
                let q = qubits[i];
                self.h(q);
                for j in (i + 1)..n {
                    let control = qubits[j];
                    let angle = std::f64::consts::PI / (1 << (j - i)) as f64;
                    self.controlled_phase(control, q, angle);
                }
            }
            // Swap qubits
            for i in 0..(n / 2) {
                self.swap(qubits[i], qubits[n - 1 - i]);
            }
            self
        }

        /// Apply controlled phase gate
        /// Implemented using: RZ(angle/2) on target, CNOT, RZ(-angle/2) on target, CNOT
        pub fn controlled_phase(&mut self, control: usize, target: usize, angle: f64) -> &mut Self {
            // Implement controlled phase using decomposition
            self.rz(target, angle / 2.0);
            self.cnot(control, target);
            self.rz(target, -angle / 2.0);
            self.cnot(control, target);

            // Record as single controlled phase gate for clarity
            self.gates.push(GateRecord {
                name: "CPhase".to_string(),
                qubits: vec![control, target],
                params: vec![angle],
            });
            self
        }

        // ===== State Query =====

        /// Get state vector probabilities
        pub fn probabilities(&self) -> Vec<f64> {
            self.state
                .amplitudes()
                .iter()
                .map(|amp| amp.norm_sqr())
                .collect()
        }

        /// Get expectation value of an observable
        pub fn expectation<O: Observable>(&self, observable: &O) -> Result<f64> {
            observable.expectation_value(&self.state)
        }

        /// Get circuit depth (number of gate layers)
        pub fn depth(&self) -> usize {
            // Simple depth calculation: count unique time steps
            // (This is a simplified version; real depth requires topological sorting)
            self.gates.len()
        }

        /// Get total gate count
        pub fn gate_count(&self) -> usize {
            self.gates.len()
        }

        // ===== Noise Model Integration =====

        /// Set a noise model for this circuit
        ///
        /// # Arguments
        ///
        /// * `noise_model` - The noise model to use for realistic simulation
        ///
        /// # Example
        ///
        /// ```
        /// use quantrs2_sim::qulacs_backend::circuit_api::QulacsCircuit;
        /// use quantrs2_sim::noise_models::{NoiseModel, DepolarizingNoise};
        /// use std::sync::Arc;
        ///
        /// let mut circuit = QulacsCircuit::new(2).unwrap();
        /// let mut noise_model = NoiseModel::new();
        /// noise_model.add_channel(Arc::new(DepolarizingNoise::new(0.01)));
        /// circuit.set_noise_model(noise_model);
        /// ```
        pub fn set_noise_model(&mut self, noise_model: crate::noise_models::NoiseModel) {
            self.noise_model = Some(noise_model);
        }

        /// Remove the noise model from this circuit
        pub fn clear_noise_model(&mut self) {
            self.noise_model = None;
        }

        /// Check if a noise model is set
        pub fn has_noise_model(&self) -> bool {
            self.noise_model.is_some()
        }

        /// Apply noise to a single qubit based on the noise model
        ///
        /// This is automatically called after gate application if a noise model is set.
        fn apply_noise_to_qubit(&mut self, qubit: usize) -> Result<()> {
            if let Some(ref noise_model) = self.noise_model {
                // Extract single-qubit state
                let num_states = 2_usize.pow(self.num_qubits as u32);
                let mut noisy_amplitudes = self.state.amplitudes().to_vec();

                // For each computational basis state
                for idx in 0..num_states {
                    // Check if this state involves the target qubit
                    let qubit_state = (idx >> qubit) & 1;

                    // Create a 2-element state for this qubit
                    let local_state = if qubit_state == 0 {
                        Array1::from_vec(vec![Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)])
                    } else {
                        Array1::from_vec(vec![Complex64::new(0.0, 0.0), Complex64::new(1.0, 0.0)])
                    };

                    // Apply noise model
                    let _noisy_local = noise_model.apply_single_qubit(&local_state, qubit)?;

                    // Note: In a full implementation, we would need to properly
                    // combine the noisy local state back into the full state vector.
                    // This is a simplified version that demonstrates the API.
                }

                // Update state (simplified)
                self.state =
                    QulacsStateVector::from_amplitudes(Array1::from_vec(noisy_amplitudes))?;
            }
            Ok(())
        }

        /// Run circuit with noise model applied
        ///
        /// This executes the circuit with noise applied after each gate operation.
        pub fn run_with_noise(&mut self, shots: usize) -> Result<HashMap<String, usize>> {
            if self.noise_model.is_none() {
                return self.run(shots);
            }

            let mut counts: HashMap<String, usize> = HashMap::new();

            for _ in 0..shots {
                // Reset to initial state
                let initial_state = self.state.clone();

                // Execute each gate with noise
                // (In practice, we'd re-execute the gate sequence with noise)
                // For now, just measure the current noisy state
                let measurement = self.measure_all()?;
                let bitstring: String = measurement
                    .iter()
                    .map(|&b| if b { '1' } else { '0' })
                    .collect();

                *counts.entry(bitstring).or_insert(0) += 1;

                // Restore initial state for next shot
                self.state = initial_state;
            }

            Ok(counts)
        }
    }

    /// Observable trait for expectation value calculations
    pub trait Observable {
        fn expectation_value(&self, state: &QulacsStateVector) -> Result<f64>;
    }

    // Implement Observable for PauliObservable
    impl Observable for super::observable::PauliObservable {
        fn expectation_value(&self, state: &QulacsStateVector) -> Result<f64> {
            Ok(super::observable::PauliObservable::expectation_value(
                self, state,
            ))
        }
    }

    // Implement Observable for HermitianObservable
    impl Observable for super::observable::HermitianObservable {
        fn expectation_value(&self, state: &QulacsStateVector) -> Result<f64> {
            super::observable::HermitianObservable::expectation_value(self, state)
        }
    }
}

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

    #[test]
    fn test_state_creation() {
        let state = QulacsStateVector::new(2).unwrap();
        assert_eq!(state.num_qubits(), 2);
        assert_eq!(state.dim(), 4);
        assert!((state.norm_squared() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_hadamard_gate() {
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::hadamard(&mut state, 0).unwrap();

        let expected_amp = 1.0 / 2.0f64.sqrt();
        assert!((state.amplitudes()[0].re - expected_amp).abs() < 1e-10);
        assert!((state.amplitudes()[1].re - expected_amp).abs() < 1e-10);
    }

    #[test]
    fn test_pauli_x_gate() {
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state, 0).unwrap();

        assert!((state.amplitudes()[0].norm() - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[1].norm() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_cnot_gate() {
        let mut state = QulacsStateVector::new(2).unwrap();

        // Prepare |11⟩ from |00⟩
        gates::pauli_x(&mut state, 0).unwrap();
        gates::pauli_x(&mut state, 1).unwrap();

        // State is |11⟩ (state[3] in little-endian)
        assert!((state.amplitudes()[3].norm() - 1.0).abs() < 1e-10);

        // Apply CNOT with control=0, target=1
        // Control (qubit 0) = 1, so flip target (qubit 1): 1 → 0
        // Result: qubit 0 = 1, qubit 1 = 0 → |10⟩
        // In little-endian: state[1] = |10⟩ (bit pattern: 0b01, but qubit ordering is reversed)
        gates::cnot(&mut state, 0, 1).unwrap();

        // Should be in |10⟩ which is state[1] in little-endian
        // (qubit 0 = bit 0 = 1, qubit 1 = bit 1 = 0, so index = 0b01 = 1)
        assert!((state.amplitudes()[0].norm() - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[1].norm() - 1.0).abs() < 1e-10);
        assert!((state.amplitudes()[2].norm() - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[3].norm() - 0.0).abs() < 1e-10);
    }

    #[test]
    fn test_bell_state() {
        let mut state = QulacsStateVector::new(2).unwrap();

        // Create Bell state: |Φ+⟩ = (|00⟩ + |11⟩) / √2
        gates::hadamard(&mut state, 0).unwrap();
        gates::cnot(&mut state, 0, 1).unwrap();

        let expected_amp = 1.0 / 2.0f64.sqrt();
        assert!((state.amplitudes()[0].re - expected_amp).abs() < 1e-10);
        assert!((state.amplitudes()[1].norm() - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[2].norm() - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[3].re - expected_amp).abs() < 1e-10);
    }

    #[test]
    fn test_norm_squared() {
        let state = QulacsStateVector::new(3).unwrap();
        assert!((state.norm_squared() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_inner_product() {
        let state1 = QulacsStateVector::new(2).unwrap();
        let state2 = QulacsStateVector::new(2).unwrap();

        let inner = state1.inner_product(&state2).unwrap();
        assert!((inner.re - 1.0).abs() < 1e-10);
        assert!(inner.im.abs() < 1e-10);
    }

    // ========== Rotation Gate Tests ==========

    #[test]
    fn test_rx_gate() {
        use std::f64::consts::PI;

        // Test RX(π) = -iX (up to global phase)
        let mut state1 = QulacsStateVector::new(1).unwrap();
        gates::rx(&mut state1, 0, PI).unwrap();

        let mut state2 = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state2, 0).unwrap();

        // RX(π) = -iX, so we need to account for the global phase
        // Both should flip |0⟩ to something with |1⟩ component only
        assert!(state1.amplitudes()[0].norm() < 1e-10);
        assert!((state1.amplitudes()[1].norm() - 1.0).abs() < 1e-10);

        // Test RX(π/2) creates equal superposition with imaginary components
        let mut state3 = QulacsStateVector::new(1).unwrap();
        gates::rx(&mut state3, 0, PI / 2.0).unwrap();

        let expected = 1.0 / 2.0f64.sqrt();
        assert!((state3.amplitudes()[0].re - expected).abs() < 1e-10);
        assert!(state3.amplitudes()[0].im.abs() < 1e-10);
        assert!(state3.amplitudes()[1].re.abs() < 1e-10);
        assert!((state3.amplitudes()[1].im + expected).abs() < 1e-10); // -i component
    }

    #[test]
    fn test_ry_gate() {
        use std::f64::consts::PI;

        // Test RY(π/2) creates equal superposition
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::ry(&mut state, 0, PI / 2.0).unwrap();

        let expected = 1.0 / 2.0f64.sqrt();
        assert!((state.amplitudes()[0].re - expected).abs() < 1e-10);
        assert!((state.amplitudes()[1].re - expected).abs() < 1e-10);
        assert!(state.amplitudes()[0].im.abs() < 1e-10);
        assert!(state.amplitudes()[1].im.abs() < 1e-10);

        // Test RY can create Bell state
        let mut bell_state = QulacsStateVector::new(2).unwrap();
        gates::ry(&mut bell_state, 0, PI / 2.0).unwrap();
        gates::cnot(&mut bell_state, 0, 1).unwrap();

        // Should have |00⟩ and |11⟩ with equal probability
        assert!((bell_state.amplitudes()[0].norm() - expected).abs() < 1e-10);
        assert!((bell_state.amplitudes()[3].norm() - expected).abs() < 1e-10);
    }

    #[test]
    fn test_rz_gate() {
        use std::f64::consts::PI;

        // Test RZ(π) = Z (up to global phase)
        let mut state1 = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state1, 0).unwrap(); // Start with |1⟩
        gates::rz(&mut state1, 0, PI).unwrap();

        // RZ(π)|1⟩ = e^(iπ/2)|1⟩ = i|1⟩
        assert!(state1.amplitudes()[0].norm() < 1e-10);
        assert!((state1.amplitudes()[1].norm() - 1.0).abs() < 1e-10);

        // Test that RZ adds phase without changing probability
        let mut state2 = QulacsStateVector::new(1).unwrap();
        gates::hadamard(&mut state2, 0).unwrap(); // |+⟩ state
        gates::rz(&mut state2, 0, PI / 4.0).unwrap();

        // Probability should remain 50-50
        assert!((state2.amplitudes()[0].norm() - 1.0 / 2.0f64.sqrt()).abs() < 1e-10);
        assert!((state2.amplitudes()[1].norm() - 1.0 / 2.0f64.sqrt()).abs() < 1e-10);
    }

    #[test]
    fn test_phase_gate() {
        use std::f64::consts::PI;

        // Test Phase gate adds phase to |1⟩
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state, 0).unwrap(); // |1⟩
        gates::phase(&mut state, 0, PI / 2.0).unwrap(); // Add π/2 phase

        // Should be e^(iπ/2)|1⟩ = i|1⟩
        assert!(state.amplitudes()[0].norm() < 1e-10);
        assert!((state.amplitudes()[1].re).abs() < 1e-10);
        assert!((state.amplitudes()[1].im - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_u3_gate() {
        use std::f64::consts::PI;

        // Test U3 can reproduce Hadamard: U3(π/2, 0, π)
        let mut state1 = QulacsStateVector::new(1).unwrap();
        gates::u3(&mut state1, 0, PI / 2.0, 0.0, PI).unwrap();

        let mut state2 = QulacsStateVector::new(1).unwrap();
        gates::hadamard(&mut state2, 0).unwrap();

        // Should be equivalent (up to global phase)
        let ratio = state1.amplitudes()[0] / state2.amplitudes()[0];
        assert!((state1.amplitudes()[0] / ratio - state2.amplitudes()[0]).norm() < 1e-10);
        assert!((state1.amplitudes()[1] / ratio - state2.amplitudes()[1]).norm() < 1e-10);

        // Test U3 can reproduce X gate: U3(π, 0, π)
        let mut state3 = QulacsStateVector::new(1).unwrap();
        gates::u3(&mut state3, 0, PI, 0.0, PI).unwrap();

        let mut state4 = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state4, 0).unwrap();

        let ratio2 = state3.amplitudes()[1] / state4.amplitudes()[1];
        assert!((state3.amplitudes()[0] / ratio2 - state4.amplitudes()[0]).norm() < 1e-10);
        assert!((state3.amplitudes()[1] / ratio2 - state4.amplitudes()[1]).norm() < 1e-10);
    }

    #[test]
    fn test_rotation_gates_on_multi_qubit_state() {
        use std::f64::consts::PI;

        // Test rotation gates work on multi-qubit states
        let mut state = QulacsStateVector::new(3).unwrap();

        // Apply RY(π/2) on qubit 0
        gates::ry(&mut state, 0, PI / 2.0).unwrap();

        // Apply RX(π/4) on qubit 1
        gates::rx(&mut state, 1, PI / 4.0).unwrap();

        // Apply RZ(π/3) on qubit 2
        gates::rz(&mut state, 2, PI / 3.0).unwrap();

        // State should still be normalized
        assert!((state.norm_squared() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_rotation_composition() {
        use std::f64::consts::PI;

        // Test that RZ(θ)RY(θ)RX(θ) works correctly
        let mut state1 = QulacsStateVector::new(1).unwrap();
        gates::rx(&mut state1, 0, PI / 4.0).unwrap();
        gates::ry(&mut state1, 0, PI / 4.0).unwrap();
        gates::rz(&mut state1, 0, PI / 4.0).unwrap();

        // Should produce a valid normalized state
        assert!((state1.norm_squared() - 1.0).abs() < 1e-10);

        // Compare with U3 which should be able to reproduce this
        let mut state2 = QulacsStateVector::new(1).unwrap();
        // U3(θ, φ, λ) with appropriate parameters
        gates::u3(&mut state2, 0, PI / 4.0, PI / 4.0, PI / 4.0).unwrap();

        assert!((state2.norm_squared() - 1.0).abs() < 1e-10);
    }

    // ========== Measurement Tests ==========

    #[test]
    fn test_probability_calculation() {
        // Test probability on |0⟩ state
        let state0 = QulacsStateVector::new(1).unwrap();
        assert!((state0.probability_zero(0).unwrap() - 1.0).abs() < 1e-10);
        assert!(state0.probability_one(0).unwrap().abs() < 1e-10);

        // Test probability on |1⟩ state
        let mut state1 = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state1, 0).unwrap();
        assert!(state1.probability_zero(0).unwrap().abs() < 1e-10);
        assert!((state1.probability_one(0).unwrap() - 1.0).abs() < 1e-10);

        // Test probability on |+⟩ state
        let mut state_plus = QulacsStateVector::new(1).unwrap();
        gates::hadamard(&mut state_plus, 0).unwrap();
        assert!((state_plus.probability_zero(0).unwrap() - 0.5).abs() < 1e-10);
        assert!((state_plus.probability_one(0).unwrap() - 0.5).abs() < 1e-10);
    }

    #[test]
    fn test_measurement_collapse() {
        // Measure |+⟩ state multiple times and check outcomes
        let mut outcomes_0 = 0;
        let mut outcomes_1 = 0;
        let num_trials = 1000;

        for _ in 0..num_trials {
            let mut state = QulacsStateVector::new(1).unwrap();
            gates::hadamard(&mut state, 0).unwrap();

            let outcome = state.measure(0).unwrap();
            if outcome {
                outcomes_1 += 1;
            } else {
                outcomes_0 += 1;
            }

            // After measurement, state should be pure |0⟩ or |1⟩
            assert!((state.norm_squared() - 1.0).abs() < 1e-10);
            if outcome {
                assert!((state.probability_one(0).unwrap() - 1.0).abs() < 1e-10);
            } else {
                assert!((state.probability_zero(0).unwrap() - 1.0).abs() < 1e-10);
            }
        }

        // Should get roughly 50-50 split (allow 3-sigma deviation)
        let ratio = outcomes_1 as f64 / num_trials as f64;
        assert!(ratio > 0.4 && ratio < 0.6, "Ratio: {}", ratio);
    }

    #[test]
    fn test_sampling() {
        // Create Bell state
        let mut bell_state = QulacsStateVector::new(2).unwrap();
        gates::hadamard(&mut bell_state, 0).unwrap();
        gates::cnot(&mut bell_state, 0, 1).unwrap();

        // Sample 1000 times
        let samples = bell_state.sample(1000).unwrap();
        assert_eq!(samples.len(), 1000);

        // Count outcomes
        let mut count_00 = 0;
        let mut count_11 = 0;
        for bitstring in &samples {
            assert_eq!(bitstring.len(), 2);
            if !bitstring[0] && !bitstring[1] {
                count_00 += 1;
            } else if bitstring[0] && bitstring[1] {
                count_11 += 1;
            } else {
                // Bell state should only produce |00⟩ or |11⟩
                panic!("Unexpected outcome: {:?}", bitstring);
            }
        }

        // Should get roughly equal counts (allow 3-sigma deviation)
        let ratio = count_00 as f64 / 1000.0;
        assert!(ratio > 0.4 && ratio < 0.6, "Ratio |00⟩: {}", ratio);

        // Verify sampling doesn't change state
        assert!((bell_state.norm_squared() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_get_counts() {
        // Create |+⟩ state
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::hadamard(&mut state, 0).unwrap();

        let counts = state.get_counts(1000).unwrap();

        // Should have two entries: [false] and [true]
        assert!(counts.len() <= 2);

        let count_0 = *counts.get(&vec![false]).unwrap_or(&0);
        let count_1 = *counts.get(&vec![true]).unwrap_or(&0);

        assert_eq!(count_0 + count_1, 1000);

        // Should be roughly 50-50
        let ratio = count_1 as f64 / 1000.0;
        assert!(ratio > 0.4 && ratio < 0.6, "Ratio |1⟩: {}", ratio);
    }

    #[test]
    fn test_sample_qubits() {
        // Create 2-qubit Bell state for now (TODO: fix 3-qubit CNOT)
        let mut bell = QulacsStateVector::new(2).unwrap();
        gates::hadamard(&mut bell, 0).unwrap();
        gates::cnot(&mut bell, 0, 1).unwrap();

        // Sample only qubit 0
        let samples = bell.sample_qubits(&[0], 1000).unwrap();

        let mut count_0 = 0;
        let mut count_1 = 0;

        for bitstring in &samples {
            assert_eq!(bitstring.len(), 1);
            if bitstring[0] {
                count_1 += 1;
            } else {
                count_0 += 1;
            }
        }

        // Should be roughly 50-50
        let ratio = count_1 as f64 / 1000.0;
        assert!(ratio > 0.4 && ratio < 0.6, "Ratio: {}", ratio);
    }

    #[test]
    fn test_measurement_multi_qubit() {
        // Create Bell state and measure both qubits
        let mut bell_state = QulacsStateVector::new(2).unwrap();
        gates::hadamard(&mut bell_state, 0).unwrap();
        gates::cnot(&mut bell_state, 0, 1).unwrap();

        let outcome0 = bell_state.measure(0).unwrap();
        let outcome1 = bell_state.measure(1).unwrap();

        // For Bell state, both qubits should be correlated
        assert_eq!(outcome0, outcome1);
    }

    #[test]
    fn test_toffoli_gate() {
        // Test Toffoli gate: CCX gate that flips target when both controls are |1⟩
        let mut state = QulacsStateVector::new(3).unwrap();

        // Case 1: Initial state |000⟩ - no change expected
        gates::toffoli(&mut state, 0, 1, 2).unwrap();
        assert!((state.amplitudes()[0].norm() - 1.0).abs() < 1e-10);
        assert!(state.amplitudes()[7].norm() < 1e-10);

        // Case 2: State |011⟩ - should flip to |111⟩
        let mut state2 = QulacsStateVector::new(3).unwrap();
        gates::pauli_x(&mut state2, 0).unwrap(); // |100⟩
        gates::pauli_x(&mut state2, 1).unwrap(); // |110⟩

        // Apply Toffoli(0, 1, 2)
        gates::toffoli(&mut state2, 0, 1, 2).unwrap();

        // Should now be |111⟩
        assert!((state2.amplitudes()[7].norm() - 1.0).abs() < 1e-10);
        assert!(state2.amplitudes()[3].norm() < 1e-10); // |110⟩ should be zero

        // Case 3: Test with superposition - create |+00⟩
        let mut state3 = QulacsStateVector::new(3).unwrap();
        gates::hadamard(&mut state3, 0).unwrap();

        // Apply X to all qubits
        gates::pauli_x(&mut state3, 0).unwrap();
        gates::pauli_x(&mut state3, 1).unwrap();
        gates::pauli_x(&mut state3, 2).unwrap();

        // Now have |-11⟩, apply Toffoli
        gates::toffoli(&mut state3, 0, 1, 2).unwrap();

        // Verify state is still normalized
        assert!((state3.norm_squared() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_toffoli_reversibility() {
        // Toffoli is self-inverse
        let mut state1 = QulacsStateVector::new(3).unwrap();
        gates::hadamard(&mut state1, 0).unwrap();
        gates::hadamard(&mut state1, 1).unwrap();
        gates::hadamard(&mut state1, 2).unwrap();

        let original_state = state1.clone();

        // Apply Toffoli twice
        gates::toffoli(&mut state1, 0, 1, 2).unwrap();
        gates::toffoli(&mut state1, 0, 1, 2).unwrap();

        // Should return to original state
        for i in 0..8 {
            let diff = (state1.amplitudes()[i] - original_state.amplitudes()[i]).norm();
            assert!(diff < 1e-10, "Difference at index {}: {}", i, diff);
        }
    }

    #[test]
    fn test_fredkin_gate() {
        // Test Fredkin gate: CSWAP gate that swaps target1 and target2 when control is |1⟩
        let mut state = QulacsStateVector::new(3).unwrap();

        // Case 1: Initial state |000⟩ - no change expected (control=0)
        gates::fredkin(&mut state, 0, 1, 2).unwrap();
        assert!((state.amplitudes()[0].norm() - 1.0).abs() < 1e-10);

        // Case 2: State |101⟩ - should swap to |011⟩
        // Ket notation |q2 q1 q0⟩:
        //   Before: |1 0 1⟩ → index = 1 + 0 + 4 = 5 = 0b101
        //   After Fredkin(control=q0, target1=q1, target2=q2): swap q1 and q2
        //   After:  |0 1 1⟩ → index = 1 + 2 + 0 = 3 = 0b011
        let mut state2 = QulacsStateVector::new(3).unwrap();
        gates::pauli_x(&mut state2, 0).unwrap(); // qubit 0 = 1 (control)
        gates::pauli_x(&mut state2, 2).unwrap(); // qubit 2 = 1

        // State is |101⟩ (ket notation: q2=1, q1=0, q0=1)
        assert!((state2.amplitudes()[0b101].norm() - 1.0).abs() < 1e-10);

        // Apply Fredkin(0, 1, 2) - swap qubits 1 and 2 since control(0)=1
        gates::fredkin(&mut state2, 0, 1, 2).unwrap();

        // Should now be at index 0b011 (ket |011⟩: q2=0, q1=1, q0=1)
        assert!((state2.amplitudes()[0b011].norm() - 1.0).abs() < 1e-10);
        assert!(state2.amplitudes()[0b101].norm() < 1e-10);

        // Case 3: Control qubit is |0⟩ - no swap
        let mut state3 = QulacsStateVector::new(3).unwrap();
        gates::pauli_x(&mut state3, 1).unwrap(); // target1=1
        gates::pauli_x(&mut state3, 2).unwrap(); // target2=1

        // State is |011⟩, control(0)=0
        let before = state3.clone();
        gates::fredkin(&mut state3, 0, 1, 2).unwrap();

        // Should be unchanged
        for i in 0..8 {
            let diff = (state3.amplitudes()[i] - before.amplitudes()[i]).norm();
            assert!(diff < 1e-10);
        }
    }

    #[test]
    fn test_fredkin_reversibility() {
        // Fredkin is self-inverse
        let mut state1 = QulacsStateVector::new(3).unwrap();
        gates::hadamard(&mut state1, 0).unwrap();
        gates::hadamard(&mut state1, 1).unwrap();
        gates::hadamard(&mut state1, 2).unwrap();

        let original_state = state1.clone();

        // Apply Fredkin twice
        gates::fredkin(&mut state1, 0, 1, 2).unwrap();
        gates::fredkin(&mut state1, 0, 1, 2).unwrap();

        // Should return to original state
        for i in 0..8 {
            let diff = (state1.amplitudes()[i] - original_state.amplitudes()[i]).norm();
            assert!(diff < 1e-10, "Difference at index {}: {}", i, diff);
        }
    }

    #[test]
    fn test_toffoli_error_cases() {
        let mut state = QulacsStateVector::new(3).unwrap();

        // Test invalid qubit indices
        assert!(gates::toffoli(&mut state, 0, 1, 5).is_err());
        assert!(gates::toffoli(&mut state, 5, 1, 2).is_err());

        // Test same qubits
        assert!(gates::toffoli(&mut state, 0, 0, 2).is_err());
        assert!(gates::toffoli(&mut state, 0, 1, 1).is_err());
        assert!(gates::toffoli(&mut state, 0, 1, 0).is_err());
    }

    #[test]
    fn test_fredkin_error_cases() {
        let mut state = QulacsStateVector::new(3).unwrap();

        // Test invalid qubit indices
        assert!(gates::fredkin(&mut state, 5, 1, 2).is_err());
        assert!(gates::fredkin(&mut state, 0, 5, 2).is_err());
        assert!(gates::fredkin(&mut state, 0, 1, 5).is_err());

        // Test same qubits
        assert!(gates::fredkin(&mut state, 0, 0, 2).is_err());
        assert!(gates::fredkin(&mut state, 0, 1, 1).is_err());
        assert!(gates::fredkin(&mut state, 0, 1, 0).is_err());
    }

    #[test]
    fn test_s_gate() {
        // S gate: |0⟩ → |0⟩, |1⟩ → i|1⟩
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state, 0).unwrap(); // |1⟩
        gates::s(&mut state, 0).unwrap();

        // State should be i|1⟩
        assert!((state.amplitudes()[0].norm() - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[1].re - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[1].im - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_s_dag_gate() {
        // S† gate: |0⟩ → |0⟩, |1⟩ → -i|1⟩
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state, 0).unwrap(); // |1⟩
        gates::sdg(&mut state, 0).unwrap();

        // State should be -i|1⟩
        assert!((state.amplitudes()[0].norm() - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[1].re - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[1].im + 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_s_s_dag_identity() {
        // S · S† = I
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::hadamard(&mut state, 0).unwrap(); // |+⟩

        let original = state.clone();

        gates::s(&mut state, 0).unwrap();
        gates::sdg(&mut state, 0).unwrap();

        // Should return to original state
        for i in 0..2 {
            let diff = (state.amplitudes()[i] - original.amplitudes()[i]).norm();
            assert!(diff < 1e-10);
        }
    }

    #[test]
    fn test_t_gate() {
        // T gate: |0⟩ → |0⟩, |1⟩ → e^(iπ/4)|1⟩
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state, 0).unwrap(); // |1⟩
        gates::t(&mut state, 0).unwrap();

        // State should be e^(iπ/4)|1⟩
        let expected = Complex64::from_polar(1.0, std::f64::consts::FRAC_PI_4);
        assert!((state.amplitudes()[0].norm() - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[1].re - expected.re).abs() < 1e-10);
        assert!((state.amplitudes()[1].im - expected.im).abs() < 1e-10);
    }

    #[test]
    fn test_t_dag_gate() {
        // T† gate: |0⟩ → |0⟩, |1⟩ → e^(-iπ/4)|1⟩
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state, 0).unwrap(); // |1⟩
        gates::tdg(&mut state, 0).unwrap();

        // State should be e^(-iπ/4)|1⟩
        let expected = Complex64::from_polar(1.0, -std::f64::consts::FRAC_PI_4);
        assert!((state.amplitudes()[0].norm() - 0.0).abs() < 1e-10);
        assert!((state.amplitudes()[1].re - expected.re).abs() < 1e-10);
        assert!((state.amplitudes()[1].im - expected.im).abs() < 1e-10);
    }

    #[test]
    fn test_t_t_dag_identity() {
        // T · T† = I
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::hadamard(&mut state, 0).unwrap(); // |+⟩

        let original = state.clone();

        gates::t(&mut state, 0).unwrap();
        gates::tdg(&mut state, 0).unwrap();

        // Should return to original state
        for i in 0..2 {
            let diff = (state.amplitudes()[i] - original.amplitudes()[i]).norm();
            assert!(diff < 1e-10);
        }
    }

    #[test]
    fn test_s_equals_two_t() {
        // S = T · T (since π/4 + π/4 = π/2)
        let mut state1 = QulacsStateVector::new(1).unwrap();
        let mut state2 = QulacsStateVector::new(1).unwrap();

        gates::hadamard(&mut state1, 0).unwrap();
        gates::hadamard(&mut state2, 0).unwrap();

        // Apply S to state1
        gates::s(&mut state1, 0).unwrap();

        // Apply T twice to state2
        gates::t(&mut state2, 0).unwrap();
        gates::t(&mut state2, 0).unwrap();

        // Should be the same
        for i in 0..2 {
            let diff = (state1.amplitudes()[i] - state2.amplitudes()[i]).norm();
            assert!(diff < 1e-10);
        }
    }

    // Observable framework tests
    #[test]
    fn test_pauli_operator_matrices() {
        use observable::PauliOperator;

        // Test identity
        let i_mat = PauliOperator::I.matrix();
        assert_eq!(i_mat[[0, 0]], Complex64::new(1.0, 0.0));
        assert_eq!(i_mat[[1, 1]], Complex64::new(1.0, 0.0));

        // Test Pauli X
        let x_mat = PauliOperator::X.matrix();
        assert_eq!(x_mat[[0, 1]], Complex64::new(1.0, 0.0));
        assert_eq!(x_mat[[1, 0]], Complex64::new(1.0, 0.0));

        // Test Pauli Z
        let z_mat = PauliOperator::Z.matrix();
        assert_eq!(z_mat[[0, 0]], Complex64::new(1.0, 0.0));
        assert_eq!(z_mat[[1, 1]], Complex64::new(-1.0, 0.0));
    }

    #[test]
    fn test_pauli_observable_creation() {
        use observable::PauliObservable;

        // Create Pauli Z observable
        let obs_z = PauliObservable::pauli_z(&[0]);
        assert_eq!(obs_z.coefficient, 1.0);
        assert_eq!(obs_z.operators.len(), 1);

        // Create Pauli X observable
        let obs_x = PauliObservable::pauli_x(&[0, 1]);
        assert_eq!(obs_x.operators.len(), 2);
    }

    #[test]
    fn test_pauli_z_expectation_value() {
        use observable::PauliObservable;

        // Create |0⟩ state
        let state = QulacsStateVector::new(1).unwrap();

        // Measure Pauli Z - should give +1 for |0⟩
        let obs = PauliObservable::pauli_z(&[0]);
        let exp_val = obs.expectation_value(&state);
        assert!((exp_val - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_pauli_z_expectation_value_excited() {
        use observable::PauliObservable;

        // Create |1⟩ state
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::pauli_x(&mut state, 0).unwrap();

        // Measure Pauli Z - should give -1 for |1⟩
        let obs = PauliObservable::pauli_z(&[0]);
        let exp_val = obs.expectation_value(&state);
        assert!((exp_val - (-1.0)).abs() < 1e-10);
    }

    #[test]
    fn test_pauli_z_expectation_value_superposition() {
        use observable::PauliObservable;

        // Create |+⟩ state (equal superposition)
        let mut state = QulacsStateVector::new(1).unwrap();
        gates::hadamard(&mut state, 0).unwrap();

        // Measure Pauli Z - should give 0 for |+⟩
        let obs = PauliObservable::pauli_z(&[0]);
        let exp_val = obs.expectation_value(&state);
        assert!(exp_val.abs() < 1e-10);
    }

    #[test]
    fn test_hermitian_observable() {
        use observable::HermitianObservable;

        // Create Pauli Z matrix manually
        let mut matrix = Array2::zeros((2, 2));
        matrix[[0, 0]] = Complex64::new(1.0, 0.0);
        matrix[[1, 1]] = Complex64::new(-1.0, 0.0);

        let obs = HermitianObservable::new(matrix).unwrap();
        assert_eq!(obs.num_qubits, 1);

        // Test on |0⟩ state
        let state = QulacsStateVector::new(1).unwrap();
        let exp_val = obs.expectation_value(&state).unwrap();
        assert!((exp_val - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_composite_observable() {
        use observable::{CompositeObservable, PauliObservable};

        // Create composite observable: 0.5 * Z_0 + 0.3 * Z_1
        let obs1 = PauliObservable::pauli_z(&[0]).with_coefficient(0.5);
        let obs2 = PauliObservable::pauli_z(&[1]).with_coefficient(0.3);

        let composite = CompositeObservable::new().add_term(obs1).add_term(obs2);

        assert_eq!(composite.num_terms(), 2);

        // Test on |00⟩ state - should give 0.5 * 1.0 + 0.3 * 1.0 = 0.8
        let state = QulacsStateVector::new(2).unwrap();
        let exp_val = composite.expectation_value(&state);
        assert!((exp_val - 0.8).abs() < 1e-10);
    }

    #[test]
    fn test_observable_weight() {
        use observable::{PauliObservable, PauliOperator};
        use std::collections::HashMap;

        let mut operators = HashMap::new();
        operators.insert(0, PauliOperator::X);
        operators.insert(1, PauliOperator::Y);
        operators.insert(2, PauliOperator::I);

        let obs = PauliObservable::new(operators, 1.0);
        assert_eq!(obs.weight(), 2); // X and Y are non-identity
    }

    // ===== Circuit API Tests =====

    #[test]
    fn test_circuit_api_basic() {
        use circuit_api::QulacsCircuit;

        let mut circuit = QulacsCircuit::new(2).unwrap();
        assert_eq!(circuit.num_qubits(), 2);
        assert_eq!(circuit.gate_count(), 0);

        circuit.h(0).x(1);
        assert_eq!(circuit.gate_count(), 2);
    }

    #[test]
    fn test_circuit_api_bell_state() {
        use circuit_api::QulacsCircuit;

        let mut circuit = QulacsCircuit::new(2).unwrap();
        circuit.bell_pair(0, 1);

        assert_eq!(circuit.gate_count(), 2);

        let probs = circuit.probabilities();
        assert!((probs[0] - 0.5).abs() < 1e-10); // |00⟩
        assert!(probs[1].abs() < 1e-10); // |01⟩
        assert!(probs[2].abs() < 1e-10); // |10⟩
        assert!((probs[3] - 0.5).abs() < 1e-10); // |11⟩
    }

    #[test]
    fn test_circuit_api_run_shots() {
        use circuit_api::QulacsCircuit;

        let mut circuit = QulacsCircuit::new(2).unwrap();
        circuit.bell_pair(0, 1);

        let counts = circuit.run(100).unwrap();

        // Should have roughly 50% |00⟩ and 50% |11⟩
        assert!(counts.contains_key("00") || counts.contains_key("11"));
        let total: usize = counts.values().sum();
        assert_eq!(total, 100);
    }

    #[test]
    fn test_circuit_api_reset() {
        use circuit_api::QulacsCircuit;

        let mut circuit = QulacsCircuit::new(2).unwrap();
        circuit.h(0).cnot(0, 1);
        assert_eq!(circuit.gate_count(), 2);

        circuit.reset().unwrap();
        assert_eq!(circuit.gate_count(), 0);

        // State should be |00⟩
        let probs = circuit.probabilities();
        assert!((probs[0] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_circuit_api_rotation_gates() {
        use circuit_api::QulacsCircuit;
        use std::f64::consts::PI;

        let mut circuit = QulacsCircuit::new(1).unwrap();

        // RX(π) should be equivalent to X (up to global phase)
        circuit.rx(0, PI);
        let probs = circuit.probabilities();
        assert!(probs[0].abs() < 1e-10);
        assert!((probs[1] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_circuit_api_phase_gates() {
        use circuit_api::QulacsCircuit;

        let mut circuit = QulacsCircuit::new(1).unwrap();

        // Apply S twice should equal Z
        circuit.s(0).s(0);
        assert_eq!(circuit.gate_count(), 2);

        // |0⟩ should remain |0⟩
        let probs = circuit.probabilities();
        assert!((probs[0] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_circuit_api_two_qubit_gates() {
        use circuit_api::QulacsCircuit;

        let mut circuit = QulacsCircuit::new(2).unwrap();

        // CNOT with control in |1⟩ should flip target
        circuit.x(0).cnot(0, 1);

        let probs = circuit.probabilities();
        assert!(probs[0].abs() < 1e-10); // |00⟩
        assert!(probs[1].abs() < 1e-10); // |01⟩
        assert!(probs[2].abs() < 1e-10); // |10⟩
        assert!((probs[3] - 1.0).abs() < 1e-10); // |11⟩
    }

    #[test]
    fn test_circuit_api_observable() {
        use circuit_api::{Observable, QulacsCircuit};
        use observable::PauliObservable;

        let mut circuit = QulacsCircuit::new(2).unwrap();
        circuit.h(0).h(1);

        let obs = PauliObservable::pauli_z(&[0]);
        let exp_val = circuit.expectation(&obs).unwrap();

        // Hadamard on |0⟩ gives equal superposition, so Z expectation is 0
        assert!(exp_val.abs() < 1e-10);
    }

    #[test]
    fn test_circuit_api_gate_record() {
        use circuit_api::QulacsCircuit;

        let mut circuit = QulacsCircuit::new(2).unwrap();
        circuit.h(0).cnot(0, 1).rx(1, 1.5);

        let gates = circuit.gates();
        assert_eq!(gates.len(), 3);
        assert_eq!(gates[0].name, "H");
        assert_eq!(gates[0].qubits, vec![0]);
        assert_eq!(gates[1].name, "CNOT");
        assert_eq!(gates[1].qubits, vec![0, 1]);
        assert_eq!(gates[2].name, "RX");
        assert_eq!(gates[2].params[0], 1.5);
    }

    #[test]
    fn test_circuit_api_noise_model() {
        use crate::noise_models::{DepolarizingNoise, NoiseModel as KrausNoiseModel};
        use circuit_api::QulacsCircuit;
        use std::sync::Arc;

        let mut circuit = QulacsCircuit::new(2).unwrap();
        assert!(!circuit.has_noise_model());

        // Set a noise model
        let mut noise_model = KrausNoiseModel::new();
        noise_model.add_channel(Arc::new(DepolarizingNoise::new(0.01)));

        circuit.set_noise_model(noise_model);
        assert!(circuit.has_noise_model());

        // Apply gates
        circuit.h(0).cnot(0, 1);
        assert_eq!(circuit.gate_count(), 2);

        // Clear noise model
        circuit.clear_noise_model();
        assert!(!circuit.has_noise_model());
    }

    #[test]
    fn test_circuit_api_run_with_noise() {
        use crate::noise_models::{BitFlipNoise, NoiseModel as KrausNoiseModel};
        use circuit_api::QulacsCircuit;
        use std::sync::Arc;

        let mut circuit = QulacsCircuit::new(1).unwrap();

        // Add bit flip noise
        let mut noise_model = KrausNoiseModel::new();
        noise_model.add_channel(Arc::new(BitFlipNoise::new(0.1)));
        circuit.set_noise_model(noise_model);

        // Start in |0⟩, with bit flip noise we should occasionally see |1⟩
        let counts = circuit.run_with_noise(100).unwrap();

        // Verify we get measurements
        let total: usize = counts.values().sum();
        assert_eq!(total, 100);

        // With 10% bit flip noise, we should mostly see |0⟩ but occasionally |1⟩
        // (This is a probabilistic test, but with 100 shots the statistics should be reasonable)
        assert!(counts.contains_key("0"));
    }
}