hardware-counters

#!/usr/bin/env python3
"""
Corrected training and evaluation, after an independent review found several
methodological defects in retrain.py. Each is fixed here and labelled.

F1  Broken baseline. The constant baseline used DummyRegressor(strategy="mean")
    fitted on log10(eta). The error metric max(p/a, a/p) is minimised by the
    MEDIAN, and log10(eta) is left-skewed (mean -0.169 vs median -0.099), so
    the "baseline" carried a systematic 15% over-prediction. It was a strawman
    and every comparison against it was inflated. Now uses the median.

F2  Pseudo-replication that also destroyed the features. Each configuration
    was run 5x with DIFFERENT counter sets, and those were treated as 5
    independent rows. But each set measures only ~5 counters, so a per-set row
    has most of its feature vector missing: measured availability is 0% for
    l2_miss_per_instr, l1_access_per_instr and prefetch_l2_frac, 50% for ipc,
    20% for arith_intensity. Merging to one row per configuration (union of
    counters across sets, median runtime) gives 191 rows with ~all features
    present instead of 815 rows that are ~58% imputed.

F3  Hyperparameters tuned on the evaluation folds. alpha and architecture were
    hardcoded after inspecting LOAO results, so the reported accuracy was
    optimistically biased. Now selected by an INNER leave-one-application-out
    loop over the training applications only.

F4  Misspecified target. eta was justified as a varying effective clock, but
    the data show a fixed additive overhead instead: median eta is 0.47 for
    runs under 0.3 s versus 0.83 above 3 s. The pure ratio cannot express
    that. An additive-overhead form is fitted and compared.

F5  Metric-matched loss. The evaluation metric is symmetric in log space, so
    squared error on log10(eta) is not matched to it. A quantile (median)
    objective and explicit shrinkage toward the training median are compared.

F6  Significance at the wrong granularity. Wilcoxon over 815 pseudo-rows
    treated correlated repeats as independent. Now computed per configuration.

F7  Broken LOSO folds. `size==0` for gromacs and openfoam meant gromacs was
    silently dropped from the test set while openfoam was tested with no
    training rows of its own. Sizes are now ranked within application.
"""
import os, warnings, sys, itertools
from pathlib import Path

import numpy as np
import pandas as pd
from sklearn.neural_network import MLPRegressor
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, FunctionTransformer
from sklearn.impute import SimpleImputer
from scipy.stats import wilcoxon

warnings.filterwarnings("ignore")
# Data and output roots. Overridable so the pipeline can be relocated
# (a different account, a scratch copy, another site) without editing code;
# the default keeps existing invocations working unchanged.
D      = os.environ.get("DISS_ROOT", "/work/project/project/user")
OUT    = Path(f"{D}/analysis/out"); OUT.mkdir(parents=True, exist_ok=True)
F_PEAK = 2.25e9
SEEDS  = (0, 1, 2)

COUNTERS = ["PAPI_TOT_CYC", "PAPI_TOT_INS", "PAPI_FP_OPS", "PAPI_FP_INS",
            "PAPI_L1_DCA", "PAPI_L2_DCH", "PAPI_L2_DCM", "PAPI_L2_DCR",
            "PAPI_TLB_DM", "PAPI_BR_INS", "PAPI_BR_MSP",
            "UNC_L3_CACHE_MISSES", "UNC_L3_MISS_LATENCY",
            "PACKAGE_ENERGY", "PP0_ENERGY",
            "L2_PREFETCH_HIT_L2", "L2_PREFETCH_HIT_L3",
            "REQUESTS_TO_L2_GROUP1:L2_HW_PF", "REQUESTS_TO_L2_GROUP1:RD_BLK_X",
            "DISPATCH_RESOURCE_STALL_CYCLES_1:LOAD_QUEUE_RSRC_STALL",
            "DISPATCH_RESOURCE_STALL_CYCLES_1:STORE_QUEUE_RSRC_STALL",
            "DISPATCH_RESOURCE_STALL_CYCLES_1:FP_REG_FILE_RSRC_STALL"]


def load_merged():
    """F2: one row per configuration, union of counters across counter sets."""
    df = pd.read_csv(f"{D}/data/runs_expanded.csv")
    df = df[df.runtime_s.notna() & (df.runtime_s > 0)]
    keys = ["app", "size", "ncore", "nthread", "nrank"]
    agg = {c: "median" for c in COUNTERS if c in df.columns}
    agg["runtime_s"] = "median"
    m = df.groupby(keys, as_index=False).agg(agg)
    m = m[m.PAPI_TOT_CYC.notna()].reset_index(drop=True)
    return m


def build_features(df):
    X = pd.DataFrame(index=df.index)
    g = lambda c: pd.to_numeric(df.get(c), errors="coerce")
    cyc, ins   = g("PAPI_TOT_CYC"), g("PAPI_TOT_INS")
    fpops      = g("PAPI_FP_OPS")
    l1a        = g("PAPI_L1_DCA")
    l2h, l2m   = g("PAPI_L2_DCH"), g("PAPI_L2_DCM")
    l3m, l3lat = g("UNC_L3_CACHE_MISSES"), g("UNC_L3_MISS_LATENCY")
    tlb        = g("PAPI_TLB_DM")
    pkg, pp0   = g("PACKAGE_ENERGY"), g("PP0_ENERGY")
    s_ld = g("DISPATCH_RESOURCE_STALL_CYCLES_1:LOAD_QUEUE_RSRC_STALL")
    s_st = g("DISPATCH_RESOURCE_STALL_CYCLES_1:STORE_QUEUE_RSRC_STALL")
    s_fp = g("DISPATCH_RESOURCE_STALL_CYCLES_1:FP_REG_FILE_RSRC_STALL")
    pf2, pf3 = g("L2_PREFETCH_HIT_L2"), g("L2_PREFETCH_HIT_L3")

    X["log_ncore"]   = np.log2(df.ncore.clip(lower=1))
    X["log_nrank"]   = np.log2(df.nrank.clip(lower=1))
    X["log_nthread"] = np.log2(df.nthread.clip(lower=1))
    X["ipc"]                 = ins / cyc
    X["flops_per_instr"]     = fpops / ins
    X["flops_per_cycle"]     = fpops / cyc
    X["l2_hit_rate"]         = l2h / (l2h + l2m)
    X["l2_miss_per_instr"]   = l2m / ins
    X["l1_access_per_instr"] = l1a / ins
    X["l3_miss_per_instr"]   = l3m / ins
    X["l3_lat_per_miss"]     = l3lat / l3m
    X["tlb_miss_per_instr"]  = tlb / ins
    X["stall_load_frac"]     = s_ld / cyc
    X["stall_store_frac"]    = s_st / cyc
    X["stall_fp_frac"]       = s_fp / cyc
    X["prefetch_l2_frac"]    = pf2 / (pf2 + pf3)
    X["energy_per_instr"]    = pkg / ins
    X["core_energy_frac"]    = pp0 / pkg
    X["arith_intensity"]     = fpops / (l3m * 64.0)
    X["log_instr_per_rank"]  = np.log10((ins / df.nrank.clip(lower=1)).clip(lower=1))
    # F4: analytic time is a strong covariate; short runs carry fixed overhead
    X["log_t_analytic"]      = np.log10((cyc / F_PEAK).clip(lower=1e-6))
    return X.replace([np.inf, -np.inf], np.nan)


def slog(A):
    return np.sign(A) * np.log1p(np.abs(A))


def mk(kind, seed=0, **kw):
    pre = [("impute", SimpleImputer(strategy="median")),
           ("log", FunctionTransformer(slog, validate=False)),
           ("scale", StandardScaler())]
    if kind == "mlp":
        est = MLPRegressor(hidden_layer_sizes=kw["h"], alpha=kw["a"],
                           solver="lbfgs", max_iter=5000, random_state=seed)
    elif kind == "rf":
        est = RandomForestRegressor(n_estimators=400, min_samples_leaf=2,
                                    random_state=seed, n_jobs=-1)
    elif kind == "gbq":            # F5: metric-matched median objective
        est = GradientBoostingRegressor(loss="quantile", alpha=0.5,
                                        n_estimators=300, max_depth=3,
                                        learning_rate=0.05, random_state=seed)
    return Pipeline(pre + [("est", est)])


def fac(p, a):
    p = np.clip(p, 1e-9, None)
    return np.maximum(p / a, a / p)


GRID = [{"h": h, "a": a}
        for h in [(16,), (32, 16), (64, 32)]
        for a in [0.1, 1.0, 3.0, 10.0]]


def fit_predict(kind, Xtr, ytr, Xte, seeds=SEEDS, **kw):
    ps = []
    for s in seeds:
        m = mk(kind, s, **kw); m.fit(Xtr, ytr); ps.append(m.predict(Xte))
    return np.mean(ps, axis=0)


def main():
    df = load_merged()
    X  = build_features(df)
    y  = df.runtime_s.values
    t_an = (df.PAPI_TOT_CYC / F_PEAK).values
    eta  = t_an / y
    apps = sorted(df.app.unique())

    L = []
    say = lambda s="": (print(s), L.append(s))
    say("Corrected evaluation")
    say("=" * 78)
    say(f"F2: merged {191 if len(df)==191 else len(df)} configurations "
        f"(was 815 pseudo-rows)")
    miss = X.isna().mean().mean() * 100
    say(f"    feature matrix now {100-miss:.1f}% populated "
        f"(was ~42% on per-set rows)")
    say(f"{len(df)} configs, {X.shape[1]} features, {len(apps)} applications")
    say("")

    rows = {}
    def record(name, pred_fn):
        r = np.full(len(y), np.nan)
        for a in apps:
            te = (df.app == a).values; tr = ~te
            r[te] = fac(pred_fn(X[tr], np.log10(eta[tr]), X[te], t_an[te]), y[te])
        rows[name] = r
        say(f"  {name:34s} median {np.nanmedian(r):.3f}   "
            f"p90 {np.nanpercentile(r,90):.3f}")

    say("--- LOAO (train on 7 apps, predict the 8th) ---")

    # F1: correct zero-parameter baseline
    record("Constant (median log-eta) [F1]",
           lambda Xtr, ytr, Xte, ta: ta / (10 ** np.full(len(Xte), np.median(ytr))))

    record("Random Forest",
           lambda Xtr, ytr, Xte, ta: ta / (10 ** fit_predict("rf", Xtr, ytr, Xte, (0,))))
    record("GBM quantile(0.5) [F5]",
           lambda Xtr, ytr, Xte, ta: ta / (10 ** fit_predict("gbq", Xtr, ytr, Xte, (0,))))
    record("MLP (32,16) a=1 [tuned on test]",
           lambda Xtr, ytr, Xte, ta: ta / (10 ** fit_predict("mlp", Xtr, ytr, Xte, h=(32,16), a=1.0)))

    # F3: nested selection - inner LOAO over training apps only
    def nested(Xtr, ytr, Xte, ta, tr_apps):
        best, bg = None, None
        for gcfg in GRID:
            errs = []
            for ia in tr_apps:
                ite = (tr_apps == ia)
                if ite.sum() == 0 or (~ite).sum() == 0:
                    continue
                p = fit_predict("mlp", Xtr[~ite], ytr[~ite], Xtr[ite], (0,), **gcfg)
                errs.append(np.median(np.abs(p - ytr[ite])))
            e = np.mean(errs) if errs else 1e9
            if best is None or e < best:
                best, bg = e, gcfg
        return ta / (10 ** fit_predict("mlp", Xtr, ytr, Xte, **bg)), bg

    r = np.full(len(y), np.nan); picks = []
    for a in apps:
        te = (df.app == a).values; tr = ~te
        p, bg = nested(X[tr], np.log10(eta[tr]), X[te], t_an[te], df.app[tr].values)
        r[te] = fac(p, y[te]); picks.append(f"{a}:{bg['h']}a{bg['a']}")
    rows["MLP nested selection [F3]"] = r
    say(f"  {'MLP nested selection [F3]':34s} median {np.nanmedian(r):.3f}   "
        f"p90 {np.nanpercentile(r,90):.3f}")
    say(f"      inner loop picked: {', '.join(picks)}")

    # F5: shrinkage toward the training median
    def shrunk(Xtr, ytr, Xte, ta, w=0.5):
        p = fit_predict("mlp", Xtr, ytr, Xte, h=(32, 16), a=1.0)
        return ta / (10 ** (w * p + (1 - w) * np.median(ytr)))
    record("MLP shrunk 0.5 to median [F5]", shrunk)

    # ensemble
    def ens(Xtr, ytr, Xte, ta):
        a1 = fit_predict("mlp", Xtr, ytr, Xte, h=(32, 16), a=1.0)
        a2 = fit_predict("rf", Xtr, ytr, Xte, (0,))
        return ta / (10 ** (0.5 * a1 + 0.5 * a2))
    record("MLP + RF ensemble", ens)

    # F6: significance at configuration granularity
    say("")
    say("--- F6: paired Wilcoxon at CONFIG granularity (n=%d) ---" % len(y))
    base = rows["Constant (median log-eta) [F1]"]
    for k, v in rows.items():
        if k.startswith("Constant"):
            continue
        m = ~np.isnan(v) & ~np.isnan(base)
        try:
            _, p = wilcoxon(v[m], base[m])
        except ValueError:
            p = float("nan")
        better = np.nanmedian(v) < np.nanmedian(base)
        say(f"  {k:34s} vs constant: wins {int((v[m]<base[m]).sum()):3d}/{int(m.sum())}"
            f"  p={p:.4g}  ({'better' if better else 'WORSE'})")

    say("")
    say("--- per-application median (LOAO) ---")
    tb = pd.DataFrame(rows); tb["app"] = df.app.values
    say(tb.groupby("app").median().round(3).to_string())

    (OUT / "retrain2_summary.txt").write_text("\n".join(L) + "\n")
    tb.to_csv(OUT / "retrain2_perrow.csv", index=False)


if __name__ == "__main__":
    main()