hardware-counters

#!/usr/bin/env python3
"""
Neural-network approach to performance prediction, compared against the
classical models from train_model.py.

The dataset is 40 configurations with 17 features. That is a difficult regime
for neural networks: a single hidden layer of 32 units already has ~600
parameters, an order of magnitude more than there are training samples, so the
network can memorise the training applications without learning anything
transferable. This script therefore does three things rather than just fitting
one MLP:

  1. sweeps architecture (depth/width) and regularisation strength, so the
     comparison is not sabotaged by an arbitrary bad choice of hyperparameters;
  2. evaluates under the same leave-one-application-out protocol as the other
     models, so the numbers are directly comparable;
  3. reports the gap between training and test error, which is the direct
     evidence of over-fitting.

As in train_model.py the target is the efficiency factor
    eta = (cycles / peak_clock) / runtime
and runtime is reconstructed as cycles / (peak_clock * eta).
"""
import warnings, json
from pathlib import Path

import numpy as np
import pandas as pd
from sklearn.neural_network import MLPRegressor
from sklearn.ensemble import RandomForestRegressor
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
D = "/work/project/project/user"
OUT = Path(f"{D}/analysis/out"); OUT.mkdir(parents=True, exist_ok=True)

import sys
sys.path.insert(0, f"{D}/analysis")
from train_model import build_features            # reuse identical features


def load():
    df = pd.read_csv(f"{D}/data/runs.csv")
    df = df[df.app != "cp2k"]
    df = df[df.runtime_s.notna()].reset_index(drop=True)
    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_an = cyc / F_PEAK
    return df, build_features(df), y, t_an, t_an / y     # df, X, y, analytic, eta


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


def mlp(hidden, alpha, lr=1e-3, iters=20000):
    return pipe(MLPRegressor(hidden_layer_sizes=hidden, alpha=alpha,
                             learning_rate_init=lr, max_iter=iters,
                             solver="lbfgs" if len(hidden) == 1 and hidden[0] <= 16
                                    else "adam",
                             random_state=RNG, early_stopping=False))


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


def loao(model_fn, X, y, t_an, eta, groups, seeds=(0, 1, 2, 3, 4)):
    """
    Leave-one-application-out, averaged over random seeds.

    Neural networks are sensitive to weight initialisation, especially with
    little data, so a single seed can be badly misleading in either direction.
    Averaging over 5 seeds gives a fairer estimate and exposes the variance.
    """
    te_f, tr_f = [], []
    for app in sorted(groups.unique()):
        te = (groups == app).values
        tr = ~te
        for s in seeds:
            m = model_fn(s)
            m.fit(X[tr], np.log10(eta[tr]))
            pt = t_an[te] / (10 ** m.predict(X[te]))
            pr = t_an[tr] / (10 ** m.predict(X[tr]))
            te_f.append(np.median(factors(pt, y[te])))
            tr_f.append(np.median(factors(pr, y[tr])))
    return float(np.median(te_f)), float(np.median(tr_f)), float(np.std(te_f))


def main():
    df, X, y, t_an, eta = load()
    groups = df.app
    print(f"{len(df)} configurations, {X.shape[1]} features, "
          f"{groups.nunique()} applications\n")

    rows = []

    # ---- reference points ------------------------------------------------
    for name, fn in [
        ("Constant efficiency (baseline)",
         lambda s: pipe(DummyRegressor(strategy="mean"))),
        ("Random Forest (best classical)",
         lambda s: pipe(RandomForestRegressor(n_estimators=500,
                                              min_samples_leaf=2,
                                              random_state=s, n_jobs=-1))),
    ]:
        t, tr, sd = loao(fn, X, y, t_an, eta, groups, seeds=(0,))
        rows.append({"model": name, "params": "-", "test_factor": t,
                     "train_factor": tr, "test_sd": sd})
        print(f"{name:34s} test {t:.3f}  train {tr:.3f}")

    print()

    # ---- neural network sweep -------------------------------------------
    archs = [(8,), (16,), (32,), (64,), (16, 8), (32, 16), (64, 32), (32, 16, 8)]
    alphas = [0.01, 0.1, 1.0, 10.0]
    best = None
    for h in archs:
        for a in alphas:
            nparam = sum(i * o for i, o in
                         zip((X.shape[1],) + h, h + (1,))) + sum(h) + 1
            t, tr, sd = loao(lambda s, h=h, a=a:
                             pipe(MLPRegressor(hidden_layer_sizes=h, alpha=a,
                                               max_iter=20000, random_state=s)),
                             X, y, t_an, eta, groups)
            tag = f"MLP {h} alpha={a}"
            rows.append({"model": tag, "params": nparam, "test_factor": t,
                         "train_factor": tr, "test_sd": sd})
            print(f"{tag:34s} test {t:.3f} +-{sd:.3f}  train {tr:.3f}  "
                  f"({nparam} params)")
            if best is None or t < best["test_factor"]:
                best = rows[-1]

    res = pd.DataFrame(rows)
    res.to_csv(OUT / "nn_sweep.csv", index=False)

    # ---- report ----------------------------------------------------------
    L = []
    L.append("Neural-network comparison")
    L.append("=" * 78)
    L.append(f"{len(df)} configurations, {X.shape[1]} features, "
             f"{groups.nunique()} applications")
    L.append("Leave-one-application-out; MLP results are the median over 5 seeds.")
    L.append("Error is the multiplicative factor; 1.00x is perfect.")
    L.append("")
    L.append(res.sort_values("test_factor").round(3).to_string(index=False))
    L.append("")

    rf   = res[res.model.str.startswith("Random Forest")].iloc[0]
    base = res[res.model.str.startswith("Constant")].iloc[0]
    L.append(f"Best MLP:        {best['model']}  test {best['test_factor']:.3f}"
             f"  (train {best['train_factor']:.3f}, {best['params']} parameters)")
    L.append(f"Random Forest:   test {rf.test_factor:.3f}  (train {rf.train_factor:.3f})")
    L.append(f"Constant eta:    test {base.test_factor:.3f}")
    L.append("")

    nn_only = res[res.model.str.startswith("MLP")]
    L.append(f"MLP test error across all {len(nn_only)} configurations swept: "
             f"{nn_only.test_factor.min():.2f}x to {nn_only.test_factor.max():.2f}x")
    L.append(f"Median seed-to-seed sd of MLP test error: "
             f"{nn_only.test_sd.median():.3f}")
    L.append("")
    gap = (nn_only.test_factor - nn_only.train_factor)
    L.append(f"Train/test gap for MLPs: median {gap.median():.3f} "
             f"(positive = fits training data better than held-out application)")
    L.append(f"Train/test gap for Random Forest: "
             f"{rf.test_factor - rf.train_factor:.3f}")

    txt = "\n".join(L)
    (OUT / "nn_summary.txt").write_text(txt + "\n")
    print("\n" + txt)


if __name__ == "__main__":
    main()