From b44ac7d61bf9d56ee2882189c8092e890cfab28b Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 4 Sep 2026 11:48:11 -0400 Subject: [PATCH 1/2] Replace the hotspots confidence ladder with boolean arithmetic (#3737) _calc_hotspots_numpy classified each Gi* z-score through a nine-branch if/elif ladder that recomputed abs(zscore) up to six times per cell. The p-value step in that ladder never changed the output: every p_value < threshold test is implied by the |z| test beside it, so the whole thing reduces to 99 for |z| > 2.58, 95 for 1.96 < |z| <= 2.58, 90 for 1.65 < |z| <= 1.96, else 0. Compute the confidence as 90*(az > 1.65) + 5*(az > 1.96) + 4*(az > 2.58) and the sign as (z > 0) - (z < 0). On a 2000x4000 float64 z-score array the classifier drops from ~18.1 ms to ~11.6 ms (random normal), ~18.3 ms to ~12.6 ms (smooth ramp) and ~18.9 ms to ~12.1 ms (10% NaN), with identical int8 results. NaN still classifies to 0 and +/-inf to +/-99. Add tests pinning the classification at each threshold and its float64 neighbours at both signs, and for NaN, +/-inf and signed zeros. The CUDA device function keeps the ladder. --- xrspatial/focal.py | 33 +++++++++-------------------- xrspatial/tests/test_focal.py | 40 +++++++++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+), 23 deletions(-) diff --git a/xrspatial/focal.py b/xrspatial/focal.py index 3863bfbe2..0bd08faae 100644 --- a/xrspatial/focal.py +++ b/xrspatial/focal.py @@ -1469,29 +1469,16 @@ def _calc_hotspots_numpy(z_array): for x in prange(cols): zscore = z_array[y, x] - # find p value - p_value = 1.0 - if abs(zscore) >= 2.33: - p_value = 0.0099 - elif abs(zscore) >= 1.65: - p_value = 0.0495 - elif abs(zscore) >= 1.29: - p_value = 0.0985 - - # confidence - confidence = 0 - if abs(zscore) > 2.58 and p_value < 0.01: - confidence = 99 - elif abs(zscore) > 1.96 and p_value < 0.05: - confidence = 95 - elif abs(zscore) > 1.65 and p_value < 0.1: - confidence = 90 - - hot_cold = 0 - if zscore > 0: - hot_cold = 1 - elif zscore < 0: - hot_cold = -1 + # Confidence is 99 for |z| > 2.58, 95 for 1.96 < |z| <= 2.58, + # 90 for 1.65 < |z| <= 1.96, else 0. The GPU twin + # (_gpu_hotspots) spells this out as a p-value / confidence + # ladder; the p-value tests there are implied by the |z| + # tests, so the ladder collapses to these three thresholds. + # NaN compares False everywhere and classifies to 0. + az = abs(zscore) + confidence = (90 * (az > 1.65) + 5 * (az > 1.96) + + 4 * (az > 2.58)) + hot_cold = (zscore > 0) - (zscore < 0) out[y, x] = hot_cold * confidence return out diff --git a/xrspatial/tests/test_focal.py b/xrspatial/tests/test_focal.py index a8216675d..a4d51b7c2 100644 --- a/xrspatial/tests/test_focal.py +++ b/xrspatial/tests/test_focal.py @@ -1541,6 +1541,46 @@ def test_hotspots_dask_cupy(): dask_cupy_hotspots.data[pad:-pad, pad:-pad].compute().get()) +def test_hotspots_classifier_thresholds_3737(): + # Regression for #3737: the z-score classifier was rewritten from a + # nine-branch threshold ladder into boolean arithmetic. Pin the + # classification at every threshold and its float64 neighbours, at + # both signs, so the open/closed side of each interval cannot drift. + from xrspatial.focal import _calc_hotspots_numpy + + thresholds = np.array([1.29, 1.65, 1.96, 2.33, 2.58]) + below = np.nextafter(thresholds, -np.inf) + above = np.nextafter(thresholds, np.inf) + z = np.stack([below, thresholds, above]) + z = np.concatenate([z, -z]) + + # Confidence is 99 for |z| > 2.58, 95 for 1.96 < |z| <= 2.58, + # 90 for 1.65 < |z| <= 1.96, else 0. 1.29 and 2.33 are p-value + # thresholds in the original ladder and never change the output. + expected = np.array([ + [0, 0, 90, 95, 95], # just below each threshold + [0, 0, 90, 95, 95], # exactly on each threshold + [0, 90, 95, 95, 99], # just above each threshold + ], dtype=np.int8) + expected = np.concatenate([expected, -expected]) + + out = _calc_hotspots_numpy(z) + assert out.dtype == np.int8 + np.testing.assert_array_equal(out, expected) + + +def test_hotspots_classifier_nonfinite_3737(): + # NaN compares False against every threshold and classifies to 0; + # +/-inf land in the top band with the matching sign; signed zeros + # are neither hot nor cold. + from xrspatial.focal import _calc_hotspots_numpy + + z = np.array([[np.nan, np.inf, -np.inf, 0.0, -0.0]]) + out = _calc_hotspots_numpy(z) + assert out.dtype == np.int8 + np.testing.assert_array_equal(out, [[0, 99, -99, 0, 0]]) + + @dask_array_available @cuda_and_cupy_available def test_hotspots_dask_cupy_matches_numpy(): From b7d0446946915d9df627b1a6ebb46faf2fd16e2f Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Fri, 4 Sep 2026 11:51:09 -0400 Subject: [PATCH 2/2] Address review: pin float32 classification, tidy the comment (#3737) The production callers hand the classifier float32 z-scores, but the new threshold test only ran on float64. Add a test that the float32 result matches the classification of the same values widened to float64, on a grid that hits every band at both signs. Put the confidence expression on one line and drop the cross-reference to the GPU kernel from the comment so it does not go stale if that kernel is collapsed later. --- xrspatial/focal.py | 12 +++++------- xrspatial/tests/test_focal.py | 25 +++++++++++++++++++++++++ 2 files changed, 30 insertions(+), 7 deletions(-) diff --git a/xrspatial/focal.py b/xrspatial/focal.py index 0bd08faae..340b3936a 100644 --- a/xrspatial/focal.py +++ b/xrspatial/focal.py @@ -1470,14 +1470,12 @@ def _calc_hotspots_numpy(z_array): zscore = z_array[y, x] # Confidence is 99 for |z| > 2.58, 95 for 1.96 < |z| <= 2.58, - # 90 for 1.65 < |z| <= 1.96, else 0. The GPU twin - # (_gpu_hotspots) spells this out as a p-value / confidence - # ladder; the p-value tests there are implied by the |z| - # tests, so the ladder collapses to these three thresholds. - # NaN compares False everywhere and classifies to 0. + # 90 for 1.65 < |z| <= 1.96, else 0. This is the p-value / + # confidence ladder collapsed to the three thresholds that + # reach the output. NaN compares False everywhere and + # classifies to 0. az = abs(zscore) - confidence = (90 * (az > 1.65) + 5 * (az > 1.96) - + 4 * (az > 2.58)) + confidence = 90 * (az > 1.65) + 5 * (az > 1.96) + 4 * (az > 2.58) hot_cold = (zscore > 0) - (zscore < 0) out[y, x] = hot_cold * confidence diff --git a/xrspatial/tests/test_focal.py b/xrspatial/tests/test_focal.py index a4d51b7c2..3670387b6 100644 --- a/xrspatial/tests/test_focal.py +++ b/xrspatial/tests/test_focal.py @@ -1569,6 +1569,31 @@ def test_hotspots_classifier_thresholds_3737(): np.testing.assert_array_equal(out, expected) +def test_hotspots_classifier_float32_3737(): + # The production callers feed the classifier float32 z-scores + # (_hotspots_numpy casts before classifying). The thresholds are + # float64 literals, so a float32 input is promoted before comparison: + # np.float32(1.65) sits just below 1.65 and np.float32(1.96) just + # above it. Pin that the float32 result equals the classification of + # the same values widened to float64, so a dtype-specific comparison + # cannot creep in. + from xrspatial.focal import _calc_hotspots_numpy + + thresholds = np.array([1.29, 1.65, 1.96, 2.33, 2.58]) + z = np.stack([np.nextafter(thresholds, -np.inf), thresholds, + np.nextafter(thresholds, np.inf)]) + z = np.concatenate([z, -z]).astype(np.float32) + z = np.concatenate([np.nextafter(z, np.float32(-np.inf)), z, + np.nextafter(z, np.float32(np.inf))]) + + out32 = _calc_hotspots_numpy(z) + out64 = _calc_hotspots_numpy(z.astype(np.float64)) + assert out32.dtype == np.int8 + np.testing.assert_array_equal(out32, out64) + # The float32 grid must still hit every band on both sides. + assert set(np.unique(out32).tolist()) == {-99, -95, -90, 0, 90, 95, 99} + + def test_hotspots_classifier_nonfinite_3737(): # NaN compares False against every threshold and classifies to 0; # +/-inf land in the top band with the matching sign; signed zeros