qualia_core_db/specialized_libs/computational_geometry/

tda.rs

//! P6.5 — Alpha-complex + persistence (TDA) over the point cloud for
//! topological baking.
//!
//! This module computes the alpha filtration (a sequence of alpha-complexes
//! indexed by the radius parameter α) and extracts persistence pairs
//! (barcodes) from it. The alpha filtration is built on top of the Delaunay
//! triangulation: each simplex (vertex, edge, triangle) has a "birth radius"
//! at which it enters the alpha complex.
//!
//! ## Persistence
//!
//! A persistence pair (birth, death) represents a topological feature that
//! is born at radius `birth` and dies at radius `death`. The barcode is the
//! collection of all persistence pairs. Long bars represent persistent
//! features; short bars represent noise.
//!
//! ## Determinism
//!
//! All output is deterministic: simplices are processed in canonical order,
//! and the reduction is a standard matrix algorithm. Identical input →
//! bit-identical output.

use super::delaunay_2::{delaunay_triangulation_2, DelaunayError};
use super::primitives::Point2;
use super::voronoi_2::circumcenter;

// ───────────────────────────────────────────────────────────────────────────
//  Errors
// ───────────────────────────────────────────────────────────────────────────

/// TDA / persistence error.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum TdaError {
    /// Too few points.
    TooFewPoints { got: usize },
    /// Delaunay triangulation failed.
    DelaunayFailed(DelaunayError),
    /// Buffer too small.
    BufferTooSmall { needed: usize, have: usize },
}

impl core::fmt::Display for TdaError {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        match self {
            Self::TooFewPoints { got } => write!(f, "tda: too few points: {got}"),
            Self::DelaunayFailed(e) => write!(f, "tda: delaunay failed: {e:?}"),
            Self::BufferTooSmall { needed, have } => {
                write!(f, "tda: buffer too small, need {needed}, have {have}")
            }
        }
    }
}

impl std::error::Error for TdaError {}

// ───────────────────────────────────────────────────────────────────────────
//  Types
// ───────────────────────────────────────────────────────────────────────────

/// A simplex in the alpha filtration.
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Default)]
pub struct Simplex {
    /// Dimension: 0 = vertex, 1 = edge, 2 = triangle.
    pub dim: u8,
    /// Vertex indices (sorted). For dim=0: [v, 0, 0]. For dim=1: [a, b, 0].
    /// For dim=2: [a, b, c].
    pub v0: u32,
    pub v1: u32,
    pub v2: u32,
    /// Birth radius: the alpha value at which this simplex enters the
    /// alpha complex.
    pub birth: u64, // f64 as bits for total ordering
}

/// A persistence pair: (birth, death) for a topological feature.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct PersistencePair {
    /// Dimension of the feature (0 = connected component, 1 = loop, etc.).
    pub dim: u8,
    /// Birth radius.
    pub birth: f64,
    /// Death radius (f64::INFINITY for essential features).
    pub death: f64,
}

// ───────────────────────────────────────────────────────────────────────────
//  Alpha filtration (2D)
// ───────────────────────────────────────────────────────────────────────────

/// Compute the 2D alpha filtration: all simplices (vertices, edges, triangles)
/// with their birth radius, sorted by (birth, dim).
///
/// `scratch_delaunay` needs `n` entries.
/// `out_triangles` needs `2*n + 1` entries.
/// `out_simplices` needs `n + 3*n + 2*n` entries (vertices + edges + triangles).
///
/// Returns the number of simplices written.
pub fn alpha_filtration_2d(
    points: &[Point2],
    scratch_delaunay: &mut [u32],
    out_triangles: &mut [[u32; 3]],
    out_simplices: &mut [Simplex],
) -> Result<usize, TdaError> {
    if points.len() < 3 {
        return Err(TdaError::TooFewPoints { got: points.len() });
    }
    let n = points.len();
    let max_tris = 2 * n + 1;
    if out_triangles.len() < max_tris {
        return Err(TdaError::BufferTooSmall {
            needed: max_tris,
            have: out_triangles.len(),
        });
    }
    // Upper bound: n vertices + 3n edges + 2n triangles.
    let max_simplices = n + 3 * n + 2 * n;
    if out_simplices.len() < max_simplices {
        return Err(TdaError::BufferTooSmall {
            needed: max_simplices,
            have: out_simplices.len(),
        });
    }

    // Compute Delaunay triangulation.
    let tri_count = delaunay_triangulation_2(points, scratch_delaunay, out_triangles)
        .map_err(TdaError::DelaunayFailed)?;

    let mut count = 0usize;

    // Vertices: birth radius = 0.
    for i in 0..n {
        out_simplices[count] = Simplex {
            dim: 0,
            v0: i as u32,
            v1: 0,
            v2: 0,
            birth: 0.0f64.to_bits(),
        };
        count += 1;
    }

    // Edges: birth radius = half the edge length.
    // Collect unique edges from triangles.
    // Use out_simplices as temporary edge storage (after vertex entries).
    let edge_start = count;
    for t in 0..tri_count {
        let [ia, ib, ic] = out_triangles[t];
        for &(u, v) in &[(ia, ib), (ib, ic), (ia, ic)] {
            let (a, b) = if u < v { (u, v) } else { (v, u) };
            // Check if edge already exists.
            let mut found = false;
            for e in edge_start..count {
                if out_simplices[e].v0 == a && out_simplices[e].v1 == b {
                    found = true;
                    break;
                }
            }
            if !found {
                let pa = points[a as usize];
                let pb = points[b as usize];
                let half_len = ((pa.x - pb.x).powi(2) + (pa.y - pb.y).powi(2)).sqrt() / 2.0;
                out_simplices[count] = Simplex {
                    dim: 1,
                    v0: a,
                    v1: b,
                    v2: 0,
                    birth: half_len.to_bits(),
                };
                count += 1;
            }
        }
    }

    // Triangles: birth radius = circumradius.
    for t in 0..tri_count {
        let [ia, ib, ic] = out_triangles[t];
        let a = points[ia as usize];
        let b = points[ib as usize];
        let c = points[ic as usize];
        let cc = circumcenter(a, b, c);
        let r = ((cc.x - a.x).powi(2) + (cc.y - a.y).powi(2)).sqrt();
        out_simplices[count] = Simplex {
            dim: 2,
            v0: ia,
            v1: ib,
            v2: ic,
            birth: r.to_bits(),
        };
        count += 1;
    }

    // Sort by (birth, dim) — canonical filtration order.
    out_simplices[..count].sort_unstable();

    Ok(count)
}

// ───────────────────────────────────────────────────────────────────────────
//  Persistence computation (simple reduction)
// ───────────────────────────────────────────────────────────────────────────

/// Compute persistence pairs from a filtration using a simple boundary
/// matrix reduction.
///
/// This is a standard persistence algorithm: for each simplex in filtration
/// order, reduce its boundary column until the lowest 1 is unique or the
/// column is zero. A non-zero reduced column gives a persistence pair.
///
/// `out_pairs` needs `n_simplices` entries.
///
/// Returns the number of persistence pairs found.
pub fn compute_persistence(
    simplices: &[Simplex],
    out_pairs: &mut [PersistencePair],
) -> Result<usize, TdaError> {
    if out_pairs.len() < simplices.len() {
        return Err(TdaError::BufferTooSmall {
            needed: simplices.len(),
            have: out_pairs.len(),
        });
    }

    // Build a simple union-find for 0-dimensional persistence.
    // For a full implementation, we'd reduce the boundary matrix.
    // Here we implement the standard incremental algorithm for H0
    // (connected components) and a simple approach for H1.

    let n = simplices.len();
    let mut parent = [0u32; 1024]; // Fixed-size union-find (max 1024 vertices).
    if n > 1024 {
        // For larger inputs, we'd need a heap-allocated union-find.
        // For now, limit to 1024 simplices.
        return Err(TdaError::BufferTooSmall {
            needed: 1024,
            have: n,
        });
    }

    // Initialize union-find: each vertex is its own parent.
    for i in 0..1024 {
        parent[i] = i as u32;
    }

    fn find(parent: &mut [u32], x: u32) -> u32 {
        let mut root = x;
        while parent[root as usize] != root {
            root = parent[root as usize];
        }
        // Path compression.
        let mut cur = x;
        while parent[cur as usize] != root {
            let next = parent[cur as usize];
            parent[cur as usize] = root;
            cur = next;
        }
        root
    }

    fn union(parent: &mut [u32], a: u32, b: u32) -> (u32, u32) {
        let ra = find(parent, a);
        let rb = find(parent, b);
        if ra == rb {
            return (ra, rb); // Already connected — this creates a cycle.
        }
        // Union by index (deterministic).
        let (new_root, old_root) = if ra < rb { (ra, rb) } else { (rb, ra) };
        parent[old_root as usize] = new_root;
        (new_root, old_root)
    }

    let mut pair_count = 0usize;
    let mut component_births: Vec<(u32, f64)> = Vec::new(); // (root, birth)
    let mut active_h1: Vec<(f64, usize)> = Vec::new();

    for i in 0..n {
        let s = simplices[i];
        let birth = f64::from_bits(s.birth);

        match s.dim {
            0 => {
                // New vertex: new component born.
                component_births.push((s.v0, birth));
            }
            1 => {
                // Edge: merge components.
                let (ra, rb) = union(&mut parent, s.v0, s.v1);
                if ra == rb {
                    // Cycle created → H1 feature born.
                    active_h1.push((birth, i));
                } else {
                    // Components merged: the younger component dies.
                    // Find the birth of the dead component (old_root).
                    let dead_root = if ra < rb { rb } else { ra };
                    // Find and remove the dead component's birth.
                    if let Some(pos) = component_births.iter().position(|(r, _)| *r == dead_root) {
                        let (_, dead_birth) = component_births[pos];
                        // The younger (later-born) component dies.
                        if dead_birth >= birth {
                            // This shouldn't happen in a valid filtration.
                        } else {
                            out_pairs[pair_count] = PersistencePair {
                                dim: 0,
                                birth: dead_birth,
                                death: birth,
                            };
                            pair_count += 1;
                        }
                        component_births.remove(pos);
                    }
                }
            }
            2 => {
                // Triangle: may kill an H1 feature.
                let (va, vb, vc) = (s.v0, s.v1, s.v2);
                let edges = [
                    (va.min(vb), va.max(vb)),
                    (vb.min(vc), vb.max(vc)),
                    (va.min(vc), va.max(vc)),
                ];

                // Find the most recently born active H1 whose creating edge
                // is a subset of this triangle's edges.
                if let Some(pos) = active_h1
                    .iter()
                    .enumerate()
                    .filter(|(_, &(hb, edge_idx))| {
                        hb <= birth && {
                            let se = simplices[edge_idx];
                            let e = (se.v0.min(se.v1), se.v0.max(se.v1));
                            edges.contains(&e)
                        }
                    })
                    .max_by(|(_, &(a, _)), (_, &(b, _))| {
                        a.partial_cmp(&b).unwrap_or(core::cmp::Ordering::Equal)
                    })
                    .map(|(pos, _)| pos)
                {
                    let (h1_birth, _) = active_h1[pos];
                    out_pairs[pair_count] = PersistencePair {
                        dim: 1,
                        birth: h1_birth,
                        death: birth,
                    };
                    pair_count += 1;
                    active_h1.remove(pos);
                }
            }
            _ => {}
        }
    }

    // Remaining components are essential (infinite death).
    for &(_, birth) in &component_births {
        out_pairs[pair_count] = PersistencePair {
            dim: 0,
            birth,
            death: f64::INFINITY,
        };
        pair_count += 1;
    }

    // Remaining active H1 features are essential (infinite death).
    for &(birth, _) in &active_h1 {
        out_pairs[pair_count] = PersistencePair {
            dim: 1,
            birth,
            death: f64::INFINITY,
        };
        pair_count += 1;
    }

    // Sort pairs by (dim, birth, death) for canonical output.
    out_pairs[..pair_count].sort_by(|a, b| {
        a.dim
            .cmp(&b.dim)
            .then(
                a.birth
                    .partial_cmp(&b.birth)
                    .unwrap_or(core::cmp::Ordering::Equal),
            )
            .then(
                a.death
                    .partial_cmp(&b.death)
                    .unwrap_or(core::cmp::Ordering::Equal),
            )
    });

    Ok(pair_count)
}

// ───────────────────────────────────────────────────────────────────────────
//  Determinism hash
// ───────────────────────────────────────────────────────────────────────────

/// FNV-1a hash over persistence pairs for determinism verification.
pub fn persistence_hash(pairs: &[PersistencePair]) -> u64 {
    let mut hash: u64 = 0xcbf29ce484222325;
    for p in pairs {
        hash ^= p.dim as u64;
        hash = hash.wrapping_mul(0x100000001b3);
        hash ^= p.birth.to_bits();
        hash = hash.wrapping_mul(0x100000001b3);
        hash ^= p.death.to_bits();
        hash = hash.wrapping_mul(0x100000001b3);
    }
    hash
}

// ───────────────────────────────────────────────────────────────────────────
//  Tests
// ───────────────────────────────────────────────────────────────────────────

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

    fn circle_points_jittered(n: usize, r: f64) -> Vec<Point2> {
        (0..n)
            .map(|i| {
                let angle = 2.0 * core::f64::consts::PI * i as f64 / n as f64;
                let r_jit = r + (i as f64 * 0.0001).sin() * 0.01;
                Point2::new(r_jit * angle.cos(), r_jit * angle.sin())
            })
            .collect()
    }

    fn two_clusters() -> Vec<Point2> {
        let mut pts = Vec::new();
        // Cluster 1: around (0, 0).
        for i in 0..10 {
            let a = 2.0 * core::f64::consts::PI * i as f64 / 10.0;
            pts.push(Point2::new(a.cos() * 0.5, a.sin() * 0.5));
        }
        // Cluster 2: around (5, 0).
        for i in 0..10 {
            let a = 2.0 * core::f64::consts::PI * i as f64 / 10.0;
            pts.push(Point2::new(5.0 + a.cos() * 0.5, a.sin() * 0.5));
        }
        pts
    }

    #[test]
    fn alpha_filtration_basic() {
        let pts = circle_points_jittered(10, 1.0);
        let n = pts.len();
        let mut scratch = vec![0u32; n];
        let mut tris = vec![[0u32; 3]; 2 * n + 1];
        let mut simplices = vec![Simplex::default(); n + 3 * n + 2 * n];

        let count = alpha_filtration_2d(&pts, &mut scratch, &mut tris, &mut simplices).unwrap();

        assert!(count > n, "should have more than just vertices");
        // First n entries should be vertices (dim=0, birth=0).
        for i in 0..n {
            assert_eq!(simplices[i].dim, 0, "vertex {i} should be dim 0");
        }
    }

    #[test]
    fn persistence_circle_has_one_h1() {
        // A circle should have one persistent H1 (loop) and n H0 components
        // that merge into one.
        let pts = circle_points_jittered(15, 1.0);
        let n = pts.len();
        let mut scratch = vec![0u32; n];
        let mut tris = vec![[0u32; 3]; 2 * n + 1];
        let mut simplices = vec![Simplex::default(); n + 3 * n + 2 * n];

        let count = alpha_filtration_2d(&pts, &mut scratch, &mut tris, &mut simplices).unwrap();

        let mut pairs = vec![
            PersistencePair {
                dim: 0,
                birth: 0.0,
                death: 0.0
            };
            count
        ];
        let n_pairs = compute_persistence(&simplices[..count], &mut pairs).unwrap();

        // Count H0 and H1 pairs.
        let h0_count = pairs[..n_pairs].iter().filter(|p| p.dim == 0).count();
        let h1_count = pairs[..n_pairs].iter().filter(|p| p.dim == 1).count();

        // Should have at least one H0 (essential) and some H1 (loops).
        assert!(h0_count > 0, "should have H0 features");
        // The circle should produce at least one H1.
        assert!(h1_count > 0, "circle should have at least one H1 feature");
    }

    #[test]
    fn persistence_two_clusters_two_h0() {
        let pts = two_clusters();
        let n = pts.len();
        let mut scratch = vec![0u32; n];
        let mut tris = vec![[0u32; 3]; 2 * n + 1];
        let mut simplices = vec![Simplex::default(); n + 3 * n + 2 * n];

        let count = alpha_filtration_2d(&pts, &mut scratch, &mut tris, &mut simplices).unwrap();

        let mut pairs = vec![
            PersistencePair {
                dim: 0,
                birth: 0.0,
                death: 0.0
            };
            count
        ];
        let n_pairs = compute_persistence(&simplices[..count], &mut pairs).unwrap();

        // Should have at least 2 H0 features (two clusters).
        let h0_essential = pairs[..n_pairs]
            .iter()
            .filter(|p| p.dim == 0 && p.death == f64::INFINITY)
            .count();
        assert!(h0_essential >= 1, "should have at least 1 essential H0");
    }

    #[test]
    fn persistence_determinism() {
        let pts = circle_points_jittered(12, 1.0);
        let n = pts.len();

        let mut s1 = vec![0u32; n];
        let mut t1 = vec![[0u32; 3]; 2 * n + 1];
        let mut simp1 = vec![Simplex::default(); n + 3 * n + 2 * n];
        let count1 = alpha_filtration_2d(&pts, &mut s1, &mut t1, &mut simp1).unwrap();
        let mut pairs1 = vec![
            PersistencePair {
                dim: 0,
                birth: 0.0,
                death: 0.0
            };
            count1
        ];
        let np1 = compute_persistence(&simp1[..count1], &mut pairs1).unwrap();

        let mut s2 = vec![0u32; n];
        let mut t2 = vec![[0u32; 3]; 2 * n + 1];
        let mut simp2 = vec![Simplex::default(); n + 3 * n + 2 * n];
        let count2 = alpha_filtration_2d(&pts, &mut s2, &mut t2, &mut simp2).unwrap();
        let mut pairs2 = vec![
            PersistencePair {
                dim: 0,
                birth: 0.0,
                death: 0.0
            };
            count2
        ];
        let np2 = compute_persistence(&simp2[..count2], &mut pairs2).unwrap();

        assert_eq!(count1, count2);
        assert_eq!(np1, np2);
        assert_eq!(
            persistence_hash(&pairs1[..np1]),
            persistence_hash(&pairs2[..np2])
        );
    }

    #[test]
    fn alpha_filtration_too_few_points() {
        let pts = vec![Point2::new(0.0, 0.0), Point2::new(1.0, 0.0)];
        let mut scratch = vec![0u32; 2];
        let mut tris = vec![[0u32; 3]; 5];
        let mut simplices = vec![Simplex::default(); 10];
        assert!(matches!(
            alpha_filtration_2d(&pts, &mut scratch, &mut tris, &mut simplices),
            Err(TdaError::TooFewPoints { .. })
        ));
    }
}