fix(rdd): exact-constancy variance rescue for constant-float outcome sides

np.var of a constant 0.7 vector is ~1.3e-32 (single-pass mean roundoff)
while R's two-pass var() is exactly 0, so the var==0 rescue was skipped
for constant non-representable-float outcomes and the resulting Inf J_MV
crashed bin construction on a zero jump. Route the side variances through
the port's _var0 exact-constancy helper (the R semantics already ported
for the fuzzy identification stop), reproducing R's J=1 rescue, its Inf
J_IMSE / NaN J_MV echoes, and the variability warning. Regression test
with a constant-0.7 side added.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VwhFMnFQGBumYfwbeQaUjm

Author: igerber

diff_diff/rdplot.py

@@ -39,6 +39,7 @@
 from ._rdrobust_port import (
     _covs_gamma,
     _normalize_kernel,
+    _var0,
     covs_drop_fun,
     qrXXinv,
     rdrobust_kweight,
@@ -854,8 +855,15 @@ def _selector_grams() -> Tuple[int, np.ndarray, np.ndarray, np.ndarray, np.ndarr
         drk_i_l = _pow_outer(x_bar_i_l, k - 1) * j_arange
         drk_i_r = _pow_outer(x_bar_i_r, k - 1) * j_arange
 
-        var_y_l = float(np.var(y_l, ddof=1))
-        var_y_r = float(np.var(y_r, ddof=1))
+        # R's two-pass var() is EXACTLY zero on a constant vector; numpy's
+        # single-pass mean leaves ~1e-32 roundoff for constants like 0.7,
+        # which would skip the var==0 rescue below and crash the bin grid on
+        # a zero jump (J_MV = Inf). Route through the port's exact-constancy
+        # helper (_var0, the #686-ported R semantics) so constant sides give
+        # var_y = 0.0 exactly, reproducing R's rescue and its 0/0 -> NaN
+        # selector echoes.
+        var_y_l = 0.0 if _var0(y_l) else float(np.var(y_l, ddof=1))
+        var_y_r = 0.0 if _var0(y_r) else float(np.var(y_r, ddof=1))
 
         mu0_i_l = rk_i_l @ gamma_k1_l
         mu0_i_r = rk_i_r @ gamma_k1_r

tests/test_rdplot.py

@@ -350,6 +350,27 @@ def test_manual_h_with_zero_effective_obs_raises(self):
         with pytest.raises(ValueError, match="zero effective observations"):
             RDPlot(h=0.1).fit(df, "y", "x")
 
+    def test_constant_float_outcome_side_rescued_exactly(self):
+        # np.var of a constant 0.7 vector is ~1.3e-32, NOT 0 (single-pass
+        # mean roundoff), while R's two-pass var() is exactly 0 - without
+        # the _var0 exact-constancy routing the var==0 rescue is skipped
+        # and the zero jump crashes bin construction (ZeroDivisionError;
+        # caught by CI codex on the rebased head). R fires the rescue:
+        # J=1 with the variability warning, Inf J_IMSE, NaN J_MV.
+        rng = np.random.default_rng(7)
+        x = np.concatenate([-rng.uniform(0.1, 1, 15), rng.uniform(0.1, 1, 15)])
+        y = np.where(x < 0, 0.7, rng.normal(size=30))
+        df = pd.DataFrame({"y": y, "x": x})
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            r = RDPlot().fit(df, "y", "x")
+        messages = [str(w.message) for w in caught]
+        assert any("variability" in m and "below" in m for m in messages)
+        assert r.J[0] == 1.0
+        assert np.isinf(r.J_IMSE[0])
+        assert np.isnan(r.J_MV[0])  # 0/0 through the MV selector, as in R
+        assert np.isfinite(r.J[1]) and r.J[1] >= 1
+
     def test_finite_singular_selector_gram_does_not_downgrade(self):
         # Companion to the overflow-ladder golden: a finite rank-deficient
         # k=4 Gram (4 distinct x values per side) is absorbed inside