From 914a24ffacaecd6e98f1713dc320e66fda41e18b Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 16 Aug 2026 12:14:00 -0400 Subject: [PATCH] Benchmark the surface_distance paths that were never timed (#3709) The benchmark file covered one of six compute paths, and the path it did cover ran on a source raster where 99.95% of pixels were sources, so the Dijkstra relaxation body never executed. - parameterize cupy and dask+cupy alongside numpy and dask, matching pathfinding.py - replace the near-uniform integer source raster with scattered point sources so the frontier crosses the grid - add time_surface_distance_bounded with a finite max_distance whose pixel radius stays inside one chunk, which routes the dask backends through map_overlap instead of the iterative tile fallback - call .compute() on dask results, following the convention in flood.py, twi.py and interpolate.py - add a numpy-only geodesic class for _dijkstra_geodesic and _precompute_dd_grid Benchmark file only, no source changes. All 50 parameter combinations were run on this host, 13s for a single pass. --- .claude/sweep-benchmarks-state.csv | 1 + benchmarks/benchmarks/surface_distance.py | 89 +++++++++++++++++++++-- 2 files changed, 85 insertions(+), 5 deletions(-) diff --git a/.claude/sweep-benchmarks-state.csv b/.claude/sweep-benchmarks-state.csv index 160a429be..53c389d5c 100644 --- a/.claude/sweep-benchmarks-state.csv +++ b/.claude/sweep-benchmarks-state.csv @@ -2,3 +2,4 @@ module,last_inspected,issue,severity_max,categories_found,notes edge_detection,2026-07-18,3672,MEDIUM,1,"No bench file existed; all 5 public funcs (sobel_x/y, prewitt_x/y, laplacian) uncovered. Compute delegates to convolve_2d (directly benchmarked, but only at 5x5/25x25 kernels, never 3x3, never via the wrappers) so MEDIUM not HIGH. Added benchmarks/benchmarks/edge_detection.py: EdgeDetection class, nx in [300,3000], numpy/cupy/dask, one timing method per public func. All 30 combos executed locally incl. cupy (GPU host). Cat 2/3/4 N/A (module implements no backends itself; no pre-existing bench to be broken)." geotiff,2026-07-02,3603,HIGH,1;2,"No benchmark existed for geotiff; open_geotiff/to_geotiff had zero asv coverage across numpy/dask/cupy. Added benchmarks/benchmarks/geotiff.py: WriteGeoTIFF (numpy/dask/cupy streaming), WriteCOG (numpy/cupy overview pyramid), ReadGeoTIFF (numpy/cupy decode), ReadGeoTIFFChunked (dask). All classes executed locally via direct call; cupy paths run on this GPU host. asv check discover fails suite-wide from an asv_runner + py3.14 metadata bug, unrelated to this file." pathfinding,2026-07-08,3645,HIGH,1;2;3,"Bench covered only numpy a_star_search at nx<=300; module also ships dask (separate sparse-Python A* + LRU chunk cache), cupy fallback, and public multi_stop_search with zero coverage. Extended AStarSearch to numpy/cupy/dask with nx up to 1000 (dask capped at 300, ~4s/call at 1000) and added MultiStopSearch (ordered + optimize_order). All combos executed locally incl. cupy (GPU host). LOW noted, not fixed: open-grid no-barrier/no-friction input is A* best case. dask+cupy not parameterized anywhere in suite (common.get_xr_dataarray has no such type). Existing bench imports/runs fine (Cat 4 clean)." +surface_distance,2026-08-16,3709,HIGH,1;2;3,"Bench (26 lines) covered 1 of 6 compute paths. Cat 2 HIGH: cupy (_sd_relax_kernel) and dask+cupy (_surface_distance_dask_cupy) never parameterized; both run fine on this GPU host. Cat 2 HIGH: default max_distance=inf routed every dask call to the iterative tile fallback, so the bounded map_overlap branch the docstring recommends was never timed, and no .compute() was called (worked only because _sd_dask_iterative computes eagerly). Cat 3 MEDIUM: get_xr_dataarray(is_int=True) made 499769/500000 pixels sources at nx=1000, so the Dijkstra relaxation body never ran; max distance 3.6 vs 160 with sparse sources, 6.4x timing gap. Cat 1 MEDIUM: method='geodesic' (_dijkstra_geodesic + _precompute_dd_grid) unbenchmarked, ~2.2x planar cost. Cat 4 clean: file imported and ran. Fixed: 4 backends, sparse point sources, bounded method, .compute(), numpy-only geodesic class. All 50 combos executed locally, 13s single pass. LOW not fixed: nx=100 gives a 50x100 grid. Note: asv discover is broken suite-wide in this conda env (asv_runner dist-metadata bug on py3.14), so verification was by direct class invocation." diff --git a/benchmarks/benchmarks/surface_distance.py b/benchmarks/benchmarks/surface_distance.py index 1cf7aad60..d2fd5fdac 100644 --- a/benchmarks/benchmarks/surface_distance.py +++ b/benchmarks/benchmarks/surface_distance.py @@ -1,26 +1,105 @@ import numpy as np +import xarray as xr from xrspatial.surface_distance import ( surface_distance, surface_allocation, surface_direction, ) +from xrspatial.utils import has_cuda_and_cupy from .common import get_xr_dataarray +def _sparse_source_raster(ny, nx, type): + """Source raster with a handful of scattered target pixels. + + ``get_xr_dataarray(is_int=True)`` draws integers over ``[-nx, nx)``, + and surface_distance treats every non-zero finite pixel as a source, + so all but roughly 1 in ``2 * nx`` pixels seed the search at distance + zero. The Dijkstra relaxation body then never runs and the benchmark + times a heap drain instead of distance propagation. Scattered point + sources make the frontier cross the whole grid, which is what these + functions are for. + """ + rng = np.random.default_rng(71942) + z = np.zeros((ny, nx), dtype=np.float32) + n_sources = max(4, (ny * nx) // 20000) + rows = rng.integers(0, ny, n_sources) + cols = rng.integers(0, nx, n_sources) + z[rows, cols] = np.arange(1, n_sources + 1, dtype=np.float32) + + chunks = (max(1, ny // 2), max(1, nx // 2)) + if type == "cupy": + if not has_cuda_and_cupy(): + raise NotImplementedError() + import cupy + z = cupy.asarray(z) + elif type == "dask": + import dask.array as da + z = da.from_array(z, chunks=chunks) + elif type == "dask+cupy": + if not has_cuda_and_cupy(): + raise NotImplementedError() + import cupy + import dask.array as da + z = da.from_array(cupy.asarray(z), chunks=chunks) + elif type != "numpy": + raise RuntimeError(f"Unrecognised type {type}") + + y = np.linspace(-90, 90, ny) + x = np.linspace(-180, 180, nx) + return xr.DataArray(z, coords=dict(y=y, x=x), dims=["y", "x"]) + + +def _compute(result): + if hasattr(result.data, "compute"): + result.data.compute() + + class SurfaceDistance: - params = ([100, 300, 1000], ["numpy", "dask"]) + params = ([100, 300, 1000], ["numpy", "cupy", "dask", "dask+cupy"]) param_names = ("nx", "type") def setup(self, nx, type): ny = nx // 2 - self.agg = get_xr_dataarray((ny, nx), type, is_int=True) + self.agg = _sparse_source_raster(ny, nx, type) self.elev = get_xr_dataarray((ny, nx), type) + # A finite max_distance whose pixel radius stays inside one chunk + # (chunks are ny//2 x nx//2) routes the dask backends through the + # bounded map_overlap branch instead of the iterative tile one. + cellsize = min(360.0 / (nx - 1), 180.0 / (ny - 1)) + self.max_distance = 20 * cellsize + def time_surface_distance(self, nx, type): - surface_distance(self.agg, self.elev) + _compute(surface_distance(self.agg, self.elev)) + + def time_surface_distance_bounded(self, nx, type): + _compute(surface_distance( + self.agg, self.elev, max_distance=self.max_distance)) def time_surface_allocation(self, nx, type): - surface_allocation(self.agg, self.elev) + _compute(surface_allocation(self.agg, self.elev)) def time_surface_direction(self, nx, type): - surface_direction(self.agg, self.elev) + _compute(surface_direction(self.agg, self.elev)) + + +class SurfaceDistanceGeodesic: + """Great-circle horizontal distances from lat/lon coordinates. + + Runs a separate numba kernel (``_dijkstra_geodesic``) behind a + precomputed per-pixel neighbour-distance grid, and costs about twice + the planar path. numpy only: the module raises NotImplementedError + for geodesic on cupy, dask, and dask+cupy. + """ + + params = [300, 1000] + param_names = ("nx",) + + def setup(self, nx): + ny = nx // 2 + self.agg = _sparse_source_raster(ny, nx, "numpy") + self.elev = get_xr_dataarray((ny, nx), "numpy") + + def time_surface_distance_geodesic(self, nx): + surface_distance(self.agg, self.elev, method="geodesic")