hardware-counters

#!/usr/bin/env python3
"""
A single, correct significance-reporting routine for the whole project.

WHY THIS FILE EXISTS
--------------------
Three defects were found in how the earlier scripts reported comparisons.

R1  Direction reported by the wrong quantity.  Every script ran a *paired*
    two-sided Wilcoxon signed-rank test on (model, baseline) error factors and
    then labelled the direction with `np.nanmedian(v) < np.nanmedian(base)`,
    i.e. by comparing two *independent* medians.  The signed-rank test is about
    the distribution of the paired DIFFERENCES, and median(a) - median(b) is
    not median(a - b).  When the two disagree you get nonsense lines such as

        GBM quantile(0.5)   wins 118/191  p=0.0004311  (WORSE)
        MLP nested selection wins  94/191  p=0.3079    (better)

    Here the direction is taken from the sign of the median of the paired
    differences, which is the quantity the test is actually about.

R2  No multiple-comparison control.  Roughly twenty model-versus-baseline
    comparisons were reported at a nominal alpha = 0.05 with no correction, so
    the family-wise error rate was of order 1 - 0.95**20 = 64%.  Holm-Bonferroni
    (step-down, uniformly more powerful than plain Bonferroni and valid under
    arbitrary dependence) is applied across each declared family.

R3  No uncertainty and no effect size.  A median error factor of 1.0548 was
    quoted to four decimals with no interval, and "p = 0.0002" was reported for
    a median improvement of 0.011x that is of no engineering consequence.  Here
    every arm carries a bootstrap 95% CI on its median (>= 10000 resamples),
    and every comparison carries the median paired difference with its own
    bootstrap CI plus the matched-pairs rank-biserial correlation.

USAGE
-----
    from stats_rigour import Family
    fam = Family("ARCHER2 LOAO")
    fam.add("Random Forest", err_rf, err_const)
    fam.add("GBM quantile", err_gbq, err_const)
    fam.finalise()          # applies Holm within the family
    print(fam.table())
"""
from __future__ import annotations

import numpy as np
import pandas as pd
from scipy.stats import wilcoxon

N_BOOT = 10000
ALPHA = 0.05


# --------------------------------------------------------------- primitives
def boot_median_ci(x, n_boot=N_BOOT, alpha=ALPHA, seed=0):
    """Percentile bootstrap CI for the median of x."""
    x = np.asarray(x, float)
    x = x[np.isfinite(x)]
    if len(x) < 3:
        return (np.nan, np.nan)
    rng = np.random.default_rng(seed)
    idx = rng.integers(0, len(x), size=(n_boot, len(x)))
    meds = np.median(x[idx], axis=1)
    lo, hi = np.percentile(meds, [100 * alpha / 2, 100 * (1 - alpha / 2)])
    return float(lo), float(hi)


def boot_paired_median_diff_ci(d, n_boot=N_BOOT, alpha=ALPHA, seed=0):
    """Percentile bootstrap CI for the median of the paired differences.

    Resampling is over PAIRS, which is what keeps the pairing intact.
    """
    d = np.asarray(d, float)
    d = d[np.isfinite(d)]
    if len(d) < 3:
        return (np.nan, np.nan)
    rng = np.random.default_rng(seed)
    idx = rng.integers(0, len(d), size=(n_boot, len(d)))
    meds = np.median(d[idx], axis=1)
    lo, hi = np.percentile(meds, [100 * alpha / 2, 100 * (1 - alpha / 2)])
    return float(lo), float(hi)


def hodges_lehmann(d):
    """Hodges-Lehmann pseudomedian: median of all Walsh averages (d_i+d_j)/2,
    i <= j.

    This is THE point estimate associated with the Wilcoxon signed-rank test,
    so its sign is the direction the test is actually testing.  It is not in
    general equal to the median of d, and it is certainly not equal to
    median(a) - median(b).
    """
    d = np.asarray(d, float)
    d = d[np.isfinite(d)]
    n = len(d)
    if n == 0:
        return np.nan
    i, j = np.triu_indices(n)
    return float(np.median((d[i] + d[j]) / 2.0))


def boot_hl_ci(d, n_boot=2000, alpha=ALPHA, seed=0):
    """Bootstrap CI for the Hodges-Lehmann estimate (O(n^2) per resample, so
    fewer resamples; 2000 is ample for a 95% interval)."""
    d = np.asarray(d, float)
    d = d[np.isfinite(d)]
    if len(d) < 3:
        return (np.nan, np.nan)
    rng = np.random.default_rng(seed)
    n = len(d)
    i, j = np.triu_indices(n)
    vals = np.empty(n_boot)
    for b in range(n_boot):
        s = d[rng.integers(0, n, n)]
        vals[b] = np.median((s[i] + s[j]) / 2.0)
    lo, hi = np.percentile(vals, [100 * alpha / 2, 100 * (1 - alpha / 2)])
    return float(lo), float(hi)


def rank_biserial(d):
    """Matched-pairs rank-biserial correlation, the effect size that belongs
    with the Wilcoxon signed-rank test.

    r = (W+ - W-) / (W+ + W-) computed on the signed ranks of |d|, so it runs
    from -1 (every pair favours the model) to +1 (every pair favours the
    baseline) given d = model - baseline.  |r| >= 0.1 small, 0.3 medium,
    0.5 large by the usual convention.
    """
    d = np.asarray(d, float)
    d = d[np.isfinite(d) & (d != 0)]
    if len(d) == 0:
        return np.nan
    from scipy.stats import rankdata
    r = rankdata(np.abs(d))
    wp, wn = r[d > 0].sum(), r[d < 0].sum()
    tot = wp + wn
    return float((wp - wn) / tot) if tot > 0 else np.nan


def holm(pvals):
    """Holm-Bonferroni step-down adjusted p-values.

    Returns adjusted p-values in the ORIGINAL order.  NaN inputs are passed
    through as NaN and excluded from the family size.
    """
    p = np.asarray(pvals, float)
    out = np.full(len(p), np.nan)
    ok = np.where(np.isfinite(p))[0]
    m = len(ok)
    if m == 0:
        return out
    order = ok[np.argsort(p[ok])]
    running = 0.0
    for i, j in enumerate(order):
        adj = (m - i) * p[j]
        running = max(running, adj)          # enforce monotonicity
        out[j] = min(1.0, running)
    return out


# --------------------------------------------------------------- comparison
def compare(name, err_model, err_base, seed=0, n_boot=N_BOOT):
    """One paired comparison, reported correctly.

    err_model / err_base are per-configuration error factors, aligned
    element-wise (same configuration in the same slot).  NaNs are dropped
    pairwise.
    """
    a = np.asarray(err_model, float)
    b = np.asarray(err_base, float)
    m = np.isfinite(a) & np.isfinite(b)
    a, b = a[m], b[m]
    n = len(a)
    if n < 3:
        return {"comparison": name, "n": n, "median_model": np.nan,
                "median_base": np.nan, "median_paired_diff": np.nan,
                "p_raw": np.nan, "direction": "not evaluable"}

    d = a - b                                  # negative => model better
    md = float(np.median(d))
    try:
        p = float(wilcoxon(a, b, zero_method="wilcox")[1])
    except ValueError:                          # all differences zero
        p = 1.0

    lo_m, hi_m = boot_median_ci(a, n_boot, seed=seed)
    lo_b, hi_b = boot_median_ci(b, n_boot, seed=seed + 1)
    lo_d, hi_d = boot_paired_median_diff_ci(d, n_boot, seed=seed + 2)
    hl = hodges_lehmann(d)
    lo_h, hi_h = boot_hl_ci(d, seed=seed + 3)

    # R1: direction from the paired test's own quantity.  The Hodges-Lehmann
    # pseudomedian is the location estimate the signed-rank test targets, so
    # it is the primary; the plain median of the differences is reported too
    # and the two agree in every case here.
    ref = hl if np.isfinite(hl) else md
    if ref < 0:
        direction = "model better"
    elif ref > 0:
        direction = "model worse"
    else:
        direction = "tied"

    med_a, med_b = float(np.median(a)), float(np.median(b))
    naive = "model better" if med_a < med_b else \
            ("model worse" if med_a > med_b else "tied")

    return {
        "comparison": name, "n": n,
        "median_model": med_a, "model_ci_lo": lo_m, "model_ci_hi": hi_m,
        "median_base": med_b, "base_ci_lo": lo_b, "base_ci_hi": hi_b,
        "median_paired_diff": md, "diff_ci_lo": lo_d, "diff_ci_hi": hi_d,
        "hodges_lehmann": hl, "hl_ci_lo": lo_h, "hl_ci_hi": hi_h,
        "wins": int((d < 0).sum()), "losses": int((d > 0).sum()),
        "rank_biserial": rank_biserial(d),
        "p_raw": p,
        "direction": direction,
        "naive_direction": naive,
        "direction_conflict": direction != naive,
    }


class Family:
    """A family of comparisons sharing one Holm-Bonferroni correction."""

    def __init__(self, label, alpha=ALPHA, n_boot=N_BOOT, seed=0):
        self.label, self.alpha, self.n_boot, self.seed = label, alpha, n_boot, seed
        self.rows = []
        self._final = False

    def add(self, name, err_model, err_base):
        self.rows.append(compare(name, err_model, err_base,
                                 seed=self.seed + 7 * len(self.rows),
                                 n_boot=self.n_boot))
        return self

    def add_not_evaluable(self, name, reason):
        self.rows.append({"comparison": name, "n": 0, "p_raw": np.nan,
                          "direction": "not evaluable", "reason": reason,
                          "median_model": np.nan, "median_base": np.nan,
                          "median_paired_diff": np.nan})
        return self

    def finalise(self):
        adj = holm([r.get("p_raw", np.nan) for r in self.rows])
        for r, pa in zip(self.rows, adj):
            r["family"] = self.label
            r["p_holm"] = pa
            r["significant_raw"] = bool(np.isfinite(r.get("p_raw", np.nan))
                                        and r["p_raw"] < self.alpha)
            r["significant_holm"] = bool(np.isfinite(pa) and pa < self.alpha)
            r["survives_holm"] = r["significant_raw"] and r["significant_holm"]
            r["lost_to_holm"] = r["significant_raw"] and not r["significant_holm"]
        self._final = True
        return self

    def frame(self):
        if not self._final:
            self.finalise()
        return pd.DataFrame(self.rows)

    def table(self, width=46):
        """Human-readable block for the text report."""
        df = self.frame()
        out = [f"### {self.label}   (family size {int(np.isfinite(df.p_raw).sum())}"
               f", Holm-Bonferroni, alpha={self.alpha})"]
        hdr = (f"  {'comparison':{width}s} {'n':>4s} {'median [95% CI]':>16s} "
               f"{'baseline [95%CI]':>16s} {'HodgesLehmann diff':>19s} "
               f"{'w/l':>8s} {'r_rb':>6s} "
               f"{'p_raw':>10s} {'p_holm':>10s}  verdict")
        out.append(hdr)
        out.append("  " + "-" * (len(hdr) - 2))
        for _, r in df.iterrows():
            if r["direction"] == "not evaluable":
                out.append(f"  {r['comparison']:{width}s} {'--':>4s} "
                           f"{'NOT EVALUABLE':>16s}   {r.get('reason','')}")
                continue
            med = f"{r.median_model:.4f} [{r.model_ci_lo:.3f},{r.model_ci_hi:.3f}]"
            bas = f"{r.median_base:.4f} [{r.base_ci_lo:.3f},{r.base_ci_hi:.3f}]"
            dif = f"{r.hodges_lehmann:+.4f} [{r.hl_ci_lo:+.3f},{r.hl_ci_hi:+.3f}]"
            wl = f"{int(r.wins)}/{int(r.losses)}"
            if r.survives_holm:
                verdict = "BETTER (Holm)" if r.direction == "model better" \
                          else "WORSE (Holm)"
            elif r.lost_to_holm:
                verdict = f"n.s. after Holm ({r.direction}, raw p<{self.alpha})"
            else:
                verdict = f"n.s. ({r.direction})"
            flag = "  [!direction conflict with naive median comparison]" \
                if r.get("direction_conflict") else ""
            out.append(f"  {r['comparison']:{width}s} {int(r.n):4d} {med:>16s} "
                       f"{bas:>16s} {dif:>19s} {wl:>8s} {r.rank_biserial:+6.3f} "
                       f"{r.p_raw:10.3g} {r.p_holm:10.3g}  {verdict}{flag}")
        out.append("")
        return "\n".join(out)


# ============================================================ project re-run
if __name__ == "__main__":
    import sys, warnings
    from pathlib import Path
    warnings.filterwarnings("ignore")

    D = "/work/project/project/user"
    OUT = Path(f"{D}/analysis/out"); OUT.mkdir(parents=True, exist_ok=True)
    sys.path.insert(0, f"{D}/analysis")

    from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
    from sklearn.neural_network import MLPRegressor
    from sklearn.linear_model import Ridge
    from sklearn.impute import SimpleImputer
    from sklearn.preprocessing import StandardScaler, FunctionTransformer
    from sklearn.pipeline import Pipeline
    from crossplatform import (load as xload, features as xfeatures,
                               F_PEAK, INTERSECTION, slog, fac)

    NJ = 8   # shared login node: do not grab all 512 cores

    def mk(kind, seed=0):
        pre = [("i", SimpleImputer(strategy="median")),
               ("l", FunctionTransformer(slog, validate=False)),
               ("s", StandardScaler())]
        est = {"rf": RandomForestRegressor(n_estimators=400, min_samples_leaf=2,
                                           random_state=seed, n_jobs=NJ),
               "gbq": GradientBoostingRegressor(loss="quantile", alpha=0.5,
                                                n_estimators=300, max_depth=3,
                                                learning_rate=0.05,
                                                random_state=seed),
               "mlp": MLPRegressor(hidden_layer_sizes=(32, 16), alpha=1.0,
                                   solver="lbfgs", max_iter=5000,
                                   random_state=seed),
               "ridge": Ridge(alpha=10.0)}[kind]
        return Pipeline(pre + [("e", est)])

    def prep(path, plat, tag):
        d = xload(path, plat)
        d["dataset"] = tag
        return d

    a2 = prep(f"{D}/data/runs_expanded.csv", "archer2", "archer2")
    cf = prep(f"{D}/data/cirrus_persets.csv", "cirrus", "cirrus-fixed")
    cw = prep(f"{D}/data/cirrus_weak.csv", "cirrus", "cirrus-weak")
    allf = pd.concat([a2, cf, cw], ignore_index=True)

    X = xfeatures(allf)
    y = allf.runtime_s.values
    fpk = allf.platform.map(F_PEAK).values
    t_an = pd.to_numeric(allf.PAPI_TOT_CYC, errors="coerce").values / fpk
    eta = t_an / y
    ok = (eta > 0) & (eta <= 1.5)
    allf = allf[ok].reset_index(drop=True)
    X = X[ok].reset_index(drop=True)
    y, t_an, eta = y[ok], t_an[ok], eta[ok]

    INTER = [c for c in INTERSECTION if c in X.columns]

    def loao(mask, cols, kind):
        """LOAO within one dataset; returns per-row error, NaN elsewhere."""
        e = np.full(len(y), np.nan)
        Xf = X[cols].values
        sub = np.where(mask)[0]
        for a in sorted(allf.app[mask].unique()):
            te = sub[allf.app.values[sub] == a]
            tr = sub[allf.app.values[sub] != a]
            if len(te) == 0 or len(tr) < 10:
                continue
            if kind == "const":
                p = t_an[te] / (10 ** np.median(np.log10(eta[tr])))
            else:
                m = mk(kind); m.fit(Xf[tr], np.log10(eta[tr]))
                p = t_an[te] / (10 ** m.predict(Xf[te]))
            e[te] = fac(p, y[te])
        return e

    def transfer(tr_mask, te_mask, cols, kind):
        e = np.full(len(y), np.nan)
        Xf = X[cols].values
        tr, te = np.where(tr_mask)[0], np.where(te_mask)[0]
        if kind == "const":
            p = t_an[te] / (10 ** np.median(np.log10(eta[tr])))
        else:
            m = mk(kind); m.fit(Xf[tr], np.log10(eta[tr]))
            p = t_an[te] / (10 ** m.predict(Xf[te]))
        e[te] = fac(p, y[te])
        return e

    L = []
    say = lambda s="": (print(s, flush=True), L.append(s))
    say("Corrected statistical reporting for all key comparisons")
    say("=" * 100)
    say("Direction is the SIGN OF THE HODGES-LEHMANN PSEUDOMEDIAN of the paired")
    say("differences (model - baseline), the location estimate the Wilcoxon")
    say("signed-rank test actually targets; NOT a comparison of two independent")
    say(f"medians.  Intervals are percentile bootstrap 95% CIs ({N_BOOT} resamples")
    say("for arm medians, 2000 for the O(n^2) Hodges-Lehmann).  p_holm is Holm-")
    say("Bonferroni adjusted within the family shown in each block heading.")
    say("r_rb is the matched-pairs rank-biserial effect size; w/l counts pairs.")
    say("Error factor = max(pred/actual, actual/pred); 1.0 is perfect.")
    say("")

    fams = []

    # ---- family 1-3: models vs constant, one family per dataset -----------
    dsets = [("archer2", "ARCHER2 (fixed size, LOAO)"),
             ("cirrus-fixed", "Cirrus fixed-size (LOAO)"),
             ("cirrus-weak", "Cirrus weak-scaled (LOAO)")]
    for tag, label in dsets:
        mask = (allf.dataset == tag).values
        base = loao(mask, INTER, "const")
        fam = Family(f"{label} vs constant baseline, intersection features")
        for kind, nm in [("rf", "Random Forest"), ("gbq", "GBM quantile(0.5)"),
                         ("mlp", "MLP (32,16) alpha=1")]:
            fam.add(nm, loao(mask, INTER, kind), base)
        fams.append(fam.finalise())
        say(fam.table())

    # ---- family 4: counter-group ablation on ARCHER2 ----------------------
    GROUPS = {
        "config": ["log_ncore", "log_nrank", "log_nthread"],
        "analytic": ["log_t_analytic"],
        "cheap": ["ipc", "log_instr_per_rank", "flops_per_instr",
                  "flops_per_cycle"],
    }
    CUM = [("config only (free)", ["config"]),
           ("+ analytic cycles (1 counter)", ["config", "analytic"]),
           ("+ instruction ratios (1 set)", ["config", "analytic", "cheap"]),
           ("+ all counters (5 sets)", None)]
    m_a2 = (allf.dataset == "archer2").values
    abl = {}
    for lab, grps in CUM:
        cols = list(X.columns) if grps is None else \
            [c for g in grps for c in GROUPS[g] if c in X.columns]
        abl[lab] = loao(m_a2, cols, "rf")
    const_a2 = loao(m_a2, INTER, "const")

    fam = Family("Counter-group ablation on ARCHER2, RF, vs CONFIG-ONLY")
    for lab in [k for k, _ in CUM][1:]:
        fam.add(lab, abl[lab], abl["config only (free)"])
    fams.append(fam.finalise())
    say(fam.table())

    fam = Family("Counter-group ablation on ARCHER2, RF, vs CONSTANT baseline")
    for lab, _ in CUM:
        fam.add(lab, abl[lab], const_a2)
    fams.append(fam.finalise())
    say(fam.table())

    # ---- family 5: cross-platform transfer -------------------------------
    m_cf = (allf.dataset == "cirrus-fixed").values
    m_cw = (allf.dataset == "cirrus-weak").values
    fam = Family("Cross-platform transfer vs constant, intersection features")
    for src, dst, sm, dm in [("ARCHER2", "Cirrus-fixed", m_a2, m_cf),
                             ("Cirrus-fixed", "ARCHER2", m_cf, m_a2),
                             ("ARCHER2", "Cirrus-weak", m_a2, m_cw),
                             ("Cirrus-weak", "ARCHER2", m_cw, m_a2)]:
        b = transfer(sm, dm, INTER, "const")
        for kind, nm in [("rf", "RF"), ("gbq", "GBQ"), ("mlp", "MLP")]:
            fam.add(f"{src} -> {dst}: {nm}", transfer(sm, dm, INTER, kind), b)
    fams.append(fam.finalise())
    say(fam.table())

    # ---- global view ------------------------------------------------------
    big = pd.concat([f.frame() for f in fams], ignore_index=True)
    big.to_csv(OUT / "stats_corrected.csv", index=False)

    say("=" * 100)
    say("WHAT THE CORRECTION COSTS")
    say("")
    lost = big[big.get("lost_to_holm", False) == True]
    say(f"  comparisons run                        : {len(big)}")
    say(f"  significant at raw alpha=0.05          : {int(big.significant_raw.sum())}")
    say(f"  still significant after Holm-Bonferroni: {int(big.survives_holm.sum())}")
    say(f"  LOST to the correction                 : {len(lost)}")
    for _, r in lost.iterrows():
        say(f"      {r['family']}  |  {r['comparison']}  "
            f"raw p={r.p_raw:.4g} -> Holm p={r.p_holm:.4g}")
    say("")
    conf = big[big.get("direction_conflict", False) == True]
    say(f"  comparisons where the OLD naive direction label (median(a) vs")
    say(f"  median(b)) disagrees with the paired direction: {len(conf)}")
    for _, r in conf.iterrows():
        say(f"      {r['family']}  |  {r['comparison']}: "
            f"paired HL={r.hodges_lehmann:+.4f} says '{r['direction']}', "
            f"naive median said '{r['naive_direction']}'")
    say("")
    say("  Comparisons that survive Holm AND favour the model:")
    win = big[(big.get("survives_holm", False) == True) &
              (big.direction == "model better")]
    if len(win) == 0:
        say("      none")
    for _, r in win.iterrows():
        say(f"      {r['family']}  |  {r['comparison']}: "
            f"{r.median_model:.4f} vs {r.median_base:.4f}, "
            f"HL diff {r.hodges_lehmann:+.4f} "
            f"[{r.hl_ci_lo:+.4f},{r.hl_ci_hi:+.4f}], "
            f"r_rb={r.rank_biserial:+.3f}, Holm p={r.p_holm:.3g}")

    (OUT / "stats_corrected.txt").write_text("\n".join(L) + "\n")
    print(f"\nwrote {OUT/'stats_corrected.txt'} and .csv")