ruvector_math/homology/
persistence.rs

1//! Persistent Homology Computation
2//!
3//! Compute birth-death pairs from a filtration using the standard algorithm.
4
5use super::{BettiNumbers, Filtration, Simplex};
6use std::collections::{HashMap, HashSet};
7
8/// Birth-death pair in persistence diagram
9#[derive(Debug, Clone, PartialEq)]
10pub struct BirthDeathPair {
11    /// Dimension of the feature (0 = component, 1 = loop, ...)
12    pub dimension: usize,
13    /// Birth time
14    pub birth: f64,
15    /// Death time (None = essential, lives forever)
16    pub death: Option<f64>,
17}
18
19impl BirthDeathPair {
20    /// Create finite interval
21    pub fn finite(dimension: usize, birth: f64, death: f64) -> Self {
22        Self {
23            dimension,
24            birth,
25            death: Some(death),
26        }
27    }
28
29    /// Create essential (infinite) interval
30    pub fn essential(dimension: usize, birth: f64) -> Self {
31        Self {
32            dimension,
33            birth,
34            death: None,
35        }
36    }
37
38    /// Persistence (lifetime) of feature
39    pub fn persistence(&self) -> f64 {
40        match self.death {
41            Some(d) => d - self.birth,
42            None => f64::INFINITY,
43        }
44    }
45
46    /// Is this an essential feature (never dies)?
47    pub fn is_essential(&self) -> bool {
48        self.death.is_none()
49    }
50
51    /// Midpoint of interval
52    pub fn midpoint(&self) -> f64 {
53        match self.death {
54            Some(d) => (self.birth + d) / 2.0,
55            None => f64::INFINITY,
56        }
57    }
58}
59
60/// Persistence diagram: collection of birth-death pairs
61#[derive(Debug, Clone)]
62pub struct PersistenceDiagram {
63    /// Birth-death pairs
64    pub pairs: Vec<BirthDeathPair>,
65    /// Maximum dimension
66    pub max_dim: usize,
67}
68
69impl PersistenceDiagram {
70    /// Create empty diagram
71    pub fn new() -> Self {
72        Self {
73            pairs: Vec::new(),
74            max_dim: 0,
75        }
76    }
77
78    /// Add a pair
79    pub fn add(&mut self, pair: BirthDeathPair) {
80        self.max_dim = self.max_dim.max(pair.dimension);
81        self.pairs.push(pair);
82    }
83
84    /// Get pairs of dimension d
85    pub fn pairs_of_dim(&self, d: usize) -> impl Iterator<Item = &BirthDeathPair> {
86        self.pairs.iter().filter(move |p| p.dimension == d)
87    }
88
89    /// Get Betti numbers at scale t
90    pub fn betti_at(&self, t: f64) -> BettiNumbers {
91        let mut b0 = 0;
92        let mut b1 = 0;
93        let mut b2 = 0;
94
95        for pair in &self.pairs {
96            let alive = pair.birth <= t && pair.death.map(|d| d > t).unwrap_or(true);
97            if alive {
98                match pair.dimension {
99                    0 => b0 += 1,
100                    1 => b1 += 1,
101                    2 => b2 += 1,
102                    _ => {}
103                }
104            }
105        }
106
107        BettiNumbers::new(b0, b1, b2)
108    }
109
110    /// Get total persistence (sum of lifetimes)
111    pub fn total_persistence(&self) -> f64 {
112        self.pairs
113            .iter()
114            .filter(|p| !p.is_essential())
115            .map(|p| p.persistence())
116            .sum()
117    }
118
119    /// Get average persistence
120    pub fn average_persistence(&self) -> f64 {
121        let finite: Vec<f64> = self
122            .pairs
123            .iter()
124            .filter(|p| !p.is_essential())
125            .map(|p| p.persistence())
126            .collect();
127
128        if finite.is_empty() {
129            0.0
130        } else {
131            finite.iter().sum::<f64>() / finite.len() as f64
132        }
133    }
134
135    /// Filter by minimum persistence
136    pub fn filter_by_persistence(&self, min_persistence: f64) -> Self {
137        Self {
138            pairs: self
139                .pairs
140                .iter()
141                .filter(|p| p.persistence() >= min_persistence)
142                .cloned()
143                .collect(),
144            max_dim: self.max_dim,
145        }
146    }
147
148    /// Number of features of each dimension
149    pub fn feature_counts(&self) -> Vec<usize> {
150        let mut counts = vec![0; self.max_dim + 1];
151        for pair in &self.pairs {
152            if pair.dimension <= self.max_dim {
153                counts[pair.dimension] += 1;
154            }
155        }
156        counts
157    }
158}
159
160impl Default for PersistenceDiagram {
161    fn default() -> Self {
162        Self::new()
163    }
164}
165
166/// Persistent homology computation
167pub struct PersistentHomology {
168    /// Working column representation (reduced boundary matrix)
169    columns: Vec<Option<HashSet<usize>>>,
170    /// Pivot to column mapping
171    pivot_to_col: HashMap<usize, usize>,
172    /// Birth times
173    birth_times: Vec<f64>,
174    /// Simplex dimensions
175    dimensions: Vec<usize>,
176}
177
178impl PersistentHomology {
179    /// Compute persistence from filtration
180    pub fn compute(filtration: &Filtration) -> PersistenceDiagram {
181        let mut ph = Self {
182            columns: Vec::new(),
183            pivot_to_col: HashMap::new(),
184            birth_times: Vec::new(),
185            dimensions: Vec::new(),
186        };
187
188        ph.run(filtration)
189    }
190
191    fn run(&mut self, filtration: &Filtration) -> PersistenceDiagram {
192        let n = filtration.simplices.len();
193        if n == 0 {
194            return PersistenceDiagram::new();
195        }
196
197        // Build simplex index mapping
198        let simplex_to_idx: HashMap<&Simplex, usize> = filtration
199            .simplices
200            .iter()
201            .enumerate()
202            .map(|(i, fs)| (&fs.simplex, i))
203            .collect();
204
205        // Initialize
206        self.columns = Vec::with_capacity(n);
207        self.birth_times = filtration.simplices.iter().map(|fs| fs.birth).collect();
208        self.dimensions = filtration
209            .simplices
210            .iter()
211            .map(|fs| fs.simplex.dim())
212            .collect();
213
214        // Build boundary matrix columns
215        for fs in &filtration.simplices {
216            let boundary = self.boundary(&fs.simplex, &simplex_to_idx);
217            self.columns.push(if boundary.is_empty() {
218                None
219            } else {
220                Some(boundary)
221            });
222        }
223
224        // Reduce matrix
225        self.reduce();
226
227        // Extract persistence pairs
228        self.extract_pairs()
229    }
230
231    /// Compute boundary of simplex as set of face indices
232    fn boundary(&self, simplex: &Simplex, idx_map: &HashMap<&Simplex, usize>) -> HashSet<usize> {
233        let mut boundary = HashSet::new();
234        for face in simplex.faces() {
235            if let Some(&idx) = idx_map.get(&face) {
236                boundary.insert(idx);
237            }
238        }
239        boundary
240    }
241
242    /// Reduce using standard persistence algorithm
243    fn reduce(&mut self) {
244        let n = self.columns.len();
245
246        for j in 0..n {
247            // Reduce column j
248            while let Some(pivot) = self.get_pivot(j) {
249                if let Some(&other) = self.pivot_to_col.get(&pivot) {
250                    // Add column 'other' to column j (mod 2)
251                    self.add_columns(j, other);
252                } else {
253                    // No collision, record pivot
254                    self.pivot_to_col.insert(pivot, j);
255                    break;
256                }
257            }
258        }
259    }
260
261    /// Get pivot (largest index) of column
262    fn get_pivot(&self, col: usize) -> Option<usize> {
263        self.columns[col]
264            .as_ref()
265            .and_then(|c| c.iter().max().copied())
266    }
267
268    /// Add column src to column dst (XOR / mod 2)
269    fn add_columns(&mut self, dst: usize, src: usize) {
270        let src_col = self.columns[src].clone();
271        if let (Some(ref mut dst_col), Some(ref src_col)) = (&mut self.columns[dst], &src_col) {
272            // Symmetric difference
273            let mut new_col = HashSet::new();
274            for &idx in dst_col.iter() {
275                if !src_col.contains(&idx) {
276                    new_col.insert(idx);
277                }
278            }
279            for &idx in src_col.iter() {
280                if !dst_col.contains(&idx) {
281                    new_col.insert(idx);
282                }
283            }
284            if new_col.is_empty() {
285                self.columns[dst] = None;
286            } else {
287                *dst_col = new_col;
288            }
289        }
290    }
291
292    /// Extract birth-death pairs from reduced matrix
293    fn extract_pairs(&self) -> PersistenceDiagram {
294        let n = self.columns.len();
295        let mut diagram = PersistenceDiagram::new();
296        let mut paired = HashSet::new();
297
298        // Process pivot pairs (death creates pair with birth)
299        for (&pivot, &col) in &self.pivot_to_col {
300            let birth = self.birth_times[pivot];
301            let death = self.birth_times[col];
302            let dim = self.dimensions[pivot];
303
304            if death > birth {
305                diagram.add(BirthDeathPair::finite(dim, birth, death));
306            }
307
308            paired.insert(pivot);
309            paired.insert(col);
310        }
311
312        // Remaining columns are essential (infinite persistence)
313        for j in 0..n {
314            if !paired.contains(&j) && self.columns[j].is_none() {
315                let dim = self.dimensions[j];
316                let birth = self.birth_times[j];
317                diagram.add(BirthDeathPair::essential(dim, birth));
318            }
319        }
320
321        diagram
322    }
323}
324
325#[cfg(test)]
326mod tests {
327    use super::*;
328    use crate::homology::{PointCloud, VietorisRips};
329
330    #[test]
331    fn test_birth_death_pair() {
332        let finite = BirthDeathPair::finite(0, 0.0, 1.0);
333        assert_eq!(finite.persistence(), 1.0);
334        assert!(!finite.is_essential());
335
336        let essential = BirthDeathPair::essential(0, 0.0);
337        assert!(essential.is_essential());
338        assert_eq!(essential.persistence(), f64::INFINITY);
339    }
340
341    #[test]
342    fn test_persistence_diagram() {
343        let mut diagram = PersistenceDiagram::new();
344        diagram.add(BirthDeathPair::essential(0, 0.0));
345        diagram.add(BirthDeathPair::finite(0, 0.0, 1.0));
346        diagram.add(BirthDeathPair::finite(1, 0.5, 2.0));
347
348        assert_eq!(diagram.pairs.len(), 3);
349
350        let betti = diagram.betti_at(0.75);
351        assert_eq!(betti.b0, 2); // Both 0-dim features alive
352        assert_eq!(betti.b1, 1); // 1-dim feature alive
353    }
354
355    #[test]
356    fn test_persistent_homology_simple() {
357        // Two points
358        let cloud = PointCloud::from_flat(&[0.0, 0.0, 1.0, 0.0], 2);
359        let rips = VietorisRips::new(1, 2.0);
360        let filtration = rips.build(&cloud);
361
362        let diagram = PersistentHomology::compute(&filtration);
363
364        // Should have:
365        // - One essential H0 (final connected component)
366        // - One finite H0 that dies when edge connects the points
367        let h0_pairs: Vec<_> = diagram.pairs_of_dim(0).collect();
368        assert!(!h0_pairs.is_empty());
369    }
370
371    #[test]
372    fn test_persistent_homology_triangle() {
373        // Three points forming triangle
374        let cloud = PointCloud::from_flat(&[0.0, 0.0, 1.0, 0.0, 0.5, 0.866], 2);
375        let rips = VietorisRips::new(2, 2.0);
376        let filtration = rips.build(&cloud);
377
378        let diagram = PersistentHomology::compute(&filtration);
379
380        // Should have H0 features (components merging)
381        let h0_count = diagram.pairs_of_dim(0).count();
382        assert!(h0_count > 0);
383    }
384
385    #[test]
386    fn test_filter_by_persistence() {
387        let mut diagram = PersistenceDiagram::new();
388        diagram.add(BirthDeathPair::finite(0, 0.0, 0.1));
389        diagram.add(BirthDeathPair::finite(0, 0.0, 1.0));
390        diagram.add(BirthDeathPair::essential(0, 0.0));
391
392        let filtered = diagram.filter_by_persistence(0.5);
393        assert_eq!(filtered.pairs.len(), 2); // Only persistence >= 0.5
394    }
395
396    #[test]
397    fn test_feature_counts() {
398        let mut diagram = PersistenceDiagram::new();
399        diagram.add(BirthDeathPair::finite(0, 0.0, 1.0));
400        diagram.add(BirthDeathPair::finite(0, 0.0, 1.0));
401        diagram.add(BirthDeathPair::finite(1, 0.0, 1.0));
402
403        let counts = diagram.feature_counts();
404        assert_eq!(counts[0], 2);
405        assert_eq!(counts[1], 1);
406    }
407}
ruvector_math/homology/persistence.rs

ruvector_math/homology/
persistence.rs
