hardware-counters

#!/usr/bin/env python3
"""
Cross-platform performance prediction: ARCHER2 (Zen 2) vs Cirrus (Zen 5).

This is the experiment the dissertation title promises. Everything before it
was single-platform, so "Modelling Hardware" was never actually tested.

THE FEATURE PROBLEM
-------------------
The two machines do not expose the same counters. ARCHER2 has cray_rapl
(PACKAGE_ENERGY, PP0_ENERGY) and cray_zenl3 (UNC_L3_*); Cirrus has neither.
Cirrus has presets ARCHER2 lacks (L1_DCM, L2_TCM, VEC_INS, FMA_INS). Any model
that must run on both is therefore restricted to the INTERSECTION, which
excludes the energy features. Those were once reported as the strongest
ARCHER2 predictors (core_energy_frac at 0.51 importance); that figure is
RETRACTED, see docs/modelling.md. Quantifying the cost of the intersection is
itself
a result.

EXPERIMENT MATRIX
-----------------
Six settings, each evaluated with the corrected protocol (merged
configurations, median-constant baseline, error = max(pred/act, act/pred)):

  1  within-A2       train and test on ARCHER2, leave-one-application-out
  2  within-CIR      train and test on Cirrus, leave-one-application-out
  3  A2 -> CIR       train on ALL of ARCHER2, predict ALL of Cirrus
  4  CIR -> A2       the reverse
  5  A2 -> CIR LOAO  train on ARCHER2 minus app X, predict Cirrus app X
                     (unseen hardware AND unseen application simultaneously)
  6  CIR -> A2 LOAO  the reverse

Settings 3-6 are the real test. 1-2 are references showing what each platform
costs on its own.

Two feature sets are compared throughout:
  full          every feature available on the training platform
  intersection  only features computable on BOTH platforms

and a `platform` indicator is tested as an extra feature in the pooled models,
to see whether the model can exploit knowing which machine it is on.
"""
import os, warnings, sys, itertools
from pathlib import Path

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.linear_model import Ridge
from sklearn.neural_network import MLPRegressor
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, FunctionTransformer
from sklearn.pipeline import Pipeline
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)

# peak clock differs per platform - this matters for the analytic term
F_PEAK = {"archer2": 2.25e9, "cirrus": 3.71e9}   # verified on compute nodes:
#   ARCHER2 EPYC 7742 max 2.25 GHz; Cirrus EPYC 9825 max 3.7146 GHz (lscpu).
#   The login node reports 2.4 GHz for a DIFFERENT part (EPYC 9745) - do not use it.

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",
            "PAPI_L1_DCM", "PAPI_L2_TCM", "PAPI_VEC_INS", "PAPI_FMA_INS",
            "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",
            "DISPATCH_RESOURCE_STALL_CYCLES_1:FP_SCHEDULER_RSRC_STALL"]

# features computable on BOTH platforms
# NOTE: branch_miss_rate was originally listed here, but PAPI_BR_INS/BR_MSP
# were never actually collected on ARCHER2 (counter set E recorded them but the
# merge dropped them), so it is 0% available there and cannot be in the shared
# set. Verified with xplat_diag.py.
INTERSECTION = ["log_ncore", "log_nrank", "log_nthread", "log_t_analytic",
                "log_instr_per_rank", "ipc", "flops_per_instr",
                "flops_per_cycle", "l2_hit_rate", "l2_miss_per_instr",
                "l1_access_per_instr", "tlb_miss_per_instr",
                "stall_load_frac", "stall_store_frac"]


def load(path, platform):
    df = pd.read_csv(path)
    if "platform" not in df.columns:
        df["platform"] = platform
    df = df[df.runtime_s.notna() & (df.runtime_s > 0)]
    keys = ["platform", "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)
    return m[m.PAPI_TOT_CYC.notna()].reset_index(drop=True)


def features(df):
    X = pd.DataFrame(index=df.index)
    g = lambda c: pd.to_numeric(df[c], errors="coerce") if c in df.columns \
                  else pd.Series(np.nan, index=df.index)
    cyc, ins = g("PAPI_TOT_CYC"), g("PAPI_TOT_INS")
    fp   = g("PAPI_FP_OPS")
    l1a  = g("PAPI_L1_DCA")
    l2h, l2m = g("PAPI_L2_DCH"), g("PAPI_L2_DCM")
    tlb  = g("PAPI_TLB_DM")
    br, brm = g("PAPI_BR_INS"), g("PAPI_BR_MSP")
    s_ld = g("DISPATCH_RESOURCE_STALL_CYCLES_1:LOAD_QUEUE_RSRC_STALL")
    s_st = g("DISPATCH_RESOURCE_STALL_CYCLES_1:STORE_QUEUE_RSRC_STALL")
    fpk  = df.platform.map(F_PEAK).values

    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["log_t_analytic"]      = np.log10((cyc / fpk).clip(lower=1e-6))
    X["log_instr_per_rank"]  = np.log10((ins / df.nrank.clip(lower=1)).clip(lower=1))
    X["ipc"]                 = ins / cyc
    X["flops_per_instr"]     = fp / ins
    X["flops_per_cycle"]     = fp / cyc
    X["l2_hit_rate"]         = l2h / (l2h + l2m)
    X["l2_miss_per_instr"]   = l2m / ins
    X["l1_access_per_instr"] = l1a / ins
    X["tlb_miss_per_instr"]  = tlb / ins
    X["branch_miss_rate"]    = brm / br
    X["stall_load_frac"]     = s_ld / cyc
    X["stall_store_frac"]    = s_st / cyc
    # ---- platform-specific extras (NaN where unavailable) ----
    X["l3_miss_per_instr"]   = g("UNC_L3_CACHE_MISSES") / ins
    X["l3_lat_per_miss"]     = g("UNC_L3_MISS_LATENCY") / g("UNC_L3_CACHE_MISSES")
    X["energy_per_instr"]    = g("PACKAGE_ENERGY") / ins
    X["core_energy_frac"]    = g("PP0_ENERGY") / g("PACKAGE_ENERGY")
    X["l1_miss_per_instr"]   = g("PAPI_L1_DCM") / ins
    X["vec_per_instr"]       = g("PAPI_VEC_INS") / ins
    X["fma_per_instr"]       = g("PAPI_FMA_INS") / ins
    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):
    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=-1),
           "gbq": GradientBoostingRegressor(loss="quantile", alpha=0.5,
                                            n_estimators=300, max_depth=3,
                                            learning_rate=0.05, random_state=seed),
           "ridge": Ridge(alpha=10.0),
           "mlp": MLPRegressor(hidden_layer_sizes=(32, 16), alpha=1.0,
                               solver="lbfgs", max_iter=5000, random_state=seed),
           }[kind]
    return Pipeline(pre + [("e", est)])


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


def run(Xtr, ytr_eta, Xte, t_an_te, y_te, kind, seeds=(0,)):
    ps = []
    for s in seeds:
        m = mk(kind, s); m.fit(Xtr, np.log10(ytr_eta))
        ps.append(m.predict(Xte))
    return fac(t_an_te / (10 ** np.mean(ps, axis=0)), y_te)


def main():
    a2  = load(f"{D}/data/runs_expanded.csv", "archer2")
    cir = load(f"{D}/data/cirrus_persets.csv", "cirrus")
    both = pd.concat([a2, cir], ignore_index=True)

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

    L = []
    say = lambda s="": (print(s, flush=True), L.append(s))
    say("Cross-platform performance prediction: ARCHER2 (Zen 2) vs Cirrus (Zen 5)")
    say("=" * 78)
    say(f"ARCHER2 {isA2.sum()} configs   Cirrus {(~isA2).sum()} configs   "
        f"total {len(both)}")
    say(f"apps: {sorted(both.app.unique())}")
    say("")
    say("feature availability (% non-null):")
    for c in X.columns:
        pa = 100 * X.loc[isA2, c].notna().mean()
        pc = 100 * X.loc[~isA2, c].notna().mean()
        mark = "  <-- intersection" if c in INTERSECTION else ""
        say(f"  {c:22s} ARCHER2 {pa:5.1f}%   Cirrus {pc:5.1f}%{mark}")
    say("")

    FEATSETS = {"intersection": [c for c in INTERSECTION if c in X.columns],
                "full":         list(X.columns)}
    rows = []

    def report(setting, featset, kind, err, n):
        med, p90 = np.nanmedian(err), np.nanpercentile(err, 90)
        rows.append({"setting": setting, "features": featset, "model": kind,
                     "n_test": n, "median": med, "p90": p90})
        say(f"  {setting:16s} {featset:13s} {kind:6s} n={n:4d}  "
            f"median {med:.4f}  p90 {p90:.4f}")

    for fs, cols in FEATSETS.items():
        Xf = X[cols].values
        say(f"--- feature set: {fs} ({len(cols)} features) ---")

        # 1/2: within-platform LOAO
        for plat, mask in [("within-A2", isA2), ("within-CIR", ~isA2)]:
            for kind in ("rf", "gbq", "mlp"):
                e = np.full(len(y), np.nan)
                sub = np.where(mask)[0]
                for a in sorted(both.app[mask].unique()):
                    te = sub[(both.app.values[sub] == a)]
                    tr = sub[(both.app.values[sub] != a)]
                    if len(te) == 0 or len(tr) < 10:
                        continue
                    e[te] = run(Xf[tr], eta[tr], Xf[te], t_an[te], y[te], kind)
                report(plat, fs, kind, e[mask], int(mask.sum()))
            # constant reference
            e = np.full(len(y), np.nan); sub = np.where(mask)[0]
            for a in sorted(both.app[mask].unique()):
                te = sub[both.app.values[sub] == a]; tr = sub[both.app.values[sub] != a]
                if len(te) == 0 or len(tr) < 10: continue
                e[te] = fac(t_an[te] / (10 ** np.median(np.log10(eta[tr]))), y[te])
            report(plat, fs, "const", e[mask], int(mask.sum()))

        # 3/4: whole-platform transfer
        for name, tr_m, te_m in [("A2->CIR", isA2, ~isA2), ("CIR->A2", ~isA2, isA2)]:
            for kind in ("rf", "gbq", "mlp", "ridge"):
                e = run(Xf[tr_m], eta[tr_m], Xf[te_m], t_an[te_m], y[te_m], kind)
                report(name, fs, kind, e, int(te_m.sum()))
            e = fac(t_an[te_m] / (10 ** np.median(np.log10(eta[tr_m]))), y[te_m])
            report(name, fs, "const", e, int(te_m.sum()))

        # 5/6: transfer AND unseen application at once
        for name, tr_m, te_m in [("A2->CIR LOAO", isA2, ~isA2),
                                 ("CIR->A2 LOAO", ~isA2, isA2)]:
            for kind in ("rf", "mlp"):
                e = np.full(len(y), np.nan)
                for a in sorted(both.app.unique()):
                    te = np.where(te_m & (both.app == a).values)[0]
                    tr = np.where(tr_m & (both.app != a).values)[0]
                    if len(te) == 0 or len(tr) < 10:
                        continue
                    e[te] = run(Xf[tr], eta[tr], Xf[te], t_an[te], y[te], kind)
                report(name, fs, kind, e[te_m], int(np.isfinite(e[te_m]).sum()))
        say("")

    res = pd.DataFrame(rows)
    res.to_csv(OUT / "crossplatform_results.csv", index=False)

    say("=== summary: best model per setting/featureset ===")
    best = (res[res.model != "const"].sort_values("median")
               .groupby(["setting", "features"]).first().reset_index())
    con  = res[res.model == "const"].set_index(["setting", "features"])["median"]
    for _, r in best.sort_values(["setting", "features"]).iterrows():
        c = con.get((r.setting, r.features), np.nan)
        say(f"  {r.setting:16s} {r.features:13s} best={r.model:5s} "
            f"{r['median']:.4f}   const={c:.4f}   "
            f"{'BEATS' if r['median'] < c else 'loses to'} constant")

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


if __name__ == "__main__":
    main()