Skip to content

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

Open
brendancol wants to merge 3 commits into
mainfrom
issue-3689
Open

Compare nearest-target distances at float32 on every backend (#3689)#3690
brendancol wants to merge 3 commits into
mainfrom
issue-3689

Conversation

@brendancol

Copy link
Copy Markdown
Contributor

Closes #3689

allocation and direction evaluated "nearest target" at a different float precision per backend: the numpy brute force compares float32-rounded distances, the CUDA kernel compared float64, and the cKDTree paths detected ties with exact float64 equality. On grids whose coordinates are not exact float64 lattice points (res 0.1, np.linspace, reprojected coords), a geometric tie shows up as float64 distances separated by ~1e-13 of rounding noise that vanishes at float32. The numpy backend called it a tie and applied the documented lowest-flat-index rule; cupy and the dask cKDTree paths each returned whichever target their own float64 arithmetic favored, so three backends could produce three different rasters. proximity itself was unaffected (the output distance is float32 either way).

Changes:

  • _proximity_cuda_kernel rounds each candidate distance to float32 inside the argmin loop, so the winning target (and the existing float32 max_distance test from proximity: max_distance comparison precision differs between CPU and GPU backends #3389) matches the CPU brute force.
  • _kdtree_query_lowest_index detects ties at float32 instead of exact float64 equality, covering both the global and tiled dask KDTree paths.
  • The Tie-breaking paragraphs in the allocation/direction docstrings now state that distances are compared at float32 precision, the precision of the output.

Backend coverage: numpy (already correct, unchanged), cupy, dask+numpy (bounded map_overlap path was already correct via the brute-force chunk kernel; unbounded KDTree path fixed), dask+cupy (bounded GPU kernel fixed; unbounded routes through the fixed KDTree path).

Also updates the accuracy-sweep state CSV row for proximity.

Test plan:

  • New test_tie_break_float32_precision_nonlattice_grid: allocation and direction on a 40x40 res-0.1 grid, all four backends, bounded and unbounded max_distance, asserting float32-tied pixels resolve to the lowest flat index everywhere. 10 of 16 parametrizations fail without the fix.
  • Full proximity suite: 595 passed locally with CUDA (cupy and dask+cupy paths executed, not skipped).
  • Reference checks from the sweep: results match scipy.ndimage.distance_transform_edt and a float64 brute force for all three metrics and modes.

The CUDA kernel picked the nearest target by comparing float64 distances
and the cKDTree paths detected ties with exact float64 equality, while
the numpy brute force compares float32-rounded distances. On grids whose
coordinates are not exact float64 lattice points (res 0.1, linspace,
reprojected coords), geometric ties differ by float64 rounding noise but
round equal at float32, so allocation and direction diverged across
backends at bisector pixels.

Round each candidate distance to float32 in the CUDA argmin loop and
detect cKDTree ties at float32, so every backend applies the documented
lowest-flat-index tie-break to the same set of ties. Document the
float32 tie precision on allocation/direction.
@github-actions github-actions Bot added the performance PR touches performance-sensitive code label Jul 22, 2026

@brendancol brendancol left a comment

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Blockers (must fix before merge)

None.

Suggestions (should fix, not blocking)

  • xrspatial/tests/test_proximity.py:705-706: the fixture sanity check computes distances with np.hypot, but the implementation computes np.sqrt(x*x + y*y) (euclidean_distance, proximity.py:106-108). The two can differ by one float64 ulp, and this test lives exactly in the regime where one ulp decides whether a pixel counts as float32-tied. Use the same sqrt(dx*dx + dy*dy) formula so the oracle and the code under test agree bit-for-bit on every platform.

Nits (optional improvements)

  • xrspatial/proximity.py:863-868: the updated 3+-way-tie caveat in _kdtree_query_lowest_index is honest about the k=2 limitation, but now that ties are float32-wide it would be worth stating the practical consequence in one sentence: a 3-way float32 tie whose lowest-flat-index member is the float64-largest resolves differently on the KDTree path than on the kernels. Rare (three targets within one float32 ulp), but the docstring is the only place a future reader will learn it.

What looks good

  • The three backends genuinely disagreed before this change (issue #3689 repro shows numpy, cupy, and dask+numpy each returning a different raster), and the fix picks the only convention all of them can implement exactly: float32, the output precision.
  • The CUDA change keeps the strict < so the lowest flat index wins a tie without any extra bookkeeping, and the float32 rounding now covers both the argmin and the existing #3389 range test with one cast.
  • isfinite still gates the KDTree tie test, so inf sentinels from bounded queries cannot alias into float32 ties.
  • The regression test self-validates its geometry (asserts at least one float32-tied pixel whose float64 noise favors the higher index) instead of hard-coding pixel positions, and 10 of its 16 parametrizations fail on the pre-fix code.
  • Bounded dask paths were left alone, correctly: their chunk function is the brute-force kernel that already compared float32.

Checklist

  • Algorithm matches reference/paper (float32 tie convention now documented in both docstrings)
  • All implemented backends produce consistent results (verified on CUDA host: numpy, cupy, dask+numpy, dask+cupy)
  • NaN handling is correct (untouched by this change; isfinite gates preserved)
  • Edge cases are covered by tests (bounded and unbounded max_distance; pre-fix failure confirmed)
  • Dask chunk boundaries handled correctly (bounded map_overlap path unchanged, chunks=(16,16) exercised)
  • No premature materialization or unnecessary copies (two .astype(np.float32) on k=2 query columns; negligible)
  • Benchmark exists or is not needed (behavior fix; one extra cast per candidate in the GPU loop is noise)
  • README feature matrix updated (not applicable, no API change)
  • Docstrings present and accurate (Tie-breaking paragraphs updated in allocation and direction)

@brendancol brendancol left a comment

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PR Review follow-up: commit eabb349

Both findings from the first pass are addressed:

  • The test oracle now computes np.sqrt(dx*dx + dy*dy), the same expression euclidean_distance uses, with a comment explaining why np.hypot was not safe here. The fixture sanity check and the code under test now agree bit-for-bit.
  • _kdtree_query_lowest_index spells out the 3-way float32 tie consequence and what it takes to hit it (three targets within one float32 ulp of the same distance).

No new findings. Tie-break tests pass (30/30 locally, CUDA backends included).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

performance PR touches performance-sensitive code

Projects

None yet

Development

Successfully merging this pull request may close these issues.

allocation/direction: nearest-target tie evaluated at different float precision per backend, diverging on non-lattice grids

1 participant