holos_tda/distances/
matrix.rs

1/// Symmetric dissimilarity matrix.
2/// No metric assumptions: entries need not satisfy the triangle inequality.
3/// Entries must be non-negative and not NaN; +inf is legal and equivalent
4/// to an absent edge.
5///
6/// The matrix has two storage forms. The compact form holds the condensed
7/// lower triangle, `n(n-1)/2` entries, and is the form every constructor
8/// builds. The full form holds both triangles row-major, `n * n` entries,
9/// so that a cofacet diameter fold reads one contiguous row per simplex
10/// vertex instead of one strided column.
11#[derive(Debug, Clone)]
12pub struct DistanceMatrix {
13    pub(super) n: usize,
14    /// Row-major n by n in the full form, the condensed lower triangle in
15    /// the compact one. Row `i` holds its entries below the diagonal first
16    /// in both, which is what [`DistanceMatrix::lower_row`] returns.
17    pub(super) data: Vec<f64>,
18    pub(super) square: bool,
19}
20
21impl DistanceMatrix {
22    /// Number of points.
23    pub fn len(&self) -> usize {
24        self.n
25    }
26
27    /// True when there are no points.
28    pub fn is_empty(&self) -> bool {
29        self.n == 0
30    }
31
32    /// Distance between points `i` and `j` (0 on the diagonal).
33    #[inline]
34    pub fn get(&self, i: usize, j: usize) -> f64 {
35        debug_assert!(i < self.n && j < self.n);
36        if self.square {
37            return self.data[i * self.n + j];
38        }
39        match i.cmp(&j) {
40            std::cmp::Ordering::Equal => 0.0,
41            std::cmp::Ordering::Greater => self.data[i * (i - 1) / 2 + j],
42            std::cmp::Ordering::Less => self.data[j * (j - 1) / 2 + i],
43        }
44    }
45
46    /// Row `i` up to the diagonal: the distances from `i` to every point
47    /// below it, in index order. Both forms store that run contiguously.
48    #[inline]
49    pub(super) fn lower_row(&self, i: usize) -> &[f64] {
50        let start = if self.square {
51            i * self.n
52        } else {
53            i * (i - 1) / 2
54        };
55        &self.data[start..start + i]
56    }
57
58    /// Count the pairs that enter the complex at `threshold`: finite and at
59    /// or below it. One pass over the condensed triangle, no allocation.
60    pub(crate) fn count_edges_at(&self, threshold: f64) -> usize {
61        (1..self.n)
62            .map(|i| {
63                self.lower_row(i)
64                    .iter()
65                    .filter(|d| d.is_finite() && **d <= threshold)
66                    .count()
67            })
68            .sum()
69    }
70    /// Default threshold: the minimum over i of the maximum over j of d(i, j).
71    /// Past that radius the complex is a cone and acquires no further homology.
72    pub fn enclosing_radius(&self) -> f64 {
73        if self.n < 2 {
74            return 0.0;
75        }
76        // Each distance folds into both endpoints' running maxima, so one
77        // pass over the lower triangle is enough. Row `i` holds its own
78        // maximum in a local until the row ends: no earlier row writes
79        // `row_max[i]`, because every column index it touches is below it.
80        let mut row_max = vec![0.0f64; self.n];
81        for i in 1..self.n {
82            let mut max_i = 0.0f64;
83            for (m, &d) in row_max[..i].iter_mut().zip(self.lower_row(i)) {
84                max_i = max_i.max(d);
85                *m = m.max(d);
86            }
87            row_max[i] = max_i;
88        }
89        row_max.into_iter().fold(f64::INFINITY, f64::min)
90    }
91}
92
93#[cfg(test)]
94thread_local! {
95    /// Conversions to the full form on this thread. Each test runs on its
96    /// own thread, so the count belongs to one test and no other test can
97    /// disturb it.
98    pub(crate) static SQUARE_BUILDS: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
99}
100
101impl DistanceMatrix {
102    /// True when this matrix holds both triangles.
103    #[cfg(test)]
104    pub(crate) fn is_square(&self) -> bool {
105        self.square
106    }
107
108    /// The same distances in the full row-major form. Every constructor
109    /// builds the compact form, so this is the only way a run reaches the
110    /// full one.
111    pub(crate) fn to_square(&self) -> Self {
112        #[cfg(test)]
113        SQUARE_BUILDS.with(|c| c.set(c.get() + 1));
114        let n = self.n;
115        let mut data = vec![0.0f64; n * n];
116        for i in 1..n {
117            let row = self.lower_row(i);
118            data[i * n..i * n + i].copy_from_slice(row);
119            for (j, &d) in row.iter().enumerate() {
120                data[j * n + i] = d;
121            }
122        }
123        Self {
124            n,
125            data,
126            square: true,
127        }
128    }
129}
holos_tda/distances/matrix.rs

holos_tda/distances/
matrix.rs
