hardware-counters

#!/usr/bin/env python3
"""
Non-circular cross-machine projection.

Every result elsewhere in this project reconstructs runtime from PAPI_TOT_CYC
measured on the configuration being predicted. That is a decomposition of a
measurement, not a projection. This script asks the question the feasibility
study actually posed: given an application profiled on a SOURCE machine
(ARCHER2), can its runtime on a TARGET machine (Cirrus) be predicted without
ever running it on the target?

Nothing measured on Cirrus is allowed into the feature vector. Only:
  - ARCHER2 runtime and ARCHER2 hardware counters for the same app/size/config
  - configuration metadata (cores, ranks, threads)
  - static, datasheet-level facts about the two machines

Evaluation is leave-one-application-out over the six applications present on
both platforms, so the test application is unseen on BOTH machines. The
reference points are deliberately trivial: assume the runtime is unchanged, and
scale it by the clock ratio. If a counter-driven model cannot beat those, the
counters are not carrying transferable information about the hardware change.

A caveat that must not be glossed: with exactly one target machine the static
target specs are constant across every row, so they cannot contribute anything
a fitted intercept does not already capture. They are included for the sake of
the interface but their coefficient is unidentifiable. Separating "machine
description" from "fitted constant" needs at least three target machines. That
is a limitation of the data, not of the formulation.
"""
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.linear_model import Ridge
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_A2 = 2.25e9       # EPYC 7742 Rome, ARCHER2 compute node
F_CIR = 3.71e9      # EPYC 9825 Turin, Cirrus compute node (NOT the 9745 login node)

# Static target-machine description. Peak clock and cores per node are from the
# system documentation. Memory bandwidth is the datasheet STREAM-triad ceiling
# per socket: Rome is 8 channels of DDR4-3200 (204.8 GB/s), Turin is 12
# channels of DDR5-6000 (576 GB/s). These are nominal peaks, not measured, and
# are used only as a machine descriptor.
SPECS = {
    "archer2": dict(f_peak=F_A2, cores_node=128, sockets=2,
                    mem_bw_socket=204.8, ddr_gen=4, zen_gen=2),
    "cirrus": dict(f_peak=F_CIR, cores_node=288, sockets=2,
                   mem_bw_socket=576.0, ddr_gen=5, zen_gen=5),
}

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

KEYS = ["app", "size", "ncore", "nthread", "nrank"]


def merge_platform(df, counters):
    """Same merge as the single-platform work: union of counters across the
    five rotating sets, median runtime. The rotation exists because the EPYC
    exposes only five counters at once."""
    df = df[df.runtime_s.notna() & (df.runtime_s > 0)]
    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 load_matched():
    a2 = merge_platform(pd.read_csv(f"{DATA}/runs_expanded.csv"), A2_COUNTERS)
    cir_raw = pd.read_csv(f"{DATA}/cirrus_persets.csv")
    cir = merge_platform(cir_raw, ["PAPI_TOT_CYC", "PAPI_TOT_INS"])
    cir = cir[KEYS + ["runtime_s", "PAPI_TOT_CYC"]].rename(
        columns={"runtime_s": "rt_cir", "PAPI_TOT_CYC": "cyc_cir"})
    m = a2.merge(cir, on=KEYS, how="inner").reset_index(drop=True)
    return a2, cir, m


def source_features(df):
    """Features from the SOURCE machine only. cyc_cir and rt_cir never appear.

    The counter ratios are dimensionless descriptions of the workload (how
    memory-bound, how branch-heavy, how much the front end stalls). The
    hypothesis under test is that these describe the application well enough to
    say how it will respond to a different microarchitecture."""
    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["src_log_runtime"] = np.log10(df.runtime_s.clip(lower=1e-6))
    X["src_log_t_analytic"] = np.log10((cyc / F_A2).clip(lower=1e-9))
    X["src_eta"] = (cyc / F_A2) / df.runtime_s
    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))

    # Static machine descriptors. Constant within this experiment (see the
    # module docstring) but they are what a multi-target version would vary.
    s, t = SPECS["archer2"], SPECS["cirrus"]
    X["clock_ratio"] = t["f_peak"] / s["f_peak"]
    X["membw_ratio"] = t["mem_bw_socket"] / s["mem_bw_socket"]
    X["cores_node_ratio"] = t["cores_node"] / s["cores_node"]
    X["zen_gen_delta"] = t["zen_gen"] - s["zen_gen"]
    return X.replace([np.inf, -np.inf], np.nan)


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


def pipe(est):
    return Pipeline([("impute", SimpleImputer(strategy="median")),
                     ("log", FunctionTransformer(slog, validate=False)),
                     ("scale", StandardScaler()),
                     ("est", est)])


def fac(pred, act):
    pred = np.clip(np.asarray(pred, float), 1e-9, None)
    return np.maximum(pred / act, act / pred)


def main():
    a2, cir, m = load_matched()
    X = source_features(m)
    y = m.rt_cir.values
    t_src = m.runtime_s.values
    apps = sorted(m.app.unique())

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

    say("Non-circular cross-machine projection: ARCHER2 -> Cirrus")
    say("=" * 78)
    say(f"ARCHER2 merged configurations : {len(a2)}")
    say(f"Cirrus  merged configurations : {len(cir)}")
    say(f"matched on {'/'.join(KEYS)}   : {len(m)}")
    say(f"applications on both machines : {len(apps)} "
        f"({', '.join(apps)})")
    say("gromacs and openfoam are ARCHER2-only and are excluded.")
    say(f"features: {X.shape[1]}, all derived from ARCHER2 plus static specs")
    say("NO Cirrus counter and NO Cirrus runtime enters any feature.")
    say("")

    say(f"{'app':10s} {'n':>4s} {'med t_A2':>10s} {'med t_CIR':>10s} "
        f"{'speedup':>9s} {'sd log10':>9s}")
    for a in apps + ["ALL"]:
        s = m if a == "ALL" else m[m.app == a]
        sp = s.runtime_s.values / s.rt_cir.values
        say(f"{a:10s} {len(s):4d} {np.median(s.runtime_s):10.3f} "
            f"{np.median(s.rt_cir):10.3f} {np.median(sp):9.3f} "
            f"{np.std(np.log10(sp)):9.3f}")
    say("")
    say(f"clock ratio f_cir/f_a2 = {F_CIR/F_A2:.3f}, so a pure clock argument")
    say(f"predicts a uniform {F_CIR/F_A2:.3f}x speedup.")
    say("")

    results = []

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

    # --- reference points that use no learning at all -------------------
    add("t_cir = t_a2 (no adjustment)", "reference",
        lambda tr, te: t_src[te])

    add("t_cir = t_a2 * f_a2/f_cir (clock ratio)", "reference",
        lambda tr, te: t_src[te] * F_A2 / F_CIR)

    # An honest floor: no source measurement at all, just the typical Cirrus
    # runtime of the other five applications.
    add("median training Cirrus runtime", "reference",
        lambda tr, te: np.full(te.sum(), np.median(y[tr])))

    # --- one fitted scalar ----------------------------------------------
    # The empirical speedup, learned from five applications and applied to the
    # sixth. This is the strongest thing that can be said without counters.
    add("t_a2 * constant speedup (fitted)", "1 parameter",
        lambda tr, te: t_src[te] / 10 ** np.median(
            np.log10(t_src[tr] / y[tr])))

    # --- learned models on ARCHER2 counters ------------------------------
    # Target the log speedup rather than the log runtime: the scale of the
    # runtime is already in t_src, so this asks the model only for the part
    # that is genuinely about the machine change.
    def rf_ratio(tr, te):
        mdl = pipe(RandomForestRegressor(n_estimators=400, min_samples_leaf=2,
                                         random_state=0, n_jobs=-1))
        mdl.fit(X[tr], np.log10(t_src[tr] / y[tr]))
        return t_src[te] / 10 ** mdl.predict(X[te])
    add("RF on A2 counters -> speedup", "learned", rf_ratio)

    def ridge_ratio(tr, te):
        mdl = pipe(Ridge(alpha=10.0))
        mdl.fit(X[tr], np.log10(t_src[tr] / y[tr]))
        return t_src[te] / 10 ** mdl.predict(X[te])
    add("Ridge on A2 counters -> speedup", "learned", ridge_ratio)

    # Direct regression on log Cirrus runtime. Less well posed (the model must
    # also reproduce the runtime scale) but it is the obvious thing to try.
    def rf_direct(tr, te):
        mdl = pipe(RandomForestRegressor(n_estimators=400, min_samples_leaf=2,
                                         random_state=0, n_jobs=-1))
        mdl.fit(X[tr], np.log10(y[tr]))
        return 10 ** mdl.predict(X[te])
    add("RF on A2 counters -> log t_cir", "learned", rf_direct)

    # Analytic transfer: assume the ARCHER2 cycle count carries over to Cirrus
    # unchanged and only the clock differs, with efficiency learned. This is
    # the classic textbook projection and is worth separating from the purely
    # empirical speedup because it uses a counter rather than a runtime.
    cyc_a2 = m.PAPI_TOT_CYC.values
    def analytic_transfer(tr, te):
        eta_tr = (cyc_a2[tr] / F_CIR) / y[tr]
        return (cyc_a2[te] / F_CIR) / 10 ** np.median(np.log10(eta_tr))
    add("A2 cycles / f_cir, constant eta", "learned", analytic_transfer)

    def analytic_rf(tr, te):
        eta_tr = (cyc_a2[tr] / F_CIR) / y[tr]
        mdl = pipe(RandomForestRegressor(n_estimators=400, min_samples_leaf=2,
                                         random_state=0, n_jobs=-1))
        mdl.fit(X[tr], np.log10(eta_tr))
        return (cyc_a2[te] / F_CIR) / 10 ** mdl.predict(X[te])
    add("A2 cycles / f_cir, RF eta", "learned", analytic_rf)

    # Fold sizes are badly unbalanced (hpcg 44 configurations, lulesh 13) and
    # hpcg is also the one application whose behaviour is atypical, so the
    # pooled median is largely an hpcg statistic. The macro-average gives each
    # application one vote and can and does reverse the ranking.
    def macro(r):
        return float(np.mean([np.nanmedian(r[(m.app == a).values]) for a in apps]))

    say("--- leave-one-application-out, error = max(pred/act, act/pred) ---")
    say(f"{'predictor':42s} {'kind':11s} {'median':>8s} {'p90':>8s} "
        f"{'macro':>8s}")
    for label, kind, r in results:
        say(f"{label:42s} {kind:11s} {np.nanmedian(r):8.3f} "
            f"{np.nanpercentile(r, 90):8.3f} {macro(r):8.3f}")
    say("")
    say("'median' and 'p90' pool all 166 configurations, so hpcg (44 of them)")
    say("dominates. 'macro' is the mean of the six per-application medians.")
    say("")

    # Significance against the strongest non-learned reference.
    best_ref = min([x for x in results if x[1] == "reference"],
                   key=lambda x: np.nanmedian(x[2]))
    say(f"--- paired Wilcoxon against the best reference "
        f"('{best_ref[0]}', n={len(y)}) ---")
    for label, kind, r in results:
        if label == best_ref[0]:
            continue
        msk = ~np.isnan(r) & ~np.isnan(best_ref[2])
        try:
            p = wilcoxon(r[msk], best_ref[2][msk]).pvalue
        except ValueError:
            p = np.nan
        wins = int((r[msk] < best_ref[2][msk]).sum())
        verdict = "better" if np.nanmedian(r) < np.nanmedian(best_ref[2]) else "WORSE"
        say(f"  {label:42s} wins {wins:3d}/{int(msk.sum())}  p={p:9.3g}  {verdict}")
    say("")

    # The circular figure, computed on exactly this subset, so the comparison
    # is like for like. This is the number the rest of the project reports.
    say("--- the circular reference, on the same 166 configurations ---")
    eta_cir = (m.cyc_cir.values / F_CIR) / y
    r_circ = np.full(len(y), np.nan)
    for a in apps:
        te = (m.app == a).values
        r_circ[te] = fac((m.cyc_cir.values[te] / F_CIR)
                         / 10 ** np.median(np.log10(eta_cir[~te])), y[te])
    say(f"  {'CIRRUS cycles / f_cir, constant eta':42s} {'circular':11s} "
        f"{np.nanmedian(r_circ):8.3f} {np.nanpercentile(r_circ, 90):8.3f} "
        f"{macro(r_circ):8.3f}")
    say("  This uses PAPI_TOT_CYC measured on Cirrus, i.e. it presupposes the")
    say("  run it claims to predict. It is the ceiling the non-circular")
    say("  predictors are being asked to approach, and the gap is the honest")
    say("  cost of not having run on the target.")
    say("")

    # Does the ARCHER2 cycle count survive the move? If cycles were invariant
    # the analytic transfer would be exact, so the spread of this ratio is a
    # direct measure of how much microarchitecture matters.
    cr = m.cyc_cir.values / cyc_a2
    say("--- how far is the cycle count from being machine invariant? ---")
    say(f"  cyc_cir / cyc_a2: median {np.median(cr):.3f}, "
        f"IQR {np.percentile(cr,25):.3f}-{np.percentile(cr,75):.3f}, "
        f"sd(log10) {np.std(np.log10(cr)):.3f}")
    for a in apps:
        s = cr[(m.app == a).values]
        say(f"    {a:10s} median {np.median(s):6.3f}  "
            f"range {s.min():6.3f}-{s.max():6.3f}")
    say("")

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

    # What would a perfect per-application constant achieve? This separates
    # "the speedup is unpredictable" from "the speedup is predictable but not
    # from an unseen application".
    say("--- oracle diagnostics ---")
    orc = np.full(len(y), np.nan)
    for a in apps:
        te = (m.app == a).values
        orc[te] = fac(t_src[te] / 10 ** np.median(np.log10(t_src[te] / y[te])),
                      y[te])
    say(f"  oracle per-application constant speedup   median "
        f"{np.nanmedian(orc):.3f}  p90 {np.nanpercentile(orc,90):.3f}")
    say("  (fitted on the test application itself, so not achievable; it")
    say("   bounds what any method that only rescales t_a2 per application")
    say("   could reach)")
    gl = fac(t_src / 10 ** np.median(np.log10(t_src / y)), y)
    say(f"  oracle single global constant speedup     median "
        f"{np.nanmedian(gl):.3f}  p90 {np.nanpercentile(gl,90):.3f}")

    tab = pd.DataFrame([dict(predictor=l, kind=k,
                             median=round(float(np.nanmedian(r)), 4),
                             p90=round(float(np.nanpercentile(r, 90)), 4))
                        for l, k, r in results])
    (OUT / "noncircular.txt").write_text("\n".join(lines) + "\n")
    tab.to_csv(OUT / "noncircular.csv", index=False)
    per.to_csv(OUT / "noncircular_perconfig.csv", index=False)
    print(f"\nwrote {OUT}/noncircular.txt")


if __name__ == "__main__":
    main()