From 9d19a8e04d0398bbcbba5e5d106ca0698c2037af Mon Sep 17 00:00:00 2001 From: Dominik Strebel Date: Fri, 6 Nov 2026 13:58:01 +0100 Subject: [PATCH 1/2] perf(data): band statistics count on several threads, summing in the band's order --- crates/data/src/raster/mod.rs | 221 ++++++++++++++++++++++++++++------ 1 file changed, 186 insertions(+), 35 deletions(-) diff --git a/crates/data/src/raster/mod.rs b/crates/data/src/raster/mod.rs index 31dc5efdb..07a516728 100644 --- a/crates/data/src/raster/mod.rs +++ b/crates/data/src/raster/mod.rs @@ -166,26 +166,41 @@ impl BandStats { BandStats::compute(data, f32::is_finite) } - fn compute(data: &[f32], valid: impl Fn(f32) -> bool) -> BandStats { + fn compute(data: &[f32], valid: impl Fn(f32) -> bool + Sync) -> BandStats { BandStats::compute_with(data, valid, None) } /// [`BandStats::compute`], taking the range, mean and spread from `stored` (statistics of the /// whole band, when `data` is a reduced overview of it). - fn compute_with(data: &[f32], valid: impl Fn(f32) -> bool, stored: Option) -> BandStats { - let (mut n, mut min, mut max, mut sum, mut sum2) = (0usize, f64::INFINITY, f64::NEG_INFINITY, 0.0f64, 0.0f64); - let mut whole = true; - for &v in data { - if !valid(v) { - continue; + /// + /// The sums are made on this thread, in the band's order, so the mean and spread are the + /// same to the bit. What doesn't depend on the order (the count, the range, whether every + /// value is whole, the histogram, the distinct values) is counted on several threads for a + /// large band ([`on_runs`]), the range while this thread sums. + fn compute_with(data: &[f32], valid: impl Fn(f32) -> bool + Sync, stored: Option) -> BandStats { + let sums = || { + let (mut sum, mut sum2) = (0.0f64, 0.0f64); + for &v in data.iter().filter(|&&v| valid(v)) { + let x = v as f64; + sum += x; + sum2 += x * x; } - let x = v as f64; - n += 1; - min = min.min(x); - max = max.max(x); - sum += x; - sum2 += x * x; - whole &= x.fract() == 0.0; + (sum, sum2) + }; + let (mut n, mut min, mut max, mut whole) = (0usize, f64::INFINITY, f64::NEG_INFINITY, true); + let (ranges, (sum, sum2)) = on_runs_beside(data, sums, |run| { + let (mut n, mut min, mut max, mut whole) = (0usize, f64::INFINITY, f64::NEG_INFINITY, true); + for &v in run.iter().filter(|&&v| valid(v)) { + let x = v as f64; + n += 1; + min = min.min(x); + max = max.max(x); + whole &= x.fract() == 0.0; + } + (n, min, max, whole) + }); + for (rn, rmin, rmax, rwhole) in ranges { + (n, min, max, whole) = (n + rn, min.min(rmin), max.max(rmax), whole && rwhole); } if n == 0 { return BandStats::default(); @@ -197,34 +212,24 @@ impl BandStats { mean = s.mean.unwrap_or(mean); std_dev = s.std_dev.unwrap_or(std_dev); } - let mut histogram = vec![0u64; HISTOGRAM_BINS]; let span = max - min; // Whole numbers within 16 bits are tracked in a bit set (fast for byte imagery); others in // a small set that gives up past 256 values. let small_ints = whole && min >= 0.0 && max <= 65_535.0; - let mut bits = if small_ints { vec![0u64; 1024] } else { Vec::new() }; - let mut distinct = std::collections::BTreeSet::new(); - let mut track = whole; let k = (HISTOGRAM_BINS as f64 - 1.0) / if span > 0.0 { span } else { 1.0 }; - for &v in data { - if !valid(v) { - continue; - } - let x = v as f64; - let bin = ((x - min) * k).round() as usize; - if let Some(h) = histogram.get_mut(bin.min(HISTOGRAM_BINS - 1)) { - *h += 1; - } - if small_ints { - let i = x as usize; - if let Some(w) = bits.get_mut(i / 64) { - *w |= 1 << (i % 64); - } - } else if track { - distinct.insert(x as i64); - track = distinct.len() <= 256; + let b = Binning { min, k, small_ints, whole }; + let mut runs = on_runs(data, |run| Tally::run(run, &valid, b)).into_iter(); + let mut tally = runs.next().unwrap_or_else(|| Tally::run(&[], &valid, b)); + for t in runs { + tally.histogram.iter_mut().zip(&t.histogram).for_each(|(a, b)| *a += b); + tally.bits.iter_mut().zip(&t.bits).for_each(|(a, b)| *a |= b); + tally.track &= t.track; + if tally.track { + tally.distinct.extend(t.distinct); + tally.track = tally.distinct.len() <= 256; } } + let Tally { histogram, bits, distinct, track } = tally; let categories = if small_ints { let n: u32 = bits.iter().map(|w| w.count_ones()).sum(); (n <= 256).then(|| { @@ -258,6 +263,91 @@ impl BandStats { } } +/// How [`Tally::run`] counts a band's values: its histogram's bins from `min`, `k` bins to a unit, +/// and how its distinct values are kept (`small_ints`: in a bit set; else, when `whole`, in a +/// set that gives up past 256 values). +#[derive(Clone, Copy)] +struct Binning { + min: f64, + k: f64, + small_ints: bool, + whole: bool, +} + +/// A band's values counted: its histogram, and its distinct values (the bit set of whole +/// numbers within 16 bits, or the set and whether it still keeps them). +struct Tally { + histogram: Vec, + bits: Vec, + distinct: std::collections::BTreeSet, + track: bool, +} + +/// Values below which a band is counted on one thread (#1252): threads pay off past about a +/// million. +const PARALLEL_VALUES: usize = 1 << 20; + +/// `f` over runs of `data` that cover it in order, its results in that order: on several +/// threads for a large band (at most 16, the calling thread among them; one in the browser), on +/// this one otherwise. A worker's panic is raised again here. +fn on_runs(data: &[f32], f: impl Fn(&[f32]) -> T + Sync) -> Vec { + on_runs_beside(data, || (), f).0 +} + +/// [`on_runs`], the calling thread running `beside` first, while the other threads start on +/// their runs. +fn on_runs_beside(data: &[f32], beside: impl FnOnce() -> U, f: impl Fn(&[f32]) -> T + Sync) -> (Vec, U) { + let threads = if cfg!(target_arch = "wasm32") || data.len() < PARALLEL_VALUES { + 1 + } else { + std::thread::available_parallelism().map_or(1, |n| n.get()).min(16).min(data.len() / (PARALLEL_VALUES / 4)) + }; + if threads <= 1 { + let beside = beside(); + return (vec![f(data)], beside); + } + let mut runs = data.chunks(data.len().div_ceil(threads)); + let first = runs.next().unwrap_or_default(); + let f = &f; + std::thread::scope(|scope| { + let handles: Vec<_> = runs.map(|run| scope.spawn(move || f(run))).collect(); + let beside = beside(); + let mut out = vec![f(first)]; + out.extend(handles.into_iter().map(|h| h.join().unwrap_or_else(|p| std::panic::resume_unwind(p)))); + (out, beside) + }) +} + +impl Tally { + /// `data`'s valid values counted on this thread. + fn run(data: &[f32], valid: &impl Fn(f32) -> bool, b: Binning) -> Tally { + let mut histogram = vec![0u64; HISTOGRAM_BINS]; + let mut bits = if b.small_ints { vec![0u64; 1024] } else { Vec::new() }; + let mut distinct = std::collections::BTreeSet::new(); + let mut track = b.whole; + for &v in data { + if !valid(v) { + continue; + } + let x = v as f64; + let bin = ((x - b.min) * b.k).round() as usize; + if let Some(h) = histogram.get_mut(bin.min(HISTOGRAM_BINS - 1)) { + *h += 1; + } + if b.small_ints { + let i = x as usize; + if let Some(w) = bits.get_mut(i / 64) { + *w |= 1 << (i % 64); + } + } else if track { + distinct.insert(x as i64); + track = distinct.len() <= 256; + } + } + Tally { histogram, bits, distinct, track } + } +} + /// The no-data value the readers leave out of [`Raster::rounded`], by one rule for every reader: /// /// - The value is the file's when it is finite, else the sidecar's (`.aux.xml`) when that is @@ -1148,6 +1238,67 @@ pub fn read_prj(path: &Path) -> Option { pub(crate) mod tests { use super::*; + /// A large band's statistics, counted on several threads, are those of one pass on one + /// thread (#1252) to the bit: the sums in the band's order, the counts, range, histogram and + /// distinct values whatever the threads. No-data, NaN, −0, byte-like values, more than 256 + /// whole values, and fractions. + #[test] + fn band_stats_on_several_threads_are_one_threads() { + // The statistics as one pass on one thread made them before #1252. + fn one_pass(data: &[f32], valid: impl Fn(f32) -> bool) -> BandStats { + let (mut n, mut min, mut max, mut sum, mut sum2) = (0usize, f64::INFINITY, f64::NEG_INFINITY, 0.0f64, 0.0f64); + let mut whole = true; + for &v in data.iter().filter(|&&v| valid(v)) { + let x = v as f64; + n += 1; + min = min.min(x); + max = max.max(x); + sum += x; + sum2 += x * x; + whole &= x.fract() == 0.0; + } + if n == 0 { + return BandStats::default(); + } + let mean = sum / n as f64; + let std_dev = (sum2 / n as f64 - mean * mean).max(0.0).sqrt(); + let span = max - min; + let k = (HISTOGRAM_BINS as f64 - 1.0) / if span > 0.0 { span } else { 1.0 }; + let mut histogram = vec![0u64; HISTOGRAM_BINS]; + let mut distinct = std::collections::BTreeSet::new(); + for &v in data.iter().filter(|&&v| valid(v)) { + let x = v as f64; + histogram[(((x - min) * k).round() as usize).min(HISTOGRAM_BINS - 1)] += 1; + if whole { + distinct.insert(x as i64); + } + } + let categories = (whole && distinct.len() <= 256).then(|| distinct.into_iter().map(|v| v as f64).collect()); + BandStats { count: n, min, max, mean, std_dev, histogram, categories } + } + let n = 2_000_003usize; + let at = |i: usize| (i * 7919 + i / 1013 * 31) % 100_003; + let bands: [(&str, Vec); 5] = [ + ("bytes", (0..n).map(|i| (at(i) % 256) as f32).collect()), + ("few whole", (0..n).map(|i| (at(i) % 200) as f32 * 3.0 - 250.0).collect()), + ("wide whole", (0..n).map(|i| at(i) as f32 - 7.0).collect()), + ("fractions", (0..n).map(|i| at(i) as f32 * 0.013 - 400.0).collect()), + ("gaps", (0..n).map(|i| [f32::NAN, -0.0, 0.0, -9999.0, 17.5][at(i) % 5]).collect()), + ]; + let bits = |s: &BandStats| (s.count, s.min.to_bits(), s.max.to_bits(), s.mean.to_bits(), s.std_dev.to_bits()); + for (name, data) in &bands { + let valid = |v: f32| v.is_finite() && v != -9999.0; + let (got, want) = (BandStats::compute(data, valid), one_pass(data, valid)); + assert_eq!(bits(&got), bits(&want), "{name}"); + assert_eq!(got.histogram, want.histogram, "{name}"); + assert_eq!(got.categories, want.categories, "{name}"); + } + // A band too small for threads, and one with no valid value. + let small: Vec = (0..1000).map(|i| (i % 37) as f32).collect(); + assert_eq!(BandStats::from_values(&small), one_pass(&small, f32::is_finite)); + assert_eq!(BandStats::from_values(&vec![f32::NAN; n]), BandStats::default()); + } + /// A band the machine can't hold is an error naming its size, not the end of the program /// (#1291), on 64-bit (the reservation fails) and wasm32 (the size overflows) alike; one it can /// is the size asked for. From b79dc30a37c33c0fc030c0f5cda305ee8e8ebeef Mon Sep 17 00:00:00 2001 From: Dominik Strebel Date: Fri, 6 Nov 2026 15:37:27 +0100 Subject: [PATCH 2/2] fix(data): one pass on one thread, threads that can't start counted here, the first equal value kept, the merge tested on any machine --- crates/data/src/raster/mod.rs | 149 ++++++++++++++++++++++------------ 1 file changed, 95 insertions(+), 54 deletions(-) diff --git a/crates/data/src/raster/mod.rs b/crates/data/src/raster/mod.rs index 07a516728..e352f9698 100644 --- a/crates/data/src/raster/mod.rs +++ b/crates/data/src/raster/mod.rs @@ -172,36 +172,41 @@ impl BandStats { /// [`BandStats::compute`], taking the range, mean and spread from `stored` (statistics of the /// whole band, when `data` is a reduced overview of it). - /// - /// The sums are made on this thread, in the band's order, so the mean and spread are the - /// same to the bit. What doesn't depend on the order (the count, the range, whether every - /// value is whole, the histogram, the distinct values) is counted on several threads for a - /// large band ([`on_runs`]), the range while this thread sums. fn compute_with(data: &[f32], valid: impl Fn(f32) -> bool + Sync, stored: Option) -> BandStats { - let sums = || { - let (mut sum, mut sum2) = (0.0f64, 0.0f64); - for &v in data.iter().filter(|&&v| valid(v)) { - let x = v as f64; - sum += x; - sum2 += x * x; - } - (sum, sum2) - }; - let (mut n, mut min, mut max, mut whole) = (0usize, f64::INFINITY, f64::NEG_INFINITY, true); - let (ranges, (sum, sum2)) = on_runs_beside(data, sums, |run| { - let (mut n, mut min, mut max, mut whole) = (0usize, f64::INFINITY, f64::NEG_INFINITY, true); + BandStats::compute_on(data, valid, stored, threads_for(data.len())) + } + + /// [`BandStats::compute_with`] on `threads` threads. On one, one pass makes the count, range, + /// sums and whether every value is whole, and a second the histogram and distinct values. On + /// several, the sums are made on this thread, in the band's order, so the mean and spread are + /// the same to the bit; what doesn't depend on the order is counted on the others + /// ([`on_runs`]), the range while this thread sums. The range keeps the first of equal values + /// (−0 or 0) in the band's order either way. + fn compute_on(data: &[f32], valid: impl Fn(f32) -> bool + Sync, stored: Option, threads: usize) -> BandStats { + let ranged = |run: &[f32], sum: bool| { + let (mut n, mut min, mut max, mut whole, mut sums) = (0usize, f64::INFINITY, f64::NEG_INFINITY, true, (0.0f64, 0.0f64)); for &v in run.iter().filter(|&&v| valid(v)) { let x = v as f64; n += 1; - min = min.min(x); - max = max.max(x); + (min, max) = (lower(min, x), higher(max, x)); whole &= x.fract() == 0.0; + if sum { + sums = (sums.0 + x, sums.1 + x * x); + } } - (n, min, max, whole) - }); - for (rn, rmin, rmax, rwhole) in ranges { - (n, min, max, whole) = (n + rn, min.min(rmin), max.max(rmax), whole && rwhole); - } + (n, min, max, whole, sums) + }; + let (n, mut min, mut max, whole, (sum, sum2)) = if threads <= 1 { + ranged(data, true) + } else { + let sums = || ranged(data, true).4; + let (runs, sums) = on_runs_beside(data, threads, sums, |run| ranged(run, false)); + let mut all = (0usize, f64::INFINITY, f64::NEG_INFINITY, true, sums); + for (rn, rmin, rmax, rwhole, _) in runs { + all = (all.0 + rn, lower(all.1, rmin), higher(all.2, rmax), all.3 && rwhole, all.4); + } + all + }; if n == 0 { return BandStats::default(); } @@ -218,7 +223,7 @@ impl BandStats { let small_ints = whole && min >= 0.0 && max <= 65_535.0; let k = (HISTOGRAM_BINS as f64 - 1.0) / if span > 0.0 { span } else { 1.0 }; let b = Binning { min, k, small_ints, whole }; - let mut runs = on_runs(data, |run| Tally::run(run, &valid, b)).into_iter(); + let mut runs = on_runs(data, threads, |run| Tally::run(run, &valid, b)).into_iter(); let mut tally = runs.next().unwrap_or_else(|| Tally::run(&[], &valid, b)); for t in runs { tally.histogram.iter_mut().zip(&t.histogram).for_each(|(a, b)| *a += b); @@ -287,22 +292,37 @@ struct Tally { /// million. const PARALLEL_VALUES: usize = 1 << 20; -/// `f` over runs of `data` that cover it in order, its results in that order: on several -/// threads for a large band (at most 16, the calling thread among them; one in the browser), on -/// this one otherwise. A worker's panic is raised again here. -fn on_runs(data: &[f32], f: impl Fn(&[f32]) -> T + Sync) -> Vec { - on_runs_beside(data, || (), f).0 +/// `x` when it is below `min`, else `min`: equal values (−0 and 0) keep the first. +fn lower(min: f64, x: f64) -> f64 { + if x < min { x } else { min } } -/// [`on_runs`], the calling thread running `beside` first, while the other threads start on -/// their runs. -fn on_runs_beside(data: &[f32], beside: impl FnOnce() -> U, f: impl Fn(&[f32]) -> T + Sync) -> (Vec, U) { - let threads = if cfg!(target_arch = "wasm32") || data.len() < PARALLEL_VALUES { +/// `x` when it is above `max`, else `max`: equal values (−0 and 0) keep the first. +fn higher(max: f64, x: f64) -> f64 { + if x > max { x } else { max } +} + +/// The threads to count a band of `len` values on: at most 16, one per 256 Ki values, and one +/// below [`PARALLEL_VALUES`] or in the browser. +fn threads_for(len: usize) -> usize { + if cfg!(target_arch = "wasm32") || len < PARALLEL_VALUES { 1 } else { - std::thread::available_parallelism().map_or(1, |n| n.get()).min(16).min(data.len() / (PARALLEL_VALUES / 4)) - }; - if threads <= 1 { + std::thread::available_parallelism().map_or(1, |n| n.get()).min(16).min(len / (PARALLEL_VALUES / 4)) + } +} + +/// `f` over `threads` runs of `data` that cover it in order, its results in that order, the +/// calling thread taking the first. A worker's panic is raised again here; a run whose thread +/// can't start is made on this one. +fn on_runs(data: &[f32], threads: usize, f: impl Fn(&[f32]) -> T + Sync) -> Vec { + on_runs_beside(data, threads, || (), f).0 +} + +/// [`on_runs`], the calling thread running `beside` first, while the other threads start on +/// their runs. +fn on_runs_beside(data: &[f32], threads: usize, beside: impl FnOnce() -> U, f: impl Fn(&[f32]) -> T + Sync) -> (Vec, U) { + if threads <= 1 || data.is_empty() { let beside = beside(); return (vec![f(data)], beside); } @@ -310,10 +330,16 @@ fn on_runs_beside(data: &[f32], beside: impl FnOnce() -> U, f: impl let first = runs.next().unwrap_or_default(); let f = &f; std::thread::scope(|scope| { - let handles: Vec<_> = runs.map(|run| scope.spawn(move || f(run))).collect(); + let started: Vec<_> = runs.map(|run| std::thread::Builder::new().spawn_scoped(scope, move || f(run)).map_err(|_| run)).collect(); let beside = beside(); let mut out = vec![f(first)]; - out.extend(handles.into_iter().map(|h| h.join().unwrap_or_else(|p| std::panic::resume_unwind(p)))); + for run in started { + out.push(match run { + Ok(h) => h.join().unwrap_or_else(|p| std::panic::resume_unwind(p)), + // The OS refused a thread: made here. + Err(run) => f(run), + }); + } (out, beside) }) } @@ -1238,21 +1264,28 @@ pub fn read_prj(path: &Path) -> Option { pub(crate) mod tests { use super::*; - /// A large band's statistics, counted on several threads, are those of one pass on one - /// thread (#1252) to the bit: the sums in the band's order, the counts, range, histogram and - /// distinct values whatever the threads. No-data, NaN, −0, byte-like values, more than 256 - /// whole values, and fractions. + /// A large band's statistics, counted on any number of threads, are those of one pass on + /// one thread (#1252) to the bit: the sums in the band's order, the counts, range, histogram + /// and distinct values whatever the runs. No-data, NaN, −0 and 0 (the first kept), byte-like + /// values, 257 small whole values, more than 256 whole values, runs that each have at most + /// 256 but more together, and fractions. The thread count is forced, so the merge is tested on + /// a machine with one core too. #[test] fn band_stats_on_several_threads_are_one_threads() { - // The statistics as one pass on one thread made them before #1252. + // The statistics as one pass on one thread made them before #1252 (the first of equal + // values kept in the range). fn one_pass(data: &[f32], valid: impl Fn(f32) -> bool) -> BandStats { let (mut n, mut min, mut max, mut sum, mut sum2) = (0usize, f64::INFINITY, f64::NEG_INFINITY, 0.0f64, 0.0f64); let mut whole = true; for &v in data.iter().filter(|&&v| valid(v)) { let x = v as f64; n += 1; - min = min.min(x); - max = max.max(x); + if x < min { + min = x; + } + if x > max { + max = x; + } sum += x; sum2 += x * x; whole &= x.fract() == 0.0; @@ -1277,26 +1310,34 @@ pub(crate) mod tests { BandStats { count: n, min, max, mean, std_dev, histogram, categories } } let n = 2_000_003usize; - let at = |i: usize| (i * 7919 + i / 1013 * 31) % 100_003; - let bands: [(&str, Vec); 5] = [ + let at = |i: usize| ((i as u64 * 7919 + i as u64 / 1013 * 31) % 100_003) as usize; + let bands: [(&str, Vec); 7] = [ ("bytes", (0..n).map(|i| (at(i) % 256) as f32).collect()), + ("257 small whole", (0..n).map(|i| (at(i) % 257) as f32).collect()), ("few whole", (0..n).map(|i| (at(i) % 200) as f32 * 3.0 - 250.0).collect()), + // Each half has 200 values, together 400: two runs keep theirs, the merge gives up. + ("halves", (0..n).map(|i| -((at(i) % 200) as f32) - if i < n / 2 { 1.0 } else { 201.0 }).collect()), ("wide whole", (0..n).map(|i| at(i) as f32 - 7.0).collect()), ("fractions", (0..n).map(|i| at(i) as f32 * 0.013 - 400.0).collect()), ("gaps", (0..n).map(|i| [f32::NAN, -0.0, 0.0, -9999.0, 17.5][at(i) % 5]).collect()), ]; let bits = |s: &BandStats| (s.count, s.min.to_bits(), s.max.to_bits(), s.mean.to_bits(), s.std_dev.to_bits()); + let valid = |v: f32| v.is_finite() && v != -9999.0; for (name, data) in &bands { - let valid = |v: f32| v.is_finite() && v != -9999.0; - let (got, want) = (BandStats::compute(data, valid), one_pass(data, valid)); - assert_eq!(bits(&got), bits(&want), "{name}"); - assert_eq!(got.histogram, want.histogram, "{name}"); - assert_eq!(got.categories, want.categories, "{name}"); + let want = one_pass(data, valid); + for threads in [1, 2, 7, 16] { + let got = BandStats::compute_on(data, valid, None, threads); + assert_eq!(bits(&got), bits(&want), "{name} on {threads}"); + assert_eq!(got.histogram, want.histogram, "{name} on {threads}"); + assert_eq!(got.categories, want.categories, "{name} on {threads}"); + } } + assert!(BandStats::compute_on(&bands[1].1, valid, None, 2).categories.is_none(), "257 small whole values"); + assert!(BandStats::compute_on(&bands[3].1, valid, None, 2).categories.is_none(), "400 in two runs of 200"); // A band too small for threads, and one with no valid value. let small: Vec = (0..1000).map(|i| (i % 37) as f32).collect(); assert_eq!(BandStats::from_values(&small), one_pass(&small, f32::is_finite)); - assert_eq!(BandStats::from_values(&vec![f32::NAN; n]), BandStats::default()); + assert_eq!(BandStats::compute_on(&vec![f32::NAN; n], valid, None, 4), BandStats::default()); } /// A band the machine can't hold is an error naming its size, not the end of the program