Opens a larger view. Escape closes it.

hardware-counters

matched_crossplatform.py

#!/usr/bin/env python3
"""
Cross-platform prediction on the EXACTLY MATCHED configurations.

WHY THIS FILE EXISTS
--------------------
P1.  dissertation/README.md states, of ARCHER2 versus Cirrus:

        "191 matched configurations each, identical codes and problem sizes."

     That is false in three ways.

       (i)  The two datasets happen to contain 191 configurations each.  That
            is a coincidence of how far each sweep got, not a matching.
       (ii) ARCHER2 has EIGHT applications (comd, gromacs, hpcg, hpl, lulesh,
            minife, openfoam, stream).  Cirrus has SIX: gromacs and openfoam
            were never built or run there.  A quarter of the ARCHER2 codes have
            no Cirrus counterpart at all.
      (iii) Even among the six shared applications the grids differ.  Cirrus
            includes core counts 125 and 288 (its node is 288-way) that ARCHER2
            (128-way) cannot reach, and ARCHER2 has an hpl size=0 entry Cirrus
            lacks.

     Joining on (app, size, ncore, nthread) gives 166 configurations present on
     BOTH machines, not 191.  Everything the earlier crossplatform.py reported
     is an UNPAIRED comparison of two differently-shaped 191-row samples, so
     any difference between the platforms is confounded with the difference in
     which configurations each sample contains.  In particular Cirrus's extra
     125- and 288-core points are its largest, most parallel, shortest-running
     runs, exactly the ones with the worst efficiency factor, so the "Cirrus is
     harder" conclusion was partly an artefact of the sample.

     This script redoes the analysis on the 166 genuinely matched
     configurations, as a paired comparison, and reports it SIDE BY SIDE with
     the unmatched 191-versus-191 version so the reader can see exactly what
     the restriction costs and what it changes.

P2.  crossplatform.py's `A2->CIR LOAO` and `CIR->A2 LOAO` settings loop over
     `both.app.unique()` (the union, eight apps) and carry

         if len(te) == 0 or len(tr) < 10:
             continue

     For app = gromacs and app = openfoam the Cirrus test mask is empty, so the
     fold is skipped in silence.  The reported n is still printed as 191 (it
     counts the mask, not the finite entries), and the constant reference for
     those settings is never computed at all, which is why the summary block
     prints `const=nan` with no explanation:

         A2->CIR LOAO  full  best=rf  1.2144  const=nan  loses to constant

     A verdict of "loses to constant" against a NaN is not a verdict.  Here
     every fold is accounted for explicitly: evaluated, or listed as NOT
     EVALUABLE with the reason, and the constant reference is computed on the
     same folds as the models so the comparison is like for like.
"""
import warnings, sys
from pathlib import Path

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.linear_model import Ridge
from sklearn.neural_network import MLPRegressor
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, FunctionTransformer
from sklearn.pipeline import Pipeline

warnings.filterwarnings("ignore")
D = "/work/project/project/user"
OUT = Path(f"{D}/analysis/out"); OUT.mkdir(parents=True, exist_ok=True)
sys.path.insert(0, f"{D}/analysis")

from crossplatform import load as xload, features as xfeatures, F_PEAK, INTERSECTION, slog, fac
from stats_rigour import Family, boot_median_ci

KEYS = ["app", "size", "ncore", "nthread"]
NJ = 8


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=NJ),
           "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),
           "ridge": Ridge(alpha=10.0)}[kind]
    return Pipeline(pre + [("e", est)])


class Bench:
    """One analysis universe (matched or unmatched), sharing all the plumbing."""

    def __init__(self, label, a2, cir):
        self.label = label
        both = pd.concat([a2, cir], ignore_index=True)
        X = xfeatures(both)
        y = both.runtime_s.values
        fpk = both.platform.map(F_PEAK).values
        t_an = pd.to_numeric(both.PAPI_TOT_CYC, errors="coerce").values / fpk
        eta = t_an / y
        ok = (eta > 0) & (eta <= 1.5)
        self.both = both[ok].reset_index(drop=True)
        self.X = X[ok].reset_index(drop=True)
        self.y, self.t_an, self.eta = y[ok], t_an[ok], eta[ok]
        self.isA2 = (self.both.platform == "archer2").values
        self.cols = [c for c in INTERSECTION if c in self.X.columns]
        self.Xf = self.X[self.cols].values
        self.n_dropped = int((~ok).sum())

    # ---- primitives ----
    def _fit(self, tr, te, kind):
        if kind == "const":
            return self.t_an[te] / (10 ** np.median(np.log10(self.eta[tr])))
        m = mk(kind); m.fit(self.Xf[tr], np.log10(self.eta[tr]))
        return self.t_an[te] / (10 ** m.predict(self.Xf[te]))

    def within_loao(self, mask, kind):
        e = np.full(len(self.y), np.nan)
        sub = np.where(mask)[0]
        for a in sorted(self.both.app[mask].unique()):
            te = sub[self.both.app.values[sub] == a]
            tr = sub[self.both.app.values[sub] != a]
            if len(te) == 0 or len(tr) < 10:
                continue
            e[te] = fac(self._fit(tr, te, kind), self.y[te])
        return e

    def transfer(self, trm, tem, kind):
        e = np.full(len(self.y), np.nan)
        e[tem] = fac(self._fit(np.where(trm)[0], np.where(tem)[0], kind),
                     self.y[tem])
        return e

    def transfer_loao(self, trm, tem, kind, audit=None):
        """P2: every fold accounted for.  Returns (errors, per-fold audit)."""
        e = np.full(len(self.y), np.nan)
        rep = []
        for a in sorted(self.both.app.unique()):
            te = np.where(tem & (self.both.app == a).values)[0]
            tr = np.where(trm & (self.both.app != a).values)[0]
            if len(te) == 0:
                rep.append((a, 0, len(tr), "NOT EVALUABLE: application absent "
                                           "on the TARGET platform"))
                continue
            if len(tr) < 10:
                rep.append((a, len(te), len(tr), f"NOT EVALUABLE: only {len(tr)}"
                            " training rows on the source platform (<10)"))
                continue
            e[te] = fac(self._fit(tr, te, kind), self.y[te])
            rep.append((a, len(te), len(tr), "evaluated"))
        if audit is not None:
            audit.extend(rep)
        return e


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

    a2 = xload(f"{D}/data/runs_expanded.csv", "archer2")
    cir = xload(f"{D}/data/cirrus_persets.csv", "cirrus")

    say("Cross-platform prediction on the EXACTLY MATCHED configurations")
    say("=" * 96)
    say("")
    say("P1: THE MATCHING CLAIM IS FALSE")
    say("")
    say(f"  README states: '191 matched configurations each, identical codes")
    say(f"  and problem sizes.'  What the data actually contain:")
    say("")
    say(f"    ARCHER2  {len(a2):3d} configurations, {a2.app.nunique()} applications: "
        f"{sorted(a2.app.unique())}")
    say(f"    Cirrus   {len(cir):3d} configurations, {cir.app.nunique()} applications: "
        f"{sorted(cir.app.unique())}")
    say("")
    say(f"    applications on ARCHER2 only : {sorted(set(a2.app) - set(cir.app))}")
    say(f"    applications on Cirrus only  : {sorted(set(cir.app) - set(a2.app))}")
    say(f"    core counts on Cirrus only   : {sorted(set(cir.ncore) - set(a2.ncore))}"
        "   (288-core node; ARCHER2 has 128)")
    say(f"    core counts on ARCHER2 only  : {sorted(set(a2.ncore) - set(cir.ncore))}")
    say("")
    key = pd.merge(a2[KEYS], cir[KEYS], on=KEYS)
    say(f"    inner join on {KEYS}")
    say(f"      EXACTLY MATCHED configurations : {len(key)}")
    say(f"      unmatched on ARCHER2           : {len(a2) - len(key)}")
    say(f"      unmatched on Cirrus            : {len(cir) - len(key)}")
    say("")
    say("    per application:")
    say(f"      {'app':10s} {'ARCHER2':>8s} {'Cirrus':>8s} {'matched':>8s}")
    for a in sorted(set(a2.app) | set(cir.app)):
        na, nc = int((a2.app == a).sum()), int((cir.app == a).sum())
        nm = int((key.app == a).sum())
        note = "  <- absent on Cirrus" if nc == 0 else ""
        say(f"      {a:10s} {na:8d} {nc:8d} {nm:8d}{note}")
    say("")
    say("  The two totals being 191 each is a coincidence of how far each sweep")
    say("  reached, not evidence of matching.  Only 166 configurations exist on")
    say("  both machines.  Comparisons over the two full 191-row samples are")
    say("  therefore UNPAIRED and confound platform with sample composition.")
    say("")
    say("  Why this biases the platform comparison in a specific direction:")
    cir_only = cir[~cir.set_index(KEYS).index.isin(key.set_index(KEYS).index)]
    cir_m = cir[cir.set_index(KEYS).index.isin(key.set_index(KEYS).index)]
    say(f"    Cirrus configurations that are matched     : n={len(cir_m):3d}  "
        f"median runtime {cir_m.runtime_s.median():7.3f} s")
    say(f"    Cirrus configurations with NO ARCHER2 twin : n={len(cir_only):3d}  "
        f"median runtime {cir_only.runtime_s.median():7.3f} s")
    say(f"    unmatched-only core counts: "
        f"{sorted(cir_only.ncore.unique())}")
    say("    The unmatched Cirrus rows are its widest and shortest runs, which")
    say("    are exactly the ones where fixed overhead dominates and the")
    say("    efficiency factor collapses.  Including them makes Cirrus look")
    say("    harder to predict for a reason that has nothing to do with Zen 5.")
    say("")

    # ---------------- build the two universes ----------------
    idx = key.set_index(KEYS).index
    a2m = a2[a2.set_index(KEYS).index.isin(idx)].reset_index(drop=True)
    cirm = cir[cir.set_index(KEYS).index.isin(idx)].reset_index(drop=True)
    a2m = a2m.sort_values(KEYS).reset_index(drop=True)
    cirm = cirm.sort_values(KEYS).reset_index(drop=True)
    assert (a2m[KEYS].values == cirm[KEYS].values).all(), "pairing broken"
    say(f"  Paired universe built and verified: {len(a2m)} ARCHER2 rows aligned")
    say(f"  row-for-row with {len(cirm)} Cirrus rows on {KEYS}.")
    say("")

    U = {"unmatched (191 v 191, as published)": Bench("unmatched", a2, cir),
         "matched (166 v 166, paired)": Bench("matched", a2m, cirm)}
    for nm, b in U.items():
        say(f"  {nm}: {len(b.y)} rows after the eta in (0,1.5] filter "
            f"({b.n_dropped} dropped), {len(b.cols)} intersection features")
    say("")

    # ---------------- descriptive, now genuinely paired ----------------
    say("=" * 96)
    say("PLATFORM DIFFERENCE ON THE MATCHED SET (a paired comparison at last)")
    say("")
    bm = U["matched (166 v 166, paired)"]
    ma, mc = bm.isA2, ~bm.isA2
    say(f"  {'quantity':28s} {'ARCHER2':>12s} {'Cirrus':>12s}")
    say("  " + "-" * 56)
    for nm, va, vc in [
        ("n", ma.sum(), mc.sum()),
        ("median runtime (s)", np.median(bm.y[ma]), np.median(bm.y[mc])),
        ("fraction under 1 s", (bm.y[ma] < 1).mean(), (bm.y[mc] < 1).mean()),
        ("median eta", np.median(bm.eta[ma]), np.median(bm.eta[mc])),
        ("sd log10(eta)", np.log10(bm.eta[ma]).std(), np.log10(bm.eta[mc]).std()),
    ]:
        say(f"  {nm:28s} {va:12.4f} {vc:12.4f}")
    bu = U["unmatched (191 v 191, as published)"]
    ua, uc = bu.isA2, ~bu.isA2
    say("")
    say("  the same on the UNMATCHED samples, for contrast:")
    say(f"  {'median runtime (s)':28s} {np.median(bu.y[ua]):12.4f} "
        f"{np.median(bu.y[uc]):12.4f}")
    say(f"  {'fraction under 1 s':28s} {(bu.y[ua]<1).mean():12.4f} "
        f"{(bu.y[uc]<1).mean():12.4f}")
    say(f"  {'sd log10(eta)':28s} {np.log10(bu.eta[ua]).std():12.4f} "
        f"{np.log10(bu.eta[uc]).std():12.4f}")
    say("")

    # ---------------- the six settings, both universes ----------------
    say("=" * 96)
    say("THE SIX SETTINGS, MATCHED VERSUS UNMATCHED, SIDE BY SIDE")
    say("")
    say("  intersection features only (the 14 computable on both machines).")
    say("  Error factor = max(pred/actual, actual/pred); [.,.] is a percentile")
    say("  bootstrap 95% CI on the median with 10000 resamples.")
    say("")

    rows = []
    audits = {}

    def emit(setting, kind, e_un, e_ma, mask_un, mask_ma):
        def cell(e, mask):
            v = e[mask]
            v = v[np.isfinite(v)]
            if len(v) == 0:
                return np.nan, np.nan, np.nan, 0
            lo, hi = boot_median_ci(v)
            return float(np.median(v)), lo, hi, len(v)
        mu, lu, hu, nu = cell(e_un, mask_un)
        mm, lm, hm, nm_ = cell(e_ma, mask_ma)
        rows.append({"setting": setting, "model": kind,
                     "n_unmatched": nu, "median_unmatched": mu,
                     "ci_lo_unmatched": lu, "ci_hi_unmatched": hu,
                     "n_matched": nm_, "median_matched": mm,
                     "ci_lo_matched": lm, "ci_hi_matched": hm,
                     "delta_matched_minus_unmatched": mm - mu})
        f = lambda m, l, h, n: (f"{'n/a':>22s}" if not np.isfinite(m)
                                else f"{m:.4f} [{l:.3f},{h:.3f}] n={n:3d}")
        say(f"  {setting:16s} {kind:6s}  unmatched {f(mu,lu,hu,nu)}   "
            f"matched {f(mm,lm,hm,nm_)}   delta {mm-mu:+.4f}"
            if np.isfinite(mm) and np.isfinite(mu) else
            f"  {setting:16s} {kind:6s}  unmatched {f(mu,lu,hu,nu)}   "
            f"matched {f(mm,lm,hm,nm_)}")

    # settings 1-2
    for setting, pick in [("within-A2", lambda b: b.isA2),
                          ("within-CIR", lambda b: ~b.isA2)]:
        for kind in ("rf", "gbq", "mlp", "const"):
            eu = bu.within_loao(pick(bu), kind)
            em = bm.within_loao(pick(bm), kind)
            emit(setting, kind, eu, em, pick(bu), pick(bm))
        say("")

    # settings 3-4
    for setting, srcf, dstf in [("A2->CIR", lambda b: b.isA2, lambda b: ~b.isA2),
                                ("CIR->A2", lambda b: ~b.isA2, lambda b: b.isA2)]:
        for kind in ("rf", "gbq", "mlp", "ridge", "const"):
            eu = bu.transfer(srcf(bu), dstf(bu), kind)
            em = bm.transfer(srcf(bm), dstf(bm), kind)
            emit(setting, kind, eu, em, dstf(bu), dstf(bm))
        say("")

    # settings 5-6, with the fold audit (P2)
    for setting, srcf, dstf in [("A2->CIR LOAO", lambda b: b.isA2, lambda b: ~b.isA2),
                                ("CIR->A2 LOAO", lambda b: ~b.isA2, lambda b: b.isA2)]:
        for kind in ("rf", "mlp", "const"):
            au, am = [], []
            eu = bu.transfer_loao(srcf(bu), dstf(bu), kind, au)
            em = bm.transfer_loao(srcf(bm), dstf(bm), kind, am)
            if kind == "rf":
                audits[setting] = (au, am)
            emit(setting, kind, eu, em, dstf(bu), dstf(bm))
        say("")

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

    # ---------------- P2: the fold audit ----------------
    say("=" * 96)
    say("P2: FOLD-BY-FOLD AUDIT OF THE TRANSFER-PLUS-UNSEEN-APPLICATION SETTINGS")
    say("")
    say("  crossplatform.py loops over the UNION of applications and skips any")
    say("  fold with an empty test set via a bare `continue`.  It then prints")
    say("  n=191 (the mask size, not the number of finite predictions) and")
    say("  never computes the constant reference for these settings, which is")
    say("  where the unexplained `const=nan` in crossplatform_summary.txt came")
    say("  from.  Every fold is now listed with its disposition.")
    say("")
    for setting, (au, am) in audits.items():
        for uname, rep in [("unmatched (191 v 191)", au), ("matched (166 v 166)", am)]:
            say(f"  {setting}, {uname}:")
            say(f"    {'app':10s} {'n_test':>7s} {'n_train':>8s}  disposition")
            for a, nt, ntr, note in rep:
                say(f"    {a:10s} {nt:7d} {ntr:8d}  {note}")
            ev = sum(1 for r in rep if r[3] == "evaluated")
            say(f"    -> {ev} of {len(rep)} folds evaluated, "
                f"{len(rep)-ev} not evaluable")
            say("")
    say("  DISPOSITION OF THE TWO MISSING APPLICATIONS")
    say("")
    say("    gromacs and openfoam are present on ARCHER2 and absent on Cirrus.")
    say("    For A2->CIR LOAO their test set is empty, so the fold cannot be")
    say("    scored: this is reported as NOT EVALUABLE (application absent on")
    say("    the target platform), not as a silent skip and not as a zero.")
    say("    For CIR->A2 LOAO their ARCHER2 test set is non-empty, but no Cirrus")
    say("    training rows for them exist either, so the fold reduces to")
    say("    predicting an application the source platform has never seen in any")
    say("    form.  That is still a legitimate (and harder) test and it IS")
    say("    scored here; the earlier code scored it too but did not say so.")
    say("")
    say("    The consequence for the headline claim: the A2->CIR LOAO figure of")
    say("    1.214 is an average over SIX applications, not eight, and the two")
    say("    that are missing (a production MD code and a production CFD code)")
    say("    are the two least like the proxy benchmarks.  The hardest setting")
    say("    in the experiment matrix is evaluated only on the easy half of the")
    say("    application space.  This must be stated wherever 1.214 is quoted.")
    say("")
    say("    The constant reference for these settings is now computed on the")
    say("    SAME folds as the models, so the 'BEATS/loses to constant' verdict")
    say("    is meaningful instead of a comparison against NaN.")
    say("")

    # ---------------- paired significance ----------------
    say("=" * 96)
    say("PAIRED SIGNIFICANCE ON THE MATCHED SET")
    say("")
    say("  Now that rows are aligned one-to-one, the model-versus-constant")
    say("  comparisons are genuinely paired.  Reported through stats_rigour.")
    say("")
    fams = []
    for setting, pick in [("within-A2 (matched)", lambda b: b.isA2),
                          ("within-CIR (matched)", lambda b: ~b.isA2)]:
        base = bm.within_loao(pick(bm), "const")
        fam = Family(f"{setting} LOAO vs constant, matched 166")
        for k in ("rf", "gbq", "mlp"):
            fam.add(k, bm.within_loao(pick(bm), k), base)
        fams.append(fam.finalise()); say(fam.table())

    fam = Family("Cross-platform transfer vs constant, matched 166")
    for setting, srcf, dstf in [("A2->CIR", lambda b: b.isA2, lambda b: ~b.isA2),
                                ("CIR->A2", lambda b: ~b.isA2, lambda b: b.isA2)]:
        base = bm.transfer(srcf(bm), dstf(bm), "const")
        for k in ("rf", "gbq", "mlp"):
            fam.add(f"{setting} {k}", bm.transfer(srcf(bm), dstf(bm), k), base)
    fams.append(fam.finalise()); say(fam.table())

    fam = Family("Transfer + unseen application vs constant, matched 166, "
                 "six evaluable folds")
    for setting, srcf, dstf in [("A2->CIR LOAO", lambda b: b.isA2, lambda b: ~b.isA2),
                                ("CIR->A2 LOAO", lambda b: ~b.isA2, lambda b: b.isA2)]:
        base = bm.transfer_loao(srcf(bm), dstf(bm), "const")
        for k in ("rf", "mlp"):
            fam.add(f"{setting} {k}", bm.transfer_loao(srcf(bm), dstf(bm), k), base)
    fams.append(fam.finalise()); say(fam.table())

    pd.concat([f.frame() for f in fams], ignore_index=True).to_csv(
        OUT / "matched_crossplatform_stats.csv", index=False)

    # ---------------- what the restriction cost ----------------
    say("=" * 96)
    say("WHAT THE RESTRICTION TO MATCHED CONFIGURATIONS COSTS AND CHANGES")
    say("")
    r = res.dropna(subset=["median_matched", "median_unmatched"])
    say(f"  {'setting':16s} {'model':6s} {'unmatched':>10s} {'matched':>10s} "
        f"{'delta':>9s}")
    say("  " + "-" * 56)
    for _, x in r.iterrows():
        say(f"  {x.setting:16s} {x.model:6s} {x.median_unmatched:10.4f} "
            f"{x.median_matched:10.4f} {x.delta_matched_minus_unmatched:+9.4f}")
    say("")
    say(f"  Sample size lost: 25 configurations per platform "
        f"(13.1%), and two of eight applications lost entirely from any")
    say("  setting that requires them on both machines.")
    say("")
    say("  Reading the deltas: where a delta is large and negative the")
    say("  unmatched result was pessimistic because of the unmatched Cirrus")
    say("  rows (its 125- and 288-core short runs), not because of the")
    say("  hardware.  Where it is near zero the published conclusion survives")
    say("  the correction.  Both are reportable, and only the matched column")
    say("  supports a statement of the form 'platform A is harder than")
    say("  platform B', because only it holds the workload fixed.")

    (OUT / "matched_crossplatform_summary.txt").write_text("\n".join(L) + "\n")
    print(f"\nwrote {OUT/'matched_crossplatform_summary.txt'}")


if __name__ == "__main__":
    main()