Address CI AI review round 4: narrow Phase 1b guards, strict boundary

Reviewer correctly flagged that the 1%-of-median rule is a Phase 2
design="auto" heuristic, not Phase 1b. Backed off that over-reach.

P1 #1: Removed the min(d)/median(d) < 0.01 check. The mass-point
guard now applies uniformly (whenever d.min() > 0 and modal fraction
at d.min() > 2%) and does not gate on boundary. This still catches
the original concern (silently routing mass-point data through the
nonparametric branch) without rejecting valid Design 1' samples like
Beta(2,2) where d.min() is strictly positive but small.

P1 #2: Tightened boundary validation. The wrapper now accepts only
boundary ~ 0 (Design 1') or boundary ~ d.min() (Design 1 continuous-
near-d_lower) within float tolerance. Off-support values -- including
the previously-allowed "boundary < d.min()" path -- are rejected with
a targeted error message.

P3: Added a public-wrapper duplicate-support regression that drives a
rank-deficient X'X through the full selector stack (boundary =
d.min(), unique minimum, only 4 distinct d values) and asserts a
specific "qrXXinv" ValueError, not LinAlgError.

Test updates:
- Removed test_boundary_zero_with_positive_d_min_rejected: the case
  it modeled is now accepted (no mass point).
- Added test_boundary_zero_thin_boundary_density_accepted: Beta(2,2)
  Design 1' with vanishing boundary density now passes.
- Added test_off_support_boundary_rejected: boundary=0.5 on U(1,2).
- Added test_negative_boundary_rejected: boundary<0 rejected.
- Updated test_nonzero_boundary: uses boundary=float(d.min()), not
  boundary=1.0 (which is off the realized support of U(1,2)).

175 tests pass (up from 172).

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

Author: igerber

diff_diff/local_linear.py

@@ -606,100 +606,57 @@ def mse_optimal_bandwidth(
         raise ValueError(f"boundary must be finite; got {boundary}")
 
     # Boundary-applicability check (Phase 1b scope).
-    # The exported wrapper is scoped to the HAD LOWER-boundary case:
-    #   - Design 1' (d_0 = 0 with support infimum at 0), or
-    #   - Design 1 continuous-near-d_lower (d_0 = d_lower = min(D_2)).
-    # In both cases ``boundary`` must lie at or below the sample minimum:
-    # no observation should fall to the left of the evaluation point,
-    # because the downstream ``local_linear_fit`` uses a one-sided kernel
-    # restricted to ``d >= d_0``. An interior or upper-boundary
-    # evaluation would silently run the boundary selector branch with a
-    # symmetric kernel (see port's ``kernel_W``) and produce a bandwidth
-    # incompatible with the fitter. Reject those inputs here until an
-    # interior-point API is documented (Phase 2+).
+    # The exported wrapper is scoped to the two documented HAD
+    # nonparametric estimands:
+    #   - Design 1' evaluates m(0) = lim_{d down 0} E[Delta Y | D_2 <= d]
+    #     with ``boundary = 0``.
+    #   - Design 1 continuous-near-d_lower evaluates m(d_lower) with
+    #     ``boundary = d_lower = min(D_2)``.
+    # Any other off-support boundary (interior, upper-boundary, or an
+    # arbitrary value between 0 and d.min()) would silently target an
+    # undocumented limit and is rejected.
     d_min = float(d.min())
-    # Allow a small numerical tolerance so boundary = min(d) passes even
-    # under floating-point noise (e.g. ``d.min()`` equals ``boundary``
-    # within 1e-12 but not bit-exact).
     _boundary_tol = 1e-12 * max(1.0, abs(d_min), abs(boundary))
-    if boundary > d_min + _boundary_tol:
+    _at_zero = abs(boundary) <= _boundary_tol
+    _at_d_min = abs(boundary - d_min) <= _boundary_tol
+    if not (_at_zero or _at_d_min):
         raise ValueError(
-            f"boundary={boundary!r} exceeds the sample minimum "
-            f"d.min()={d_min!r}. The Phase 1b wrapper is scoped to the "
-            f"HAD lower-boundary case (d_0 <= min(d)). For interior or "
-            f"upper-boundary evaluation, use "
-            f"diff_diff._nprobust_port.lpbwselect_mse_dpi directly with "
-            f"interior=True; note that the non-boundary branches are not "
-            f"separately parity-tested against nprobust."
+            f"boundary={boundary!r} is not at a supported HAD estimand. "
+            f"The Phase 1b public wrapper accepts only boundary ~ 0 "
+            f"(Design 1') or boundary ~ d.min()={d_min!r} (Design 1 "
+            f"continuous-near-d_lower). Off-support values would "
+            f"silently target an undocumented limit. For interior or "
+            f"other boundary points, use "
+            f"diff_diff._nprobust_port.lpbwselect_mse_dpi directly and "
+            f"note that those paths are not separately parity-tested."
         )
 
-    # Design classification (REGISTRY.md Phase 2 design="auto" rule
-    # plus paper Section 3.2.4 mass-point redirection):
-    #
-    #   Design 1'        <=> d_lower = inf supp(D_2) = 0
-    #   Design 1 masspt  <=> d_lower > 0, P(D_2 = d_lower) > 0
-    #   Design 1 cont    <=> d_lower > 0, D_2 continuous near d_lower
-    #
-    # The CCF nonparametric selector is appropriate for Design 1' and
-    # Design 1 continuous-near-d_lower. Mass-point Design 1 requires a
-    # 2SLS sample-average estimator (Phase 2, not Phase 1b).
+    # Mass-point design check (paper Section 3.2.4, REGISTRY 2% rule).
+    # When d_min > 0 and there is bunching at d_min, Design 1 requires
+    # the 2SLS sample-average path (Phase 2), not the CCF nonparametric
+    # selector. The check applies independently of the boundary the
+    # user supplied: mass-point data is never appropriate for this
+    # wrapper. The check explicitly excludes d_min ~ 0, which is the
+    # Design 1' "untreated units present" subcase that the paper's
+    # simulations and the Garrett et al. (2020) application accept.
     _MASS_POINT_THRESHOLD = 0.02  # REGISTRY rule: > 2% modal-min
-    eps_eq = 1e-12 * max(1.0, abs(d_min))
-    at_d_min_mask = np.abs(d - d_min) <= eps_eq
-    modal_fraction = float(np.mean(at_d_min_mask))
-
-    if boundary > _boundary_tol:
-        # User supplied boundary = d_min > 0 explicitly. This is
-        # Design 1; distinguish mass-point vs continuous-near-d_lower.
+    if d_min > _boundary_tol:
+        eps_eq = 1e-12 * max(1.0, abs(d_min))
+        at_d_min_mask = np.abs(d - d_min) <= eps_eq
+        modal_fraction = float(np.mean(at_d_min_mask))
         if modal_fraction > _MASS_POINT_THRESHOLD:
             raise NotImplementedError(
-                f"Detected mass-point design: the lower boundary "
-                f"d_lower={d_min!r} has modal fraction "
-                f"{modal_fraction:.4f} > {_MASS_POINT_THRESHOLD:.2f}. "
-                f"Per de Chaisemartin et al. (2026) Section 3.2.4 and "
-                f"the methodology registry, this case requires the 2SLS "
-                f"sample-average estimator with instrument 1{{D_2 > "
-                f"d_lower}}, not the nonparametric CCF local-polynomial "
-                f"bandwidth selector. That estimator is queued for "
+                f"Detected mass-point design at d.min()={d_min!r} "
+                f"(modal fraction {modal_fraction:.4f} > "
+                f"{_MASS_POINT_THRESHOLD:.2f}). Per de Chaisemartin et "
+                f"al. (2026) Section 3.2.4, Design 1 mass-point cases "
+                f"require the 2SLS sample-average estimator with "
+                f"instrument 1{{D_2 > d_lower}}, not the CCF "
+                f"nonparametric selector. That path is queued for "
                 f"Phase 2 (HeterogeneousAdoptionDiD). For continuous "
                 f"near-d_lower designs (modal fraction <= "
-                f"{_MASS_POINT_THRESHOLD:.2f}), this wrapper is applicable."
-            )
-    else:
-        # User passed boundary <= 0 (default). Intent is Design 1'.
-        # Apply the REGISTRY's `min(d) < 0.01 * median(d)` rule: only
-        # accept when the support infimum is effectively 0 relative to
-        # the data scale. Otherwise force the user to disambiguate.
-        d_median = float(np.median(d))
-        _DESIGN_1_PRIME_RATIO = 0.01  # REGISTRY: min(d)/median(d) < 1%
-        # Guard against pathological all-zero or all-negative median.
-        effective_threshold = _DESIGN_1_PRIME_RATIO * max(abs(d_median), 1e-12)
-        if d_min > effective_threshold:
-            if modal_fraction > _MASS_POINT_THRESHOLD:
-                # Mass-point design dressed as Design 1': same
-                # NotImplementedError as the boundary>0 branch.
-                raise NotImplementedError(
-                    f"Detected mass-point design at d.min()={d_min!r} "
-                    f"(modal fraction {modal_fraction:.4f} > "
-                    f"{_MASS_POINT_THRESHOLD:.2f}), but boundary=0 "
-                    f"implies Design 1' intent. Either: (a) pass "
-                    f"boundary=d.min() to make the mass-point "
-                    f"classification explicit (which will then raise "
-                    f"NotImplementedError directing to the 2SLS path "
-                    f"per de Chaisemartin et al. 2026 Section 3.2.4), "
-                    f"or (b) reconsider whether the data is truly "
-                    f"Design 1' (support at 0)."
-                )
-            raise ValueError(
-                f"Ambiguous design: boundary=0 but d.min()={d_min!r} "
-                f"exceeds the Design 1' threshold of "
-                f"0.01 * median(d) = {effective_threshold!r}. This "
-                f"dataset does not satisfy Design 1' (support infimum "
-                f"at 0). Either: (a) pass boundary=d.min() for the "
-                f"Design 1 continuous-near-d_lower path, or (b) verify "
-                f"the data truly has support at 0 (in which case "
-                f"d.min() would be much closer to zero relative to the "
-                f"data scale)."
+                f"{_MASS_POINT_THRESHOLD:.2f}), this wrapper is "
+                f"applicable."
             )
 
     # Defer heavy import to call time to avoid import-cycle risk.

tests/test_bandwidth_selector.py

@@ -313,30 +313,40 @@ def test_boundary_equal_to_min_d_accepted(self):
         assert h > 0.0
 
     def test_boundary_zero_design_1_prime_accepted(self):
-        """Design 1' data: d.min() effectively 0 relative to median,
-        so boundary=0 passes the REGISTRY 1%-of-median rule."""
+        """Design 1' with support at 0: boundary=0 passes."""
         rng = np.random.default_rng(20260419)
-        # d.min() ~ 5e-4 << 0.01 * median(d) ~ 5e-3
         d = rng.uniform(0.0, 1.0, size=3000)
         y = d + d**2 + rng.normal(0, 0.5, size=3000)
         h = mse_optimal_bandwidth(d, y, boundary=0.0)
         assert np.isfinite(h)
         assert h > 0.0
 
-    def test_boundary_zero_with_positive_d_min_rejected(self):
-        """Default boundary=0 with d.min() substantially positive (>1%
-        of median) is not Design 1'; force the caller to pass
-        boundary=d.min() instead of silently misclassifying."""
+    def test_boundary_zero_thin_boundary_density_accepted(self):
+        """Beta(2,2) Design 1' case: boundary density vanishes at 0
+        but the estimand is well-defined. Must not be mistakenly
+        rejected by any design heuristic."""
+        rng = np.random.default_rng(20260419)
+        d = rng.beta(2.0, 2.0, size=2000)  # f_D(0) = 0
+        y = d + d**2 + rng.normal(0, 0.5, size=2000)
+        h = mse_optimal_bandwidth(d, y, boundary=0.0)
+        assert np.isfinite(h)
+        assert h > 0.0
+
+    def test_boundary_zero_with_data_far_from_zero_fails_gracefully(self):
+        """boundary=0 passes the boundary-validation and mass-point
+        checks but then hits the per-stage count guard deeper in the
+        selector because the kernel window is empty (d ~ U(0.5, 1.0)
+        has no data near 0). Must surface a clear ValueError, not an
+        opaque failure."""
         rng = np.random.default_rng(2026)
-        # d.min() ~ 0.5, median ~ 0.75 -> ratio 0.67 >> 0.01.
-        d = rng.uniform(0.5, 1.0, size=1500)
+        d = rng.uniform(0.5, 1.0, size=1500)  # d.min() ~ 0.5, no mass
         y = d + rng.normal(0, 0.3, size=1500)
-        with pytest.raises(ValueError, match="Ambiguous design"):
+        with pytest.raises(ValueError, match="lprobust_bw"):
             mse_optimal_bandwidth(d, y, boundary=0.0)
 
     def test_boundary_zero_with_d_min_mass_point_rejected(self):
-        """boundary=0 default path with d.min() > 0 AND mass at d.min()
-        must also be rejected, pointing to the 2SLS redirection."""
+        """boundary=0 with d.min() > 0 AND mass at d.min() is a
+        Design 1 mass-point design and must be redirected to 2SLS."""
         rng = np.random.default_rng(2026)
         n_mass = 300  # 15% at 0.1
         n_cont = 1700
@@ -347,6 +357,25 @@ def test_boundary_zero_with_d_min_mass_point_rejected(self):
         with pytest.raises(NotImplementedError, match="mass-point"):
             mse_optimal_bandwidth(d, y, boundary=0.0)
 
+    def test_off_support_boundary_rejected(self):
+        """boundary must equal 0 or d.min() within tolerance; any
+        other lower off-support value must be rejected."""
+        rng = np.random.default_rng(2026)
+        d = rng.uniform(1.0, 2.0, size=1500)  # d.min() ~ 1.0
+        y = d + rng.normal(0, 0.3, size=1500)
+        # boundary = 0.5 is between 0 and d.min(); neither documented
+        # estimand.
+        with pytest.raises(ValueError, match="not at a supported HAD estimand"):
+            mse_optimal_bandwidth(d, y, boundary=0.5)
+
+    def test_negative_boundary_rejected(self):
+        """boundary < 0 is off-support and rejected."""
+        rng = np.random.default_rng(2026)
+        d = rng.uniform(0.0, 1.0, size=1500)
+        y = d + rng.normal(0, 0.3, size=1500)
+        with pytest.raises(ValueError, match="not at a supported HAD estimand"):
+            mse_optimal_bandwidth(d, y, boundary=-0.1)
+
     def test_mass_point_design_rejected(self):
         """Design 1 mass-point case (boundary > 0, modal fraction > 2%)
         must be rejected with NotImplementedError pointing to 2SLS."""
@@ -394,30 +423,33 @@ def test_rank_deficient_design_raises_valueerror(self):
         with pytest.raises(ValueError, match="qrXXinv"):
             qrXXinv(X)
 
-    def test_full_stack_rank_deficient_raises_valueerror(self):
-        """End-to-end rank-deficiency integration test: a sample with
-        enough rows but only a handful of distinct in-window d values
-        drives rank-deficient X'X at some DPI stage. The public wrapper
-        must surface a clear ValueError (either from qrXXinv's Cholesky
-        guard or from the stage-count guards), not IndexError or
-        LinAlgError.
+    def test_wrapper_rank_deficient_raises_valueerror(self):
+        """Public-wrapper regression: a continuous-near-d_lower sample
+        whose kernel window contains too few DISTINCT d values drives
+        a rank-deficient X'X in one of the DPI stages. The wrapper
+        must surface a clear ValueError from qrXXinv's Cholesky guard,
+        not an opaque LinAlgError.
+
+        Construction: d.min() is unique (modal_fraction = 1/G < 2% so
+        mass-point check passes), but the remaining data concentrates
+        on a single value so the kernel window has only 2 distinct d
+        values and the design-matrix columns become linearly
+        dependent at higher polynomial orders.
         """
-        from diff_diff._nprobust_port import lpbwselect_mse_dpi
-
-        # Only 3 distinct d values; each pilot stage's design matrix
-        # will have at most 3 independent rows, but the B2 fit needs
-        # o_B+2 = 5 or 6 columns -> rank deficient.
-        d = np.repeat([0.1, 0.2, 0.3], 50).astype(np.float64)
         rng = np.random.default_rng(2026)
-        y = d + rng.normal(0, 0.1, size=d.size)
-        with pytest.raises(ValueError):
-            lpbwselect_mse_dpi(y, d, eval_point=0.0, bwcheck=None)
-        # Same through the public wrapper (boundary=0 is fine since
-        # d.min()=0.1 triggers either ambiguous-design or rank-deficient
-        # rejection; we only care that it is a ValueError and not an
-        # opaque linear-algebra crash).
-        with pytest.raises((ValueError, NotImplementedError)):
-            mse_optimal_bandwidth(d, y)
+        # G = 151: d.min=0.1 unique, 50 obs each at 0.15 / 0.3 / 0.4.
+        # Modal fraction = 1/151 < 2% passes mass-point check.
+        # The B1 / B2 auxiliary fits at stage d1 use h_B1 = h_B2 =
+        # range = 0.3, which retains all 4 distinct values (0.1, 0.15,
+        # 0.3, 0.4). The B1 design matrix has o_B+1 = 5 columns but
+        # only 4 independent rows -> rank-deficient X'X -> Cholesky
+        # fails in qrXXinv.
+        d = np.concatenate(
+            [[0.1], np.full(50, 0.15), np.full(50, 0.3), np.full(50, 0.4)]
+        )
+        y = d + rng.normal(0, 0.01, size=d.size)
+        with pytest.raises(ValueError, match="qrXXinv"):
+            mse_optimal_bandwidth(d, y, boundary=float(d.min()))
 
 
 class TestKernelDispatch:
@@ -446,12 +478,15 @@ def test_kernel_epa_vs_tri_differ(self):
 
 class TestBoundary:
     def test_nonzero_boundary(self):
-        """Design 1 continuous-near-d_lower case: evaluate at a non-zero
-        boundary d_0 = d_lower."""
+        """Design 1 continuous-near-d_lower case: boundary = d.min()
+        (not the theoretical infimum of the support). Under the
+        strict boundary-applicability check, the user must pass the
+        sample minimum, not a known theoretical lower bound like 1.0
+        on U(1, 2) data."""
         rng = np.random.default_rng(2026)
         d = rng.uniform(1.0, 2.0, size=1500)
         y = 3.0 + 0.5 * (d - 1.0) ** 2 + rng.normal(0, 0.3, size=1500)
-        h = mse_optimal_bandwidth(d, y, boundary=1.0)
+        h = mse_optimal_bandwidth(d, y, boundary=float(d.min()))
         assert np.isfinite(h)
         assert h > 0.0