"""Analysis of the ECDLP studies: statistics with bootstrap CIs, empirical complexity fits, figures and LaTeX tables.

Inputs : results/synthetic_ecdlp/{scaling_generic,scaling_j0,sweep_r,parallel,tmto}_<host>.json
Outputs: results/synthetic_ecdlp/analysis.json, figures/ecdlp_*.pdf (+ copies in paper/figures), paper/tables/ecdlp_*.tex
"""
import glob, json, math, pathlib, re, sys
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
ROOT = pathlib.Path(__file__).resolve().parents[2]
R = ROOT / "results/synthetic_ecdlp"; F = ROOT / "figures"; PF = ROOT / "paper/figures"; T = ROOT / "paper/tables"
for d in (F, PF, T):
    d.mkdir(parents=True, exist_ok=True)
plt.rcParams.update({"font.family": "serif", "font.size": 9, "mathtext.fontset": "stix", "axes.grid": True, "grid.alpha": 0.3})
rng = np.random.default_rng(34)
ANALYSIS = {}
MACROS = {}


def load(study):
    files = sorted(glob.glob(str(R / f"{study}_*.json")))
    files = [f for f in files if ".chunk" not in f and ".raw" not in f]
    if not files:
        return None, None
    d = json.loads(pathlib.Path(files[0]).read_text())
    runs = [r for r in d["measurements"]["runs"] if r]
    for r in runs:
        r["config"] = r["label"].split("/")[1]; r["inst"] = r["label"].split("/")[0]
        r["bits"] = int(re.match(r"[gj](\d+)b", r["inst"]).group(1))
        r["n"] = int(r["params"]["instance"]["n"]) if "instance" in r["params"] else int(r["params"]["n"])
    return d, runs


def boot_mean(x, B=2000):
    x = np.asarray(x, float)
    if len(x) < 2:
        return float(x.mean()), float(x.mean()), float(x.mean())
    means = rng.choice(x, size=(B, len(x)), replace=True).mean(axis=1)
    return float(x.mean()), float(np.percentile(means, 2.5)), float(np.percentile(means, 97.5))


def boot_ratio(a, b, B=2000):
    """Bootstrap CI for mean(a)/mean(b) with paired resampling when lengths match."""
    a, b = np.asarray(a, float), np.asarray(b, float)
    if len(a) == len(b):
        idx = rng.integers(0, len(a), size=(B, len(a)))
        r = a[idx].mean(axis=1) / b[idx].mean(axis=1)
    else:
        r = rng.choice(a, size=(B, len(a))).mean(axis=1) / rng.choice(b, size=(B, len(b))).mean(axis=1)
    return float(a.mean() / b.mean()), float(np.percentile(r, 2.5)), float(np.percentile(r, 97.5))


def fit_power(ns, ys, with_log=False):
    """Least squares on log y = log a + α log n (+ β log log2 n). Returns (a, α, β, r²) and bootstrap CI for α."""
    ns, ys = np.asarray(ns, float), np.asarray(ys, float)
    X = np.column_stack([np.ones_like(ns), np.log(ns)] + ([np.log(np.log2(ns))] if with_log else []))
    y = np.log(ys)
    coef, *_ = np.linalg.lstsq(X, y, rcond=None)
    pred = X @ coef; r2 = 1 - ((y - pred) ** 2).sum() / ((y - y.mean()) ** 2).sum()
    alphas = []
    for _ in range(1000):
        idx = rng.integers(0, len(ns), len(ns))
        c, *_ = np.linalg.lstsq(X[idx], y[idx], rcond=None); alphas.append(c[1])
    return {"a": float(math.exp(coef[0])), "alpha": float(coef[1]), "beta": float(coef[2]) if with_log else 0.0, "r2": float(r2),
            "alpha_ci": [float(np.percentile(alphas, 2.5)), float(np.percentile(alphas, 97.5))]}


def summarise(runs, key="group_ops"):
    out = {}
    for cfg in sorted({r["config"] for r in runs}):
        out[cfg] = {}
        for bits in sorted({r["bits"] for r in runs if r["config"] == cfg}):
            rs = [r for r in runs if r["config"] == cfg and r["bits"] == bits]
            n = float(np.mean([r["n"] for r in rs]))
            ops = [r[key] for r in rs]; mean, lo, hi = boot_mean(ops)
            its = [r.get("iterations", 0) for r in rs]
            out[cfg][bits] = {"n_mean": n, "count": len(rs), "ops_mean": mean, "ops_ci": [lo, hi], "ops_median": float(np.median(ops)),
                              "ops_sd": float(np.std(ops)), "ops_over_sqrt_n": mean / math.sqrt(n), "iter_mean": float(np.mean(its)),
                              "iter_over_sqrt_n": float(np.mean(its)) / math.sqrt(n), "wall_mean": float(np.mean([r["wall_s"] for r in rs])),
                              "cpu_mean": float(np.mean([r["cpu_s"] for r in rs])), "rss_max_mb": max(r["peak_rss_bytes"] for r in rs) / 2**20,
                              "table_entries_mean": float(np.mean([r.get("table_entries", 0) or 0 for r in rs])),
                              "dp_count_mean": float(np.mean([r.get("dp_count", 0) or 0 for r in rs])),
                              "fruitless_mean": float(np.mean([r.get("fruitless_cycles", 0) or 0 for r in rs])),
                              "ok": all(r["ok"] for r in rs)}
    return out


LABELS = {"bsgs": "BSGS", "rho_floyd": "rho, Floyd", "rho_brent": "rho, Brent", "rho_dp": "rho, DP ($r=20$)", "rho_dp_r128": "rho, DP ($r=128$)",
          "rho_dp_neg": "rho, DP + negation ($r=128$)", "kangaroo": "kangaroo (full interval)",
          "rho_dp_aut3": r"rho, DP + $\phi$-orbits", "rho_dp_neg_aut3": r"rho, DP + negation + $\phi$ ($r=128$)", "rho_dp_neg_aut3_r20": r"rho, DP + negation + $\phi$ ($r=20$)"}
THEORY = {"bsgs": lambda n: 2 * math.sqrt(n), "rho_floyd": lambda n: 3 * math.sqrt(math.pi * n / 2), "rho_brent": lambda n: 1.0 * math.sqrt(math.pi * n / 2) * 1.15,
          "rho_dp": lambda n: math.sqrt(math.pi * n / 2), "rho_dp_r128": lambda n: math.sqrt(math.pi * n / 2), "rho_dp_neg": lambda n: math.sqrt(math.pi * n / 4),
          "kangaroo": lambda n: 2 * math.sqrt(n), "rho_dp_aut3": lambda n: math.sqrt(math.pi * n / 6), "rho_dp_neg_aut3": lambda n: math.sqrt(math.pi * n / 12),
          "rho_dp_neg_aut3_r20": lambda n: math.sqrt(math.pi * n / 12)}


def analyse_scaling():
    d, runs = load("scaling_generic")
    if not runs:
        return
    S = summarise(runs); ANALYSIS["scaling_generic"] = {"host": d["host"], "summary": S}
    # fits per config on group_ops
    fits = {}
    for cfg, per in S.items():
        ns = [v["n_mean"] for v in per.values()]; ys = [v["ops_mean"] for v in per.values()]
        if len(ns) >= 3:
            big = [(v["n_mean"], v["ops_mean"], v["iter_mean"]) for b, v in per.items() if b >= 32]
            fits[cfg] = {"power": fit_power(ns, ys), "power_log": fit_power(ns, ys, True),
                         "power_ge32": fit_power([t[0] for t in big], [t[1] for t in big]) if len(big) >= 3 else None,
                         "power_iter_ge32": fit_power([t[0] for t in big], [t[2] for t in big]) if len(big) >= 3 and all(t[2] > 0 for t in big) else None}
    ANALYSIS["scaling_generic"]["fits"] = fits
    # figure: runtime scaling (ops and wall) log-log with theory lines
    fig, axes = plt.subplots(1, 2, figsize=(7.2, 3.3))
    for cfg, per in S.items():
        bits = sorted(per); ns = [per[b]["n_mean"] for b in bits]
        axes[0].errorbar(ns, [per[b]["ops_mean"] for b in bits], yerr=[[per[b]["ops_mean"] - per[b]["ops_ci"][0] for b in bits], [per[b]["ops_ci"][1] - per[b]["ops_mean"] for b in bits]], marker="o", ms=3, lw=1, capsize=2, label=LABELS.get(cfg, cfg))
        axes[1].plot(ns, [per[b]["wall_mean"] for b in bits], marker="o", ms=3, lw=1, label=LABELS.get(cfg, cfg))
    nn = np.logspace(math.log10(2**20), math.log10(2**56), 50)
    axes[0].plot(nn, np.sqrt(np.pi * nn / 2), "k--", lw=0.8, label=r"$\sqrt{\pi n/2}$")
    axes[0].plot(nn, 2 * np.sqrt(nn), "k:", lw=0.8, label=r"$2\sqrt{n}$")
    for ax in axes:
        ax.set_xscale("log", base=2); ax.set_yscale("log"); ax.set_xlabel("group order $n$")
    axes[0].set_ylabel("group operations (mean, 95% CI)"); axes[1].set_ylabel("wall-clock seconds (1 thread)")
    axes[0].legend(fontsize=6, ncol=1, loc="lower right"); fig.tight_layout()
    fig.savefig(F / "ecdlp_scaling.pdf"); fig.savefig(PF / "ecdlp_scaling.pdf"); plt.close(fig)
    # figure: normalised ops/√n with theory constants
    fig, ax = plt.subplots(figsize=(6.4, 3.2))
    for cfg, per in S.items():
        bits = sorted(per)
        ax.plot(bits, [per[b]["ops_over_sqrt_n"] for b in bits], marker="o", ms=3, lw=1, label=LABELS.get(cfg, cfg))
    for c, lab in [(2, r"$2$ (BSGS, kangaroo)"), (math.sqrt(math.pi / 2), r"$\sqrt{\pi/2}$ (rho)"), (math.sqrt(math.pi / 4), r"$\sqrt{\pi/4}$ (negation)")]:
        ax.axhline(c, color="gray", lw=0.6, ls="--"); ax.text(bits[-1] + 0.3, c, lab, fontsize=6, va="center")
    ax.set_xlabel("group size (bits)"); ax.set_ylabel(r"group operations / $\sqrt{n}$"); ax.set_yscale("log"); ax.legend(fontsize=6, ncol=2)
    fig.tight_layout(); fig.savefig(F / "ecdlp_constants.pdf"); fig.savefig(PF / "ecdlp_constants.pdf"); plt.close(fig)
    # table
    rows = [r"\begin{tabular}{@{}llrrrrrr@{}}", r"\toprule Algorithm & Bits & Inst. & Group ops (mean) & 95\% CI & ops$/\sqrt n$ & Theory & Wall (s) \\ \midrule"]
    for cfg, per in S.items():
        for i, (b, v) in enumerate(sorted(per.items())):
            if b not in (20, 28, 36, 44, 48, 52, 56):
                continue
            th = THEORY[cfg](v["n_mean"]) / math.sqrt(v["n_mean"])
            rows.append(f"{LABELS.get(cfg, cfg) if i == 0 else ''} & {b} & {v['count']} & {v['ops_mean']:,.0f} & [{v['ops_ci'][0]:,.0f}, {v['ops_ci'][1]:,.0f}] & {v['ops_over_sqrt_n']:.2f} & {th:.2f} & {v['wall_mean']:.3f} \\\\".replace(",", r"\,"))
    rows += [r"\bottomrule", r"\end{tabular}"]
    (T / "ecdlp_scaling.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")
    rows = [r"\begin{tabular}{@{}lrrrrrr@{}}", r"\toprule Algorithm & $\hat\alpha$ (all sizes, ops) & 95\% CI & $R^2$ & $\hat\alpha$ ($\ge 32$ bits, ops) & 95\% CI & $\hat\alpha$ ($\ge 32$ bits, iterations) \\ \midrule"]
    for cfg, ff in fits.items():
        p, pg, pi = ff["power"], ff["power_ge32"], ff["power_iter_ge32"]
        g = f"{pg['alpha']:.3f} & [{pg['alpha_ci'][0]:.3f}, {pg['alpha_ci'][1]:.3f}]" if pg else "--- & ---"
        it = f"{pi['alpha']:.3f} [{pi['alpha_ci'][0]:.3f}, {pi['alpha_ci'][1]:.3f}]" if pi else "---"
        rows.append(f"{LABELS.get(cfg, cfg)} & {p['alpha']:.3f} & [{p['alpha_ci'][0]:.3f}, {p['alpha_ci'][1]:.3f}] & {p['r2']:.4f} & {g} & {it} \\\\")
    rows += [r"\bottomrule", r"\end{tabular}"]
    (T / "ecdlp_fits.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")
    for cfg, ff in fits.items():
        k = cfg.replace('_', '').replace('128', 'big')
        MACROS[f"alpha{k}"] = f"{ff['power']['alpha']:.3f}"
        MACROS[f"alphaLo{k}"] = f"{ff['power']['alpha_ci'][0]:.3f}"; MACROS[f"alphaHi{k}"] = f"{ff['power']['alpha_ci'][1]:.3f}"
        if ff["power_ge32"]:
            MACROS[f"alphaG{k}"] = f"{ff['power_ge32']['alpha']:.3f}"; MACROS[f"alphaGLo{k}"] = f"{ff['power_ge32']['alpha_ci'][0]:.3f}"; MACROS[f"alphaGHi{k}"] = f"{ff['power_ge32']['alpha_ci'][1]:.3f}"
    alphas = [ff["power_ge32"]["alpha"] for ff in fits.values() if ff["power_ge32"]]
    MACROS["alphaGMin"] = f"{min(alphas):.3f}"; MACROS["alphaGMax"] = f"{max(alphas):.3f}"
    big = max(S["rho_dp"]); MACROS["scalingMaxBits"] = str(big); MACROS["scalingInstances"] = str(S["rho_dp"][big]["count"])
    MACROS["rhoDpConstBig"] = f"{S['rho_dp'][big]['ops_over_sqrt_n']:.2f}"; MACROS["rhoDpIterConstBig"] = f"{S['rho_dp'][big]['iter_over_sqrt_n']:.2f}"
    MACROS["rhoNegIterConstBig"] = f"{S['rho_dp_neg'][big]['iter_over_sqrt_n']:.2f}"
    MACROS["bsgsConstBig"] = f"{S['bsgs'][max(S['bsgs'])]['ops_over_sqrt_n']:.2f}"; MACROS["kangarooConstBig"] = f"{S['kangaroo'][max(S['kangaroo'])]['ops_over_sqrt_n']:.2f}"
    MACROS["floydConstBig"] = f"{S['rho_floyd'][max(S['rho_floyd'])]['ops_over_sqrt_n']:.2f}"; MACROS["brentConstBig"] = f"{S['rho_brent'][max(S['rho_brent'])]['ops_over_sqrt_n']:.2f}"
    MACROS["scalingAllOk"] = "all" if all(v["ok"] for per in S.values() for v in per.values()) else "NOT ALL"
    MACROS["scalingHost"] = d["host"]["cpu"]
    # speedup neg vs plain (paired by instance)
    sp = {}
    for b in sorted(S["rho_dp_neg"]):
        a = [r["iterations"] for r in runs if r["config"] == "rho_dp_r128" and r["bits"] == b]
        c = [r["iterations"] for r in runs if r["config"] == "rho_dp_neg" and r["bits"] == b]
        sp[b] = boot_ratio(a, c)
    ANALYSIS["scaling_generic"]["neg_speedup_iterations"] = sp
    MACROS["negSpeedupBig"] = f"{sp[big][0]:.2f}"; MACROS["negSpeedupBigLo"] = f"{sp[big][1]:.2f}"; MACROS["negSpeedupBigHi"] = f"{sp[big][2]:.2f}"
    PB = 44
    allsp = [r["iterations"] for r in runs if r["config"] == "rho_dp_r128" and r["bits"] >= PB]; allc = [r["iterations"] for r in runs if r["config"] == "rho_dp_neg" and r["bits"] >= PB]
    pooled = boot_ratio(np.array(allsp) / np.array([math.sqrt(r["n"]) for r in runs if r["config"] == "rho_dp_r128" and r["bits"] >= PB]),
                        np.array(allc) / np.array([math.sqrt(r["n"]) for r in runs if r["config"] == "rho_dp_neg" and r["bits"] >= PB]))
    MACROS["pooledMinBits"] = str(PB); MACROS["pooledN"] = str(len(allsp))
    # plain rho constant pooled over ≥ PB bits, with CI, against √(π/2)
    c = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] in ("rho_dp", "rho_dp_r128") and r["bits"] >= PB]
    m_, lo_, hi_ = boot_mean(c); MACROS["rhoConstPooled"] = f"{m_:.3f}"; MACROS["rhoConstPooledLo"] = f"{lo_:.3f}"; MACROS["rhoConstPooledHi"] = f"{hi_:.3f}"; MACROS["rhoConstPooledN"] = str(len(c))
    c = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == "rho_dp_neg" and r["bits"] >= PB]
    m_, lo_, hi_ = boot_mean(c); MACROS["negConstPooled"] = f"{m_:.3f}"; MACROS["negConstPooledLo"] = f"{lo_:.3f}"; MACROS["negConstPooledHi"] = f"{hi_:.3f}"
    c = [r["group_ops"] / math.sqrt(r["n"]) for r in runs if r["config"] == "bsgs" and r["bits"] >= 36]
    m_, lo_, hi_ = boot_mean(c); MACROS["bsgsConstPooled"] = f"{m_:.3f}"; MACROS["bsgsConstPooledLo"] = f"{lo_:.3f}"; MACROS["bsgsConstPooledHi"] = f"{hi_:.3f}"
    c = [r["group_ops"] / math.sqrt(r["n"]) for r in runs if r["config"] == "kangaroo" and r["bits"] >= PB]
    m_, lo_, hi_ = boot_mean(c); MACROS["kangConstPooled"] = f"{m_:.3f}"; MACROS["kangConstPooledLo"] = f"{lo_:.3f}"; MACROS["kangConstPooledHi"] = f"{hi_:.3f}"
    ANALYSIS["scaling_generic"]["pooled_constants"] = {k: v for k, v in MACROS.items() if "Pooled" in k}
    MACROS["negSpeedupPooled"] = f"{pooled[0]:.2f}"; MACROS["negSpeedupPooledLo"] = f"{pooled[1]:.2f}"; MACROS["negSpeedupPooledHi"] = f"{pooled[2]:.2f}"
    ANALYSIS["scaling_generic"]["neg_speedup_pooled_ge36"] = pooled


def analyse_j0():
    d, runs = load("scaling_j0")
    if not runs:
        return
    S = summarise(runs); ANALYSIS["scaling_j0"] = {"host": d["host"], "summary": S}
    base = "rho_dp"
    ratios = {}
    for cfg in ("rho_dp_neg", "rho_dp_aut3", "rho_dp_neg_aut3", "rho_dp_neg_aut3_r20"):
        ratios[cfg] = {}
        for b in sorted(S[cfg]):
            a = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == base and r["bits"] == b]
            c = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == cfg and r["bits"] == b]
            ratios[cfg][b] = boot_ratio(a, c)
        PB = 44
        a = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == base and r["bits"] >= PB]
        c = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == cfg and r["bits"] >= PB]
        ratios[cfg]["pooled_ge36"] = boot_ratio(a, c)  # key kept for compatibility; pooled over bits ≥ PB
        ratios[cfg]["pooled_min_bits"] = PB; ratios[cfg]["pooled_n"] = len(a)
        m_, lo_, hi_ = boot_mean(c); ratios[cfg]["const_pooled"] = [m_, lo_, hi_]
        ratios[cfg]["cpu_ns_per_iter_ge40"] = float(np.mean([1e9 * r["cpu_s"] / max(r["iterations"], 1) for r in runs if r["config"] == cfg and r["bits"] >= 40]))
        # wall-clock ratio (what the user actually gains)
        a = [r["wall_s"] for r in runs if r["config"] == base and r["bits"] >= 40]; c = [r["wall_s"] for r in runs if r["config"] == cfg and r["bits"] >= 40]
        ratios[cfg]["wall_pooled_ge40"] = boot_ratio(a, c)
        a = [r["group_ops"] / math.sqrt(r["n"]) for r in runs if r["config"] == base and r["bits"] >= PB]; c = [r["group_ops"] / math.sqrt(r["n"]) for r in runs if r["config"] == cfg and r["bits"] >= PB]
        ratios[cfg]["ops_pooled_ge36"] = boot_ratio(a, c)
    a = [r["iterations"] / math.sqrt(r["n"]) for r in runs if r["config"] == base and r["bits"] >= PB]
    m_, lo_, hi_ = boot_mean(a); ratios["_plain_const_pooled"] = [m_, lo_, hi_]
    ratios["_plain_cpu_ns_per_iter_ge40"] = float(np.mean([1e9 * r["cpu_s"] / max(r["iterations"], 1) for r in runs if r["config"] == base and r["bits"] >= 40]))
    ANALYSIS["scaling_j0"]["iteration_speedup_vs_plain"] = ratios
    fig, axes = plt.subplots(1, 2, figsize=(7.2, 3.2))
    for cfg, lab, th in [("rho_dp_neg", "negation", math.sqrt(2)), ("rho_dp_aut3", r"$\phi$-orbits", math.sqrt(3)), ("rho_dp_neg_aut3", r"negation + $\phi$ ($r=128$)", math.sqrt(6)), ("rho_dp_neg_aut3_r20", r"negation + $\phi$ ($r=20$)", math.sqrt(6))]:
        bits = sorted(b for b in ratios[cfg] if isinstance(b, int))
        m = [ratios[cfg][b][0] for b in bits]; lo = [ratios[cfg][b][0] - ratios[cfg][b][1] for b in bits]; hi = [ratios[cfg][b][2] - ratios[cfg][b][0] for b in bits]
        l = axes[0].errorbar(bits, m, yerr=[lo, hi], marker="o", ms=3, lw=1, capsize=2, label=lab)
        axes[0].axhline(th, color=l[0].get_color(), lw=0.6, ls="--")
    axes[0].set_xlabel("group size (bits)"); axes[0].set_ylabel("iteration speed-up vs plain rho (paired, 95% CI)"); axes[0].legend(fontsize=6)
    axes[0].set_title(r"dashed: $\sqrt{2}$, $\sqrt{3}$, $\sqrt{6}$", fontsize=8)
    for cfg in ("rho_dp_neg", "rho_dp_aut3", "rho_dp_neg_aut3", "rho_dp_neg_aut3_r20"):
        bits = sorted(S[cfg]); axes[1].plot(bits, [S[cfg][b]["fruitless_mean"] / max(S[cfg][b]["iter_mean"], 1) * 1e3 for b in bits], marker="o", ms=3, lw=1, label=LABELS.get(cfg, cfg))
    axes[1].set_xlabel("group size (bits)"); axes[1].set_ylabel("fruitless cycles per 1000 iterations"); axes[1].legend(fontsize=6)
    fig.tight_layout(); fig.savefig(F / "ecdlp_endomorphism.pdf"); fig.savefig(PF / "ecdlp_endomorphism.pdf"); plt.close(fig)
    rows = [r"\begin{tabular}{@{}lrrrrr@{}}", r"\toprule Classes & Theory & Iterations$/\sqrt n$ (95\% CI) & Iteration speed-up (paired, 95\% CI) & ns per iteration & Wall-clock speed-up \\ \midrule",
            f"none (plain rho, $r=128$) & $\\sqrt{{\\pi/2}} = 1.253$ & {ratios['_plain_const_pooled'][0]:.3f} [{ratios['_plain_const_pooled'][1]:.3f}, {ratios['_plain_const_pooled'][2]:.3f}] & 1 & {ratios['_plain_cpu_ns_per_iter_ge40']:.0f} & 1 \\\\"]
    for cfg, lab, th in [("rho_dp_neg", r"$\{\pm P\}$ (negation)", r"$\sqrt2 = 1.414$"), ("rho_dp_aut3", r"$\{P,\phi P,\phi^2P\}$", r"$\sqrt3 = 1.732$"), ("rho_dp_neg_aut3", r"order-6 orbits, $r=128$", r"$\sqrt6 = 2.449$"), ("rho_dp_neg_aut3_r20", r"order-6 orbits, $r=20$", r"$\sqrt6 = 2.449$")]:
        p, o, w, cp = ratios[cfg]["pooled_ge36"], ratios[cfg]["ops_pooled_ge36"], ratios[cfg]["wall_pooled_ge40"], ratios[cfg]["const_pooled"]
        rows.append(f"{lab} & {th} & {cp[0]:.3f} [{cp[1]:.3f}, {cp[2]:.3f}] & {p[0]:.3f} [{p[1]:.3f}, {p[2]:.3f}] & {ratios[cfg]['cpu_ns_per_iter_ge40']:.0f} & {w[0]:.2f} [{w[1]:.2f}, {w[2]:.2f}] \\\\")
    rows += [r"\bottomrule", r"\end{tabular}"]
    (T / "ecdlp_j0.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")
    for cfg, key in [("rho_dp_neg", "Neg"), ("rho_dp_aut3", "Aut"), ("rho_dp_neg_aut3", "NegAut"), ("rho_dp_neg_aut3_r20", "NegAutRtwenty")]:
        p = ratios[cfg]["pooled_ge36"]; MACROS[f"jzero{key}Speedup"] = f"{p[0]:.2f}"; MACROS[f"jzero{key}SpeedupLo"] = f"{p[1]:.2f}"; MACROS[f"jzero{key}SpeedupHi"] = f"{p[2]:.2f}"
        w = ratios[cfg]["wall_pooled_ge40"]; MACROS[f"jzero{key}Wall"] = f"{w[0]:.2f}"; MACROS[f"jzero{key}WallLo"] = f"{w[1]:.2f}"; MACROS[f"jzero{key}WallHi"] = f"{w[2]:.2f}"
        MACROS[f"jzero{key}NsIter"] = f"{ratios[cfg]['cpu_ns_per_iter_ge40']:.0f}"; MACROS[f"jzero{key}Const"] = f"{ratios[cfg]['const_pooled'][0]:.3f}"
    MACROS["jzeroPlainNsIter"] = f"{ratios['_plain_cpu_ns_per_iter_ge40']:.0f}"; MACROS["jzeroPlainConst"] = f"{ratios['_plain_const_pooled'][0]:.3f}"
    MACROS["jzeroPooledN"] = str(ratios["rho_dp_neg"]["pooled_n"]); MACROS["jzeroPooledMinBits"] = str(ratios["rho_dp_neg"]["pooled_min_bits"])
    MACROS["jzeroMaxBits"] = str(max(S["rho_dp"])); MACROS["jzeroHost"] = d["host"]["cpu"]
    MACROS["jzeroAllOk"] = "all" if all(v["ok"] for per in S.values() for v in per.values()) else "NOT ALL"
    b = max(S["rho_dp_neg_aut3"]); MACROS["jzeroFruitlessPerK"] = f"{S['rho_dp_neg_aut3'][b]['fruitless_mean'] / max(S['rho_dp_neg_aut3'][b]['iter_mean'], 1) * 1e3:.2f}"
    MACROS["jzeroFruitlessPerKRtwenty"] = f"{S['rho_dp_neg_aut3_r20'][b]['fruitless_mean'] / max(S['rho_dp_neg_aut3_r20'][b]['iter_mean'], 1) * 1e3:.2f}"


def analyse_sweep():
    d, runs = load("sweep_r")
    if not runs:
        return
    res = {}
    for bits in sorted({r["bits"] for r in runs}):
        res[bits] = {}
        for cfg in sorted({r["config"] for r in runs}):
            rs = [r for r in runs if r["bits"] == bits and r["config"] == cfg]
            if not rs:
                continue
            x = [r["iterations"] / math.sqrt(r["n"]) for r in rs]; m, lo, hi = boot_mean(x)
            res[bits][cfg] = {"iter_over_sqrt_n": m, "ci": [lo, hi], "count": len(rs), "ops_over_sqrt_n": float(np.mean([r["group_ops"] / math.sqrt(r["n"]) for r in rs]))}
    ANALYSIS["sweep_r"] = {"host": d["host"], "results": res}
    fig, ax = plt.subplots(figsize=(6.4, 3.2))
    rs_ = [2, 3, 4, 8, 16, 20, 32, 64, 128]
    for bits, mk in zip(sorted(res), ("o", "s")):
        for mixed, ls in ((0, "-"), (2, "--"), (4, ":")):
            xs, ys, lo, hi = [], [], [], []
            for r in rs_:
                k = f"r{r}_m{mixed}"
                if k in res[bits]:
                    xs.append(r); ys.append(res[bits][k]["iter_over_sqrt_n"]); lo.append(ys[-1] - res[bits][k]["ci"][0]); hi.append(res[bits][k]["ci"][1] - ys[-1])
            if xs:
                ax.errorbar(xs, ys, yerr=[lo, hi], marker=mk, ms=3, lw=1, ls=ls, capsize=2, label=f"{bits} bits, mixed={mixed or 'none'}")
    ax.axhline(math.sqrt(math.pi / 2), color="k", lw=0.8, ls="--"); ax.text(2.1, math.sqrt(math.pi / 2) * 1.02, r"random mapping $\sqrt{\pi/2}$", fontsize=7)
    ax.set_xscale("log", base=2); ax.set_xlabel("number of partitions $r$ of the adding walk"); ax.set_ylabel(r"iterations / $\sqrt{n}$ (mean, 95% CI)"); ax.legend(fontsize=6, ncol=2)
    fig.tight_layout(); fig.savefig(F / "ecdlp_sweep_r.pdf"); fig.savefig(PF / "ecdlp_sweep_r.pdf"); plt.close(fig)
    rows = [r"\begin{tabular}{@{}lrrrrrr@{}}", r"\toprule $r$ & \multicolumn{3}{c}{36 bits: iterations$/\sqrt n$} & \multicolumn{3}{c}{40 bits: iterations$/\sqrt n$} \\", r" & plain & mixed 2 & mixed 4 & plain & mixed 2 & mixed 4 \\ \midrule"]
    for r in rs_:
        cells = []
        for bits in (36, 40):
            for mixed in (0, 2, 4):
                k = f"r{r}_m{mixed}"; v = res.get(bits, {}).get(k)
                cells.append(f"{v['iter_over_sqrt_n']:.3f}" if v else "---")
        rows.append(f"{r} & " + " & ".join(cells) + r" \\")
    rows += [r"\midrule", r"random-mapping constant $\sqrt{\pi/2}$ & 1.253 & & & 1.253 & & \\", r"\bottomrule", r"\end{tabular}"]
    (T / "ecdlp_sweep_r.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")
    best = min(((bits, cfg, v["iter_over_sqrt_n"], v["ci"]) for bits in res for cfg, v in res[bits].items() if cfg != "r2_m0"), key=lambda t: t[2])
    MACROS["sweepBestCfg"] = best[1].replace("_", r"\_"); MACROS["sweepBestConst"] = f"{best[2]:.3f}"; MACROS["sweepBestLo"] = f"{best[3][0]:.3f}"; MACROS["sweepBestHi"] = f"{best[3][1]:.3f}"
    MACROS["sweepRtwo"] = f"{res[40]['r2_m0']['iter_over_sqrt_n']:.2f}" if "r2_m0" in res.get(40, {}) else "n/a"
    MACROS["sweepRthree"] = f"{res[40]['r3_m0']['iter_over_sqrt_n']:.2f}" if "r3_m0" in res.get(40, {}) else "n/a"
    MACROS["sweepRtwenty"] = f"{res[40]['r20_m0']['iter_over_sqrt_n']:.3f}"; MACROS["sweepRhundred"] = f"{res[40]['r128_m0']['iter_over_sqrt_n']:.3f}"
    MACROS["sweepConfigs"] = str(len({r["config"] for r in runs})); MACROS["sweepInstances"] = str(len({r["inst"] for r in runs}))


def analyse_parallel():
    d, runs = load("parallel")
    if not runs:
        return
    res = {}
    for t in sorted({r["threads"] for r in runs}):
        rs = [r for r in runs if r["threads"] == t]
        res[t] = {"wall_mean": float(np.mean([r["wall_s"] for r in rs])), "ops_mean": float(np.mean([r["group_ops"] for r in rs])),
                  "ops_per_wall_s": float(np.mean([r["group_ops"] / r["wall_s"] for r in rs])), "count": len(rs), "cpu_mean": float(np.mean([r["cpu_s"] for r in rs]))}
    base = res[1]
    for t, v in res.items():
        v["speedup_throughput"] = v["ops_per_wall_s"] / base["ops_per_wall_s"]; v["efficiency"] = v["speedup_throughput"] / t
        v["ops_ratio_to_1thread"] = v["ops_mean"] / base["ops_mean"]
    ANALYSIS["parallel"] = {"host": d["host"], "results": res}
    fig, axes = plt.subplots(1, 2, figsize=(7.0, 3.0))
    ts = sorted(res); axes[0].plot(ts, [res[t]["speedup_throughput"] for t in ts], marker="o", ms=3, label="measured throughput speed-up"); axes[0].plot(ts, ts, "k--", lw=0.8, label="linear")
    axes[0].set_xlabel("threads"); axes[0].set_ylabel("throughput speed-up vs 1 thread"); axes[0].legend(fontsize=7)
    axes[1].plot(ts, [res[t]["ops_ratio_to_1thread"] for t in ts], marker="s", ms=3, color="C3"); axes[1].axhline(1, color="k", lw=0.8, ls="--")
    axes[1].set_xlabel("threads"); axes[1].set_ylabel("total group ops, relative to 1 thread")
    fig.tight_layout(); fig.savefig(F / "ecdlp_parallel.pdf"); fig.savefig(PF / "ecdlp_parallel.pdf"); plt.close(fig)
    rows = [r"\begin{tabular}{@{}rrrrrr@{}}", r"\toprule Threads & Wall (s) & Group ops ($\times 10^6$) & Mops/s & Speed-up & Efficiency \\ \midrule"]
    for t in ts:
        v = res[t]; rows.append(f"{t} & {v['wall_mean']:.2f} & {v['ops_mean']/1e6:.1f} & {v['ops_per_wall_s']/1e6:.0f} & {v['speedup_throughput']:.2f} & {v['efficiency']:.2f} \\\\")
    rows += [r"\bottomrule", r"\end{tabular}"]
    (T / "ecdlp_parallel.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")
    tmax = max(ts); MACROS["parMaxThreads"] = str(tmax); MACROS["parSpeedupMax"] = f"{res[tmax]['speedup_throughput']:.1f}"; MACROS["parEffMax"] = f"{res[tmax]['efficiency']:.2f}"
    t16 = 16 if 16 in res else tmax; MACROS["parSpeedupSixteen"] = f"{res[t16]['speedup_throughput']:.1f}"; MACROS["parEffSixteen"] = f"{res[t16]['efficiency']:.2f}"
    MACROS["parOpsRatioMax"] = f"{res[tmax]['ops_ratio_to_1thread']:.2f}"; MACROS["parHost"] = d["host"]["cpu"]; MACROS["parCores"] = str(d["host"]["cores"])


def analyse_tmto():
    d, runs = load("tmto")
    if not runs:
        return
    bs = {}
    for r in runs:
        if r["config"].startswith("bsgs_m"):
            m = int(r["config"][6:]); bs.setdefault(m, []).append(r)
    rh = {}
    for r in runs:
        if r["config"].startswith("rho_dp"):
            dpb = int(r["config"][6:]); rh.setdefault(dpb, []).append(r)
    B = {m: {"ops_mean": float(np.mean([r["group_ops"] for r in rs])), "table_entries": float(np.mean([r["table_entries"] for r in rs])), "table_mb": float(np.mean([r.get("table_bytes", 0) for r in rs])) / 2**20,
             "rss_mb": float(np.mean([r["peak_rss_bytes"] for r in rs])) / 2**20, "wall": float(np.mean([r["wall_s"] for r in rs])), "n": float(np.mean([r["n"] for r in rs])), "pred": float(np.mean([m + r["n"] / (2 * m) for r in rs]))} for m, rs in bs.items()}
    Rr = {dpb: {"ops_mean": float(np.mean([r["group_ops"] for r in rs])), "ops_ci": boot_mean([r["group_ops"] for r in rs])[1:], "dp_count": float(np.mean([r["dp_count"] for r in rs])), "dp_mb": float(np.mean([r.get("dp_bytes", 0) for r in rs])) / 2**20,
                "rss_mb": float(np.mean([r["peak_rss_bytes"] for r in rs])) / 2**20, "wall": float(np.mean([r["wall_s"] for r in rs])), "n": float(np.mean([r["n"] for r in rs])), "iter": float(np.mean([r["iterations"] for r in rs]))} for dpb, rs in rh.items()}
    ANALYSIS["tmto"] = {"host": d["host"], "bsgs": B, "rho": Rr}
    fig, axes = plt.subplots(1, 2, figsize=(7.2, 3.2))
    ms = sorted(B); axes[0].plot([B[m]["table_entries"] for m in ms], [B[m]["ops_mean"] for m in ms], "o-", ms=3, label="BSGS measured (36 bits)")
    axes[0].plot([B[m]["table_entries"] for m in ms], [B[m]["pred"] for m in ms], "k--", lw=0.8, label=r"$m + n/(2m)$ (table keyed by $x$: $\pm j$ resolved)")
    axes[0].set_xscale("log", base=2); axes[0].set_yscale("log"); axes[0].set_xlabel("table entries (memory)"); axes[0].set_ylabel("group operations"); axes[0].legend(fontsize=7)
    ds = sorted(Rr); axes[1].errorbar([max(Rr[k]["dp_count"], 1) for k in ds], [Rr[k]["ops_mean"] for k in ds], yerr=[[Rr[k]["ops_mean"] - Rr[k]["ops_ci"][0] for k in ds], [Rr[k]["ops_ci"][1] - Rr[k]["ops_mean"] for k in ds]], fmt="s-", ms=3, capsize=2, label="rho, DP (40 bits)")
    for k in ds:
        axes[1].annotate(f"$2^{{-{k}}}$", (max(Rr[k]["dp_count"], 1), Rr[k]["ops_mean"]), xytext=(3, 3), textcoords="offset points", fontsize=6)
    axes[1].axhline(math.sqrt(math.pi * Rr[ds[0]]["n"] / 2), color="k", lw=0.8, ls="--", label=r"$\sqrt{\pi n/2}$")
    axes[1].set_xscale("log", base=2); axes[1].set_xlabel("stored distinguished points (memory)"); axes[1].set_ylabel("group operations"); axes[1].legend(fontsize=7)
    fig.tight_layout(); fig.savefig(F / "ecdlp_tmto.pdf"); fig.savefig(PF / "ecdlp_tmto.pdf"); plt.close(fig)
    rows = [r"\begin{tabular}{@{}lrrrrr@{}}", r"\toprule Algorithm / parameter & Memory (entries) & Memory (MB) & Group ops & Predicted & Wall (s) \\ \midrule"]
    for m in ms:
        v = B[m]; rows.append(f"BSGS, $m=2^{{{int(math.log2(m))}}}$ (36 bits) & {v['table_entries']:,.0f} & {v['table_mb']:.1f} & {v['ops_mean']:,.0f} & {v['pred']:,.0f} & {v['wall']:.3f} \\\\".replace(",", r"\,"))
    rows.append(r"\midrule")
    for k in ds:
        v = Rr[k]; rows.append(f"rho, DP density $2^{{-{k}}}$ (40 bits) & {v['dp_count']:,.0f} & {v['dp_mb']:.2f} & {v['ops_mean']:,.0f} & {math.sqrt(math.pi * v['n'] / 2):,.0f} & {v['wall']:.3f} \\\\".replace(",", r"\,"))
    rows += [r"\bottomrule", r"\end{tabular}"]
    (T / "ecdlp_tmto.tex").write_text("\\resizebox{\\textwidth}{!}{%\n" + "\n".join(rows) + "}\n")
    MACROS["tmtoBsgsMinM"] = f"2^{{{int(math.log2(min(ms)))}}}"; MACROS["tmtoBsgsMaxM"] = f"2^{{{int(math.log2(max(ms)))}}}"
    dev = max(abs(B[m]["ops_mean"] / B[m]["pred"] - 1) for m in ms); MACROS["tmtoBsgsMaxDev"] = f"{100*dev:.1f}"
    MACROS["tmtoRhoDpMin"] = str(min(ds)); MACROS["tmtoRhoDpMax"] = str(max(ds))
    MACROS["tmtoRhoOpsRange"] = f"{min(Rr[k]['ops_mean'] for k in ds)/1e6:.2f}--{max(Rr[k]['ops_mean'] for k in ds)/1e6:.2f}"


def main():
    analyse_scaling(); analyse_j0(); analyse_sweep(); analyse_parallel(); analyse_tmto()
    (R / "analysis.json").write_text(json.dumps(ANALYSIS, indent=1, default=str))
    (T / "numbers_ecdlp.tex").write_text("% generated by experiments/synthetic_ecdlp/analyze_ecdlp.py\n" + "\n".join(f"\\newcommand{{\\{k}}}{{{v}}}" for k, v in MACROS.items()) + "\n")
    print("studies analysed:", list(ANALYSIS), "| macros:", len(MACROS))


if __name__ == "__main__":
    main()