Merge pull request #2686 from vangel76/perf/stat-select-2679

perf: select instead of sort in MEDIAN, PERCENTILE, QUARTILE, LARGE, SMALL

crates/functions/src/stat/descriptive.rs

@@ -33,19 +33,40 @@ pub(super) fn var_p(x: &[f64]) -> R<f64> {
     Ok(devsq(x) / x.len() as f64)
 }
 
-fn sorted(mut v: Vec<f64>) -> Vec<f64> {
+pub(super) fn sorted(mut v: Vec<f64>) -> Vec<f64> {
     v.sort_by(f64::total_cmp);
     v
 }
 
+/// Sorts `v` if it is already monotone in `total_cmp` order (a descending slice is reversed) and
+/// says whether it did. At most two passes, each stopping at the first out-of-order pair, so
+/// unordered data costs next to nothing. Needed because the std sort recognises ascending and
+/// descending runs in O(n) with a small constant; selection is O(n) too but measured 1.2-1.8x
+/// slower on such input.
+fn sorted_if_monotone(v: &mut [f64]) -> bool {
+    if v.is_sorted_by(|a, b| a.total_cmp(b).is_le()) {
+        return true;
+    }
+    if v.is_sorted_by(|a, b| a.total_cmp(b).is_ge()) {
+        v.reverse();
+        return true;
+    }
+    false
+}
+
 pub(super) fn median(v: Vec<f64>) -> R<f64> {
     if v.is_empty() {
         return Err(CellError::Num);
     }
-    let s = sorted(v);
-    let n = s.len();
-    let hi = s.get(n / 2).copied().unwrap_or(0.0);
-    if n % 2 == 1 { Ok(hi) } else { Ok(s.get(n / 2 - 1).copied().unwrap_or(0.0).midpoint(hi)) }
+    let mut v = v;
+    let n = v.len();
+    if sorted_if_monotone(&mut v) {
+        let hi = v.get(n / 2).copied().unwrap_or(0.0);
+        return if n % 2 == 1 { Ok(hi) } else { Ok(v.get(n / 2 - 1).copied().unwrap_or(0.0).midpoint(hi)) };
+    }
+    let (lower, hi, _) = v.select_nth_unstable_by(n / 2, f64::total_cmp);
+    let hi = *hi;
+    if n % 2 == 1 { Ok(hi) } else { Ok(lower.iter().copied().max_by(f64::total_cmp).unwrap_or(0.0).midpoint(hi)) }
 }
 
 /// Values occurring at least twice, ordered by first occurrence, with their (equal, maximal) count.
@@ -126,34 +147,54 @@ pub(super) fn max_of(v: &[f64]) -> f64 {
 }
 
 /// k-th largest (k is rounded up like Excel).
-fn kth(v: Vec<f64>, k: f64, largest: bool) -> R<f64> {
+pub(super) fn kth(v: Vec<f64>, k: f64, largest: bool) -> R<f64> {
     let k = k.ceil();
     if k.is_nan() || k < 1.0 || k > v.len() as f64 {
         return Err(CellError::Num);
     }
-    let s = sorted(v);
+    let mut v = v;
     let k = k as usize;
-    let idx = if largest { s.len() - k } else { k - 1 };
-    s.get(idx).copied().ok_or(CellError::Num)
+    let idx = if largest { v.len() - k } else { k - 1 };
+    // Selection instead of a full sort (monotone input is read directly); same total order, so the
+    // same bit pattern as `sorted(v)[idx]`.
+    if sorted_if_monotone(&mut v) {
+        return v.get(idx).copied().ok_or(CellError::Num);
+    }
+    Ok(*v.select_nth_unstable_by(idx, f64::total_cmp).1)
 }
 
-fn interp(s: &[f64], pos: f64) -> R<f64> {
+/// Value at fractional rank `pos` of the sorted data, linearly interpolated, found by selection.
+fn interp(mut v: Vec<f64>, pos: f64) -> R<f64> {
     let i = pos.floor();
     let frac = pos - i;
     let i = i as usize;
-    let lo = *s.get(i).ok_or(CellError::Num)?;
-    match s.get(i + 1) {
-        Some(hi) if frac > 0.0 => Ok(lo + frac * (hi - lo)),
-        _ => Ok(lo),
+    if i >= v.len() {
+        return Err(CellError::Num);
     }
+    if sorted_if_monotone(&mut v) {
+        let lo = v.get(i).copied().ok_or(CellError::Num)?;
+        return match v.get(i + 1) {
+            Some(hi) if frac > 0.0 => Ok(lo + frac * (hi - lo)),
+            _ => Ok(lo),
+        };
+    }
+    let (_, lo, upper) = v.select_nth_unstable_by(i, f64::total_cmp);
+    let lo = *lo;
+    // The next order statistic is the minimum of the upper partition.
+    if frac > 0.0
+        && let Some(hi) = upper.iter().copied().min_by(f64::total_cmp)
+    {
+        return Ok(lo + frac * (hi - lo));
+    }
+    Ok(lo)
 }
 
 pub(super) fn percentile_inc(v: Vec<f64>, k: f64) -> R<f64> {
     if v.is_empty() || !(0.0..=1.0).contains(&k) {
         return Err(CellError::Num);
     }
-    let s = sorted(v);
-    interp(&s, k * (s.len() - 1) as f64)
+    let pos = k * (v.len() - 1) as f64;
+    interp(v, pos)
 }
 
 pub(super) fn percentile_exc(v: Vec<f64>, k: f64) -> R<f64> {
@@ -165,8 +206,7 @@ pub(super) fn percentile_exc(v: Vec<f64>, k: f64) -> R<f64> {
     if pos < 1.0 || pos > n {
         return Err(CellError::Num);
     }
-    let s = sorted(v);
-    interp(&s, pos - 1.0)
+    interp(v, pos - 1.0)
 }
 
 pub(super) fn quartile_inc(v: Vec<f64>, q: f64) -> R<f64> {

crates/functions/src/stat/tests.rs

@@ -1,4 +1,6 @@
+use super::descriptive::{kth, median, percentile_exc, percentile_inc, sorted};
 use super::{MAX_REGRESSORS, aggregate_values, aggregate_values_k, check_regression_size};
+use crate::util::R;
 use crate::util::testutil::*;
 use gridcraft_core::{CellError, Value};
 
@@ -200,6 +202,126 @@ fn min_max_order() {
     close(ev("RANK.AVG", vec![n(3.5), r2, n(1.0)]), 3.5);
 }
 
+// Reference implementations: the previous full-sort versions, as the oracle for the selection ones.
+fn old_interp(s: &[f64], pos: f64) -> R<f64> {
+    let i = pos.floor();
+    let frac = pos - i;
+    let i = i as usize;
+    let lo = *s.get(i).ok_or(CellError::Num)?;
+    match s.get(i + 1) {
+        Some(hi) if frac > 0.0 => Ok(lo + frac * (hi - lo)),
+        _ => Ok(lo),
+    }
+}
+
+fn old_median(v: Vec<f64>) -> R<f64> {
+    if v.is_empty() {
+        return Err(CellError::Num);
+    }
+    let s = sorted(v);
+    let n = s.len();
+    let hi = s.get(n / 2).copied().unwrap_or(0.0);
+    if n % 2 == 1 { Ok(hi) } else { Ok(s.get(n / 2 - 1).copied().unwrap_or(0.0).midpoint(hi)) }
+}
+
+fn old_kth(v: Vec<f64>, k: f64, largest: bool) -> R<f64> {
+    let k = k.ceil();
+    if k.is_nan() || k < 1.0 || k > v.len() as f64 {
+        return Err(CellError::Num);
+    }
+    let s = sorted(v);
+    let k = k as usize;
+    let idx = if largest { s.len() - k } else { k - 1 };
+    s.get(idx).copied().ok_or(CellError::Num)
+}
+
+fn old_percentile_inc(v: Vec<f64>, k: f64) -> R<f64> {
+    if v.is_empty() || !(0.0..=1.0).contains(&k) {
+        return Err(CellError::Num);
+    }
+    let s = sorted(v);
+    old_interp(&s, k * (s.len() - 1) as f64)
+}
+
+fn old_percentile_exc(v: Vec<f64>, k: f64) -> R<f64> {
+    let n = v.len() as f64;
+    if v.is_empty() || !(k > 0.0 && k < 1.0) {
+        return Err(CellError::Num);
+    }
+    let pos = k * (n + 1.0);
+    if pos < 1.0 || pos > n {
+        return Err(CellError::Num);
+    }
+    let s = sorted(v);
+    old_interp(&s, pos - 1.0)
+}
+
+/// Same Ok bit pattern or same error.
+fn same(a: R<f64>, b: R<f64>) -> bool {
+    match (a, b) {
+        (Ok(x), Ok(y)) => x.to_bits() == y.to_bits(),
+        (Err(x), Err(y)) => x == y,
+        _ => false,
+    }
+}
+
+#[test]
+fn selection_matches_full_sort_bit_for_bit() {
+    let mut state = 0x9E37_79B9_7F4A_7C15_u64;
+    let mut next = move || {
+        state ^= state << 13;
+        state ^= state >> 7;
+        state ^= state << 17;
+        state
+    };
+    let specials = [0.0, -0.0, f64::INFINITY, f64::NEG_INFINITY, f64::NAN, -f64::NAN, f64::MAX, f64::MIN_POSITIVE, 1.0, -1.0];
+    for len in (0..=40).chain([63, 64, 65, 100, 257, 1000, 4096]) {
+        for round in 0..24 {
+            let mut v: Vec<f64> = (0..len)
+                .map(|_| match (round % 4, next() % 8) {
+                    // Few distinct values: many ties.
+                    (0, r) => (r % 3) as f64,
+                    // Specials, including both zeros and NaNs.
+                    (1, r) => specials[(next() as usize).wrapping_add(r as usize) % specials.len()],
+                    // Wide random range.
+                    (2, _) => (next() as i64 as f64) / 1e3,
+                    // Integers with some specials mixed in.
+                    (_, 0) => specials[next() as usize % specials.len()],
+                    _ => (next() % 50) as f64,
+                })
+                .collect();
+            // Rounds 12..24 repeat the data as already ascending, already descending and nearly
+            // ascending (ends swapped): the monotone shortcut and its fallback to selection.
+            match round / 4 {
+                3 => v.sort_by(f64::total_cmp),
+                4 => {
+                    v.sort_by(f64::total_cmp);
+                    v.reverse();
+                }
+                5 => {
+                    v.sort_by(f64::total_cmp);
+                    if let Some(last) = v.len().checked_sub(1) {
+                        v.swap(0, last);
+                    }
+                }
+                _ => {}
+            }
+            assert!(same(median(v.clone()), old_median(v.clone())), "median len={len} round={round}");
+            let ks = [-1.0, 0.0, 0.2, 1.0, 1.5, 2.0, len as f64 / 2.0, len as f64 - 0.5, len as f64, len as f64 + 1.0, f64::NAN];
+            for &k in &ks {
+                for largest in [true, false] {
+                    assert!(same(kth(v.clone(), k, largest), old_kth(v.clone(), k, largest)), "kth len={len} k={k}");
+                }
+            }
+            let ps = [-0.1, 0.0, 0.01, 0.1, 0.25, 0.3, 0.5, 0.75, 0.9, 0.99, 1.0, 1.1, f64::NAN, (next() % 1001) as f64 / 1000.0];
+            for &p in &ps {
+                assert!(same(percentile_inc(v.clone(), p), old_percentile_inc(v.clone(), p)), "pinc len={len} p={p}");
+                assert!(same(percentile_exc(v.clone(), p), old_percentile_exc(v.clone(), p)), "pexc len={len} p={p}");
+            }
+        }
+    }
+}
+
 #[test]
 fn percentiles() {
     close(ev("PERCENTILE.INC", vec![c(&[1.0, 3.0, 2.0, 4.0]), n(0.3)]), 1.9);