1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136

use std::{collections::HashSet, ops::ControlFlow};

use kiddo::{ImmutableKdTree, SquaredEuclidean};
use nalgebra::Point3;

use crate::geometry::{approx_cmp_f64, approx_eq_point_f64, FloatEq, FloatOrdering};

#[derive(Debug, Clone, PartialEq)]
pub struct MultiCoordPoint {
    pub(crate) point: Point3<f64>,
    pub(crate) atom_ids: Vec<usize>,
}

impl MultiCoordPoint {
    pub fn new(point: Point3<f64>, atom_ids: Vec<usize>) -> Self {
        Self { point, atom_ids }
    }
    pub(crate) fn cn(&self) -> usize {
        self.atom_ids.len()
    }
    pub(crate) fn merge_with(&self, rhs: &Self) -> Option<MultiCoordPoint> {
        if let FloatEq::Eq = approx_eq_point_f64(self.point, rhs.point) {
            let new_connecting_atoms = [self.atom_ids.clone(), rhs.atom_ids.clone()].to_vec();
            let new_connecting_atoms_set: HashSet<usize> =
                new_connecting_atoms.concat().into_iter().collect();
            let mut new_connecting_atoms_array: Vec<usize> =
                new_connecting_atoms_set.into_iter().collect();
            new_connecting_atoms_array.sort();
            Some(MultiCoordPoint::new(self.point, new_connecting_atoms_array))
        } else {
            None
        }
    }

    pub fn no_closer_atoms(
        self,
        kdtree: &ImmutableKdTree<f64, 3>,
        dist: f64,
    ) -> Option<MultiCoordPoint> {
        let query: [f64; 3] = self.point.into();
        let closer_than_dist = kdtree
            .within::<SquaredEuclidean>(&query, dist.powi(2))
            .into_iter()
            .try_for_each(|nb| {
                if let FloatOrdering::Less = approx_cmp_f64(nb.distance, dist.powi(2)) {
                    ControlFlow::Break(nb)
                } else {
                    ControlFlow::Continue(())
                }
            });
        if !matches!(closer_than_dist, ControlFlow::Break(_)) {
            Some(self)
        } else {
            None
        }
    }
    pub fn dedup_points(
        points: &[MultiCoordPoint],
        kdtree: &ImmutableKdTree<f64, 3>,
        dist: f64,
    ) -> Vec<MultiCoordPoint> {
        let mut visited = vec![false; points.len()];
        points
            .iter()
            .enumerate()
            .filter_map(|(now, curr_p)| {
                Self::look_for_same_points((now, curr_p), points, &mut visited)
            })
            .filter_map(|p| p.no_closer_atoms(kdtree, dist))
            .collect()
    }
    fn look_for_same_points(
        curr: (usize, &MultiCoordPoint),
        points: &[MultiCoordPoint],
        visited: &mut [bool],
    ) -> Option<MultiCoordPoint> {
        let (now, curr_p) = curr;
        if !visited[now] {
            visited[now] = true;
            let same_points: Vec<MultiCoordPoint> = points
                .iter()
                .enumerate()
                .filter_map(|(to_check, p)| {
                    if visited[to_check] {
                        None
                    } else if let Some(merged) = curr_p.merge_with(p) {
                        visited[to_check] = true;
                        Some(merged)
                    } else {
                        None
                    }
                })
                .collect();
            if same_points.is_empty() {
                Some(curr_p.clone())
            } else {
                Some(
                    same_points
                        .into_iter()
                        .reduce(|acc, x| acc.merge_with(&x).unwrap())
                        .unwrap(),
                )
            }
        } else {
            None
        }
    }

    pub fn atom_ids(&self) -> &[usize] {
        self.atom_ids.as_ref()
    }

    pub fn point(&self) -> Point3<f64> {
        self.point
    }
}

#[derive(Debug, Clone, PartialEq)]
pub struct DelegatePoint<const N: usize> {
    pub(crate) point: Point3<f64>,
    pub(crate) atom_ids: [usize; N],
}

impl<const N: usize> DelegatePoint<N> {
    pub fn new(point: Point3<f64>, atom_ids: [usize; N]) -> Self {
        Self { point, atom_ids }
    }

    pub fn point(&self) -> Point3<f64> {
        self.point
    }

    pub fn atom_ids(&self) -> &[usize] {
        &self.atom_ids
    }
}