diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index d819f284c..c4e4ccc3a 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -27,7 +27,7 @@ normalize,2026-05-01,,,,rescale and standardize across all 4 backends. NaN/inf f perlin,2026-04-10T12:00:00Z,,,,Improved Perlin noise implementation correct. Fade/gradient functions verified. Backend-consistent. Continuous at cell boundaries. polygon_clip,2026-06-10,3186,HIGH,5,"Cat5 backend inconsistency: dask+cupy clip_polygon rasterizes the mask with a uniform chunk size from the raster's first chunk, then feeds raster+mask to da.map_blocks (positional block pairing). Non-uniform raster chunks gave the mask a different block layout -> IndexError/ValueError (or silent mis-stamp). Repro (8,6) rechunk ((3,5),(6,)) on dask+cupy raised ValueError Shapes do not align; dask+numpy was fine via xarray.where rechunk. Fix #3186/PR: rechunk cond to raster.data.chunks[-2:] before map_blocks; added non-uniform regression tests for dask+numpy and dask+cupy. use_cuda->gpu migration in that branch was already landed by #3089/#3122. CUDA available; cupy+dask+cupy verified, 25 tests pass. Cats 1-4 clean: numpy path uses raster.where, cupy path operates on raw arrays, NaN inputs preserved, no neighborhood ops/curvature. Prior fix #1197/#1200 (crop+all_touched) merged and unrelated." polygonize,2026-05-29,2606,HIGH,5,"Cat 5 HIGH: dask connectivity=8 cross-chunk merge filled diagonal notch where same-value regions meet only at a corner across a chunk boundary; total area exceeded raster. Hole ring was dropped because containment tested hole[0] (on exterior at pinch). Fixed via _ring_interior_point in PR for #2606. numpy, dask+numpy, dask+cupy area parity now holds; 4-conn was already correct. cupy + dask+cupy paths validated on GPU host. Other cats clean: NaN masked on numpy/cupy float paths (tested), _is_close handles +/-inf via exact-equality short-circuit, atol/rtol/simplify_tolerance reject NaN/inf, integer GPU CCL matches numpy." -proximity,2026-06-09,3108,HIGH,4;5,"Cat5/Cat4: bounded GREAT_CIRCLE dask (numpy+cupy) missed targets across the +/-180 antimeridian seam: _halo_depth sized x-halo as linear parallel-arc sum, but haversine is periodic in lon and chords shorten near poles, so array-space adjacency is no lower bound on spherical distance; numpy/cupy (brute force) found the wrap target (~111 km), dask returned NaN. Fixed in #3108 via chord bound 2R asin(cos(lat_max)|sin(dlon/2)|) + x-axis fold when seam/180-deg chord within max_distance (covers over-pole too). CUDA host: cupy + dask+cupy executed, 417+ tests pass. Cat1-3 clean (float32 output documented; NaN via isfinite consistent; bounds guards correct; tie-break unified in #2881). LOW (not fixed): great_circle_distance uses WGS84 equatorial radius 6378137 as sphere radius (~0.1% vs mean-radius convention) but documented and exposed as param." +proximity,2026-06-18,3389,MEDIUM,1;5,"Cat1/Cat5 (MEDIUM, fixed #3389): CPU brute-force _distance casts to float32 before best_dist<=max_distance, but the CUDA kernel compared the float64 distance; a max_distance in a target's float32 ulp gap qualified it on numpy/dask+numpy yet rejected it on cupy/dask+cupy (GREAT_CIRCLE proximity + EUCLIDEAN/MANHATTAN allocation/direction). Fixed by rounding best_dist to float32 in the kernel before the range test. CUDA host: 518 tests pass incl cupy + dask+cupy. Out of scope (separate pre-existing edge, predates fix): bounded dask EUCLIDEAN allocation can still NaN a target exactly at a float32-ulp max_distance via the halo-sizing path. Prior #3108 fix (antimeridian halo) still verified. LOW (unfixed): great_circle_distance uses WGS84 equatorial radius 6378137 as sphere radius (~0.1% vs mean-radius), but documented and param-exposed." rasterize,2026-06-18,3384,HIGH,1;5,dask all_touched polygon-boundary supercover walk re-extracted per tile in tile-local float pixel coords; floor() tie on an on-grid boundary segment diverged from eager (Cat1 precision -> Cat5 backend split). Fixed by extracting boundary float segments in the global grid frame and shifting by integer tile offset (#3384). Verified numpy/cupy/dask+numpy/dask+cupy. reproject,2026-06-12,3274,HIGH,1;4,"3 confirmed bugs, all kernel-vs-PROJ parity: #3274 HIGH LAEA inverse spurious /rq (2.6 km err for 3035) + _authalic_apa inverse series wrong (4.8 m in AEA/CEA inverses; PROJ 3-term = 1.6 mm), CPU+CUDA kernels both; #3275 HIGH _is_wgs84_compatible_ellipsoid passes R-defined spheres (MODIS sinusoidal 18.9 km err) and _aea_params/_cea_params lack the guard entirely (23.8 km on spherical aea/cea); #3276 MEDIUM itrf helmert scale 1e-9 but PROJ +s is ppm (1e-6), ~23 mm err. Verified clean: merc/emerc/UTM/tmerc/LCC/polar stere (incl lat_ts akm1) forward+inverse <=1e-5 m vs pyproj; resampling kernels NaN handling and GDAL renorm match across numpy/cupy/dask (CUDA run, gpu-vs-cpu 1.3e-7); dask footprint chunk-skip bbox is a superset in all probed cases (no holes). LOW (documented only): _source_footprint_in_target probe array typo uses x-midpoint mx as a latitude in last 3 ys entries (bbox superset, correctness unaffected)." resample,2026-05-29,2610,HIGH,3;5,"dask interp (nearest/bilinear) overlap depth=1 too small on downsample; block-centered source coord landed past chunk, map_coordinates clamped to edge -> wrong seam rows. Fixed PR #2627 via per-axis _downsample_radius. cupy+dask+cupy verified." diff --git a/xrspatial/proximity.py b/xrspatial/proximity.py index bbcdf05f9..e383a17f9 100644 --- a/xrspatial/proximity.py +++ b/xrspatial/proximity.py @@ -634,6 +634,12 @@ def _proximity_cuda_kernel(target_xs, target_ys, target_vals, n_targets, best_dist = d best_idx = k + # Round the winning distance to float32 before the range test so the + # in-range decision matches the CPU brute-force path, where _distance + # returns np.float32(d). Comparing the float64 distance here let a target + # whose distance rounds down across a float32 ulp pass on the CPU but fail + # on the GPU (or the reverse) when max_distance sits in that ulp gap. + best_dist = np.float32(best_dist) if best_idx >= 0 and best_dist <= max_distance: if process_mode == PROXIMITY: out[iy, ix] = best_dist diff --git a/xrspatial/tests/test_proximity.py b/xrspatial/tests/test_proximity.py index b7f254dd0..c364e8939 100644 --- a/xrspatial/tests/test_proximity.py +++ b/xrspatial/tests/test_proximity.py @@ -2499,3 +2499,75 @@ def test_proximity_numpy_no_scipy_fallback_is_exact(): expected = _exact_planar_proximity( raster.data, raster['x'].data, raster['y'].data, 'EUCLIDEAN') np.testing.assert_allclose(result.data, expected, rtol=1e-5, atol=1e-7) + + +# --------------------------------------------------------------------------- +# Issue #3389: the max_distance comparison must use the same precision on the +# CPU brute-force and the CUDA kernel. _distance casts to float32 before the +# CPU compares best_dist <= max_distance, while the kernel used to compare the +# float64 distance. A max_distance landing inside the float32 ulp gap of a +# target's distance then qualified the target on numpy but rejected it on cupy +# (or the reverse). +# --------------------------------------------------------------------------- + + +def test_max_distance_float32_boundary_gpu_matches_cpu(): + """GREAT_CIRCLE proximity at a float32-ulp max_distance agrees CPU/GPU.""" + if not has_cuda_and_cupy(): + pytest.skip("Requires CUDA and CuPy") + import cupy + + lon = np.linspace(0, 39, 40) + lat = np.array([0.0]) + data = np.zeros((1, 40), dtype=np.float64) + data[0, 0] = 1.0 + j = 20 + d64 = great_circle_distance(lon[0], lon[j], lat[0], lat[0]) + # max_distance between the float32-rounded and the true float64 distance. + md = (float(np.float32(d64)) + d64) / 2.0 + assert float(np.float32(d64)) < md < d64, "needs a real float32 ulp gap" + + kwargs = dict(x='lon', y='lat', distance_metric='GREAT_CIRCLE', + max_distance=md) + numpy_raster = xr.DataArray( + data, dims=['lat', 'lon'], coords={'lon': lon, 'lat': lat}) + numpy_val = proximity(numpy_raster, **kwargs).data[0, j] + + cupy_raster = xr.DataArray( + cupy.asarray(data), dims=['lat', 'lon'], + coords={'lon': lon, 'lat': lat}) + cupy_val = proximity(cupy_raster, **kwargs).data.get()[0, j] + + # Both keep the target: its float32 distance is <= max_distance. + assert np.isnan(numpy_val) == np.isnan(cupy_val) + np.testing.assert_allclose(numpy_val, cupy_val, rtol=1e-6) + + +def test_max_distance_float32_boundary_allocation_gpu_matches_cpu(): + """EUCLIDEAN allocation (brute-force on both CPU and GPU) agrees at the + float32 ulp boundary.""" + if not has_cuda_and_cupy(): + pytest.skip("Requires CUDA and CuPy") + import cupy + + lon = np.array([0.0, 0.1, 0.7]) + lat = np.array([0.0]) + data = np.zeros((1, 3), dtype=np.float64) + data[0, 0] = 9.0 + j = 2 + d64 = euclidean_distance(lon[0], lon[j], lat[0], lat[0]) + md = (float(np.float32(d64)) + d64) / 2.0 + assert float(np.float32(d64)) < md < d64, "needs a real float32 ulp gap" + + kwargs = dict(x='lon', y='lat', max_distance=md) + numpy_raster = xr.DataArray( + data, dims=['lat', 'lon'], coords={'lon': lon, 'lat': lat}) + numpy_val = allocation(numpy_raster, **kwargs).data[0, j] + + cupy_raster = xr.DataArray( + cupy.asarray(data), dims=['lat', 'lon'], + coords={'lon': lon, 'lat': lat}) + cupy_val = allocation(cupy_raster, **kwargs).data.get()[0, j] + + assert np.isnan(numpy_val) == np.isnan(cupy_val) + np.testing.assert_allclose(numpy_val, cupy_val, rtol=1e-6)