hardware-counters

#!/usr/bin/env python3
"""
Retrain on the expanded dataset.

The expansion changed the problem in a way that matters for interpretation:
previously each application had exactly ONE problem size, so "application" and
"problem size" were perfectly confounded and a model could score well by
recognising five fixed operating points. Size now varies within each
application, which makes the task genuinely harder but the result meaningful.

Two validation protocols are therefore reported:

  LOAO   leave-one-application-out. Train on 7 codes, predict the 8th.
         Directly measures transfer to an unseen workload. Comparable to the
         earlier results.

  LOSO   leave-one-size-out (within application). Train on all applications
         but hold out the largest problem size of each. Measures whether the
         model extrapolates to a bigger problem than it has seen -- arguably
         the more useful capability in practice, since one usually has small
         runs of a code and wants to predict the large one.

Target remains the efficiency factor eta = (cycles / peak_clock) / runtime.
"""
import warnings, sys
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.linear_model import Ridge
from sklearn.dummy import DummyRegressor
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, FunctionTransformer
from sklearn.impute import SimpleImputer
from scipy.stats import wilcoxon

warnings.filterwarnings("ignore")
D      = "/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, 3, 4)


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


def build_features(df):
    """Scale-free behavioural ratios + configuration descriptors."""
    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")

    # configuration (free to know in advance, no profiling needed)
    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))

    # behavioural ratios
    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)

    # work per rank: instructions is a size proxy that transfers across codes
    X["log_instr_per_rank"]  = np.log10((ins / df["nrank"].clip(lower=1)).clip(lower=1))
    return X.replace([np.inf, -np.inf], np.nan)


def mlp(hidden=(32, 16), alpha=1.0):
    def fn(seed):
        return Pipeline([
            ("impute", SimpleImputer(strategy="median")),
            ("log",    FunctionTransformer(signed_log1p, validate=False)),
            ("scale",  StandardScaler()),
            ("est",    MLPRegressor(hidden_layer_sizes=hidden, alpha=alpha,
                                    solver="lbfgs", max_iter=5000,
                                    random_state=seed)),
        ])
    return fn


def plain(f):
    return lambda s: Pipeline([("impute", SimpleImputer(strategy="median")),
                               ("scale", StandardScaler()), ("est", f(s))])


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


def cv(fn, X, y, t_an, eta, folds, seeds=SEEDS):
    """folds: iterable of boolean test masks."""
    acc, cnt = np.zeros(len(y)), np.zeros(len(y))
    for s in seeds:
        for te in folds:
            if te.sum() == 0 or (~te).sum() == 0:
                continue
            m = fn(s)
            m.fit(X[~te], np.log10(eta[~te]))
            p = t_an[te] / (10 ** m.predict(X[te]))
            acc[te] += factors(p, y[te]); cnt[te] += 1
    return np.where(cnt > 0, acc / np.maximum(cnt, 1), np.nan)


def main():
    df = pd.read_csv(f"{D}/data/runs_expanded.csv")
    df = df[df.PAPI_TOT_CYC.notna() & df.runtime_s.notna()]
    df = df[df.runtime_s > 0].reset_index(drop=True)

    y    = df.runtime_s.values
    cyc  = pd.to_numeric(df.PAPI_TOT_CYC, errors="coerce").values
    t_an = cyc / F_PEAK
    eta  = t_an / y
    # drop physically impossible rows (eta > 1 means faster than peak clock)
    ok = (eta > 0) & (eta <= 1.5)
    df, y, t_an, eta = df[ok].reset_index(drop=True), y[ok], t_an[ok], eta[ok]
    X = build_features(df)

    L = []
    say = lambda s="": (print(s), L.append(s))
    say("Retraining on the expanded dataset")
    say("=" * 78)
    say(f"{len(df)} samples, {X.shape[1]} features, {df.app.nunique()} applications")
    say(f"runtime {y.min():.1f}-{y.max():.1f} s")
    say("")
    say("samples per application:")
    for a, n in df.app.value_counts().sort_index().items():
        sizes = df[df.app == a]["size"].nunique()
        thr   = df[df.app == a]["nthread"].nunique()
        say(f"  {a:9s} {n:4d}   ({sizes} problem size(s), {thr} thread config(s))")
    say("")

    models = {
        "Constant baseline": plain(lambda s: DummyRegressor(strategy="mean")),
        "Ridge":             plain(lambda s: Ridge(alpha=10.0)),
        "Random Forest":     plain(lambda s: RandomForestRegressor(
                                 n_estimators=500, min_samples_leaf=2,
                                 random_state=s, n_jobs=-1)),
        "Gradient Boosting": plain(lambda s: GradientBoostingRegressor(
                                 n_estimators=300, max_depth=3,
                                 learning_rate=0.05, random_state=s)),
        "MLP (32,16)":       mlp((32, 16), 1.0),
        "MLP (64,32)":       mlp((64, 32), 1.0),
        "MLP (64,32,16)":    mlp((64, 32, 16), 1.0),
    }

    # ---- protocol 1: leave-one-application-out --------------------------
    apps  = sorted(df.app.unique())
    loao  = [(df.app == a).values for a in apps]
    say("--- LOAO: leave one APPLICATION out (train on n-1 codes) ---")
    res_a = {}
    for name, fn in models.items():
        sd = (0,) if name in ("Random Forest", "Constant baseline",
                              "Ridge", "Gradient Boosting") else SEEDS
        r = cv(fn, X, y, t_an, eta, loao, sd)
        res_a[name] = r
        say(f"  {name:20s} median {np.nanmedian(r):.3f}   "
            f"p90 {np.nanpercentile(r,90):.3f}")
    say("")

    # ---- protocol 2: hold out the largest size of every application -----
    big = df.groupby("app")["size"].transform("max")
    te_big = ((df["size"] == big) & (df["size"] > 0)).values
    say(f"--- LOSO: hold out the LARGEST problem size of each app "
        f"({te_big.sum()} test rows) ---")
    say("    tests extrapolation to a bigger problem than seen in training")
    res_s = {}
    for name, fn in models.items():
        sd = (0,) if name in ("Random Forest", "Constant baseline",
                              "Ridge", "Gradient Boosting") else SEEDS
        r = cv(fn, X, y, t_an, eta, [te_big], sd)
        res_s[name] = r
        say(f"  {name:20s} median {np.nanmedian(r[te_big]):.3f}")
    say("")

    # ---- per-application detail for the two best ------------------------
    say("--- LOAO median error factor by application ---")
    tb = pd.DataFrame({k: v for k, v in res_a.items()})
    tb["app"] = df.app.values
    say(tb.groupby("app").median().round(3).to_string())
    say("")

    # ---- significance ----------------------------------------------------
    say("--- paired significance vs Random Forest (LOAO) ---")
    rf = res_a["Random Forest"]
    for name in ["MLP (32,16)", "MLP (64,32)", "MLP (64,32,16)",
                 "Gradient Boosting", "Constant baseline"]:
        v = res_a[name]
        m = ~np.isnan(v) & ~np.isnan(rf)
        try:
            _, p = wilcoxon(v[m], rf[m])
        except ValueError:
            p = float("nan")
        better = np.nanmedian(v) < np.nanmedian(rf)
        say(f"  {name:20s} wins {int((v[m]<rf[m]).sum()):4d}/{int(m.sum())}  "
            f"p={p:.4f}  ({'better' if better else 'worse'})")

    pd.DataFrame({**res_a, "app": df.app.values, "ncore": df.ncore.values,
                  "size": df["size"].values, "nthread": df.nthread.values}
                 ).to_csv(OUT / "retrain_perrow.csv", index=False)
    (OUT / "retrain_summary.txt").write_text("\n".join(L) + "\n")


if __name__ == "__main__":
    main()