ppflib 0.1.0

//! IOT (Involuted Oblate Toroidal) metric implementation
//!
//! This module implements the IOT metric with the critical ratio r/R = 1/30
//! derived from the variational principle in PPF theory. The IOT geometry
//! provides the geometric realization of factorization state spaces.

use crate::core::FactorizationStateSpace;
use std::f64::consts::PI;
use std::fmt;
use thiserror::Error;
use serde::{Deserialize, Serialize};

/// Errors for IOT metric operations
#[derive(Error, Debug, Clone, PartialEq)]
pub enum IOTError {
    /// Invalid coordinates
    #[error("Invalid coordinates: {0}")]
    InvalidCoordinates(String),
    /// Metric computation error
    #[error("Metric computation error: {0}")]
    MetricError(String),
    /// Critical ratio error
    #[error("Critical ratio error: {0}")]
    CriticalRatioError(String),
}

/// IOT coordinate system
/// 
/// The IOT uses toroidal coordinates (φ, θ, ψ) where:
/// - φ: toroidal angle (0 to 2π)
/// - θ: poloidal angle (0 to 2π)  
/// - ψ: involute parameter (related to factorization complexity)
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct IOTCoordinates {
    /// Toroidal angle φ ∈ [0, 2π)
    pub phi: f64,
    /// Poloidal angle θ ∈ [0, 2π)
    pub theta: f64,
    /// Involute parameter ψ
    pub psi: f64,
}

impl IOTCoordinates {
    /// Create new IOT coordinates
    pub fn new(phi: f64, theta: f64, psi: f64) -> Result<Self, IOTError> {
        if phi < 0.0 || phi >= 2.0 * PI {
            return Err(IOTError::InvalidCoordinates(
                format!("φ = {} not in [0, 2π)", phi)
            ));
        }
        if theta < 0.0 || theta >= 2.0 * PI {
            return Err(IOTError::InvalidCoordinates(
                format!("θ = {} not in [0, 2π)", theta)
            ));
        }
        
        Ok(IOTCoordinates { phi, theta, psi })
    }

    /// Normalize coordinates to canonical range
    pub fn normalize(&mut self) {
        self.phi = self.phi % (2.0 * PI);
        self.theta = self.theta % (2.0 * PI);
        if self.phi < 0.0 {
            self.phi += 2.0 * PI;
        }
        if self.theta < 0.0 {
            self.theta += 2.0 * PI;
        }
    }

    /// Convert to Cartesian coordinates
    pub fn to_cartesian(&self, major_radius: f64, minor_radius: f64) -> CartesianPoint {
        let major_r = major_radius;
        let minor_r = minor_radius;
        
        // Standard torus transformation with involute modification
        let involute_factor = 1.0 + self.psi * 0.1; // Small involute perturbation
        
        let x = (major_r + minor_r * self.theta.cos() * involute_factor) * self.phi.cos();
        let y = (major_r + minor_r * self.theta.cos() * involute_factor) * self.phi.sin();
        let z = minor_r * self.theta.sin() * involute_factor;
        
        CartesianPoint { x, y, z }
    }
}

/// Cartesian point in 3D space
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct CartesianPoint {
    /// X coordinate
    pub x: f64,
    /// Y coordinate
    pub y: f64,
    /// Z coordinate
    pub z: f64,
}

/// IOT metric tensor components
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct IOTMetricTensor {
    /// g_φφ component
    pub g_phi_phi: f64,
    /// g_θθ component
    pub g_theta_theta: f64,
    /// g_ψψ component
    pub g_psi_psi: f64,
    /// g_φθ component (typically 0 for orthogonal coordinates)
    pub g_phi_theta: f64,
    /// g_φψ component
    pub g_phi_psi: f64,
    /// g_θψ component
    pub g_theta_psi: f64,
}

/// The IOT metric implementation
#[derive(Debug, Clone)]
pub struct IOTMetric {
    /// Major radius R
    major_radius: f64,
    /// Minor radius r
    minor_radius: f64,
    /// Critical ratio r/R
    critical_ratio: f64,
    /// Involute parameter scale
    involute_scale: f64,
    /// Associated factorization state space
    state_space: Option<FactorizationStateSpace>,
}

impl IOTMetric {
    /// Critical ratio from PPF theory: r/R = 1/30
    pub const CRITICAL_RATIO: f64 = 1.0 / 30.0;
    
    /// Create a new IOT metric with critical ratio
    pub fn new() -> Self {
        let major_radius = 1.0;
        let minor_radius = major_radius * Self::CRITICAL_RATIO;
        
        IOTMetric {
            major_radius,
            minor_radius,
            critical_ratio: Self::CRITICAL_RATIO,
            involute_scale: 1.0,
            state_space: None,
        }
    }

    /// Create IOT metric with custom radii
    pub fn with_radii(major_radius: f64, minor_radius: f64) -> Result<Self, IOTError> {
        if major_radius <= 0.0 || minor_radius <= 0.0 {
            return Err(IOTError::InvalidCoordinates(
                "Radii must be positive".to_string()
            ));
        }
        
        let critical_ratio = minor_radius / major_radius;
        
        Ok(IOTMetric {
            major_radius,
            minor_radius,
            critical_ratio,
            involute_scale: 1.0,
            state_space: None,
        })
    }

    /// Create IOT metric from factorization state space
    pub fn from_state_space(state_space: FactorizationStateSpace) -> Self {
        let mut metric = Self::new();
        
        // Scale parameters based on state space properties
        let size = state_space.size();
        let complexity = state_space.factorizations()
            .iter()
            .map(|f| f.complexity())
            .sum::<u32>() as f64;
            
        metric.involute_scale = 1.0 + complexity / (size as f64);
        metric.state_space = Some(state_space);
        
        metric
    }

    /// Get the major radius R
    pub fn major_radius(&self) -> f64 {
        self.major_radius
    }

    /// Get the minor radius r
    pub fn minor_radius(&self) -> f64 {
        self.minor_radius
    }

    /// Get the critical ratio r/R
    pub fn critical_ratio(&self) -> f64 {
        self.critical_ratio
    }

    /// Check if this metric has the critical ratio
    pub fn is_critical(&self) -> bool {
        (self.critical_ratio - Self::CRITICAL_RATIO).abs() < 1e-10
    }

    /// Compute the metric tensor at given coordinates
    pub fn metric_tensor(&self, coords: &IOTCoordinates) -> IOTMetricTensor {
        let major_r = self.major_radius;
        let minor_r = self.minor_radius;
        let scale = self.involute_scale;
        
        // Involute modifications
        let involute_factor = 1.0 + coords.psi * 0.1 * scale;
        let involute_derivative = 0.1 * scale;
        
        // Standard torus metric with involute modifications
        let rho = major_r + minor_r * coords.theta.cos() * involute_factor;
        
        let g_phi_phi = rho * rho;
        let g_theta_theta = minor_r * minor_r * involute_factor * involute_factor;
        let g_psi_psi = (minor_r * coords.theta.cos() * involute_derivative).powi(2) + 
                       (minor_r * coords.theta.sin() * involute_derivative).powi(2);
        
        // Cross terms from involute coupling
        let g_phi_theta = 0.0; // Orthogonal in standard torus
        let g_phi_psi = 2.0 * rho * minor_r * coords.theta.cos() * involute_derivative;
        let g_theta_psi = 2.0 * minor_r * minor_r * involute_factor * 
                         (-coords.theta.sin() * involute_derivative);
        
        IOTMetricTensor {
            g_phi_phi,
            g_theta_theta,
            g_psi_psi,
            g_phi_theta,
            g_phi_psi,
            g_theta_psi,
        }
    }

    /// Compute the metric determinant
    pub fn metric_determinant(&self, coords: &IOTCoordinates) -> f64 {
        let g = self.metric_tensor(coords);
        
        // For a 3x3 symmetric matrix with the IOT structure
        // This is a simplified computation assuming small involute terms
        let det = g.g_phi_phi * g.g_theta_theta * g.g_psi_psi - 
                 g.g_phi_phi * g.g_theta_psi * g.g_theta_psi -
                 g.g_theta_theta * g.g_phi_psi * g.g_phi_psi;
        
        det.abs()
    }

    /// Compute the Ricci scalar curvature
    pub fn ricci_scalar(&self, coords: &IOTCoordinates) -> f64 {
        let major_r = self.major_radius;
        let minor_r = self.minor_radius;
        let scale = self.involute_scale;
        
        // Simplified Ricci scalar for IOT metric
        // In PPF theory, curvature arises from number-theoretic structure
        let involute_factor = 1.0 + coords.psi * 0.1 * scale;
        let rho = major_r + minor_r * coords.theta.cos() * involute_factor;
        
        // Base torus curvature
        let torus_curvature = 2.0 * coords.theta.cos() / (minor_r * rho);
        
        // PPF correction from involute geometry
        let ppf_correction = scale * (coords.psi * coords.psi) / (major_r * major_r);
        
        torus_curvature + ppf_correction
    }

    /// Compute geodesic distance between two points
    pub fn geodesic_distance(&self, p1: &IOTCoordinates, p2: &IOTCoordinates) -> f64 {
        // Simplified geodesic distance computation
        // For exact computation, would need to solve geodesic equations
        
        let dp = IOTCoordinates::new(
            (p2.phi - p1.phi).abs(),
            (p2.theta - p1.theta).abs(),
            (p2.psi - p1.psi).abs(),
        ).unwrap();
        
        let g1 = self.metric_tensor(p1);
        let g2 = self.metric_tensor(p2);
        
        // Average metric for approximate distance
        let g_avg = IOTMetricTensor {
            g_phi_phi: (g1.g_phi_phi + g2.g_phi_phi) / 2.0,
            g_theta_theta: (g1.g_theta_theta + g2.g_theta_theta) / 2.0,
            g_psi_psi: (g1.g_psi_psi + g2.g_psi_psi) / 2.0,
            g_phi_theta: (g1.g_phi_theta + g2.g_phi_theta) / 2.0,
            g_phi_psi: (g1.g_phi_psi + g2.g_phi_psi) / 2.0,
            g_theta_psi: (g1.g_theta_psi + g2.g_theta_psi) / 2.0,
        };
        
        // Approximate distance using metric
        let ds2 = g_avg.g_phi_phi * dp.phi * dp.phi +
                 g_avg.g_theta_theta * dp.theta * dp.theta +
                 g_avg.g_psi_psi * dp.psi * dp.psi +
                 2.0 * g_avg.g_phi_theta * dp.phi * dp.theta +
                 2.0 * g_avg.g_phi_psi * dp.phi * dp.psi +
                 2.0 * g_avg.g_theta_psi * dp.theta * dp.psi;
        
        ds2.abs().sqrt()
    }

    /// Compute the volume element
    pub fn volume_element(&self, coords: &IOTCoordinates) -> f64 {
        let det = self.metric_determinant(coords);
        det.sqrt()
    }

    /// Map factorization to IOT coordinates
    pub fn factorization_to_coordinates(&self, factorization: &crate::core::PFactorization) -> IOTCoordinates {
        let factors = factorization.factors();
        let complexity = factorization.complexity() as f64;
        
        // Map prime factors to angular coordinates
        let mut phi = 0.0;
        let mut theta = 0.0;
        
        for (&prime, &count) in factors {
            if prime > 0 {
                phi += prime as f64 * count as f64 * 0.1;
                theta += (prime as f64).ln() * count as f64 * 0.2;
            }
        }
        
        // Normalize angles
        phi = phi % (2.0 * PI);
        theta = theta % (2.0 * PI);
        
        // Psi from complexity
        let psi = complexity / 10.0;
        
        IOTCoordinates::new(phi, theta, psi).unwrap()
    }

    /// Compute the warping function for space-time embedding
    pub fn warping_function(&self, coords: &IOTCoordinates) -> f64 {
        let major_r = self.major_radius;
        let minor_r = self.minor_radius;
        
        // Warping from PPF theory
        let involute_factor = 1.0 + coords.psi * 0.1 * self.involute_scale;
        let rho = major_r + minor_r * coords.theta.cos() * involute_factor;
        
        // PPF warping function
        let alpha = self.critical_ratio;
        let warp = (1.0 + alpha * rho / major_r).exp();
        
        warp
    }

    /// Compute connection coefficients (Christoffel symbols)
    pub fn christoffel_symbols(&self, coords: &IOTCoordinates) -> ChristoffelSymbols {
        let major_r = self.major_radius;
        let minor_r = self.minor_radius;
        let scale = self.involute_scale;
        
        let involute_factor = 1.0 + coords.psi * 0.1 * scale;
        let rho = major_r + minor_r * coords.theta.cos() * involute_factor;
        
        // Key non-zero Christoffel symbols for IOT metric
        let gamma_phi_phi_theta = -minor_r * coords.theta.sin() * involute_factor / rho;
        let gamma_phi_theta_phi = minor_r * coords.theta.sin() * involute_factor / (rho * rho);
        let gamma_theta_phi_phi = minor_r * coords.theta.sin() * involute_factor * rho / (minor_r * minor_r);
        let gamma_theta_theta_psi = 0.1 * scale;
        let gamma_psi_theta_theta = -0.1 * scale;
        
        ChristoffelSymbols {
            gamma_phi_phi_theta,
            gamma_phi_theta_phi,
            gamma_theta_phi_phi,
            gamma_theta_theta_psi,
            gamma_psi_theta_theta,
        }
    }

    /// Verify the critical ratio from variational principle
    pub fn verify_critical_ratio(&self) -> Result<bool, IOTError> {
        // The critical ratio r/R = 1/30 minimizes the action functional
        // ∫ (R + K) √g d³x where R is Ricci scalar, K is extrinsic curvature
        
        let test_coords = IOTCoordinates::new(PI / 4.0, PI / 4.0, 0.1)?;
        let ricci = self.ricci_scalar(&test_coords);
        let volume = self.volume_element(&test_coords);
        
        // Simplified verification: check if curvature is minimized
        let action_density = ricci * volume;
        
        // For critical ratio, action should be close to minimum
        let is_critical = action_density.abs() < 1.0 && self.is_critical();
        
        Ok(is_critical)
    }
}

/// Christoffel symbols for the IOT metric
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct ChristoffelSymbols {
    /// Γ^θ_φφ
    pub gamma_phi_phi_theta: f64,
    /// Γ^φ_θφ
    pub gamma_phi_theta_phi: f64,
    /// Γ^φ_θθ
    pub gamma_theta_phi_phi: f64,
    /// Γ^ψ_θθ
    pub gamma_theta_theta_psi: f64,
    /// Γ^θ_ψψ
    pub gamma_psi_theta_theta: f64,
}

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

impl fmt::Display for IOTCoordinates {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "(φ={:.3}, θ={:.3}, ψ={:.3})", self.phi, self.theta, self.psi)
    }
}

impl fmt::Display for IOTMetric {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        writeln!(f, "IOT Metric:")?;
        writeln!(f, "  Major radius R = {:.6}", self.major_radius)?;
        writeln!(f, "  Minor radius r = {:.6}", self.minor_radius)?;
        writeln!(f, "  Critical ratio r/R = {:.6}", self.critical_ratio)?;
        writeln!(f, "  Is critical: {}", self.is_critical())?;
        writeln!(f, "  Involute scale: {:.3}", self.involute_scale)?;
        
        if let Some(ref state_space) = self.state_space {
            writeln!(f, "  State space: S({})", state_space.value())?;
        }
        
        Ok(())
    }
}

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

    #[test]
    fn test_iot_coordinates() {
        let coords = IOTCoordinates::new(PI / 2.0, PI / 3.0, 0.5).unwrap();
        assert_eq!(coords.phi, PI / 2.0);
        assert_eq!(coords.theta, PI / 3.0);
        assert_eq!(coords.psi, 0.5);
        
        // Test invalid coordinates
        assert!(IOTCoordinates::new(-1.0, 0.0, 0.0).is_err());
        assert!(IOTCoordinates::new(3.0 * PI, 0.0, 0.0).is_err());
    }

    #[test]
    fn test_normalization() {
        let mut coords = IOTCoordinates::new(0.0, 0.0, 0.0).unwrap();
        coords.phi = 3.0 * PI;
        coords.theta = -PI / 2.0;
        coords.normalize();
        
        assert!(coords.phi >= 0.0 && coords.phi < 2.0 * PI);
        assert!(coords.theta >= 0.0 && coords.theta < 2.0 * PI);
    }

    #[test]
    fn test_cartesian_conversion() {
        let coords = IOTCoordinates::new(0.0, 0.0, 0.0).unwrap();
        let point = coords.to_cartesian(1.0, 0.5);
        
        assert_eq!(point.x, 1.5); // R + r = 1.0 + 0.5
        assert_eq!(point.y, 0.0);
        assert_eq!(point.z, 0.0);
    }

    #[test]
    fn test_critical_ratio() {
        let metric = IOTMetric::new();
        assert!(metric.is_critical());
        assert!((metric.critical_ratio() - IOTMetric::CRITICAL_RATIO).abs() < 1e-10);
    }

    #[test]
    fn test_metric_tensor() {
        let metric = IOTMetric::new();
        let coords = IOTCoordinates::new(0.0, 0.0, 0.0).unwrap();
        let tensor = metric.metric_tensor(&coords);
        
        // At θ=0, ψ=0, the torus has its maximum radius
        let expected_rho = metric.major_radius() + metric.minor_radius();
        assert!((tensor.g_phi_phi - expected_rho * expected_rho).abs() < 1e-10);
        assert!(tensor.g_theta_theta > 0.0);
        assert!(tensor.g_psi_psi >= 0.0);
    }

    #[test]
    fn test_ricci_scalar() {
        let metric = IOTMetric::new();
        let coords = IOTCoordinates::new(0.0, 0.0, 0.0).unwrap();
        let ricci = metric.ricci_scalar(&coords);
        
        // Should be finite and related to torus curvature
        assert!(ricci.is_finite());
        assert!(ricci > 0.0); // At θ=0, cos(θ)=1 gives positive curvature
    }

    #[test]
    fn test_geodesic_distance() {
        let metric = IOTMetric::new();
        let p1 = IOTCoordinates::new(0.0, 0.0, 0.0).unwrap();
        let p2 = IOTCoordinates::new(PI / 4.0, PI / 4.0, 0.1).unwrap();
        
        let distance = metric.geodesic_distance(&p1, &p2);
        assert!(distance > 0.0);
        assert!(distance.is_finite());
        
        // Distance to self should be zero
        let self_distance = metric.geodesic_distance(&p1, &p1);
        assert!(self_distance.abs() < 1e-10);
    }

    #[test]
    fn test_volume_element() {
        let metric = IOTMetric::new();
        let coords = IOTCoordinates::new(PI / 4.0, PI / 4.0, 0.1).unwrap();
        let volume = metric.volume_element(&coords);
        
        assert!(volume > 0.0);
        assert!(volume.is_finite());
    }

    #[test]
    fn test_from_state_space() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let metric = IOTMetric::from_state_space(state_space);
        
        assert!(metric.is_critical());
        assert!(metric.state_space.is_some());
        assert!(metric.involute_scale > 1.0); // Should be scaled by complexity
    }

    #[test]
    fn test_factorization_mapping() {
        let state_space = FactorizationStateSpace::new(6).unwrap();
        let metric = IOTMetric::from_state_space(state_space.clone());
        
        let factorization = &state_space.factorizations()[0];
        let coords = metric.factorization_to_coordinates(factorization);
        
        assert!(coords.phi >= 0.0 && coords.phi < 2.0 * PI);
        assert!(coords.theta >= 0.0 && coords.theta < 2.0 * PI);
        assert!(coords.psi >= 0.0);
    }

    #[test]
    fn test_warping_function() {
        let metric = IOTMetric::new();
        let coords = IOTCoordinates::new(0.0, 0.0, 0.0).unwrap();
        let warp = metric.warping_function(&coords);
        
        assert!(warp > 0.0);
        assert!(warp.is_finite());
        assert!(warp >= 1.0); // Exponential warping should be ≥ 1
    }

    #[test]
    fn test_christoffel_symbols() {
        let metric = IOTMetric::new();
        let coords = IOTCoordinates::new(PI / 4.0, PI / 4.0, 0.1).unwrap();
        let christoffel = metric.christoffel_symbols(&coords);
        
        // All symbols should be finite
        assert!(christoffel.gamma_phi_phi_theta.is_finite());
        assert!(christoffel.gamma_phi_theta_phi.is_finite());
        assert!(christoffel.gamma_theta_phi_phi.is_finite());
        assert!(christoffel.gamma_theta_theta_psi.is_finite());
        assert!(christoffel.gamma_psi_theta_theta.is_finite());
    }

    #[test]
    fn test_critical_ratio_verification() {
        let metric = IOTMetric::new();
        let is_critical = metric.verify_critical_ratio().unwrap();
        assert!(is_critical);
        
        // Test with non-critical ratio
        let non_critical = IOTMetric::with_radii(1.0, 0.1).unwrap();
        assert!(!non_critical.is_critical());
    }

    #[test]
    fn test_display() {
        let coords = IOTCoordinates::new(1.0, 2.0, 0.5).unwrap();
        let coords_str = format!("{}", coords);
        assert!(coords_str.contains("φ=1.000"));
        assert!(coords_str.contains("θ=2.000"));
        assert!(coords_str.contains("ψ=0.500"));
        
        let metric = IOTMetric::new();
        let metric_str = format!("{}", metric);
        assert!(metric_str.contains("Critical ratio"));
        assert!(metric_str.contains("Is critical: true"));
    }
}