mzpeaks/feature/
simple.rs

1use std::{
2    cmp::Ordering,
3    marker::PhantomData,
4    ops::{Bound, RangeBounds},
5};
6
7#[cfg(feature = "serde")]
8use serde::{Deserialize, Serialize};
9
10use crate::{coordinate::CoordinateLike, IntensityMeasurement};
11
12use super::util::{NonNan, EMPTY_Y, EMPTY_Z};
13use super::{
14    traits::{CoArrayOps, FeatureLike, SplittableFeatureLike, TimeInterval},
15    FeatureLikeMut, TimeArray,
16};
17
18/// A feature-like type that doesn't have a variable first dimension, instead
19/// using a constant value
20#[derive(Debug, Default, Clone)]
21#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
22pub struct SimpleFeature<X, Y> {
23    /// A constant value for the first dimension
24    pub label: f64,
25    y: Vec<f64>,
26    z: Vec<f32>,
27    _x: PhantomData<X>,
28    _y: PhantomData<Y>,
29}
30
31impl<X, Y> CoArrayOps for SimpleFeature<X, Y> {}
32
33impl<X, Y> SimpleFeature<X, Y> {
34    pub fn new(label: f64, y: Vec<f64>, z: Vec<f32>) -> Self {
35        Self {
36            label,
37            y,
38            z,
39            _x: PhantomData,
40            _y: PhantomData,
41        }
42    }
43
44    pub fn as_view(&self) -> SimpleFeatureView<'_, X, Y> {
45        SimpleFeatureView::new(self.label, &self.y, &self.z)
46    }
47
48    pub fn into_inner(self) -> (f64, Vec<f64>, Vec<f32>) {
49        (self.label, self.y, self.z)
50    }
51
52    /// Create an empty [`SimpleFeature`]
53    pub fn empty(label: f64) -> Self {
54        Self::with_capacity(0, label)
55    }
56
57    /// Createa a new empty [`SimpleFeature`] with pre-allocated capacity
58    pub fn with_capacity(capacity: usize, label: f64) -> Self {
59        Self {
60            label,
61            y: Vec::with_capacity(capacity),
62            z: Vec::with_capacity(capacity),
63            _x: PhantomData,
64            _y: PhantomData,
65        }
66    }
67
68    fn sort_by_y(&mut self) {
69        let n = self.len();
70        let mut mask: Vec<_> = (0..n).collect();
71        mask.sort_by_key(|i| NonNan::new(self.y[*i]));
72
73        const TOMBSTONE: usize = usize::MAX;
74
75        for idx in 0..n {
76            if mask[idx] != TOMBSTONE {
77                let mut current_idx = idx;
78                loop {
79                    let next_idx = mask[current_idx];
80                    mask[current_idx] = TOMBSTONE;
81                    if mask[next_idx] == TOMBSTONE {
82                        break;
83                    }
84                    self.y.swap(current_idx, next_idx);
85                    self.z.swap(current_idx, next_idx);
86                    current_idx = next_idx;
87                }
88            }
89        }
90    }
91
92    /// Push a new data point onto the feature and ensure the time ordering invariant is satisfied.
93    pub fn push<T: CoordinateLike<X> + IntensityMeasurement>(&mut self, pt: &T, time: f64) {
94        self.push_raw(pt.coordinate(), time, pt.intensity())
95    }
96
97    /// Push a new data point onto the feature and ensure the time ordering invariant is satisfied.
98    pub fn push_raw(&mut self, _x: f64, y: f64, z: f32) {
99        let needs_sort = self.len() > 0 && self.y.last().copied().unwrap() > y;
100        if !self.is_empty() && y == *self.y.last().unwrap() {
101            *self.z.last_mut().unwrap() += z;
102        } else {
103            self.y.push(y);
104            self.z.push(z);
105        }
106        if needs_sort {
107            self.sort_by_y();
108        }
109    }
110
111    fn apex_y(&self) -> Option<f64> {
112        self.apex_of(&self.y, &self.z)
113    }
114
115    fn find_y(&self, y: f64) -> (Option<usize>, f64) {
116        if self.is_empty() {
117            return (None, y);
118        }
119        match self.y.binary_search_by(|yi| y.total_cmp(yi).reverse()) {
120            Ok(i) => {
121                let low = i.saturating_sub(3);
122                (low..(low + 6).min(self.len()))
123                    .map(|i| (Some(i), (self.y[i] - y).abs()))
124                    .min_by(|(_, e), (_, d)| e.total_cmp(d))
125                    .unwrap()
126            }
127            Err(i) => {
128                let low = i.saturating_sub(3);
129                (low..(low + 6).min(self.len()))
130                    .map(|i| (Some(i), (self.y[i] - y).abs()))
131                    .min_by(|(_, e), (_, d)| e.total_cmp(d))
132                    .unwrap()
133            }
134        }
135    }
136}
137
138impl<X, Y> TimeArray<Y> for SimpleFeature<X, Y> {
139    fn time_view(&self) -> &[f64] {
140        &self.y
141    }
142
143    fn intensity_view(&self) -> &[f32] {
144        &self.z
145    }
146}
147
148impl<X, Y> PartialEq for SimpleFeature<X, Y> {
149    fn eq(&self, other: &Self) -> bool {
150        self.label == other.label
151            && self.y == other.y
152            && self.z == other.z
153            && self._x == other._x
154            && self._y == other._y
155    }
156}
157
158impl<X, Y> PartialOrd for SimpleFeature<X, Y> {
159    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
160        match self.label.partial_cmp(&other.label) {
161            Some(core::cmp::Ordering::Equal) => {}
162            ord => return ord,
163        }
164        match self.start_time().partial_cmp(&other.start_time()) {
165            Some(core::cmp::Ordering::Equal) => {}
166            ord => return ord,
167        }
168        self.z.partial_cmp(&other.z)
169    }
170}
171
172impl<X, Y> CoordinateLike<X> for SimpleFeature<X, Y> {
173    fn coordinate(&self) -> f64 {
174        self.label
175    }
176}
177
178impl<X, Y> TimeInterval<Y> for SimpleFeature<X, Y> {
179    fn start_time(&self) -> Option<f64> {
180        self.y.first().copied()
181    }
182
183    fn end_time(&self) -> Option<f64> {
184        self.y.last().copied()
185    }
186
187    fn apex_time(&self) -> Option<f64> {
188        self.apex_y()
189    }
190
191    fn area(&self) -> f32 {
192        let mut it = self.y.iter().zip(self.z.iter());
193        if let Some((first_y, first_z)) = it.next() {
194            let (_y, _z, acc) =
195                it.fold((first_y, first_z, 0.0), |(last_y, last_z, acc), (y, z)| {
196                    let step = (last_z + z) / 2.0;
197                    let dy = y - last_y;
198                    (y, z, acc + (step as f64 * dy))
199                });
200            acc as f32
201        } else {
202            0.0
203        }
204    }
205
206    fn iter_time(&self) -> impl Iterator<Item = f64> {
207        self.y.iter().copied()
208    }
209
210    fn find_time(&self, time: f64) -> (Option<usize>, f64) {
211        self.find_y(time)
212    }
213}
214
215impl<X, Y> IntensityMeasurement for SimpleFeature<X, Y> {
216    fn intensity(&self) -> f32 {
217        self.z.iter().sum()
218    }
219}
220
221impl<X, Y> FeatureLike<X, Y> for SimpleFeature<X, Y> {
222    fn len(&self) -> usize {
223        self.y.len()
224    }
225
226    fn iter(&self) -> impl Iterator<Item = (f64, f64, f32)> {
227        self.y
228            .iter()
229            .zip(self.z.iter())
230            .map(|(y, z)| (self.label, *y, *z))
231    }
232}
233
234pub struct IntoIter<X, Y> {
235    xval: f64,
236    yiter: std::vec::IntoIter<f64>,
237    ziter: std::vec::IntoIter<f32>,
238    _x: PhantomData<X>,
239    _y: PhantomData<Y>,
240}
241
242impl<X, Y> Iterator for IntoIter<X, Y> {
243    type Item = (f64, f64, f32);
244
245    fn next(&mut self) -> Option<Self::Item> {
246        let y = self.yiter.next();
247        let z = self.ziter.next();
248        match (y, z) {
249            (Some(y), Some(z)) => Some((self.xval, y, z)),
250            _ => None,
251        }
252    }
253}
254
255impl<X, Y> IntoIterator for SimpleFeature<X, Y> {
256    type Item = (f64, f64, f32);
257
258    type IntoIter = IntoIter<X, Y>;
259
260    fn into_iter(self) -> Self::IntoIter {
261        let xval = self.coordinate();
262        IntoIter {
263            xval,
264            yiter: self.y.into_iter(),
265            ziter: self.z.into_iter(),
266            _x: PhantomData::<X>,
267            _y: PhantomData::<Y>,
268        }
269    }
270}
271
272impl<X, Y> Extend<(f64, f64, f32)> for SimpleFeature<X, Y> {
273    fn extend<T: IntoIterator<Item = (f64, f64, f32)>>(&mut self, iter: T) {
274        for (_x, y, z) in iter {
275            self.y.push(y);
276            self.z.push(z);
277        }
278    }
279}
280
281impl<X, Y, P: CoordinateLike<X> + IntensityMeasurement> Extend<(P, f64)> for SimpleFeature<X, Y> {
282    fn extend<T: IntoIterator<Item = (P, f64)>>(&mut self, iter: T) {
283        for (x, t) in iter {
284            self.push(&x, t)
285        }
286    }
287}
288
289impl<X, Y, P: CoordinateLike<X> + IntensityMeasurement> FromIterator<(P, f64)>
290    for SimpleFeature<X, Y>
291{
292    fn from_iter<T: IntoIterator<Item = (P, f64)>>(iter: T) -> Self {
293        let mut it = iter.into_iter();
294        if let Some((peak, time)) = it.next() {
295            let mut this = Self::empty(peak.coordinate());
296            this.push(&peak, time);
297            this.extend(it);
298            this
299        } else {
300            Self::empty(0.0)
301        }
302    }
303}
304
305impl<X, Y> FromIterator<(f64, f64, f32)> for SimpleFeature<X, Y> {
306    fn from_iter<T: IntoIterator<Item = (f64, f64, f32)>>(iter: T) -> Self {
307        let mut it = iter.into_iter();
308        if let Some((x, time, intensity)) = it.next() {
309            let mut this = Self::empty(x);
310            this.push_raw(x, time, intensity);
311            this.extend(it);
312            this
313        } else {
314            Self::empty(0.0)
315        }
316    }
317}
318
319#[derive(Debug)]
320pub struct IterMut<'a, X, Y> {
321    yiter: std::slice::IterMut<'a, f64>,
322    ziter: std::slice::IterMut<'a, f32>,
323    value: f64,
324    _x: PhantomData<X>,
325    _y: PhantomData<Y>,
326}
327
328impl<'a, X, Y> Iterator for IterMut<'a, X, Y> {
329    type Item = (&'a mut f64, &'a mut f64, &'a mut f32);
330
331    fn next(&mut self) -> Option<Self::Item> {
332        let x = std::ptr::addr_of_mut!(self.value);
333        let y = self.yiter.next();
334        let z = self.ziter.next();
335
336        match (y, z) {
337            (Some(y), Some(z)) => Some((unsafe { &mut *x }, y, z)),
338            _ => None,
339        }
340    }
341}
342
343impl<X, Y> FeatureLikeMut<X, Y> for SimpleFeature<X, Y> {
344    fn iter_mut(&mut self) -> impl Iterator<Item = (&mut f64, &mut f64, &mut f32)> {
345        IterMut {
346            yiter: self.y.iter_mut(),
347            ziter: self.z.iter_mut(),
348            value: self.label,
349            _x: PhantomData::<X>,
350            _y: PhantomData::<Y>,
351        }
352    }
353
354    fn push<T: CoordinateLike<X> + IntensityMeasurement>(&mut self, pt: &T, time: f64) {
355        self.push(pt, time);
356    }
357
358    fn push_raw(&mut self, x: f64, y: f64, z: f32) {
359        self.push_raw(x, y, z);
360    }
361
362    fn clear(&mut self) {
363        self.y.clear();
364        self.z.clear();
365    }
366
367    fn reserve(&mut self, capacity: usize) {
368        self.y.reserve(capacity);
369        self.z.reserve(capacity);
370    }
371}
372
373/// A non-owning version of [`SimpleFeature`]
374#[derive(Debug, Clone, Copy)]
375#[cfg_attr(feature = "serde", derive(serde::Serialize))]
376pub struct SimpleFeatureView<'a, X, Y> {
377    pub label: f64,
378    y: &'a [f64],
379    z: &'a [f32],
380    _x: PhantomData<X>,
381    _y: PhantomData<Y>,
382}
383
384impl<X, Y> CoArrayOps for SimpleFeatureView<'_, X, Y> {}
385
386impl<'a, X, Y> SimpleFeatureView<'a, X, Y> {
387    pub fn new(label: f64, y: &'a [f64], z: &'a [f32]) -> Self {
388        Self {
389            label,
390            y,
391            z,
392            _x: PhantomData,
393            _y: PhantomData,
394        }
395    }
396
397    pub fn empty(label: f64) -> Self {
398        Self {
399            label,
400            y: EMPTY_Y,
401            z: EMPTY_Z,
402            _x: PhantomData,
403            _y: PhantomData,
404        }
405    }
406
407    fn apex_y(&self) -> Option<f64> {
408        self.apex_of(self.y, self.z)
409    }
410
411    pub fn to_owned(&self) -> SimpleFeature<X, Y> {
412        SimpleFeature::new(self.label, self.y.to_vec(), self.z.to_vec())
413    }
414
415    pub fn into_inner(self) -> (f64, &'a [f64], &'a [f32]) {
416        (self.label, self.y, self.z)
417    }
418
419    fn find_y(&self, y: f64) -> (Option<usize>, f64) {
420        if self.is_empty() {
421            return (None, y);
422        }
423        match self.y.binary_search_by(|yi| y.total_cmp(yi).reverse()) {
424            Ok(i) => {
425                let low = i.saturating_sub(3);
426                (low..(low + 6).min(self.len()))
427                    .map(|i| (Some(i), (self.y[i] - y).abs()))
428                    .min_by(|(_, e), (_, d)| e.total_cmp(d))
429                    .unwrap()
430            }
431            Err(i) => {
432                let low = i.saturating_sub(3);
433                (low..(low + 6).min(self.len()))
434                    .map(|i| (Some(i), (self.y[i] - y).abs()))
435                    .min_by(|(_, e), (_, d)| e.total_cmp(d))
436                    .unwrap()
437            }
438        }
439    }
440}
441
442impl<X, Y> PartialEq for SimpleFeatureView<'_, X, Y> {
443    fn eq(&self, other: &Self) -> bool {
444        self.label == other.label
445            && self.y == other.y
446            && self.z == other.z
447            && self._x == other._x
448            && self._y == other._y
449    }
450}
451
452impl<X, Y> PartialOrd for SimpleFeatureView<'_, X, Y> {
453    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
454        match self.label.partial_cmp(&other.label) {
455            Some(core::cmp::Ordering::Equal) => {}
456            ord => return ord,
457        }
458        match self.start_time().partial_cmp(&other.start_time()) {
459            Some(core::cmp::Ordering::Equal) => {}
460            ord => return ord,
461        }
462        self.z.partial_cmp(other.z)
463    }
464}
465
466impl<X, Y> CoordinateLike<X> for SimpleFeatureView<'_, X, Y> {
467    fn coordinate(&self) -> f64 {
468        self.label
469    }
470}
471
472impl<X, Y> TimeInterval<Y> for SimpleFeatureView<'_, X, Y> {
473    fn start_time(&self) -> Option<f64> {
474        self.y.first().copied()
475    }
476
477    fn end_time(&self) -> Option<f64> {
478        self.y.last().copied()
479    }
480
481    fn apex_time(&self) -> Option<f64> {
482        self.apex_y()
483    }
484
485    fn area(&self) -> f32 {
486        let mut it = self.y.iter().zip(self.z.iter());
487        if let Some((first_y, first_z)) = it.next() {
488            let (_y, _z, acc) =
489                it.fold((first_y, first_z, 0.0), |(last_y, last_z, acc), (y, z)| {
490                    let step = (last_z + z) / 2.0;
491                    let dy = y - last_y;
492                    (y, z, acc + (step as f64 * dy))
493                });
494            acc as f32
495        } else {
496            0.0
497        }
498    }
499
500    fn iter_time(&self) -> impl Iterator<Item = f64> {
501        self.y.iter().copied()
502    }
503
504    fn find_time(&self, time: f64) -> (Option<usize>, f64) {
505        self.find_y(time)
506    }
507}
508
509impl<X, Y> IntensityMeasurement for SimpleFeatureView<'_, X, Y> {
510    fn intensity(&self) -> f32 {
511        self.z.iter().sum()
512    }
513}
514
515impl<X, Y> FeatureLike<X, Y> for SimpleFeatureView<'_, X, Y> {
516    fn len(&self) -> usize {
517        self.y.len()
518    }
519
520    fn iter(&self) -> impl Iterator<Item = (f64, f64, f32)> {
521        self.y
522            .iter()
523            .zip(self.z.iter())
524            .map(|(y, z)| (self.label, *y, *z))
525    }
526}
527
528impl<'a, X, Y> SplittableFeatureLike<'a, X, Y> for SimpleFeature<X, Y> {
529    type ViewType = SimpleFeatureView<'a, X, Y>;
530
531    fn split_at(&'a self, index: usize) -> (Self::ViewType, Self::ViewType) {
532        if !self.is_empty() {
533            let before = Self::ViewType::new(self.label, &self.y[..index], &self.z[..index]);
534            let after = Self::ViewType::new(self.label, &self.y[index..], &self.z[index..]);
535            (before, after)
536        } else {
537            let before = Self::ViewType::new(self.label, EMPTY_Y, EMPTY_Z);
538            let after = Self::ViewType::new(self.label, EMPTY_Y, EMPTY_Z);
539            (before, after)
540        }
541    }
542
543    fn split_at_time(&'a self, point: f64) -> (Self::ViewType, Self::ViewType) {
544        if let Some(point) = self.find_time(point).0 {
545            self.split_at(point)
546        } else {
547            let before = Self::ViewType::new(self.label, EMPTY_Y, EMPTY_Z);
548            let after = Self::ViewType::new(self.label, EMPTY_Y, EMPTY_Z);
549            (before, after)
550        }
551    }
552
553    fn slice<I: RangeBounds<usize> + Clone>(&'a self, bounds: I) -> Self::ViewType {
554        let start = match bounds.start_bound() {
555            Bound::Included(i) => *i,
556            Bound::Excluded(i) => *i,
557            Bound::Unbounded => 0,
558        };
559        let end = match bounds.end_bound() {
560            Bound::Included(i) => *i + 1,
561            Bound::Excluded(i) => *i,
562            Bound::Unbounded => self.len(),
563        };
564        let y = &self.y[start..end];
565        let z = &self.z[start..end];
566        Self::ViewType::new(self.label, y, z)
567    }
568}
569
570impl<'a, X, Y> SplittableFeatureLike<'a, X, Y> for SimpleFeatureView<'a, X, Y> {
571    type ViewType = Self;
572
573    fn split_at_time(&'a self, point: f64) -> (Self::ViewType, Self::ViewType) {
574        if let Some(point) = self.find_time(point).0 {
575            self.split_at(point)
576        } else {
577            let before = Self::ViewType::new(self.label, EMPTY_Y, EMPTY_Z);
578            let after = Self::ViewType::new(self.label, EMPTY_Y, EMPTY_Z);
579            (before, after)
580        }
581    }
582
583    fn slice<I: RangeBounds<usize> + Clone>(&'a self, bounds: I) -> Self::ViewType {
584        let start = bounds.start_bound();
585        let end = bounds.end_bound();
586        match (start, end) {
587            (Bound::Included(i), Bound::Included(j)) => {
588                Self::ViewType::new(self.label, &self.y[*i..=*j], &self.z[*i..=*j])
589            }
590            (Bound::Included(i), Bound::Excluded(j)) => {
591                Self::ViewType::new(self.label, &self.y[*i..*j], &self.z[*i..*j])
592            }
593            (Bound::Included(i), Bound::Unbounded) => {
594                Self::ViewType::new(self.label, &self.y[*i..], &self.z[*i..])
595            }
596            (Bound::Excluded(i), Bound::Included(j)) => {
597                Self::ViewType::new(self.label, &self.y[*i..=*j], &self.z[*i..=*j])
598            }
599            (Bound::Excluded(i), Bound::Excluded(j)) => {
600                Self::ViewType::new(self.label, &self.y[*i..*j], &self.z[*i..*j])
601            }
602            (Bound::Excluded(i), Bound::Unbounded) => {
603                Self::ViewType::new(self.label, &self.y[*i..], &self.z[*i..])
604            }
605            (Bound::Unbounded, Bound::Included(j)) => {
606                Self::ViewType::new(self.label, &self.y[..=*j], &self.z[..=*j])
607            }
608            (Bound::Unbounded, Bound::Excluded(j)) => {
609                Self::ViewType::new(self.label, &self.y[..*j], &self.z[..*j])
610            }
611            (Bound::Unbounded, Bound::Unbounded) => Self::ViewType::new(self.label, self.y, self.z),
612        }
613    }
614
615    fn split_at(&'a self, index: usize) -> (Self::ViewType, Self::ViewType) {
616        if !self.is_empty() {
617            let before = Self::ViewType::new(self.label, &self.y[..index], &self.z[..index]);
618            let after = Self::ViewType::new(self.label, &self.y[index..], &self.z[index..]);
619            (before, after)
620        } else {
621            let before = Self::ViewType::new(self.label, EMPTY_Y, EMPTY_Z);
622            let after = Self::ViewType::new(self.label, EMPTY_Y, EMPTY_Z);
623            (before, after)
624        }
625    }
626}
627
628impl<X, Y> TimeArray<Y> for SimpleFeatureView<'_, X, Y> {
629    fn time_view(&self) -> &[f64] {
630        self.y
631    }
632
633    fn intensity_view(&self) -> &[f32] {
634        self.z
635    }
636}
mzpeaks/feature/simple.rs

mzpeaks/feature/
simple.rs
