Honor the boundary parameter on the cupy backend in bilateral()

.claude/sweep-accuracy-state.csv

@@ -1,7 +1,7 @@
 module,last_inspected,issue,severity_max,categories_found,notes
 aspect,2026-06-02,2827,MEDIUM,5,"Cat5 backend divergence: planar cupy _gpu snapped aspect>359.999 to 0 (no such clamp in numpy _cpu, whose range is [0,360) and never reaches 360), so cupy/dask+cupy disagreed with numpy by ~360 on near-degenerate gradients (gx~0+, gy>0). Removing the clamp exposed a 2nd divergence: GPU used coarse 57.29578 vs numpy 180/pi, flipping the >90 compass branch and yielding exact 360 vs 0 on uint32/uint64 random data. Fix #2827/PR #2833: GPU reuses RADIAN and wraps >=360 back to [0,360). Cats 1-4 clean; geodesic path canonicalizes consistently on CPU+GPU and was left untouched. CUDA available; cupy+dask+cupy verified (235 tests pass, numpy-vs-cupy max abs diff 0 over 360 rasters). Dedup: prior aspect fixes #2780 (cellsize)/#2774 (dask mem guard)/#2781 (oracle) all merged and unrelated. Note: PR review COMMENT could not be posted to GitHub (auto-mode permission denial); findings recorded in PR run instead."
 balanced_allocation,2026-04-14T12:00:00Z,1203,,,float32 allocation array caused source ID mismatch for non-integer IDs. Fix in PR #1205.
-bilateral,2026-05-01,,,,"No CRIT/HIGH/MEDIUM. Sigma underflow validated via sqrt(tiny) bound; oversize sigma clamped. float64 throughout numpy/cupy. NaN center returns NaN; NaN neighbors skipped (denom not incremented). w_sum>0 guard avoids div-by-zero. map_overlap depth==kernel radius. CUDA bounds correct. Inf input could yield 0*inf=NaN in v_sum but unvalidated input is general xrspatial pattern, not bilateral-specific."
+bilateral,2026-07-03,#3625,HIGH,backend-inconsistency,"HIGH fixed: cupy backend silently ignored boundary param (nearest/reflect/wrap behaved as nan), diverging from numpy/dask+numpy/dask+cupy; issue #3625, fixed via _bilateral_cupy_boundary wrapper + GPU boundary parity test. assert_boundary_mode_correctness only covers numpy vs dask+numpy (coverage gap noted for test-coverage sweep). Reference checks: matches textbook Tomasi-Manduchi brute force to 3e-16 and scipy gaussian_filter Gaussian limit to 6e-13; skimage delta traces to its asymmetric spatial LUT window convention, not a divergence. Invariants (constant passthrough, rot90, convex hull) pass. LOW: docstring claims same dtype but float32 input gives float64 output (all backends). LOW: Inf-input divergence numpy(NaN) vs cupy(inf) at inf pixel, unvalidated-input pattern per previous audit."
 contour,2026-05-01,,,,"Marching squares correct: NaN check uses self-inequality, loop bounds (ny-1,nx-1) cover all quads, dask overlap depth=1 matches 2x2 stencil, float64 cast consistent across backends, saddle disambiguation via center value. No CRIT/HIGH issues; minor LOW (Inf inputs not specifically rejected) not flagged."
 corridor,2026-05-01,,LOW,1,"LOW: corridor inherits float32 from cost_distance; for very large accumulated costs, normalized = corridor - corridor_min loses precision near min (intrinsic to upstream dtype, not corridor itself). NaN handling correct (skipna min, np.isfinite check before normalize). All 4 backends route through pure xarray arithmetic; threshold uses dask/cupy/numpy where with try/except dispatch. No CRIT/HIGH issues."
 cost_distance,2026-06-16,3369,CRITICAL,5,"CRITICAL heap overflow (#3369/this PR): numba Dijkstra kernels _cost_distance_kernel + _cost_distance_tile_kernel sized the binary min-heap at height*width, but a lazy-deletion heap pushes a pixel on every improving relaxation, so push count exceeds h*w on non-uniform friction. _heap_push then writes OOB -> heap corruption, SIGABRT (exit 134, 'corrupted size vs prev_size') on iterative dask path; UB on numpy path. Reference heapq Dijkstra hits 44 pushes on a 6x6=36 grid. Fix: max_heap = h*w*(n_neighbors+1), tile kernel adds +2*(w+h)+4 for phase-2 boundary seeds. Verified: cupy relax kernel (parallel Bellman-Ford) does NOT use this heap, GPU path unaffected. CUDA available; numpy/cupy/dask+numpy/dask+cupy all agree post-fix over 30+40 random adversarial grids; 88 module tests pass (4 new regression tests). Cats 1-4 clean: dist float64 / out float32 fine; inf/nan/zero friction all impassable (tested); bounds guards use >=h/>=w; planar algorithm, no curvature expected. Supersedes prior #1191 (cupy max_iterations h+w->h*w, fixed PR #1192)."

xrspatial/bilateral.py

@@ -185,6 +185,15 @@ def _bilateral_cupy(data, radius, sigma_spatial, sigma_range):
     return out
 
 
+def _bilateral_cupy_boundary(data, radius, sigma_spatial, sigma_range,
+                             boundary='nan'):
+    if boundary == 'nan':
+        return _bilateral_cupy(data, radius, sigma_spatial, sigma_range)
+    padded = _pad_array(data, radius, boundary)
+    result = _bilateral_cupy(padded, radius, sigma_spatial, sigma_range)
+    return result[radius:-radius, radius:-radius]
+
+
 # ---------------------------------------------------------------------------
 # Dask + CuPy backend
 # ---------------------------------------------------------------------------
@@ -216,7 +225,10 @@ def _bilateral(data, radius, sigma_spatial, sigma_range, boundary='nan'):
             _bilateral_numpy_boundary,
             boundary=boundary,
         ),
-        cupy_func=_bilateral_cupy,
+        cupy_func=partial(
+            _bilateral_cupy_boundary,
+            boundary=boundary,
+        ),
         dask_func=partial(
             _bilateral_dask_numpy,
             boundary=boundary,

xrspatial/tests/test_bilateral.py

@@ -9,6 +9,7 @@
 import xarray as xr
 
 from xrspatial.bilateral import bilateral, _bilateral_numpy, _kernel_radius
+from xrspatial.utils import VALID_BOUNDARY_MODES
 from xrspatial.tests.general_checks import (
     assert_boundary_mode_correctness,
     create_test_raster,
@@ -334,6 +335,40 @@ def test_bilateral_boundary_modes(random_data_969):
     )
 
 
+@cuda_and_cupy_available
+@pytest.mark.parametrize('mode', VALID_BOUNDARY_MODES)
+def test_bilateral_boundary_modes_gpu(random_data_969, mode):
+    """cupy must honor every boundary mode, matching numpy (issue #3625)."""
+    numpy_agg = create_test_raster(random_data_969)
+    numpy_result = bilateral(numpy_agg, boundary=mode)
+
+    cupy_agg = create_test_raster(random_data_969, backend='cupy')
+    cupy_result = bilateral(cupy_agg, boundary=mode)
+    np.testing.assert_allclose(
+        numpy_result.data, cupy_result.data.get(),
+        equal_nan=True, rtol=1e-6,
+        err_msg=f'cupy diverges from numpy for boundary={mode!r}',
+    )
+
+
+@dask_array_available
+@cuda_and_cupy_available
+@pytest.mark.parametrize('mode', VALID_BOUNDARY_MODES)
+def test_bilateral_boundary_modes_dask_gpu(random_data_969, mode):
+    """dask+cupy must honor every boundary mode, matching numpy."""
+    numpy_agg = create_test_raster(random_data_969)
+    numpy_result = bilateral(numpy_agg, boundary=mode)
+
+    dask_cupy_agg = create_test_raster(random_data_969, backend='dask+cupy',
+                                       chunks=(15, 15))
+    dask_cupy_result = bilateral(dask_cupy_agg, boundary=mode)
+    np.testing.assert_allclose(
+        numpy_result.data, dask_cupy_result.data.compute().get(),
+        equal_nan=True, rtol=1e-4,
+        err_msg=f'dask+cupy diverges from numpy for boundary={mode!r}',
+    )
+
+
 # ---- 3D multi-band test ----
 
 def test_bilateral_3d():