Reason or Problem
_calc_hotspots_numpy in xrspatial/focal.py classifies each Gi* z-score with a nine-branch if/elif ladder. It first maps abs(zscore) to one of four p-values, then tests abs(zscore) again against three confidence thresholds combined with p_value < ..., and finally picks the sign in a third pair of branches. abs(zscore) gets evaluated up to six times per cell, and the data-dependent branches keep the loop from vectorizing.
The p-value step does nothing. Each p_value < threshold test is already implied by the abs(zscore) test sitting next to it:
|z| > 2.58 implies |z| >= 2.33, so p_value = 0.0099 < 0.01|z| > 1.96 implies |z| >= 1.65, so p_value <= 0.0495 < 0.05|z| > 1.65 implies |z| >= 1.65, so p_value <= 0.0495 < 0.1
The 1.29 threshold never reaches the output. The whole ladder reduces to: confidence 99 when |z| > 2.58, 95 when 1.96 < |z| <= 2.58, 90 when 1.65 < |z| <= 1.96, otherwise 0.
Proposal
Replace the ladder with boolean arithmetic:
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
NaN z-scores compare False everywhere and classify to 0, same as today. +inf and -inf classify to +99 and -99, same as today. Output stays int8 and the signature does not change.
A spike measured the rewrite at roughly 0.65x to 0.72x of the current wall time on a 2000x4000 float64 z-score array (random normal, smooth ramp, and 10% NaN inputs), with bitwise identical int8 output.
Design: Put a short comment above the new code stating the collapsed thresholds so nobody has to re-derive them from the old ladder. The CUDA device function _gpu_hotspots carries the same ladder. Leave it alone for now: branch cost on the GPU is a different question and that kernel is not on the timed CPU path.
Value: The classification runs once per cell after the convolution on the numpy and dask+numpy backends, so a cheaper classification shows up directly in hotspots() wall time on large rasters.
Stakeholders and Impacts
Users of hotspots() on numpy and dask+numpy. No API or output change.
Drawbacks
The arithmetic form is less literal than the ladder. That is the reason for the comment documenting the thresholds.
Alternatives
Keep the ladder but hoist abs(zscore) into a local. That removes the repeated calls but keeps the branches, so it recovers only part of the speedup.
Additional Notes or Context
Verification plan: extract the current function from origin/main into a scratch module and assert np.array_equal against the new one on a large random array, a dense linspace(-4, 4) sweep that crosses every threshold at both signs, the exact threshold values and their float64 neighbours, and an array with NaN, +inf, -inf and signed zeros. Add the boundary and non-finite cases as a pytest test. The asv class FocalHotspots already benchmarks this path.
Reason or Problem
_calc_hotspots_numpyinxrspatial/focal.pyclassifies each Gi* z-score with a nine-branchif/elifladder. It first mapsabs(zscore)to one of four p-values, then testsabs(zscore)again against three confidence thresholds combined withp_value < ..., and finally picks the sign in a third pair of branches.abs(zscore)gets evaluated up to six times per cell, and the data-dependent branches keep the loop from vectorizing.The p-value step does nothing. Each
p_value < thresholdtest is already implied by theabs(zscore)test sitting next to it:|z| > 2.58implies|z| >= 2.33, sop_value = 0.0099 < 0.01|z| > 1.96implies|z| >= 1.65, sop_value <= 0.0495 < 0.05|z| > 1.65implies|z| >= 1.65, sop_value <= 0.0495 < 0.1The 1.29 threshold never reaches the output. The whole ladder reduces to: confidence 99 when
|z| > 2.58, 95 when1.96 < |z| <= 2.58, 90 when1.65 < |z| <= 1.96, otherwise 0.Proposal
Replace the ladder with boolean arithmetic:
NaN z-scores compare False everywhere and classify to 0, same as today.
+infand-infclassify to+99and-99, same as today. Output staysint8and the signature does not change.A spike measured the rewrite at roughly 0.65x to 0.72x of the current wall time on a 2000x4000 float64 z-score array (random normal, smooth ramp, and 10% NaN inputs), with bitwise identical int8 output.
Design: Put a short comment above the new code stating the collapsed thresholds so nobody has to re-derive them from the old ladder. The CUDA device function
_gpu_hotspotscarries the same ladder. Leave it alone for now: branch cost on the GPU is a different question and that kernel is not on the timed CPU path.Value: The classification runs once per cell after the convolution on the numpy and dask+numpy backends, so a cheaper classification shows up directly in
hotspots()wall time on large rasters.Stakeholders and Impacts
Users of
hotspots()on numpy and dask+numpy. No API or output change.Drawbacks
The arithmetic form is less literal than the ladder. That is the reason for the comment documenting the thresholds.
Alternatives
Keep the ladder but hoist
abs(zscore)into a local. That removes the repeated calls but keeps the branches, so it recovers only part of the speedup.Additional Notes or Context
Verification plan: extract the current function from
origin/maininto a scratch module and assertnp.array_equalagainst the new one on a large random array, a denselinspace(-4, 4)sweep that crosses every threshold at both signs, the exact threshold values and their float64 neighbours, and an array with NaN, +inf, -inf and signed zeros. Add the boundary and non-finite cases as a pytest test. The asv classFocalHotspotsalready benchmarks this path.