persistence_agent/
point_cloud.rs1use crate::error::PersistenceError;
2
3#[derive(Debug, Clone)]
5pub struct PointCloud {
6 pub points: Vec<Vec<f64>>,
7 pub distance_matrix: Vec<Vec<f64>>,
8}
9
10impl PointCloud {
11 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 pub fn distance(&self, i: usize, j: usize) -> f64 {
50 self.distance_matrix[i][j]
51 }
52
53 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 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 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}