persistence_agent/
boundary.rs

1use crate::error::PersistenceError;
2use crate::vietoris_rips::VietorisRipsComplex;
3
4/// Boundary matrix over Z/2Z with column reduction for persistence computation.
5#[derive(Debug, Clone)]
6pub struct BoundaryMatrix {
7    pub matrix: Vec<Vec<u8>>,
8    pub n_simplices: usize,
9}
10
11impl BoundaryMatrix {
12    /// Build the mod-2 boundary matrix from a Vietoris-Rips complex.
13    ///
14    /// The columns are indexed by simplices (in the same order as `vr.simplices`),
15    /// and rows are indexed by the *simplices one dimension lower*.
16    ///
17    /// We use a dense representation for clarity: `matrix[col][row]` ∈ {0, 1}.
18    /// For each column (a k-simplex), the rows with value 1 correspond to the
19    /// (k−1)-faces of that simplex.
20    pub fn build(vr: &VietorisRipsComplex) -> Result<Self, PersistenceError> {
21        let n = vr.n_simplices();
22        if n == 0 {
23            return Err(PersistenceError::EmptyFiltration);
24        }
25
26        // Map each simplex to its index for fast face lookup
27        let mut simplex_index = std::collections::HashMap::new();
28        for (i, s) in vr.simplices.iter().enumerate() {
29            simplex_index.insert(s.clone(), i);
30        }
31
32        let mut matrix = vec![vec![0u8; n]; n];
33
34        for (col_idx, simplex) in vr.simplices.iter().enumerate() {
35            if simplex.len() <= 1 {
36                // Vertices have no boundary
37                continue;
38            }
39            // All (k−1)-faces: remove one vertex at a time
40            for skip in 0..simplex.len() {
41                let face: Vec<usize> = simplex
42                    .iter()
43                    .enumerate()
44                    .filter(|(i, _)| *i != skip)
45                    .map(|(_, &v)| v)
46                    .collect();
47                if let Some(&row_idx) = simplex_index.get(&face) {
48                    matrix[col_idx][row_idx] ^= 1;
49                }
50            }
51        }
52
53        Ok(Self {
54            matrix,
55            n_simplices: n,
56        })
57    }
58
59    /// Reduce the matrix using the standard algorithm.
60    /// After reduction, column `j` is either zero or its lowest nonzero row `low(j)`
61    /// is unique. Returns a mapping `low -> j` for non-zero reduced columns.
62    pub fn reduce(&mut self) -> std::collections::HashMap<usize, usize> {
63        let n = self.n_simplices;
64        let mut low_to_col: std::collections::HashMap<usize, usize> =
65            std::collections::HashMap::new();
66
67        for j in 0..n {
68            while let Some(lj) = self.low(j) {
69                if let Some(&k) = low_to_col.get(&lj) {
70                    self.add_columns(k, j);
71                } else {
72                    low_to_col.insert(lj, j);
73                    break;
74                }
75            }
76        }
77
78        low_to_col
79    }
80
81    /// Lowest nonzero row in column `j`, or None if the column is zero.
82    pub fn low(&self, j: usize) -> Option<usize> {
83        self.matrix[j]
84            .iter()
85            .enumerate()
86            .rev()
87            .find(|(_, &v)| v == 1)
88            .map(|(i, _)| i)
89    }
90
91    /// Add column `src` to column `dst` (mod 2).
92    fn add_columns(&mut self, src: usize, dst: usize) {
93        for i in 0..self.n_simplices {
94            self.matrix[dst][i] ^= self.matrix[src][i];
95        }
96    }
97}
persistence_agent/boundary.rs

persistence_agent/
boundary.rs
