Address CI review P0 + P1 on Phase 1c

P0: bias_corrected_local_linear now routes the CI through
`safe_inference()` so degenerate cases with `se_robust <= 0` or
non-finite `se_robust` (e.g., exact-fit / constant-y) return
`(NaN, NaN)` rather than a misleading zero-width or infinite CI.
Matches the repo-wide inference contract (CLAUDE.md Key Design
Pattern #6).

P1: Auto-bandwidth path now calls `lpbwselect_mse_dpi` directly
with `cluster`, `vce`, and `nnmatch` forwarded. Previously it went
through `mse_optimal_bandwidth` which hard-codes unclustered / nn /
nnmatch=3, silently mismatching the downstream `lprobust` fit's
reported estimator.

Tests added: TestNaNSafeCI (constant-y + near-zero-SE) and
TestAutoBandwidthForwardsParameters (cluster+auto, vce='hc1'+auto,
nnmatch=5+auto), all asserting the selected bandwidth changes when
the corresponding parameter changes (catches silent fallback).

Also: suppress spurious BLAS FPE warnings in lprobust_bw's hc/hc2/hc3
branch (numpy issue #21432 pattern), newly reachable via the
wired-through vce='hc1' auto-bandwidth path.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: igerber

diff_diff/_nprobust_port.py

@@ -519,12 +519,18 @@ def lprobust_bw(
             dups_B = dups[ind1] if dups is not None else None
             dupsid_B = dupsid[ind1] if dupsid is not None else None
         if vce in ("hc0", "hc1", "hc2", "hc3"):
-            predicts_B = R_B1 @ beta_B1
-            if vce in ("hc2", "hc3"):
-                hii_B = np.empty(n_B1, dtype=np.float64)
-                RW1 = R_B1 * eW1[:, None]
-                for i in range(n_B1):
-                    hii_B[i] = R_B1[i, :] @ invG_B1 @ RW1[i, :]
+            # Suppress spurious BLAS FPE warnings (numpy issue #21432
+            # pattern); matmul on some platforms (Accelerate / OpenBLAS)
+            # sets divide/overflow flags on SIMD intermediates even when
+            # input and output are finite.
+            with np.errstate(divide="ignore", over="ignore",
+                             invalid="ignore", under="ignore"):
+                predicts_B = R_B1 @ beta_B1
+                if vce in ("hc2", "hc3"):
+                    hii_B = np.empty(n_B1, dtype=np.float64)
+                    RW1 = R_B1 * eW1[:, None]
+                    for i in range(n_B1):
+                        hii_B[i] = R_B1[i, :] @ invG_B1 @ RW1[i, :]
         res_B = lprobust_res(
             eX1,
             eY1,

diff_diff/local_linear.py

@@ -1073,16 +1073,59 @@ def bias_corrected_local_linear(
     # bandwidth_diagnostics.b_mse for callers that want to inspect or
     # override. The paper (de Chaisemartin et al. 2026) likewise uses a
     # single h*_G throughout Equation 8.
+    #
+    # In auto mode, cluster / vce / nnmatch are forwarded to
+    # ``lpbwselect_mse_dpi`` so bandwidth selection reflects the same
+    # estimator the final ``lprobust`` call will use. Calling
+    # ``mse_optimal_bandwidth`` (the public wrapper) instead would hard-code
+    # ``cluster=None, vce="nn", nnmatch=3`` and silently mismatch the
+    # downstream fit — a methodology bug (CI review PR #340 P1).
     bw_source: str
     bw_diag: Optional[BandwidthResult] = None
     if h is None and b is None:
-        bw_diag = mse_optimal_bandwidth(
-            d=d,
+        # Defer heavy import to call time to avoid import-cycle risk.
+        from diff_diff._nprobust_port import lpbwselect_mse_dpi
+
+        stages = lpbwselect_mse_dpi(
             y=y,
-            boundary=boundary,
+            x=d,
+            cluster=cluster_arr,
+            eval_point=float(boundary),
+            p=1,
+            q=2,
+            deriv=0,
+            kernel=nprobust_kernel,
+            bwcheck=21,
+            bwregul=1.0,
+            vce=vce,
+            nnmatch=nnmatch,
+            interior=False,
+        )
+        bw_diag = BandwidthResult(
+            h_mse=stages.h_mse_dpi,
+            b_mse=stages.b_mse_dpi,
+            c_bw=stages.c_bw,
+            bw_mp2=stages.bw_mp2,
+            bw_mp3=stages.bw_mp3,
+            stage_d1_V=stages.stage_d1.V,
+            stage_d1_B1=stages.stage_d1.B1,
+            stage_d1_B2=stages.stage_d1.B2,
+            stage_d1_R=stages.stage_d1.R,
+            stage_d2_V=stages.stage_d2.V,
+            stage_d2_B1=stages.stage_d2.B1,
+            stage_d2_B2=stages.stage_d2.B2,
+            stage_d2_R=stages.stage_d2.R,
+            stage_b_V=stages.stage_b.V,
+            stage_b_B1=stages.stage_b.B1,
+            stage_b_B2=stages.stage_b.B2,
+            stage_b_R=stages.stage_b.R,
+            stage_h_V=stages.stage_h.V,
+            stage_h_B1=stages.stage_h.B1,
+            stage_h_B2=stages.stage_h.B2,
+            stage_h_R=stages.stage_h.R,
+            n=int(d.shape[0]),
             kernel=kernel,
-            weights=None,  # already validated
-            return_diagnostics=True,
+            boundary=float(boundary),
         )
         h_val = float(bw_diag.h_mse)
         b_val = h_val  # rho=1 default to match nprobust
@@ -1125,14 +1168,21 @@ def bias_corrected_local_linear(
         bwcheck=21,
     )
 
-    # --- Bias-corrected CI (lprobust summary.lprobust:420-421) ---
-    # z_{1 - alpha/2}; Python uses scipy.stats.norm.ppf. For parity with R
-    # on the golden tests, the golden JSON stores R's qnorm value.
-    from scipy.stats import norm as _norm
-
-    z = float(_norm.ppf(1.0 - alpha / 2.0))
-    ci_low = result.tau_bc - z * result.se_rb
-    ci_high = result.tau_bc + z * result.se_rb
+    # --- Bias-corrected CI via safe_inference (NaN-safe gate) ---
+    # When se_robust is zero, negative, or non-finite (e.g., exact-fit
+    # cases where the residual vector collapses), ALL inference fields —
+    # including the CI — must return NaN. This enforces the repo-wide
+    # inference contract (CLAUDE.md Key Design Pattern #6; CI review
+    # PR #340 P0) rather than returning a misleading zero-width or infinite
+    # CI. safe_inference also handles the R z = qnorm(1 - alpha/2) critical
+    # value via scipy.stats.norm.ppf (the golden JSON stores R's z so
+    # parity tests consume R's value directly and drift is pure
+    # tau.bc + z * se.rb arithmetic).
+    from diff_diff.utils import safe_inference
+
+    _, _, (ci_low, ci_high) = safe_inference(
+        result.tau_bc, result.se_rb, alpha=float(alpha)
+    )
 
     return BiasCorrectedFit(
         estimate_classical=result.tau_cl,

tests/test_bias_corrected_lprobust.py

@@ -387,6 +387,114 @@ def test_auto_bandwidth_covers_truth_at_95(self):
         assert fit.ci_low <= beta0_true <= fit.ci_high
 
 
+# =============================================================================
+# NaN-safe CI (CI review PR #340 P0)
+# =============================================================================
+
+
+class TestNaNSafeCI:
+    """``bias_corrected_local_linear`` must route the CI through
+    ``safe_inference`` so degenerate cases with ``se_robust <= 0`` or
+    non-finite ``se_robust`` surface as ``ci_low = ci_high = NaN`` rather
+    than a misleading zero-width or infinite CI."""
+
+    def test_constant_y_returns_nan_ci(self):
+        """Constant y makes residuals zero; se_robust collapses to 0. CI
+        must be (NaN, NaN), not a finite zero-width CI."""
+        d = np.linspace(0.0, 1.0, 200)
+        y = np.full_like(d, 1.5)  # zero residuals everywhere
+        fit = bias_corrected_local_linear(d, y, h=0.3, b=0.3)
+        assert fit.se_robust == 0.0 or not np.isfinite(fit.se_robust)
+        assert np.isnan(fit.ci_low)
+        assert np.isnan(fit.ci_high)
+
+    def test_near_zero_se_returns_nan_ci(self):
+        """Near-constant y produces a tiny se_robust; the NaN-safe gate
+        fires when it hits zero exactly (covers the exact-fit edge case
+        the CI review flagged without tripping a pre-existing Phase 1b
+        ZeroDivisionError in the auto-bandwidth selector on truly
+        constant y, which is tracked separately)."""
+        rng = np.random.default_rng(0)
+        d = np.linspace(0.0, 1.0, 500)
+        # Residuals near machine epsilon; tau_bc stays finite.
+        y = 0.1 * np.ones_like(d) + rng.normal(0, 1e-300, size=d.shape)
+        fit = bias_corrected_local_linear(d, y, h=0.3, b=0.3)
+        # Either the inference should be fully valid, OR the CI gate has
+        # correctly fired. The contract is: se_rb <= 0 / non-finite =>
+        # NaN CI.
+        if not (np.isfinite(fit.se_robust) and fit.se_robust > 0):
+            assert np.isnan(fit.ci_low)
+            assert np.isnan(fit.ci_high)
+
+
+# =============================================================================
+# Auto-bandwidth parameter forwarding (CI review PR #340 P1)
+# =============================================================================
+
+
+class TestAutoBandwidthForwardsParameters:
+    """Auto-bandwidth must forward ``cluster``, ``vce``, and ``nnmatch`` to
+    the bandwidth selector. Calling the public ``mse_optimal_bandwidth``
+    wrapper would hard-code ``cluster=None, vce="nn", nnmatch=3`` and
+    silently mismatch the downstream ``lprobust`` fit — a methodology
+    bug. These tests pin the correct wiring."""
+
+    def _smoke_data(self, seed=33):
+        rng = np.random.default_rng(seed)
+        G = 600
+        d = rng.uniform(0.0, 1.0, G)
+        y = d + d ** 2 + rng.normal(0, 0.3, G)
+        return d, y
+
+    def test_auto_cluster_returns_finite(self):
+        """Auto-bandwidth with cluster produces a finite BiasCorrectedFit.
+
+        No R parity anchor (nprobust's internal lpbwselect has a
+        singleton-cluster bug on the pilot fits); this test pins that the
+        Python path completes and uses the clustered bandwidth downstream,
+        not the unclustered one.
+        """
+        d, y = self._smoke_data()
+        cluster = np.repeat(np.arange(30), 20)
+        fit_cluster = bias_corrected_local_linear(d, y, cluster=cluster)
+        fit_uncluster = bias_corrected_local_linear(d, y)
+        assert fit_cluster.bandwidth_source == "auto"
+        assert np.isfinite(fit_cluster.estimate_bias_corrected)
+        assert np.isfinite(fit_cluster.se_robust)
+        # The clustered bandwidth should differ from the unclustered one
+        # (different residual meat feeds into Stage-2/3 AMSE minimization).
+        # If the wrapper were silently passing cluster=None, these would
+        # be identical to bit-parity.
+        assert fit_cluster.h != fit_uncluster.h
+
+    def test_auto_vce_hc1_returns_finite(self):
+        """Auto-bandwidth with non-default vce must use the requested vce
+        during bandwidth selection, not silently fall back to nn."""
+        d, y = self._smoke_data()
+        fit_hc1 = bias_corrected_local_linear(d, y, vce="hc1")
+        fit_nn = bias_corrected_local_linear(d, y, vce="nn")
+        assert fit_hc1.bandwidth_source == "auto"
+        assert np.isfinite(fit_hc1.estimate_bias_corrected)
+        assert np.isfinite(fit_hc1.se_robust)
+        # Different residual definitions yield different stage-2/3 AMSE
+        # and therefore different bandwidths. Bit-identity would indicate
+        # the selector silently ignored vce.
+        assert fit_hc1.h != fit_nn.h
+
+    def test_auto_nnmatch_non_default_returns_finite(self):
+        """Auto-bandwidth with non-default nnmatch must forward it to the
+        selector, not silently use the hard-coded default of 3."""
+        d, y = self._smoke_data()
+        fit_nn5 = bias_corrected_local_linear(d, y, nnmatch=5)
+        fit_nn3 = bias_corrected_local_linear(d, y, nnmatch=3)
+        assert fit_nn5.bandwidth_source == "auto"
+        assert np.isfinite(fit_nn5.estimate_bias_corrected)
+        # nnmatch controls the NN residual construction; different values
+        # give different meat matrices and therefore different stage
+        # bandwidths.
+        assert fit_nn5.h != fit_nn3.h
+
+
 # =============================================================================
 # Validator idempotence (regression gate for the Phase 1b extraction)
 # =============================================================================