sciencehub-team/constellarationPublic

# /// script
# requires-python = ">=3.11"
# dependencies = ["marimo>=0.25.1", "numpy", "pandas", "pyarrow", "matplotlib"]
# ///

import marimo

__generated_with = "0.25.1"
app = marimo.App(width="medium", app_title="The ConStellaration leaderboard, scored again")

with app.setup:
    import numpy as np
    import pandas as pd

    ENTRIES = "data/entries.parquet"
    RESCORED = "data/rescored.parquet"
    TOLERANCE = 0.01
    FLUX_COMPRESSION_BOUND = 0.9
    TITLES = {"geometrical": "Geometrical", "simple_to_build": "Simple-to-build QI", "mhd_stable": "MHD-stable QI"}


@app.function
def hypervolume(gradient_scale_length, aspect_ratio, reference=(1.0, 20.0)) -> float:
    """Area dominated by the points in the plane of (−L∇B, A), up to the reference point.

    The MHD-stable score of constellaration 0.2.6 for feasible points; equal
    to pymoo's Hypervolume to 1e-12 on all 58 leaderboard submissions.
    """
    x = -np.asarray(gradient_scale_length, dtype=float)
    y = np.asarray(aspect_ratio, dtype=float)
    keep = (x < reference[0]) & (y < reference[1])
    x, y = x[keep], y[keep]
    order = np.argsort(x, kind="stable")
    x, y = x[order], y[order]
    edges = np.append(x, reference[0])
    lowest = np.minimum.accumulate(y) if len(y) else y
    return float(np.sum(np.diff(edges) * (reference[1] - lowest)))


@app.function
def load_entries(path: str = ENTRIES) -> pd.DataFrame:
    """The leaderboard as recorded on 25 October 2026: one row per single-objective
    submission and one per point of each MHD-stable submission."""
    return pd.read_parquet(path)


@app.function
def mhd_stable_scores(entries: pd.DataFrame, tolerance: float = TOLERANCE, top: int = 10) -> pd.DataFrame:
    """The `top` MHD-stable submissions by hypervolume when only points within `tolerance` count.

    All submissions are ranked; a resubmission with the same user and
    recorded score counts once. "recorded rank" is the order of the
    recorded scores.
    """
    points = entries[entries["problem"] == "mhd_stable"]
    rows = []
    for (user, submitted), group in points.groupby(["user", "submitted"]):
        kept = group[group["feasibility"] <= tolerance]
        rows.append({
            "user": user,
            "submitted": submitted[:10],
            "points": len(group),
            "recorded score": group["entry_score"].iloc[0],
            "score at this tolerance": hypervolume(kept["gradient_scale_length"], kept["aspect_ratio"]),
            "points in the last tenth of the band": float(np.mean((group["feasibility"] > 0.9 * TOLERANCE) & (group["feasibility"] <= TOLERANCE))),
        })
    table = pd.DataFrame(rows).sort_values("recorded score", ascending=False).drop_duplicates(["user", "recorded score"])
    table.insert(0, "recorded rank", range(1, len(table) + 1))
    table.insert(0, "rank at this tolerance", table["score at this tolerance"].rank(ascending=False, method="min").astype(int))
    return table.sort_values(["rank at this tolerance", "recorded rank"]).head(top).reset_index(drop=True)


@app.cell
def _():
    import marimo as mo
    import matplotlib.pyplot as plt

    return mo, plt


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    # The ConStellaration leaderboard, scored again

    The ConStellaration leaderboard on Hugging Face ranks stellarator boundaries
    with Proxima Fusion's evaluator, constellaration 0.2.6
    ([Cadena et al. 2025](https://arxiv.org/abs/2506.19583)). By 25 October
    2026 it held 174 evaluated submissions. The best entries had passed the
    paper's baselines on all three problems: a maximum elongation of 1.18
    against 1.27 on the geometrical problem, a normalised gradient scale
    length of 13.4 against 8.6 on the simple-to-build problem, and a
    hypervolume of 249 against 130 on the MHD-stable problem.

    We check whether the recorded scores reproduce on another machine, how
    much of each leading score depends on the 1 % tolerance the evaluator
    grants every constraint, and how the leading MHD-stable fronts change
    when the flux-compression proxy is averaged the way Proxima corrected it
    in version 0.3.0. The forward-model runs were made with `scripts/rescore.py`
    of this project; this notebook only reads their results.
    """)
    return


@app.cell
def _():
    entries = load_entries()
    rescored = pd.read_parquet(RESCORED)
    return entries, rescored


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    ## The recorded scores reproduce

    We evaluated the paper's baselines (the first submissions of the paper's
    first author), the three best single-objective entries per problem, and a
    random sample of points from the best MHD-stable front, in a fresh
    environment built like the leaderboard's (Python 3.10, constellaration
    0.2.6, booz-xform 0.0.9 compiled from source) on an x86-64 server. The
    reference boundary and the paper's baselines also ran in ScienceHub's
    kernel image (Python 3.13, booz-xform 0.1.0), which the evaluate
    notebook of this project uses.
    """)
    return


@app.cell
def _(mo, rescored):
    _checked = rescored[rescored["recorded_feasibility"].notna() & rescored["flux_compression_pest"].isna() & (rescored["run"] == "leaderboard environment")]
    _single = _checked[_checked["problem"] != "mhd_stable"]
    _multi = _checked[_checked["problem"] == "mhd_stable"]
    _table = pd.DataFrame({
        "entry": _single["user"] + ", " + _single["submitted"].str[:10],
        "problem": _single["problem"].map(TITLES),
        "recorded score": _single["recorded_entry_score"],
        "our score": _single["score"],
        "recorded feasibility": _single["recorded_feasibility"],
        "our feasibility": _single["feasibility"],
    })
    _gap_single = float((_single["score"] - _single["recorded_entry_score"]).abs().max())
    _gap_points = float((_multi["feasibility"] - _multi["recorded_feasibility"]).abs().max())
    _gap_objective = float((_multi["minimum_normalized_magnetic_gradient_scale_length"] - _multi["recorded_gradient_scale_length"]).abs().max())
    _metrics = list(rescored.loc[:, "feasibility":"flux_compression_in_regions_of_bad_curvature"].columns.drop(["flux_compression_pest", "flux_compression_vmec"]))
    _pairs = rescored[rescored["run"] == "kernel image"].merge(
        rescored[rescored["run"] == "leaderboard environment"], on=["user", "submitted", "problem", "index"], suffixes=("", "_leaderboard")
    )
    _kernel = _pairs[_metrics].to_numpy(dtype=float)
    _leaderboard = _pairs[[f"{_name}_leaderboard" for _name in _metrics]].to_numpy(dtype=float)
    assert np.array_equal(np.isnan(_kernel), np.isnan(_leaderboard))
    _values = int(np.sum(~np.isnan(_kernel)))
    _gap_kernel = float(np.nanmax(np.abs(_kernel - _leaderboard)))
    mo.vstack([
        mo.ui.table(_table, selection=None),
        mo.md(
            f"Largest difference in a single-objective score: **{_gap_single:.1e}**. "
            f"Over {len(_multi)} MHD-stable points, the largest difference in feasibility is **{_gap_points:.1e}** "
            f"and in L∇B **{_gap_objective:.1e}**. "
            f"For the {len(_pairs)} boundaries evaluated in both environments, the largest difference over "
            f"{_values} metric values, feasibility and score included, is **{_gap_kernel:g}**."
        ),
    ])
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    The evaluator is deterministic: VMEC++ runs single-threaded, and the
    same boundary gives the same numbers on another machine to the precision
    shown. The kernel image, with Python 3.13, booz-xform 0.1.0, simsopt
    1.9.0, jax 0.5.3 and scipy 1.17.1 instead of the leaderboard's versions,
    gives the same metrics bit for bit. A score on
    the leaderboard can therefore be checked by anyone with a minute of CPU
    time per boundary and about 3 GiB of memory.

    ## How much of each score lies in the tolerance

    A boundary counts as feasible when every normalised constraint violation
    is at most 0.01. Optimisers use this margin: the best entry on each
    single-objective problem sits at feasibility 0.0090 (geometrical) and
    0.0097 (simple-to-build), and on the
    ten leading MHD-stable fronts between 17 % and 62 % of the points lie in
    the last tenth of the band, above 0.009. The table ranks the 55 distinct
    MHD-stable submissions (58 before removing resubmissions) by the hypervolume of their points within the
    tolerance you choose and shows the ten best.
    """)
    return


@app.cell
def _(mo):
    tolerance = mo.ui.slider(0.0, 0.01, step=0.0005, value=0.001, label="Tolerance", show_value=True, full_width=True)
    tolerance
    return (tolerance,)


@app.cell
def _(entries, mo, tolerance):
    _table = mhd_stable_scores(entries, tolerance.value)
    mo.ui.table(_table.round({"recorded score": 2, "score at this tolerance": 2, "points in the last tenth of the band": 2}), selection=None)
    return


@app.cell
def _(entries, mo, plt):
    _grid = np.linspace(0, TOLERANCE, 41)
    _leaders = pd.concat([mhd_stable_scores(entries, TOLERANCE, top=4), mhd_stable_scores(entries, 0.001, top=1)])
    _fig, _ax = plt.subplots(figsize=(7, 3.6))
    _points = entries[entries["problem"] == "mhd_stable"]
    for (_, _row), _color in zip(_leaders.iterrows(), ("#1f2328", "#0969da", "#9a6700", "#57606a", "#cf222e")):
        _group = _points[(_points["user"] == _row["user"]) & (_points["submitted"].str[:10] == _row["submitted"]) & (_points["entry_score"] == _row["recorded score"])]
        _curve = [hypervolume(*_group.loc[_group["feasibility"] <= _t, ["gradient_scale_length", "aspect_ratio"]].to_numpy().T) for _t in _grid]
        _ax.plot(_grid, _curve, color=_color, lw=1.6, label=f"{_row['user']} ({_row['recorded score']:.1f})")
    _ax.set_xlabel("tolerance on the normalised constraint violation")
    _ax.set_ylabel("hypervolume")
    _ax.legend(frameon=False, fontsize=8, loc="lower right")
    for _side in ("top", "right"):
        _ax.spines[_side].set_visible(False)
    _fig.tight_layout()
    _fig
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    At a tolerance of 0.001 the four leading hypervolumes fall from 243–249
    to 167–186: between a quarter and a third of each leading score is
    earned by points that break a constraint by more than 0.1 %. Most of
    their points sit in the upper part of the band: 80 % of Eppy's 688
    points and 86 % of Wenyang Li's 783 lie above 0.005. Other fronts lose
    less. The front canyon-ant submitted on 9 September, 16th on the
    leaderboard with 43 points, keeps 94 % of its hypervolume and scores
    192.5, the highest at this tolerance. The front phanerozoic submitted on
    15 September, 15th, keeps 93 % with only 15 of its 668 points and scores
    192.2. Wenyang Li follows with 185.8 and the recorded leader, Eppy, with
    184.1; at a tolerance of 0 Wenyang Li leads with 185.7 ahead of Eppy
    with 183.9. A front scores well at a tight tolerance when its few points
    with margin are spread along it. The leaderboard rewards only the
    official tolerance, so the order at 0.001 shows how much of each score
    rests on the band, not which method would win a benchmark with a
    tighter one. (The third and fourth recorded submissions tie at tight
    tolerances because 507 of the fourth's 512 points also appear in the
    third.) The tolerance is part of the benchmark's definition, so none of
    this changes the leaderboard. It does mean that a comparison which cares
    about designs with margin should report the score at a tighter
    tolerance beside the official one.
    """)
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    ## The flux-compression correction

    The MHD-stable problem bounds a geometric proxy for turbulent transport:
    the flux-surface average, at ρ = 0.7, of $|\nabla r|^2$ where the
    curvature drift is bad and zero elsewhere. Version 0.2.6 averages with the VMEC Jacobian; in
    version 0.3.0 Proxima changed the weight to the field-aligned (PEST)
    Jacobian, calling the old one a bug
    ([constellaration#108](https://github.com/proximafusion/constellaration/issues/108),
    [#109](https://github.com/proximafusion/constellaration/pull/109)). The
    leaderboard stays on 0.2.6.

    To see what the correction alone changes, we solved each boundary's
    equilibrium once with 0.2.6 and computed the proxy with both weights from
    that same equilibrium. A point that was feasible under 0.2.6 stays
    feasible exactly when the corrected proxy is at most 0.9 × 1.01; no other
    constraint is affected. The two leading fronts were evaluated in full
    (688 and 783 points). The third and fourth recorded fronts, 1,364 and 512
    points, were evaluated at 120 points each, drawn at random.
    """)
    return


@app.cell
def _(rescored):
    flux = rescored[rescored["run"] == "flux compression only"].copy()
    _solved = flux["error"].isna()
    # The VMEC-weighted value recomputed here must be the one 0.2.6 reports.
    assert np.allclose(
        flux.loc[_solved, "flux_compression_vmec"].astype(float),
        flux.loc[_solved, "flux_compression_in_regions_of_bad_curvature"].astype(float),
        rtol=1e-12,
    )
    # A failed equilibrium counts as not feasible and stays in the point counts.
    flux["feasible_after_correction"] = (
        _solved
        & (flux["recorded_feasibility"] <= TOLERANCE)
        & (flux["flux_compression_pest"].astype(float) <= FLUX_COMPRESSION_BOUND * (1 + TOLERANCE))
    )
    return (flux,)


@app.cell
def _(flux, mo, plt):
    mo.stop(flux.empty, mo.md("The flux-compression runs are not in `data/rescored.parquet` yet."))
    _fig, _ax = plt.subplots(figsize=(5.5, 4.5))
    _ax.scatter(flux["flux_compression_vmec"], flux["flux_compression_pest"], s=6, color="#1f2328", alpha=0.6, lw=0)
    _limit = FLUX_COMPRESSION_BOUND * (1 + TOLERANCE)
    _ax.axhline(_limit, color="#cf222e", lw=0.8, ls="--")
    _ax.axvline(_limit, color="#cf222e", lw=0.8, ls="--")
    _lo, _hi = float(flux[["flux_compression_vmec", "flux_compression_pest"]].min().min()), float(flux[["flux_compression_vmec", "flux_compression_pest"]].max().max())
    _ax.plot([_lo, _hi], [_lo, _hi], color="#57606a", lw=0.8)
    _ax.set_xlabel("flux compression, VMEC weight (0.2.6)")
    _ax.set_ylabel("flux compression, PEST weight (0.3.0)")
    for _side in ("top", "right"):
        _ax.spines[_side].set_visible(False)
    _fig.tight_layout()
    _fig
    return


@app.cell
def _(entries, flux, mo):
    mo.stop(flux.empty)
    _points = entries[entries["problem"] == "mhd_stable"]
    _rows = []
    for (_user, _submitted), _group in flux.groupby(["user", "submitted"]):
        _all = _points[(_points["user"] == _user) & (_points["submitted"] == _submitted)]
        _complete = len(_group) == len(_all)
        _n = len(_group)
        _share = float(_group["feasible_after_correction"].mean())
        # Wilson score interval, 95 %.
        _z = 1.96
        _centre = (_share + _z**2 / (2 * _n)) / (1 + _z**2 / _n)
        _half_width = _z * np.sqrt(_share * (1 - _share) / _n + _z**2 / (4 * _n**2)) / (1 + _z**2 / _n)
        _row = {
            "user": _user,
            "submitted": _submitted[:10],
            "points evaluated": f"{_n} of {len(_all)}",
            "failed runs": int(_group["error"].notna().sum()),
            "still feasible": f"{_share:.0%}" + ("" if _complete else f" (95 % interval {_centre - _half_width:.0%} to {_centre + _half_width:.0%})"),
            "recorded score": round(float(_all["entry_score"].iloc[0]), 2),
        }
        _ok = _group[_group["feasible_after_correction"]]
        _recorded = _all.set_index("index").loc[_ok["index"]]
        _score = hypervolume(_recorded["gradient_scale_length"], _recorded["aspect_ratio"])
        # Adding points never lowers a hypervolume, so a sample's score is a lower bound.
        _row["score after correction"] = f"{_score:.2f}" if _complete else f"at least {_score:.2f}"
        _rows.append(_row)
    mo.ui.table(pd.DataFrame(_rows).sort_values("recorded score", ascending=False), selection=None)
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    The corrected weight raises the proxy for most points. The median ratio
    of the PEST-weighted to the VMEC-weighted value is 1.09 on Eppy's front
    and 1.10 on Wenyang Li's, against 1.04 on phanerozoic's and 1.03 on
    mmdeuss's; over all 1,711 points it lies between 0.96 and 1.20. Every point
    of the two leading fronts was feasible under 0.2.6 and every equilibrium
    solved, but only 181 of Eppy's 688 points (26 %) and 140 of Wenyang Li's
    783 (18 %) stay within the bound. Their scores fall from 249.14 to
    208.83 and from 248.23 to 208.05, so the two keep their order, 0.78
    apart instead of 0.91.

    The sampled fronts lose fewer points: 47 of the 120 sampled from
    phanerozoic's front stay within the bound (39 %, 95 % interval
    31–48 %), and 58 of the 120 from mmdeuss's (48 %, 40–57 %). Their
    corrected scores are not known, but the surviving sample points alone
    score 220.47 and 223.84, and adding the rest of a front can only raise
    that. Both lower bounds are above the corrected
    scores of the first two fronts. Under the corrected proxy, the third and
    fourth recorded fronts would therefore both rank above the first and
    second; the samples do not settle the order between those two, and
    fronts below the fourth were not evaluated, so we cannot say which
    submission would lead.

    This is a statement about the benchmark, not about the designs. The
    proxy is a geometric stand-in for turbulent transport, not a transport
    calculation, and each submission was optimised against the 0.2.6
    version of it, which the leaderboard still uses. A front optimised
    against the corrected proxy could win back much of what it loses here.
    What the result shows is that the order at the top of the MHD-stable
    leaderboard depends on how one constraint is averaged, and that the
    correction raises the proxy more on the two leading fronts than on the
    next two.
    """)
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    ## Limits and references

    **What was run.** Scores in the first section come from a fresh
    installation of constellaration 0.2.6 on Python 3.10 with booz-xform
    0.0.9 built from source, in the Docker image `python:3.10-slim`, on a
    Hetzner CPX42 server; the flux-compression runs use the same environment
    with the Boozer transform and the QI residual switched off, which the
    proxy does not depend on. The tolerance analysis needs no new runs: it
    uses the per-point feasibility the leaderboard recorded.

    **What was not run.** The third and fourth fronts are sampled, so their
    share of surviving points carries a 95 % Wilson interval and their
    corrected score is a lower bound. A boundary whose equilibrium failed
    would count as not feasible; none of the 1,711 runs failed. The
    corrected score keeps every other metric from 0.2.6 and so isolates the
    change of weight; it is not an evaluation with constellaration 0.3.0. A
    front that loses points under the correction could recover some
    hypervolume by re-optimisation, which we did not attempt.

    **Data.** `data/entries.parquet` is derived from the results dataset of
    the leaderboard Space,
    [proxima-fusion/constellaration-bench-results](https://huggingface.co/datasets/proxima-fusion/constellaration-bench-results)
    (178 rows, 174 of them evaluated submissions from users other than
    `test`, last modified 25 October 2026, no licence stated). It
    keeps the recorded objectives, feasibilities and scores with the
    submitters' public names and dates, and none of the submitted boundaries.
    `data/rescored.parquet` holds our own evaluations; `build_data.py` writes
    both files.

    Cadena, S. A. et al. (2025). ConStellaration: A dataset of QI-like
    stellarator plasma boundaries and optimization benchmarks. *NeurIPS 2025
    Datasets and Benchmarks Track*. [arXiv:2506.19583](https://arxiv.org/abs/2506.19583).
    """)
    return


if __name__ == "__main__":
    app.run()