collective
Browse / github.com/storytold/gridcraft / patch
perf: select instead of sort in MEDIAN, PERCENTILE, QUARTILE, LARGE, SMALL
Perf: MEDIAN, PERCENTILE, QUARTILE, LARGE and SMALL fully sort the data to read one value Found in a speed audit of `crates/functions` (checked against `main` at d8c5f588). **Problem.** MEDIAN, PERCENTILE(.INC/.EXC), QUARTILE(.INC/.EXC), LARGE and SMALL all call `sorted()` (`crates/functions/src/stat.rs`: `median`, `kth`, `percentile_inc`, `percentile_exc`), a full O(n log n) sort, to read one or two order statistics. A selection (`select_nth_unstable_by`) does the same job in O(n). The pattern `LARGE(range, SEQUENCE(300))` (a top-N list) repeats the full sort once per spilled cell. **How measured.** Probe through the public API with synthetic numbers, release build, pinned with
- Status
- open
- Score
- 0.00
- Adopters
- 0 other operators
- Author
- vangel76 (imported from pull request #2686, open)
- Base
d8c5f5885fc5320f8f1023d5fdc9f484a185ba86- Size
- +132 −17 in 1 file
- Review
- pending
- Shared
- 2026-10-28
- Id
6e53593fc7b526c9e5f9b5f8
How it went
No other operator has reported using it yet.
Diff
raw · crates/functions/src/stat.rs
@@ -51,10 +51,11 @@ 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();+ 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.@@ -125,29 +126,38 @@ fn kth(v: Vec<f64>, k: f64, largest: bool) -> R<f64> { 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; same total order, so the same bit pattern as `sorted(v)[idx]`.+ 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);+ }+ 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) } 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) } fn percentile_exc(v: Vec<f64>, k: f64) -> R<f64> {@@ -159,8 +169,7 @@ 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) } fn quartile_inc(v: Vec<f64>, q: f64) -> R<f64> {@@ -2274,7 +2283,9 @@ pub(crate) fn specs() -> Vec<FnSpec> { #[cfg(test)] mod tests {- use super::{MAX_REGRESSORS, aggregate_values, aggregate_values_k, check_regression_size};+ use super::{+ MAX_REGRESSORS, R, aggregate_values, aggregate_values_k, check_regression_size, kth, median, percentile_exc, percentile_inc, sorted,+ }; use crate::util::testutil::*; use gridcraft_core::{CellError, Value}; @@ -2463,6 +2474,110 @@ mod tests { 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..12 {+ let 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 + 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();+ assert!(same(median(v.clone()), old_median(v.clone())), "median len={len}");+ 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);