hardware-counters

#!/usr/bin/env python3
"""
Improved neural-network approach.

The first MLP attempt (train_nn.py) was significantly worse than Random Forest
and no better than a no-counter baseline. Diagnosing why identified five
distinct problems, each addressed here:

  H1  Dataset size. Each configuration was collapsed into a single row,
      discarding four of the five counter-set runs. But each set is an
      INDEPENDENT execution with its own Thread Time, so the usable dataset is
      ~190 samples, not 40. This is the single biggest fix and costs no extra
      compute -- the data was already on disk.

  H2  Feature skew. Four features had |skew| > 2 with max/min ratios up to
      5.7e4 (arith_intensity). StandardScaler centres these but leaves the
      tail, so a handful of configurations dominate the squared-error loss.
      Fixed with a signed log1p transform before scaling.

  H3  Structural missingness. 68 of 680 cells are missing, and the pattern is
      per-application (GROMACS 30, OpenFOAM 28, HPL/STREAM 0) because
      pat_report omits counters reading zero. Median imputation therefore
      fills a held-out application with the *other* applications' typical
      values -- systematically wrong exactly where it matters. Fixed by adding
      binary missingness indicators so the network can distinguish "absent"
      from "average".

  H4  Over-parameterisation. With 32 training rows even an 8-unit layer has
      153 parameters. Partly fixed by H1; also addressed with a bagged
      ensemble of small networks rather than one larger network.

  H5  Training protocol. adam with max_iter=20000 and no validation split runs
      to convergence on a tiny sample every time. lbfgs is deterministic and
      better suited to small n; seed-ensembling replaces early stopping.

Evaluation is unchanged: leave-one-application-out, error reported as the
multiplicative factor max(pred/actual, actual/pred).
"""
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, BaggingRegressor
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)


# ---------------------------------------------------------------- features
def features(df):
    """Scale-free behavioural ratios + missingness indicators (H3)."""
    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)
    X = X.replace([np.inf, -np.inf], np.nan)

    # H3: tell the model which counters were absent rather than imputed
    for c in list(X.columns):
        if X[c].isna().any():
            X[f"{c}__missing"] = X[c].isna().astype(float)
    return X


def signed_log1p(A):
    """H2: compress heavy tails while preserving sign and zero."""
    return np.sign(A) * np.log1p(np.abs(A))


def nn_pipe(hidden, alpha, seed, bag=1):
    """H2 + H5: log1p -> scale -> (bagged) lbfgs MLP."""
    base = MLPRegressor(hidden_layer_sizes=hidden, alpha=alpha,
                        solver="lbfgs", max_iter=5000, random_state=seed)
    est = base if bag == 1 else BaggingRegressor(
        estimator=base, n_estimators=bag, max_samples=0.8,
        random_state=seed, n_jobs=-1)
    return Pipeline([
        ("impute", SimpleImputer(strategy="median")),
        ("log",    FunctionTransformer(signed_log1p, validate=False)),
        ("scale",  StandardScaler()),
        ("est",    est),
    ])


def simple_pipe(est):
    return Pipeline([("impute", SimpleImputer(strategy="median")),
                     ("scale", StandardScaler()), ("est", est)])


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


# ---------------------------------------------------------------- data
def load(path, per_set):
    df = pd.read_csv(path)
    df = df[df.app != "cp2k"]
    if per_set:
        # set C has no PAPI_TOT_CYC, so no analytic term -> cannot form eta
        df = df[df.PAPI_TOT_CYC.notna()]
        y = df.runtime_s.values
    else:
        wc = [c for c in ["wall_B", "wall_C", "wall_D", "wall_E"] if c in df]
        df = df[df.runtime_s.notna()]
        y = df[wc].median(axis=1).fillna(df.runtime_s).values
    df = df.reset_index(drop=True)
    cyc = pd.to_numeric(df.PAPI_TOT_CYC, errors="coerce").values
    t_an = cyc / F_PEAK
    return df, features(df), y, t_an, t_an / y


def loao(fn, X, y, t_an, eta, groups, seeds=SEEDS):
    """Return per-row error factor averaged over seeds."""
    acc = np.zeros(len(y))
    for s in seeds:
        for app in sorted(groups.unique()):
            te = (groups == app).values
            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])
    return acc / len(seeds)


def main():
    L = []
    say = lambda s="": (print(s), L.append(s))

    say("Improved neural network: fixing the identified holes")
    say("=" * 78)

    results = {}

    # ---------- baseline: original 40-row dataset -------------------------
    df, X, y, t_an, eta = load(f"{D}/data/runs.csv", per_set=False)
    say(f"\n[1] ORIGINAL DATASET  {len(df)} rows, {X.shape[1]} features")
    g = df.app
    r_rf  = loao(lambda s: simple_pipe(RandomForestRegressor(
                 n_estimators=500, min_samples_leaf=2, random_state=s, n_jobs=-1)),
                 X, y, t_an, eta, g, seeds=(0,))
    r_bl  = loao(lambda s: simple_pipe(DummyRegressor(strategy="mean")),
                 X, y, t_an, eta, g, seeds=(0,))
    r_old = loao(lambda s: Pipeline([
                 ("i", SimpleImputer(strategy="median")),
                 ("sc", StandardScaler()),
                 ("e", MLPRegressor(hidden_layer_sizes=(32, 16, 8), alpha=10.0,
                                    max_iter=20000, random_state=s))]),
                 X, y, t_an, eta, g)
    r_new = loao(lambda s: nn_pipe((16, 8), 1.0, s, bag=10), X, y, t_an, eta, g)
    for k, v in [("Random Forest", r_rf), ("constant baseline", r_bl),
                 ("MLP original", r_old), ("MLP improved", r_new)]:
        say(f"    {k:22s} median {np.median(v):.3f}")
    results["orig"] = dict(rf=r_rf, bl=r_bl, old=r_old, new=r_new)

    # ---------- H1: per-counter-set dataset -------------------------------
    p2 = Path(f"{D}/data/runs_persets.csv")
    if not p2.exists():
        say("\n[2] per-set dataset not available yet - run parse_persets.py")
        return

    df2, X2, y2, t2, eta2 = load(p2, per_set=True)
    g2 = df2.app
    say(f"\n[2] PER-COUNTER-SET DATASET  {len(df2)} rows, {X2.shape[1]} features "
        f"({len(df2)/len(df):.1f}x more data, same experiments)")

    r2_rf  = loao(lambda s: simple_pipe(RandomForestRegressor(
                  n_estimators=500, min_samples_leaf=2, random_state=s, n_jobs=-1)),
                  X2, y2, t2, eta2, g2, seeds=(0,))
    r2_bl  = loao(lambda s: simple_pipe(DummyRegressor(strategy="mean")),
                  X2, y2, t2, eta2, g2, seeds=(0,))
    r2_old = loao(lambda s: Pipeline([
                  ("i", SimpleImputer(strategy="median")),
                  ("sc", StandardScaler()),
                  ("e", MLPRegressor(hidden_layer_sizes=(32, 16, 8), alpha=10.0,
                                     max_iter=20000, random_state=s))]),
                  X2, y2, t2, eta2, g2)

    say("\n    architecture / regularisation sweep for the improved MLP:")
    best, best_r = None, None
    for hid in [(16,), (32,), (16, 8), (32, 16), (64, 32)]:
        for al in [0.1, 1.0, 10.0]:
            r = loao(lambda s, h=hid, a=al: nn_pipe(h, a, s, bag=10),
                     X2, y2, t2, eta2, g2)
            med = np.median(r)
            say(f"      MLP {str(hid):10s} alpha={al:<5} bag=10   median {med:.3f}")
            if best is None or med < best:
                best, best_r, best_cfg = med, r, (hid, al)
    say(f"\n    best: MLP {best_cfg[0]} alpha={best_cfg[1]}  median {best:.3f}")

    for k, v in [("Random Forest", r2_rf), ("constant baseline", r2_bl),
                 ("MLP original", r2_old), ("MLP improved", best_r)]:
        say(f"    {k:22s} median {np.median(v):.3f}")

    # ---------- significance ---------------------------------------------
    say("\n[3] PAIRED SIGNIFICANCE (per-set dataset, Wilcoxon)")
    for name, a, b in [("improved MLP vs original MLP", best_r, r2_old),
                       ("improved MLP vs baseline",     best_r, r2_bl),
                       ("improved MLP vs RandomForest", best_r, r2_rf)]:
        try:
            _, p = wilcoxon(a, b)
        except ValueError:
            p = float("nan")
        say(f"    {name:32s} wins {int((a<b).sum()):3d}/{len(a)}  p = {p:.4f}")

    say("\n[4] PER-APPLICATION median error factor (per-set dataset)")
    tb = pd.DataFrame({"app": g2.values, "RandomForest": r2_rf,
                       "MLP_original": r2_old, "MLP_improved": best_r,
                       "baseline": r2_bl}).groupby("app").median().round(3)
    say(tb.to_string())

    pd.DataFrame({"app": g2.values, "ncore": df2.ncore.values,
                  "cset": df2.cset.values, "rf": r2_rf, "mlp_old": r2_old,
                  "mlp_new": best_r, "baseline": r2_bl}
                 ).to_csv(OUT / "nn2_perrow.csv", index=False)
    (OUT / "nn2_summary.txt").write_text("\n".join(L) + "\n")


if __name__ == "__main__":
    main()