Skip to main content

holos_tda/filtration/
flag.rs

1use std::collections::BTreeMap;
2
3use crate::SparseDistanceMatrix;
4
5use super::grade::{FiltrationError, ScalarGrade};
6use super::simplex::{FilteredSimplex, FilteredSimplicialComplex};
7
8/// Per-dimension bounds for an explicit complex materialization.
9#[derive(Debug, Clone, Copy, PartialEq, Eq)]
10pub struct ComplexLimits {
11    /// Largest accepted vertex count.
12    pub max_vertices: usize,
13    /// Largest accepted edge count.
14    pub max_edges: usize,
15    /// Largest accepted triangle count.
16    pub max_triangles: usize,
17    /// Largest accepted simplex count in each dimension above two.
18    pub max_higher_simplices: usize,
19}
20
21impl Default for ComplexLimits {
22    fn default() -> Self {
23        Self {
24            max_vertices: 1_000_000,
25            max_edges: 20_000_000,
26            max_triangles: 100_000_000,
27            max_higher_simplices: 100_000_000,
28        }
29    }
30}
31
32impl ComplexLimits {
33    pub(crate) fn for_dimension(self, dimension: usize) -> usize {
34        match dimension {
35            0 => self.max_vertices,
36            1 => self.max_edges,
37            2 => self.max_triangles,
38            _ => self.max_higher_simplices,
39        }
40    }
41}
42
43/// Limits and filtration choices for a flag-complex materialization.
44#[derive(Debug, Clone, Copy, PartialEq)]
45pub struct FlagComplexParams {
46    /// Highest simplex dimension to materialize.
47    pub max_dimension: usize,
48    /// Largest included edge weight. `None` includes every finite listed edge.
49    pub threshold: Option<f64>,
50    /// Per-dimension resource bounds.
51    pub limits: ComplexLimits,
52}
53
54impl FilteredSimplicialComplex<ScalarGrade> {
55    /// Materialize the flag complex of a sparse graph with stable labels.
56    pub fn from_flag_graph(
57        input: &SparseDistanceMatrix,
58        labels: &[usize],
59        params: FlagComplexParams,
60    ) -> Result<Self, FiltrationError> {
61        let threshold = validate_flag_params(input, labels, params)?;
62        let zero = ScalarGrade::new(0.0)?;
63        let vertices = labels
64            .iter()
65            .map(|&vertex| FilteredSimplex::new(vec![vertex], zero))
66            .collect::<Vec<_>>();
67        let mut simplices = vec![vertices];
68        let label_to_local: BTreeMap<_, _> = labels
69            .iter()
70            .copied()
71            .enumerate()
72            .map(|(local, label)| (label, local))
73            .collect();
74        for dimension in 1..=params.max_dimension {
75            let next = extend_flag_dimension(
76                input,
77                labels,
78                &label_to_local,
79                &simplices[dimension - 1],
80                dimension,
81                threshold,
82                params.limits.for_dimension(dimension),
83            )?;
84            simplices.push(next);
85        }
86        Self::new(labels.to_vec(), simplices)
87    }
88}
89
90fn validate_flag_params(
91    input: &SparseDistanceMatrix,
92    labels: &[usize],
93    params: FlagComplexParams,
94) -> Result<f64, FiltrationError> {
95    if labels.len() != input.len() || labels.windows(2).any(|pair| pair[0] >= pair[1]) {
96        return Err(FiltrationError::new(
97            "flag-complex labels must match the graph and increase strictly",
98        ));
99    }
100    let threshold = params.threshold.unwrap_or(f64::INFINITY);
101    if threshold.is_nan() || threshold < 0.0 {
102        return Err(FiltrationError::new(
103            "a flag-complex threshold must be non-negative",
104        ));
105    }
106    if labels.len() > params.limits.max_vertices {
107        return Err(FiltrationError::new(
108            "flag-complex vertices exceed the per-dimension limit",
109        ));
110    }
111    Ok(threshold)
112}
113
114fn extend_flag_dimension(
115    input: &SparseDistanceMatrix,
116    labels: &[usize],
117    label_to_local: &BTreeMap<usize, usize>,
118    faces: &[FilteredSimplex<ScalarGrade>],
119    dimension: usize,
120    threshold: f64,
121    limit: usize,
122) -> Result<Vec<FilteredSimplex<ScalarGrade>>, FiltrationError> {
123    let mut simplices = Vec::new();
124    for simplex in faces {
125        let last_local = label_to_local[&simplex.vertices()[dimension - 1]];
126        for (local_vertex, &label) in labels.iter().enumerate().skip(last_local + 1) {
127            if let Some(coface) = extend_flag_simplex(
128                input,
129                label_to_local,
130                simplex,
131                local_vertex,
132                label,
133                threshold,
134            )? {
135                simplices.push(coface);
136                if simplices.len() > limit {
137                    return Err(FiltrationError::new(format!(
138                        "flag-complex dimension {dimension} exceeds the simplex limit"
139                    )));
140                }
141            }
142        }
143    }
144    Ok(simplices)
145}
146
147fn extend_flag_simplex(
148    input: &SparseDistanceMatrix,
149    label_to_local: &BTreeMap<usize, usize>,
150    simplex: &FilteredSimplex<ScalarGrade>,
151    local_vertex: usize,
152    label: usize,
153    threshold: f64,
154) -> Result<Option<FilteredSimplex<ScalarGrade>>, FiltrationError> {
155    let mut value = simplex.grade().value();
156    for member in simplex.vertices() {
157        let edge = input.get(label_to_local[member], local_vertex);
158        if !edge.is_finite() || edge > threshold {
159            return Ok(None);
160        }
161        value = value.max(edge);
162    }
163    let mut vertices = simplex.vertices().to_vec();
164    vertices.push(label);
165    Ok(Some(FilteredSimplex::new(
166        vertices,
167        ScalarGrade::new(value)?,
168    )))
169}