Compare nearest-target distances at float32 on every backend (#3689)

.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."

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

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."""