holos_tda/distances/
construction.rs

1use super::matrix::DistanceMatrix;
2use super::sparse::SparseDistanceMatrix;
3use crate::{Error, Result};
4
5impl DistanceMatrix {
6    /// Euclidean distances of a point cloud. Coordinates must be finite.
7    pub fn from_points(points: &[Vec<f64>]) -> Result<Self> {
8        let n = points.len();
9        validate_points(points)?;
10        let mut data = Vec::with_capacity(n.saturating_sub(1) * n / 2);
11        for i in 1..n {
12            for j in 0..i {
13                data.push(euclidean(&points[i], &points[j]));
14            }
15        }
16        Ok(Self {
17            n,
18            data,
19            square: false,
20        })
21    }
22
23    /// Build from the condensed lower triangle, row by row: d(1,0), d(2,0),
24    /// d(2,1), d(3,0), and so on. An empty vector means one point (n = 1).
25    /// Only [`DistanceMatrix::from_points`] can build an empty *space*
26    /// (n = 0).
27    pub fn from_condensed(mut condensed: Vec<f64>) -> Result<Self> {
28        let m = condensed.len();
29        let n = ((1.0 + 8.0 * m as f64).sqrt() as usize).div_ceil(2);
30        if n * (n - 1) / 2 != m {
31            return Err(Error::InvalidInput(format!(
32                "condensed length {m} is not n(n-1)/2 for any n"
33            )));
34        }
35        for (i, d) in condensed.iter_mut().enumerate() {
36            if d.is_nan() {
37                return Err(Error::InvalidDistance(format!(
38                    "NaN at condensed index {i}"
39                )));
40            }
41            if *d < 0.0 {
42                return Err(Error::InvalidDistance(format!(
43                    "negative entry {d} at condensed index {i}"
44                )));
45            }
46            if *d == 0.0 {
47                *d = 0.0;
48            }
49        }
50        Ok(Self {
51            n,
52            data: condensed,
53            square: false,
54        })
55    }
56}
57
58impl DistanceMatrix {
59    /// The thresholded graph: the pairs [`DistanceMatrix::count_edges_at`]
60    /// counts, over the same vertex set. A vertex with no edge keeps its
61    /// place and its essential H0 bar.
62    ///
63    /// One pass over the lower triangle counts the degrees, and a second
64    /// pass files each kept pair under both of its endpoints. Row `i`
65    /// reaches vertex `v` before any later row does, and it lists the
66    /// neighbors below `v` in ascending order, so every list comes out
67    /// sorted.
68    pub(crate) fn to_sparse_at(&self, threshold: f64) -> Result<SparseDistanceMatrix> {
69        let n = self.n;
70        if n > u32::MAX as usize {
71            return Err(Error::InvalidInput(format!(
72                "sparse matrix holds at most {} points, got {n}",
73                u32::MAX
74            )));
75        }
76        let keep = |d: f64| d.is_finite() && d <= threshold;
77        let mut degree = vec![0usize; n];
78        for i in 1..n {
79            for (j, &d) in self.lower_row(i).iter().enumerate() {
80                if keep(d) {
81                    degree[i] += 1;
82                    degree[j] += 1;
83                }
84            }
85        }
86        let mut offsets = vec![0usize; n + 1];
87        let mut total = 0usize;
88        for (v, &deg) in degree.iter().enumerate() {
89            offsets[v] = total;
90            total += deg;
91        }
92        offsets[n] = total;
93
94        let mut indices = vec![0u32; total];
95        let mut values = vec![0.0f64; total];
96        let mut cursor = offsets[..n].to_vec();
97        let mut max_distance = 0.0f64;
98        for i in 1..n {
99            for (j, &d) in self.lower_row(i).iter().enumerate() {
100                if !keep(d) {
101                    continue;
102                }
103                max_distance = max_distance.max(d);
104                indices[cursor[i]] = j as u32;
105                values[cursor[i]] = d;
106                cursor[i] += 1;
107                indices[cursor[j]] = i as u32;
108                values[cursor[j]] = d;
109                cursor[j] += 1;
110            }
111        }
112        Ok(SparseDistanceMatrix {
113            n,
114            offsets,
115            indices,
116            values,
117            max_distance,
118        })
119    }
120}
121
122pub(super) fn validate_points(points: &[Vec<f64>]) -> Result<usize> {
123    let dimensions = points.first().map_or(0, Vec::len);
124    if let Some(point) = points.iter().find(|point| point.len() != dimensions) {
125        return Err(Error::InvalidInput(format!(
126            "inconsistent point dimensions: {} vs {}",
127            dimensions,
128            point.len()
129        )));
130    }
131    if points
132        .iter()
133        .flatten()
134        .any(|coordinate| !coordinate.is_finite())
135    {
136        return Err(Error::InvalidInput("non-finite coordinate".into()));
137    }
138    Ok(dimensions)
139}
140
141/// Scaled two-norm: exact where the naive sum of squares would overflow or
142/// underflow. Finite coordinates whose difference still overflows f64 give
143/// +inf.
144pub(super) fn euclidean(a: &[f64], b: &[f64]) -> f64 {
145    let m = a
146        .iter()
147        .zip(b)
148        .map(|(x, y)| (x - y).abs())
149        .fold(0.0f64, f64::max);
150    if m == 0.0 {
151        return 0.0;
152    }
153    if m.is_infinite() {
154        return f64::INFINITY;
155    }
156    let s: f64 = a
157        .iter()
158        .zip(b)
159        .map(|(x, y)| {
160            let r = (x - y) / m;
161            r * r
162        })
163        .sum();
164    m * s.sqrt()
165}
holos_tda/distances/construction.rs

holos_tda/distances/
construction.rs
