feat(analysis): Viewshed, IDW and Kriging take as many times the work as the threads they run on (#1291 step 2)

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));

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<T: Sync, S, R: Send>(
         items: &[T],

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<usize>, 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,19 +468,23 @@ 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());
     }
     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 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]

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`).