diff --git a/crates/analysis/src/kriging.rs b/crates/analysis/src/kriging.rs index c175ba1a4..271dc9bcc 100644 --- a/crates/analysis/src/kriging.rs +++ b/crates/analysis/src/kriging.rs @@ -14,7 +14,8 @@ use rustgis_geom::{Coord, Extent}; use crate::rasterops::{NODATA, max_pixels}; -/// Maximum cell-point pairs Kriging computes (about a second's work). +/// Maximum cell-point pairs Kriging computes on one thread (about a second's work); a run on several +/// threads may take as many times more ([`crate::parallel::work_cap`]). pub const MAX_KRIGING_WORK: usize = 40_000_000; /// Theoretical semivariogram models for spatial autocorrelation. @@ -343,15 +344,19 @@ fn kriging_on( let k = if neighbours == 0 { 12.min(points.len()) } else { neighbours.clamp(2, points.len().min(32)) }; let k3 = (k * k * k) / 3; let work_per_cell = points.len().saturating_add(k3); - if w.saturating_mul(h).saturating_mul(work_per_cell) > MAX_KRIGING_WORK { - return Err(format!("{w} × {h} cells from {} points is too much work; use a larger cell size", points.len())); - } + let work = w.saturating_mul(h).saturating_mul(work_per_cell); - // Each row on its own, on several threads when the work is large (#1252); a thread has its - // own distances and system buffers. + // Each row on its own, on several threads when the work is large (#1252), and as many times + // the work (#1291); a thread has its own distances and system buffers. 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 threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, work, scratch_bytes)); + let cap = crate::parallel::work_cap(MAX_KRIGING_WORK, threads); + if work > cap { + return Err(format!( + "{w} × {h} cells from {} points is too much work, more than the {cap} this computer takes on {threads} thread(s); use a larger cell size", + points.len() + )); + } let scratch = || (Vec::<(f64, f64, Coord)>::with_capacity(points.len()), KrigingBuffers::new()); let (mut pred_data, mut var_data) = (crate::rasterops::new_band(w, h, NODATA)?, crate::rasterops::new_band(w, h, NODATA)?); let rows = pred_data.chunks_mut(w).zip(var_data.chunks_mut(w)); diff --git a/crates/analysis/src/parallel.rs b/crates/analysis/src/parallel.rs index 85408220b..b86a37cb5 100644 --- a/crates/analysis/src/parallel.rs +++ b/crates/analysis/src/parallel.rs @@ -42,6 +42,13 @@ pub(crate) fn threads_for_rows(rows: usize, work: usize, scratch_bytes: usize) - cores().min(rows).min(affordable).max(1) } +/// A work cap tuned for one thread (a second or two of its work) for a run on `threads` threads: +/// as many times larger, since the threads share the work and the run takes about as long +/// (#1291 step 2). On one thread, and so in the browser, the cap as it is. +pub(crate) fn work_cap(one_thread: usize, threads: usize) -> usize { + one_thread.saturating_mul(threads.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) @@ -225,6 +232,14 @@ impl Drop for StopOnPanic<'_> { mod tests { use super::*; + #[test] + fn a_work_cap_grows_with_the_threads_a_run_takes() { + assert_eq!(work_cap(400, 1), 400); + assert_eq!(work_cap(400, 16), 6400); + assert_eq!(work_cap(400, 0), 400, "no threads counts as one"); + assert_eq!(work_cap(usize::MAX, 2), usize::MAX, "no overflow"); + } + /// The infallible case, for the tests. fn map_on( items: &[T], diff --git a/crates/analysis/src/rasterops.rs b/crates/analysis/src/rasterops.rs index 6764dae91..ffe5d9965 100644 --- a/crates/analysis/src/rasterops.rs +++ b/crates/analysis/src/rasterops.rs @@ -350,9 +350,10 @@ pub fn line_of_sight( Some((target_visible, out)) } -/// Most ray steps a viewshed takes, and cell-point pairs an interpolation weighs (each about a -/// second or two). -const MAX_VIEWSHED_WORK: i64 = 400_000_000; +/// Most ray steps a viewshed takes, and cell-point pairs an interpolation weighs, on one thread +/// (each about a second or two); a run on several threads may take as many times more +/// ([`crate::parallel::work_cap`]). +const MAX_VIEWSHED_WORK: usize = 400_000_000; pub(crate) const MAX_IDW_WORK: usize = 400_000_000; /// Cells visible from an observer (1) or not (0) within `radius_m`; outside the radius no-data. @@ -376,18 +377,24 @@ fn viewshed_on(threads: Option, r: &Raster, b: usize, observer: Coord, he let rad_r = (radius_m / cy).ceil() as i64; // Each cell walks a ray back to the observer: bound cells × ray length, not just cells. let cells = (2 * rad_c + 1).saturating_mul(2 * rad_r + 1); - let ray_work = cells.saturating_mul(rad_c.max(rad_r) / 2 + 1); - if cells > 25_000_000 || ray_work > MAX_VIEWSHED_WORK { + let ray_work = usize::try_from(cells.saturating_mul(rad_c.max(rad_r) / 2 + 1)).unwrap_or(usize::MAX); + if cells > 25_000_000 { return Err("the radius covers too many cells; use a smaller radius".into()); } - let w = r.width as usize; - let mut out = new_band(w, r.height as usize, NODATA)?; - const R: f64 = 6_371_000.0 * 4.0 / 3.0; // Each cell walks its own ray and writes only itself, so the rows within the radius are made - // on several threads when the work is large (#1252). + // on several threads when the work is large (#1252), and may take as many times the work. let (first, last) = ((or - rad_r).max(0), (or + rad_r).min(r.height as i64 - 1)); let rows = usize::try_from(last - first + 1).unwrap_or(0); - let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(rows, usize::try_from(ray_work).unwrap_or(usize::MAX), 0)); + let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(rows, ray_work, 0)); + let cap = crate::parallel::work_cap(MAX_VIEWSHED_WORK, threads); + if ray_work > cap { + return Err(format!( + "the radius takes {ray_work} ray steps, more than the {cap} this computer takes on {threads} thread(s); use a smaller radius" + )); + } + let w = r.width as usize; + let mut out = new_band(w, r.height as usize, NODATA)?; + const R: f64 = 6_371_000.0 * 4.0 / 3.0; let window = out.chunks_mut(w).skip(usize::try_from(first).unwrap_or(0)).take(rows); crate::control::begin(rows); crate::parallel::for_each_row( @@ -461,10 +468,6 @@ fn idw_on( if w.saturating_mul(h) > max_pixels() / 4 { return Err(format!("{w} × {h} cells is too many, more than this computer's {}; use a larger cell size", max_pixels() / 4)); } - // Every cell weighs every point. - if w.saturating_mul(h).saturating_mul(points.len()) > MAX_IDW_WORK { - return Err(format!("{w} × {h} cells from {} points is too much work; use a larger cell size", points.len())); - } if !(power.is_finite() && power >= 0.0) { return Err("the distance coefficient must be zero or more".into()); } @@ -472,8 +475,16 @@ fn idw_on( // Each row on its own, on several threads when the work is large (#1252); a thread's // distances go into its own buffer. 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)); + // Every cell weighs every point: as many times the work as there are threads (#1291). + let work = w.saturating_mul(h).saturating_mul(points.len()); + let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_rows(h, work, scratch_bytes)); + let cap = crate::parallel::work_cap(MAX_IDW_WORK, threads); + if work > cap { + return Err(format!( + "{w} × {h} cells from {} points is too much work, more than the {cap} cell-point pairs this computer takes on {threads} thread(s); use a larger cell size", + points.len() + )); + } let mut data = new_band(w, h, NODATA)?; crate::control::begin(h); crate::parallel::for_each_row( @@ -1212,6 +1223,19 @@ mod tests { } } + /// A run past one thread's work, on one thread, is refused, naming the cap and the threads. + #[test] + fn idw_past_one_threads_work_is_refused_on_one_naming_the_cap() { + // 1000 × 1000 cells from 500 points: 500 million pairs, past one thread's 400 million. + let pts: Vec<(Coord, f64)> = (0..500).map(|i| (Coord::new((i % 25) as f64 / 2.5, (i / 25) as f64 / 2.0), i as f64)).collect(); + let extent = Extent::new(0.0, 0.0, 10.0, 10.0); + let e = idw_on(Some(1), &pts, &extent, 0.01, 2.0, 0, "idw").unwrap_err(); + assert!( + e.contains("too much work") && e.contains(&format!("{MAX_IDW_WORK} cell-point pairs")) && e.contains("on 1 thread(s)"), + "{e}" + ); + } + /// 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] diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index dd130de68..b8320385f 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -70,6 +70,15 @@ switching to another map, a result that adds to or changes the map is refused (o tables or a file is still made). Opening another project drops the run. Scripts, models, the AI Assistant, the command line and the web app wait for each tool to finish. +**Viewshed**, **IDW Interpolation** and **Kriging** refuse a run that would take more than a second +or two of work on one thread (a ray step, a cell-point pair). On the desktop they run on several +cores (up to 16), and they take as many times more work as the threads they run on. Threads share +cores and memory, so a run at that limit can take up to about twice as long as a one-thread run at +its own; in the web app, on one thread, the limits stay as they are. The error names the +limit on this computer and the threads it counted. A tool's own default cell size still follows the +one-thread limit, so it makes the same grid on any computer; a finer cell size given is what the +larger limit lets through. + From the command line or a script (`gp.run`), a built-in tool's parameter with a list of choices, such as Buffer's units, takes one of them in any letter case (`"meters"` for `Meters`); any other value is an error that names the choices (Buffer's units no longer take the abbreviations `m`, `mi` or `metres`).