holos_tda/bipersistence/
module.rs

1//! Construction and accessors for finite bipersistence modules.
2
3use std::collections::BTreeMap;
4use std::sync::Arc;
5
6use crate::bifiltration::DegreeRipsBifiltration;
7use crate::cohomology::{CohomologySpace, cohomology_restriction, cohomology_space};
8use crate::{Error, Result, SparseDistanceMatrix};
9
10use super::linear::{LinearMap, linear_from_restriction, public_terms, rank};
11use super::{
12    Bigrade, BipersistenceLimits, BipersistenceMap, BipersistenceMapColumn, BipersistenceNode,
13};
14
15/// Exact finite `H¹` module of a degree-Rips bifiltration.
16///
17/// Cover maps are cohomology restrictions. Their ranks equal the ranks of the
18/// dual homology inclusion maps over the same field.
19#[derive(Debug, Clone)]
20pub struct BipersistenceModule {
21    pub(super) degree_rips: DegreeRipsBifiltration,
22    pub(super) modulus: u32,
23    pub(super) scales: Vec<f64>,
24    pub(super) minimum_degrees: Vec<usize>,
25    pub(super) nodes: Vec<BipersistenceNode>,
26    pub(super) cover_maps: Vec<BipersistenceMap>,
27    pub(super) cover_positions: BTreeMap<(Bigrade, Bigrade), usize>,
28    pub(super) graphs: Vec<Arc<SparseDistanceMatrix>>,
29    pub(super) spaces: Vec<Arc<CohomologySpace>>,
30    pub(super) limits: BipersistenceLimits,
31}
32
33impl BipersistenceModule {
34    /// Construct the complete finite `H¹` module on a degree-Rips grid.
35    ///
36    /// Construction checks every cover restriction and every commutative
37    /// square. The degree-Rips input must include homology dimension one.
38    pub fn from_degree_rips(
39        degree_rips: &DegreeRipsBifiltration,
40        modulus: u32,
41        limits: BipersistenceLimits,
42    ) -> Result<Self> {
43        validate_h1_dimension(degree_rips)?;
44        let bifiltration = degree_rips.bifiltration();
45        let scales = bifiltration.scales().collect::<Vec<_>>();
46        let minimum_degrees = bifiltration.minimum_degrees().to_vec();
47        let (node_count, cover_count) = grid_sizes(&scales, &minimum_degrees, limits)?;
48        let BuiltNodes {
49            nodes,
50            graphs,
51            spaces,
52            space_positions,
53        } = build_nodes(
54            bifiltration,
55            &scales,
56            &minimum_degrees,
57            modulus,
58            limits,
59            node_count,
60        )?;
61
62        let mut module = Self {
63            degree_rips: degree_rips.clone(),
64            modulus,
65            scales,
66            minimum_degrees,
67            nodes,
68            cover_maps: Vec::with_capacity(cover_count),
69            cover_positions: BTreeMap::new(),
70            graphs,
71            spaces,
72            limits,
73        };
74        module.populate_covers(&space_positions)?;
75        module.validate_map_terms()?;
76        module.check_squares()?;
77        Ok(module)
78    }
79
80    /// Prime coefficient modulus.
81    pub fn modulus(&self) -> u32 {
82        self.modulus
83    }
84
85    /// Scale values in strict ascending order.
86    pub fn scales(&self) -> &[f64] {
87        &self.scales
88    }
89
90    /// Minimum degrees in strict descending order.
91    pub fn minimum_degrees(&self) -> &[usize] {
92        &self.minimum_degrees
93    }
94
95    /// Nodes in scale-major, then density-major order.
96    pub fn nodes(&self) -> &[BipersistenceNode] {
97        &self.nodes
98    }
99
100    /// Checked horizontal and vertical cover maps.
101    pub fn cover_maps(&self) -> &[BipersistenceMap] {
102        &self.cover_maps
103    }
104
105    /// Source degree-Rips bifiltration.
106    pub fn degree_rips(&self) -> &DegreeRipsBifiltration {
107        &self.degree_rips
108    }
109
110    /// Return one module node.
111    pub fn node(&self, grade: Bigrade) -> Result<&BipersistenceNode> {
112        self.validate_grade(grade)?;
113        Ok(&self.nodes[self.node_index(grade)])
114    }
115
116    /// Return the zero-weight graph used for `H¹` at one grid node.
117    ///
118    /// The graph retains inactive vertices as isolated labels. It represents
119    /// positive-dimensional flag cohomology, not `H⁰` of the degree-Rips
120    /// slice.
121    pub fn h1_graph(&self, grade: Bigrade) -> Result<&SparseDistanceMatrix> {
122        self.validate_grade(grade)?;
123        Ok(self.graphs[self.node_index(grade)].as_ref())
124    }
125
126    /// Return the canonical `H¹` space at one grid node.
127    pub fn cohomology_space(&self, grade: Bigrade) -> Result<&CohomologySpace> {
128        self.validate_grade(grade)?;
129        Ok(self.spaces[self.node_index(grade)].as_ref())
130    }
131
132    /// Return one checked cover map.
133    pub fn cover_map(&self, lower: Bigrade, upper: Bigrade) -> Result<&BipersistenceMap> {
134        let position = self.cover_positions.get(&(lower, upper)).ok_or_else(|| {
135            Error::InvalidInput("the requested grades do not form a grid cover".into())
136        })?;
137        Ok(&self.cover_maps[*position])
138    }
139
140    pub(super) fn scale_count(&self) -> usize {
141        self.scales.len()
142    }
143
144    pub(super) fn density_count(&self) -> usize {
145        self.minimum_degrees.len()
146    }
147
148    pub(super) fn node_index(&self, grade: Bigrade) -> usize {
149        grade.scale() * self.density_count() + grade.density()
150    }
151
152    pub(super) fn validate_grade(&self, grade: Bigrade) -> Result<()> {
153        if grade.scale() >= self.scale_count() || grade.density() >= self.density_count() {
154            return Err(Error::InvalidInput(
155                "bipersistence grade is outside the finite parameter grid".into(),
156            ));
157        }
158        Ok(())
159    }
160
161    pub(super) fn validate_comparable(&self, lower: Bigrade, upper: Bigrade) -> Result<()> {
162        self.validate_grade(lower)?;
163        self.validate_grade(upper)?;
164        if !lower.precedes(upper) {
165            return Err(Error::InvalidInput(
166                "bipersistence map grades are not comparable".into(),
167            ));
168        }
169        Ok(())
170    }
171
172    fn push_cover(
173        &mut self,
174        lower: Bigrade,
175        upper: Bigrade,
176        space_positions: &[usize],
177        map_cache: &mut BTreeMap<(usize, usize), Arc<LinearMap>>,
178    ) -> Result<()> {
179        let lower_position = self.node_index(lower);
180        let upper_position = self.node_index(upper);
181        let key = (
182            space_positions[upper_position],
183            space_positions[lower_position],
184        );
185        let linear = if let Some(linear) = map_cache.get(&key) {
186            Arc::clone(linear)
187        } else {
188            let restriction = cohomology_restriction(
189                self.graphs[upper_position].as_ref(),
190                self.spaces[upper_position].as_ref(),
191                self.graphs[lower_position].as_ref(),
192                self.spaces[lower_position].as_ref(),
193            )?;
194            let linear = Arc::new(linear_from_restriction(
195                &restriction,
196                self.spaces[upper_position].as_ref(),
197                self.spaces[lower_position].as_ref(),
198            )?);
199            map_cache.insert(key, Arc::clone(&linear));
200            linear
201        };
202        let map = self.public_map(lower, upper, linear.as_ref());
203        let position = self.cover_maps.len();
204        self.cover_positions.insert((lower, upper), position);
205        self.cover_maps.push(map);
206        Ok(())
207    }
208
209    pub(super) fn cover_linear(&self, lower: Bigrade, upper: Bigrade) -> Result<LinearMap> {
210        self.linear_from_public(self.cover_map(lower, upper)?)
211    }
212
213    pub(super) fn linear_from_public(&self, map: &BipersistenceMap) -> Result<LinearMap> {
214        let source_rank = self.node(map.upper_grade)?.rank;
215        let target_rank = self.node(map.lower_grade)?.rank;
216        if map.columns.len() != source_rank {
217            return Err(Error::InvalidInput(
218                "bipersistence map has the wrong source rank".into(),
219            ));
220        }
221        let mut columns = Vec::with_capacity(source_rank);
222        for (source, column) in map.columns.iter().enumerate() {
223            if column.source_basis_index != source {
224                return Err(Error::InvalidInput(
225                    "bipersistence map columns are not canonical".into(),
226                ));
227            }
228            columns.push(super::linear::coordinate_vector(
229                &column.image,
230                target_rank,
231                self.modulus,
232                "bipersistence map image is not canonical",
233            )?);
234        }
235        Ok(LinearMap {
236            source_rank,
237            target_rank,
238            columns,
239        })
240    }
241
242    pub(super) fn public_map(
243        &self,
244        lower: Bigrade,
245        upper: Bigrade,
246        linear: &LinearMap,
247    ) -> BipersistenceMap {
248        BipersistenceMap {
249            lower_grade: lower,
250            upper_grade: upper,
251            source_space: self.nodes[self.node_index(upper)].space,
252            target_space: self.nodes[self.node_index(lower)].space,
253            rank: rank(linear.columns.clone(), linear.target_rank, self.modulus),
254            columns: linear
255                .columns
256                .iter()
257                .enumerate()
258                .map(|(source_basis_index, image)| BipersistenceMapColumn {
259                    source_basis_index,
260                    image: public_terms(image),
261                })
262                .collect(),
263        }
264    }
265
266    fn check_squares(&self) -> Result<()> {
267        for scale in 0..self.scale_count().saturating_sub(1) {
268            for density in 0..self.density_count().saturating_sub(1) {
269                let lower = Bigrade::new(scale, density);
270                let upper = Bigrade::new(scale + 1, density + 1);
271                let scale_first = LinearMap::compose(
272                    &self.cover_linear(lower, Bigrade::new(scale, density + 1))?,
273                    &self.cover_linear(Bigrade::new(scale, density + 1), upper)?,
274                    self.modulus,
275                )?;
276                let density_first = LinearMap::compose(
277                    &self.cover_linear(lower, Bigrade::new(scale + 1, density))?,
278                    &self.cover_linear(Bigrade::new(scale + 1, density), upper)?,
279                    self.modulus,
280                )?;
281                if scale_first != density_first {
282                    return Err(Error::InvalidInput(format!(
283                        "bipersistence square at ({scale}, {density}) does not commute"
284                    )));
285                }
286            }
287        }
288        Ok(())
289    }
290}
291
292fn validate_h1_dimension(degree_rips: &DegreeRipsBifiltration) -> Result<()> {
293    if degree_rips.max_homology_dimension() < 1 {
294        return Err(Error::InvalidInput(
295            "an H1 bipersistence module needs degree-Rips dimension one".into(),
296        ));
297    }
298    Ok(())
299}
300
301fn grid_sizes(
302    scales: &[f64],
303    minimum_degrees: &[usize],
304    limits: BipersistenceLimits,
305) -> Result<(usize, usize)> {
306    let node_count = scales
307        .len()
308        .checked_mul(minimum_degrees.len())
309        .ok_or_else(|| Error::InvalidInput("bipersistence grid size overflows".into()))?;
310    if node_count > limits.max_nodes {
311        return Err(Error::InvalidInput(format!(
312            "bipersistence node count exceeds the limit {}",
313            limits.max_nodes
314        )));
315    }
316    let horizontal = scales
317        .len()
318        .saturating_sub(1)
319        .checked_mul(minimum_degrees.len())
320        .ok_or_else(|| Error::InvalidInput("bipersistence cover count overflows".into()))?;
321    let vertical = minimum_degrees
322        .len()
323        .saturating_sub(1)
324        .checked_mul(scales.len())
325        .ok_or_else(|| Error::InvalidInput("bipersistence cover count overflows".into()))?;
326    let cover_count = horizontal
327        .checked_add(vertical)
328        .ok_or_else(|| Error::InvalidInput("bipersistence cover count overflows".into()))?;
329    if cover_count > limits.max_cover_maps {
330        return Err(Error::InvalidInput(format!(
331            "bipersistence cover count exceeds the limit {}",
332            limits.max_cover_maps
333        )));
334    }
335    Ok((node_count, cover_count))
336}
337
338fn build_nodes(
339    bifiltration: &crate::bifiltration::MulticriticalBifiltration,
340    scales: &[f64],
341    minimum_degrees: &[usize],
342    modulus: u32,
343    limits: BipersistenceLimits,
344    node_count: usize,
345) -> Result<BuiltNodes> {
346    let mut nodes = Vec::with_capacity(node_count);
347    let mut graphs = Vec::with_capacity(node_count);
348    let mut spaces = Vec::with_capacity(node_count);
349    let mut space_positions = Vec::with_capacity(node_count);
350    let mut graph_cache = BTreeMap::new();
351    let mut total_rank = 0usize;
352    for scale in 0..scales.len() {
353        for density in 0..minimum_degrees.len() {
354            let grade = Bigrade::new(scale, density);
355            let slice = bifiltration.slice(grade)?;
356            let graph = slice.h1_graph(bifiltration.vertex_count())?;
357            let key = ActiveGraphKey::from_graph(&graph);
358            let (graph, space, space_position) =
359                if let Some((cached_graph, cached_space, space_position)) = graph_cache.get(&key) {
360                    (
361                        Arc::clone(cached_graph),
362                        Arc::clone(cached_space),
363                        *space_position,
364                    )
365                } else {
366                    let graph = Arc::new(graph);
367                    let space = Arc::new(cohomology_space(
368                        graph.as_ref(),
369                        1,
370                        0.0,
371                        modulus,
372                        limits.cohomology,
373                    )?);
374                    let space_position = graph_cache.len();
375                    graph_cache.insert(
376                        key,
377                        (Arc::clone(&graph), Arc::clone(&space), space_position),
378                    );
379                    (graph, space, space_position)
380                };
381            total_rank = total_rank
382                .checked_add(space.rank())
383                .ok_or_else(|| Error::InvalidInput("bipersistence total rank overflows".into()))?;
384            if total_rank > limits.max_total_rank {
385                return Err(Error::InvalidInput(format!(
386                    "bipersistence total rank exceeds the limit {}",
387                    limits.max_total_rank
388                )));
389            }
390            nodes.push(BipersistenceNode {
391                grade,
392                space: space.id(),
393                rank: space.rank(),
394            });
395            graphs.push(graph);
396            spaces.push(space);
397            space_positions.push(space_position);
398        }
399    }
400    Ok(BuiltNodes {
401        nodes,
402        graphs,
403        spaces,
404        space_positions,
405    })
406}
407
408struct BuiltNodes {
409    nodes: Vec<BipersistenceNode>,
410    graphs: Vec<Arc<SparseDistanceMatrix>>,
411    spaces: Vec<Arc<CohomologySpace>>,
412    space_positions: Vec<usize>,
413}
414
415#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord)]
416struct ActiveGraphKey {
417    vertex_count: usize,
418    edges: Vec<(usize, usize, u64)>,
419}
420
421impl ActiveGraphKey {
422    fn from_graph(graph: &SparseDistanceMatrix) -> Self {
423        Self {
424            vertex_count: graph.len(),
425            edges: graph
426                .edges()
427                .map(|(u, v, value)| (u, v, value.to_bits()))
428                .collect(),
429        }
430    }
431}
432
433impl BipersistenceModule {
434    fn populate_covers(&mut self, space_positions: &[usize]) -> Result<()> {
435        let mut map_cache = BTreeMap::new();
436        for scale in 0..self.scale_count() {
437            for density in 0..self.density_count() {
438                let lower = Bigrade::new(scale, density);
439                if scale + 1 < self.scale_count() {
440                    self.push_cover(
441                        lower,
442                        Bigrade::new(scale + 1, density),
443                        space_positions,
444                        &mut map_cache,
445                    )?;
446                }
447                if density + 1 < self.density_count() {
448                    self.push_cover(
449                        lower,
450                        Bigrade::new(scale, density + 1),
451                        space_positions,
452                        &mut map_cache,
453                    )?;
454                }
455            }
456        }
457        Ok(())
458    }
459
460    fn validate_map_terms(&self) -> Result<()> {
461        let map_terms = self
462            .cover_maps
463            .iter()
464            .flat_map(|map| &map.columns)
465            .map(|column| column.image.len())
466            .try_fold(0usize, |total, count| total.checked_add(count))
467            .ok_or_else(|| Error::InvalidInput("bipersistence map term count overflows".into()))?;
468        if map_terms > self.limits.max_map_terms {
469            return Err(Error::InvalidInput(format!(
470                "bipersistence map term count exceeds the limit {}",
471                self.limits.max_map_terms
472            )));
473        }
474        Ok(())
475    }
476}
holos_tda/bipersistence/module.rs

holos_tda/bipersistence/
module.rs
