"""Statistical testing (spec §10.1): bootstrap 95% CIs at the pair level,
paired comparisons across seeds (Wilcoxon signed-rank + Cohen's d), Holm
correction across RQ families.
"""
from __future__ import annotations

import numpy as np
from scipy import stats as scipy_stats
from typing import Callable, TypeVar

Pair = TypeVar("Pair")


def bootstrap_metric_ci(
    pairs: list[Pair], metric_fn: Callable[[list[Pair]], float],
    n_resamples: int = 10_000, ci: float = 0.95, seed: int = 42,
) -> tuple[float, float, float]:
    """Pair-level CI for any metric, including non-additive F1/MCC/PR-AUC.

    The same sampled pair indices must be used for both halves of each pair.
    """
    if not pairs:
        raise ValueError("cannot bootstrap an empty pair set")
    if n_resamples <= 0 or not 0 < ci < 1:
        raise ValueError("n_resamples and ci must be positive and valid")
    rng = np.random.default_rng(seed)
    estimate = float(metric_fn(pairs))
    n = len(pairs)
    values = np.empty(n_resamples)
    for i in range(n_resamples):
        indices = rng.integers(0, n, size=n)
        values[i] = float(metric_fn([pairs[j] for j in indices]))
    alpha = (1 - ci) / 2
    return estimate, float(np.quantile(values, alpha)), float(np.quantile(values, 1 - alpha))


def mcnemar_label_correct(a: list[bool], b: list[bool]) -> dict:
    """Exact two-sided McNemar test for paired label correctness (OOD, n=134).

    Each position must refer to the same pair in both model outputs. Exact
    binomial inference remains valid when discordant counts are small.
    """
    if len(a) != len(b) or not a:
        raise ValueError("paired correctness vectors must have equal nonzero length")
    a_only = sum(bool(x) and not bool(y) for x, y in zip(a, b))
    b_only = sum(not bool(x) and bool(y) for x, y in zip(a, b))
    both = sum(bool(x) and bool(y) for x, y in zip(a, b))
    neither = len(a) - a_only - b_only - both
    discordant = a_only + b_only
    p_value = float(scipy_stats.binomtest(min(a_only, b_only), discordant, 0.5).pvalue) if discordant else 1.0
    return {"n_pairs": len(a), "both_correct": both, "a_only_correct": a_only,
            "b_only_correct": b_only, "neither_correct": neither,
            "p_value": p_value, "method": "exact_two_sided_mcnemar"}


def bootstrap_ci(
    values: list[float], n_resamples: int = 10_000, ci: float = 0.95, seed: int = 42
) -> tuple[float, float, float]:
    """Returns (point_estimate, ci_low, ci_high) via the percentile bootstrap
    over the pair-level metric values (spec: "bootstrap 95% CIs (10k
    resamples) at the pair level")."""
    values = np.asarray(values, dtype=float)
    if len(values) == 0:
        return (0.0, 0.0, 0.0)
    rng = np.random.default_rng(seed)
    point = float(np.mean(values))
    n = len(values)
    resample_means = np.empty(n_resamples)
    for i in range(n_resamples):
        sample = rng.choice(values, size=n, replace=True)
        resample_means[i] = sample.mean()
    alpha = (1 - ci) / 2
    lo = float(np.quantile(resample_means, alpha))
    hi = float(np.quantile(resample_means, 1 - alpha))
    return (point, lo, hi)


def cohens_d_paired(a: list[float], b: list[float]) -> float:
    """Cohen's d for paired samples (a - b), using the pooled SD of the
    differences — the standard formulation for a paired design."""
    a, b = np.asarray(a, dtype=float), np.asarray(b, dtype=float)
    diffs = a - b
    sd = diffs.std(ddof=1)
    if sd == 0:
        return 0.0
    return float(diffs.mean() / sd)


def wilcoxon_signed_rank(a: list[float], b: list[float]) -> dict:
    """Paired comparison across seeds (spec §10.1). Returns the statistic,
    p-value, and Cohen's d. Falls back gracefully when all paired differences
    are zero (scipy raises in that degenerate case)."""
    a, b = np.asarray(a, dtype=float), np.asarray(b, dtype=float)
    diffs = a - b
    if np.allclose(diffs, 0):
        return {"statistic": 0.0, "p_value": 1.0, "cohens_d": 0.0}
    stat, p = scipy_stats.wilcoxon(a, b)
    return {"statistic": float(stat), "p_value": float(p), "cohens_d": cohens_d_paired(a, b)}


def holm_correction(p_values: dict[str, float], alpha: float = 0.05) -> dict[str, dict]:
    """Holm-Bonferroni step-down correction across a named family of
    hypotheses (spec: "Holm correction across RQ families"). Returns, per
    name, the adjusted alpha threshold used at its rank and whether it's
    significant under the family-wise correction."""
    ranked = sorted(p_values.items(), key=lambda kv: kv[1])
    m = len(ranked)
    out: dict[str, dict] = {}
    reject_all_remaining = False
    for i, (name, p) in enumerate(ranked):
        threshold = alpha / (m - i)
        significant = (not reject_all_remaining) and (p <= threshold)
        if not significant:
            reject_all_remaining = True  # Holm: once one fails, all later (larger p) fail too
        out[name] = {"p_value": p, "threshold": threshold, "significant": significant, "rank": i + 1}
    return out