smoothed_z_score/
lib.rs

1#[derive(Clone, Copy, PartialEq, Debug)]
2pub enum Peak {
3    Low,
4    High,
5}
6
7pub struct PeaksDetector {
8    threshold: f64,
9    influence: f64,
10    window: Vec<f64>,
11}
12
13impl PeaksDetector {
14    pub fn new(lag: usize, threshold: f64, influence: f64) -> PeaksDetector {
15        PeaksDetector {
16            threshold,
17            influence,
18            window: Vec::with_capacity(lag),
19        }
20    }
21
22    pub fn signal(&mut self, value: f64) -> Option<Peak> {
23        if self.window.len() < self.window.capacity() {
24            self.window.push(value);
25            None
26        } else if let (Some((mean, stddev)), Some(&window_last)) = (self.stats(), self.window.last()) {
27            self.window.remove(0);
28            if (value - mean).abs() > (self.threshold * stddev) {
29                let next_value =
30                    (value * self.influence) + ((1. - self.influence) * window_last);
31                self.window.push(next_value);
32                Some(if value > mean { Peak::High } else { Peak::Low })
33            } else {
34                self.window.push(value);
35                None
36            }
37        } else {
38            None
39        }
40    }
41
42    pub fn stats(&self) -> Option<(f64, f64)> {
43        if self.window.is_empty() {
44            None
45        } else {
46            let window_len = self.window.len() as f64;
47            let mean = self.window.iter().fold(0., |a, v| a + v) / window_len;
48            let sq_sum = self.window.iter().fold(0., |a, v| a + ((v - mean) * (v - mean)));
49            let stddev = (sq_sum / window_len).sqrt();
50            Some((mean, stddev))
51        }
52    }
53}
54
55pub struct PeaksIter<I, F> {
56    source: I,
57    signal: F,
58    detector: PeaksDetector,
59}
60
61pub trait PeaksFilter<I> where I: Iterator {
62    fn peaks<F>(self, detector: PeaksDetector, signal: F) -> PeaksIter<I, F>
63        where F: FnMut(&I::Item) -> f64;
64}
65
66impl<I> PeaksFilter<I> for I where I: Iterator {
67    fn peaks<F>(self, detector: PeaksDetector, signal: F) -> PeaksIter<I, F>
68        where F: FnMut(&I::Item) -> f64
69    {
70        PeaksIter { source: self, signal, detector, }
71    }
72}
73
74impl<I, F> Iterator for PeaksIter<I, F> where I: Iterator, F: FnMut(&I::Item) -> f64 {
75    type Item = (I::Item, Peak);
76
77    fn next(&mut self) -> Option<Self::Item> {
78        while let Some(item) = self.source.next() {
79            let value = (self.signal)(&item);
80            if let Some(peak) = self.detector.signal(value) {
81                return Some((item, peak));
82            }
83        }
84        None
85    }
86}
87
88#[cfg(test)]
89mod tests {
90    use super::{Peak, PeaksDetector, PeaksFilter};
91
92    #[test]
93    fn sample_data() {
94        let input = vec![
95            1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0,
96            1.0, 1.0, 1.0, 1.1, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 1.1, 1.0, 0.8, 0.9, 1.0,
97            1.2, 0.9, 1.0, 1.0, 1.1, 1.2, 1.0, 1.5, 1.0, 3.0, 2.0, 5.0, 3.0, 2.0, 1.0, 1.0, 1.0, 0.9, 1.0,
98            1.0, 3.0, 2.6, 4.0, 3.0, 3.2, 2.0, 1.0, 1.0, 0.8, 4.0, 4.0, 2.0, 2.5, 1.0, 1.0, 1.0
99        ];
100        let output: Vec<_> = input
101            .into_iter()
102            .enumerate()
103            .peaks(PeaksDetector::new(30, 5.0, 0.0), |e| e.1)
104            .map(|((i, _), p)| (i, p))
105            .collect();
106        assert_eq!(output, vec![
107            (45, Peak::High), (47, Peak::High), (48, Peak::High), (49, Peak::High),
108            (50, Peak::High), (51, Peak::High), (58, Peak::High), (59, Peak::High),
109            (60, Peak::High), (61, Peak::High), (62, Peak::High), (63, Peak::High),
110            (67, Peak::High), (68, Peak::High), (69, Peak::High), (70, Peak::High),
111        ]);
112    }
113}