diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 5f6f509f2..6d2b39673 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.claude/sweep-accuracy-state.csv @@ -28,7 +28,7 @@ pathfinding,2026-07-08,3646;3647,HIGH,2,"Re-audit after 2026-07-03 sweep fixes ( 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-22,3442,MEDIUM,5,"Cat5 backend divergence: EUCLIDEAN proximity/allocation/direction on numpy and dask+numpy use cKDTree, which dropped a target sitting exactly at max_distance to NaN (clearest at max_distance=0: a target pixel is distance 0 from itself), while cupy/dask+cupy brute-force kept it (0.0 / target value / 0). Root cause: cKDTree squares its exclusive distance_upper_bound for p=2, and the old np.nextafter(0,inf) widen underflows when squared, collapsing the bound back to max_distance; _kdtree_query_lowest_index passed the raw exclusive bound with no widen at all. Fix #3442/PR #3443: route all 3 cKDTree call sites through one _inclusive_upper_bound helper (relative eps + absolute floor, robust for p=1 and p=2). MANHATTAN (p=1) was unaffected. Verified all 4 backends agree; 543 proximity tests pass incl. 25 new cross-backend regression cases. CUDA available; cupy + dask+cupy executed. Cats 1-4 clean: float32/float64 boundary already handled (#3389), antimeridian halo (#3108), irregular-coord halo (#2908), non-finite target_values rejected (#2868), non-monotonic coords rejected (#2875), tie-break row-major across backends (#3090/#2881)." +proximity,2026-07-22,3689,MEDIUM,5,"Cat5 backend divergence: allocation/direction argmin compared distances at different precision per backend: numpy brute force float32 (_distance rounds), CUDA kernel float64, cKDTree tie detection exact float64. On non-lattice coord grids (res 0.1/linspace/reprojected) geometric ties differ by ~1e-13 float64 noise but round equal at float32, so numpy applied the lowest-flat-index rule while cupy and dask kdtree each picked their own float64-closer target; 92/300 random 2-target layouts on a 40x40 res-0.1 grid had >=1 diverging pixel, and cupy vs cKDTree also disagreed with each other (squared vs rooted float64 orderings). proximity output unaffected (float32 either way). Fix #3689: round candidates to float32 in the CUDA argmin loop, detect cKDTree ties at float32, document float32 tie precision in both docstrings; 16-case cross-backend regression test on a res-0.1 grid (10 failed pre-fix). CUDA available; cupy + dask+cupy executed, 595 proximity tests pass. Cats 1-4,6 clean: matches scipy.ndimage.distance_transform_edt exactly and float64 brute force for all 3 metrics/modes/bounded+unbounded; target-pixel/translation invariants hold. Dedup honored: #3443 inclusive bound, #3392/#3555 f32-ulp bound widen, #3389 range-test precision, #3090 tie order. Test-coverage note: all prior tie-break tests used integer-lattice coords, which is why this class never fired." 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 c560b776d..bd54a5696 100644 --- a/xrspatial/proximity.py +++ b/xrspatial/proximity.py @@ -628,18 +628,21 @@ def _proximity_cuda_kernel(target_xs, target_ys, target_vals, n_targets, best_dist = 1.0e38 best_idx = -1 + # Round each candidate distance to float32 before comparing, matching the + # CPU brute-force path where _distance returns np.float32(d). This makes + # both the argmin and the max_distance range test float32 on every + # backend: two targets whose float64 distances differ only past the + # float32 mantissa are a tie, and the strict < keeps the first (lowest + # flat-index) target, the tie-break documented on allocation/direction. + # Comparing float64 here instead let the float64-closer target win on the + # GPU while the CPU called it a tie and picked the lowest index. for k in range(n_targets): - d = _gpu_distance(px, target_xs[k], py, target_ys[k], distance_metric) + d = np.float32(_gpu_distance( + px, target_xs[k], py, target_ys[k], distance_metric)) if d < best_dist: 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 @@ -852,10 +855,21 @@ def _kdtree_query_lowest_index(tree, query_pts, p, max_distance): the tie-break policy documented on ``allocation``/``direction``. Query the two nearest targets; wherever they are equidistant, keep the one - with the smaller index. This resolves 2-way ties, which is what grid - geometry produces in practice. A pixel equidistant to three or more targets - relies on cKDTree returning the lower index among the rest, which it does - for the row-major target order used here but does not strictly promise. + with the smaller index. "Equidistant" is evaluated at float32 precision, + matching the brute-force and CUDA kernels, which compare float32-rounded + distances (see ``_distance``). Detecting ties with exact float64 equality + instead let this path pick the float64-closer target on grids whose + coordinates are not exact float64 lattice points (e.g. res 0.1), where a + geometric tie shows up as float64 distances separated by rounding noise + but identical in float32 -- diverging from the kernels' lowest-index pick. + This resolves 2-way ties, which is what grid geometry produces in + practice. A pixel tied with three or more targets relies on the two + float64-smallest of them (the pair cKDTree returns) including the + lowest-index one, which holds for exact ties under the row-major target + order used here but is not strictly promised. Concretely: a 3-way float32 + tie whose lowest-flat-index member is the float64-largest of the three + resolves to a different target here than on the brute-force/CUDA kernels. + That takes three targets within one float32 ulp of the same distance. """ # Match the inclusive dist <= max_distance check of the brute-force/CUDA # backends; cKDTree's bound is exclusive (see _inclusive_upper_bound). @@ -869,10 +883,12 @@ def _kdtree_query_lowest_index(tree, query_pts, p, max_distance): distance_upper_bound=upper) dists = dists2[:, 0] indices = idx2[:, 0] - # A tie exists where both neighbours are finite and equidistant. Prefer the - # smaller index in that case so the result is independent of cKDTree's - # internal traversal order. - tied = np.isfinite(dists2[:, 1]) & (dists2[:, 1] == dists) + # A tie exists where both neighbours are finite and equidistant at float32 + # precision. Prefer the smaller index in that case so the result is + # independent of cKDTree's internal traversal order and of float64 + # rounding noise the kernels' float32 comparison cannot see. + tied = (np.isfinite(dists2[:, 1]) + & (dists2[:, 1].astype(np.float32) == dists.astype(np.float32))) if tied.any(): indices = np.where(tied, np.minimum(idx2[:, 0], idx2[:, 1]), indices) return dists, indices @@ -1777,12 +1793,14 @@ def allocation( expanding each chunk's borders to cover `max_distance`; otherwise the nearest targets are found with a KDTree query over all target pixels. - Tie-breaking: when two or more targets are exactly equidistant from a - pixel, the target with the lowest flat (row-major) index wins, i.e. the - first target encountered when scanning the raster top-to-bottom and - left-to-right. This policy is identical across all backends (numpy, cupy, - dask+numpy, dask+cupy), so the allocated value is deterministic regardless - of which backend computes it. + Tie-breaking: when two or more targets are equidistant from a pixel, the + target with the lowest flat (row-major) index wins, i.e. the first target + encountered when scanning the raster top-to-bottom and left-to-right. + Distances are compared at float32 precision (the precision of the + output), so targets whose distances differ only beyond the float32 + mantissa count as equidistant. This policy is identical across all + backends (numpy, cupy, dask+numpy, dask+cupy), so the allocated value is + deterministic regardless of which backend computes it. Parameters ---------- @@ -1939,12 +1957,15 @@ def direction( expanding each chunk's borders to cover `max_distance`; otherwise the nearest targets are found with a KDTree query over all target pixels. - Tie-breaking: when two or more targets are exactly equidistant from a - pixel, the direction is computed toward the target with the lowest flat - (row-major) index, i.e. the first target encountered when scanning the - raster top-to-bottom and left-to-right. This policy is identical across - all backends (numpy, cupy, dask+numpy, dask+cupy), so the reported - direction is deterministic regardless of which backend computes it. + Tie-breaking: when two or more targets are equidistant from a pixel, the + direction is computed toward the target with the lowest flat (row-major) + index, i.e. the first target encountered when scanning the raster + top-to-bottom and left-to-right. Distances are compared at float32 + precision (the precision of the output), so targets whose distances + differ only beyond the float32 mantissa count as equidistant. This policy + is identical across all backends (numpy, cupy, dask+numpy, dask+cupy), so + the reported direction is deterministic regardless of which backend + computes it. Parameters ---------- diff --git a/xrspatial/tests/test_proximity.py b/xrspatial/tests/test_proximity.py index 3bb31502f..d065b4047 100644 --- a/xrspatial/tests/test_proximity.py +++ b/xrspatial/tests/test_proximity.py @@ -670,6 +670,65 @@ def test_direction_tie_break_chunk_column_order( general_output_checks(raster, result, expected) +@pytest.fixture +def tie_break_nonlattice_raster_data(): + # Regression data for issue #3689. On a res-0.1 grid the linspace + # coordinates are not exact float64 lattice points, so a pixel that is + # geometrically equidistant to two targets sees float64 distances + # separated by rounding noise (~1e-13) that vanishes at float32, the + # precision of the output. The tie must be detected at float32 and + # resolved to the lowest flat (row-major) index on every backend. The + # numpy brute force already compared float32 distances, but the CUDA + # kernel compared float64 and the cKDTree paths detected ties with exact + # float64 equality, so each backend allocated such pixels to whichever + # target its own float64 arithmetic favoured. + data = np.zeros((40, 40), dtype=np.float64) + data[27, 35] = 1.0 # target A, flat index 1115 + data[37, 25] = 2.0 # target B, flat index 1505 + return data + + +@pytest.mark.parametrize("backend", ['numpy', 'dask+numpy', 'cupy', 'dask+cupy']) +@pytest.mark.parametrize("max_distance", [np.inf, 3.5]) +@pytest.mark.parametrize("func", [allocation, direction]) +def test_tie_break_float32_precision_nonlattice_grid( + backend, max_distance, func, tie_break_nonlattice_raster_data): + data = tie_break_nonlattice_raster_data + res_attrs = {'res': (0.1, 0.1)} + + numpy_raster = create_test_raster(data, backend='numpy', attrs=res_attrs) + + # Sanity-check the fixture geometry: at least one pixel must be tied at + # float32 while its float64 distances favour target B (the higher flat + # index). Those are the pixels where a float64 comparison diverges from + # the documented float32 tie-break. + # Use the implementation's exact formula (sqrt of the sum of squares, see + # euclidean_distance), not np.hypot: the two can differ by one float64 + # ulp, which is enough to flip a float32 tie in this regime. + xs = numpy_raster['x'].data + ys = numpy_raster['y'].data + dx_a, dy_a = xs[None, :] - xs[35], ys[:, None] - ys[27] + dx_b, dy_b = xs[None, :] - xs[25], ys[:, None] - ys[37] + d_a = np.sqrt(dx_a * dx_a + dy_a * dy_a) + d_b = np.sqrt(dx_b * dx_b + dy_b * dy_b) + tied_f32 = np.float32(d_a) == np.float32(d_b) + assert (tied_f32 & (d_b < d_a) & (d_a <= max_distance)).any() + + expected = func(numpy_raster, max_distance=max_distance).data + + # Every float32-tied pixel within range resolves to target A, the lowest + # flat index: allocation returns A's value; direction returns the angle + # toward A, which differs from the angle toward B by tens of degrees. + in_range = tied_f32 & (np.float32(d_a) <= max_distance) + if func is allocation: + assert (expected[in_range] == 1.0).all() + + raster = create_test_raster( + data, backend=backend, attrs=res_attrs, chunks=(16, 16)) + result = func(raster, max_distance=max_distance) + general_output_checks(raster, result, expected) + + @pytest.mark.skipif(da is None, reason="dask is not installed") def test_proximity_dask_kdtree_no_targets(): """No target pixels found → result is all NaN."""