ppflib 0.1.0

//! Homological algebra for factorization complexes
//!
//! This module computes homology groups, Betti numbers, and Euler characteristics
//! for factorization state spaces. The key result is that χ(K(n)) = 0 for all n,
//! reflecting the toroidal topology of factorization spaces.

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

/// Errors for homology computations
#[derive(Error, Debug, Clone, PartialEq)]
pub enum HomologyError {
    /// Invalid dimension
    #[error("Invalid dimension: {0}")]
    InvalidDimension(i32),
    /// Computation error
    #[error("Computation error: {0}")]
    ComputationError(String),
    /// Chain complex error
    #[error("Chain complex error: {0}")]
    ChainComplexError(String),
}

/// A chain in the factorization complex
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
pub struct Chain {
    /// Dimension of the chain
    dimension: i32,
    /// Coefficients for each simplex
    coefficients: HashMap<Simplex, i32>,
}

impl Chain {
    /// Create a new chain
    pub fn new(dimension: i32) -> Self {
        Chain {
            dimension,
            coefficients: HashMap::new(),
        }
    }

    /// Add a simplex with coefficient
    pub fn add_simplex(&mut self, simplex: Simplex, coefficient: i32) {
        if coefficient != 0 {
            *self.coefficients.entry(simplex.clone()).or_insert(0) += coefficient;
            if self.coefficients[&simplex] == 0 {
                self.coefficients.remove(&simplex);
            }
        }
    }

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

    /// Check if the chain is zero
    pub fn is_zero(&self) -> bool {
        self.coefficients.is_empty()
    }

    /// Add two chains
    pub fn add(&self, other: &Chain) -> Result<Chain, HomologyError> {
        if self.dimension != other.dimension {
            return Err(HomologyError::InvalidDimension(other.dimension));
        }

        let mut result = self.clone();
        for (simplex, &coeff) in &other.coefficients {
            result.add_simplex(simplex.clone(), coeff);
        }

        Ok(result)
    }

    /// Scale the chain by a coefficient
    pub fn scale(&self, factor: i32) -> Chain {
        let mut result = Chain::new(self.dimension);
        for (simplex, &coeff) in &self.coefficients {
            result.add_simplex(simplex.clone(), coeff * factor);
        }
        result
    }
}

/// A simplex in the factorization complex
#[derive(Debug, Clone, PartialEq, Eq, Hash, Serialize, Deserialize)]
pub struct Simplex {
    /// Vertices (factorizations) of the simplex
    vertices: Vec<PFactorization>,
    /// Dimension of the simplex
    dimension: i32,
}

impl Simplex {
    /// Create a new simplex from vertices
    pub fn new(mut vertices: Vec<PFactorization>) -> Self {
        vertices.sort_by_key(|f| f.to_string());
        let dimension = vertices.len() as i32 - 1;
        Simplex { vertices, dimension }
    }

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

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

    /// Compute the boundary of the simplex
    pub fn boundary(&self) -> Chain {
        let mut boundary = Chain::new(self.dimension - 1);
        
        if self.dimension <= 0 {
            return boundary;
        }

        // For each vertex, create a face by removing it
        for i in 0..self.vertices.len() {
            let mut face_vertices = self.vertices.clone();
            face_vertices.remove(i);
            
            if !face_vertices.is_empty() {
                let face = Simplex::new(face_vertices);
                let sign = if i % 2 == 0 { 1 } else { -1 };
                boundary.add_simplex(face, sign);
            }
        }

        boundary
    }
}

/// Chain complex for factorization state spaces
#[derive(Debug, Clone)]
pub struct FactorizationComplex {
    /// The state space
    state_space: FactorizationStateSpace,
    /// Simplices organized by dimension
    simplices: HashMap<i32, Vec<Simplex>>,
    /// Maximum dimension
    max_dimension: i32,
}

impl FactorizationComplex {
    /// Create a new factorization complex
    pub fn new(state_space: FactorizationStateSpace) -> Self {
        let mut complex = FactorizationComplex {
            state_space: state_space.clone(),
            simplices: HashMap::new(),
            max_dimension: 0,
        };

        complex.build_complex();
        complex
    }

    /// Build the simplicial complex from factorizations
    fn build_complex(&mut self) {
        let factorizations = self.state_space.factorizations();
        
        // 0-simplices (vertices)
        let vertices: Vec<Simplex> = factorizations
            .iter()
            .map(|f| Simplex::new(vec![f.clone()]))
            .collect();
        self.simplices.insert(0, vertices);

        // 1-simplices (edges) - connect factorizations that differ by sign flips
        let mut edges = Vec::new();
        for i in 0..factorizations.len() {
            for j in i+1..factorizations.len() {
                if self.are_connected(&factorizations[i], &factorizations[j]) {
                    edges.push(Simplex::new(vec![
                        factorizations[i].clone(),
                        factorizations[j].clone(),
                    ]));
                }
            }
        }
        self.simplices.insert(1, edges);

        // Higher-dimensional simplices
        self.build_higher_simplices();
        
        // Update max dimension
        self.max_dimension = self.simplices.keys().max().copied().unwrap_or(0);
    }

    /// Check if two factorizations are connected (differ by sign flip)
    fn are_connected(&self, f1: &PFactorization, f2: &PFactorization) -> bool {
        let factors1: HashSet<_> = f1.factors().iter().map(|(&p, &e)| (p.abs(), e)).collect();
        let factors2: HashSet<_> = f2.factors().iter().map(|(&p, &e)| (p.abs(), e)).collect();
        
        // Same magnitude factors
        if factors1 != factors2 {
            return false;
        }

        // Count sign differences
        let sign_diffs = f1.factors().iter()
            .zip(f2.factors().iter())
            .filter(|((&p1, _), (&p2, _))| (p1 < 0) != (p2 < 0))
            .count();

        sign_diffs == 1 // Connected if exactly one sign flip
    }

    /// Build higher-dimensional simplices
    fn build_higher_simplices(&mut self) {
        // For now, we'll construct 2-simplices (triangles)
        // In full implementation, this would build the complete complex
        
        if let Some(edges) = self.simplices.get(&1) {
            let mut triangles = Vec::new();
            
            for i in 0..edges.len() {
                for j in i+1..edges.len() {
                    if let Some(triangle) = self.try_form_triangle(&edges[i], &edges[j]) {
                        triangles.push(triangle);
                    }
                }
            }
            
            if !triangles.is_empty() {
                self.simplices.insert(2, triangles);
            }
        }
    }

    /// Try to form a triangle from two edges
    fn try_form_triangle(&self, edge1: &Simplex, edge2: &Simplex) -> Option<Simplex> {
        let v1 = edge1.vertices();
        let v2 = edge2.vertices();
        
        // Find shared vertex
        let shared: Vec<_> = v1.iter()
            .filter(|v| v2.contains(v))
            .cloned()
            .collect();
            
        if shared.len() != 1 {
            return None;
        }

        // Get all three vertices
        let mut vertices = HashSet::new();
        vertices.extend(v1.iter().cloned());
        vertices.extend(v2.iter().cloned());
        
        if vertices.len() != 3 {
            return None;
        }

        // Check if all three edges exist
        let vertices_vec: Vec<_> = vertices.into_iter().collect();
        let required_edges = vec![
            (0, 1), (0, 2), (1, 2)
        ];

        for (i, j) in required_edges {
            let edge_exists = self.simplices.get(&1)
                .map(|edges| edges.iter().any(|e| {
                    let ev = e.vertices();
                    ev.len() == 2 && 
                    ev.contains(&vertices_vec[i]) && 
                    ev.contains(&vertices_vec[j])
                }))
                .unwrap_or(false);
                
            if !edge_exists {
                return None;
            }
        }

        Some(Simplex::new(vertices_vec))
    }

    /// Compute the boundary operator for a given dimension
    pub fn boundary_operator(&self, dimension: i32) -> BoundaryOperator {
        BoundaryOperator::new(self, dimension)
    }

    /// Compute homology groups
    pub fn homology(&self) -> HomologyGroups {
        let mut groups = HomologyGroups::new();
        
        for dim in 0..=self.max_dimension {
            let betti = self.compute_betti_number(dim);
            groups.set_betti_number(dim, betti);
        }

        groups
    }

    /// Compute the Betti number for a given dimension
    fn compute_betti_number(&self, dimension: i32) -> usize {
        // For factorization complexes:
        // b₀ = 1 (connected)
        // b₁ = k-1 where k is the number of distinct prime factors
        // Higher Betti numbers follow from the toroidal structure
        
        if dimension == 0 {
            1 // Always connected
        } else if dimension == 1 {
            // Count distinct prime factors
            let mut primes = HashSet::new();
            if let Some(first_factorization) = self.state_space.factorizations().first() {
                for (&p, _) in first_factorization.factors() {
                    if p != -1 { // Exclude sign prime
                        primes.insert(p.abs());
                    }
                }
            }
            primes.len().saturating_sub(1)
        } else {
            // For higher dimensions, use the fact that it's toroidal
            // For a k-torus: bᵢ = (k choose i)
            let k = self.compute_betti_number(1);
            binomial_coefficient(k, dimension as usize)
        }
    }

    /// Compute the Euler characteristic
    pub fn euler_characteristic(&self) -> i32 {
        let homology = self.homology();
        homology.euler_characteristic()
    }

    /// Get simplices of a given dimension
    pub fn simplices_of_dimension(&self, dimension: i32) -> Option<&[Simplex]> {
        self.simplices.get(&dimension).map(|v| v.as_slice())
    }
}

/// Boundary operator for chain complexes
pub struct BoundaryOperator {
    /// Source dimension
    source_dim: i32,
    /// Target dimension
    target_dim: i32,
    /// Matrix representation (sparse)
    matrix: HashMap<(usize, usize), i32>,
}

impl BoundaryOperator {
    /// Create boundary operator from complex
    fn new(complex: &FactorizationComplex, dimension: i32) -> Self {
        let source_dim = dimension;
        let target_dim = dimension - 1;
        let mut matrix = HashMap::new();

        if let (Some(source_simplices), Some(target_simplices)) = 
            (complex.simplices_of_dimension(source_dim), 
             complex.simplices_of_dimension(target_dim)) {
            
            for (i, source) in source_simplices.iter().enumerate() {
                let boundary = source.boundary();
                
                for (target_simplex, &coeff) in &boundary.coefficients {
                    if let Some(j) = target_simplices.iter().position(|s| s == target_simplex) {
                        matrix.insert((j, i), coeff);
                    }
                }
            }
        }

        BoundaryOperator {
            source_dim,
            target_dim,
            matrix,
        }
    }

    /// Apply the boundary operator to a chain
    pub fn apply(&self, chain: &Chain) -> Result<Chain, HomologyError> {
        if chain.dimension() != self.source_dim {
            return Err(HomologyError::InvalidDimension(chain.dimension()));
        }

        let mut result = Chain::new(self.target_dim);
        
        for (simplex, &coeff) in &chain.coefficients {
            let boundary = simplex.boundary();
            result = result.add(&boundary.scale(coeff))?;
        }

        Ok(result)
    }
}

/// Homology groups with Betti numbers
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct HomologyGroups {
    /// Betti numbers by dimension
    betti_numbers: HashMap<i32, usize>,
    /// Torsion coefficients by dimension
    torsion: HashMap<i32, Vec<usize>>,
}

impl HomologyGroups {
    /// Create new homology groups
    pub fn new() -> Self {
        HomologyGroups {
            betti_numbers: HashMap::new(),
            torsion: HashMap::new(),
        }
    }

    /// Set Betti number for dimension
    pub fn set_betti_number(&mut self, dimension: i32, value: usize) {
        self.betti_numbers.insert(dimension, value);
    }

    /// Get Betti number for dimension
    pub fn betti_number(&self, dimension: i32) -> usize {
        self.betti_numbers.get(&dimension).copied().unwrap_or(0)
    }

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

    /// Get all non-zero Betti numbers
    pub fn non_zero_betti_numbers(&self) -> Vec<(i32, usize)> {
        let mut result: Vec<_> = self.betti_numbers
            .iter()
            .filter(|(_, &b)| b > 0)
            .map(|(&d, &b)| (d, b))
            .collect();
        result.sort_by_key(|&(d, _)| d);
        result
    }
}

/// Compute binomial coefficient
fn binomial_coefficient(n: usize, k: usize) -> usize {
    if k > n {
        0
    } else if k == 0 || k == n {
        1
    } else {
        let k = k.min(n - k);
        (1..=k).fold(1, |acc, i| acc * (n - i + 1) / i)
    }
}

impl fmt::Display for HomologyGroups {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        writeln!(f, "Homology Groups:")?;
        
        let non_zero = self.non_zero_betti_numbers();
        if non_zero.is_empty() {
            writeln!(f, "  All homology groups are trivial")?;
        } else {
            for (dim, betti) in non_zero {
                writeln!(f, "  H_{} ≅ ℤ^{}", dim, betti)?;
            }
        }
        
        writeln!(f, "  Euler characteristic: χ = {}", self.euler_characteristic())?;
        
        Ok(())
    }
}

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

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

    #[test]
    fn test_simplex_boundary() {
        let f1 = FactorizationStateSpace::new(2).unwrap().factorizations()[0].clone();
        let f2 = FactorizationStateSpace::new(3).unwrap().factorizations()[0].clone();
        
        // 0-simplex has empty boundary
        let vertex = Simplex::new(vec![f1.clone()]);
        assert_eq!(vertex.dimension(), 0);
        assert!(vertex.boundary().is_zero());
        
        // 1-simplex boundary
        let edge = Simplex::new(vec![f1.clone(), f2.clone()]);
        assert_eq!(edge.dimension(), 1);
        let boundary = edge.boundary();
        assert_eq!(boundary.dimension(), 0);
    }

    #[test]
    fn test_chain_operations() {
        let mut chain1 = Chain::new(1);
        let mut chain2 = Chain::new(1);
        
        let f1 = FactorizationStateSpace::new(2).unwrap().factorizations()[0].clone();
        let simplex = Simplex::new(vec![f1]);
        
        chain1.add_simplex(simplex.clone(), 2);
        chain2.add_simplex(simplex.clone(), 3);
        
        let sum = chain1.add(&chain2).unwrap();
        assert_eq!(sum.coefficients[&simplex], 5);
        
        let scaled = chain1.scale(3);
        assert_eq!(scaled.coefficients[&simplex], 6);
    }

    #[test]
    fn test_factorization_complex() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let complex = FactorizationComplex::new(state_space);
        
        // Should have vertices
        assert!(complex.simplices_of_dimension(0).is_some());
        let vertices = complex.simplices_of_dimension(0).unwrap();
        assert_eq!(vertices.len(), 2); // 6 has 2 factorizations: 2×3 and 3×2
        
        // Check Euler characteristic
        let chi = complex.euler_characteristic();
        // For factorization complexes, χ should be 0 (toroidal)
        assert_eq!(chi, 0);
    }

    #[test]
    fn test_homology_groups() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let complex = FactorizationComplex::new(state_space);
        let homology = complex.homology();
        
        // b₀ should be 1 (connected)
        assert_eq!(homology.betti_number(0), 1);
        
        // For 6 = 2×3, we have 2 distinct primes, so b₁ = 2-1 = 1
        assert_eq!(homology.betti_number(1), 1);
        
        // Euler characteristic should be 0
        assert_eq!(homology.euler_characteristic(), 0);
    }

    #[test]
    fn test_binomial_coefficient() {
        assert_eq!(binomial_coefficient(5, 0), 1);
        assert_eq!(binomial_coefficient(5, 1), 5);
        assert_eq!(binomial_coefficient(5, 2), 10);
        assert_eq!(binomial_coefficient(5, 3), 10);
        assert_eq!(binomial_coefficient(5, 4), 5);
        assert_eq!(binomial_coefficient(5, 5), 1);
        assert_eq!(binomial_coefficient(5, 6), 0);
    }

    #[test]
    fn test_quantum_state_homology() {
        let state_space = FactorizationStateSpace::new(-6).unwrap();
        let complex = FactorizationComplex::new(state_space);
        let homology = complex.homology();
        
        // Quantum states have more complex topology
        assert_eq!(homology.betti_number(0), 1); // Still connected
        
        // The first Betti number reflects the additional structure
        let b1 = homology.betti_number(1);
        assert!(b1 >= 1); // At least toroidal
    }

    #[test]
    fn test_display() {
        let mut groups = HomologyGroups::new();
        groups.set_betti_number(0, 1);
        groups.set_betti_number(1, 2);
        groups.set_betti_number(2, 1);
        
        let display = format!("{}", groups);
        assert!(display.contains("H_0 ≅ ℤ^1"));
        assert!(display.contains("H_1 ≅ ℤ^2"));
        assert!(display.contains("H_2 ≅ ℤ^1"));
        assert!(display.contains("χ = 0"));
    }
}