Christmas eve checkpoint

Well, it's not pretty, but it's Christmas eve and it's time to ship.

1 files changed, 145 insertions, 0 deletions

diff --git a/src/maxmin.rs b/src/maxmin.rs
new file mode 100644
index 0000000..6bea2d7
--- /dev/null
+++ b/src/maxmin.rs
@@ -0,0 +1,145 @@
+use std::cmp::Ordering;
+
+pub fn normalize(data: &mut Vec<f64>) {
+    if data.is_empty() {
+        return;
+    }
+
+    let max = *data
+        .iter()
+        .max_by(|a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal))
+        .unwrap();
+    let min = *data
+        .iter()
+        .min_by(|a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal))
+        .unwrap();
+    for x in data.iter_mut() {
+        *x = (*x - min) / (max - min);
+    }
+}
+
+pub fn smooth(data: &mut Vec<f64>, sigma: f64) {
+    const WINDOW_SIZE: usize = 7;
+
+    if data.is_empty() {
+        return;
+    }
+
+    // Build Gaussian filter
+    let sigsig = 2.0 * sigma.powi(2);
+    let mut g: Vec<f64> = (0..WINDOW_SIZE)
+        .map(|i| {
+            let x2 = (i as f64).powi(2);
+            std::f64::consts::E.powf(-x2 / sigsig)
+        })
+        .collect();
+
+    // Normalize Gaussian filter
+    let filter_sum: f64 = g.iter().sum();
+    g.iter_mut().for_each(|x| *x /= filter_sum);
+
+    // Smooth data
+    let first = data.first().unwrap();
+    let mut buffer = [*first; WINDOW_SIZE];
+    let mut bi = 0;
+    for (_loc, x) in data.iter_mut().enumerate() {
+        buffer[bi] = *x;
+        *x = 0.0;
+        for (i, b) in g.iter().enumerate() {
+            *x += b * buffer[(bi + i) % WINDOW_SIZE];
+        }
+        bi = (bi + 1) % WINDOW_SIZE;
+    }
+}
+
+pub fn derive(data: &mut Vec<f64>) {
+    const WINDOW_SIZE: usize = 3;
+
+    if data.is_empty() {
+        return;
+    }
+
+    // Calculate the derivative.
+    let first = data.first().unwrap();
+    let mut buffer = [*first; WINDOW_SIZE];
+    let mut bi = 0;
+    for (_loc, x) in data.iter_mut().enumerate() {
+        buffer[bi] = *x;
+        *x = (buffer[bi] - buffer[(bi + WINDOW_SIZE - 1) % (WINDOW_SIZE)])
+            / (WINDOW_SIZE as f64 - 1.0);
+        bi = (bi + 1) % WINDOW_SIZE;
+    }
+}
+
+enum State {
+    Start,
+    NoClue(usize, f64),
+    SeekingMinimum(usize, f64),
+    SeekingMaximum(usize, f64),
+}
+
+#[derive(Debug)]
+pub enum Extremum {
+    Maximum(usize),
+    Minimum(usize),
+}
+
+pub fn find_extrema(data: &[f64]) -> Vec<Extremum> {
+    let mut result = Vec::new();
+    let mut state = State::Start;
+    for (i, v) in data.iter().enumerate() {
+        let datum = *v;
+        state = match state {
+            State::Start => State::NoClue(i, datum),
+            State::NoClue(pi, prior) => {
+                if datum > prior {
+                    result.push(Extremum::Minimum(pi));
+                    State::SeekingMaximum(i, datum)
+                } else {
+                    result.push(Extremum::Maximum(pi));
+                    State::SeekingMinimum(i, datum)
+                }
+            }
+            State::SeekingMinimum(pi, prior) => {
+                if datum <= prior {
+                    State::SeekingMinimum(i, datum)
+                } else {
+                    result.push(Extremum::Minimum(pi));
+                    State::SeekingMaximum(i, datum)
+                }
+            }
+            State::SeekingMaximum(pi, prior) => {
+                if datum >= prior {
+                    State::SeekingMaximum(i, datum)
+                } else {
+                    result.push(Extremum::Maximum(pi));
+                    State::SeekingMinimum(i, datum)
+                }
+            }
+        };
+    }
+
+    match state {
+        State::Start => (),
+        State::NoClue(..) => (),
+        State::SeekingMinimum(i, _) => result.push(Extremum::Minimum(i)),
+        State::SeekingMaximum(i, _) => result.push(Extremum::Maximum(i)),
+    }
+
+    result
+}
+
+#[cfg(test)]
+mod tests {
+    #[test]
+    fn normalize_normalizes() {
+        let mut x = vec![1.0, 2.0, 3.0, 4.0];
+
+        super::normalize(&mut x);
+
+        assert!((0.0 - x[0]).abs() < 0.001);
+        assert!((0.333 - x[1]).abs() < 0.001);
+        assert!((0.667 - x[2]).abs() < 0.001);
+        assert!((1.0 - x[3]).abs() < 0.001);
+    }
+}