From 8a3a5d8fc75af3e5e1e43013c80edba5761662bb Mon Sep 17 00:00:00 2001 From: Dominik Strebel Date: Sun, 1 Nov 2026 15:49:41 +0100 Subject: [PATCH 1/3] perf(analysis): IDW and kriging rows on several cores IDW and Ordinary Kriging make each raster row on its own, so their rows now go through parallel::map_rows (#1252): rows taken one at a time by up to 16 threads, the calling thread among them, above MIN_ROW_WORK cell-point pairs, and on one thread in the browser. Each thread has its own distance buffer, and for kriging its own KrigingBuffers. The surfaces are the same, bit for bit, as on one thread. The helper now runs over indices (try_each_on), which try_map_on wraps, and a drop guard sets the stop flag when a worker panics, so the other threads stop after the items they are making. Select by Location's notice for selecting features without geometry has a test. On 40 threads, IDW of 10,000 points onto 230 x 100 cells takes 0.55 s instead of 6.1 s, and of 1,000 points onto 920 x 400 cells 0.88 s instead of 9.9 s. --- crates/analysis/src/kriging.rs | 53 +++++++++++--- crates/analysis/src/parallel.rs | 119 ++++++++++++++++++++++++++++--- crates/analysis/src/rasterops.rs | 103 ++++++++++++++++++-------- crates/analysis/src/tools.rs | 18 +++++ 4 files changed, 245 insertions(+), 48 deletions(-) diff --git a/crates/analysis/src/kriging.rs b/crates/analysis/src/kriging.rs index cec80fab5..3f8ba75ad 100644 --- a/crates/analysis/src/kriging.rs +++ b/crates/analysis/src/kriging.rs @@ -313,6 +313,19 @@ pub fn kriging( params: &VariogramParams, neighbours: usize, name: &str, +) -> Result { + kriging_on(None, points, extent, cell_deg, params, neighbours, name) +} + +/// [`kriging`] on `threads` threads, or as many as its work calls for. +fn kriging_on( + threads: Option, + points: &[(Coord, f64)], + extent: &Extent, + cell_deg: f64, + params: &VariogramParams, + neighbours: usize, + name: &str, ) -> Result { if points.len() < 2 { return Err("kriging needs at least two points with values".into()); @@ -334,12 +347,13 @@ pub fn kriging( return Err(format!("{w} × {h} cells from {} points is too much work; use a larger cell size", points.len())); } - let mut pred_data = Vec::with_capacity(w * h); - let mut var_data = Vec::with_capacity(w * h); - let mut near: Vec<(f64, f64, Coord)> = Vec::with_capacity(points.len()); - let mut buffers = KrigingBuffers::new(); - - for row in 0..h { + // Each row on its own, on several threads when the work is large (#1252); a thread has its + // own distances and system buffers. + let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, w.saturating_mul(h).saturating_mul(work_per_cell))); + let scratch = || (Vec::<(f64, f64, Coord)>::with_capacity(points.len()), KrigingBuffers::new()); + let rows = crate::parallel::map_rows(h, threads, scratch, |(near, buffers), row| { + let mut pred_data = Vec::with_capacity(w); + let mut var_data = Vec::with_capacity(w); let kx = (extent.max_y - (row as f64 + 0.5) * cell_deg).to_radians().cos().abs().max(1e-6); for col in 0..w { let c = Coord::new(extent.min_x + (col as f64 + 0.5) * cell_deg, extent.max_y - (row as f64 + 0.5) * cell_deg); @@ -356,11 +370,14 @@ pub fn kriging( near.truncate(k); } - let (val, var) = ordinary_kriging_point(c, &near, params, kx, &mut buffers); + let (val, var) = ordinary_kriging_point(c, near, params, kx, buffers); pred_data.push(if val.is_finite() { val as f32 } else { NODATA }); var_data.push(if var.is_finite() { var as f32 } else { NODATA }); } - } + (pred_data, var_data) + }); + let (pred, var): (Vec>, Vec>) = rows.into_iter().unzip(); + let (pred_data, var_data) = (pred.concat(), var.concat()); let gt = GeoTransform { origin_x: extent.min_x, origin_y: extent.max_y, pixel_width: cell_deg, pixel_height: -cell_deg }; Raster::new(name, w as u32, h as u32, vec![pred_data, var_data], gt, Crs::LonLat, Some(NODATA as f64), ColorHint::Gray, "Float32") @@ -436,6 +453,26 @@ mod tests { } } + /// Rows made on several threads give the same surfaces, bit for bit, as one thread (#1252). + #[test] + fn kriging_on_several_threads_is_kriging_on_one() { + let pts: Vec<(Coord, f64)> = (0..300) + .map(|i| { + let (x, y) = ((i * 37 % 300) as f64 / 30.0, (i * 91 % 300) as f64 / 30.0); + (Coord::new(x, 40.0 + y), (x * 1.3).sin() * 100.0 + y * 7.0) + }) + .collect(); + let extent = Extent::new(0.0, 40.0, 10.0, 50.0); + let params = auto_variogram(&pts, &extent, VariogramModel::Spherical); + let one = kriging_on(Some(1), &pts, &extent, 0.2, ¶ms, 12, "k").unwrap(); + for threads in [2, 7, 16] { + let many = kriging_on(Some(threads), &pts, &extent, 0.2, ¶ms, 12, "k").unwrap(); + for (a, b) in one.bands.iter().zip(&many.bands) { + assert!(a.data.iter().map(|v| v.to_bits()).eq(b.data.iter().map(|v| v.to_bits())), "{threads} threads"); + } + } + } + #[test] fn test_high_latitude_metric_consistency() { let lat = 60.0; diff --git a/crates/analysis/src/parallel.rs b/crates/analysis/src/parallel.rs index 16215eefe..42e3852c7 100644 --- a/crates/analysis/src/parallel.rs +++ b/crates/analysis/src/parallel.rs @@ -12,13 +12,46 @@ pub(crate) const BATCH: usize = 16; /// Each worker's stack, as large as the main thread's usually is. const STACK_BYTES: usize = 8 << 20; +/// The work below which a raster tool's rows are made on one thread (cell-point pairs, say: a +/// few milliseconds of it). +pub(crate) const MIN_ROW_WORK: usize = 2_000_000; + /// How many threads to make `items` items on: one below `min_items`, always one in the browser, /// and no more than there are batches. pub(crate) fn threads_for(items: usize, min_items: usize) -> usize { if cfg!(target_arch = "wasm32") || items < min_items { return 1; } - std::thread::available_parallelism().map_or(1, |n| n.get()).min(MAX_THREADS).min(items.div_ceil(BATCH)).max(1) + cores().min(items.div_ceil(BATCH)).max(1) +} + +/// How many threads to make `rows` raster rows on, `work` in all: one below [`MIN_ROW_WORK`] and +/// in the browser, and no more than there are rows. +pub(crate) fn threads_for_rows(rows: usize, work: usize) -> usize { + if cfg!(target_arch = "wasm32") || work < MIN_ROW_WORK { + return 1; + } + cores().min(rows).max(1) +} + +/// The cores to use, at most [`MAX_THREADS`]. +fn cores() -> usize { + std::thread::available_parallelism().map_or(1, |n| n.get()).min(MAX_THREADS) +} + +/// `f` on each of `rows` raster rows by its number, the rows in their order, on `threads` +/// threads ([`threads_for_rows`]). A row is costly enough to be taken one at a time, so rows +/// that cost more than others are shared out too. Each thread makes its own `scratch()` once. +pub(crate) fn map_rows( + rows: usize, + threads: usize, + scratch: impl Fn() -> S + Sync, + f: impl Fn(&mut S, usize) -> R + Sync, +) -> Vec { + match try_each_on(rows, 1, threads, scratch, |s, i| Ok::(f(s, i))) { + Ok(made) => made, + Err(never) => match never {}, + } } /// `f` on each of `items` with its index, the results in their order: on several threads when @@ -45,26 +78,53 @@ pub(crate) fn try_map_on( threads: usize, scratch: impl Fn() -> S + Sync, f: impl Fn(&mut S, usize, &T) -> Result + Sync, +) -> Result, E> { + try_each_on(items.len(), BATCH, threads, scratch, |s, i| f(s, i, &items[i])) +} + +/// Sets the stop flag when the thread it lives on panics, so the others stop after the items +/// they are making rather than make the rest for nothing. +struct StopOnPanic<'a>(&'a std::sync::atomic::AtomicBool); + +impl Drop for StopOnPanic<'_> { + fn drop(&mut self) { + if std::thread::panicking() { + self.0.store(true, std::sync::atomic::Ordering::Relaxed); + } + } +} + +/// [`try_map_on`] over the indices `0..len`, taken `batch` at a time. +fn try_each_on( + len: usize, + batch: usize, + threads: usize, + scratch: impl Fn() -> S + Sync, + f: impl Fn(&mut S, usize) -> Result + Sync, ) -> Result, E> { use std::sync::atomic::{AtomicBool, AtomicUsize, Ordering}; - if threads <= 1 || items.len() <= 1 { + if threads <= 1 || len <= 1 { let mut s = scratch(); - return items.iter().enumerate().map(|(i, item)| f(&mut s, i, item)).collect(); + return (0..len).map(|i| f(&mut s, i)).collect(); } + let batch = batch.max(1); let next = AtomicUsize::new(0); let failed = AtomicBool::new(false); let work = || { + let _stop = StopOnPanic(&failed); let mut s = scratch(); let mut made = Vec::new(); while !failed.load(Ordering::Relaxed) { - let from = next.fetch_add(BATCH, Ordering::Relaxed); - let Some(batch) = items.get(from..items.len().min(from + BATCH)).filter(|b| !b.is_empty()) else { break }; - for (k, item) in batch.iter().enumerate() { - match f(&mut s, from + k, item) { - Ok(r) => made.push((from + k, r)), + let from = next.fetch_add(batch, Ordering::Relaxed); + if from >= len { + break; + } + for i in from..len.min(from.saturating_add(batch)) { + match f(&mut s, i) { + Ok(r) => made.push((i, r)), Err(e) => { failed.store(true, Ordering::Relaxed); - return Err((from + k, e)); + return Err((i, e)); } } } @@ -82,7 +142,7 @@ pub(crate) fn try_map_on( parts.push(h.join().unwrap_or_else(|e| std::panic::resume_unwind(e))); } }); - let mut made: Vec<(usize, R)> = Vec::with_capacity(items.len()); + let mut made: Vec<(usize, R)> = Vec::with_capacity(len); let mut first: Option<(usize, E)> = None; for part in parts { match part { @@ -217,6 +277,45 @@ mod tests { assert!(ran.load(std::sync::atomic::Ordering::Relaxed) < items.len()); } + #[test] + fn rows_come_back_in_their_order_one_at_a_time() { + let alone: Vec = (0..37).map(|r| r * 3).collect(); + for threads in [1, 2, 5, 16, 64] { + assert_eq!(map_rows(37, threads, || (), |_, r| r * 3), alone, "{threads} threads"); + } + assert!(map_rows(0, 4, || (), |_, r| r).is_empty()); + // Small work, and the browser: one thread; never more threads than rows. + assert_eq!(threads_for_rows(100, MIN_ROW_WORK - 1), 1); + assert!(threads_for_rows(3, usize::MAX) <= 3); + } + + #[test] + fn a_workers_panic_stops_the_other_threads() { + // A worker panics on its first item; the others stop after the items they are making, + // far short of all of them. Unless no worker could start: then nothing panics. + let items: Vec = (0..2_000).collect(); + let caller = std::thread::current().id(); + let ran = std::sync::atomic::AtomicUsize::new(0); + let r = std::panic::catch_unwind(|| { + map_on( + &items, + 4, + || (), + |_, _, x| { + ran.fetch_add(1, std::sync::atomic::Ordering::Relaxed); + if std::thread::current().id() != caller { + panic!("a worker's panic"); + } + std::thread::sleep(std::time::Duration::from_millis(1)); + *x + }, + ) + }); + if r.is_err() { + assert!(ran.load(std::sync::atomic::Ordering::Relaxed) < items.len() / 2, "{}", ran.load(std::sync::atomic::Ordering::Relaxed)); + } + } + #[test] fn no_more_threads_than_batches() { assert_eq!(threads_for(10, 256), 1); diff --git a/crates/analysis/src/rasterops.rs b/crates/analysis/src/rasterops.rs index e5879536c..fad4f502c 100644 --- a/crates/analysis/src/rasterops.rs +++ b/crates/analysis/src/rasterops.rs @@ -374,6 +374,19 @@ pub fn viewshed(r: &Raster, b: usize, observer: Coord, height: f64, radius_m: f6 /// Inverse distance weighted interpolation of point values onto a longitude/latitude grid, from /// each cell's `neighbours` nearest points, or from every point (as QGIS does) when it is 0. pub fn idw(points: &[(Coord, f64)], extent: &Extent, cell_deg: f64, power: f64, neighbours: usize, name: &str) -> Result { + idw_on(None, points, extent, cell_deg, power, neighbours, name) +} + +/// [`idw`] on `threads` threads, or as many as its work calls for. +fn idw_on( + threads: Option, + points: &[(Coord, f64)], + extent: &Extent, + cell_deg: f64, + power: f64, + neighbours: usize, + name: &str, +) -> Result { if points.len() < 2 { return Err("interpolation needs at least two points with values".into()); } @@ -395,37 +408,46 @@ pub fn idw(points: &[(Coord, f64)], extent: &Extent, cell_deg: f64, power: f64, return Err("the distance coefficient must be zero or more".into()); } let k = if neighbours == 0 { points.len() } else { neighbours.clamp(1, points.len()) }; - let mut data = Vec::with_capacity(w * h); - let mut near: Vec<(f64, f64)> = Vec::with_capacity(points.len()); - for row in 0..h { - // East-west degrees shrink with latitude, which only changes by row. - let kx = (extent.max_y - (row as f64 + 0.5) * cell_deg).to_radians().cos(); - for col in 0..w { - let c = Coord::new(extent.min_x + (col as f64 + 0.5) * cell_deg, extent.max_y - (row as f64 + 0.5) * cell_deg); - near.clear(); - near.extend(points.iter().map(|(p, v)| { - let dx = (p.x - c.x) * kx; - ((dx * dx + (p.y - c.y).powi(2)).sqrt(), *v) - })); - if k < near.len() { - near.select_nth_unstable_by(k - 1, |a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); - near.truncate(k); + // Each row on its own, on several threads when the work is large (#1252); a thread's + // distances go into its own buffer. + let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, w.saturating_mul(h).saturating_mul(points.len()))); + let rows = crate::parallel::map_rows( + h, + threads, + || Vec::with_capacity(points.len()), + |near: &mut Vec<(f64, f64)>, row| { + let mut data = Vec::with_capacity(w); + // East-west degrees shrink with latitude, which only changes by row. + let kx = (extent.max_y - (row as f64 + 0.5) * cell_deg).to_radians().cos(); + for col in 0..w { + let c = Coord::new(extent.min_x + (col as f64 + 0.5) * cell_deg, extent.max_y - (row as f64 + 0.5) * cell_deg); + near.clear(); + near.extend(points.iter().map(|(p, v)| { + let dx = (p.x - c.x) * kx; + ((dx * dx + (p.y - c.y).powi(2)).sqrt(), *v) + })); + if k < near.len() { + near.select_nth_unstable_by(k - 1, |a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); + near.truncate(k); + } + let v = if let Some(exact) = near.iter().find(|(d, _)| *d < 1e-12) { + exact.1 + } else { + // Weights relative to the nearest point's: the same ratios as 1/dᵖ, without + // overflowing for a high power. + let nearest = near.iter().map(|(d, _)| *d).fold(f64::INFINITY, f64::min); + let (num, den) = near.iter().fold((0.0, 0.0), |(n, d), (dist, v)| { + let wgt = (nearest / dist).powf(power); + (n + wgt * v, d + wgt) + }); + num / den + }; + data.push(if v.is_finite() { v as f32 } else { NODATA }); } - let v = if let Some(exact) = near.iter().find(|(d, _)| *d < 1e-12) { - exact.1 - } else { - // Weights relative to the nearest point's: the same ratios as 1/dᵖ, without - // overflowing for a high power. - let nearest = near.iter().map(|(d, _)| *d).fold(f64::INFINITY, f64::min); - let (num, den) = near.iter().fold((0.0, 0.0), |(n, d), (dist, v)| { - let wgt = (nearest / dist).powf(power); - (n + wgt * v, d + wgt) - }); - num / den - }; - data.push(if v.is_finite() { v as f32 } else { NODATA }); - } - } + data + }, + ); + let data = rows.concat(); let gt = GeoTransform { origin_x: extent.min_x, origin_y: extent.max_y, pixel_width: cell_deg, pixel_height: -cell_deg }; Raster::new(name, w as u32, h as u32, vec![data], gt, Crs::LonLat, Some(NODATA as f64), ColorHint::Gray, "Float32") .map_err(|e| e.to_string()) @@ -1084,6 +1106,27 @@ impl Expr { mod tests { use super::*; + /// Rows made on several threads give the same surface, bit for bit, as one thread (#1252), + /// from every point and from the nearest few. + #[test] + fn idw_on_several_threads_is_idw_on_one() { + let pts: Vec<(Coord, f64)> = (0..500) + .map(|i| { + let (x, y) = ((i * 37 % 500) as f64 / 50.0, (i * 91 % 500) as f64 / 50.0); + (Coord::new(x, 40.0 + y), (x * 1.3).sin() * 100.0 + y * 7.0) + }) + .collect(); + let extent = Extent::new(0.0, 40.0, 10.0, 50.0); + for neighbours in [0, 8] { + let one = idw_on(Some(1), &pts, &extent, 0.1, 2.0, neighbours, "idw").unwrap_or_else(|e| panic!("{e}")); + for threads in [2, 7, 16] { + let many = idw_on(Some(threads), &pts, &extent, 0.1, 2.0, neighbours, "idw").unwrap_or_else(|e| panic!("{e}")); + let bits = |r: &Raster| r.bands[0].data.iter().map(|v| v.to_bits()).collect::>(); + assert_eq!(bits(&one), bits(&many), "{threads} threads, {neighbours} neighbours"); + } + } + } + #[test] fn terrain_rows_match_scalar_edges_nodata_and_orientations() { for (w, h) in [(1, 1), (1, 5), (2, 3), (3, 2), (3, 3), (9, 7)] { diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index 875781d3e..08cd2a77f 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -7295,6 +7295,24 @@ mod tests { } } + #[test] + fn select_by_location_says_when_selecting_features_lack_geometry() { + let mut polys = Dataset::new("polys", GeometryType::Polygon, vec![]); + polys.push(Some(sq(0.0, 0.0, 2.0)), vec![]); + polys.push(None, vec![]); // a selector without geometry + let mut pts = Dataset::new("pts", GeometryType::Point, vec![]); + pts.push(Some(Geometry::Point(Coord::new(1.0, 1.0))), vec![]); + pts.push(Some(Geometry::Point(Coord::new(10.0, 10.0))), vec![]); + let m = Map(HashMap::from([("polys".to_string(), polys), ("pts".to_string(), pts)])); + let res = crate::run("selectByLocation", &m, &serde_json::json!({"input": "pts", "select": "polys"})).unwrap(); + // One notice for it, and the selector with geometry still selects. + assert_eq!(res.messages, vec!["Selected 1 features in pts", "1 of 2 selecting features lack geometry, so they select nothing"]); + match res.output { + Output::Selection { rows, .. } => assert_eq!(rows, vec![0]), + _ => panic!("expected selection"), + } + } + #[test] fn select_by_location_invert_excludes_features_without_geometry() { let mut polys = Dataset::new("polys", GeometryType::Polygon, vec![]); From d11e2fcb93123f657fe284e04a0349e7864b5211 Mon Sep 17 00:00:00 2001 From: Dominik Strebel Date: Sun, 1 Nov 2026 16:09:03 +0100 Subject: [PATCH 2/3] perf(analysis): bound row threads by their scratch, and lower their threshold Threads for a raster tool's rows are capped so that their scratch (a buffer per point for IDW and kriging) stays within 256 MiB beyond one thread's, so millions of points don't reserve a copy on each of 16 threads. MIN_ROW_WORK drops to 200,000 cell-point pairs, about 5 ms; the old comment understated it tenfold. The rows go into the surface one at a time instead of through concat, so it isn't held twice. The panic test asserts the panic unconditionally. --- crates/analysis/src/kriging.rs | 12 +++++++++--- crates/analysis/src/parallel.rs | 33 +++++++++++++++++++------------- crates/analysis/src/rasterops.rs | 8 ++++++-- 3 files changed, 35 insertions(+), 18 deletions(-) diff --git a/crates/analysis/src/kriging.rs b/crates/analysis/src/kriging.rs index 3f8ba75ad..c01cd8c7c 100644 --- a/crates/analysis/src/kriging.rs +++ b/crates/analysis/src/kriging.rs @@ -349,7 +349,9 @@ fn kriging_on( // Each row on its own, on several threads when the work is large (#1252); a thread has its // own distances and system buffers. - let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, w.saturating_mul(h).saturating_mul(work_per_cell))); + let scratch_bytes = points.len().saturating_mul(std::mem::size_of::<(f64, f64, Coord)>()); + let threads = + threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, w.saturating_mul(h).saturating_mul(work_per_cell), scratch_bytes)); let scratch = || (Vec::<(f64, f64, Coord)>::with_capacity(points.len()), KrigingBuffers::new()); let rows = crate::parallel::map_rows(h, threads, scratch, |(near, buffers), row| { let mut pred_data = Vec::with_capacity(w); @@ -376,8 +378,12 @@ fn kriging_on( } (pred_data, var_data) }); - let (pred, var): (Vec>, Vec>) = rows.into_iter().unzip(); - let (pred_data, var_data) = (pred.concat(), var.concat()); + // Each row dropped as it is copied, so the surfaces aren't held twice. + let (mut pred_data, mut var_data) = (Vec::with_capacity(w * h), Vec::with_capacity(w * h)); + for (p, v) in rows { + pred_data.extend(p); + var_data.extend(v); + } let gt = GeoTransform { origin_x: extent.min_x, origin_y: extent.max_y, pixel_width: cell_deg, pixel_height: -cell_deg }; Raster::new(name, w as u32, h as u32, vec![pred_data, var_data], gt, Crs::LonLat, Some(NODATA as f64), ColorHint::Gray, "Float32") diff --git a/crates/analysis/src/parallel.rs b/crates/analysis/src/parallel.rs index 42e3852c7..728fe64f4 100644 --- a/crates/analysis/src/parallel.rs +++ b/crates/analysis/src/parallel.rs @@ -12,9 +12,11 @@ pub(crate) const BATCH: usize = 16; /// Each worker's stack, as large as the main thread's usually is. const STACK_BYTES: usize = 8 << 20; -/// The work below which a raster tool's rows are made on one thread (cell-point pairs, say: a -/// few milliseconds of it). -pub(crate) const MIN_ROW_WORK: usize = 2_000_000; +/// The work below which a raster tool's rows are made on one thread: cell-point pairs, say, +/// about 30 ns each for IDW, so some 5 ms of it, against tens of microseconds to start a thread. +pub(crate) const MIN_ROW_WORK: usize = 200_000; +/// The scratch memory a raster tool's threads may take together, beyond one thread's. +const SCRATCH_BUDGET: usize = 256 << 20; /// How many threads to make `items` items on: one below `min_items`, always one in the browser, /// and no more than there are batches. @@ -25,13 +27,16 @@ pub(crate) fn threads_for(items: usize, min_items: usize) -> usize { cores().min(items.div_ceil(BATCH)).max(1) } -/// How many threads to make `rows` raster rows on, `work` in all: one below [`MIN_ROW_WORK`] and -/// in the browser, and no more than there are rows. -pub(crate) fn threads_for_rows(rows: usize, work: usize) -> usize { +/// How many threads to make `rows` raster rows on, `work` in all, each thread taking +/// `scratch_bytes` of scratch: one below [`MIN_ROW_WORK`] and in the browser, no more than there +/// are rows, and no more than [`SCRATCH_BUDGET`] holds (a scratch per point, for millions of +/// points, is large). +pub(crate) fn threads_for_rows(rows: usize, work: usize, scratch_bytes: usize) -> usize { if cfg!(target_arch = "wasm32") || work < MIN_ROW_WORK { return 1; } - cores().min(rows).max(1) + let affordable = 1 + SCRATCH_BUDGET / scratch_bytes.max(1); + cores().min(rows).min(affordable).max(1) } /// The cores to use, at most [`MAX_THREADS`]. @@ -285,14 +290,17 @@ mod tests { } assert!(map_rows(0, 4, || (), |_, r| r).is_empty()); // Small work, and the browser: one thread; never more threads than rows. - assert_eq!(threads_for_rows(100, MIN_ROW_WORK - 1), 1); - assert!(threads_for_rows(3, usize::MAX) <= 3); + assert_eq!(threads_for_rows(100, MIN_ROW_WORK - 1, 0), 1); + assert!(threads_for_rows(3, usize::MAX, 0) <= 3); + // Scratch larger than the budget: one thread more than it holds, none beyond. + assert!(threads_for_rows(100, usize::MAX, SCRATCH_BUDGET) <= 2); + assert_eq!(threads_for_rows(100, usize::MAX, usize::MAX), 1); } #[test] fn a_workers_panic_stops_the_other_threads() { // A worker panics on its first item; the others stop after the items they are making, - // far short of all of them. Unless no worker could start: then nothing panics. + // far short of all of them. let items: Vec = (0..2_000).collect(); let caller = std::thread::current().id(); let ran = std::sync::atomic::AtomicUsize::new(0); @@ -311,9 +319,8 @@ mod tests { }, ) }); - if r.is_err() { - assert!(ran.load(std::sync::atomic::Ordering::Relaxed) < items.len() / 2, "{}", ran.load(std::sync::atomic::Ordering::Relaxed)); - } + assert!(r.is_err(), "no worker ran an item"); + assert!(ran.load(std::sync::atomic::Ordering::Relaxed) < items.len() / 2, "{}", ran.load(std::sync::atomic::Ordering::Relaxed)); } #[test] diff --git a/crates/analysis/src/rasterops.rs b/crates/analysis/src/rasterops.rs index fad4f502c..40627d54e 100644 --- a/crates/analysis/src/rasterops.rs +++ b/crates/analysis/src/rasterops.rs @@ -410,7 +410,9 @@ fn idw_on( let k = if neighbours == 0 { points.len() } else { neighbours.clamp(1, points.len()) }; // Each row on its own, on several threads when the work is large (#1252); a thread's // distances go into its own buffer. - let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, w.saturating_mul(h).saturating_mul(points.len()))); + let scratch_bytes = points.len().saturating_mul(std::mem::size_of::<(f64, f64)>()); + let threads = + threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, w.saturating_mul(h).saturating_mul(points.len()), scratch_bytes)); let rows = crate::parallel::map_rows( h, threads, @@ -447,7 +449,9 @@ fn idw_on( data }, ); - let data = rows.concat(); + // Each row dropped as it is copied, so the surface isn't held twice. + let mut data = Vec::with_capacity(w * h); + rows.into_iter().for_each(|r| data.extend(r)); let gt = GeoTransform { origin_x: extent.min_x, origin_y: extent.max_y, pixel_width: cell_deg, pixel_height: -cell_deg }; Raster::new(name, w as u32, h as u32, vec![data], gt, Crs::LonLat, Some(NODATA as f64), ColorHint::Gray, "Float32") .map_err(|e| e.to_string()) From 6b4afd5b108b299d447a1a7b9d57d258a704fc10 Mon Sep 17 00:00:00 2001 From: Dominik Strebel Date: Sun, 1 Nov 2026 16:31:52 +0100 Subject: [PATCH 3/3] perf(analysis): write IDW and kriging rows where they lie map_rows collected every row before they were copied into the surface, so a raster was held about twice over, and a tall grid of one-cell rows made a vector per row. parallel::for_each_row instead hands each thread its next output row (chunks_mut of the bands, taken under a lock), so the surfaces are held once, as on one thread. A worker's panic stops the others and reaches the caller, as in try_map. try_each_on takes its batch of 16 again, since only try_map_on uses it. --- crates/analysis/src/kriging.rs | 21 +++---- crates/analysis/src/parallel.rs | 97 +++++++++++++++++++++++++------- crates/analysis/src/rasterops.rs | 16 ++---- 3 files changed, 90 insertions(+), 44 deletions(-) diff --git a/crates/analysis/src/kriging.rs b/crates/analysis/src/kriging.rs index c01cd8c7c..d2790e899 100644 --- a/crates/analysis/src/kriging.rs +++ b/crates/analysis/src/kriging.rs @@ -353,11 +353,11 @@ fn kriging_on( let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, w.saturating_mul(h).saturating_mul(work_per_cell), scratch_bytes)); let scratch = || (Vec::<(f64, f64, Coord)>::with_capacity(points.len()), KrigingBuffers::new()); - let rows = crate::parallel::map_rows(h, threads, scratch, |(near, buffers), row| { - let mut pred_data = Vec::with_capacity(w); - let mut var_data = Vec::with_capacity(w); + let (mut pred_data, mut var_data) = (vec![NODATA; w * h], vec![NODATA; w * h]); + let rows = pred_data.chunks_mut(w).zip(var_data.chunks_mut(w)); + crate::parallel::for_each_row(rows, threads, scratch, |(near, buffers), row, (pred, var): (&mut [f32], &mut [f32])| { let kx = (extent.max_y - (row as f64 + 0.5) * cell_deg).to_radians().cos().abs().max(1e-6); - for col in 0..w { + for (col, (pred, var)) in pred.iter_mut().zip(var.iter_mut()).enumerate() { let c = Coord::new(extent.min_x + (col as f64 + 0.5) * cell_deg, extent.max_y - (row as f64 + 0.5) * cell_deg); near.clear(); near.extend(points.iter().map(|(p, v)| { @@ -372,18 +372,11 @@ fn kriging_on( near.truncate(k); } - let (val, var) = ordinary_kriging_point(c, near, params, kx, buffers); - pred_data.push(if val.is_finite() { val as f32 } else { NODATA }); - var_data.push(if var.is_finite() { var as f32 } else { NODATA }); + let (val, var_est) = ordinary_kriging_point(c, near, params, kx, buffers); + *pred = if val.is_finite() { val as f32 } else { NODATA }; + *var = if var_est.is_finite() { var_est as f32 } else { NODATA }; } - (pred_data, var_data) }); - // Each row dropped as it is copied, so the surfaces aren't held twice. - let (mut pred_data, mut var_data) = (Vec::with_capacity(w * h), Vec::with_capacity(w * h)); - for (p, v) in rows { - pred_data.extend(p); - var_data.extend(v); - } let gt = GeoTransform { origin_x: extent.min_x, origin_y: extent.max_y, pixel_width: cell_deg, pixel_height: -cell_deg }; Raster::new(name, w as u32, h as u32, vec![pred_data, var_data], gt, Crs::LonLat, Some(NODATA as f64), ColorHint::Gray, "Float32") diff --git a/crates/analysis/src/parallel.rs b/crates/analysis/src/parallel.rs index 728fe64f4..dd93197c3 100644 --- a/crates/analysis/src/parallel.rs +++ b/crates/analysis/src/parallel.rs @@ -44,19 +44,45 @@ fn cores() -> usize { std::thread::available_parallelism().map_or(1, |n| n.get()).min(MAX_THREADS) } -/// `f` on each of `rows` raster rows by its number, the rows in their order, on `threads` -/// threads ([`threads_for_rows`]). A row is costly enough to be taken one at a time, so rows -/// that cost more than others are shared out too. Each thread makes its own `scratch()` once. -pub(crate) fn map_rows( - rows: usize, +/// `f` on each of a raster's `rows` (its output rows, `chunks_mut` of the bands it writes) with +/// the row's number, on `threads` threads ([`threads_for_rows`]), this one among them. Each row +/// is written where it lies, so the output is held once, as on one thread. A row is costly +/// enough to be taken one at a time, so rows that cost more than others are shared out too. +/// Each thread makes its own `scratch()` once; one that panics stops the others after the rows +/// they are making, and its panic is raised again here. +pub(crate) fn for_each_row( + rows: impl Iterator + Send, threads: usize, scratch: impl Fn() -> S + Sync, - f: impl Fn(&mut S, usize) -> R + Sync, -) -> Vec { - match try_each_on(rows, 1, threads, scratch, |s, i| Ok::(f(s, i))) { - Ok(made) => made, - Err(never) => match never {}, + f: impl Fn(&mut S, usize, T) + Sync, +) { + use std::sync::atomic::{AtomicBool, Ordering}; + if threads <= 1 { + let mut s = scratch(); + rows.enumerate().for_each(|(i, row)| f(&mut s, i, row)); + return; } + let rows = std::sync::Mutex::new(rows.enumerate()); + let failed = AtomicBool::new(false); + let work = || { + let _stop = StopOnPanic(&failed); + let mut s = scratch(); + while !failed.load(Ordering::Relaxed) { + // The lock is held only to take a row; a poisoned one means a thread panicked. + let Some((i, row)) = rows.lock().ok().and_then(|mut r| r.next()) else { break }; + f(&mut s, i, row); + } + }; + std::thread::scope(|scope| { + let others: Vec<_> = + (1..threads).filter_map(|_| std::thread::Builder::new().stack_size(STACK_BYTES).spawn_scoped(scope, work).ok()).collect(); + work(); + for h in others { + if let Err(e) = h.join() { + std::panic::resume_unwind(e); + } + } + }); } /// `f` on each of `items` with its index, the results in their order: on several threads when @@ -84,7 +110,7 @@ pub(crate) fn try_map_on( scratch: impl Fn() -> S + Sync, f: impl Fn(&mut S, usize, &T) -> Result + Sync, ) -> Result, E> { - try_each_on(items.len(), BATCH, threads, scratch, |s, i| f(s, i, &items[i])) + try_each_on(items.len(), threads, scratch, |s, i| f(s, i, &items[i])) } /// Sets the stop flag when the thread it lives on panics, so the others stop after the items @@ -99,10 +125,9 @@ impl Drop for StopOnPanic<'_> { } } -/// [`try_map_on`] over the indices `0..len`, taken `batch` at a time. +/// [`try_map_on`] over the indices `0..len`. fn try_each_on( len: usize, - batch: usize, threads: usize, scratch: impl Fn() -> S + Sync, f: impl Fn(&mut S, usize) -> Result + Sync, @@ -112,7 +137,6 @@ fn try_each_on( let mut s = scratch(); return (0..len).map(|i| f(&mut s, i)).collect(); } - let batch = batch.max(1); let next = AtomicUsize::new(0); let failed = AtomicBool::new(false); let work = || { @@ -120,11 +144,11 @@ fn try_each_on( let mut s = scratch(); let mut made = Vec::new(); while !failed.load(Ordering::Relaxed) { - let from = next.fetch_add(batch, Ordering::Relaxed); + let from = next.fetch_add(BATCH, Ordering::Relaxed); if from >= len { break; } - for i in from..len.min(from.saturating_add(batch)) { + for i in from..len.min(from.saturating_add(BATCH)) { match f(&mut s, i) { Ok(r) => made.push((i, r)), Err(e) => { @@ -283,12 +307,22 @@ mod tests { } #[test] - fn rows_come_back_in_their_order_one_at_a_time() { - let alone: Vec = (0..37).map(|r| r * 3).collect(); + fn each_row_is_written_in_its_place() { + let alone: Vec = (0..37 * 5).map(|i| i / 5 * 3 + i % 5).collect(); for threads in [1, 2, 5, 16, 64] { - assert_eq!(map_rows(37, threads, || (), |_, r| r * 3), alone, "{threads} threads"); + let mut out = vec![0usize; 37 * 5]; + for_each_row( + out.chunks_mut(5), + threads, + || (), + |_, r, row: &mut [usize]| { + row.iter_mut().enumerate().for_each(|(c, o)| *o = r * 3 + c); + }, + ); + assert_eq!(out, alone, "{threads} threads"); } - assert!(map_rows(0, 4, || (), |_, r| r).is_empty()); + let mut none: Vec = Vec::new(); + for_each_row(none.chunks_mut(5), 4, || (), |_, _, _: &mut [u8]| panic!("no rows")); // Small work, and the browser: one thread; never more threads than rows. assert_eq!(threads_for_rows(100, MIN_ROW_WORK - 1, 0), 1); assert!(threads_for_rows(3, usize::MAX, 0) <= 3); @@ -323,6 +357,29 @@ mod tests { assert!(ran.load(std::sync::atomic::Ordering::Relaxed) < items.len() / 2, "{}", ran.load(std::sync::atomic::Ordering::Relaxed)); } + #[test] + fn a_workers_panic_on_a_row_reaches_the_caller_and_stops_the_others() { + let caller = std::thread::current().id(); + let ran = std::sync::atomic::AtomicUsize::new(0); + let mut out = vec![0u8; 2_000]; + let r = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { + for_each_row( + out.chunks_mut(1), + 4, + || (), + |_, _, _: &mut [u8]| { + ran.fetch_add(1, std::sync::atomic::Ordering::Relaxed); + if std::thread::current().id() != caller { + panic!("a worker's panic"); + } + std::thread::sleep(std::time::Duration::from_millis(1)); + }, + ) + })); + assert!(r.is_err(), "no worker took a row"); + assert!(ran.load(std::sync::atomic::Ordering::Relaxed) < 1_000, "{}", ran.load(std::sync::atomic::Ordering::Relaxed)); + } + #[test] fn no_more_threads_than_batches() { assert_eq!(threads_for(10, 256), 1); diff --git a/crates/analysis/src/rasterops.rs b/crates/analysis/src/rasterops.rs index 40627d54e..4357e60fe 100644 --- a/crates/analysis/src/rasterops.rs +++ b/crates/analysis/src/rasterops.rs @@ -413,15 +413,15 @@ fn idw_on( let scratch_bytes = points.len().saturating_mul(std::mem::size_of::<(f64, f64)>()); let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, w.saturating_mul(h).saturating_mul(points.len()), scratch_bytes)); - let rows = crate::parallel::map_rows( - h, + let mut data = vec![NODATA; w * h]; + crate::parallel::for_each_row( + data.chunks_mut(w), threads, || Vec::with_capacity(points.len()), - |near: &mut Vec<(f64, f64)>, row| { - let mut data = Vec::with_capacity(w); + |near: &mut Vec<(f64, f64)>, row, out: &mut [f32]| { // East-west degrees shrink with latitude, which only changes by row. let kx = (extent.max_y - (row as f64 + 0.5) * cell_deg).to_radians().cos(); - for col in 0..w { + for (col, o) in out.iter_mut().enumerate() { let c = Coord::new(extent.min_x + (col as f64 + 0.5) * cell_deg, extent.max_y - (row as f64 + 0.5) * cell_deg); near.clear(); near.extend(points.iter().map(|(p, v)| { @@ -444,14 +444,10 @@ fn idw_on( }); num / den }; - data.push(if v.is_finite() { v as f32 } else { NODATA }); + *o = if v.is_finite() { v as f32 } else { NODATA }; } - data }, ); - // Each row dropped as it is copied, so the surface isn't held twice. - let mut data = Vec::with_capacity(w * h); - rows.into_iter().for_each(|r| data.extend(r)); let gt = GeoTransform { origin_x: extent.min_x, origin_y: extent.max_y, pixel_width: cell_deg, pixel_height: -cell_deg }; Raster::new(name, w as u32, h as u32, vec![data], gt, Crs::LonLat, Some(NODATA as f64), ColorHint::Gray, "Float32") .map_err(|e| e.to_string())