1#!/usr/bin/env python32"""Generate paper/tables/*.tex (LaTeX fragments), paper/tables/numbers.tex (macros) and experiments/EXPERIMENT_MANIFEST.csv3from results/processed/summary.csv, results/processed/boundary_summary.csv, results/raw/ablation_*.csv, results/raw/skew_*.csv4and results/raw/benchmarks.jsonl. Degrades gracefully when a result set is absent (the fragment then says so)."""5import glob, json, pathlib, sys6import numpy as np, pandas as pd7ROOT = pathlib.Path(__file__).resolve().parents[1]8PROC, RAW, TAB = ROOT / "results" / "processed", ROOT / "results" / "raw", ROOT / "paper" / "tables"9TAB.mkdir(parents=True, exist_ok=True)10LAB = {"r1_q0": "(1,0) plain average", "r1_q1": "(1,1) average + 1 Newton", "r1_q2": "(1,2) average + 2 Newton", "r2_q0": "(2,0) half-block jackknife",11 "r2_q1": "(2,1) half-block + 1 Newton", "r2_q2": "(2,2) half-block + 2 Newton", "r3_q0": "(3,0) quarter-block", "r3_q1": "(3,1) quarter-block + 1 Newton",12 "r3_q2": "(3,2) quarter-block + 2 Newton", "cs_q0": "analytic Cox--Snell, $q=0$", "cs_q1": "analytic Cox--Snell + 1 Newton", "csl1": "CSL, 1 round", "csl2": "CSL, 2 rounds",13 "aee": "AEE (Lin--Xi)", "stale_r1_q1": "(1,1) stale Hessian", "stale_r2_q1": "(2,1) stale Hessian", "stale_r1_q2": "(1,2) stale Hessian", "central": "centralized MLE"}14MAIN = ["r1_q0", "r1_q1", "r1_q2", "r2_q0", "r2_q1", "r2_q2", "r3_q0", "r3_q1", "cs_q1", "csl1", "csl2", "aee", "stale_r2_q1", "central"]151617def fmt(x, d=3):18 if x is None or (isinstance(x, float) and not np.isfinite(x)):19 return "--"20 a = abs(x)21 if a == 0:22 return "0"23 if a < 1e-3 or a >= 1e4:24 return f"{x:.2e}".replace("e-0", "e-").replace("e+0", "e+")25 return f"{x:.{d}f}"262728def table_cell(s, label, p, N, ms, fname, caption, ests=MAIN):29 d = s[(s.label == label) & (s.p == p) & (s.N == N) & (s.m.isin(ms))]30 if d.empty:31 (TAB / fname).write_text(f"% no data for {label} p={p} N={N}\n"); return False32 ms = sorted(d.m.unique())33 if len(ms) > 4: # keep at most four local sizes: the two extremes and two interior ones34 idx = sorted(set([0, len(ms) - 1, len(ms) // 3, (2 * len(ms)) // 3]))35 ms = [ms[i] for i in idx]36 lines = [r"\begin{table}[p]", r"\centering\scriptsize", r"\begin{threeparttable}", f"\\caption{{{caption}}}\\label{{tab:{fname[:-4]}}}",37 r"\begin{tabular}{@{}l" + "r" * (3 * len(ms)) + "@{}}", r"\toprule",38 "Estimator & " + " & ".join(f"\\multicolumn{{3}}{{c}}{{$m = {m}$, $K = {int(N//m):,}$}}" for m in ms) + r" \\",39 " & " + " & ".join("gap & RMSE & cov" for _ in ms) + r" \\", r"\midrule"]40 for e in ests:41 row = [LAB.get(e, e)]42 for m in ms:43 dd = d[(d.estimator == e) & (d.m == m)]44 if dd.empty:45 row += ["--", "--", "--"]; continue46 dd = dd.iloc[0]47 row += [fmt(dd["rms_gap"]), fmt(dd["rmse"]), (f"{dd['cov95']:.3f}" if "cov95" in dd and np.isfinite(dd["cov95"]) else "--")]48 lines.append(" & ".join(row) + r" \\")49 reps = int(d.reps.iloc[0]); cf = d.drop_duplicates("m")[["m", "central_fail", "valid_frac_r2", "valid_frac_r3"]].sort_values("m") if "valid_frac_r3" in d else None50 lines += [r"\bottomrule", r"\end{tabular}", r"\begin{tablenotes}\footnotesize",51 f"\\item {reps} paired replications per column. gap $=$ root mean square of $\\sqrt N\\,\\lVert\\tilde\\theta - \\widehat\\theta_N\\rVert$; RMSE $=$ root mean square of $\\sqrt N\\,\\lVert\\tilde\\theta - \\theta_0\\rVert$; cov $=$ empirical coverage of nominal 95\\% Wald intervals (observed information at the estimate; `--' when not computed). Monte Carlo standard errors of the gaps and RMSEs are in \\texttt{{results/processed/summary.csv}}; coverage standard error $\\approx {np.sqrt(0.95*0.05/reps):.3f}$."]52 if cf is not None:53 lines.append("\\item Fraction of workers with all local fits valid (orders 2 / 3): " + "; ".join(f"$m={int(r.m)}$: {r.valid_frac_r2:.3f} / {r.valid_frac_r3:.3f}" for r in cf.itertuples()) + ".")54 lines += [r"\end{tablenotes}", r"\end{threeparttable}", r"\end{table}"]55 (TAB / fname).write_text("\n".join(lines) + "\n"); return True565758def table_rates(s, fname="mc_rates.tex"):59 """Log-log slopes of rms_gap vs m at fixed N (largest N per model) for the main estimators."""60 rows = []61 for (label, p), d in s.groupby(["label", "p"]):62 N = d.N.max()63 for e in ["r1_q0", "r1_q1", "r1_q2", "r2_q0", "r2_q1", "r3_q0", "r3_q1", "cs_q1"]:64 dd = d[(d.N == N) & (d.estimator == e) & (d.rms_gap > 1e-9)].sort_values("m")65 if len(dd) >= 3:66 slope = np.polyfit(np.log(dd.m), np.log(dd.rms_gap), 1)[0]67 rows.append((label, p, N, e, slope, len(dd)))68 if not rows:69 (TAB / fname).write_text("% no data\n"); return70 th = {"r1_q0": -1, "r1_q1": -2, "r1_q2": -4, "r2_q0": -2, "r2_q1": -4, "r3_q0": -3, "r3_q1": -6, "cs_q1": -4}71 df = pd.DataFrame(rows, columns=["label", "p", "N", "estimator", "slope", "n"])72 piv = df.pivot_table(index=["label", "p", "N"], columns="estimator", values="slope")73 cols = [c for c in ["r1_q0", "r1_q1", "r1_q2", "r2_q0", "r2_q1", "r3_q0", "r3_q1", "cs_q1"] if c in piv.columns]74 lines = [r"\begin{table}[t]", r"\centering\small", r"\begin{threeparttable}", r"\caption{Empirical log-log slopes of the RMS scaled gap to $\widehat\theta_N$ against $m$ at fixed $N$ (largest $N$ of each design), with the slopes predicted by the leading bias term.}\label{tab:mc_rates}",75 r"\begin{tabular}{@{}llr" + "r" * len(cols) + "@{}}", r"\toprule", "Design & $p$ & $N$ & " + " & ".join(c.replace("_", "\\_") for c in cols) + r" \\",76 "predicted & & & " + " & ".join(str(th[c]) for c in cols) + r" \\", r"\midrule"]77 for idx, r in piv.iterrows():78 lines.append(f"{idx[0].replace('_', ' ')} & {idx[1]} & {int(idx[2]):,} & " + " & ".join(fmt(r[c], 2) for c in cols) + r" \\")79 lines += [r"\bottomrule", r"\end{tabular}", r"\begin{tablenotes}\footnotesize", r"\item Slopes are least-squares fits over the available $m$ grid; they approach the predicted exponent only where the bias term dominates the random remainder and the sample is in the asymptotic regime. The predicted slope for the Cox--Snell + 1 Newton estimator assumes an exact first-order bias formula ($-4$). Source: \texttt{results/processed/summary.csv}.", r"\end{tablenotes}", r"\end{threeparttable}", r"\end{table}"]80 (TAB / fname).write_text("\n".join(lines) + "\n")818283def table_boundary(fname="boundary.tex"):84 f = PROC / "boundary_summary.csv"85 if not f.exists():86 (TAB / fname).write_text("% no boundary data\n"); return87 b = pd.read_csv(f)88 sys.path.insert(0, str(ROOT / "experiments")); from boundary_pred import predicted_mean_gap89 b["pred_finite_m"] = [predicted_mean_gap(int(r.r), int(r.q), int(r.m), int(r.K), int(r.N)) for r in b.itertuples()]90 b.to_csv(PROC / "boundary_summary_with_prediction.csv", index=False)91 lines = [r"\begin{table}[t]", r"\centering\scriptsize", r"\begin{threeparttable}", r"\caption{Boundary study in the exponential model (Theorem~\ref{thm:exponential}(v)): mean of $\sqrt N(\theta^{(q)}_r - \widehat\theta_N)$ when $K \approx cN^{\alpha(r,q)}$, against the limit $-\theta_0(\beta_r c_{\mathrm{eff}}^r)^{2^q}$ ($q\ge1$) or $\beta_r\theta_0 c_{\mathrm{eff}}^r$ ($q=0$), with $c_{\mathrm{eff}} = K/N^{\alpha}$ after rounding $m$ to a multiple of the finest split.}\label{tab:boundary}",92 r"\begin{tabular}{@{}ccrrrrrrrr@{}}", r"\toprule", r"$(r,q)$ & $c$ & $N$ & $K$ & $m$ & reps & mean gap (s.e.) & finite-$m$ prediction & limit at $c_{\mathrm{eff}}$ & gap/limit \\", r"\midrule"]93 for (r, q, c), d in b.groupby(["r", "q", "c"]):94 for row in d.sort_values("N").itertuples():95 ratio = row.mean_gap / row.limit_c_eff if row.limit_c_eff != 0 else np.nan96 lines.append(f"({r},{q}) & {c:g} & $2^{{{int(np.log2(row.N))}}}$ & {row.K:,} & {row.m} & {row.reps} & {row.mean_gap:+.4g} ({row.se:.2g}) & {(f'{row.pred_finite_m:+.4g}' if np.isfinite(row.pred_finite_m) else '--')} & {row.limit_c_eff:+.4g} & {ratio:.2f} \\\\")97 lines += [r"\bottomrule", r"\end{tabular}", r"\begin{tablenotes}\footnotesize", r"\item Exact-Gamma simulation of block sufficient statistics; seed 20261028 plus a cell hash. The ratio converges to one only as $m\to\infty$; with $m \asymp N^{1/(2^{q+1}r)}$ the local sample sizes in the table are small and the finite-$m$ bias $\theta_0\sum_j w_js_j/(m-s_j)$ exceeds its leading term. Source: \texttt{results/processed/boundary\_summary.csv}.", r"\end{tablenotes}", r"\end{threeparttable}", r"\end{table}"]98 (TAB / fname).write_text("\n".join(lines) + "\n")99100101def table_ablation(model, fname):102 files = sorted(glob.glob(str(RAW / f"ablation_{model}_*.csv")))103 if not files:104 (TAB / fname).write_text(f"% no ablation data for {model}\n"); return None105 d = pd.concat([pd.read_csv(f) for f in files], ignore_index=True)106 d = d[np.isfinite(d.gap)]107 g = d.groupby(["N", "m", "variant", "pattern", "r", "q", "hessian", "bits"]).agg(rms=("gap", lambda x: np.sqrt(np.mean(x ** 2))), mean0=("gap0", "mean"), n=("gap", "size"), valid=("valid_frac", "mean")).reset_index()108 g.to_csv(PROC / f"ablation_{model}_summary.csv", index=False)109 # (a) patterns x q at the two largest m of the largest N110 N = g.N.max(); ms = sorted(g[g.N == N].m.unique())111 lines = [r"\begin{table}[p]", r"\centering\scriptsize", r"\begin{threeparttable}", f"\\caption{{Ablation ({model}, $N = {N:,}$): RMS of $\\sqrt N\\,\\lVert\\tilde\\theta - \\widehat\\theta_N\\rVert$ by split pattern, order $r$, Newton rounds $q$ and curvature used (exact / stale averaged local / single machine); quantized exact curvature with $b$ mantissa bits for $q = 1$.}}\\label{{tab:abl_{model}}}",112 r"\begin{tabular}{@{}llr" + "r" * len(ms) + "@{}}", r"\toprule", "Pattern & curvature & $q$ & " + " & ".join(f"$m = {m}$" for m in ms) + r" \\", r"\midrule"]113 for pat in ["1", "1-2", "1-3", "1-2-4", "1-2-3", "1-2-4-8", "1-2-3-4"]:114 for hess in ["exact", "stale", "single"]:115 for q in range(0, 4):116 if hess != "exact" and q == 0:117 continue118 vals = []119 for m in ms:120 dd = g[(g.N == N) & (g.m == m) & (g.pattern == pat) & (g.q == q) & (g.hessian == hess) & (g.bits == 64)]121 vals.append(fmt(float(dd.rms.iloc[0])) if len(dd) else "--")122 if all(v == "--" for v in vals):123 continue124 lines.append(f"({pat.replace('-', ',')}) & {hess} & {q} & " + " & ".join(vals) + r" \\")125 lines.append(r"\midrule")126 for pat in ["1", "1-2", "1-2-4"]:127 for hess, lab in [("quantizedH", "$b$-bit Hessian"), ("quantizedg", "$b$-bit score"), ("quantizedBoth", "$b$-bit both")]:128 for bits in sorted(g.bits.unique()):129 if bits == 64:130 continue131 vals = []132 for m in ms:133 dd = g[(g.N == N) & (g.m == m) & (g.pattern == pat) & (g.q == 1) & (g.hessian == hess) & (g.bits == bits)]134 vals.append(fmt(float(dd.rms.iloc[0])) if len(dd) else "--")135 if all(v == "--" for v in vals):136 continue137 lines.append(f"({pat.replace('-', ',')}) & {lab}, $b={bits}$ & 1 & " + " & ".join(vals) + r" \\")138 rc = d.groupby(["N", "m"]).rep.nunique(); reps = f"{int(rc.min())}--{int(rc.max())}" if rc.min() != rc.max() else str(int(rc.max()))139 lines += [r"\bottomrule", r"\end{tabular}", r"\begin{tablenotes}\footnotesize", f"\\item {reps} paired replications per column (exact count per cell in the processed CSV); $K = N/m$. Rows are omitted when the pattern does not divide $m$. Source: \\texttt{{results/processed/ablation\\_{model}\\_summary.csv}}.", r"\end{tablenotes}", r"\end{threeparttable}", r"\end{table}"]140 (TAB / fname).write_text("\n".join(lines) + "\n"); return g141142143def table_skew(fname="skew.tex"):144 files = sorted(glob.glob(str(RAW / "skew_*.csv")))145 if not files:146 (TAB / fname).write_text("% no skew data\n"); return147 d = pd.concat([pd.read_csv(f) for f in files], ignore_index=True)148 g = d.groupby(["model", "K", "ratio", "agg", "pilot"]).agg(gap0=("gap_q0", lambda x: np.sqrt(np.mean(x ** 2))), gap1=("gap_q1", lambda x: np.sqrt(np.mean(x ** 2))), gap2=("gap_q2", lambda x: np.sqrt(np.mean(x ** 2))), n=("rep", "size"), N=("N", "first")).reset_index()149 g.to_csv(PROC / "skew_summary.csv", index=False)150 lines = [r"\begin{table}[t]", r"\centering\scriptsize", r"\begin{threeparttable}", r"\caption{Unequal block sizes (Proposition~\ref{prop:robust}): RMS of $\sqrt N\,\lVert\tilde\theta - \widehat\theta_N\rVert$ for size-weighted and equal-weight aggregation of the raw and half-block-jackknife local estimates, with $q = 0, 1, 2$ Newton corrections. Sizes grow geometrically with max/min ratio as shown; $N \approx 2^{18}$.}\label{tab:skew}",151 r"\begin{tabular}{@{}llrrllrrr@{}}", r"\toprule", r"Model & $K$ & ratio & reps & weights & pilot & $q=0$ & $q=1$ & $q=2$ \\", r"\midrule"]152 for r in g.sort_values(["model", "K", "ratio", "agg", "pilot"]).itertuples():153 lines.append(f"{r.model} & {r.K} & {r.ratio:g} & {r.n} & {r.agg} & {r.pilot} & {fmt(r.gap0)} & {fmt(r.gap1)} & {fmt(r.gap2)} \\\\")154 lines += [r"\bottomrule", r"\end{tabular}", r"\begin{tablenotes}\footnotesize", r"\item Source: \texttt{results/processed/skew\_summary.csv}.", r"\end{tablenotes}", r"\end{threeparttable}", r"\end{table}"]155 (TAB / fname).write_text("\n".join(lines) + "\n")156157158def table_bench(fname="bench.tex"):159 f = RAW / "benchmarks.jsonl"160 if not f.exists():161 (TAB / fname).write_text("% no benchmark data\n"); return None162 rows = [json.loads(l) for l in open(f) if l.strip()]163 recs = []164 for r in rows:165 st = r["stats"]; A = st["stages"]["A_B_local_fit_and_upload"]; rounds = [v for k, v in st["stages"].items() if k.startswith("D_F")]166 recs.append({"label": r["label"], "K": r["K"], "hosts": r["n_hosts"], "p": r["p"], "m": r["m"], "r": len(r["splits"]), "q": r["q"], "hessian": r["hessian_mode"],167 "t_total": r["t_total_after_connect"], "t_local": A["wall"], "fit_mean": A["worker_fit_s_mean"], "fit_max": A["worker_fit_s_max"], "lat99_A": A["latency_p99"],168 "t_rounds": sum(v["wall"] for v in rounds), "round_wall": np.mean([v["wall"] for v in rounds]) if rounds else 0.0, "der_max": max([v["worker_der_s_max"] for v in rounds]) if rounds else 0.0,169 "lat99_round": max([v["latency_p99"] for v in rounds]) if rounds else 0.0, "up_B": st["bytes_up"], "down_B": st["bytes_down"], "msgs": st["messages_up"] + st["messages_down"], "grad_final": r["newton"][-1]["grad_norm"] if r["newton"] else np.nan, "t_connect": r["t_connect"]})170 d = pd.DataFrame(recs)171 g = d.groupby(["label", "K", "hosts", "p", "r", "q", "hessian"]).agg(reps=("t_total", "size"), t_total=("t_total", "median"), t_local=("t_local", "median"), fit_mean=("fit_mean", "median"), fit_max=("fit_max", "median"),172 t_rounds=("t_rounds", "median"), round_wall=("round_wall", "median"), der_max=("der_max", "median"), lat99=("lat99_round", "median"), up_B=("up_B", "first"), down_B=("down_B", "first"), msgs=("msgs", "first"), grad=("grad_final", "median"), t_connect=("t_connect", "median")).reset_index()173 g.to_csv(PROC / "benchmarks_summary.csv", index=False)174 lines = [r"\begin{table}[p]", r"\centering\scriptsize", r"\begin{threeparttable}", r"\caption{Measured multi-node benchmarks of the protocol (logistic regression, $N = 2^{22}$, medians over repeats): wall-clock of the local phase (slowest worker) and of the synchronized derivative phases, worker compute times, 99th-percentile reply latency, and bytes moved.}\label{tab:bench}",175 r"\begin{tabular}{@{}lrrrrrlrrrrrrrr@{}}", r"\toprule", r"Design & $K$ & hosts & $p$ & $r$ & $q$ & curvature & local (s) & fit max (s) & rounds (s) & per round (s) & der.\ max (s) & lat.\ p99 (s) & up (kB) & down (kB) \\", r"\midrule"]176 for r in g.itertuples():177 lines.append(f"{r.label.replace('_', ' ')} & {r.K} & {r.hosts} & {r.p} & {r.r} & {r.q} & {r.hessian} & {r.t_local:.3f} & {r.fit_max:.3f} & {r.t_rounds:.3f} & {r.round_wall:.3f} & {r.der_max:.3f} & {r.lat99:.3f} & {r.up_B/1e3:.1f} & {r.down_B/1e3:.1f} \\\\")178 lines += [r"\bottomrule", r"\end{tabular}", r"\begin{tablenotes}\footnotesize", r"\item Python 3.9/3.14 + NumPy over TCP (length-prefixed pickle frames) on the MacLustr LAN (1 GbE/10 GbE mixed); coordinator on a Mac Studio M4 Max. `local' is the wall-clock from the first `fit' broadcast to the last upload, i.e.\ the slowest worker; `rounds' sums the synchronized derivative phases. Bytes include frame headers and pickle overhead. Repeats and raw logs: \texttt{results/raw/benchmarks.jsonl}.", r"\end{tablenotes}", r"\end{threeparttable}", r"\end{table}"]179 (TAB / fname).write_text("\n".join(lines) + "\n"); return g180181182def manifest():183 rows = []184 for f in sorted(glob.glob(str(RAW / "*.manifest.json"))):185 m = json.load(open(f)); c = m["cell"]186 rows.append({"run_id": m["run_id"], "git_commit": m["git_commit"], "python_version": m["python"], "numpy": m["numpy"], "dependency_lock_sha256": "", "seed": m["seed"], "N": c["N"], "K": c["K"], "m": c["m"], "p": c["p"],187 "DGP": json.dumps(c["dgp"]), "method": "estimate_all", "optimizer": "batched Newton, step halving, tol 1e-9 (central 1e-12)", "tolerance": "1e-9", "replicates": m["reps"], "parallel_backend": "xargs -P / nohup", "node_count": 1,188 "hardware": m["host"] + " " + m["machine"], "OS": "macOS", "wall_seconds": round(m["wall_seconds"], 1), "communicated_bytes": 0, "rounds": "in-memory (no network)", "status": m["status"], "raw_output_sha256": m["raw_output_sha256"]})189 pd.DataFrame(rows).to_csv(ROOT / "experiments" / "EXPERIMENT_MANIFEST.csv", index=False)190 return len(rows)191192193def main():194 s = pd.read_csv(PROC / "summary.csv") if (PROC / "summary.csv").exists() else pd.DataFrame()195 made = {}196 if not s.empty:197 for (label, p), d in s.groupby(["label", "p"]):198 N = d.N.max()199 made[(label, p)] = table_cell(s, label, p, N, sorted(d.m.unique()), f"mc_{label}_p{p}.tex", f"Monte Carlo, {label.replace('_', ' ')} design, $p = {p}$, $N = {N:,}$: scaled gap to $\\widehat\\theta_N$, scaled RMSE and 95\\% coverage by local sample size $m$.")200 table_rates(s)201 table_boundary(); table_ablation("exponential", "abl_exponential.tex"); table_ablation("logistic", "abl_logistic.tex"); table_skew(); table_bench()202 n = manifest()203 # numbers macros204 macros = {"numMCcells": str(s.groupby(["label", "p", "N", "K"]).ngroups if not s.empty else 0), "numManifest": str(n)}205 if not s.empty:206 macros["numMCrepsExp"] = str(int(s[s.label == "exp"].reps.max())) if (s.label == "exp").any() else "0"207 macros["numMCrepsLogit"] = str(int(s[s.label == "logit"].reps.max())) if (s.label == "logit").any() else "0"208 bf = PROC / "benchmarks_summary.csv"209 macros["benchHosts"] = str(int(pd.read_csv(bf).hosts.max())) if bf.exists() else "0"210 with open(TAB / "numbers.tex", "w") as f:211 f.write("% Auto-generated by experiments/make_tables.py\n")212 for k, v in macros.items():213 f.write(f"\\newcommand{{\\{k}}}{{{v}}}\n")214 print("tables written:", sorted(p.name for p in TAB.glob("*.tex")), "manifest rows:", n)215216217if __name__ == "__main__":218 main()219