hardware-counters

#!/usr/bin/env python3
"""
Is Cirrus intrinsically hard to predict, or was that an artefact of holding
problem sizes fixed?

The cross-platform experiment kept problem sizes identical to ARCHER2 so the
two machines ran the same work. On a 288-core Zen 5 node that made many runs
sub-second (miniFE 90.5%, HPL 65%, CoMD 64%), where fixed overheads dominate
and the efficiency factor collapses. No model beat a constant within Cirrus.

This compares three datasets under one protocol:
    archer2        191 configs, fixed sizes
    cirrus-fixed   191 configs, same sizes as ARCHER2  (comparability)
    cirrus-weak    ~165 configs, work per rank constant (predictability)

If cirrus-weak becomes predictable, the earlier result was a design artefact
and the lesson is that comparability and predictability trade off. If it does
not, Zen 5 is genuinely harder to model.
"""
import os, warnings, sys
from pathlib import Path
import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, FunctionTransformer
from sklearn.pipeline import Pipeline
from scipy.stats import wilcoxon

warnings.filterwarnings("ignore")
# Data and output roots. Overridable so the pipeline can be relocated
# (a different account, a scratch copy, another site) without editing code;
# the default keeps existing invocations working unchanged.
D   = os.environ.get("DISS_ROOT", "/work/project/project/user")
OUT = Path(f"{D}/analysis/out"); OUT.mkdir(parents=True, exist_ok=True)
F_PEAK = {"archer2": 2.25e9, "cirrus": 3.71e9}

sys.path.insert(0, f"{D}/analysis")
from crossplatform import COUNTERS, INTERSECTION, features, slog, fac


def load(path, platform, scaling):
    df = pd.read_csv(path)
    df["platform"] = platform
    df = df[df.runtime_s.notna() & (df.runtime_s > 0)]
    keys = ["platform", "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)
    m = m[m.PAPI_TOT_CYC.notna()].reset_index(drop=True)
    m["dataset"] = scaling
    return m


def mk(kind, seed=0):
    pre = [("i", SimpleImputer(strategy="median")),
           ("l", FunctionTransformer(slog, validate=False)),
           ("s", StandardScaler())]
    est = {"rf": RandomForestRegressor(n_estimators=400, min_samples_leaf=2,
                                       random_state=seed, n_jobs=-1),
           "gbq": GradientBoostingRegressor(loss="quantile", alpha=0.5,
                                            n_estimators=300, max_depth=3,
                                            learning_rate=0.05, random_state=seed),
           "mlp": MLPRegressor(hidden_layer_sizes=(32, 16), alpha=1.0,
                               solver="lbfgs", max_iter=5000, random_state=seed),
           }[kind]
    return Pipeline(pre + [("e", est)])


def loao(sub, X, kind):
    """Leave-one-application-out within one dataset."""
    y, t_an, eta = sub.y.values, sub.t_an.values, sub.eta.values
    err = np.full(len(sub), np.nan)
    for a in sorted(sub.app.unique()):
        te = (sub.app == a).values
        tr = ~te
        if te.sum() == 0 or tr.sum() < 10:
            continue
        if kind == "const":
            p = t_an[te] / (10 ** np.median(np.log10(eta[tr])))
        else:
            m = mk(kind); m.fit(X[tr], np.log10(eta[tr]))
            p = t_an[te] / (10 ** m.predict(X[te]))
        err[te] = fac(p, y[te])
    return err


def main():
    a2 = load(f"{D}/data/runs_expanded.csv",  "archer2", "archer2")
    cf = load(f"{D}/data/cirrus_persets.csv", "cirrus",  "cirrus-fixed")
    cw = load(f"{D}/data/cirrus_weak.csv",    "cirrus",  "cirrus-weak")
    allf = pd.concat([a2, cf, cw], ignore_index=True)

    X = features(allf)
    y = allf.runtime_s.values
    fpk = allf.platform.map(F_PEAK).values
    t_an = pd.to_numeric(allf.PAPI_TOT_CYC, errors="coerce").values / fpk
    eta = t_an / y
    ok = (eta > 0) & (eta <= 1.5)
    allf = allf[ok].reset_index(drop=True)
    X = X[ok].reset_index(drop=True)
    allf["y"], allf["t_an"], allf["eta"] = y[ok], t_an[ok], eta[ok]

    L = []
    say = lambda s="": (print(s, flush=True), L.append(s))
    say("Weak scaling: is Cirrus intrinsically hard, or was it the problem size?")
    say("=" * 78)

    say("--- dataset characteristics ---")
    hdr = f"  {'dataset':14s} {'n':>4s} {'med rt':>8s} {'<1s':>6s} {'sd log10(eta)':>14s} {'med eta':>8s}"
    say(hdr); say("  " + "-" * (len(hdr) - 2))
    for ds in ["archer2", "cirrus-fixed", "cirrus-weak"]:
        s = allf[allf.dataset == ds]
        if len(s) == 0:
            continue
        say(f"  {ds:14s} {len(s):4d} {s.y.median():8.2f} "
            f"{100*(s.y<1).mean():5.1f}% {np.log10(s.eta).std():14.4f} "
            f"{s.eta.median():8.3f}")
    say("")

    say("--- leave-one-application-out, within each dataset ---")
    rows = []
    for ds in ["archer2", "cirrus-fixed", "cirrus-weak"]:
        mask = (allf.dataset == ds).values
        if mask.sum() == 0:
            continue
        sub = allf[mask].reset_index(drop=True)
        for fs, cols in [("intersection", [c for c in INTERSECTION if c in X.columns]),
                         ("full", list(X.columns))]:
            Xs = X.loc[mask, cols].values
            res = {}
            for kind in ("const", "rf", "gbq", "mlp"):
                e = loao(sub, Xs, kind)
                res[kind] = e
                rows.append({"dataset": ds, "features": fs, "model": kind,
                             "n": int(mask.sum()),
                             "median": np.nanmedian(e),
                             "p90": np.nanpercentile(e, 90)})
            best = min(("rf", "gbq", "mlp"), key=lambda k: np.nanmedian(res[k]))
            c, b = np.nanmedian(res["const"]), np.nanmedian(res[best])
            m = ~np.isnan(res[best]) & ~np.isnan(res["const"])
            try:
                _, p = wilcoxon(res[best][m], res["const"][m])
            except ValueError:
                p = float("nan")
            say(f"  {ds:14s} {fs:13s} best={best:4s} {b:.4f}   const={c:.4f}   "
                f"p={p:.4g}  {'BEATS' if b < c else 'loses to'} constant")
    say("")

    res = pd.DataFrame(rows)
    res.to_csv(OUT / "weak_results.csv", index=False)
    say("--- all models ---")
    say(res.pivot_table(index=["dataset", "features"], columns="model",
                        values="median").round(4).to_string())
    say("")

    # transfer with the weak-scaled Cirrus data
    say("--- cross-platform transfer using WEAK-scaled Cirrus ---")
    cols = [c for c in INTERSECTION if c in X.columns]
    for src, dst in [("archer2", "cirrus-weak"), ("cirrus-weak", "archer2")]:
        trm = (allf.dataset == src).values
        tem = (allf.dataset == dst).values
        Xtr, Xte = X.loc[trm, cols].values, X.loc[tem, cols].values
        e_c = fac(allf.t_an.values[tem] /
                  (10 ** np.median(np.log10(allf.eta.values[trm]))),
                  allf.y.values[tem])
        line = f"  {src} -> {dst}:  const {np.nanmedian(e_c):.4f}"
        for kind in ("rf", "gbq", "mlp"):
            m = mk(kind); m.fit(Xtr, np.log10(allf.eta.values[trm]))
            e = fac(allf.t_an.values[tem] / (10 ** m.predict(Xte)),
                    allf.y.values[tem])
            line += f"   {kind} {np.nanmedian(e):.4f}"
        say(line)

    (OUT / "weak_summary.txt").write_text("\n".join(L) + "\n")


if __name__ == "__main__":
    main()