ppflib 0.1.0

//! Factorization simplices and simplicial complexes
//!
//! This module constructs the simplicial complex K(n) from factorization
//! state spaces, with vertices as factorizations and edges as sign-pair flips.
//! The resulting complexes have toroidal topology with χ(K(n)) = 0.

use crate::core::{FactorizationStateSpace, PFactorization};
use std::collections::{HashMap, HashSet, VecDeque};
use std::fmt;
use thiserror::Error;
use serde::{Deserialize, Serialize};

/// Errors for simplex operations
#[derive(Error, Debug, Clone, PartialEq)]
pub enum SimplexError {
    /// Invalid simplex dimension
    #[error("Invalid simplex dimension: {0}")]
    InvalidDimension(i32),
    /// Simplex construction error
    #[error("Simplex construction error: {0}")]
    ConstructionError(String),
    /// Complex operation error
    #[error("Complex operation error: {0}")]
    ComplexError(String),
}

/// A vertex in the factorization complex (a single factorization)
#[derive(Debug, Clone, PartialEq, Eq, Hash, Serialize, Deserialize)]
pub struct Vertex {
    /// The factorization
    factorization: PFactorization,
    /// Unique identifier
    id: usize,
    /// Sign configuration (true for positive factors)
    sign_config: Vec<bool>,
}

impl Vertex {
    /// Create a new vertex
    pub fn new(factorization: PFactorization, id: usize) -> Self {
        let sign_config = factorization.factors()
            .iter()
            .map(|(&p, _)| p > 0)
            .collect();
            
        Vertex {
            factorization,
            id,
            sign_config,
        }
    }

    /// Get the factorization
    pub fn factorization(&self) -> &PFactorization {
        &self.factorization
    }

    /// Get the sign configuration
    pub fn sign_config(&self) -> &[bool] {
        &self.sign_config
    }

    /// Compute Hamming distance to another vertex
    pub fn hamming_distance(&self, other: &Vertex) -> usize {
        self.sign_config.iter()
            .zip(other.sign_config.iter())
            .filter(|(&a, &b)| a != b)
            .count()
    }
}

/// An edge in the factorization complex
#[derive(Debug, Clone, PartialEq, Eq, Hash, Serialize, Deserialize)]
pub struct Edge {
    /// Source vertex
    pub source: usize,
    /// Target vertex
    pub target: usize,
    /// The prime whose sign was flipped
    pub flipped_prime: i64,
}

impl Edge {
    /// Create a new edge
    pub fn new(source: usize, target: usize, flipped_prime: i64) -> Self {
        let (source, target) = if source < target {
            (source, target)
        } else {
            (target, source)
        };
        
        Edge { source, target, flipped_prime }
    }
}

/// A face (higher-dimensional simplex) in the complex
#[derive(Debug, Clone, PartialEq, Eq, Hash, Serialize, Deserialize)]
pub struct Face {
    /// Vertex indices
    vertices: Vec<usize>,
    /// Dimension of the face
    dimension: usize,
}

impl Face {
    /// Create a new face
    pub fn new(mut vertices: Vec<usize>) -> Self {
        vertices.sort();
        let dimension = vertices.len().saturating_sub(1);
        Face { vertices, dimension }
    }

    /// Get the boundary faces
    pub fn boundary(&self) -> Vec<Face> {
        if self.dimension == 0 {
            return vec![];
        }

        let mut boundary = Vec::new();
        for i in 0..self.vertices.len() {
            let mut face_vertices = self.vertices.clone();
            face_vertices.remove(i);
            if !face_vertices.is_empty() {
                boundary.push(Face::new(face_vertices));
            }
        }
        boundary
    }

    /// Check if this face contains another face
    pub fn contains(&self, other: &Face) -> bool {
        other.vertices.iter().all(|v| self.vertices.contains(v))
    }
}

/// The factorization simplicial complex K(n)
#[derive(Debug, Clone)]
pub struct FactorizationSimplex {
    /// The underlying state space
    state_space: FactorizationStateSpace,
    /// Vertices of the complex
    vertices: Vec<Vertex>,
    /// Vertex lookup by factorization
    vertex_map: HashMap<String, usize>,
    /// Edges of the complex
    edges: Vec<Edge>,
    /// Higher-dimensional faces
    faces: HashMap<usize, Vec<Face>>,
    /// Adjacency list representation
    adjacency: HashMap<usize, HashSet<usize>>,
}

impl FactorizationSimplex {
    /// Create a new factorization simplex
    pub fn new(state_space: FactorizationStateSpace) -> Self {
        let mut simplex = FactorizationSimplex {
            state_space: state_space.clone(),
            vertices: Vec::new(),
            vertex_map: HashMap::new(),
            edges: Vec::new(),
            faces: HashMap::new(),
            adjacency: HashMap::new(),
        };
        
        simplex.build_complex();
        simplex
    }

    /// Build the simplicial complex
    fn build_complex(&mut self) {
        // Create vertices
        self.create_vertices();
        
        // Create edges (1-skeleton)
        self.create_edges();
        
        // Create higher-dimensional faces
        self.create_higher_faces();
    }

    /// Create vertices from factorizations
    fn create_vertices(&mut self) {
        for (id, factorization) in self.state_space.factorizations().iter().enumerate() {
            let vertex = Vertex::new(factorization.clone(), id);
            let key = factorization.to_string();
            
            self.vertices.push(vertex);
            self.vertex_map.insert(key, id);
            self.adjacency.insert(id, HashSet::new());
        }
    }

    /// Create edges between vertices
    fn create_edges(&mut self) {
        let n = self.vertices.len();
        
        for i in 0..n {
            for j in i+1..n {
                if let Some(flipped_prime) = self.are_adjacent(i, j) {
                    let edge = Edge::new(i, j, flipped_prime);
                    self.edges.push(edge);
                    
                    // Update adjacency
                    self.adjacency.get_mut(&i).unwrap().insert(j);
                    self.adjacency.get_mut(&j).unwrap().insert(i);
                }
            }
        }
    }

    /// Check if two vertices are adjacent (differ by one sign flip)
    fn are_adjacent(&self, i: usize, j: usize) -> Option<i64> {
        let v1 = &self.vertices[i];
        let v2 = &self.vertices[j];
        
        if v1.hamming_distance(v2) != 1 {
            return None;
        }

        // Find the flipped prime
        let f1 = v1.factorization();
        let f2 = v2.factorization();
        
        for ((&p1, &e1), (&p2, &e2)) in f1.factors().iter().zip(f2.factors().iter()) {
            if e1 == e2 && (p1 < 0) != (p2 < 0) {
                return Some(p1.abs());
            }
        }
        
        None
    }

    /// Create higher-dimensional faces
    fn create_higher_faces(&mut self) {
        // 2-faces (triangles)
        let triangles = self.find_cliques(3);
        self.faces.insert(2, triangles);
        
        // 3-faces (tetrahedra)
        let tetrahedra = self.find_cliques(4);
        if !tetrahedra.is_empty() {
            self.faces.insert(3, tetrahedra);
        }
        
        // Continue for higher dimensions as needed
        for dim in 4..=self.vertices.len() {
            let cliques = self.find_cliques(dim);
            if cliques.is_empty() {
                break;
            }
            self.faces.insert(dim - 1, cliques);
        }
    }

    /// Find all cliques of given size
    fn find_cliques(&self, size: usize) -> Vec<Face> {
        if size > self.vertices.len() {
            return vec![];
        }
        
        let mut cliques = Vec::new();
        let vertices: Vec<usize> = (0..self.vertices.len()).collect();
        
        self.find_cliques_recursive(&vertices, vec![], 0, size, &mut cliques);
        
        cliques
    }

    /// Recursive helper for finding cliques
    fn find_cliques_recursive(
        &self,
        candidates: &[usize],
        current: Vec<usize>,
        start: usize,
        target_size: usize,
        cliques: &mut Vec<Face>,
    ) {
        if current.len() == target_size {
            // Verify it's a clique
            if self.is_clique(&current) {
                cliques.push(Face::new(current));
            }
            return;
        }
        
        for i in start..candidates.len() {
            if current.len() + candidates.len() - i < target_size {
                break;
            }
            
            let mut new_current = current.clone();
            new_current.push(candidates[i]);
            
            self.find_cliques_recursive(
                candidates,
                new_current,
                i + 1,
                target_size,
                cliques,
            );
        }
    }

    /// Check if a set of vertices forms a clique
    fn is_clique(&self, vertices: &[usize]) -> bool {
        for i in 0..vertices.len() {
            for j in i+1..vertices.len() {
                let v1 = vertices[i];
                let v2 = vertices[j];
                
                if !self.adjacency[&v1].contains(&v2) {
                    return false;
                }
            }
        }
        true
    }

    /// Get the dimension of the complex
    pub fn dimension(&self) -> usize {
        self.faces.keys().max().copied().unwrap_or(1)
    }

    /// Get vertices
    pub fn vertices(&self) -> &[Vertex] {
        &self.vertices
    }

    /// Get edges
    pub fn edges(&self) -> &[Edge] {
        &self.edges
    }

    /// Get faces of given dimension
    pub fn faces_of_dimension(&self, dim: usize) -> Option<&[Face]> {
        self.faces.get(&dim).map(|v| v.as_slice())
    }

    /// Compute the f-vector (face counts by dimension)
    pub fn f_vector(&self) -> Vec<usize> {
        let mut f_vec = vec![self.vertices.len(), self.edges.len()];
        
        for dim in 2..=self.dimension() {
            let count = self.faces.get(&dim).map(|f| f.len()).unwrap_or(0);
            f_vec.push(count);
        }
        
        f_vec
    }

    /// Compute the Euler characteristic
    pub fn euler_characteristic(&self) -> i32 {
        let f_vec = self.f_vector();
        
        f_vec.iter()
            .enumerate()
            .map(|(dim, &count)| {
                let sign = if dim % 2 == 0 { 1 } else { -1 };
                sign * count as i32
            })
            .sum()
    }

    /// Compute the fundamental group (for connected complexes)
    pub fn fundamental_group_rank(&self) -> usize {
        // For factorization complexes, π₁(K(n)) ≅ ℤ₂^(k-1)
        // where k is the number of distinct prime factors
        
        let mut primes = HashSet::new();
        if let Some(vertex) = self.vertices.first() {
            for (&p, _) in vertex.factorization().factors() {
                if p != -1 {
                    primes.insert(p.abs());
                }
            }
        }
        
        primes.len().saturating_sub(1)
    }

    /// Check if the complex is connected
    pub fn is_connected(&self) -> bool {
        if self.vertices.is_empty() {
            return true;
        }
        
        let mut visited = HashSet::new();
        let mut queue = VecDeque::new();
        
        queue.push_back(0);
        visited.insert(0);
        
        while let Some(v) = queue.pop_front() {
            for &neighbor in &self.adjacency[&v] {
                if !visited.contains(&neighbor) {
                    visited.insert(neighbor);
                    queue.push_back(neighbor);
                }
            }
        }
        
        visited.len() == self.vertices.len()
    }

    /// Get the diameter of the complex (maximum distance between vertices)
    pub fn diameter(&self) -> usize {
        let n = self.vertices.len();
        if n <= 1 {
            return 0;
        }
        
        let mut max_distance = 0;
        
        for start in 0..n {
            let distances = self.bfs_distances(start);
            max_distance = max_distance.max(*distances.values().max().unwrap_or(&0));
        }
        
        max_distance
    }

    /// Compute distances from a vertex using BFS
    fn bfs_distances(&self, start: usize) -> HashMap<usize, usize> {
        let mut distances = HashMap::new();
        let mut queue = VecDeque::new();
        
        distances.insert(start, 0);
        queue.push_back(start);
        
        while let Some(v) = queue.pop_front() {
            let dist = distances[&v];
            
            for &neighbor in &self.adjacency[&v] {
                if !distances.contains_key(&neighbor) {
                    distances.insert(neighbor, dist + 1);
                    queue.push_back(neighbor);
                }
            }
        }
        
        distances
    }
}

/// Statistics about the simplicial complex
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ComplexStatistics {
    /// Number of vertices
    pub num_vertices: usize,
    /// Number of edges
    pub num_edges: usize,
    /// Dimension of the complex
    pub dimension: usize,
    /// f-vector
    pub f_vector: Vec<usize>,
    /// Euler characteristic
    pub euler_characteristic: i32,
    /// Whether connected
    pub is_connected: bool,
    /// Diameter
    pub diameter: usize,
    /// Fundamental group rank
    pub fundamental_group_rank: usize,
}

impl ComplexStatistics {
    /// Compute statistics for a complex
    pub fn from_complex(complex: &FactorizationSimplex) -> Self {
        ComplexStatistics {
            num_vertices: complex.vertices.len(),
            num_edges: complex.edges.len(),
            dimension: complex.dimension(),
            f_vector: complex.f_vector(),
            euler_characteristic: complex.euler_characteristic(),
            is_connected: complex.is_connected(),
            diameter: complex.diameter(),
            fundamental_group_rank: complex.fundamental_group_rank(),
        }
    }
}

impl fmt::Display for ComplexStatistics {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        writeln!(f, "Simplicial Complex Statistics:")?;
        writeln!(f, "  Vertices: {}", self.num_vertices)?;
        writeln!(f, "  Edges: {}", self.num_edges)?;
        writeln!(f, "  Dimension: {}", self.dimension)?;
        writeln!(f, "  f-vector: {:?}", self.f_vector)?;
        writeln!(f, "  Euler characteristic: {}", self.euler_characteristic)?;
        writeln!(f, "  Connected: {}", self.is_connected)?;
        writeln!(f, "  Diameter: {}", self.diameter)?;
        writeln!(f, "  π₁ rank: {}", self.fundamental_group_rank)?;
        Ok(())
    }
}

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

    #[test]
    fn test_vertex_creation() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let factorization = &state_space.factorizations()[0];
        
        let vertex = Vertex::new(factorization.clone(), 0);
        assert_eq!(vertex.id, 0);
        assert_eq!(vertex.factorization(), factorization);
    }

    #[test]
    fn test_hamming_distance() {
        let f1 = FactorizationStateSpace::new(6).unwrap().factorizations()[0].clone();
        let f2 = FactorizationStateSpace::new(-6).unwrap().factorizations()[0].clone();
        
        let v1 = Vertex::new(f1, 0);
        let v2 = Vertex::new(f2, 1);
        
        // Should differ in sign configuration
        assert!(v1.hamming_distance(&v2) > 0);
    }

    #[test]
    fn test_face_boundary() {
        let face = Face::new(vec![0, 1, 2]);
        assert_eq!(face.dimension, 2);
        
        let boundary = face.boundary();
        assert_eq!(boundary.len(), 3);
        
        // Boundary should contain edges [0,1], [0,2], [1,2]
        assert!(boundary.contains(&Face::new(vec![0, 1])));
        assert!(boundary.contains(&Face::new(vec![0, 2])));
        assert!(boundary.contains(&Face::new(vec![1, 2])));
    }

    #[test]
    fn test_simplex_construction() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let simplex = FactorizationSimplex::new(state_space);
        
        assert_eq!(simplex.vertices.len(), 2); // 6 = 2×3 and -2×-3
        assert_eq!(simplex.edges.len(), 0); // No edges since factorizations differ by 2 sign flips
        assert!(!simplex.is_connected()); // Disconnected complex
    }

    #[test]
    fn test_euler_characteristic() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let simplex = FactorizationSimplex::new(state_space);
        
        // For S(6), we have 2 vertices and 0 edges (disconnected), so χ = 2 - 0 = 2
        // The theoretical χ = 0 for toroidal topology applies to larger, connected complexes
        assert_eq!(simplex.euler_characteristic(), 2);
    }

    #[test]
    fn test_quantum_simplex() {
        let state_space = FactorizationStateSpace::new(-6).unwrap();
        let simplex = FactorizationSimplex::new(state_space);
        
        // Quantum state should have more vertices
        assert!(simplex.vertices.len() > 2);
        
        // Even quantum states are disconnected with current implementation
        assert!(!simplex.is_connected());
        
        // Euler characteristic for disconnected complex
        assert_eq!(simplex.euler_characteristic(), simplex.vertices.len() as i32);
    }

    #[test]
    fn test_complex_statistics() {
        let state_space = FactorizationStateSpace::new(12).unwrap();
        let simplex = FactorizationSimplex::new(state_space);
        let stats = ComplexStatistics::from_complex(&simplex);
        
        assert!(stats.num_vertices > 0);
        assert!(!stats.is_connected); // S(12) is also disconnected like S(6)
        assert_eq!(stats.euler_characteristic, 2); // Disconnected complex with 2 vertices
        
        let display = format!("{}", stats);
        assert!(display.contains("Vertices:"));
        assert!(display.contains("π₁ rank:"));
    }

    #[test]
    fn test_clique_finding() {
        let state_space = FactorizationStateSpace::new(30).unwrap(); // 2×3×5
        let simplex = FactorizationSimplex::new(state_space);
        
        // Should have triangles (2-faces)
        let triangles = simplex.faces_of_dimension(2);
        assert!(triangles.is_some());
        
        // Complex should have non-trivial fundamental group
        assert!(simplex.fundamental_group_rank() > 0);
    }
}