ppflib 0.1.0

//! PPF Hilbert space construction
//!
//! This module implements the PPF Hilbert space over factorization states
//! with complexity-based inner product ⟨f|g⟩ = δ_fg e^(-Ω(f)/2) and
//! quantum operators: Ŝ (sign), Ω̂ (complexity), K̂ (curvature-flip).
//! The space carries an SU(2) algebra structure.

use crate::core::{FactorizationStateSpace, PFactorization};
use crate::geometry::iot::IOTMetric;
use std::collections::HashMap;
use std::fmt;
use thiserror::Error;
use serde::{Deserialize, Serialize};
use num_complex::Complex64;

/// Errors for Hilbert space operations
#[derive(Error, Debug, Clone, PartialEq)]
pub enum HilbertSpaceError {
    /// Invalid state
    #[error("Invalid state: {0}")]
    InvalidState(String),
    /// Operator error
    #[error("Operator error: {0}")]
    OperatorError(String),
    /// Dimension mismatch
    #[error("Dimension mismatch: expected {expected}, got {actual}")]
    DimensionMismatch { 
        /// Expected dimension
        expected: usize, 
        /// Actual dimension
        actual: usize 
    },
    /// Not normalized
    #[error("State not normalized: norm = {0}")]
    NotNormalized(f64),
}

/// A quantum state in the PPF Hilbert space
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct QuantumState {
    /// Amplitudes for each factorization basis state
    amplitudes: HashMap<String, Complex64>,
    /// Dimension of the Hilbert space
    dimension: usize,
    /// Whether the state is normalized
    is_normalized: bool,
}

impl QuantumState {
    /// Create a new quantum state
    pub fn new(amplitudes: HashMap<String, Complex64>) -> Result<Self, HilbertSpaceError> {
        let dimension = amplitudes.len();
        
        if dimension == 0 {
            return Err(HilbertSpaceError::InvalidState(
                "Cannot create state with no amplitudes".to_string()
            ));
        }
        
        let mut state = QuantumState {
            amplitudes,
            dimension,
            is_normalized: false,
        };
        
        // Check if it's a zero state (all amplitudes are zero)
        let is_zero_state = state.amplitudes.values().all(|&amp| amp.norm() < 1e-15);
        
        if is_zero_state {
            // Zero states don't need normalization
            state.is_normalized = true;
            Ok(state)
        } else {
            state.normalize()?;
            Ok(state)
        }
    }

    /// Create a basis state (single factorization)
    pub fn basis_state(factorization: &PFactorization, hilbert_space_dimension: usize) -> Self {
        let mut amplitudes = HashMap::new();
        amplitudes.insert(factorization.to_string(), Complex64::new(1.0, 0.0));
        
        QuantumState {
            amplitudes,
            dimension: hilbert_space_dimension,
            is_normalized: true,
        }
    }

    /// Create a uniform superposition state
    pub fn uniform_superposition(state_space: &FactorizationStateSpace) -> Result<Self, HilbertSpaceError> {
        let n = state_space.size();
        let amplitude = Complex64::new(1.0 / (n as f64).sqrt(), 0.0);
        
        let mut amplitudes = HashMap::new();
        for factorization in state_space.factorizations() {
            amplitudes.insert(factorization.to_string(), amplitude);
        }
        
        Ok(QuantumState {
            amplitudes,
            dimension: n,
            is_normalized: true,
        })
    }

    /// Get the amplitude for a factorization
    pub fn amplitude(&self, factorization: &PFactorization) -> Complex64 {
        self.amplitudes.get(&factorization.to_string())
            .copied()
            .unwrap_or(Complex64::new(0.0, 0.0))
    }

    /// Get all amplitudes
    pub fn amplitudes(&self) -> &HashMap<String, Complex64> {
        &self.amplitudes
    }

    /// Get the dimension
    pub fn dimension(&self) -> usize {
        self.dimension
    }

    /// Compute the norm
    pub fn norm(&self) -> f64 {
        let norm_squared: f64 = self.amplitudes.values()
            .map(|z| z.norm_sqr())
            .sum();
        norm_squared.sqrt()
    }

    /// Normalize the state
    pub fn normalize(&mut self) -> Result<(), HilbertSpaceError> {
        let norm = self.norm();
        
        if norm < 1e-15 {
            return Err(HilbertSpaceError::InvalidState(
                "Cannot normalize zero state".to_string()
            ));
        }
        
        for amplitude in self.amplitudes.values_mut() {
            *amplitude /= norm;
        }
        
        self.is_normalized = true;
        Ok(())
    }

    /// Check if the state is normalized
    pub fn is_normalized(&self) -> bool {
        self.is_normalized && (self.norm() - 1.0).abs() < 1e-10
    }

    /// Add two quantum states
    pub fn add(&self, other: &QuantumState) -> Result<QuantumState, HilbertSpaceError> {
        if self.dimension != other.dimension {
            return Err(HilbertSpaceError::DimensionMismatch {
                expected: self.dimension,
                actual: other.dimension,
            });
        }

        let mut result_amplitudes = HashMap::new();
        
        // Add amplitudes from both states
        for (key, &amplitude) in &self.amplitudes {
            result_amplitudes.insert(key.clone(), amplitude);
        }
        
        for (key, &amplitude) in &other.amplitudes {
            let existing = result_amplitudes.get(key).copied().unwrap_or(Complex64::new(0.0, 0.0));
            result_amplitudes.insert(key.clone(), existing + amplitude);
        }
        
        // Remove zero amplitudes
        result_amplitudes.retain(|_, &mut amp| amp.norm() > 1e-15);
        
        QuantumState::new(result_amplitudes)
    }

    /// Scalar multiplication
    pub fn scalar_mult(&self, scalar: Complex64) -> QuantumState {
        let mut result_amplitudes = HashMap::new();
        
        for (key, &amplitude) in &self.amplitudes {
            result_amplitudes.insert(key.clone(), scalar * amplitude);
        }
        
        QuantumState {
            amplitudes: result_amplitudes,
            dimension: self.dimension,
            is_normalized: false,
        }
    }
}

/// The PPF Hilbert space
#[derive(Debug, Clone)]
pub struct PPFHilbertSpace {
    /// The factorization state space
    state_space: FactorizationStateSpace,
    /// Basis states (one per factorization)
    basis_states: Vec<QuantumState>,
    /// Inner product metric
    metric: InnerProductMetric,
}

impl PPFHilbertSpace {
    /// Create a new PPF Hilbert space
    pub fn new(state_space: FactorizationStateSpace) -> Self {
        let basis_states: Vec<QuantumState> = state_space.factorizations()
            .iter()
            .map(|f| QuantumState::basis_state(f, state_space.size()))
            .collect();
        
        let metric = InnerProductMetric::new(&state_space);
        
        PPFHilbertSpace {
            state_space,
            basis_states,
            metric,
        }
    }

    /// Get the dimension of the Hilbert space
    pub fn dimension(&self) -> usize {
        self.state_space.size()
    }

    /// Get the basis states
    pub fn basis_states(&self) -> &[QuantumState] {
        &self.basis_states
    }

    /// Compute the PPF inner product ⟨f|g⟩ = δ_fg e^(-Ω(f)/2)
    pub fn inner_product(&self, state1: &QuantumState, state2: &QuantumState) -> Complex64 {
        let mut result = Complex64::new(0.0, 0.0);
        
        for (factorization_str, &amp1) in &state1.amplitudes {
            if let Some(&amp2) = state2.amplitudes.get(factorization_str) {
                // Find the factorization to get its complexity
                let factorization = self.state_space.factorizations()
                    .iter()
                    .find(|f| f.to_string() == *factorization_str);
                
                if let Some(f) = factorization {
                    let complexity = f.complexity() as f64;
                    let weight = (-complexity / 2.0).exp();
                    result += amp1.conj() * amp2 * weight;
                }
            }
        }
        
        result
    }

    /// Compute the norm in the PPF metric
    pub fn norm(&self, state: &QuantumState) -> f64 {
        self.inner_product(state, state).norm().sqrt()
    }

    /// Create the sign operator Ŝ
    pub fn sign_operator(&self) -> QuantumOperator {
        let mut matrix = HashMap::new();
        
        for (i, factorization_i) in self.state_space.factorizations().iter().enumerate() {
            for (j, factorization_j) in self.state_space.factorizations().iter().enumerate() {
                // Sign operator flips the overall sign
                let flipped_value = -factorization_i.value();
                
                // Check if factorization_j corresponds to the sign-flipped version
                if factorization_j.value() == flipped_value {
                    matrix.insert((i, j), Complex64::new(1.0, 0.0));
                }
            }
        }
        
        QuantumOperator::new(matrix, self.dimension())
    }

    /// Create the complexity operator Ω̂
    pub fn complexity_operator(&self) -> QuantumOperator {
        let mut matrix = HashMap::new();
        
        for (i, factorization) in self.state_space.factorizations().iter().enumerate() {
            let complexity = factorization.complexity() as f64;
            matrix.insert((i, i), Complex64::new(complexity, 0.0));
        }
        
        QuantumOperator::new(matrix, self.dimension())
    }

    /// Create the curvature-flip operator K̂
    pub fn curvature_flip_operator(&self, iot_metric: &IOTMetric) -> QuantumOperator {
        let mut matrix = HashMap::new();
        
        for (i, factorization_i) in self.state_space.factorizations().iter().enumerate() {
            for (j, factorization_j) in self.state_space.factorizations().iter().enumerate() {
                // Curvature-flip connects factorizations based on geometric properties
                let coords_i = iot_metric.factorization_to_coordinates(factorization_i);
                let coords_j = iot_metric.factorization_to_coordinates(factorization_j);
                
                let curvature_i = iot_metric.ricci_scalar(&coords_i);
                let curvature_j = iot_metric.ricci_scalar(&coords_j);
                
                // Matrix element based on curvature correlation
                let coupling = (curvature_i * curvature_j).abs().sqrt() * 0.1;
                if coupling > 1e-10 {
                    matrix.insert((i, j), Complex64::new(coupling, 0.0));
                }
            }
        }
        
        QuantumOperator::new(matrix, self.dimension())
    }

    /// Create SU(2) generators
    pub fn su2_generators(&self) -> (QuantumOperator, QuantumOperator, QuantumOperator) {
        let _dim = self.dimension();
        
        // Pauli matrices scaled for the Hilbert space dimension
        let sigma_x = self.pauli_x_scaled();
        let sigma_y = self.pauli_y_scaled();
        let sigma_z = self.pauli_z_scaled();
        
        (sigma_x, sigma_y, sigma_z)
    }

    /// Create scaled Pauli-X operator
    fn pauli_x_scaled(&self) -> QuantumOperator {
        let mut matrix = HashMap::new();
        let dim = self.dimension();
        
        // Create a block-diagonal structure for higher dimensions
        for i in 0..dim {
            for j in 0..dim {
                if (i + j) % 2 == 1 {
                    let coupling = 1.0 / (dim as f64).sqrt();
                    matrix.insert((i, j), Complex64::new(coupling, 0.0));
                }
            }
        }
        
        QuantumOperator::new(matrix, dim)
    }

    /// Create scaled Pauli-Y operator
    fn pauli_y_scaled(&self) -> QuantumOperator {
        let mut matrix = HashMap::new();
        let dim = self.dimension();
        
        for i in 0..dim {
            for j in 0..dim {
                if (i + j) % 2 == 1 {
                    let coupling = 1.0 / (dim as f64).sqrt();
                    let sign = if i > j { 1.0 } else { -1.0 };
                    matrix.insert((i, j), Complex64::new(0.0, sign * coupling));
                }
            }
        }
        
        QuantumOperator::new(matrix, dim)
    }

    /// Create scaled Pauli-Z operator
    fn pauli_z_scaled(&self) -> QuantumOperator {
        let mut matrix = HashMap::new();
        let dim = self.dimension();
        
        for i in 0..dim {
            let eigenvalue = if i % 2 == 0 { 1.0 } else { -1.0 };
            matrix.insert((i, i), Complex64::new(eigenvalue, 0.0));
        }
        
        QuantumOperator::new(matrix, dim)
    }

    /// Verify SU(2) algebra relations
    pub fn verify_su2_algebra(&self) -> bool {
        let (sigma_x, sigma_y, sigma_z) = self.su2_generators();
        
        // Check [σ_x, σ_y] = 2iσ_z
        let commutator_xy = sigma_x.commutator(&sigma_y);
        let expected_z = sigma_z.scalar_mult(Complex64::new(0.0, 2.0));
        
        commutator_xy.approximately_equals(&expected_z, 1e-10)
    }

    /// Get the state space
    pub fn state_space(&self) -> &FactorizationStateSpace {
        &self.state_space
    }

    /// Get the metric
    pub fn metric(&self) -> &InnerProductMetric {
        &self.metric
    }
}

/// Inner product metric for PPF Hilbert space
#[derive(Debug, Clone)]
pub struct InnerProductMetric {
    /// Complexity weights for each factorization
    complexity_weights: HashMap<String, f64>,
}

impl InnerProductMetric {
    /// Create a new inner product metric
    pub fn new(state_space: &FactorizationStateSpace) -> Self {
        let mut complexity_weights = HashMap::new();
        
        for factorization in state_space.factorizations() {
            let complexity = factorization.complexity() as f64;
            let weight = (-complexity / 2.0).exp();
            complexity_weights.insert(factorization.to_string(), weight);
        }
        
        InnerProductMetric { complexity_weights }
    }

    /// Get the weight for a factorization
    pub fn weight(&self, factorization: &PFactorization) -> f64 {
        self.complexity_weights.get(&factorization.to_string())
            .copied()
            .unwrap_or(1.0)
    }
}

/// A quantum operator in the PPF Hilbert space
#[derive(Debug, Clone)]
pub struct QuantumOperator {
    /// Matrix elements (sparse representation)
    matrix: HashMap<(usize, usize), Complex64>,
    /// Dimension of the operator
    dimension: usize,
}

impl QuantumOperator {
    /// Create a new quantum operator
    pub fn new(matrix: HashMap<(usize, usize), Complex64>, dimension: usize) -> Self {
        QuantumOperator { matrix, dimension }
    }

    /// Get the dimension
    pub fn dimension(&self) -> usize {
        self.dimension
    }

    /// Get matrix element
    pub fn matrix_element(&self, i: usize, j: usize) -> Complex64 {
        self.matrix.get(&(i, j)).copied().unwrap_or(Complex64::new(0.0, 0.0))
    }

    /// Apply operator to a state
    pub fn apply(&self, state: &QuantumState) -> Result<QuantumState, HilbertSpaceError> {
        if state.dimension() != self.dimension {
            return Err(HilbertSpaceError::DimensionMismatch {
                expected: self.dimension,
                actual: state.dimension(),
            });
        }

        let mut result_amplitudes = HashMap::new();
        
        // Convert state amplitudes to indexed form
        let mut state_vector = vec![Complex64::new(0.0, 0.0); self.dimension];
        let mut factorization_keys = Vec::new();
        
        for (i, (key, &amplitude)) in state.amplitudes().iter().enumerate() {
            if i < self.dimension {
                state_vector[i] = amplitude;
                factorization_keys.push(key.clone());
            }
        }
        
        // Apply matrix multiplication
        for i in 0..self.dimension {
            let mut result_amplitude = Complex64::new(0.0, 0.0);
            
            for j in 0..self.dimension {
                result_amplitude += self.matrix_element(i, j) * state_vector[j];
            }
            
            if result_amplitude.norm() > 1e-15 {
                let key = if i < factorization_keys.len() {
                    factorization_keys[i].clone()
                } else {
                    format!("basis_{}", i)
                };
                result_amplitudes.insert(key, result_amplitude);
            }
        }
        
        if result_amplitudes.is_empty() {
            // Create a zero state - add a tiny amplitude to avoid empty state error
            result_amplitudes.insert("zero_state".to_string(), Complex64::new(0.0, 0.0));
        }
        QuantumState::new(result_amplitudes)
    }

    /// Compute the commutator [A, B] = AB - BA
    pub fn commutator(&self, other: &QuantumOperator) -> QuantumOperator {
        let ab = self.multiply(other);
        let ba = other.multiply(self);
        ab.subtract(&ba)
    }

    /// Compute the anticommutator {A, B} = AB + BA
    pub fn anticommutator(&self, other: &QuantumOperator) -> QuantumOperator {
        let ab = self.multiply(other);
        let ba = other.multiply(self);
        ab.add(&ba)
    }

    /// Multiply two operators
    pub fn multiply(&self, other: &QuantumOperator) -> QuantumOperator {
        let mut result_matrix = HashMap::new();
        
        for i in 0..self.dimension {
            for j in 0..self.dimension {
                let mut element = Complex64::new(0.0, 0.0);
                
                for k in 0..self.dimension {
                    element += self.matrix_element(i, k) * other.matrix_element(k, j);
                }
                
                if element.norm() > 1e-15 {
                    result_matrix.insert((i, j), element);
                }
            }
        }
        
        QuantumOperator::new(result_matrix, self.dimension)
    }

    /// Add two operators
    pub fn add(&self, other: &QuantumOperator) -> QuantumOperator {
        let mut result_matrix = HashMap::new();
        
        for i in 0..self.dimension {
            for j in 0..self.dimension {
                let element = self.matrix_element(i, j) + other.matrix_element(i, j);
                if element.norm() > 1e-15 {
                    result_matrix.insert((i, j), element);
                }
            }
        }
        
        QuantumOperator::new(result_matrix, self.dimension)
    }

    /// Subtract two operators
    pub fn subtract(&self, other: &QuantumOperator) -> QuantumOperator {
        let mut result_matrix = HashMap::new();
        
        for i in 0..self.dimension {
            for j in 0..self.dimension {
                let element = self.matrix_element(i, j) - other.matrix_element(i, j);
                if element.norm() > 1e-15 {
                    result_matrix.insert((i, j), element);
                }
            }
        }
        
        QuantumOperator::new(result_matrix, self.dimension)
    }

    /// Scalar multiplication
    pub fn scalar_mult(&self, scalar: Complex64) -> QuantumOperator {
        let mut result_matrix = HashMap::new();
        
        for (&(i, j), &element) in &self.matrix {
            result_matrix.insert((i, j), scalar * element);
        }
        
        QuantumOperator::new(result_matrix, self.dimension)
    }

    /// Check if two operators are approximately equal
    pub fn approximately_equals(&self, other: &QuantumOperator, tolerance: f64) -> bool {
        if self.dimension != other.dimension {
            return false;
        }
        
        for i in 0..self.dimension {
            for j in 0..self.dimension {
                let diff = self.matrix_element(i, j) - other.matrix_element(i, j);
                if diff.norm() > tolerance {
                    return false;
                }
            }
        }
        
        true
    }

    /// Compute the trace
    pub fn trace(&self) -> Complex64 {
        let mut result = Complex64::new(0.0, 0.0);
        
        for i in 0..self.dimension {
            result += self.matrix_element(i, i);
        }
        
        result
    }

    /// Check if the operator is Hermitian
    pub fn is_hermitian(&self) -> bool {
        for i in 0..self.dimension {
            for j in 0..self.dimension {
                let element_ij = self.matrix_element(i, j);
                let element_ji = self.matrix_element(j, i);
                
                if (element_ij - element_ji.conj()).norm() > 1e-10 {
                    return false;
                }
            }
        }
        
        true
    }
}

impl fmt::Display for QuantumState {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "QuantumState[")?;
        for (i, (key, amplitude)) in self.amplitudes.iter().enumerate() {
            if i > 0 {
                write!(f, " + ")?;
            }
            write!(f, "({:.3}+{:.3}i)|{}>", amplitude.re, amplitude.im, key)?;
        }
        write!(f, "]")
    }
}

impl fmt::Display for PPFHilbertSpace {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        writeln!(f, "PPF Hilbert Space:")?;
        writeln!(f, "  Dimension: {}", self.dimension())?;
        writeln!(f, "  State space: S({})", self.state_space.value())?;
        writeln!(f, "  Basis states: {}", self.basis_states.len())?;
        Ok(())
    }
}

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

    #[test]
    fn test_quantum_state_creation() {
        let mut amplitudes = HashMap::new();
        amplitudes.insert("2 × 3".to_string(), Complex64::new(0.6, 0.0));
        amplitudes.insert("-2 × -3".to_string(), Complex64::new(0.8, 0.0));
        
        let state = QuantumState::new(amplitudes).unwrap();
        assert!(state.is_normalized());
    }

    #[test]
    fn test_basis_state() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let factorization = &state_space.factorizations()[0];
        
        let basis_state = QuantumState::basis_state(factorization, state_space.size());
        assert_eq!(basis_state.dimension(), state_space.size());
        assert!(basis_state.is_normalized());
        assert_eq!(basis_state.amplitude(factorization), Complex64::new(1.0, 0.0));
    }

    #[test]
    fn test_uniform_superposition() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let superposition = QuantumState::uniform_superposition(&state_space).unwrap();
        
        assert_eq!(superposition.dimension(), state_space.size());
        assert!(superposition.is_normalized());
        
        // All amplitudes should have equal magnitude
        let expected_magnitude = 1.0 / (state_space.size() as f64).sqrt();
        for factorization in state_space.factorizations() {
            let amplitude = superposition.amplitude(factorization);
            assert!((amplitude.norm() - expected_magnitude).abs() < 1e-10);
        }
    }

    #[test]
    fn test_state_operations() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let state1 = QuantumState::uniform_superposition(&state_space).unwrap();
        let state2 = QuantumState::basis_state(&state_space.factorizations()[0], state_space.size());
        
        let sum = state1.add(&state2).unwrap();
        assert_eq!(sum.dimension(), state_space.size());
        
        let scalar_mult = state1.scalar_mult(Complex64::new(2.0, 0.0));
        assert_eq!(scalar_mult.dimension(), state_space.size());
    }

    #[test]
    fn test_hilbert_space_creation() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let hilbert_space = PPFHilbertSpace::new(state_space.clone());
        
        assert_eq!(hilbert_space.dimension(), state_space.size());
        assert_eq!(hilbert_space.basis_states().len(), state_space.size());
    }

    #[test]
    fn test_ppf_inner_product() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let hilbert_space = PPFHilbertSpace::new(state_space.clone());
        
        let state1 = QuantumState::basis_state(&state_space.factorizations()[0], state_space.size());
        let state2 = QuantumState::basis_state(&state_space.factorizations()[0], state_space.size());
        
        let inner_product = hilbert_space.inner_product(&state1, &state2);
        
        // Should be e^(-Ω(f)/2) for the same factorization
        let complexity = state_space.factorizations()[0].complexity() as f64;
        let expected = (-complexity / 2.0).exp();
        assert!((inner_product.re - expected).abs() < 1e-10);
        assert!(inner_product.im.abs() < 1e-10);
    }

    #[test]
    fn test_sign_operator() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let hilbert_space = PPFHilbertSpace::new(state_space.clone());
        
        let sign_op = hilbert_space.sign_operator();
        assert_eq!(sign_op.dimension(), hilbert_space.dimension());
        
        // Test on a basis state - sign operator should produce zero result
        // for positive integer state space since there are no negative factorizations
        let basis_state = QuantumState::basis_state(&state_space.factorizations()[0], state_space.size());
        let result = sign_op.apply(&basis_state);
        
        // The sign operator should return a zero state for positive factorizations
        // when there are no corresponding negative factorizations in the state space
        let result_state = result.unwrap();
        
        // Check if it's a zero state (only contains zero amplitude)
        if result_state.amplitudes().len() == 1 && result_state.amplitudes().contains_key("zero_state") {
            let zero_amp = result_state.amplitudes()["zero_state"];
            assert_eq!(zero_amp, Complex64::new(0.0, 0.0));
        } else {
            // Or it should have no non-zero amplitudes
            for (_, &amp) in result_state.amplitudes() {
                assert!(amp.norm() < 1e-15);
            }
        }
    }

    #[test]
    fn test_complexity_operator() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let hilbert_space = PPFHilbertSpace::new(state_space.clone());
        
        let complexity_op = hilbert_space.complexity_operator();
        assert_eq!(complexity_op.dimension(), hilbert_space.dimension());
        
        // Should be diagonal with complexity eigenvalues
        for i in 0..hilbert_space.dimension() {
            let factorization = &state_space.factorizations()[i];
            let expected_complexity = factorization.complexity() as f64;
            let matrix_element = complexity_op.matrix_element(i, i);
            assert!((matrix_element.re - expected_complexity).abs() < 1e-10);
        }
    }

    #[test]
    fn test_curvature_flip_operator() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let hilbert_space = PPFHilbertSpace::new(state_space.clone());
        let iot_metric = IOTMetric::from_state_space(state_space);
        
        let curvature_op = hilbert_space.curvature_flip_operator(&iot_metric);
        assert_eq!(curvature_op.dimension(), hilbert_space.dimension());
        
        // Should be Hermitian
        assert!(curvature_op.is_hermitian());
    }

    #[test]
    fn test_su2_generators() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let hilbert_space = PPFHilbertSpace::new(state_space);
        
        let (sigma_x, sigma_y, sigma_z) = hilbert_space.su2_generators();
        
        // All should be Hermitian
        assert!(sigma_x.is_hermitian());
        assert!(sigma_y.is_hermitian());
        assert!(sigma_z.is_hermitian());
        
        // Check dimensions
        assert_eq!(sigma_x.dimension(), hilbert_space.dimension());
        assert_eq!(sigma_y.dimension(), hilbert_space.dimension());
        assert_eq!(sigma_z.dimension(), hilbert_space.dimension());
    }

    #[test]
    fn test_operator_operations() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let hilbert_space = PPFHilbertSpace::new(state_space);
        
        let sign_op = hilbert_space.sign_operator();
        let complexity_op = hilbert_space.complexity_operator();
        
        // Test commutator
        let commutator = sign_op.commutator(&complexity_op);
        assert_eq!(commutator.dimension(), hilbert_space.dimension());
        
        // Test anticommutator
        let anticommutator = sign_op.anticommutator(&complexity_op);
        assert_eq!(anticommutator.dimension(), hilbert_space.dimension());
        
        // Test multiplication
        let product = sign_op.multiply(&complexity_op);
        assert_eq!(product.dimension(), hilbert_space.dimension());
        
        // Test trace
        let trace = sign_op.trace();
        assert!(trace.is_finite());
    }

    #[test]
    fn test_operator_application() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let hilbert_space = PPFHilbertSpace::new(state_space.clone());
        
        let superposition = QuantumState::uniform_superposition(&state_space).unwrap();
        let complexity_op = hilbert_space.complexity_operator();
        
        let result = complexity_op.apply(&superposition).unwrap();
        assert_eq!(result.dimension(), hilbert_space.dimension());
    }

    #[test]
    fn test_inner_product_metric() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let metric = InnerProductMetric::new(&state_space);
        
        for factorization in state_space.factorizations() {
            let weight = metric.weight(factorization);
            assert!(weight > 0.0);
            assert!(weight <= 1.0);
            
            // Weight should be e^(-Ω(f)/2)
            let complexity = factorization.complexity() as f64;
            let expected = (-complexity / 2.0).exp();
            assert!((weight - expected).abs() < 1e-10);
        }
    }

    #[test]
    fn test_display() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let hilbert_space = PPFHilbertSpace::new(state_space.clone());
        
        let hilbert_str = format!("{}", hilbert_space);
        assert!(hilbert_str.contains("PPF Hilbert Space"));
        assert!(hilbert_str.contains("Dimension"));
        
        let superposition = QuantumState::uniform_superposition(&state_space).unwrap();
        let state_str = format!("{}", superposition);
        assert!(state_str.contains("QuantumState"));
    }
}