hardware-counters

#!/usr/bin/env python3
"""
Retraction analysis: why `core_energy_frac` looked like the strongest predictor.

The project previously reported `core_energy_frac` (the PP0/PACKAGE energy
ratio) as the single strongest predictor of the efficiency factor eta, at 0.509
random forest impurity importance, and described it as the clearest novelty
claim. This script establishes where that number came from, decomposes the
three separate defects that produced it, and reports the corrected ranking.

PROVENANCE. The 0.509 figure is in results/model/feature_importance.csv, which
was written by scripts/analysis/train_model.py. That script read data/runs.csv,
which is NOT the 815-row per-counter-set table: it is 40 rows (5 applications
by 8 core counts, one problem size each) built by an earlier merge that joined
counter sets on (app, ncore). It carried 17 features and no analytic time term.
Establishing this matters, because the correction has to name the right defect.

THREE DEFECTS, each demonstrated separately below.

  A  Missing scale covariates. train_model.py had no `log_t_analytic` and no
     `log_instr_per_rank`. eta is strongly driven by run length (short runs
     carry fixed MPI_Init/IO/instrumentation overhead), so with no term for
     run length the forest had to reach for whatever correlated with it.
     Adding just those two columns to the SAME 40 rows drops core_energy_frac
     from 0.526 to 0.148.

  B  One problem size per application, at n = 40. In runs.csv each application
     appears at exactly one problem size, so application identity, problem size
     and counter signature are perfectly confounded, and the sample is small
     enough for a spurious association to dominate a greedy split. The raw
     correlation between core_energy_frac and log10(eta) is 0.44 on those 40
     rows but only 0.088 on the SAME FIVE applications once several problem
     sizes per application are present.

  C  Structural missingness. On the 815-row per-counter-set table the energy
     pair sat on separate registers and was recorded in every set, while every
     other counter appeared in only one or two sets. Median imputation
     collapses a mostly-missing column to a near-constant, a constant cannot
     reduce impurity, and the importance mass therefore migrates to whatever
     happens to be observed. This script shows impurity importance on that
     table correlates with per-feature OBSERVATION FREQUENCY rather than with
     signal, and confirms it with a control that re-imposes the same
     missingness pattern on the corrected data.

Robustness: permutation importance on held-out application folds is computed as
a model-agnostic alternative to impurity importance, and the two rankings are
compared.

Outputs: out/importance_artefact.txt, out/fig_importance_artefact.png
"""
import os, sys, warnings
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.stats import spearmanr, pearsonr
from sklearn.ensemble import RandomForestRegressor
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.inspection import permutation_importance
from sklearn.model_selection import GroupKFold

warnings.filterwarnings("ignore")

D = os.environ.get("DISS_ROOT", "/work/project/project/user")
sys.path.insert(0, f"{D}/analysis")
sys.path.insert(0, f"{D}/repo/user/dissertation/scripts/analysis")
from retrain2 import build_features, load_merged, F_PEAK  # noqa: E402
import train_model as legacy                              # noqa: E402

OUT = Path(f"{D}/analysis/out"); OUT.mkdir(parents=True, exist_ok=True)
SEED = 0
SCALE_TERMS = ["log_t_analytic", "log_instr_per_rank", "log_nrank", "log_nthread"]

L = []
def say(s=""):
    print(s, flush=True); L.append(s)


def median_impute(X):
    """Median imputation that KEEPS all-NaN columns.

    sklearn's SimpleImputer silently drops a column that is entirely missing.
    Three features are 100% missing on the per-set table, and dropping them
    would hide exactly the effect under study, so an all-NaN column is filled
    with zero (a constant, which is the degenerate case of median imputation)
    and retained, keeping the importance vectors comparable across tables.
    """
    Xi = X.copy()
    for c in Xi.columns:
        med = Xi[c].median()
        Xi[c] = Xi[c].fillna(0.0 if not np.isfinite(med) else med)
    return Xi


def rf_importance(X, y, seed=SEED, n=500):
    rf = RandomForestRegressor(n_estimators=n, min_samples_leaf=2,
                               random_state=seed, n_jobs=-1)
    rf.fit(median_impute(X).values, y)
    return pd.Series(rf.feature_importances_, index=X.columns)


def rankof(s, f):
    return int(s.rank(ascending=False)[f])


def eta_sq(v, groups):
    """Share of a feature's variance that lies BETWEEN applications.

    A feature with eta^2 near 1 is close to an application label. On a table
    with one problem size per application such a feature cannot be separated
    from application identity, which is defect B.
    """
    m = np.isfinite(v)
    v = np.asarray(v)[m]; g = np.asarray(groups)[m]
    if len(v) < 3:
        return np.nan
    s = pd.Series(v); gg = pd.Series(g)
    gm, n = s.groupby(gg).mean(), s.groupby(gg).size()
    ssb = (n * (gm - s.mean()) ** 2).sum(); sst = ((s - s.mean()) ** 2).sum()
    return ssb / sst if sst > 0 else np.nan


def load_legacy40():
    """data/runs.csv exactly as train_model.py consumed it: the source of 0.509."""
    df = pd.read_csv(f"{D}/data/runs.csv")
    df = df[df.app != "cp2k"]
    df = df[df.runtime_s.notna()].reset_index(drop=True)
    wc = [c for c in ["wall_B", "wall_C", "wall_D", "wall_E"] if c in df.columns]
    y = (df[wc].median(axis=1).fillna(df.runtime_s)).values
    t_an = pd.to_numeric(df.PAPI_TOT_CYC, errors="coerce").values / F_PEAK
    return df, legacy.build_features(df), np.log10(t_an / y), t_an


def load_persets():
    """The 815-row per-counter-set table, rows that can form a target."""
    df = pd.read_csv(f"{D}/data/runs_expanded.csv")
    df = df[df.runtime_s.notna() & (df.runtime_s > 0)]
    return df[df.PAPI_TOT_CYC.notna()].reset_index(drop=True)


def target(df):
    t_an = (pd.to_numeric(df.PAPI_TOT_CYC, errors="coerce") / F_PEAK).values
    return np.log10(t_an / df.runtime_s.values)


def main():
    say("Retraction analysis: the core_energy_frac importance claim")
    say("=" * 78)
    say("")
    say("Claim under test: `core_energy_frac` (PP0/PACKAGE) is the strongest")
    say("predictor of the efficiency factor eta, at 0.509 RF importance, and is")
    say("the project's clearest novelty claim.")
    say("")

    # ------------------------------------------------------------ 0. provenance
    d40, X40, y40, ta40 = load_legacy40()
    i40 = rf_importance(X40, y40)
    say("--- 0. Provenance of the 0.509 figure ---")
    say(f"results/model/feature_importance.csv was written by train_model.py from")
    say(f"data/runs.csv: {len(d40)} rows, {X40.shape[1]} features, "
        f"{d40.app.nunique()} applications, one problem size each.")
    say(f"Reproduced here: core_energy_frac = {i40['core_energy_frac']:.3f} "
        f"(file records 0.509), rank {rankof(i40,'core_energy_frac')}.")
    say("It was NOT computed on the 815-row per-counter-set table. That matters,")
    say("because it means three defects contributed, not one.")
    say("")

    # -------------------------------------------- data for the corrected view
    ps = load_persets(); yp = target(ps)
    mg = load_merged();  ym = target(mg)
    Xm = build_features(mg)
    feats = list(Xm.columns)
    Xp = build_features(ps)[feats]

    obs_p = 1.0 - Xp.isna().mean()
    obs_m = 1.0 - Xm.isna().mean()

    ip = rf_importance(Xp, yp)
    im = rf_importance(Xm, ym)

    # ------------------------------------------------ A. missing scale terms
    say("--- A. Defect one: no analytic time term ---")
    X40s = X40.copy()
    X40s["log_t_analytic"] = np.log10(np.clip(ta40, 1e-6, None))
    X40s["log_instr_per_rank"] = np.log10(
        np.clip(pd.to_numeric(d40.PAPI_TOT_INS, errors="coerce"), 1, None))
    i40s = rf_importance(X40s, y40)
    say("eta is dominated by run length: short runs carry a fixed MPI_Init, I/O")
    say("and instrumentation overhead. train_model.py gave the forest no term")
    say("for run length, so the split had to be made on a correlate of it.")
    say(f"  same 40 rows, 17 features            core_energy_frac "
        f"{i40['core_energy_frac']:.3f}  rank {rankof(i40,'core_energy_frac')}")
    say(f"  + log_t_analytic, log_instr_per_rank core_energy_frac "
        f"{i40s['core_energy_frac']:.3f}  rank {rankof(i40s,'core_energy_frac')}")
    say(f"  the added term takes the top slot: log_t_analytic "
        f"{i40s['log_t_analytic']:.3f}")
    say(f"  adding two columns and changing nothing else removes "
        f"{(1-i40s['core_energy_frac']/i40['core_energy_frac'])*100:.0f}% of the claimed importance.")
    say("")

    # ------------------------------- B. one problem size per app, and small n
    say("--- B. Defect two: one problem size per application, n = 40 ---")
    five = sorted(d40.app.unique())
    sub = mg.app.isin(five).values
    r40 = pd.Series(X40.core_energy_frac.values).corr(pd.Series(y40))
    r5 = pd.Series(Xm.core_energy_frac.values[sub]).corr(pd.Series(ym[sub]))
    rall = pd.Series(Xm.core_energy_frac.values).corr(pd.Series(ym))
    say("In runs.csv every application appears at exactly one problem size, so")
    say("application identity, problem size and counter signature cannot be")
    say("separated. The association does not survive when sizes vary.")
    say(f"  corr(core_energy_frac, log10 eta), 40 rows, one size per app : {r40:+.3f}")
    say(f"  corr, SAME five applications, merged, several sizes each      : {r5:+.3f}")
    say(f"  corr, all eight applications, merged                          : {rall:+.3f}")
    e40 = eta_sq(X40.core_energy_frac.values, d40.app.values)
    em = eta_sq(Xm.core_energy_frac.values, mg.app.values)
    et40 = eta_sq(y40, d40.app.values)
    say(f"  between-application variance share of core_energy_frac: "
        f"{e40:.3f} on the 40 rows, {em:.3f} merged")
    say(f"  between-application variance share of the target itself: {et40:.3f} (40 rows)")
    say("  With five applications, one size each and eight rows apiece, a feature")
    say("  that half tracks application identity can split the target almost as")
    say("  well as a label would, and impurity importance rewards exactly that.")
    i5 = rf_importance(Xm[sub].drop(columns=SCALE_TERMS), ym[sub])
    say(f"  merged, same five applications, same 17 features: core_energy_frac "
        f"{i5['core_energy_frac']:.3f}, rank {rankof(i5,'core_energy_frac')} "
        f"(n = {int(sub.sum())})")
    acc = []
    for s in range(30):
        idx = np.random.default_rng(s).choice(len(mg), 40, replace=False)
        acc.append(rf_importance(Xm.iloc[idx].drop(columns=SCALE_TERMS),
                                 ym[idx], seed=s, n=300))
    isub = pd.concat(acc, axis=1).mean(axis=1)
    say(f"  merged, subsampled to n=40 (30 draws), 17 features: core_energy_frac "
        f"{isub['core_energy_frac']:.3f}, rank {rankof(isub,'core_energy_frac')}")
    say("  Small n alone inflates it a little; small n plus a single problem size")
    say("  per application is what inflates it to the top of the table.")
    say("")

    # ------------------------------------------- C. structural missingness
    say("--- C. Defect three: structural missingness on the 815-row table ---")
    say(f"  per-counter-set table : {len(ps):4d} rows, "
        f"{obs_p.mean()*100:.1f}% of feature cells observed")
    say(f"  merged configurations : {len(mg):4d} rows, "
        f"{obs_m.mean()*100:.1f}% of feature cells observed")
    say("")
    say(f"{'feature':24s} {'obs815':>7s} {'imp815':>8s} {'rk':>3s}   "
        f"{'obs191':>7s} {'imp191':>8s} {'rk':>3s}")
    for f in ip.sort_values(ascending=False).index:
        say(f"{f:24s} {obs_p[f]*100:6.1f}% {ip[f]:8.3f} {rankof(ip,f):3d}   "
            f"{obs_m[f]*100:6.1f}% {im[f]:8.3f} {rankof(im,f):3d}")
    say("")
    say("If impurity importance measured signal it would be unrelated to how")
    say("often a feature happens to be recorded. It is not:")
    for label, imp_, obs in (("815 pseudo-rows", ip, obs_p),
                             ("191 merged configs", im, obs_m)):
        rs, ps_ = spearmanr(obs.values, imp_.values)
        rr, pr_ = pearsonr(obs.values, imp_.values)
        say(f"  {label:22s} Spearman rho = {rs:+.3f} (p = {ps_:.3g}), "
            f"Pearson r = {rr:+.3f} (p = {pr_:.3g})")
    say("  The association is present on the structurally-missing table and")
    say("  absent on the merged one, which is the signature of the artefact.")
    fully = obs_p[obs_p > 0.99].index.tolist()
    say(f"  fully observed on the 815-row table ({len(fully)}): {', '.join(fully)}")
    say(f"    summed importance there {ip[fully].sum():.3f} "
        f"({ip[fully].sum()*100:.0f}% of the total, from "
        f"{len(fully)}/{len(feats)} features)")
    say(f"    the same features on the merged table {im[fully].sum():.3f} "
        f"({im[fully].sum()*100:.0f}%)")
    zero = obs_p[obs_p < 0.001].index.tolist()
    say(f"  never observed there ({len(zero)}): {', '.join(zero)}")
    say(f"    importance there {ip[zero].sum():.3f}, on the merged table "
        f"{im[zero].sum():.3f}, so real signal was being suppressed as well as")
    say("    misattributed.")
    say("")
    say("  Control: re-impose the 815-row missingness on the merged data. Same")
    say("  191 configurations, same signal, only the observation pattern copied")
    say("  across, averaged over 20 random masks.")
    accm = []
    for rep in range(20):
        rng = np.random.default_rng(rep)
        Xmask = Xm.copy()
        for f in feats:
            keep = obs_p[f]
            if keep < 0.999:
                Xmask.loc[rng.random(len(Xmask)) > keep, f] = np.nan
        accm.append(rf_importance(Xmask, ym, seed=rep, n=300))
    imask = pd.concat(accm, axis=1).mean(axis=1)
    say(f"    masked-merged core_energy_frac {imask['core_energy_frac']:.3f}, "
        f"rank {rankof(imask,'core_energy_frac')} "
        f"(unmasked {im['core_energy_frac']:.3f}, rank {rankof(im,'core_energy_frac')})")
    say(f"    masked-merged importance vs 815-row observation frequency: "
        f"Spearman rho = {spearmanr(obs_p.values, imask.values)[0]:+.3f} "
        f"(p = {spearmanr(obs_p.values, imask.values)[1]:.3g})")
    say(f"    masked-merged ranking vs the 815-row ranking:   rho = "
        f"{spearmanr(ip.values, imask.values)[0]:+.3f}")
    say(f"    unmasked-merged ranking vs the 815-row ranking: rho = "
        f"{spearmanr(ip.values, im.values)[0]:+.3f}")
    say("  The mask alone recreates most of the old ranking from data known to")
    say("  carry the corrected signal, so the pattern, not the physics, produced it.")
    say("")

    # -------------------------------------------- the corrected ranking
    say("--- The corrected result ---")
    say(f"  core_energy_frac, 40-row runs.csv (retracted) : "
        f"{i40['core_energy_frac']:.3f}, rank {rankof(i40,'core_energy_frac')} of {X40.shape[1]}")
    say(f"  core_energy_frac, 815 per-set rows            : "
        f"{ip['core_energy_frac']:.3f}, rank {rankof(ip,'core_energy_frac')} of {len(feats)}")
    say(f"  core_energy_frac, 191 merged configurations   : "
        f"{im['core_energy_frac']:.3f}, rank {rankof(im,'core_energy_frac')} of {len(feats)}")
    say("  corrected top five (merged, impurity):")
    for f, v in im.sort_values(ascending=False).head(5).items():
        say(f"    {f:24s} {v:.3f}")
    say("")

    # -------------------------------- robustness: permutation importance
    say("--- Robustness: permutation importance on held-out folds ---")
    say("Impurity importance is measured on training data and is biased towards")
    say("high-cardinality features. Permutation importance is measured on held-")
    say("out folds grouped by application, matching the LOAO protocol, so it")
    say("cannot reward a feature the model merely memorised.")
    pi_acc = []
    gkf = GroupKFold(n_splits=min(5, mg.app.nunique()))
    for tr, te in gkf.split(Xm, ym, groups=mg.app.values):
        pipe = Pipeline([("impute", SimpleImputer(strategy="median")),
                         ("rf", RandomForestRegressor(n_estimators=400,
                                                      min_samples_leaf=2,
                                                      random_state=SEED,
                                                      n_jobs=-1))])
        pipe.fit(Xm.iloc[tr], ym[tr])
        r = permutation_importance(pipe, Xm.iloc[te], ym[te], n_repeats=20,
                                   random_state=SEED,
                                   scoring="neg_mean_absolute_error")
        pi_acc.append(pd.Series(r.importances_mean, index=feats))
    pi = pd.concat(pi_acc, axis=1).mean(axis=1)

    say("")
    say(f"{'feature':24s} {'perm':>9s} {'rk':>3s}  {'impurity':>9s} {'rk':>3s}")
    for f in pi.sort_values(ascending=False).index:
        say(f"{f:24s} {pi[f]:9.4f} {rankof(pi,f):3d}  {im[f]:9.4f} {rankof(im,f):3d}")
    rho_agree, p_agree = spearmanr(pi.values, im.values)
    say("")
    say(f"  permutation vs impurity ranking agreement: Spearman rho = "
        f"{rho_agree:+.3f} (p = {p_agree:.3g})")
    say("  The two measures agree on direction but not in detail, which is the")
    say("  expected outcome and a reason not to over-read any single ranking.")
    say(f"  Both place core_energy_frac in the bottom half: permutation "
        f"{pi['core_energy_frac']:.4f}, rank {rankof(pi,'core_energy_frac')} of {len(feats)}.")
    say(f"  permutation top four: "
        f"{', '.join(pi.sort_values(ascending=False).head(4).index)}")
    say(f"  impurity top four   : "
        f"{', '.join(im.sort_values(ascending=False).head(4).index)}")
    say("  Both agree that log_t_analytic leads and that the remaining signal is")
    say("  in instruction mix and cache behaviour, not in energy. The features")
    say("  they disagree about (stall_fp_frac, log_instr_per_rank) are the ones")
    say("  where impurity importance is known to be optimistic.")
    say("")

    # ---------------------------------------------------------- conclusion
    say("--- Conclusion ---")
    say("The claim is retracted. `core_energy_frac` is not the strongest")
    say(f"predictor of eta; on the corrected 191-configuration table it is worth")
    say(f"{im['core_energy_frac']:.3f} (rank {rankof(im,'core_energy_frac')} of {len(feats)}) by impurity and rank "
        f"{rankof(pi,'core_energy_frac')} by held-out permutation")
    say("importance. The leading terms are the analytic time scale, the")
    say("instruction mix and cache behaviour.")
    say("")
    say("The 0.509 was produced by three compounding defects, none of them about")
    say("energy: a missing run-length covariate (A), a design in which each")
    say("application appeared at a single problem size at n = 40 (B), and, on the")
    say("later per-counter-set table, structural missingness that made the energy")
    say("pair the only fully observed feature (C). Each is demonstrated above by")
    say("changing one thing at a time.")
    say("")
    say("What replaces the claim is a methodological result, and it generalises")
    say("beyond this study. Impurity importance is not interpretable on a feature")
    say("matrix assembled by counter multiplexing. Rotating counter sets is the")
    say("standard way to widen a feature space past a small hardware register")
    say("budget, and it produces structurally missing columns by construction;")
    say("imputation then makes those columns uninformative, and the importance")
    say("mass has to go somewhere. Any counter that sits on a separate register")
    say("and is therefore recorded in every set, RAPL energy being the usual")
    say("case, will absorb it and look like a discovery. The safeguards are to")
    say("merge to one row per configuration before interpreting importance, to")
    say("include an explicit scale covariate so run length is not attributed to a")
    say("correlate, to vary problem size within application so importance cannot")
    say("proxy application identity, and to confirm with held-out permutation")
    say("importance rather than training-set impurity.")

    (OUT / "importance_artefact.txt").write_text("\n".join(L) + "\n")

    # ------------------------------------------------------------- figure
    fig, ax = plt.subplots(2, 2, figsize=(15, 11))

    a = ax[0, 0]
    steps = ["runs.csv\n40 rows, 17 feats\n(as published)",
             "+ scale covariates\n(defect A)",
             "merged 191\n17 feats\n(defects A+B)",
             "merged 191\n21 feats\n(corrected)"]
    vals = [i40["core_energy_frac"], i40s["core_energy_frac"],
            rf_importance(Xm.drop(columns=SCALE_TERMS), ym)["core_energy_frac"],
            im["core_energy_frac"]]
    a.bar(range(4), vals, color=["#c44", "#d86", "#8ac", "#48c"])
    for i, v in enumerate(vals):
        a.text(i, v + 0.012, f"{v:.3f}", ha="center", fontsize=9)
    a.set_xticks(range(4)); a.set_xticklabels(steps, fontsize=7)
    a.set_ylabel("core_energy_frac impurity importance")
    a.set_title("(a) The retracted claim, defect by defect")
    a.set_ylim(0, max(vals) * 1.2)

    a = ax[0, 1]
    a.scatter(obs_p.values * 100, ip.values, c="#c44", s=34, label="815 per-set rows")
    a.scatter(obs_m.values * 100, im.values, c="#48c", s=34, marker="^",
              label="191 merged configs")
    for f in ip.sort_values(ascending=False).head(5).index:
        a.annotate(f, (obs_p[f] * 100, ip[f]), fontsize=7,
                   xytext=(3, 3), textcoords="offset points")
    rs_p = spearmanr(obs_p.values, ip.values)[0]
    rs_m = spearmanr(obs_m.values, im.values)[0]
    a.set_xlabel("percentage of rows in which the feature is observed")
    a.set_ylabel("impurity importance")
    a.set_title(f"(b) Defect C: importance tracks observability\n"
                f"per-set rho = {rs_p:+.2f}, merged rho = {rs_m:+.2f}")
    a.legend(fontsize=8)

    a = ax[1, 0]
    order = ip.sort_values(ascending=False).index[:12]
    yy = np.arange(len(order))
    a.barh(yy - 0.27, ip[order].values, 0.27, label="815 per-set rows", color="#c44")
    a.barh(yy, imask[order].values, 0.27,
           label="191 merged, per-set mask re-imposed", color="#e93")
    a.barh(yy + 0.27, im[order].values, 0.27, label="191 merged, unmasked", color="#48c")
    a.set_yticks(yy); a.set_yticklabels(order, fontsize=8); a.invert_yaxis()
    a.set_xlabel("random forest impurity importance")
    a.set_title("(c) Control: the missingness pattern alone recreates the ranking")
    a.legend(fontsize=8)

    a = ax[1, 1]
    o2 = im.sort_values(ascending=False).index[:12]
    y2 = np.arange(len(o2))
    a.barh(y2 - 0.2, im[o2].values / max(im.max(), 1e-12), 0.4,
           label="impurity, training set", color="#48c")
    a.barh(y2 + 0.2, np.clip(pi[o2].values, 0, None) / max(pi.max(), 1e-12), 0.4,
           label="permutation, held-out folds", color="#3a7")
    a.set_yticks(y2); a.set_yticklabels(o2, fontsize=8); a.invert_yaxis()
    a.set_xlabel("importance, normalised to the maximum of each measure")
    a.set_title(f"(d) Corrected ranking, two measures, rho = {rho_agree:+.2f}")
    a.legend(fontsize=8)

    fig.suptitle("Feature importance under structural missingness: "
                 "retraction of the core_energy_frac claim", fontsize=13)
    fig.tight_layout(rect=[0, 0, 1, 0.97])
    fig.savefig(OUT / "fig_importance_artefact.png", dpi=150)
    print(f"\nwrote {OUT}/importance_artefact.txt and fig_importance_artefact.png")


if __name__ == "__main__":
    main()