holos_tda/filtration/
flag.rs1use std::collections::BTreeMap;
2
3use crate::SparseDistanceMatrix;
4
5use super::grade::{FiltrationError, ScalarGrade};
6use super::simplex::{FilteredSimplex, FilteredSimplicialComplex};
7
8#[derive(Debug, Clone, Copy, PartialEq, Eq)]
10pub struct ComplexLimits {
11 pub max_vertices: usize,
13 pub max_edges: usize,
15 pub max_triangles: usize,
17 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#[derive(Debug, Clone, Copy, PartialEq)]
45pub struct FlagComplexParams {
46 pub max_dimension: usize,
48 pub threshold: Option<f64>,
50 pub limits: ComplexLimits,
52}
53
54impl FilteredSimplicialComplex<ScalarGrade> {
55 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}