hardware-counters

#!/usr/bin/env python3
"""
Baseline ladder: how much of the prediction is done by each successive layer of
information, under one identical protocol.

The motivation is a scoping problem rather than a modelling one. The pipeline
predicts an efficiency factor eta and reconstructs runtime as
t_analytic / eta, with t_analytic = PAPI_TOT_CYC / f_peak. The cycle count is
measured on the configuration being "predicted", so t_analytic is not a
prediction at all, it is a measurement. Any headline accuracy figure is
therefore dominated by information that was handed to the model for free. This
script quantifies exactly how much, by evaluating five rungs of increasing
measurement cost under the same leave-one-application-out folds, the same
merged 191-configuration dataset and the same multiplicative error metric.

Rungs:
  0  median training runtime            no measurement of the target at all
  1  configuration only                 free metadata, still no measurement
  2  t_analytic with eta = 1            one counter, measured on the target run
  3  constant eta                       rung 2 plus one fitted scalar
  4  Random Forest on all counters      five counter sets, measured on target

Rungs 0 and 1 are the only ones that are genuinely predictive in the sense the
feasibility study meant. Rungs 2 to 4 all presuppose an instrumented run of the
application on the target machine.
"""
import warnings
from pathlib import Path

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

warnings.filterwarnings("ignore")

D = "/work/project/project/user"
DATA = f"{D}/repo/user/dissertation/data"
OUT = Path(f"{D}/analysis/out")
OUT.mkdir(parents=True, exist_ok=True)
F_PEAK = 2.25e9

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"]

CFG_FEATURES = ["log_ncore", "log_nrank", "log_nthread"]


def load_merged():
    """One row per configuration: union of counters over the five sets, median
    runtime. Per-set rows are ~58% imputed because each set measures only its
    own five counters, so merging is not merely a de-duplication step."""
    df = pd.read_csv(f"{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)
    return m[m.PAPI_TOT_CYC.notna()].reset_index(drop=True)


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))
    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 rf_pipeline(seed=0):
    return Pipeline([("impute", SimpleImputer(strategy="median")),
                     ("log", FunctionTransformer(slog, validate=False)),
                     ("scale", StandardScaler()),
                     ("est", RandomForestRegressor(n_estimators=400,
                                                   min_samples_leaf=2,
                                                   random_state=seed,
                                                   n_jobs=-1))])


def fac(pred, act):
    """Symmetric multiplicative error. Runtime spans 0.1 s to 215 s, so an
    absolute metric would be decided entirely by the slowest applications."""
    pred = np.clip(np.asarray(pred, float), 1e-9, None)
    return np.maximum(pred / act, act / pred)


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())

    lines = []
    say = lambda s="": (print(s), lines.append(s))

    say("Baseline ladder: what each layer of measurement actually buys")
    say("=" * 78)
    say(f"{len(df)} merged configurations, {len(apps)} applications: "
        f"{', '.join(apps)}")
    say(f"f_peak = {F_PEAK:.3g} Hz (ARCHER2 EPYC 7742)")
    say("Protocol: leave-one-application-out, error = max(pred/act, act/pred)")
    say("")

    # How much of the answer is already contained in the measured cycle count,
    # before any model is fitted. If this is near 1 the modelling contribution
    # is cosmetic.
    lr = np.corrcoef(np.log10(df.PAPI_TOT_CYC.values), np.log10(y))[0, 1]
    say(f"log-log correlation r(PAPI_TOT_CYC, runtime) = {lr:.4f} "
        f"(r^2 = {lr**2:.4f})")
    say(f"eta = t_analytic/runtime: median {np.median(eta):.3f}, "
        f"IQR {np.percentile(eta,25):.3f}-{np.percentile(eta,75):.3f}, "
        f"sd(log10 eta) {np.std(np.log10(eta)):.3f}")
    say("")

    rungs = []          # (label, cost, per-config error array)

    def add(label, cost, pred_fn):
        r = np.full(len(y), np.nan)
        for a in apps:
            te = (df.app == a).values
            r[te] = fac(pred_fn(~te, te), y[te])
        rungs.append((label, cost, r))
        return r

    # Rung 0: no measurement of the target configuration whatsoever. The median
    # is the correct constant because the metric is symmetric in log space.
    add("0  median training runtime", "nothing",
        lambda tr, te: np.full(te.sum(), np.median(y[tr])))

    # Rung 1: free configuration metadata. Note this is still a real regression
    # problem, it just has no counters. Under LOAO the model must extrapolate a
    # scaling curve for an application it has never seen, which is why it can
    # do worse than the constant.
    def cfg_only(tr, te):
        m = rf_pipeline()
        m.fit(X.loc[tr, CFG_FEATURES], np.log10(y[tr]))
        return 10 ** m.predict(X.loc[te, CFG_FEATURES])
    add("1  configuration only (RF)", "free metadata", cfg_only)

    # Rung 2: the analytic term alone, no fitting at all. This is where the
    # circularity enters: PAPI_TOT_CYC comes from an instrumented run of this
    # exact configuration on this exact machine.
    add("2  t_analytic, eta = 1", "1 counter, on target",
        lambda tr, te: t_an[te])

    # Rung 3: rung 2 corrected by a single fitted scalar, the median training
    # efficiency. Zero features, one parameter.
    add("3  t_analytic / constant eta", "1 counter, on target",
        lambda tr, te: t_an[te] / 10 ** np.median(np.log10(eta[tr])))

    # Rung 4: the full learned model. Same analytic backbone, eta now a
    # function of 21 counter-derived features.
    def rf_eta(tr, te):
        m = rf_pipeline()
        m.fit(X[tr], np.log10(eta[tr]))
        return t_an[te] / 10 ** m.predict(X[te])
    add("4  t_analytic / RF(counters)", "5 counter sets, on target", rf_eta)

    say("--- the ladder ---")
    say(f"{'rung':32s} {'measurement cost':22s} {'median':>8s} {'p90':>8s} "
        f"{'d.med':>8s} {'d.p90':>8s} {'p':>10s}")
    tab = []
    prev = None
    for label, cost, r in rungs:
        med, p90 = np.nanmedian(r), np.nanpercentile(r, 90)
        if prev is None:
            dmed = dp90 = np.nan
            pval = np.nan
        else:
            dmed = np.nanmedian(prev) - med
            dp90 = np.nanpercentile(prev, 90) - p90
            msk = ~np.isnan(r) & ~np.isnan(prev)
            try:
                pval = wilcoxon(r[msk], prev[msk]).pvalue
            except ValueError:
                pval = np.nan
        say(f"{label:32s} {cost:22s} {med:8.3f} {p90:8.3f} "
            f"{dmed:8.3f} {dp90:8.3f} {pval:10.2}")
        tab.append(dict(rung=label, cost=cost, median=round(med, 4),
                        p90=round(p90, 4),
                        delta_median_vs_prev=None if np.isnan(dmed) else round(dmed, 4),
                        delta_p90_vs_prev=None if np.isnan(dp90) else round(dp90, 4),
                        wilcoxon_p_vs_prev=None if np.isnan(pval) else float(f"{pval:.3g}")))
        prev = r

    say("")
    say("d.med and d.p90 are the reduction in error contributed by that rung")
    say("relative to the rung above it, so positive is an improvement.")
    say("")

    # The single number the scoping argument turns on.
    r0 = rungs[0][2]
    r2 = rungs[2][2]
    r4 = rungs[4][2]
    span = np.nanmedian(r0) - np.nanmedian(r4)
    say("--- attribution of the total gain ---")
    say(f"total gain, rung 0 to rung 4: {span:.3f} error factor")
    for i in range(1, len(rungs)):
        d = np.nanmedian(rungs[i - 1][2]) - np.nanmedian(rungs[i][2])
        say(f"  {rungs[i][0]:32s} contributes {d:7.3f}  "
            f"({100*d/span:5.1f}% of the total)")
    say("")
    say("Rung 2 is a measurement, not a prediction. Everything from rung 2")
    say("downwards requires an instrumented run of the target configuration.")

    say("")
    say("--- per-application median error ---")
    per = pd.DataFrame({lab: r for lab, _, r in rungs})
    per["app"] = df.app.values
    say(per.groupby("app").median().round(3).to_string())

    (OUT / "baseline_ladder.txt").write_text("\n".join(lines) + "\n")
    pd.DataFrame(tab).to_csv(OUT / "baseline_ladder.csv", index=False)
    per.to_csv(OUT / "baseline_ladder_perconfig.csv", index=False)
    print(f"\nwrote {OUT}/baseline_ladder.txt and .csv")


if __name__ == "__main__":
    main()