filtration_domination/points/
mod.rs

1//! Point clouds: create and modify them.
2use num::Float;
3use ordered_float::OrderedFloat;
4use rand::distributions::Distribution;
5use rand::Rng;
6use std::fmt::Formatter;
7
8use crate::distance_matrix::DistanceMatrix;
9
10pub mod input;
11pub mod output;
12
13/// A point in `R^N`.
14#[derive(Debug, Copy, Clone, Eq, PartialEq)]
15pub struct Point<T, const N: usize>(pub [T; N]);
16
17impl<T: Float, const N: usize> Point<T, N> {
18    /// Computes the Euclidean distance between the given points.
19    pub fn euclidean_distance(&self, other: &Point<T, N>) -> T {
20        let mut d = T::zero();
21        for i in 0..N {
22            d = d + (self.0[i] - other.0[i]).powi(2);
23        }
24        d.sqrt()
25    }
26
27    /// Computes the norm of the point.
28    pub fn norm(&self) -> T {
29        let mut d = T::zero();
30        for i in 0..N {
31            d = d + (self.0[i]).powi(2);
32        }
33        d.sqrt()
34    }
35
36    /// Compute the unit vector of the same direction.
37    pub fn normalize(&mut self) {
38        let norm = self.norm();
39        for i in 0..N {
40            self.0[i] = self.0[i] / norm;
41        }
42    }
43
44    /// Sample a random point, coordinate by coordinate.
45    pub fn random<D: Distribution<T>, R: Rng>(distribution: &D, rng: &mut R) -> Point<T, N> {
46        let mut p = Point([T::zero(); N]);
47        for x in p.0.iter_mut() {
48            *x = rng.sample(distribution);
49        }
50        p
51    }
52}
53
54impl<T: Float + std::fmt::Display, const N: usize> std::fmt::Display for Point<T, N> {
55    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
56        write!(f, "(")?;
57        for i in 0..N {
58            write!(f, "{}", self.0[i])?;
59            if i != N - 1 {
60                write!(f, ", ")?;
61            }
62        }
63        write!(f, ")")
64    }
65}
66
67impl<T: Float, const N: usize> std::ops::Sub for Point<T, N> {
68    type Output = Point<T, N>;
69
70    fn sub(self, rhs: Self) -> Self::Output {
71        let mut point = Point([T::zero(); N]);
72        for i in 0..N {
73            point.0[i] = self.0[i] - rhs.0[i];
74        }
75        point
76    }
77}
78
79impl<T: Float, S: Copy, const N: usize> From<[S; N]> for Point<T, N>
80where
81    S: Into<T>,
82{
83    fn from(x: [S; N]) -> Self {
84        let mut point = Point([T::zero(); N]);
85        for (i, &x_i) in x.iter().enumerate() {
86            point.0[i] = x_i.into();
87        }
88        point
89    }
90}
91
92/// A collection of points.
93pub struct PointCloud<T: Float, const N: usize>(pub Vec<Point<T, N>>);
94
95impl<T: Float, const N: usize> Default for PointCloud<T, N> {
96    fn default() -> Self {
97        Self::new()
98    }
99}
100
101impl<T: Float, const N: usize> PointCloud<T, N> {
102    /// Create a new empty point cloud.
103    pub fn new() -> Self {
104        Self(Vec::new())
105    }
106
107    /// Add a new point.
108    pub fn push_point(&mut self, p: Point<T, N>) {
109        self.0.push(p)
110    }
111
112    /// Return the distance matrix of the point cloud, where the order is the order in which the
113    /// points where added.
114    pub fn distance_matrix(&self) -> DistanceMatrix<T> {
115        let n = self.len();
116        let mut matrix = DistanceMatrix::new(n);
117        for u in 0..n {
118            for v in (u + 1)..n {
119                matrix.set(u, v, self.0[u].euclidean_distance(&self.0[v]))
120            }
121        }
122        matrix
123    }
124
125    /// Returns the number of points in the point cloud.
126    pub fn len(&self) -> usize {
127        self.0.len()
128    }
129
130    /// Returns whether the point cloud has no points.
131    pub fn is_empty(&self) -> bool {
132        self.0.is_empty()
133    }
134}
135
136impl<const N: usize> From<PointCloud<f64, N>> for PointCloud<OrderedFloat<f64>, N> {
137    fn from(points: PointCloud<f64, N>) -> Self {
138        let mut result: PointCloud<OrderedFloat<f64>, N> = PointCloud::new();
139        for p in points.0.into_iter() {
140            result.push_point(p.0.try_into().unwrap());
141        }
142        result
143    }
144}