parse_monitors/pressure/
telescope.rs

1use super::{Record, Result};
2use flate2::read::GzDecoder;
3use itertools::{Itertools, MinMaxResult::MinMax};
4use std::{cmp::Ordering, fmt::Display, fs::File, io::Read, path::Path};
5
6fn partition(data: &[f64]) -> Option<(Vec<f64>, f64, Vec<f64>)> {
7    match data.len() {
8        0 => None,
9        _ => {
10            let (pivot_slice, tail) = data.split_at(1);
11            let pivot = pivot_slice[0];
12            let (left, right) = tail.iter().fold((vec![], vec![]), |mut splits, next| {
13                {
14                    let (ref mut left, ref mut right) = &mut splits;
15                    if next < &pivot {
16                        left.push(*next);
17                    } else {
18                        right.push(*next);
19                    }
20                }
21                splits
22            });
23
24            Some((left, pivot, right))
25        }
26    }
27}
28
29fn select(data: &[f64], k: usize) -> Option<f64> {
30    let part = partition(data);
31
32    match part {
33        None => None,
34        Some((left, pivot, right)) => {
35            let pivot_idx = left.len();
36
37            match pivot_idx.cmp(&k) {
38                Ordering::Equal => Some(pivot),
39                Ordering::Greater => select(&left, k),
40                Ordering::Less => select(&right, k - (pivot_idx + 1)),
41            }
42        }
43    }
44}
45
46fn median(data: &[f64]) -> Option<f64> {
47    let size = data.len();
48
49    match size {
50        even if even % 2 == 0 => {
51            let fst_med = select(data, (even / 2) - 1);
52            let snd_med = select(data, even / 2);
53
54            match (fst_med, snd_med) {
55                (Some(fst), Some(snd)) => Some((fst + snd) / 2.0),
56                _ => None,
57            }
58        }
59        odd => select(data, odd / 2),
60    }
61}
62
63/// Telescope mount surface pressure
64#[derive(Default, Debug)]
65pub struct Telescope {
66    // The pressure file
67    pub filename: String,
68    // the segment surface pressure [Pa]
69    pub pressure: Vec<f64>,
70    // the area vector along the surface normal
71    pub area_ijk: Vec<[f64; 3]>,
72    // the (x,y,z) coordinate where the pressure is applied
73    pub xyz: Vec<[f64; 3]>,
74}
75#[cfg(feature = "rstar")]
76pub mod rtree {
77    use rstar::{PointDistance, RTreeObject, AABB};
78    #[allow(dead_code)]
79    #[derive(Default, Debug, PartialEq, Clone, serde::Serialize)]
80    pub struct Node {
81        pub pressure: f64,
82        pub area_ijk: [f64; 3],
83        pub xyz: [f64; 3],
84    }
85    impl From<Vec<Node>> for super::Telescope {
86        fn from(nodes: Vec<Node>) -> Self {
87            let mut telescope: super::Telescope = Default::default();
88            nodes.into_iter().for_each(|node| {
89                telescope.pressure.push(node.pressure);
90                telescope.area_ijk.push(node.area_ijk);
91                telescope.xyz.push(node.xyz);
92            });
93            telescope
94        }
95    }
96    impl RTreeObject for Node {
97        type Envelope = AABB<[f64; 3]>;
98
99        fn envelope(&self) -> Self::Envelope {
100            AABB::from_point(self.xyz)
101        }
102    }
103    impl PointDistance for Node {
104        fn distance_2(&self, point: &[f64; 3]) -> f64 {
105            self.xyz
106                .iter()
107                .zip(point)
108                .map(|(x, x0)| x - x0)
109                .map(|x| x * x)
110                .sum()
111        }
112    }
113    impl super::Telescope {
114        /// Returns the [r-tree](https://docs.rs/rstar/latest/rstar/) of the all the pressure [Node]s
115        pub fn to_rtree(self) -> rstar::RTree<Node> {
116            let mut tree = rstar::RTree::new();
117            for i in 0..self.len() {
118                tree.insert(Node {
119                    pressure: self.pressure[i],
120                    area_ijk: self.area_ijk[i],
121                    xyz: self.xyz[i],
122                });
123            }
124            tree
125        }
126    }
127}
128impl Telescope {
129    /// Loads the pressure data
130    pub fn from_path<P: AsRef<Path>>(path: P) -> Result<Self> {
131        let data_path = Path::new(path.as_ref());
132        let mut contents = String::new();
133        let csv_file = File::open(&data_path)?;
134        let mut gz = GzDecoder::new(csv_file);
135        gz.read_to_string(&mut contents)?;
136        let mut rdr = csv::Reader::from_reader(contents.as_bytes());
137        let mut telescope: Telescope = Telescope {
138            filename: String::from(data_path.file_name().map(|x| x.to_str()).flatten().unwrap()),
139            ..Default::default()
140        };
141        for result in rdr.deserialize() {
142            let row: Record = result?;
143            telescope.pressure.push(row.pressure);
144            telescope
145                .area_ijk
146                .push([row.area_i, row.area_j, row.area_k]);
147            telescope.xyz.push([row.x, row.y, row.z]);
148        }
149        Ok(telescope)
150    }
151    pub fn len(&self) -> usize {
152        self.pressure.len()
153    }
154    pub fn is_empty(&self) -> bool {
155        self.len() == 0
156    }
157    /// Ierator over pressure
158    pub fn pressure_iter(&self) -> impl Iterator<Item = &f64> {
159        self.pressure.iter()
160    }
161    /// Returns the mean pressure
162    pub fn mean_pressure(&self) -> f64 {
163        self.pressure_iter().sum::<f64>() / self.len() as f64
164    }
165    /// Returns the median pressure
166    pub fn median_pressure(&self) -> Option<f64> {
167        median(&self.pressure)
168    }
169    /// Returns the pressure minimum and maximum
170    pub fn minmax_pressure(&self) -> Option<(f64, f64)> {
171        match self.pressure_iter().minmax() {
172            MinMax(x, y) => Some((*x, *y)),
173            _ => None,
174        }
175    }
176    /// Iterator over pressure area projected onto normal to the surface a mode location
177    pub fn area_ijk_iter(&self) -> impl Iterator<Item = &[f64; 3]> {
178        self.area_ijk.iter()
179    }
180    /// Iterator over pressure node coordinates
181    pub fn xyz_iter(&self) -> impl Iterator<Item = &[f64; 3]> {
182        self.xyz.iter()
183    }
184    /// Iterator over pressure node x coordinates
185    pub fn x_iter(&self) -> impl Iterator<Item = f64> + '_ {
186        self.xyz.iter().map(|xyz| xyz[0])
187    }
188    /// Iterator over pressure node y coordinates
189    pub fn y_iter(&self) -> impl Iterator<Item = f64> + '_ {
190        self.xyz.iter().map(|xyz| xyz[1])
191    }
192    /// Iterator over pressure node z coordinates
193    pub fn z_iter(&self) -> impl Iterator<Item = f64> + '_ {
194        self.xyz.iter().map(|xyz| xyz[2])
195    }
196    /// Returns the range of the x coordinates
197    pub fn minmax_x(&self) -> Option<(f64, f64)> {
198        match self.x_iter().minmax() {
199            MinMax(x, y) => Some((x, y)),
200            _ => None,
201        }
202    }
203    /// Returns the range of the y coordinates
204    pub fn minmax_y(&self) -> Option<(f64, f64)> {
205        match self.y_iter().minmax() {
206            MinMax(x, y) => Some((x, y)),
207            _ => None,
208        }
209    }
210    /// Returns the range of the z coordinates
211    pub fn minmax_z(&self) -> Option<(f64, f64)> {
212        match self.z_iter().minmax() {
213            MinMax(x, y) => Some((x, y)),
214            _ => None,
215        }
216    }
217    /// Returns the pressure areas
218    pub fn area_mag(&self) -> Vec<f64> {
219        self.area_ijk
220            .iter()
221            .map(|ijk| ijk.iter().map(|x| x * x).sum::<f64>().sqrt())
222            .collect()
223    }
224    /// Returns the total area
225    pub fn total_area(&self) -> f64 {
226        self.area_mag().into_iter().sum::<f64>()
227    }
228}
229impl Display for Telescope {
230    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
231        writeln!(f, "{} [{:#}]:", self.filename, self.len())?;
232        writeln!(f, " - pressure:")?;
233        writeln!(f, "  - mean  : {:.3}pa", self.mean_pressure())?;
234        writeln!(f, "  - median: {:.3}pa", self.median_pressure().unwrap())?;
235        self.minmax_pressure()
236            .map(|x| writeln!(f, "  - minmax: {:.3?}pa", x))
237            .transpose()?;
238        writeln!(f, " - total area: {:.3}m^2", self.total_area())?;
239        writeln!(f, " - volume:")?;
240        self.minmax_x()
241            .map(|x| writeln!(f, "  - x minmax: {:.3?}m", x))
242            .transpose()?;
243        self.minmax_y()
244            .map(|x| writeln!(f, "  - y minmax: {:.3?}m", x))
245            .transpose()?;
246        self.minmax_z()
247            .ok_or(std::fmt::Error)
248            .map(|x| write!(f, "  - z minmax: {:.3?}m", x))?
249    }
250}
251
252#[cfg(test)]
253mod tests {
254    use super::Telescope;
255
256    #[test]
257    fn loading() {
258        let telescope =
259            Telescope::from_path("data/Telescope_p_telescope_7.000000e+02.csv.z").unwrap();
260        println!("{telescope}");
261    }
262
263    #[cfg(feature = "rstar")]
264    #[test]
265    fn rtree() {
266        use super::rtree::Node;
267        let telescope =
268            Telescope::from_path("data/Telescope_p_telescope_7.000000e+02.csv.z").unwrap();
269        println!("{telescope}");
270        let rtree = telescope.to_rtree();
271        let node = rtree.locate_at_point(&[-6.71743755562523, 1.18707466192993, -4.88465284781676]);
272        assert_eq!(
273            node.unwrap().clone(),
274            Node {
275                pressure: -7.29796930557952,
276                area_ijk: [
277                    1.81847333261638e-06,
278                    1.84800982152254e-06,
279                    0.0126174546658134,
280                ],
281                xyz: [-6.71743755562523, 1.18707466192993, -4.88465284781676],
282            }
283        );
284    }
285}
parse_monitors/pressure/telescope.rs

parse_monitors/pressure/
telescope.rs
