hardware-counters

#!/usr/bin/env python3
"""
Performance prediction from CrayPat hardware counters (ARCHER2, AMD EPYC 7742).

--------------------------------------------------------------------------
Problem framing
--------------------------------------------------------------------------
CrayPat reports hardware counters for rank 0 only. Rank-0 cycle count is
strongly related to wall time (log-log r = 0.97), but it is NOT simply
cycles/peak-clock: the *effective* clock implied by cycles/runtime ranges from
1.98 GHz for GROMACS (compute-bound, near the 2.25 GHz peak) down to 0.23 GHz
for STREAM at 128 cores, where rank 0 spends most of its time stalled on
memory or waiting at MPI barriers.

That gap is the interesting part. We therefore decompose

        runtime  =  cycles / (f_peak * eta)

where eta in (0, 1] is a dimensionless *efficiency factor* absorbing stalls,
imbalance and communication. Cycles are measured; eta is what the model has
to learn, and it is exactly what the behavioural counters describe.

Predicting log10(eta) rather than log10(runtime) has two benefits:
  * it removes the trivially-predictable magnitude term, so reported skill
    reflects real learning rather than "more cycles = longer runtime";
  * the target is bounded and comparable across applications spanning 4 s to
    215 s, so no single slow application dominates the fit.

Three approaches are compared:
  analytic   runtime = cycles / peak clock. No fitting, no training data.
  direct     regress log10(runtime) on counters (the naive framing).
  residual   regress log10(eta) on counters, then reconstruct runtime.

--------------------------------------------------------------------------
Validation
--------------------------------------------------------------------------
Leave-one-application-out: train on four applications, predict the fifth.
A random split would place e.g. STREAM at 32 cores in train and STREAM at 64
in test -- near-duplicates -- and would badly overstate accuracy. LOAO
measures the quantity that matters: transfer to an unseen workload.

Outputs (in --outdir): metrics.csv, predictions.csv,
feature_importance.csv, summary.txt
"""
import argparse, warnings
from pathlib import Path

import numpy as np
import pandas as pd
from sklearn.linear_model import Ridge, Lasso
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.dummy import DummyRegressor
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer

warnings.filterwarnings("ignore")
RNG = 42
F_PEAK = 2.25e9          # AMD EPYC 7742 max clock (Hz)


def build_features(df):
    """Scale-free behavioural ratios describing *how* the machine is used."""
    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"])
    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)      # Roofline-style
    return X.replace([np.inf, -np.inf], np.nan)


def make_models():
    def pipe(est):
        return Pipeline([("impute", SimpleImputer(strategy="median")),
                         ("scale",  StandardScaler()),
                         ("est",    est)])
    return {
        "Mean baseline":     pipe(DummyRegressor(strategy="mean")),
        "Ridge":             pipe(Ridge(alpha=10.0)),
        "Lasso":             pipe(Lasso(alpha=0.05, max_iter=50000)),
        "Random Forest":     pipe(RandomForestRegressor(
                                 n_estimators=500, min_samples_leaf=2,
                                 random_state=RNG, n_jobs=-1)),
        "Gradient Boosting": pipe(GradientBoostingRegressor(
                                 n_estimators=300, max_depth=2,
                                 learning_rate=0.05, random_state=RNG)),
    }


def factors(pred, actual):
    """Multiplicative error: 1.0 is perfect, 2.0 means out by 2x either way."""
    pred = np.clip(pred, 1e-9, None)
    return np.maximum(pred / actual, actual / pred)


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--data",   default="/work/project/project/user/data/runs.csv")
    ap.add_argument("--outdir", default="/work/project/project/user/analysis/out")
    a = ap.parse_args()
    out = Path(a.outdir); out.mkdir(parents=True, exist_ok=True)

    df = pd.read_csv(a.data)
    df = df[df.app != "cp2k"]                       # no runtime recorded
    df = df[df.runtime_s.notna()].reset_index(drop=True)

    # Target: median wall time of counter sets B-E. Set A is always the first
    # run in a job and absorbs cold-cache / filesystem warm-up, so it is a
    # systematic outlier (HPL c1: 8.9 s for A vs ~4.6 s for B-E).
    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) if wc else df.runtime_s).values

    cyc      = pd.to_numeric(df.PAPI_TOT_CYC, errors="coerce").values
    t_analytic = cyc / F_PEAK                       # zero-parameter prediction
    eta        = t_analytic / y                     # efficiency factor in (0,1]
    X, groups  = build_features(df), df.app
    apps       = sorted(groups.unique())

    rows, preds = [], []

    # ---------------- analytic baseline (no training at all) --------------
    for app in apps:
        m = (groups == app).values
        f = factors(t_analytic[m], y[m])
        rows.append({"approach": "analytic", "model": "cycles / peak clock",
                     "held_out_app": app, "median_factor": float(np.median(f)),
                     "max_factor": float(f.max()),
                     "MAPE_pct": float(np.mean(np.abs(t_analytic[m]-y[m])/y[m])*100)})
        for nc, aa, pp in zip(df.ncore[m], y[m], t_analytic[m]):
            preds.append({"approach": "analytic", "model": "cycles / peak clock",
                          "app": app, "ncore": int(nc),
                          "actual_s": float(aa), "predicted_s": float(pp)})

    # ---------------- learned models --------------------------------------
    for approach in ("direct", "residual"):
        target = np.log10(y) if approach == "direct" else np.log10(eta)
        for name, model in make_models().items():
            for app in apps:
                te = (groups == app).values
                tr = ~te
                model.fit(X[tr], target[tr])
                p = model.predict(X[te])
                pred = 10 ** p if approach == "direct" else t_analytic[te] / (10 ** p)
                f = factors(pred, y[te])
                rows.append({"approach": approach, "model": name,
                             "held_out_app": app,
                             "median_factor": float(np.median(f)),
                             "max_factor": float(f.max()),
                             "MAPE_pct": float(np.mean(np.abs(pred-y[te])/y[te])*100)})
                for nc, aa, pp in zip(df.ncore[te], y[te], pred):
                    preds.append({"approach": approach, "model": name, "app": app,
                                  "ncore": int(nc), "actual_s": float(aa),
                                  "predicted_s": float(pp)})

    metrics = pd.DataFrame(rows)
    metrics.to_csv(out / "metrics.csv", index=False)
    pd.DataFrame(preds).to_csv(out / "predictions.csv", index=False)

    # ---------------- report ---------------------------------------------
    L = []
    L.append("Performance prediction from CrayPat hardware counters")
    L.append("=" * 78)
    L.append(f"{len(df)} configurations, {len(apps)} applications, "
             f"runtime {y.min():.1f}-{y.max():.1f} s")
    L.append("Validation: leave-one-application-out (train on 4, predict the 5th)")
    L.append("Error = multiplicative factor max(pred/act, act/pred); 1.00x is perfect")
    L.append("")
    L.append("Effective clock implied by cycles/runtime varies from 1.98 GHz")
    L.append("(GROMACS, compute-bound) to 0.23 GHz (STREAM at 128 cores, memory")
    L.append("stalled) against a 2.25 GHz peak. Learning that efficiency factor")
    L.append("is the actual prediction problem.")
    L.append("")

    agg = (metrics.groupby(["approach", "model"])
                  .agg(median_factor=("median_factor", "median"),
                       worst=("max_factor", "max"),
                       MAPE_pct=("MAPE_pct", "mean"))
                  .sort_values("median_factor"))
    L.append("All approaches ranked by typical error:")
    L.append(agg.round(2).to_string())
    L.append("")

    for ap_ in ("analytic", "direct", "residual"):
        sub = metrics[metrics.approach == ap_]
        piv = sub.pivot(index="model", columns="held_out_app", values="median_factor")
        L.append(f"--- {ap_}: median factor per held-out application ---")
        L.append(piv.round(2).to_string())
        L.append("")

    best = agg.index[0]
    L.append(f"Best: {best[1]} ({best[0]}) -- typical {agg.iloc[0]['median_factor']:.2f}x, "
             f"worst {agg.iloc[0]['worst']:.2f}x")
    bl = metrics[(metrics.model == "Mean baseline")].median_factor.median()
    L.append(f"Mean baseline (no counters, predict training mean): {bl:.2f}x")
    L.append("")

    rf = make_models()["Random Forest"]
    rf.fit(X, np.log10(eta))
    imp = (pd.Series(rf.named_steps["est"].feature_importances_, index=X.columns)
             .sort_values(ascending=False))
    imp.to_csv(out / "feature_importance.csv", header=["importance"])
    L.append("Which counters explain the efficiency factor (Random Forest):")
    for k, v in imp.head(10).items():
        L.append(f"  {k:22s} {v:.3f}")

    txt = "\n".join(L)
    (out / "summary.txt").write_text(txt + "\n")
    print(txt)


if __name__ == "__main__":
    main()