Split out statistics module

src/benchmark/measurement.rs

@@ -2,10 +2,8 @@ use std::process::ExitStatus;
 
 use serde::Serialize;
 
-use crate::quantity::{
+use crate::quantity::statistics::{max, mean, median, min, modified_zscores, standard_deviation};
-    max, mean, median, min, modified_zscores, serialize_information, serialize_time,
+use crate::quantity::{serialize_information, serialize_time, Information, Time};
-    standard_deviation, Information, Time,
-};
 use crate::util::exit_code::extract_exit_code;
 
 fn serialize_exit_status<S>(exit_status: &ExitStatus, serializer: S) -> Result<S::Ok, S::Error>

src/quantity/mod.rs

@@ -1,13 +1,9 @@
 use core::f64;
 use std::marker::PhantomData;
-use std::ops::Add;
-use std::ops::AddAssign;
-use std::ops::Div;
 
 use serde::ser::SerializeStruct;
 use serde::Serializer;
 
-use uom::num_traits;
 use uom::si;
 
 pub use si::f64::{Information, Ratio, Time};
@@ -18,6 +14,7 @@ pub use uom::num_traits::Zero;
 
 pub use units::{InformationUnit, IsUnit, TimeUnit};
 
+pub mod statistics;
 mod units;
 
 pub trait FormatQuantity {
@@ -83,8 +80,12 @@ impl FormatQuantity for Information {
     fn suitable_unit(&self) -> InformationUnit {
         if *self < Information::new::<kibibyte>(1.0) {
             InformationUnit::Byte
-        } else {
+        } else if *self < Information::new::<mebibyte>(1.0) {
             InformationUnit::KibiByte
+        } else if *self < Information::new::<gibibyte>(1.0) {
+            InformationUnit::MebiByte
+        } else {
+            InformationUnit::GibiByte
         }
     }
 
@@ -132,113 +133,6 @@ where
     state.end()
 }
 
-pub trait UnsafeRawValue {
-    fn unsafe_raw_value(&self) -> f64;
-    fn unsafe_from_raw_value(value: f64) -> Self;
-}
-
-impl UnsafeRawValue for Time {
-    fn unsafe_raw_value(&self) -> f64 {
-        self.get::<second>()
-    }
-
-    fn unsafe_from_raw_value(value: f64) -> Self {
-        Time::new::<second>(value)
-    }
-}
-
-impl UnsafeRawValue for Information {
-    fn unsafe_raw_value(&self) -> f64 {
-        self.get::<byte>()
-    }
-
-    fn unsafe_from_raw_value(value: f64) -> Self {
-        Information::new::<byte>(value)
-    }
-}
-
-macro_rules! quantity_fn {
-    ($name:ident, $unwrapped_values:ident, $body:expr) => {
-        pub fn $name<Q: UnsafeRawValue>(values: &[Q]) -> Q {
-            let $unwrapped_values: Vec<_> = values.iter().map(|q| q.unsafe_raw_value()).collect();
-            let result_value = $body;
-
-            Q::unsafe_from_raw_value(result_value)
-        }
-    };
-}
-
-/// A max function that assumes no NaNs and at least one element
-pub fn max<V: PartialOrd>(values: impl IntoIterator<Item = V>) -> V {
-    values
-        .into_iter()
-        .max_by(|a, b| a.partial_cmp(b).expect("No NaN values"))
-        .expect("'max' requires at least one element")
-}
-
-/// A min function that assumes no NaNs and at least one element
-pub fn min<V: PartialOrd>(values: impl IntoIterator<Item = V>) -> V {
-    values
-        .into_iter()
-        .min_by(|a, b| a.partial_cmp(b).expect("No NaN values"))
-        .expect("'min' requires at least one element")
-}
-
-/// A mean function that assumes at least one element
-pub fn mean<Q, P>(values: impl IntoIterator<Item = Q>) -> Q
-where
-    Q: AddAssign + num_traits::Zero + Div<Ratio, Output = P>,
-    P: Into<Q>,
-{
-    let mut sum = Q::zero();
-    let mut count = 0;
-    for value in values {
-        sum += value;
-        count += 1;
-    }
-
-    let count = Ratio::new::<ratio>(count as f64);
-    (sum / count).into()
-}
-
-pub fn median<Q, P>(values: impl IntoIterator<Item = Q>) -> Q
-where
-    Q: Copy + PartialOrd + Add<Output = Q> + Div<Ratio, Output = P>,
-    P: Into<Q>,
-{
-    let mut values = values.into_iter().collect::<Vec<_>>();
-    values.sort_by(|a, b| a.partial_cmp(b).expect("No NaN values"));
-
-    let len = values.len();
-    if len % 2 == 0 {
-        let mid = len / 2;
-        let a = &values[mid - 1];
-        let b = &values[mid];
-        ((*a + *b) / Ratio::new::<ratio>(2.)).into()
-    } else {
-        values[len / 2]
-    }
-}
-
-quantity_fn!(standard_deviation, values, {
-    let mean_value = statistical::mean(&values);
-    statistical::standard_deviation(&values, Some(mean_value))
-});
-
-pub fn modified_zscores<Q: UnsafeRawValue>(values: &[Q]) -> Vec<f64> {
-    let values: Vec<_> = values.iter().map(|q| q.unsafe_raw_value()).collect();
-    crate::outlier_detection::modified_zscores(&values)
-}
-
-#[test]
-fn test_max() {
-    assert_eq!(1.0, max([1.0]));
-    assert_eq!(-1.0, max([-1.0]));
-    assert_eq!(-1.0, max([-2.0, -1.0]));
-    assert_eq!(1.0, max([-1.0, 1.0]));
-    assert_eq!(1.0, max([-1.0, 1.0, 0.0]));
-}
-
 #[test]
 fn test_time() {
     let time = Time::new::<millisecond>(123.4);
@@ -266,32 +160,7 @@ fn test_information() {
 }
 
 #[test]
-fn test_format() {
+fn test_suiteable_unit_time() {
-    let time = Time::new::<millisecond>(123.4);
-    assert_eq!(time.format(TimeUnit::Second), "0.123 s");
-    assert_eq!(time.format(TimeUnit::MilliSecond), "123.4 ms");
-    assert_eq!(time.format(TimeUnit::MicroSecond), "123400.0 µs");
-
-    let peak_memory_usage = Information::new::<kibibyte>(8.);
-    assert_eq!(peak_memory_usage.format(InformationUnit::Byte), "8192 B");
-    assert_eq!(
-        peak_memory_usage.format(InformationUnit::KibiByte),
-        "8.0 KiB"
-    );
-}
-
-#[test]
-fn test_mean() {
-    let values = vec![
-        Time::new::<millisecond>(123.4),
-        Time::new::<millisecond>(234.5),
-    ];
-    let result = mean(values.into_iter());
-    assert_eq!(result.format(TimeUnit::MilliSecond), "178.9 ms");
-}
-
-#[test]
-fn test_suiteable_unit() {
     assert_eq!(Time::new::<second>(1.3).suitable_unit(), TimeUnit::Second);
     assert_eq!(Time::new::<second>(1.0).suitable_unit(), TimeUnit::Second);
     assert_eq!(
@@ -312,6 +181,41 @@ fn test_suiteable_unit() {
     );
 }
 
+#[test]
+fn test_suitable_unit_information() {
+    assert_eq!(
+        Information::new::<byte>(512.0).suitable_unit(),
+        InformationUnit::Byte
+    );
+    assert_eq!(
+        Information::new::<byte>(2048.0).suitable_unit(),
+        InformationUnit::KibiByte
+    );
+    assert_eq!(
+        Information::new::<byte>(2_097_152.0).suitable_unit(),
+        InformationUnit::MebiByte
+    );
+    assert_eq!(
+        Information::new::<byte>(2_147_483_648.0).suitable_unit(),
+        InformationUnit::GibiByte
+    );
+}
+
+#[test]
+fn test_format() {
+    let time = Time::new::<millisecond>(123.4);
+    assert_eq!(time.format(TimeUnit::Second), "0.123 s");
+    assert_eq!(time.format(TimeUnit::MilliSecond), "123.4 ms");
+    assert_eq!(time.format(TimeUnit::MicroSecond), "123400.0 µs");
+
+    let peak_memory_usage = Information::new::<kibibyte>(8.);
+    assert_eq!(peak_memory_usage.format(InformationUnit::Byte), "8192 B");
+    assert_eq!(
+        peak_memory_usage.format(InformationUnit::KibiByte),
+        "8.0 KiB"
+    );
+}
+
 #[test]
 fn test_format_duration_unit_with_unit() {
     let out = Time::new::<second>(1.3).format(TimeUnit::Second);
@@ -323,34 +227,3 @@ fn test_format_duration_unit_with_unit() {
     let out = Time::new::<second>(1.3).format(TimeUnit::MicroSecond);
     assert_eq!("1300000.0 µs", out);
 }
-
-#[test]
-fn statistics() {
-    let values = vec![
-        Time::new::<second>(1.0),
-        Time::new::<second>(2.0),
-        Time::new::<second>(3.0),
-    ];
-
-    assert_eq!(
-        mean(values.iter().copied()).format(TimeUnit::Second),
-        "2.000 s"
-    );
-    assert_eq!(
-        median(values.iter().copied()).format(TimeUnit::Second),
-        "2.000 s"
-    );
-    assert_eq!(min(&values).format(TimeUnit::Second), "1.000 s");
-    assert_eq!(max(&values).format(TimeUnit::Second), "3.000 s");
-    assert_eq!(
-        standard_deviation(&values).format(TimeUnit::Second),
-        "1.000 s"
-    );
-
-    let values = vec![
-        Information::new::<byte>(1.0),
-        Information::new::<byte>(2.0),
-        Information::new::<byte>(3.0),
-    ];
-    mean(values.iter().copied()).format(InformationUnit::Byte);
-}

src/quantity/statistics.rs

@@ -8,64 +8,6 @@ use uom::si::ratio::ratio;
 
 use super::{byte, second, Information, Time};
 
-/// A min function that assumes no NaNs and at least one element
-pub fn min<Q: PartialOrd>(values: impl IntoIterator<Item = Q>) -> Q {
-    values
-        .into_iter()
-        .min_by(|a, b| a.partial_cmp(b).expect("No NaN values"))
-        .expect("'min' requires at least one element")
-}
-
-/// A max function that assumes no NaNs and at least one element
-pub fn max<Q: PartialOrd>(values: impl IntoIterator<Item = Q>) -> Q {
-    values
-        .into_iter()
-        .max_by(|a, b| a.partial_cmp(b).expect("No NaN values"))
-        .expect("'max' requires at least one element")
-}
-
-/// A mean function that assumes at least one element
-pub fn mean<Q, P>(values: impl IntoIterator<Item = Q>) -> Q
-where
-    Q: AddAssign + num_traits::Zero + Div<Ratio, Output = P>,
-    P: Into<Q>,
-{
-    let mut sum = Q::zero();
-    let mut count = 0;
-    for value in values {
-        sum += value;
-        count += 1;
-    }
-
-    let count = Ratio::new::<ratio>(count as f64);
-    (sum / count).into()
-}
-
-pub fn median<Q, P>(values: impl IntoIterator<Item = Q>) -> Q
-where
-    Q: Copy + PartialOrd + Add<Output = Q> + Div<Ratio, Output = P>,
-    P: Into<Q>,
-{
-    let mut values = values.into_iter().collect::<Vec<_>>();
-    assert!(
-        !values.is_empty(),
-        "median requires at least one measurement"
-    );
-    values.sort_by(|a, b| a.partial_cmp(b).expect("No NaN values"));
-
-    let len = values.len();
-    if len % 2 == 0 {
-        let mid = len / 2;
-        let a = &values[mid - 1];
-        let b = &values[mid];
-        ((*a + *b) / Ratio::new::<ratio>(2.)).into()
-    } else {
-        values[len / 2]
-    }
-}
-
-/// This is just sad, but after several hours trying to figure out how to write
-/// a generic `standard_deviation` function using `uom` quantities, I gave up.
 pub trait UnsafeRawValue {
     fn unsafe_raw_value(&self) -> f64;
     fn unsafe_from_raw_value(value: f64) -> Self;
@@ -91,140 +33,123 @@ impl UnsafeRawValue for Information {
     }
 }
 
-fn standard_deviation_f64(values: &[f64]) -> f64 {
+/// A max function that assumes no NaNs and at least one element
-    let mean_value = mean(values.iter().copied());
+pub fn max<Q: PartialOrd>(values: impl IntoIterator<Item = Q>) -> Q {
+    values
+        .into_iter()
+        .max_by(|a, b| a.partial_cmp(b).expect("No NaN values"))
+        .expect("'max' requires at least one element")
+}
 
-    let mut squared_deviations = 0.;
+/// A min function that assumes no NaNs and at least one element
-    let mut n = 0;
+pub fn min<Q: PartialOrd>(values: impl IntoIterator<Item = Q>) -> Q {
+    values
+        .into_iter()
+        .min_by(|a, b| a.partial_cmp(b).expect("No NaN values"))
+        .expect("'min' requires at least one element")
+}
+
+/// A mean function that assumes at least one element
+pub fn mean<Q, P>(values: impl IntoIterator<Item = Q>) -> Q
+where
+    Q: AddAssign + num_traits::Zero + Div<Ratio, Output = P>,
+    P: Into<Q>,
+{
+    let mut sum = Q::zero();
+    let mut count = 0;
     for value in values {
-        let deviation = value - mean_value;
+        sum += value;
-        squared_deviations += deviation * deviation;
+        count += 1;
-        n += 1;
     }
 
-    (1. / ((n - 1) as f64) * squared_deviations).sqrt()
+    let count = Ratio::new::<ratio>(count as f64);
+    (sum / count).into()
 }
 
-pub fn standard_deviation<Q: UnsafeRawValue>(values: impl IntoIterator<Item = Q> + Clone) -> Q {
+pub fn median<Q, P>(values: impl IntoIterator<Item = Q>) -> Q
-    let values: Vec<_> = values.into_iter().map(|q| q.unsafe_raw_value()).collect();
+where
-    let result = standard_deviation_f64(&values);
+    Q: Copy + PartialOrd + Add<Output = Q> + Div<Ratio, Output = P>,
-    Q::unsafe_from_raw_value(result)
+    P: Into<Q>,
+{
+    let mut values = values.into_iter().collect::<Vec<_>>();
+    values.sort_by(|a, b| a.partial_cmp(b).expect("No NaN values"));
+
+    let len = values.len();
+    if len % 2 == 0 {
+        let mid = len / 2;
+        let a = &values[mid - 1];
+        let b = &values[mid];
+        ((*a + *b) / Ratio::new::<ratio>(2.)).into()
+    } else {
+        values[len / 2]
+    }
 }
 
-/// Compute modified Z-scores for a given sample. A (unmodified) Z-score is defined by
+pub fn standard_deviation<Q: UnsafeRawValue>(values: &[Q]) -> Q {
-/// `(x_i - x_mean)/x_stddev` whereas the modified Z-score is defined by `(x_i - x_median)/MAD`
+    let values: Vec<_> = values.iter().map(|q| q.unsafe_raw_value()).collect();
-/// where MAD is the median absolute deviation.
+    let result_value = {
-///
+        let mean_value = statistical::mean(&values);
-/// References:
+        statistical::standard_deviation(&values, Some(mean_value))
-/// - <https://en.wikipedia.org/wiki/Median_absolute_deviation>
+    };
-pub fn modified_zscores_f64(xs: &[f64]) -> Vec<f64> {
+    Q::unsafe_from_raw_value(result_value)
-    assert!(!xs.is_empty());
-
-    // Compute sample median:
-    let x_median = median(xs.iter().copied());
-
-    // Compute the absolute deviations from the median:
-    let deviations: Vec<f64> = xs.iter().map(|x| (x - x_median).abs()).collect();
-
-    // Compute median absolute deviation:
-    let mad = median(deviations.iter().copied());
-
-    // Handle MAD == 0 case
-    let mad = if mad > 0.0 { mad } else { f64::EPSILON };
-
-    // Compute modified Z-scores (x_i - x_median) / MAD
-    xs.iter().map(|&x| (x - x_median) / mad).collect()
 }
 
 pub fn modified_zscores<Q: UnsafeRawValue>(values: &[Q]) -> Vec<f64> {
     let values: Vec<_> = values.iter().map(|q| q.unsafe_raw_value()).collect();
-    modified_zscores_f64(&values)
+    crate::outlier_detection::modified_zscores(&values)
 }
 
-#[cfg(test)]
+#[test]
-mod tests {
+fn test_max() {
-    use super::*;
+    assert_eq!(1.0, max([1.0]));
-    use approx::assert_relative_eq;
+    assert_eq!(-1.0, max([-1.0]));
-    use uom::si::information::kibibyte;
+    assert_eq!(-1.0, max([-2.0, -1.0]));
-    use uom::si::time::{microsecond, millisecond};
+    assert_eq!(1.0, max([-1.0, 1.0]));
-
+    assert_eq!(1.0, max([-1.0, 1.0, 0.0]));
-    #[test]
+}
-    fn test_min() {
+
-        assert_eq!(1.0, min([1.0]));
+#[test]
-        assert_eq!(-1.0, min([-1.0]));
+fn test_mean() {
-        assert_eq!(-2.0, min([-2.0, -1.0]));
+    use crate::quantity::{FormatQuantity, TimeUnit};
-        assert_eq!(-1.0, min([-1.0, 1.0]));
+    use uom::si::time::millisecond;
-        assert_eq!(-1.0, min([1.0, -1.0, 0.0]));
+
-
+    let values = vec![
-        let values = vec![
+        Time::new::<millisecond>(123.4),
-            Information::new::<kibibyte>(1.0),
+        Time::new::<millisecond>(234.5),
-            Information::new::<byte>(2.0),
+    ];
-            Information::new::<byte>(3.0),
+    let result = mean(values.into_iter());
-        ];
+    assert_eq!(result.format(TimeUnit::MilliSecond), "178.9 ms");
-        assert_eq!(min(&values).get::<byte>(), 2.0);
+}
-    }
+
-
+#[test]
-    #[test]
+fn statistics() {
-    fn test_max() {
+    use crate::quantity::{FormatQuantity, InformationUnit, TimeUnit};
-        assert_eq!(1.0, max([1.0]));
+
-        assert_eq!(-1.0, max([-1.0]));
+    let values = vec![
-        assert_eq!(-1.0, max([-2.0, -1.0]));
+        Time::new::<second>(1.0),
-        assert_eq!(1.0, max([-1.0, 1.0]));
+        Time::new::<second>(2.0),
-        assert_eq!(1.0, max([-1.0, 1.0, 0.0]));
+        Time::new::<second>(3.0),
-
+    ];
-        let values = vec![
+
-            Information::new::<byte>(1.0),
+    assert_eq!(
-            Information::new::<kibibyte>(2.0),
+        mean(values.iter().copied()).format(TimeUnit::Second),
-            Information::new::<byte>(3.0),
+        "2.000 s"
-        ];
+    );
-        assert_eq!(max(&values).get::<kibibyte>(), 2.0);
+    assert_eq!(
-    }
+        median(values.iter().copied()).format(TimeUnit::Second),
-
+        "2.000 s"
-    #[test]
+    );
-    fn test_mean() {
+    assert_eq!(min(&values).format(TimeUnit::Second), "1.000 s");
-        assert_eq!(1.0, mean([1.0]));
+    assert_eq!(max(&values).format(TimeUnit::Second), "3.000 s");
-        assert_relative_eq!(2.0, mean([1.0, 3.0]));
+    assert_eq!(
-
+        standard_deviation(&values).format(TimeUnit::Second),
-        let values = [
+        "1.000 s"
-            Time::new::<millisecond>(100.0),
+    );
-            Time::new::<millisecond>(200.0),
+
-            Time::new::<microsecond>(600_000.0),
+    let values = vec![
-        ];
+        Information::new::<byte>(1.0),
-        let result = mean(values);
+        Information::new::<byte>(2.0),
-        assert_relative_eq!(result.get::<millisecond>(), 300.0);
+        Information::new::<byte>(3.0),
-    }
+    ];
-
+    mean(values.iter().copied()).format(InformationUnit::Byte);
-    #[test]
-    fn test_median() {
-        assert_eq!(1.0, median([1.0]));
-        assert_relative_eq!(2.0, median([1.0, 3.0]));
-
-        let values = [
-            Time::new::<millisecond>(100.0),
-            Time::new::<millisecond>(200.0),
-            Time::new::<microsecond>(600_000.0),
-        ];
-        let result = median(values);
-        assert_relative_eq!(result.get::<millisecond>(), 200.0);
-
-        let values = [
-            Time::new::<millisecond>(100.0),
-            Time::new::<millisecond>(200.0),
-            Time::new::<microsecond>(300_000.0),
-            Time::new::<microsecond>(600_000.0),
-        ];
-        let result = median(values);
-        assert_relative_eq!(result.get::<millisecond>(), 250.0);
-    }
-
-    #[test]
-    fn test_standard_deviation() {
-        let values = [
-            Time::new::<millisecond>(100.0),
-            Time::new::<millisecond>(200.0),
-            Time::new::<microsecond>(300_000.0),
-        ];
-        let result = standard_deviation(values);
-        assert_relative_eq!(result.get::<millisecond>(), 100.0);
-    }
 }