persistence_agent/
vietoris_rips.rs

1use crate::error::PersistenceError;
2use crate::point_cloud::PointCloud;
3
4/// A Vietoris-Rips complex built from a point cloud.
5/// Simplices are added as epsilon grows — a k-simplex [v0,…,vk] appears at the
6/// maximum pairwise distance among its vertices.
7#[derive(Debug, Clone)]
8pub struct VietorisRipsComplex {
9    pub simplices: Vec<Vec<usize>>,
10    pub filtration_values: Vec<f64>,
11    pub max_dimension: usize,
12}
13
14impl VietorisRipsComplex {
15    /// Build the Vietoris-Rips complex up to `max_dimension`.
16    /// For large point clouds you can pass `max_eps` to prune simplices whose
17    /// filtration value exceeds it; pass `f64::INFINITY` to include everything.
18    pub fn build(
19        cloud: &PointCloud,
20        max_dimension: usize,
21        max_eps: f64,
22    ) -> Result<Self, PersistenceError> {
23        let n = cloud.n_points();
24        if n == 0 {
25            return Err(PersistenceError::EmptyCloud);
26        }
27
28        let mut simplices: Vec<Vec<usize>> = Vec::new();
29        let mut filtration_values: Vec<f64> = Vec::new();
30
31        // Dimension 0: each vertex appears at epsilon = 0
32        for i in 0..n {
33            simplices.push(vec![i]);
34            filtration_values.push(0.0);
35        }
36
37        // Iteratively build higher-dimensional simplices
38        let mut current_simplices: Vec<Vec<usize>> = (0..n).map(|i| vec![i]).collect();
39
40        for _dim in 0..max_dimension {
41            let mut candidates: Vec<(f64, Vec<usize>)> = Vec::new();
42            for simplex in &current_simplices {
43                // Try to extend by adding any vertex with index > max(simplex)
44                let max_v = *simplex.iter().max().unwrap();
45                for v in (max_v + 1)..n {
46                    let mut new_simplex = simplex.clone();
47                    new_simplex.push(v);
48                    new_simplex.sort();
49
50                    // Filtration value = max pairwise distance
51                    let filt = max_pairwise(cloud, &new_simplex);
52                    if filt <= max_eps {
53                        candidates.push((filt, new_simplex));
54                    }
55                }
56            }
57
58            // Sort by filtration value, break ties by lexicographic order
59            candidates.sort_by(|a, b| {
60                a.0.partial_cmp(&b.0)
61                    .unwrap_or(std::cmp::Ordering::Equal)
62                    .then_with(|| a.1.cmp(&b.1))
63            });
64
65            current_simplices = candidates.iter().map(|(_, s)| s.clone()).collect();
66            for (filt, s) in &candidates {
67                simplices.push(s.clone());
68                filtration_values.push(*filt);
69            }
70
71            if current_simplices.is_empty() {
72                break;
73            }
74        }
75
76        Ok(Self {
77            simplices,
78            filtration_values,
79            max_dimension,
80        })
81    }
82
83    pub fn n_simplices(&self) -> usize {
84        self.simplices.len()
85    }
86
87    /// Return the dimension of a simplex at given index.
88    pub fn simplex_dimension(&self, idx: usize) -> usize {
89        self.simplices[idx].len().saturating_sub(1)
90    }
91
92    /// Indices of all simplices of a given dimension.
93    pub fn simplices_of_dimension(&self, dim: usize) -> Vec<usize> {
94        self.simplices
95            .iter()
96            .enumerate()
97            .filter(|(_, s)| s.len().saturating_sub(1) == dim)
98            .map(|(i, _)| i)
99            .collect()
100    }
101
102    /// All (dimension, simplex-index) pairs sorted by filtration value then dimension.
103    pub fn sorted_filtration(&self) -> Vec<(usize, usize)> {
104        let mut indexed: Vec<(usize, usize, f64)> = self
105            .simplices
106            .iter()
107            .enumerate()
108            .map(|(i, s)| (s.len().saturating_sub(1), i, self.filtration_values[i]))
109            .collect();
110        indexed.sort_by(|a, b| {
111            a.2.partial_cmp(&b.2)
112                .unwrap_or(std::cmp::Ordering::Equal)
113                .then_with(|| a.0.cmp(&b.0))
114        });
115        indexed.into_iter().map(|(d, i, _)| (d, i)).collect()
116    }
117}
118
119fn max_pairwise(cloud: &PointCloud, simplex: &[usize]) -> f64 {
120    let mut mx = 0.0_f64;
121    for i in 0..simplex.len() {
122        for j in (i + 1)..simplex.len() {
123            mx = mx.max(cloud.distance(simplex[i], simplex[j]));
124        }
125    }
126    mx
127}
persistence_agent/vietoris_rips.rs

persistence_agent/
vietoris_rips.rs
