"""
R DIDmultiplegtDYN parity tests for the dCDH estimator.

Loads pre-computed golden values from
``benchmarks/data/dcdh_dynr_golden_values.json`` (generated by the R
script at ``benchmarks/R/generate_dcdh_dynr_test_values.R``) and asserts
that the Python ``ChaisemartinDHaultfoeuille`` implementation matches R
``did_multiplegt_dyn`` at l = 1 (which is numerically identical to the
AER 2020 paper's ``DID_M``).

**R is only needed to regenerate the JSON file**, not to run these tests.
The committed JSON is loaded directly and the assertions run on any
Python-only environment. Tests skip ONLY if the JSON file is absent
(run the R script to populate it). The ``require_r_dcdh`` fixture is
intentionally NOT used here so the parity regression suite remains
active in CI even when DIDmultiplegtDYN is unavailable.

Tolerances follow the existing ``test_csdid_ported.py`` convention:
``rtol=1e-4`` for pure-direction point estimates, looser tolerances for
SEs and mixed-direction scenarios (see class docstring).
"""

import json
from pathlib import Path
from typing import List

import pandas as pd
import pytest

from diff_diff import ChaisemartinDHaultfoeuille

# ---------------------------------------------------------------------------
# Fixtures
# ---------------------------------------------------------------------------

GOLDEN_VALUES_PATH = (
    Path(__file__).parents[1] / "benchmarks" / "data" / "dcdh_dynr_golden_values.json"
)


def _as_dict(slot):
    """Coerce a fixture slot to a dict for type-stable iteration.

    R's ``jsonlite::toJSON`` serializes empty named lists (``list()``)
    as JSON arrays (``[]``); populated named lists serialize as JSON
    objects (``{}``). The dCDH `predict_het` extractors return empty
    lists for the no-rows case (e.g., scenarios without ``placebo`` set
    have empty ``placebo_predict_het`` slots), so consumers iterating
    ``.items()`` need to handle the array form. This helper coerces
    empty list ``[]``, ``None`` (missing slot), or any non-dict value
    to ``{}``; populated dicts pass through unchanged. Used by all
    parity-test classes that read the optional ``placebo_predict_het``
    / ``placebo_horizons`` slots.
    """
    return slot if isinstance(slot, dict) else {}


@pytest.fixture(scope="module")
def golden_values():
    """
    Load R DIDmultiplegtDYN golden values from the committed JSON file.

    Skips ONLY if the JSON file is absent (R is not needed to run this
    fixture — the committed JSON is the source of truth). To regenerate
    the JSON, run::

        Rscript benchmarks/R/generate_dcdh_dynr_test_values.R

    which requires R + the `DIDmultiplegtDYN` and `jsonlite` packages.
    """
    if not GOLDEN_VALUES_PATH.exists():
        pytest.skip(
            f"dCDH golden values file not found at {GOLDEN_VALUES_PATH}; "
            "run: Rscript benchmarks/R/generate_dcdh_dynr_test_values.R"
        )
    with open(GOLDEN_VALUES_PATH) as f:
        return json.load(f)["scenarios"]


def _golden_to_df(data_dict: dict) -> pd.DataFrame:
    """Reconstruct a panel DataFrame from a golden-values data export."""
    return pd.DataFrame(
        {
            "group": data_dict["group"],
            "period": data_dict["period"],
            "treatment": data_dict["treatment"],
            "outcome": data_dict["outcome"],
        }
    )


def _fit_dcdh_l1(df: pd.DataFrame):
    """Fit ChaisemartinDHaultfoeuille at l=1 and return the results object."""
    est = ChaisemartinDHaultfoeuille()
    return est.fit(
        df,
        outcome="outcome",
        group="group",
        time="period",
        treatment="treatment",
    )


# ---------------------------------------------------------------------------
# Parity tests — one per scenario in the golden values JSON
# ---------------------------------------------------------------------------


class TestDCDHDynRParity:
    """
    Each scenario asserts |python - r| within tolerance for DID_1.

    Tolerances by scenario type:

    - **Pure-direction scenarios** (``joiners_only``, ``leavers_only``):
      tight ``POINT_RTOL = 1e-4`` (0.01%). Python and R agree exactly when
      every switching cohort uses controls of a single baseline type.
    - **Hand-calculable worked example**: exact arithmetic agreement
      (``abs=1e-6``).
    - **Mixed-direction scenarios** (``single_switch_mixed``,
      ``mixed_single_switch``): looser ``MIXED_POINT_RTOL = 0.025`` (2.5%).
      Python implements the AER 2020 Theorem 3 formula with period-based
      stable-control sets (``stable_0(t)`` is any cell with
      ``D_{,t-1}=D_{,t}=0``); R ``DIDmultiplegtDYN`` implements the
      dynamic-companion-paper formula with cohort-based control sets
      (additionally requires ``D_{g,1}`` to match the side). The two
      formulas agree exactly when no group's baseline is on the opposite
      side of where it currently sits — i.e., always for pure-direction
      panels and for the worked example. They disagree by O(1%) when both
      joiners and leavers exist AND some joiners' post-switch cells could
      serve as leavers' controls (or vice versa). This is documented in
      ``REGISTRY.md`` under ``ChaisemartinDHaultfoeuille`` as a
      ``**Note (deviation from R DIDmultiplegtDYN):**``.
    """

    POINT_RTOL = 1e-4
    MIXED_POINT_RTOL = 0.025  # See class docstring for justification.
    SE_RTOL = 1e-3
    MIXED_SE_RTOL = 0.10  # SE is also affected by the cohort/period control set choice.

    def test_parity_single_switch_mixed(self, golden_values):
        scenario = golden_values.get("single_switch_mixed")
        if scenario is None:
            pytest.skip("scenario 'single_switch_mixed' not in golden values")
        df = _golden_to_df(scenario["data"])
        results = _fit_dcdh_l1(df)
        r_results = scenario["results"]
        # Mixed scenario: see class docstring re. cohort vs. period control sets.
        # Point estimate is checked at 2.5% rtol; SE comparison is intentionally
        # omitted because the SE depends on the same cohort-vs-period design
        # choice and is fundamentally apples-to-oranges between the two
        # implementations on mixed-direction panels.
        assert results.overall_att == pytest.approx(
            r_results["overall_att"], rel=self.MIXED_POINT_RTOL
        )

    PURE_DIRECTION_SE_RTOL = 0.05  # 5% rtol on pure-direction scenarios after the full IF fix

    def test_parity_joiners_only(self, golden_values):
        scenario = golden_values.get("joiners_only")
        if scenario is None:
            pytest.skip("scenario 'joiners_only' not in golden values")
        df = _golden_to_df(scenario["data"])
        results = _fit_dcdh_l1(df)
        r_results = scenario["results"]
        # Pure-direction scenario: tight point-estimate parity expected
        # (cell counts and the full IF formula produce identical answers
        # to R DIDmultiplegtDYN under pure-direction conditions).
        assert results.overall_att == pytest.approx(r_results["overall_att"], rel=self.POINT_RTOL)
        # SE parity: looser 5% rtol because Python's cohort-recentered
        # plug-in formula and R's may use slightly different small-sample
        # corrections. The full Lambda^G_{g,l=1} influence function fix
        # narrowed this gap from ~18% to ~3% on the test scenarios.
        assert results.overall_se == pytest.approx(
            r_results["overall_se"], rel=self.PURE_DIRECTION_SE_RTOL
        )

    def test_parity_leavers_only(self, golden_values):
        scenario = golden_values.get("leavers_only")
        if scenario is None:
            pytest.skip("scenario 'leavers_only' not in golden values")
        df = _golden_to_df(scenario["data"])
        results = _fit_dcdh_l1(df)
        r_results = scenario["results"]
        # Pure-direction scenario: tight point-estimate + 5% SE parity
        assert results.overall_att == pytest.approx(r_results["overall_att"], rel=self.POINT_RTOL)
        assert results.overall_se == pytest.approx(
            r_results["overall_se"], rel=self.PURE_DIRECTION_SE_RTOL
        )

    def test_parity_mixed_single_switch(self, golden_values):
        scenario = golden_values.get("mixed_single_switch")
        if scenario is None:
            pytest.skip("scenario 'mixed_single_switch' not in golden values")
        df = _golden_to_df(scenario["data"])
        results = _fit_dcdh_l1(df)
        r_results = scenario["results"]
        # Mixed scenario: see class docstring re. cohort vs. period control sets.
        assert results.overall_att == pytest.approx(
            r_results["overall_att"], rel=self.MIXED_POINT_RTOL
        )

    def test_parity_hand_calculable_worked_example(self, golden_values):
        """
        Cross-check the 4-group worked example panel against R.

        This is the panel where the Python implementation produces
        DID_M = 2.5 exactly. R DIDmultiplegtDYN should agree (modulo
        small numerical differences in how SE is computed).
        """
        scenario = golden_values.get("hand_calculable_worked_example")
        if scenario is None:
            pytest.skip("scenario 'hand_calculable_worked_example' not in golden values")
        df = _golden_to_df(scenario["data"])
        results = _fit_dcdh_l1(df)
        r_results = scenario["results"]
        # Tight tolerance for this exact-arithmetic case
        assert results.overall_att == pytest.approx(r_results["overall_att"], abs=1e-6)


# ---------------------------------------------------------------------------
# Phase 2: Multi-horizon parity tests
# ---------------------------------------------------------------------------


def _fit_dcdh_multi(df: pd.DataFrame, L_max: int):
    """Fit ChaisemartinDHaultfoeuille with L_max and return results."""
    est = ChaisemartinDHaultfoeuille()
    return est.fit(
        df,
        outcome="outcome",
        group="group",
        time="period",
        treatment="treatment",
        L_max=L_max,
    )


class TestDCDHDynRParityMultiHorizon:
    """
    Multi-horizon parity tests against R DIDmultiplegtDYN ``effects > 1``.

    Each scenario asserts per-horizon |python - r| within tolerance for
    DID_l at each horizon l = 1..L_max. The multi-horizon scenarios
    are generated by the R script with ``effects = 3`` or ``effects = 5``
    and ``placebo > 0``.

    Tolerances follow the same convention as the l=1 parity tests:
    pure-direction at 1e-4 for point estimates, looser for SEs and
    mixed-direction panels.
    """

    POINT_RTOL = 1e-4
    MIXED_POINT_RTOL = 0.025
    # SE tolerance: the cell-count vs obs-count weighting deviation
    # (documented in REGISTRY.md) compounds across horizons and panels.
    # 10% covers the observed gap on all pure-direction multi-horizon
    # scenarios.
    SE_RTOL = 0.10

    def _check_multi_horizon(self, golden_values, scenario_name, L_max, rtol):
        scenario = golden_values.get(scenario_name)
        if scenario is None:
            pytest.skip(f"scenario '{scenario_name}' not in golden values")
        r_effects = scenario["results"].get("effects", {})
        if not r_effects:
            pytest.skip(f"scenario '{scenario_name}' has no multi-horizon effects")
        df = _golden_to_df(scenario["data"])
        results = _fit_dcdh_multi(df, L_max=L_max)

        for horizon_str, r_data in r_effects.items():
            horizon = int(horizon_str)
            if horizon not in results.event_study_effects:
                continue
            py_eff = results.event_study_effects[horizon]["effect"]
            r_eff = r_data["overall_att"]
            assert py_eff == pytest.approx(r_eff, rel=rtol), (
                f"Horizon {horizon}: Python DID_{horizon}={py_eff:.6f} vs "
                f"R DID_{horizon}={r_eff:.6f} (rtol={rtol})"
            )

    def _check_multi_horizon_se(self, golden_values, scenario_name, L_max, se_rtol):
        """Check per-horizon SE and placebo SE against R golden values."""
        scenario = golden_values.get(scenario_name)
        if scenario is None:
            pytest.skip(f"scenario '{scenario_name}' not in golden values")
        df = _golden_to_df(scenario["data"])
        results = _fit_dcdh_multi(df, L_max=L_max)

        # Per-horizon SE
        r_effects = scenario["results"].get("effects", {})
        for horizon_str, r_data in r_effects.items():
            horizon = int(horizon_str)
            if horizon not in results.event_study_effects:
                continue
            py_se = results.event_study_effects[horizon]["se"]
            r_se = r_data.get("overall_se")
            if r_se is not None and r_se > 0:
                assert py_se == pytest.approx(r_se, rel=se_rtol), (
                    f"Horizon {horizon}: Python SE={py_se:.6f} vs "
                    f"R SE={r_se:.6f} (rtol={se_rtol})"
                )

        # Placebo SE
        r_placebos = scenario["results"].get("placebos", {})
        if r_placebos and results.placebo_event_study:
            for lag_str, r_pl_data in r_placebos.items():
                lag = int(lag_str)
                neg_key = -lag
                if neg_key not in results.placebo_event_study:
                    continue
                py_pl_se = results.placebo_event_study[neg_key]["se"]
                r_pl_se = r_pl_data.get("se")
                if r_pl_se is not None and r_pl_se > 0:
                    assert py_pl_se == pytest.approx(r_pl_se, rel=se_rtol), (
                        f"Placebo lag {lag}: Python SE={py_pl_se:.6f} vs "
                        f"R SE={r_pl_se:.6f} (rtol={se_rtol})"
                    )

    def test_parity_joiners_only_multi_horizon(self, golden_values):
        self._check_multi_horizon(
            golden_values, "joiners_only_multi_horizon", L_max=3, rtol=self.POINT_RTOL
        )

    def test_parity_joiners_only_multi_horizon_se(self, golden_values):
        self._check_multi_horizon_se(
            golden_values, "joiners_only_multi_horizon", L_max=3, se_rtol=self.SE_RTOL
        )

    def test_parity_leavers_only_multi_horizon(self, golden_values):
        self._check_multi_horizon(
            golden_values, "leavers_only_multi_horizon", L_max=3, rtol=self.POINT_RTOL
        )

    def test_parity_leavers_only_multi_horizon_se(self, golden_values):
        self._check_multi_horizon_se(
            golden_values, "leavers_only_multi_horizon", L_max=3, se_rtol=self.SE_RTOL
        )

    def test_parity_mixed_single_switch_multi_horizon(self, golden_values):
        self._check_multi_horizon(
            golden_values,
            "mixed_single_switch_multi_horizon",
            L_max=5,
            rtol=self.MIXED_POINT_RTOL,
        )

    def test_parity_joiners_only_long_multi_horizon(self, golden_values):
        self._check_multi_horizon(
            golden_values,
            "joiners_only_long_multi_horizon",
            L_max=5,
            rtol=self.POINT_RTOL,
        )

    def test_parity_joiners_only_long_multi_horizon_se(self, golden_values):
        # Long panel (L_max=5): the cell-count vs obs-count weighting
        # deviation compounds at far horizons where the panel thins.
        # Use wider tolerance than the shorter (L_max=3) scenarios.
        self._check_multi_horizon_se(
            golden_values, "joiners_only_long_multi_horizon", L_max=5, se_rtol=0.15
        )


# ---------------------------------------------------------------------------
# Phase 3: Covariate and trend parity tests
# ---------------------------------------------------------------------------


def _golden_to_df_with_covariates(data_dict: dict) -> pd.DataFrame:
    """Reconstruct a panel DataFrame including covariate columns."""
    cols = {
        "group": data_dict["group"],
        "period": data_dict["period"],
        "treatment": data_dict["treatment"],
        "outcome": data_dict["outcome"],
    }
    if "X1" in data_dict:
        cols["X1"] = data_dict["X1"]
    return pd.DataFrame(cols)


def _golden_to_df_with_extra(
    data_dict: dict, extra_cols: List[str]
) -> pd.DataFrame:
    """Reconstruct a panel DataFrame including arbitrary extra columns.

    Used for scenarios that ship non-covariate side columns (e.g.,
    ``state`` for ``trends_nonparam``).
    """
    cols = {
        "group": data_dict["group"],
        "period": data_dict["period"],
        "treatment": data_dict["treatment"],
        "outcome": data_dict["outcome"],
    }
    for c in extra_cols:
        if c in data_dict:
            cols[c] = data_dict[c]
    return pd.DataFrame(cols)


class TestDCDHDynRParityPhase3:
    """
    Phase 3 parity tests: covariates (DID^X) and linear trends (DID^{fd}).

    Tests that the Python implementation matches R ``did_multiplegt_dyn``
    with ``controls`` and ``trends_lin`` options on identical data.

    Tolerances are wider than Phase 1/2 because the covariate and trend
    adjustments involve additional OLS steps that may amplify the
    cell-count vs obs-count weighting deviation documented in REGISTRY.md.
    """

    # Controls-only: observed gap 0.15%-0.26% (from OLS residualization).
    # Trends-only: exact (0.0000%) at both horizons after cumulation fix.
    # Combined: observed gap 0.30%-0.59% (from OLS residualization only).
    # SE: 3-5% from cell-count weighting; 12-18% for cumulated SEs.
    POINT_RTOL = 0.01    # 1% for controls (observed: 0.26%)
    SE_RTOL = 0.20       # 20% for SE (cell-count weighting + cumulation)

    def _check_phase3_scenario(
        self, golden_values, scenario_name, L_max, controls=None,
        trends_linear=None, point_rtol=None, se_rtol=None,
    ):
        scenario = golden_values.get(scenario_name)
        if scenario is None:
            pytest.skip(f"scenario {scenario_name!r} not in golden values")

        df = _golden_to_df_with_covariates(scenario["data"])
        est = ChaisemartinDHaultfoeuille()
        results = est.fit(
            df, outcome="outcome", group="group", time="period",
            treatment="treatment", L_max=L_max,
            controls=controls, trends_linear=trends_linear,
        )
        r_results = scenario["results"]
        rtol = point_rtol or self.POINT_RTOL
        se_tol = se_rtol or self.SE_RTOL

        # When trends_linear is active, R returns cumulated level effects
        # (delta^{fd}_l), not second-differences (DID^{fd}_l). Compare
        # against linear_trends_effects (cumulated) instead of
        # event_study_effects (second-differences).
        if trends_linear:
            py_effects = results.linear_trends_effects
            assert py_effects is not None, "linear_trends_effects is None"
        else:
            py_effects = results.event_study_effects

        # Check per-horizon effects
        for h_str, r_eff in r_results.get("effects", {}).items():
            h = int(h_str)
            assert h in py_effects, (
                f"Horizon {h} missing from Python results"
            )
            py_eff = py_effects[h]["effect"]
            assert py_eff == pytest.approx(
                r_eff["overall_att"], rel=rtol
            ), f"h={h}: Python={py_eff:.4f} vs R={r_eff['overall_att']:.4f}"

            # SE comparison (wider tolerance)
            py_se = py_effects[h]["se"]
            r_se = r_eff["overall_se"]
            if py_se > 0 and r_se > 0:
                assert py_se == pytest.approx(
                    r_se, rel=se_tol
                ), f"h={h} SE: Python={py_se:.4f} vs R={r_se:.4f}"

    def test_parity_joiners_only_controls(self, golden_values):
        """DID^X with controls vs R did_multiplegt_dyn(..., controls='X1').

        Observed gap: 0.15% at h=1, 0.26% at h=2. Deterministic on
        identical data - the small gap is from the documented cell-count
        vs obs-count weighting deviation in REGISTRY.md.
        """
        self._check_phase3_scenario(
            golden_values, "joiners_only_controls", L_max=2,
            controls=["X1"],
            point_rtol=self.POINT_RTOL,
        )

    def test_parity_joiners_only_trends_lin(self, golden_values):
        """DID^{fd} with trends_linear vs R did_multiplegt_dyn(..., trends_lin=TRUE).

        Exact match (0.0000%) at both horizons after per-group cumulation fix.
        """
        self._check_phase3_scenario(
            golden_values, "joiners_only_trends_lin", L_max=2,
            trends_linear=True,
            point_rtol=1e-4,  # exact match
        )

    def test_parity_joiners_only_controls_trends_lin(self, golden_values):
        """DID^{X,fd} with controls + trends vs R.

        Observed gap: 0.30%-0.59% (from OLS residualization step only;
        the trends cumulation is now exact after per-group cumulation fix).
        """
        self._check_phase3_scenario(
            golden_values, "joiners_only_controls_trends_lin", L_max=2,
            controls=["X1"], trends_linear=True,
            point_rtol=self.POINT_RTOL,
        )


# ---------------------------------------------------------------------------
# Phase 3 extension: by_path per-path event-study disaggregation parity
# ---------------------------------------------------------------------------


class TestDCDHDynRParityByPath:
    """
    Parity tests for ``by_path`` against R DIDmultiplegtDYN.

    R's ``did_multiplegt_dyn(..., by_path=k)`` returns ``k`` per-path
    result slots (``res$by_level_1``, ``res$by_level_2``, ...) keyed by
    frequency rank with path labels in ``res$by_levels``. The generator
    script captures these under ``results.by_path`` as an ordered list.

    Paths are matched by **tuple label via set equality**, not by
    frequency rank: R's tied-frequency tiebreak convention is
    undocumented, while Python's is lexicographic. Set-equality catches
    silent divergence where Python and R select non-identical top-k
    path sets under equal frequencies, and the symmetric-difference
    diagnostic in the assertion message makes the failure actionable.

    Per-path switcher counts are hard-asserted **before** SE
    comparison so that window-eligibility divergences surface
    explicitly rather than as apples-to-oranges SE gaps.
    """

    # Point-estimate tolerance: 1e-9 locks in the observed exact-match
    # property on the committed goldens (measured rtol ~1e-11 across
    # all (path, horizon) cells in both scenarios, well inside R's
    # 10-digit JSON rounding envelope). This matches the claim in
    # REGISTRY.md and CHANGELOG.md that per-path point estimates agree
    # with R exactly. Any material regression on point parity will
    # trip this; loosening to 2.5% would silently accept a ~6 orders
    # of magnitude regression.
    POINT_RTOL = 1e-9
    # SE tolerance: 12% sits between Phase 2 multi-horizon SE_RTOL
    # (0.10 for short panels) and the long-panel widening (0.15). Per-
    # path slices have ~10-40 switchers - comparable to Phase 2's long-
    # panel pure-direction slices where the cell-count vs obs-count
    # weighting envelope (documented in REGISTRY.md) widens from 0.10
    # to 0.15 as horizon length compounds the small-sample correction.
    SE_RTOL = 0.12

    def _path_key_from_r_label(self, r_label: str):
        return tuple(int(x) for x in r_label.split(","))

    def _compare_by_path(self, scenario, by_path, L_max, point_rtol, se_rtol):
        import warnings

        df = _golden_to_df(scenario["data"])
        est = ChaisemartinDHaultfoeuille(drop_larger_lower=False, by_path=by_path)
        # Suppress the per-(path, horizon) degenerate-cohort UserWarning
        # emitted by the foundation when a path subset produces an
        # identically-zero centered IF.
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", UserWarning)
            results = est.fit(
                df,
                outcome="outcome",
                group="group",
                time="period",
                treatment="treatment",
                L_max=L_max,
            )

        r_by_path = scenario["results"]["by_path"]
        assert results.path_effects is not None

        py_keys = set(results.path_effects.keys())
        r_keys = {self._path_key_from_r_label(e["path"]) for e in r_by_path}
        assert py_keys == r_keys, (
            f"Path-set mismatch between Python and R.\n"
            f"  Python only: {py_keys - r_keys}\n"
            f"  R only:      {r_keys - py_keys}"
        )

        for r_path_entry in r_by_path:
            path_key = self._path_key_from_r_label(r_path_entry["path"])
            py_path = results.path_effects[path_key]

            # Assert the public frequency_rank contract matches R. Both
            # committed scenarios are constructed with unique path
            # frequencies (scenario 13 via mixed_single_switch pattern,
            # scenario 14 via deterministic counts 40/25/10/5) so rank
            # ordering is unambiguous and must agree; a regression in
            # path ranking or top-k tiebreak handling should fail here
            # even if the selected path set and per-path effects remain
            # correct.
            assert py_path["frequency_rank"] == r_path_entry["frequency_rank"], (
                f"path={path_key}: frequency_rank mismatch "
                f"py={py_path['frequency_rank']} vs r={r_path_entry['frequency_rank']}"
            )

            for h_str, r_h in r_path_entry["horizons"].items():
                h = int(h_str)
                assert (
                    h in py_path["horizons"]
                ), f"path={path_key}: horizon {h} missing from Python path_effects"
                py_h = py_path["horizons"][h]

                # Per-path switcher count sanity check BEFORE SE comparison:
                # if R and Python disagree on which groups land in this
                # path's window-eligible set at horizon h, SE gaps become
                # apples-to-oranges and point-estimate gaps are masked.
                assert py_h["n_obs"] == int(r_h["n_switchers"]), (
                    f"path={path_key} h={h}: switcher-count mismatch "
                    f"py={py_h['n_obs']} vs r={int(r_h['n_switchers'])} "
                    f"- window-eligibility divergence; investigate before "
                    f"comparing SE."
                )

                assert py_h["effect"] == pytest.approx(r_h["effect"], rel=point_rtol), (
                    f"path={path_key} h={h}: "
                    f"py={py_h['effect']:.4f} vs r={r_h['effect']:.4f}"
                )

                # Assert matching finite/missing state BEFORE the numeric
                # tolerance check so that a silent regression to 0/NaN on
                # one side (while R stays finite, or vice versa) fails the
                # test rather than skipping this cell.
                import math

                py_se = py_h["se"]
                r_se = r_h["se"]
                py_finite_positive = math.isfinite(py_se) and py_se > 0.0
                r_finite_positive = math.isfinite(r_se) and r_se > 0.0
                assert py_finite_positive == r_finite_positive, (
                    f"path={path_key} h={h} SE state mismatch "
                    f"(py_se={py_se}, r_se={r_se}): one side is "
                    f"finite+positive while the other is 0/NaN/inf. "
                    f"Verify the path subset produces matching "
                    f"variance-identifiability on both sides."
                )
                if py_finite_positive and r_finite_positive:
                    assert py_se == pytest.approx(r_se, rel=se_rtol), (
                        f"path={path_key} h={h} SE: "
                        f"py={py_se:.4f} vs r={r_se:.4f}"
                    )

    def test_parity_mixed_single_switch_by_path(self, golden_values):
        """2-path basic case: mixed_single_switch at L_max=3, by_path=2."""
        scenario = golden_values.get("mixed_single_switch_by_path")
        if scenario is None:
            pytest.skip("scenario 'mixed_single_switch_by_path' not in golden values")
        self._compare_by_path(
            scenario,
            by_path=2,
            L_max=3,
            point_rtol=self.POINT_RTOL,
            se_rtol=self.SE_RTOL,
        )

    def test_parity_multi_path_reversible_by_path(self, golden_values):
        """Top-k ranking case: 4 observed paths, by_path=3 selects top-3."""
        scenario = golden_values.get("multi_path_reversible_by_path")
        if scenario is None:
            pytest.skip("scenario 'multi_path_reversible_by_path' not in golden values")
        self._compare_by_path(
            scenario,
            by_path=3,
            L_max=3,
            point_rtol=self.POINT_RTOL,
            se_rtol=self.SE_RTOL,
        )


class TestDCDHDynRParityByPathPlacebo:
    """
    Parity tests for ``by_path + placebo`` against R DIDmultiplegtDYN.

    R's ``did_multiplegt_dyn(..., by_path=k, placebo=N)`` re-runs the
    estimator per path; each ``res$by_level_i$results$Placebos`` row
    holds the per-path backward-horizon placebo. The R generator
    captures these under the same ``horizons`` dict on each per-path
    entry, with negative-int string keys ("-1", "-2", ...) parallel
    to the existing positive-keyed event-study horizons.

    Per-path placebos inherit the same cross-path cohort-sharing SE
    deviation from R that ``path_effects`` shows (full-panel cohort-
    centered plug-in vs R's per-path re-run); the
    ``multi_path_reversible_by_path_placebo`` scenario is constructed
    on a deterministic, single-path-per-cohort DGP so analytical SE
    tracks R within the same Phase-2 envelope used by
    ``TestDCDHDynRParityByPath`` for positive horizons.

    R's placebo iteration is conditional on per-path cohort
    eligibility: paths whose smallest F_g cohort has backward index
    ``F_g - 1 - lag < 0`` produce fewer placebo rows than the user
    requested. The test iterates over the rows R actually produced
    (negative-string keys present in ``r_path_entry["horizons"]``)
    rather than over the requested ``placebo`` parameter.
    """

    POINT_RTOL = 1e-9
    SE_RTOL = 0.12

    def _path_key_from_r_label(self, r_label: str):
        return tuple(int(x) for x in r_label.split(","))

    def test_parity_multi_path_reversible_by_path_placebo(self, golden_values):
        """Per-path placebos R-parity on the cohort-clean deterministic
        DGP. Bundle of (path, lag) cells produced by R is compared row-
        by-row to Python's ``path_placebo_event_study``.
        """
        import math
        import warnings

        scenario = golden_values.get("multi_path_reversible_by_path_placebo")
        if scenario is None:
            pytest.skip(
                "scenario 'multi_path_reversible_by_path_placebo' not in golden values"
            )

        df = _golden_to_df(scenario["data"])
        est = ChaisemartinDHaultfoeuille(
            drop_larger_lower=False, by_path=3, placebo=True
        )
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", UserWarning)
            results = est.fit(
                df,
                outcome="outcome",
                group="group",
                time="period",
                treatment="treatment",
                L_max=3,
            )

        r_by_path = scenario["results"]["by_path"]
        assert results.path_placebo_event_study is not None

        # Same path-set equality check as positive horizons; R-parity
        # over a different surface but the path enumeration must match.
        py_keys = set(results.path_placebo_event_study.keys())
        r_keys = {self._path_key_from_r_label(e["path"]) for e in r_by_path}
        assert py_keys == r_keys, (
            f"Path-set mismatch between Python and R placebos.\n"
            f"  Python only: {py_keys - r_keys}\n"
            f"  R only:      {r_keys - py_keys}"
        )

        for r_path_entry in r_by_path:
            path_key = self._path_key_from_r_label(r_path_entry["path"])
            py_lag_dict = results.path_placebo_event_study[path_key]

            # Iterate over R's negative-keyed horizons. R generator emits
            # `as.character(-h)` keys ("-1", "-2", ...); convert via int().
            for h_str, r_h in r_path_entry["horizons"].items():
                h = int(h_str)
                if h >= 0:
                    continue  # positive horizons covered by TestDCDHDynRParityByPath

                assert h in py_lag_dict, (
                    f"path={path_key}: placebo lag {h} present in R goldens "
                    f"but missing from Python path_placebo_event_study"
                )
                py_h = py_lag_dict[h]

                assert py_h["n_obs"] == int(r_h["n_switchers"]), (
                    f"path={path_key} lag={h}: switcher-count mismatch "
                    f"py={py_h['n_obs']} vs r={int(r_h['n_switchers'])} "
                    f"- per-path placebo eligibility divergence; investigate "
                    f"before comparing SE."
                )

                assert py_h["effect"] == pytest.approx(
                    r_h["effect"], rel=self.POINT_RTOL
                ), (
                    f"path={path_key} lag={h}: "
                    f"py={py_h['effect']:.4f} vs r={r_h['effect']:.4f}"
                )

                py_se = py_h["se"]
                r_se = r_h["se"]
                py_finite_positive = math.isfinite(py_se) and py_se > 0.0
                r_finite_positive = math.isfinite(r_se) and r_se > 0.0
                assert py_finite_positive == r_finite_positive, (
                    f"path={path_key} lag={h} placebo SE state mismatch "
                    f"(py_se={py_se}, r_se={r_se})"
                )
                if py_finite_positive and r_finite_positive:
                    assert py_se == pytest.approx(r_se, rel=self.SE_RTOL), (
                        f"path={path_key} lag={h} placebo SE: "
                        f"py={py_se:.4f} vs r={r_se:.4f}"
                    )


class TestDCDHDynRParityByPathControls:
    """
    Parity tests for ``by_path + controls`` (DID^X residualization)
    against R DIDmultiplegtDYN 2.3.3.

    R's ``did_multiplegt_dyn(..., by_path=k, controls="X1")`` re-runs
    the estimator per path with a path-restricted subsample (path's
    switchers + same-baseline not-yet-treated controls); see the
    canonical R source at ``R/R/did_multiplegt_dyn.R`` lines 393-411
    where the per-path dispatch loop calls ``did_multiplegt_main()``
    once per path with ``df_main`` filtered to that path's groups
    plus same-baseline not-yet-treated controls. Our architecture
    residualizes once on the full panel before path enumeration.

    On the ``multi_path_reversible`` DGP, every switcher shares
    baseline ``D_{g,1}=0``, so R's per-path control pool reduces to
    the global control pool we use — and the per-baseline OLS
    residualization coefficients coincide on the one-observation-per-
    ``(g, t)`` panel this fixture produces. Per-path point estimates
    then match R bit-exactly (measured rtol ~1e-11). On
    multi-observation-per-cell panels the library's equal-cell-
    weighting first stage diverges from R's ``N_gt``-weighted first
    stage per the existing DID^X cell-weighting deviation documented
    in ``docs/methodology/REGISTRY.md`` "Note (Phase 3 DID^X covariate
    adjustment)" — that deviation is independent of the by_path lift.
    Per-path SE inherits the documented cross-path cohort-sharing
    deviation (Phase 2 envelope, ~6.5% rtol on this scenario).

    **On multi-baseline switcher panels** the residualization
    coefficients vary per path under R's per-path call, producing a
    point-estimate deviation between Python and R. The fixture
    intentionally restricts to the single-baseline scenario to keep
    the parity claim tight, and the production estimator emits a
    ``UserWarning`` whenever ``by_path + controls`` is fit on a
    multi-baseline panel so practitioners do not silently consume
    estimates that disagree with R.
    """

    POINT_RTOL = 1e-9
    SE_RTOL = 0.12

    def _path_key_from_r_label(self, r_label: str):
        return tuple(int(x) for x in r_label.split(","))

    def test_parity_multi_path_reversible_by_path_controls(self, golden_values):
        """3-path case with covariate residualization: by_path=3, controls=X1."""
        import math
        import warnings

        scenario = golden_values.get("multi_path_reversible_by_path_controls")
        if scenario is None:
            pytest.skip(
                "scenario 'multi_path_reversible_by_path_controls' not in golden values"
            )

        df = _golden_to_df_with_covariates(scenario["data"])
        est = ChaisemartinDHaultfoeuille(drop_larger_lower=False, by_path=3)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", UserWarning)
            results = est.fit(
                df,
                outcome="outcome",
                group="group",
                time="period",
                treatment="treatment",
                controls=["X1"],
                L_max=3,
            )

        r_by_path = scenario["results"]["by_path"]
        assert results.path_effects is not None

        py_keys = set(results.path_effects.keys())
        r_keys = {self._path_key_from_r_label(e["path"]) for e in r_by_path}
        assert py_keys == r_keys, (
            f"Path-set mismatch.\n"
            f"  Python only: {py_keys - r_keys}\n"
            f"  R only:      {r_keys - py_keys}"
        )

        for r_path_entry in r_by_path:
            path_key = self._path_key_from_r_label(r_path_entry["path"])
            py_path = results.path_effects[path_key]

            assert py_path["frequency_rank"] == r_path_entry["frequency_rank"], (
                f"path={path_key}: frequency_rank mismatch "
                f"py={py_path['frequency_rank']} vs r={r_path_entry['frequency_rank']}"
            )

            for h_str, r_h in r_path_entry["horizons"].items():
                h = int(h_str)
                assert (
                    h in py_path["horizons"]
                ), f"path={path_key}: horizon {h} missing from Python path_effects"
                py_h = py_path["horizons"][h]

                assert py_h["n_obs"] == int(r_h["n_switchers"]), (
                    f"path={path_key} h={h}: switcher-count mismatch "
                    f"py={py_h['n_obs']} vs r={int(r_h['n_switchers'])}"
                )

                assert py_h["effect"] == pytest.approx(
                    r_h["effect"], rel=self.POINT_RTOL
                ), (
                    f"path={path_key} h={h}: "
                    f"py={py_h['effect']:.4f} vs r={r_h['effect']:.4f}"
                )

                py_se = py_h["se"]
                r_se = r_h["se"]
                py_finite_positive = math.isfinite(py_se) and py_se > 0.0
                r_finite_positive = math.isfinite(r_se) and r_se > 0.0
                assert py_finite_positive == r_finite_positive, (
                    f"path={path_key} h={h} SE state mismatch "
                    f"(py_se={py_se}, r_se={r_se})"
                )
                if py_finite_positive and r_finite_positive:
                    assert py_se == pytest.approx(r_se, rel=self.SE_RTOL), (
                        f"path={path_key} h={h} SE: "
                        f"py={py_se:.4f} vs r={r_se:.4f}"
                    )


class TestDCDHDynRParityByPathTrendsLinear:
    """
    Parity tests for ``by_path + trends_linear`` (DID^{fd} group-specific
    linear trends) against R DIDmultiplegtDYN 2.3.3.

    R's ``did_multiplegt_dyn(..., by_path=k, trends_lin=TRUE)`` re-runs
    the estimator per path with a path-restricted subsample
    (``R/R/did_multiplegt_dyn.R`` lines 226-257 dispatcher → per-path
    ``did_multiplegt_main()`` call with ``trends_lin=TRUE``). Inside
    ``did_multiplegt_main``, R first-differences the outcome and drops
    ``F_g==2`` switchers and ``time==1`` rows, then loops through
    horizons computing per-horizon DID values and **cumulating** them
    (R's ``Effect_l`` under ``trends_lin=TRUE`` is the cumulated level
    effect ``delta_l``, NOT the raw second-difference DID^{fd}_l —
    verified by the existing parity test at
    ``test_parity_joiners_only_trends_lin`` lines 403-409).

    Python's architecture first-differences once globally before path
    enumeration, then disaggregates per path on the global
    ``multi_horizon_dids`` (raw DID^{fd}_l surfaces on
    ``path_effects``). The new ``path_cumulated_event_study`` field
    surfaces the cumulated ``delta_l`` per path (sum of per-horizon
    DID^{fd} values for the path's switchers, eligibility-restricted at
    each horizon — mirrors the global ``linear_trends_effects``
    cumulation at ``chaisemartin_dhaultfoeuille.py:3340-3398``).

    The fixture ``single_baseline_multi_path_by_path_trends_lin`` is
    designed for clean parity:
    - Single-baseline (all switchers ``D_{g,1}=0``) eliminates the
      multi-baseline divergence pattern (same architectural family as
      ``controls`` — see ``TestDCDHDynRParityByPathControls``)
    - Cohort-single-path (each ``F_g`` maps to exactly one path)
      eliminates the cross-path cohort-sharing deviation that PR #360
      documented for ``path_effects``
    - All ``F_g >= 4`` keeps trends_lin's ``F_g==2`` filter a no-op
      and gives both Python and R 2+ valid pre-window Z values
      (avoids the boundary-case divergence at ``F_g=3``)

    Per-path cumulated point estimates match R bit-exactly (rtol ~1e-9)
    on event horizons under these conditions. Cumulated SE inherits the
    cross-path cohort-sharing deviation amplified by summation across
    horizons; ``CUM_SE_RTOL=0.20`` is widened from the 0.12 used for
    non-cumulated by_path parity (per the conservative upper-bound SE
    formula at ``_compute_path_cumulated_event_study``).

    **Placebo divergence (documented):** under ``trends_linear=True`` +
    ``by_path``, the per-path placebo at lag 1 diverges between Python
    and R for paths whose switchers have minimal pre-window depth
    (e.g., ``F_g=4`` switchers in path 1 of this fixture). R's per-path
    placebo computation re-runs on the path-restricted subsample with
    different control eligibility than Python's global-then-disaggregate
    architecture surfaces. The placebo parity assertion is therefore
    skipped here; placebo + ``trends_linear`` is exercised via internal
    regression in ``TestByPathTrendsLinear`` (finite values, bootstrap
    inheritance) but not pinned to R bit-by-bit. See REGISTRY.md
    "Note (Phase 3 by_path)" sub-paragraph "Per-path linear-trends
    DID^{fd}" for the full deviation documentation.
    """

    POINT_RTOL = 1e-9
    POINT_ATOL = 1e-9
    CUM_SE_RTOL = 0.20

    def _path_key_from_r_label(self, r_label: str):
        return tuple(int(x) for x in r_label.split(","))

    def test_parity_single_baseline_multi_path_trends_lin(self, golden_values):
        """3-path single-baseline trends_lin: by_path=3, trends_linear=True."""
        import math
        import warnings

        scenario = golden_values.get(
            "single_baseline_multi_path_by_path_trends_lin"
        )
        if scenario is None:
            pytest.skip(
                "scenario 'single_baseline_multi_path_by_path_trends_lin' "
                "not in golden values"
            )

        df = _golden_to_df(scenario["data"])
        est = ChaisemartinDHaultfoeuille(
            drop_larger_lower=False, by_path=3, placebo=True
        )
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", UserWarning)
            results = est.fit(
                df,
                outcome="outcome",
                group="group",
                time="period",
                treatment="treatment",
                trends_linear=True,
                L_max=3,
            )

        r_by_path = scenario["results"]["by_path"]
        assert results.path_cumulated_event_study is not None, (
            "path_cumulated_event_study should be populated under "
            "by_path + trends_linear=True"
        )
        py_cum = results.path_cumulated_event_study

        py_keys = set(py_cum.keys())
        r_keys = {self._path_key_from_r_label(e["path"]) for e in r_by_path}
        assert py_keys == r_keys, (
            f"Path-set mismatch.\n"
            f"  Python only: {py_keys - r_keys}\n"
            f"  R only:      {r_keys - py_keys}"
        )

        for r_path_entry in r_by_path:
            path_key = self._path_key_from_r_label(r_path_entry["path"])
            py_path_cum = py_cum[path_key]

            for h_str, r_h in r_path_entry["horizons"].items():
                h = int(h_str)
                # Skip placebo horizons (negative keys) — see class
                # docstring for the documented deviation.
                if h <= 0:
                    continue
                assert h in py_path_cum, (
                    f"path={path_key}: horizon {h} missing from "
                    f"path_cumulated_event_study"
                )
                py_h = py_path_cum[h]

                assert py_h["n_obs"] == int(r_h["n_switchers"]), (
                    f"path={path_key} h={h}: switcher-count mismatch "
                    f"py={py_h['n_obs']} vs r={int(r_h['n_switchers'])}"
                )

                assert py_h["effect"] == pytest.approx(
                    r_h["effect"], rel=self.POINT_RTOL, abs=self.POINT_ATOL
                ), (
                    f"path={path_key} h={h}: cumulated effect mismatch "
                    f"py={py_h['effect']:.6f} vs r={r_h['effect']:.6f}"
                )

                py_se = py_h["se"]
                r_se = r_h["se"]
                py_finite_positive = math.isfinite(py_se) and py_se > 0.0
                r_finite_positive = math.isfinite(r_se) and r_se > 0.0
                assert py_finite_positive == r_finite_positive, (
                    f"path={path_key} h={h} SE state mismatch "
                    f"(py_se={py_se}, r_se={r_se})"
                )
                if py_finite_positive and r_finite_positive:
                    assert py_se == pytest.approx(
                        r_se, rel=self.CUM_SE_RTOL
                    ), (
                        f"path={path_key} h={h} cumulated SE: "
                        f"py={py_se:.4f} vs r={r_se:.4f}"
                    )


class TestDCDHDynRParityByPathTrendsNonparam:
    """
    Parity tests for ``by_path + trends_nonparam`` (state-set trends)
    against R DIDmultiplegtDYN 2.3.3.

    R's ``did_multiplegt_dyn(..., by_path=k, trends_nonparam="state")``
    re-runs the estimator per path with a path-restricted subsample
    (same dispatch loop as other ``by_path`` modes). The
    ``trends_nonparam`` parameter is passed through to
    ``did_multiplegt_main()`` which uses the column as a grouping key
    for control-pool restriction (``did_multiplegt_main`` lines
    414-415: ``grp_cols <- c(grp_cols, trends_nonparam)``). R does
    NOT first-difference and does NOT cumulate under
    ``trends_nonparam``; ``Effect_l`` is a normal per-horizon DID
    with set-restricted controls.

    Python's architecture validates the set-membership column once
    globally and stores ``set_ids_arr`` (aligned with ``all_groups``);
    the ``set_ids`` parameter is threaded through the four per-path IF
    helpers (``_compute_path_effects``,
    ``_compute_path_placebos``, ``_collect_path_bootstrap_inputs``,
    ``_collect_path_placebo_bootstrap_inputs``) so per-path analytical
    SE, bootstrap, placebos, and sup-t bands all consume the
    set-restricted control pool. Per-path point estimates and placebos
    match R bit-exactly under the cohort-single-path
    ``multi_path_reversible`` DGP. Per-path SE inherits the documented
    cross-path cohort-sharing deviation (Phase 2 envelope, ~13% rtol
    on this scenario).
    """

    POINT_RTOL = 1e-9
    SE_RTOL = 0.15

    def _path_key_from_r_label(self, r_label: str):
        return tuple(int(x) for x in r_label.split(","))

    def test_parity_multi_path_reversible_by_path_trends_nonparam(
        self, golden_values
    ):
        """3-path case with state-set trends: by_path=3, trends_nonparam='state'."""
        import math
        import warnings

        scenario = golden_values.get(
            "multi_path_reversible_by_path_trends_nonparam"
        )
        if scenario is None:
            pytest.skip(
                "scenario 'multi_path_reversible_by_path_trends_nonparam' "
                "not in golden values"
            )

        df = _golden_to_df_with_extra(scenario["data"], extra_cols=["state"])
        est = ChaisemartinDHaultfoeuille(
            drop_larger_lower=False, by_path=3, placebo=True
        )
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", UserWarning)
            results = est.fit(
                df,
                outcome="outcome",
                group="group",
                time="period",
                treatment="treatment",
                trends_nonparam="state",
                L_max=3,
            )

        r_by_path = scenario["results"]["by_path"]
        assert results.path_effects is not None

        py_keys = set(results.path_effects.keys())
        r_keys = {self._path_key_from_r_label(e["path"]) for e in r_by_path}
        assert py_keys == r_keys, (
            f"Path-set mismatch.\n"
            f"  Python only: {py_keys - r_keys}\n"
            f"  R only:      {r_keys - py_keys}"
        )

        for r_path_entry in r_by_path:
            path_key = self._path_key_from_r_label(r_path_entry["path"])
            py_path = results.path_effects[path_key]
            py_placebos = (
                results.path_placebo_event_study.get(path_key, {})
                if results.path_placebo_event_study is not None
                else {}
            )

            assert py_path["frequency_rank"] == r_path_entry["frequency_rank"], (
                f"path={path_key}: frequency_rank mismatch "
                f"py={py_path['frequency_rank']} vs r={r_path_entry['frequency_rank']}"
            )

            for h_str, r_h in r_path_entry["horizons"].items():
                h = int(h_str)
                if h > 0:
                    assert h in py_path["horizons"], (
                        f"path={path_key}: horizon {h} missing from "
                        f"Python path_effects"
                    )
                    py_h = py_path["horizons"][h]
                else:
                    assert h in py_placebos, (
                        f"path={path_key}: placebo {h} missing from "
                        f"Python path_placebo_event_study"
                    )
                    py_h = py_placebos[h]

                assert py_h["n_obs"] == int(r_h["n_switchers"]), (
                    f"path={path_key} h={h}: switcher-count mismatch "
                    f"py={py_h['n_obs']} vs r={int(r_h['n_switchers'])}"
                )

                assert py_h["effect"] == pytest.approx(
                    r_h["effect"], rel=self.POINT_RTOL
                ), (
                    f"path={path_key} h={h}: "
                    f"py={py_h['effect']:.6f} vs r={r_h['effect']:.6f}"
                )

                py_se = py_h["se"]
                r_se = r_h["se"]
                py_finite_positive = math.isfinite(py_se) and py_se > 0.0
                r_finite_positive = math.isfinite(r_se) and r_se > 0.0
                assert py_finite_positive == r_finite_positive, (
                    f"path={path_key} h={h} SE state mismatch "
                    f"(py_se={py_se}, r_se={r_se})"
                )
                if py_finite_positive and r_finite_positive:
                    assert py_se == pytest.approx(r_se, rel=self.SE_RTOL), (
                        f"path={path_key} h={h} SE: "
                        f"py={py_se:.4f} vs r={r_se:.4f}"
                    )


class TestDCDHDynRParityByPathNonBinary:
    """
    Parity tests for ``by_path`` + non-binary integer treatment against
    R DIDmultiplegtDYN 2.3.3.

    Wave 3 #8 lift. R's ``did_multiplegt_dyn(..., by_path=k)`` accepts
    discrete non-binary treatment naturally; the per-path string
    encoding (``"0,2,2,2"``) maps bit-for-bit to Python's tuple key
    ``(0, 2, 2, 2)`` whenever every D value renders as a single
    character (nonnegative single-digit states). R's
    ``substr(path_index$path, 1, 1)`` baseline-derivation (in
    ``did_multiplegt_by_path``) breaks for multi-character baselines
    (``D >= 10`` or negative ``D``): it captures only the first
    character of the comma-separated path string. Python's tuple-key
    matching is correct in both regimes — the shipped R-parity
    scenario stays in ``{0, 1, 2}`` so the R bug does not interfere.

    The fixture mirrors Scenario 17's single-baseline custom DGP:
    78 switchers across 3 non-binary paths, 20 never-treated and 20
    always-treated controls, F_g >= 4. Per-path point estimates
    match R bit-exactly. Per-path SE inherits the documented
    cross-path cohort-sharing deviation (Phase 2 envelope, ~5% rtol
    on this scenario; SE_RTOL=0.15 for headroom).
    """

    POINT_RTOL = 1e-9
    # Near-zero placebo (e.g. -0.003) amplifies relative difference;
    # use atol=1e-9 as a fallback so the cumulated relative+absolute
    # `pytest.approx` envelope matches bit-exact event horizons and
    # small-value placebo points consistently.
    POINT_ATOL = 1e-9
    SE_RTOL = 0.15

    def _path_key_from_r_label(self, r_label: str):
        return tuple(int(x) for x in r_label.split(","))

    def test_parity_multi_path_reversible_by_path_non_binary(
        self, golden_values
    ):
        """3-path case with D in {0, 1, 2}: by_path=3, placebo=1."""
        import math
        import warnings

        scenario = golden_values.get(
            "multi_path_reversible_by_path_non_binary"
        )
        if scenario is None:
            pytest.skip(
                "scenario 'multi_path_reversible_by_path_non_binary' "
                "not in golden values"
            )

        df = _golden_to_df(scenario["data"])
        est = ChaisemartinDHaultfoeuille(
            drop_larger_lower=False, by_path=3, placebo=True
        )
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", UserWarning)
            results = est.fit(
                df,
                outcome="outcome",
                group="group",
                time="period",
                treatment="treatment",
                L_max=3,
            )

        r_by_path = scenario["results"]["by_path"]
        assert results.path_effects is not None

        py_keys = set(results.path_effects.keys())
        r_keys = {self._path_key_from_r_label(e["path"]) for e in r_by_path}
        assert py_keys == r_keys, (
            f"Path-set mismatch.\n"
            f"  Python only: {py_keys - r_keys}\n"
            f"  R only:      {r_keys - py_keys}"
        )

        for r_path_entry in r_by_path:
            path_key = self._path_key_from_r_label(r_path_entry["path"])
            py_path = results.path_effects[path_key]
            py_placebos = (
                results.path_placebo_event_study.get(path_key, {})
                if results.path_placebo_event_study is not None
                else {}
            )

            assert py_path["frequency_rank"] == r_path_entry["frequency_rank"], (
                f"path={path_key}: frequency_rank mismatch "
                f"py={py_path['frequency_rank']} vs r={r_path_entry['frequency_rank']}"
            )

            for h_str, r_h in r_path_entry["horizons"].items():
                h = int(h_str)
                if h > 0:
                    assert h in py_path["horizons"], (
                        f"path={path_key}: horizon {h} missing from "
                        f"Python path_effects"
                    )
                    py_h = py_path["horizons"][h]
                else:
                    assert h in py_placebos, (
                        f"path={path_key}: placebo {h} missing from "
                        f"Python path_placebo_event_study"
                    )
                    py_h = py_placebos[h]

                assert py_h["n_obs"] == int(r_h["n_switchers"]), (
                    f"path={path_key} h={h}: switcher-count mismatch "
                    f"py={py_h['n_obs']} vs r={int(r_h['n_switchers'])}"
                )

                assert py_h["effect"] == pytest.approx(
                    r_h["effect"], rel=self.POINT_RTOL, abs=self.POINT_ATOL
                ), (
                    f"path={path_key} h={h}: "
                    f"py={py_h['effect']:.6f} vs r={r_h['effect']:.6f}"
                )

                py_se = py_h["se"]
                r_se = r_h["se"]
                py_finite_positive = math.isfinite(py_se) and py_se > 0.0
                r_finite_positive = math.isfinite(r_se) and r_se > 0.0
                assert py_finite_positive == r_finite_positive, (
                    f"path={path_key} h={h} SE state mismatch "
                    f"(py_se={py_se}, r_se={r_se})"
                )
                if py_finite_positive and r_finite_positive:
                    assert py_se == pytest.approx(r_se, rel=self.SE_RTOL), (
                        f"path={path_key} h={h} SE: "
                        f"py={py_se:.4f} vs r={r_se:.4f}"
                    )


class TestDCDHDynRParityHeterogeneity:
    """Parity tests for global ``predict_het`` against R DIDmultiplegtDYN.

    Wave 5 #11 introduces this as the FIRST ``predict_het`` parity
    baseline in the repo. Anchors the per-path tolerance for
    ``TestDCDHDynRParityByPathHeterogeneity`` since the per-path R call
    is ``did_multiplegt_main(..., predict_het=...)`` per path-restricted
    subsample with no additional numerical loss.

    Fixture: scenario 20 (``multi_path_reversible_predict_het``) — 90
    switchers across 3 reversal paths + 30 never-treated controls,
    binary ``het_x`` with effect = 5.0 + 3.0 * het_x. R is invoked with
    ``dont_drop_larger_lower=TRUE`` so the reversal-path eligibility
    matches Python's ``drop_larger_lower=False`` requirement.
    """

    BETA_RTOL = 1e-6
    BETA_ATOL = 1e-9
    SE_RTOL = 1e-5
    INFERENCE_RTOL = 1e-4  # p_value, conf_int — t-distribution implementation tol

    def test_parity_multi_path_reversible_predict_het(self, golden_values):
        """Global heterogeneity coefficient parity per horizon."""
        import warnings

        scenario = golden_values.get("multi_path_reversible_predict_het")
        if scenario is None:
            pytest.skip(
                "scenario 'multi_path_reversible_predict_het' not in "
                "golden values"
            )

        df = _golden_to_df_with_extra(scenario["data"], extra_cols=["het_x"])
        est = ChaisemartinDHaultfoeuille(drop_larger_lower=False)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", UserWarning)
            results = est.fit(
                df,
                outcome="outcome",
                group="group",
                time="period",
                treatment="treatment",
                L_max=3,
                heterogeneity="het_x",
            )

        assert results.heterogeneity_effects is not None
        r_predict_het = scenario["results"]["predict_het"]

        for h_str, r_h in r_predict_het.items():
            h = int(h_str)
            assert h in results.heterogeneity_effects, (
                f"horizon {h} missing from Python heterogeneity_effects"
            )
            py_h = results.heterogeneity_effects[h]
            assert py_h["beta"] == pytest.approx(
                r_h["beta"], rel=self.BETA_RTOL, abs=self.BETA_ATOL
            ), f"h={h} beta: py={py_h['beta']:.6f} vs r={r_h['beta']:.6f}"
            assert py_h["se"] == pytest.approx(r_h["se"], rel=self.SE_RTOL), (
                f"h={h} se: py={py_h['se']:.6f} vs r={r_h['se']:.6f}"
            )
            # `t_stat = beta / se` is invariant to the Wald-test
            # critical-value distribution; pin it at SE_RTOL so a
            # regression in beta or se surfaces here too.
            assert py_h["t_stat"] == pytest.approx(r_h["t"], rel=self.SE_RTOL), (
                f"h={h} t_stat: py={py_h['t_stat']:.6f} vs r={r_h['t']:.6f}"
            )
            assert int(py_h["n_obs"]) == int(r_h["n_obs"]), (
                f"h={h} n_obs: py={py_h['n_obs']} vs r={r_h['n_obs']}"
            )
            # `p_value` and `conf_int` parity (post-2026-05-15 df threading).
            # `_compute_heterogeneity_test` now passes
            # `df = n_obs - rank(design)` to `safe_inference`, matching
            # R's t-distribution with OLS df (full-rank designs have
            # `rank == n_params`). Pinned at INFERENCE_RTOL = 1e-4
            # because Wald-test critical values come from
            # `scipy.stats.t.ppf` and `t.sf` which are implementation-
            # aligned with R's `qt`/`pt` to ~6 sig figs.
            assert py_h["p_value"] == pytest.approx(
                r_h["p_value"], rel=self.INFERENCE_RTOL
            ), f"h={h} p_value: py={py_h['p_value']:.6e} vs r={r_h['p_value']:.6e}"
            assert py_h["conf_int"][0] == pytest.approx(
                r_h["ci_lo"], rel=self.INFERENCE_RTOL
            ), (
                f"h={h} ci_lo: py={py_h['conf_int'][0]:.6f} "
                f"vs r={r_h['ci_lo']:.6f}"
            )
            assert py_h["conf_int"][1] == pytest.approx(
                r_h["ci_hi"], rel=self.INFERENCE_RTOL
            ), (
                f"h={h} ci_hi: py={py_h['conf_int'][1]:.6f} "
                f"vs r={r_h['ci_hi']:.6f}"
            )


class TestDCDHDynRParityByPathHeterogeneity:
    """Parity tests for ``by_path`` + ``predict_het`` against R DIDmultiplegtDYN.

    Wave 5 #11 lift. R's per-path dispatcher
    (``did_multiplegt_dyn.R:226-257``) re-runs
    ``did_multiplegt_main(..., predict_het=...)`` on each path-restricted
    subsample. Python mirrors via ``_compute_path_heterogeneity_test``,
    which loops ``_compute_heterogeneity_test`` over selected paths
    with a ``path_groups`` filter restricting eligibility to the path's
    switchers. Cohort dummies absorb baseline by construction, so
    multi-baseline switcher panels do not produce R-divergence.

    Fixture: scenario 21 (``multi_path_reversible_by_path_predict_het``)
    — same DGP as scenario 20 so per-path tolerances inherit from
    ``TestDCDHDynRParityHeterogeneity``.
    """

    BETA_RTOL = 1e-6
    BETA_ATOL = 1e-9
    SE_RTOL = 1e-5
    INFERENCE_RTOL = 1e-4

    def _path_key_from_r_label(self, r_label: str):
        return tuple(int(x) for x in r_label.split(","))

    def test_parity_multi_path_reversible_by_path_predict_het(
        self, golden_values
    ):
        """Per-path heterogeneity coefficient parity per (path, horizon)."""
        import warnings

        scenario = golden_values.get(
            "multi_path_reversible_by_path_predict_het"
        )
        if scenario is None:
            pytest.skip(
                "scenario 'multi_path_reversible_by_path_predict_het' "
                "not in golden values"
            )

        df = _golden_to_df_with_extra(scenario["data"], extra_cols=["het_x"])
        est = ChaisemartinDHaultfoeuille(drop_larger_lower=False, by_path=3)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", UserWarning)
            results = est.fit(
                df,
                outcome="outcome",
                group="group",
                time="period",
                treatment="treatment",
                L_max=3,
                heterogeneity="het_x",
            )

        assert results.path_heterogeneity_effects is not None
        r_by_path = scenario["results"]["by_path_predict_het"]

        py_keys = set(results.path_heterogeneity_effects.keys())
        r_keys = {self._path_key_from_r_label(e["path"]) for e in r_by_path}
        assert py_keys == r_keys, (
            f"Path-set mismatch.\n"
            f"  Python only: {py_keys - r_keys}\n"
            f"  R only:      {r_keys - py_keys}"
        )

        for r_path_entry in r_by_path:
            path_key = self._path_key_from_r_label(r_path_entry["path"])
            py_horizons = results.path_heterogeneity_effects[path_key]

            for h_str, r_h in r_path_entry["horizons"].items():
                h = int(h_str)
                assert h in py_horizons, (
                    f"path={path_key}: horizon {h} missing from Python "
                    f"path_heterogeneity_effects"
                )
                py_h = py_horizons[h]
                assert py_h["beta"] == pytest.approx(
                    r_h["beta"], rel=self.BETA_RTOL, abs=self.BETA_ATOL
                ), (
                    f"path={path_key} h={h} beta: "
                    f"py={py_h['beta']:.6f} vs r={r_h['beta']:.6f}"
                )
                assert py_h["se"] == pytest.approx(
                    r_h["se"], rel=self.SE_RTOL
                ), (
                    f"path={path_key} h={h} se: "
                    f"py={py_h['se']:.6f} vs r={r_h['se']:.6f}"
                )
                # `t_stat = beta / se` is invariant to the Wald-test
                # critical-value distribution; pin at SE_RTOL.
                assert py_h["t_stat"] == pytest.approx(
                    r_h["t"], rel=self.SE_RTOL
                ), (
                    f"path={path_key} h={h} t_stat: "
                    f"py={py_h['t_stat']:.6f} vs r={r_h['t']:.6f}"
                )
                assert int(py_h["n_obs"]) == int(r_h["n_obs"]), (
                    f"path={path_key} h={h} n_obs: "
                    f"py={py_h['n_obs']} vs r={r_h['n_obs']}"
                )
                # `p_value` and `conf_int` parity — same df-threading
                # rationale as the global heterogeneity class.
                assert py_h["p_value"] == pytest.approx(
                    r_h["p_value"], rel=self.INFERENCE_RTOL
                ), (
                    f"path={path_key} h={h} p_value: "
                    f"py={py_h['p_value']:.6e} vs r={r_h['p_value']:.6e}"
                )
                assert py_h["conf_int"][0] == pytest.approx(
                    r_h["ci_lo"], rel=self.INFERENCE_RTOL
                ), (
                    f"path={path_key} h={h} ci_lo: "
                    f"py={py_h['conf_int'][0]:.6f} vs r={r_h['ci_lo']:.6f}"
                )
                assert py_h["conf_int"][1] == pytest.approx(
                    r_h["ci_hi"], rel=self.INFERENCE_RTOL
                ), (
                    f"path={path_key} h={h} ci_hi: "
                    f"py={py_h['conf_int'][1]:.6f} vs r={r_h['ci_hi']:.6f}"
                )


class TestDCDHDynRParityByPathHeterogeneityWithPlacebo:
    """Parity tests for ``by_path`` + ``predict_het`` + ``placebo``.

    R-verified: ``did_multiplegt_dyn(by_path, predict_het, placebo)``
    emits per-path heterogeneity OLS results on backward (placebo)
    horizons via R's per-by_level dispatcher
    (``DIDmultiplegtDYN:::did_multiplegt_main`` placebo block at the
    ``effect = matrix(-i, ...)`` rbind site). Python mirrors via
    ``_compute_heterogeneity_test(..., placebo=N)`` which iterates
    forward (1..L_max) and backward (-1..-N) horizons in a single loop.

    Fixture: scenario 22 (``multi_path_reversible_predict_het_with_placebo``)
    — same DGP as scenario 21 (reuses ``d20``) so per-path tolerances
    inherit from ``TestDCDHDynRParityByPathHeterogeneity``.

    R syntax note: scenarios 20/21 use ``predict_het = list("X", c(1,2,3))``
    which is rejected by R when ``placebo > 0`` (forward indices > placebo
    count error). Scenario 22 uses the ``c(-1)`` sentinel that triggers
    "compute het for ALL forward (1..effects) AND ALL placebo
    (1..placebo) positions" — see the R script comment at
    ``benchmarks/R/generate_dcdh_dynr_test_values.R`` scenario 22.
    """

    BETA_RTOL = 1e-6
    BETA_ATOL = 1e-9
    SE_RTOL = 1e-5
    INFERENCE_RTOL = 1e-4

    def _path_key_from_r_label(self, r_label: str):
        return tuple(int(x) for x in r_label.split(","))

    def test_parity_multi_path_reversible_predict_het_with_placebo(
        self, golden_values
    ):
        """Per-path heterogeneity on forward + backward horizons."""
        import warnings

        scenario = golden_values.get(
            "multi_path_reversible_predict_het_with_placebo"
        )
        if scenario is None:
            pytest.skip(
                "scenario 'multi_path_reversible_predict_het_with_placebo' "
                "not in golden values"
            )

        df = _golden_to_df_with_extra(scenario["data"], extra_cols=["het_x"])
        # Python's `placebo=True` triggers all-backward-horizons (1..L_max)
        # placebo computation; the parity check iterates only R's emitted
        # subset (placebo=2 in scenario 22 -> backward horizons -1, -2).
        # Python's extra backward horizons (e.g., -3 if eligible) are
        # allowed and not parity-checked.
        est = ChaisemartinDHaultfoeuille(
            drop_larger_lower=False, by_path=3, placebo=True
        )
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", UserWarning)
            results = est.fit(
                df,
                outcome="outcome",
                group="group",
                time="period",
                treatment="treatment",
                L_max=3,
                heterogeneity="het_x",
            )

        assert results.path_heterogeneity_effects is not None
        r_by_path = scenario["results"]["by_path_predict_het"]

        py_keys = set(results.path_heterogeneity_effects.keys())
        r_keys = {self._path_key_from_r_label(e["path"]) for e in r_by_path}
        assert py_keys == r_keys, (
            f"Path-set mismatch.\n"
            f"  Python only: {py_keys - r_keys}\n"
            f"  R only:      {r_keys - py_keys}"
        )

        for r_path_entry in r_by_path:
            path_key = self._path_key_from_r_label(r_path_entry["path"])
            py_horizons = results.path_heterogeneity_effects[path_key]

            # Forward horizons (positive int keys).
            for h_str, r_h in r_path_entry["horizons"].items():
                h = int(h_str)
                assert h > 0
                assert h in py_horizons, (
                    f"path={path_key}: forward horizon {h} missing"
                )
                py_h = py_horizons[h]
                self._assert_horizon_parity(path_key, h, py_h, r_h)

            # Placebo horizons (negative int keys). R-verified: scenario 22
            # has 2 placebo rows per path (-1, -2); Python mirrors with
            # negative-int keys in path_heterogeneity_effects.
            for h_str, r_h in _as_dict(
                r_path_entry.get("placebo_horizons")
            ).items():
                h = int(h_str)
                assert h < 0
                assert h in py_horizons, (
                    f"path={path_key}: placebo horizon {h} missing "
                    f"from Python path_heterogeneity_effects"
                )
                py_h = py_horizons[h]
                self._assert_horizon_parity(path_key, h, py_h, r_h)

    def _assert_horizon_parity(self, path_key, h, py_h, r_h):
        """Pin all 6 inference fields against R."""
        assert py_h["beta"] == pytest.approx(
            r_h["beta"], rel=self.BETA_RTOL, abs=self.BETA_ATOL
        ), (
            f"path={path_key} h={h} beta: "
            f"py={py_h['beta']:.6f} vs r={r_h['beta']:.6f}"
        )
        assert py_h["se"] == pytest.approx(r_h["se"], rel=self.SE_RTOL), (
            f"path={path_key} h={h} se: "
            f"py={py_h['se']:.6f} vs r={r_h['se']:.6f}"
        )
        assert py_h["t_stat"] == pytest.approx(
            r_h["t"], rel=self.SE_RTOL
        ), (
            f"path={path_key} h={h} t_stat: "
            f"py={py_h['t_stat']:.6f} vs r={r_h['t']:.6f}"
        )
        assert int(py_h["n_obs"]) == int(r_h["n_obs"]), (
            f"path={path_key} h={h} n_obs: "
            f"py={py_h['n_obs']} vs r={r_h['n_obs']}"
        )
        assert py_h["p_value"] == pytest.approx(
            r_h["p_value"], rel=self.INFERENCE_RTOL
        ), (
            f"path={path_key} h={h} p_value: "
            f"py={py_h['p_value']:.6e} vs r={r_h['p_value']:.6e}"
        )
        assert py_h["conf_int"][0] == pytest.approx(
            r_h["ci_lo"], rel=self.INFERENCE_RTOL
        ), (
            f"path={path_key} h={h} ci_lo: "
            f"py={py_h['conf_int'][0]:.6f} vs r={r_h['ci_lo']:.6f}"
        )
        assert py_h["conf_int"][1] == pytest.approx(
            r_h["ci_hi"], rel=self.INFERENCE_RTOL
        ), (
            f"path={path_key} h={h} ci_hi: "
            f"py={py_h['conf_int'][1]:.6f} vs r={r_h['ci_hi']:.6f}"
        )


class TestDCDHDynRParityHeterogeneityWithPlacebo:
    """Parity tests for global ``predict_het`` + ``placebo`` (no by_path).

    Per codex R1 P1 #2: the Phase 1A change extended the global
    `_compute_heterogeneity_test` loop to cover backward horizons, so
    `results.heterogeneity_effects` now includes negative-int keys when
    ``placebo=True`` and ``heterogeneity=`` are co-set. This class pins
    the global surface independently of the per-path version
    (``TestDCDHDynRParityByPathHeterogeneityWithPlacebo``).

    Fixture: scenario 23 (``multi_path_reversible_predict_het_with_placebo_global``)
    — same DGP as scenarios 20-22 (reuses ``d20``) so per-horizon
    tolerances inherit from ``TestDCDHDynRParityHeterogeneity``.

    R syntax note: ``predict_het = list("X", c(-1))`` triggers
    "compute heterogeneity for ALL forward (1..effects) AND ALL placebo
    (1..placebo) positions" — same sentinel as scenario 22.
    """

    BETA_RTOL = 1e-6
    BETA_ATOL = 1e-9
    SE_RTOL = 1e-5
    INFERENCE_RTOL = 1e-4

    def test_parity_multi_path_reversible_predict_het_with_placebo_global(
        self, golden_values
    ):
        """Global heterogeneity on forward + backward horizons."""
        import warnings

        scenario = golden_values.get(
            "multi_path_reversible_predict_het_with_placebo_global"
        )
        if scenario is None:
            pytest.skip(
                "scenario 'multi_path_reversible_predict_het_with_placebo_global' "
                "not in golden values"
            )

        df = _golden_to_df_with_extra(scenario["data"], extra_cols=["het_x"])
        # Python's `placebo=True` triggers all-backward-horizons
        # (1..L_max) placebo computation; the parity check iterates only
        # R's emitted subset (placebo=2 in scenario 23 -> backward
        # horizons -1, -2). Python's extra backward horizons (e.g., -3
        # if eligible) are allowed and not parity-checked.
        est = ChaisemartinDHaultfoeuille(
            drop_larger_lower=False, placebo=True
        )
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", UserWarning)
            results = est.fit(
                df,
                outcome="outcome",
                group="group",
                time="period",
                treatment="treatment",
                L_max=3,
                heterogeneity="het_x",
            )

        assert results.heterogeneity_effects is not None
        py_het = results.heterogeneity_effects

        # Forward horizons (positive int keys).
        r_predict_het = scenario["results"]["predict_het"]
        for h_str, r_h in r_predict_het.items():
            h = int(h_str)
            assert h > 0
            assert h in py_het, f"forward horizon {h} missing"
            self._assert_horizon_parity(h, py_het[h], r_h)

        # Placebo horizons (negative int keys).
        r_placebo_het = _as_dict(scenario["results"].get("placebo_predict_het"))
        for h_str, r_h in r_placebo_het.items():
            h = int(h_str)
            assert h < 0
            assert h in py_het, (
                f"placebo horizon {h} missing from Python heterogeneity_effects"
            )
            self._assert_horizon_parity(h, py_het[h], r_h)

    def _assert_horizon_parity(self, h, py_h, r_h):
        """Pin all 6 inference fields against R."""
        assert py_h["beta"] == pytest.approx(
            r_h["beta"], rel=self.BETA_RTOL, abs=self.BETA_ATOL
        ), f"h={h} beta: py={py_h['beta']:.6f} vs r={r_h['beta']:.6f}"
        assert py_h["se"] == pytest.approx(r_h["se"], rel=self.SE_RTOL), (
            f"h={h} se: py={py_h['se']:.6f} vs r={r_h['se']:.6f}"
        )
        assert py_h["t_stat"] == pytest.approx(r_h["t"], rel=self.SE_RTOL), (
            f"h={h} t_stat: py={py_h['t_stat']:.6f} vs r={r_h['t']:.6f}"
        )
        assert int(py_h["n_obs"]) == int(r_h["n_obs"]), (
            f"h={h} n_obs: py={py_h['n_obs']} vs r={r_h['n_obs']}"
        )
        assert py_h["p_value"] == pytest.approx(
            r_h["p_value"], rel=self.INFERENCE_RTOL
        ), (
            f"h={h} p_value: py={py_h['p_value']:.6e} vs r={r_h['p_value']:.6e}"
        )
        assert py_h["conf_int"][0] == pytest.approx(
            r_h["ci_lo"], rel=self.INFERENCE_RTOL
        ), f"h={h} ci_lo: py={py_h['conf_int'][0]:.6f} vs r={r_h['ci_lo']:.6f}"
        assert py_h["conf_int"][1] == pytest.approx(
            r_h["ci_hi"], rel=self.INFERENCE_RTOL
        ), f"h={h} ci_hi: py={py_h['conf_int'][1]:.6f} vs r={r_h['ci_hi']:.6f}"
