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.

Author: brendancol

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

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():