diff options
Diffstat (limited to 'src/maxmin.rs')
-rw-r--r-- | src/maxmin.rs | 145 |
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); + } +} |