Skip to main content

persistence_agent/
point_cloud.rs

1use crate::error::PersistenceError;
2
3/// A set of points in Euclidean space with a precomputed distance matrix.
4#[derive(Debug, Clone)]
5pub struct PointCloud {
6    pub points: Vec<Vec<f64>>,
7    pub distance_matrix: Vec<Vec<f64>>,
8}
9
10impl PointCloud {
11    /// Build a point cloud, computing the pairwise Euclidean distance matrix.
12    pub fn new(points: Vec<Vec<f64>>) -> Result<Self, PersistenceError> {
13        if points.is_empty() {
14            return Err(PersistenceError::EmptyCloud);
15        }
16        let dim = points[0].len();
17        for p in points.iter() {
18            if p.len() != dim {
19                return Err(PersistenceError::DimensionMismatch {
20                    expected: dim,
21                    actual: p.len(),
22                });
23            }
24        }
25        let n = points.len();
26        let mut dm = vec![vec![0.0; n]; n];
27        for i in 0..n {
28            for j in (i + 1)..n {
29                let d = euclidean(&points[i], &points[j]);
30                dm[i][j] = d;
31                dm[j][i] = d;
32            }
33        }
34        Ok(Self {
35            points,
36            distance_matrix: dm,
37        })
38    }
39
40    pub fn n_points(&self) -> usize {
41        self.points.len()
42    }
43
44    pub fn dimension(&self) -> usize {
45        self.points[0].len()
46    }
47
48    /// Distance between points at indices `i` and `j`.
49    pub fn distance(&self, i: usize, j: usize) -> f64 {
50        self.distance_matrix[i][j]
51    }
52
53    /// k-nearest-neighbor graph: returns for each point the indices of its k nearest
54    /// neighbors (excluding itself), sorted by distance ascending.
55    pub fn knn(&self, k: usize) -> Result<Vec<Vec<usize>>, PersistenceError> {
56        let n = self.n_points();
57        if k == 0 || k >= n {
58            return Err(PersistenceError::InvalidK { k, n });
59        }
60        let mut result = Vec::with_capacity(n);
61        for i in 0..n {
62            let mut pairs: Vec<(f64, usize)> = (0..n)
63                .filter(|&j| j != i)
64                .map(|j| (self.distance_matrix[i][j], j))
65                .collect();
66            pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
67            result.push(pairs[..k].iter().map(|(_, j)| *j).collect());
68        }
69        Ok(result)
70    }
71
72    /// Maximum pairwise distance in the cloud.
73    pub fn max_distance(&self) -> f64 {
74        let n = self.n_points();
75        let mut mx: f64 = 0.0;
76        for i in 0..n {
77            for j in (i + 1)..n {
78                mx = mx.max(self.distance_matrix[i][j]);
79            }
80        }
81        mx
82    }
83
84    /// Return the sorted list of unique pairwise distances (filtration thresholds).
85    pub fn unique_distances(&self) -> Vec<f64> {
86        let n = self.n_points();
87        let mut dists = Vec::new();
88        for i in 0..n {
89            for j in (i + 1)..n {
90                dists.push(self.distance_matrix[i][j]);
91            }
92        }
93        dists.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
94        dists.dedup_by(|a, b| (*a - *b).abs() < 1e-12);
95        dists
96    }
97}
98
99fn euclidean(a: &[f64], b: &[f64]) -> f64 {
100    a.iter()
101        .zip(b.iter())
102        .map(|(x, y)| (x - y) * (x - y))
103        .sum::<f64>()
104        .sqrt()
105}