Historical Continuity as
Justification for Dark Matter

def initialize(root, output_dir=None, data_dir=None):
    root = Path(root).resolve()
    expected = json.loads((root / "checks/expected-results.json").read_text())
    out = Path(output_dir).resolve() if output_dir else root / "reproduced"
    data = Path(data_dir).resolve() if data_dir else root / "inputs"
    protected = [
        root / name
        for name in [
            "analysis",
            "inputs",
            "checks",
            "scripts",
            "config",
            "resources",
            "figures",
            "manuscript-figures",
            ".git",
        ]
    ]
    if (
        out == root
        or any(out == p or p in out.parents for p in protected)
        or out == data
        or data in out.parents
    ):
        raise ValueError(
            "Choose a separate reconstruction directory; supplied evidence is protected"
        )
    from materialize_inputs import materialize

    materialize(data)
    checked = {}
    for name, entry in expected["inputs"].items():
        path = data / (name + ".parquet")
        actual = sha256(path)
        if actual != entry["sha256"]:
            raise ValueError("Input checksum differs: " + path.name)
        checked[name] = actual
    (out / "analysis").mkdir(parents=True, exist_ok=True)
    (out / "figures").mkdir(exist_ok=True)
    return {
        "root": root,
        "data": data,
        "output_root": out,
        "analysis_dir": out / "analysis",
        "figure_dir": out / "figures",
        "expected": expected,
        "checked_inputs": checked,
    }

"""Restore the largest input table byte-for-byte from its transport parts."""

from pathlib import Path
import hashlib
import json
import os
import shutil

ROOT = Path(__file__).resolve().parents[1]


def sha256(path):
    h = hashlib.sha256()
    with path.open("rb") as f:
        for block in iter(lambda: f.read(1024 * 1024), b""):
            h.update(block)
    return h.hexdigest()


def materialize(data_dir=None):
    data = Path(data_dir) if data_dir is not None else ROOT / "inputs"
    spec = json.loads((ROOT / "INPUT-STORAGE.json").read_text())
    target = data / spec["table"]
    if target.is_file():
        if sha256(target) != spec["sha256"]:
            raise ValueError("Existing input table has a different checksum.")
        return target
    missing = [p["file"] for p in spec["parts"] if not (data / p["file"]).is_file()]
    if missing:
        raise FileNotFoundError(
            "Full input package required; missing transport parts: "
            + ", ".join(missing)
        )
    temporary = target.with_suffix(".assembling")
    try:
        with temporary.open("wb") as output:
            for part in spec["parts"]:
                path = data / part["file"]
                if (
                    path.stat().st_size != part["bytes"]
                    or sha256(path) != part["sha256"]
                ):
                    raise ValueError(
                        "Transport-part checksum differs: " + part["file"]
                    )
                with path.open("rb") as source:
                    shutil.copyfileobj(source, output)
        if (
            temporary.stat().st_size != spec["bytes"]
            or sha256(temporary) != spec["sha256"]
        ):
            raise ValueError("Assembled input checksum differs.")
        os.replace(temporary, target)
    finally:
        if temporary.exists():
            temporary.unlink()
    return target


if __name__ == "__main__":
    print("Verified input:", materialize().name)

from pathlib import Path
import sys, os, json, subprocess

start = Path.cwd().resolve()
root = next(
    (p for p in [start, *start.parents] if (p / "config/workspace.json").is_file()),
    None,
)
if root is None:
    raise FileNotFoundError("Start this notebook inside its computational package.")
(root / "verification").mkdir(exist_ok=True)
cache = root / ".runtime/matplotlib"
cache.mkdir(parents=True, exist_ok=True)
os.environ["MPLCONFIGDIR"] = str(cache)
sys.path.insert(0, str(root / "scripts"))
from modules import analysis, plotting

state = analysis.initialize(root)
print("Verified five local input hashes:")
for name, checksum in state["checked_inputs"].items():
    print(name, checksum)

def build_model_sample(ctx):
    OUT = ctx["analysis_dir"]
    DATA = ctx["data"]

    def read(name, **kwargs):
        return pd.read_parquet(DATA / (name + ".parquet"), **kwargs)

    papers = read("papers", columns=["bibcode", "year", "doctype", "title"])
    classes = (
        read("paper_arxiv_classes")
        .sort_values(["bibcode", "class_pos"])
        .drop_duplicates("bibcode")
    )
    meta = papers.merge(
        classes[["bibcode", "arxiv_class", "arxiv_category"]],
        on="bibcode",
        how="left",
        validate="one_to_one",
    )
    cites = read(
        "paper_citations", columns=["bibcode", "citing_bibcode"]
    ).drop_duplicates()
    citers = set(cites.loc[cites.bibcode.eq(BH), "citing_bibcode"])
    refyears = pd.to_numeric(
        cites.loc[cites.citing_bibcode.eq(BH), "bibcode"].str[:4], errors="coerce"
    )
    bh_refs = cites.loc[cites.citing_bibcode.eq(BH)].copy()
    bh_refs["ref_year"] = pd.to_numeric(bh_refs.bibcode.str[:4], errors="coerce")
    early = set(bh_refs.loc[bh_refs.ref_year.between(1800, 1959), "bibcode"])
    pd.DataFrame({"bibcode": sorted(early)}).merge(
        papers[["bibcode", "title"]], how="left"
    ).to_csv(OUT / "bh-linked-early-references.csv", index=False)
    cohort = meta.loc[
        meta.year.between(2018, 2025)
        & meta.doctype.isin(["article", "eprint"])
        & meta.arxiv_category.notna()
    ].copy()
    cohort["cites_bh"] = cohort.bibcode.isin(citers).astype(int)
    refs = cites[cites.citing_bibcode.isin(cohort.bibcode)].copy()
    refs["ref_year"] = pd.to_numeric(refs.bibcode.str[:4], errors="coerce")
    refs = refs.loc[refs.ref_year.between(1800, 2025)]
    refs["in_early"] = refs.bibcode.isin(early)
    refs["other_early"] = refs.ref_year.lt(1960) & ~refs.in_early
    signals = (
        refs.groupby("citing_bibcode")
        .agg(
            n_linked_refs=("bibcode", "size"),
            bh_early=("in_early", "max"),
            other_early=("other_early", "max"),
        )
        .reset_index()
    )
    cohort = cohort.merge(
        signals,
        left_on="bibcode",
        right_on="citing_bibcode",
        how="inner",
        validate="one_to_one",
    )
    keep = [
        "astrophysics",
        "high-energy physics",
        "general relativity and quantum cosmology",
        "physics",
        "nuclear physics",
        "quantum physics",
    ]
    cohort = cohort[cohort.arxiv_category.isin(keep)].copy()
    cohort["year_bin"] = pd.cut(
        cohort.year,
        [2018, 2020, 2022, 2024, 2026],
        right=False,
        labels=["2018–2019", "2020–2021", "2022–2023", "2024–2025"],
    ).astype(str)
    cohort["stratum"] = cohort.year_bin + " / " + cohort.arxiv_category
    cohort["log_refs"] = np.log1p(cohort.n_linked_refs)
    cohort[
        [
            "bibcode",
            "year",
            "arxiv_category",
            "cites_bh",
            "n_linked_refs",
            "bh_early",
            "other_early",
            "stratum",
        ]
    ].to_csv(OUT / "co-citation-analysis.csv", index=False)
    ctx["papers"] = papers
    ctx["meta"] = meta
    ctx["cites"] = cites
    ctx["citers"] = citers
    ctx["early"] = early
    ctx["cohort"] = cohort
    ctx["keep"] = keep
    return ctx

analysis.build_model_sample(state)
print("Eligible model records:", len(state["cohort"]))
print("Model-sample review citers:", int(state["cohort"].cites_bh.sum()))
print("Recoverable early references:", len(state["early"]))
assert len(state["cohort"]) == 55140
state["cohort"].groupby("cites_bh")[["bh_early", "other_early"]].agg(
    ["sum", "mean"]
)

def fit_models_and_strata(ctx):
    OUT = ctx["analysis_dir"]
    DATA = ctx["data"]

    def read(name, **kwargs):
        return pd.read_parquet(DATA / (name + ".parquet"), **kwargs)

    cohort = ctx["cohort"]
    results = {}
    models = {}
    for outcome in ["bh_early", "other_early"]:
        cohort[outcome] = cohort[outcome].astype(int)
        with warnings.catch_warnings(record=True) as caught:
            model = smf.glm(
                outcome + " ~ cites_bh + log_refs + C(stratum)",
                data=cohort,
                family=sm.families.Binomial(),
            ).fit(maxiter=100)
        models[outcome] = model
        beta = float(model.params["cites_bh"])
        ci = model.conf_int().loc["cites_bh"]
        a = cohort.assign(cites_bh=1)
        b = cohort.assign(cites_bh=0)
        diff = model.predict(a) - model.predict(b)
        results[outcome] = {
            "odds_ratio": float(np.exp(beta)),
            "wald_ci95": [float(np.exp(ci.iloc[0])), float(np.exp(ci.iloc[1]))],
            "mean_standardized_difference": float(diff.mean()),
            "median_standardized_difference": float(diff.median()),
            "converged": bool(model.converged),
            "warnings": sorted(set(str(x.message) for x in caught)),
        }
        model.summary2().tables[1].to_csv(OUT / (outcome + "-coefficients.csv"))
    strata = (
        cohort.groupby(["year_bin", "arxiv_category", "stratum"])
        .agg(
            mean_log_refs=("log_refs", "mean"),
            n_bh=("cites_bh", "sum"),
            n_total=("cites_bh", "size"),
        )
        .reset_index()
    )
    strata["n_nonbh"] = strata.n_total - strata.n_bh
    strata = strata[(strata.n_bh >= 10) & (strata.n_nonbh >= 100)].copy()
    strata["difference"] = models["bh_early"].predict(
        strata.assign(cites_bh=1, log_refs=strata.mean_log_refs)
    ) - models["bh_early"].predict(
        strata.assign(cites_bh=0, log_refs=strata.mean_log_refs)
    )
    strata.to_csv(OUT / "stratum-differences.csv", index=False)
    ctx["cohort"] = cohort
    ctx["results"] = results
    ctx["models"] = models
    ctx["strata"] = strata
    return ctx

analysis.fit_models_and_strata(state)
print(json.dumps(state["results"], indent=2))
state["strata"][["year_bin", "arxiv_category", "n_bh", "n_nonbh", "difference"]]

def verify_citation_rates(ctx):
    OUT = ctx["analysis_dir"]
    DATA = ctx["data"]

    def read(name, **kwargs):
        return pd.read_parquet(DATA / (name + ".parquet"), **kwargs)

    papers = ctx["papers"]
    rates = read("avg_citations_result")
    assert rates.bibcode.is_unique
    # Independently verify the stored rates against annual citation metrics, using
    # the original notebook's inclusive publication-to-latest-metric convention.
    annual = pd.read_parquet(
        DATA / "paper_metrics_long.parquet", filters=[("metric", "==", "citations")]
    )
    annual = annual.merge(
        papers[["bibcode", "year"]],
        on="bibcode",
        how="left",
        validate="many_to_one",
    )
    annual = annual[annual.metric_year.ge(annual.year)]
    rate_check = annual.groupby("bibcode").agg(
        publication_year=("year", "first"),
        latest_metric_year=("metric_year", "max"),
        total_citations=("value", "sum"),
    )
    rate_check["years_exposed"] = (
        rate_check.latest_metric_year - rate_check.publication_year + 1
    )
    rate_check["recomputed_rate"] = (
        rate_check.total_citations / rate_check.years_exposed
    )
    assert set(rate_check.index) == set(rates.bibcode)
    rate_check = rate_check.join(
        rates.set_index("bibcode")[["avg_citations_per_year"]],
        validate="one_to_one",
    )
    assert np.allclose(
        rate_check.recomputed_rate, rate_check.avg_citations_per_year
    )
    rate_check.to_csv(OUT / "citation-rate-verification.csv")
    del annual
    expected = (
        rankdata(rates.avg_citations_per_year.to_numpy(), method="average")
        / len(rates)
        * 100
    )
    assert np.allclose(expected, rates.percentile)
    target = rates.loc[rates.bibcode.eq(BH)].iloc[0]
    pos = rates[rates.avg_citations_per_year.gt(0)].sort_values(
        "avg_citations_per_year"
    )
    ctx["rates"] = rates
    ctx["rate_check"] = rate_check
    ctx["target"] = target
    ctx["pos"] = pos
    return ctx

analysis.verify_citation_rates(state)
assert len(state["rates"]) == 158367
print("Verified rate records:", len(state["rate_check"]))
print("Positive-rate records:", len(state["pos"]))
print("Review annual rate:", float(state["target"].avg_citations_per_year))
print("Review percentile:", float(state["target"].percentile))

def write_counts_and_summary(ctx):
    OUT = ctx["analysis_dir"]
    DATA = ctx["data"]

    def read(name, **kwargs):
        return pd.read_parquet(DATA / (name + ".parquet"), **kwargs)

    papers = ctx["papers"]
    meta = ctx["meta"]
    citers = ctx["citers"]
    early = ctx["early"]
    cohort = ctx["cohort"]
    results = ctx["results"]
    rate_check = ctx["rate_check"]
    rates = ctx["rates"]
    target = ctx["target"]
    pos = ctx["pos"]
    circulation = (
        meta[meta.bibcode.isin(citers) & meta.arxiv_category.notna()]
        .groupby("arxiv_category")
        .size()
        .sort_values(ascending=False)
    )
    circulation.to_csv(OUT / "bh-citers-by-field.csv", header=["papers"])
    assert circulation.loc[
        [
            "high-energy physics",
            "astrophysics",
            "general relativity and quantum cosmology",
            "physics",
        ]
    ].tolist() == [233, 180, 76, 48]
    # Comparators are named historical treatments, not an exhaustive genre population.
    comparator_ids = [
        BH,
        "2017NatAs...1E..59D",
        "2010dmp..book.....S",
        "1999PASP..111..657V",
    ]
    comparators = meta[meta.bibcode.isin(comparator_ids)].merge(
        rates, on="bibcode", how="left"
    )
    comparators[
        ["bibcode", "title", "year", "avg_citations_per_year", "percentile"]
    ].to_csv(OUT / "historical-comparators.csv", index=False)
    inputs = {
        n: {
            "path": "inputs/" + n + ".parquet",
            "sha256": hashlib.sha256(
                (DATA / (n + ".parquet")).read_bytes()
            ).hexdigest(),
        }
        for n in [
            "papers",
            "paper_citations",
            "paper_arxiv_classes",
            "paper_metrics_long",
            "avg_citations_result",
        ]
    }
    summary = {
        "full_corpus": len(papers),
        "year_range": [int(papers.year.min()), int(papers.year.max())],
        "rate_records": len(rates),
        "rate_records_verified_from_annual_metrics": len(rate_check),
        "rate_max_absolute_error": float(
            (rate_check.recomputed_rate - rate_check.avg_citations_per_year)
            .abs()
            .max()
        ),
        "positive_rate_records": len(pos),
        "bh_rate": float(target.avg_citations_per_year),
        "bh_percentile": float(target.percentile),
        "bh_graph_citers": len(citers),
        "classified_bh_citers": int(circulation.sum()),
        "linked_early_reference_count": len(early),
        "analysis_n": len(cohort),
        "analysis_bh_citers": int(cohort.cites_bh.sum()),
        "analysis_outcome_counts": cohort.groupby("cites_bh")[
            ["bh_early", "other_early"]
        ]
        .agg(["sum", "mean"])
        .to_json(),
        "results": results,
        "inputs": inputs,
        "limits": [
            "Retained linked references are incomplete relative to the full bibliographies.",
            "Reference-list length means linked, date-eligible references, not the full reference list.",
            "Citation overlap is an association, not interpretive adherence or causal impact.",
            "No original ADS retrieval timestamp is embedded in these parquet inputs.",
        ],
    }
    (OUT / "bibliometric-summary.json").write_text(
        json.dumps(summary, indent=2) + "\n"
    )
    ctx["circulation"] = circulation
    ctx["comparators"] = comparators
    ctx["summary"] = summary
    return ctx

analysis.write_counts_and_summary(state)
assert int(state["circulation"].sum()) == 551
print(state["circulation"].to_string())
print("Graph review citers:", len(state["citers"]))
state["comparators"][["bibcode", "avg_citations_per_year", "percentile"]]

"""Refit both reported models and verify the numerical evidence for the revision."""

from pathlib import Path
import argparse
import json
import numpy as np
import pandas as pd
from scipy.stats import rankdata
import statsmodels.api as sm
import statsmodels.formula.api as smf

ROOT = Path(__file__).resolve().parents[1]
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument(
    "--reproduced",
    type=Path,
    help="Also compare tables rebuilt from retained inputs.",
)
parser.add_argument("--report", type=Path, help="Write a JSON verification report.")
args = parser.parse_args()
expected = json.loads((ROOT / "checks/expected-results.json").read_text())
a = ROOT / "analysis"
c = pd.read_csv(a / "co-citation-analysis.csv")
assert c.bibcode.is_unique and len(c) == expected["analysis_n"]
assert int(c.cites_bh.sum()) == expected["analysis_bh_citers"]
assert c.year.between(2018, 2025).all() and (c.n_linked_refs >= 1).all()
assert set(c.cites_bh.unique()) == {0, 1}
c["log_refs"] = np.log1p(c.n_linked_refs)
models = {}
for outcome in ["bh_early", "other_early"]:
    c[outcome] = c[outcome].astype(int)
    fit = smf.glm(
        outcome + " ~ cites_bh + log_refs + C(stratum)",
        data=c,
        family=sm.families.Binomial(),
    ).fit(maxiter=100)
    assert fit.converged
    saved = pd.read_csv(a / (outcome + "-coefficients.csv"), index_col=0)
    assert list(saved.index) == list(fit.params.index)
    np.testing.assert_allclose(fit.params, saved["Coef."], rtol=1e-9, atol=1e-10)
    np.testing.assert_allclose(fit.bse, saved["Std.Err."], rtol=1e-9, atol=1e-10)
    np.testing.assert_allclose(
        fit.conf_int(), saved[["[0.025", "0.975]"]], rtol=1e-9, atol=1e-10
    )
    estimate = float(np.exp(fit.params["cites_bh"]))
    interval = np.exp(fit.conf_int().loc["cites_bh"]).tolist()
    delta = float(
        (
            fit.predict(c.assign(cites_bh=1)) - fit.predict(c.assign(cites_bh=0))
        ).mean()
    )
    target = expected["results"][outcome]
    np.testing.assert_allclose(estimate, target["odds_ratio"], rtol=1e-9)
    np.testing.assert_allclose(interval, target["wald_ci95"], rtol=1e-9)
    np.testing.assert_allclose(
        delta, target["mean_standardized_difference"], rtol=1e-9
    )
    models[outcome] = {
        "odds_ratio": estimate,
        "wald_ci95": interval,
        "mean_standardized_difference": delta,
        "converged": True,
    }
    if outcome == "bh_early":
        strata = pd.read_csv(a / "stratum-differences.csv")
        counts = c.groupby("stratum").agg(
            mean_log_refs=("log_refs", "mean"),
            n_bh=("cites_bh", "sum"),
            n_total=("bibcode", "size"),
        )
        for row in strata.itertuples():
            values = counts.loc[row.stratum]
            np.testing.assert_allclose(row.mean_log_refs, values.mean_log_refs)
            assert row.n_bh == values.n_bh and row.n_total == values.n_total
            assert (
                row.n_nonbh == row.n_total - row.n_bh
                and row.n_bh >= 10
                and row.n_nonbh >= 100
            )
        differences = fit.predict(
            strata.assign(cites_bh=1, log_refs=strata.mean_log_refs)
        ) - fit.predict(strata.assign(cites_bh=0, log_refs=strata.mean_log_refs))
        np.testing.assert_allclose(differences, strata.difference, rtol=1e-9)

rates = pd.read_csv(a / "citation-rate-verification.csv")
assert rates.bibcode.is_unique and len(rates) == expected["rate_records"]
assert (
    rates.years_exposed == rates.latest_metric_year - rates.publication_year + 1
).all()
np.testing.assert_allclose(
    rates.total_citations / rates.years_exposed,
    rates.avg_citations_per_year,
    rtol=1e-9,
)
np.testing.assert_allclose(
    rates.recomputed_rate, rates.avg_citations_per_year, rtol=1e-9
)
rates["independent_percentile"] = (
    rankdata(rates.avg_citations_per_year, method="average") / len(rates) * 100
)
bh = rates.set_index("bibcode").loc["2018RvMP...90d5002B"]
np.testing.assert_allclose(bh.avg_citations_per_year, expected["bh_rate"])
np.testing.assert_allclose(bh.independent_percentile, expected["bh_percentile"])
comparators = pd.read_csv(a / "historical-comparators.csv")
for row in comparators.itertuples():
    matching = rates.set_index("bibcode").loc[row.bibcode]
    np.testing.assert_allclose(
        row.avg_citations_per_year, matching.avg_citations_per_year
    )
    np.testing.assert_allclose(row.percentile, matching.independent_percentile)
fields = (
    pd.read_csv(a / "bh-citers-by-field.csv").set_index("arxiv_category").papers
)
assert fields.sum() == expected["classified_bh_citers"] == 551
assert fields.loc[
    [
        "high-energy physics",
        "astrophysics",
        "general relativity and quantum cosmology",
        "physics",
    ]
].tolist() == [233, 180, 76, 48]
assert (
    len(pd.read_csv(a / "bh-linked-early-references.csv"))
    == expected["linked_early_reference_count"]
    == 7
)

compared = []
if args.reproduced:
    for p in sorted(a.glob("*.csv")):
        rebuilt = args.reproduced / "analysis" / p.name
        pd.testing.assert_frame_equal(
            pd.read_csv(p),
            pd.read_csv(rebuilt),
            check_exact=False,
            rtol=1e-9,
            atol=1e-10,
        )
        compared.append(p.name)
    summary = json.loads(
        (args.reproduced / "analysis/bibliometric-summary.json").read_text()
    )
    assert (
        summary == expected
    ), "Rebuilt numerical summary differs from the revision targets."
report = {
    "passed": True,
    "model_records": len(c),
    "bh_model_citers": int(c.cites_bh.sum()),
    "rate_records_verified": len(rates),
    "bh_rate": float(bh.avg_citations_per_year),
    "bh_percentile": float(bh.independent_percentile),
    "models": models,
    "rebuilt_tables_compared": compared,
}
if args.report:
    args.report.write_text(json.dumps(report, indent=2) + "\n")
print(json.dumps(report, indent=2))

check = subprocess.run(
    [
        sys.executable,
        str(root / "scripts/check_results.py"),
        "--reproduced",
        "reproduced",
        "--report",
        "verification/results.json",
    ],
    cwd=root,
    capture_output=True,
    text=True,
)
if check.returncode:
    raise RuntimeError(check.stdout + check.stderr)
report = json.loads((root / "verification/results.json").read_text())
assert report["passed"] and len(report["rebuilt_tables_compared"]) == 8
print(
    "Both models, rate identities, summary, and all eight CSV comparisons passed."
)
print("Compared tables:", *report["rebuilt_tables_compared"], sep="\n")

def export_all(root):
    root = Path(root).resolve()
    analysis.initialize(root)  # Verify input identities before any export.
    out = root / "reproduced/figures"
    out.mkdir(parents=True, exist_ok=True)
    configure_fonts(root / "resources/fonts")
    records = [export_impact(root, out)]
    result = subprocess.run(
        [
            "Rscript",
            "--vanilla",
            str(root / "scripts/02-export-waffle.R"),
            str(root),
        ],
        cwd=root,
        capture_output=True,
        text=True,
    )
    (root / "verification/waffle-export.log").write_text(
        result.stdout + result.stderr
    )
    if result.returncode:
        raise RuntimeError(
            "R field-display export failed; see verification/waffle-export.log"
        )
    validate_waffle_pdf(root)
    records.append(export_strata(root, out))
    records.append(
        {
            "figure": "fig-arxiv-citing-patterns",
            "pdf": "fig-arxiv-citing-patterns.pdf",
            "source": "Author stepwise R waffle/bar construction, using the verified frozen primary-class counts",
            "displayed_counts": {
                "high-energy physics": 233,
                "astrophysics": 180,
                "general relativity and quantum cosmology": 76,
                "physics": 48,
            },
            "classified_records": 551,
            "outside_display": 14,
            **PROFILES["fig-arxiv-citing-patterns"],
        }
    )
    for row in records:
        row.update(
            presentation_revision="0.1.6",
            text_font="Bundled STIX Two Text",
            emphasis="Medium",
            background="white",
        )
        row.setdefault("minimum_normal_text_pt", 9)
    (root / "reproduced/figure-registry.json").write_text(
        json.dumps(records, indent=2) + "\n"
    )
    return records

def export_impact(root, out):
    rates = pd.read_parquet(root / "inputs/avg_citations_result.parquet")
    positive = rates[rates.avg_citations_per_year.gt(0)].copy()
    target = positive.loc[positive.bibcode.eq(analysis.BH)].iloc[0]
    rng = np.random.default_rng(42)
    y_jitter = rng.uniform(-1.5, 1.5, len(positive))
    x_jitter = rng.uniform(-1.5, 1.5, len(positive))
    fig, ax = plt.subplots(figsize=(156 / 25.4, 90 / 25.4))
    ax.scatter(
        positive.avg_citations_per_year + x_jitter,
        positive.percentile + y_jitter,
        alpha=0.15,
        color="#FFA239",
        marker="2",
        linewidths=0.2,
        edgecolors="black",
        s=10,
    )
    ax.scatter(
        [target.avg_citations_per_year],
        [target.percentile],
        color="#E37434",
        s=50,
        marker="o",
        edgecolors="black",
        linewidths=0.5,
        label=analysis.BH,
        zorder=5,
    )
    ax.axvline(
        positive.avg_citations_per_year.mean(),
        color="#2cb869",
        linestyle="--",
        linewidth=1,
        label="Overall mean",
    )
    ax.axvline(
        positive.avg_citations_per_year.median(),
        color="#333333",
        linestyle="--",
        linewidth=1,
        label="Overall median",
    )
    ax.annotate(
        f"Bibcode: {analysis.BH}\nAvg cit/year: {target.avg_citations_per_year:.2f}\nPercentile: {target.percentile:.1f}",
        xy=(target.avg_citations_per_year, target.percentile),
        xytext=(0.51, 0.57),
        textcoords="axes fraction",
        fontsize=9,
        color="#333333",
        bbox={
            "boxstyle": "round",
            "facecolor": "white",
            "edgecolor": "#777777",
            "linewidth": 0.6,
        },
        arrowprops={"color": "#777777", "arrowstyle": "->", "linewidth": 0.6},
    )
    ax.set_xscale("symlog")
    ax.set_xlim(1, 550)
    ax.xaxis.set_major_formatter(
        ticker.FuncFormatter(lambda v, p: f"{int(v):,}" if v >= 1 else "")
    )
    ax.set_xticks([1, 2, 5, 10, 25, 50, 100, 200, 500])
    ax.tick_params(labelsize=9)
    ax.grid(True, linestyle="--", alpha=0.3)
    ax.set_title("Percentile vs. average citations per year", fontsize=11)
    ax.set_xlabel("Average citations per year (symlog)", fontsize=10)
    ax.set_ylabel("Percentile", fontsize=10)
    ax.legend(fontsize=9, loc="lower right", framealpha=0.95)
    fig.tight_layout(pad=0.7)
    path = save_pdf(fig, out / "fig-citation-impact.pdf", 156, 90)
    plt.close(fig)
    return {
        "figure": "fig-citation-impact",
        "pdf": path.name,
        "source": "Retained annual-rate table; original notebook scatter construction",
        "seed": 42,
        "jitter_order": "y followed by x, uniform -1.5 to 1.5; unchanged from retained notebook",
        "positive_rate_records": len(positive),
        "percentile_denominator": len(rates),
        "target_rate": float(target.avg_citations_per_year),
        "target_percentile": float(target.percentile),
        **PROFILES["fig-citation-impact"],
    }

# absolute_tiles
absolute_tiles <- make_waffle_tiles(arxiv_plot_data, n_cols = 13, normalized = FALSE)

# absolute_waffle_plot
absolute_waffle_plot <- plot_arxiv_waffle(absolute_tiles, title = NULL, subtitle = NULL, caption = NULL)

# normalized_tiles
normalized_tiles <- make_total_normalized_waffle_tiles(arxiv_plot_data, n_cols = 25, max_total_squares = 500)

# normalized_waffle_plot
normalized_waffle_plot <- plot_arxiv_single_waffle(normalized_tiles, title = "arXiv class resolution", subtitle = NULL, 
    caption = NULL)

# waffle_design
waffle_design <- "\nAAAAA\nBBCCC\n"

# arxiv_waffle_plot
arxiv_waffle_plot <- absolute_waffle_plot/two_panel_left + two_panel_right + plot_layout(design = waffle_design, 
    heights = c(1, 1), guides = "keep") + plot_annotation(title = "Citing patterns across arXiv classes", 
    subtitle = "Absolute counts and within-category shares of BH citers by arXiv category and subclass.") & theme(plot.title = element_text(size = 11, 
    face = "bold", family = "DAMADI STIX Two Text"), plot.subtitle = element_text(size = 9, family = "DAMADI STIX Two Text"))

# 3. Check field counts against the accepted revision, then export at final size.
expected <- c("high-energy physics"=233L,"astrophysics"=180L,"general relativity and quantum cosmology"=76L,"physics"=48L)
actual <- setNames(category_totals$total_n,as.character(category_totals$arxiv_category))
stopifnot(identical(as.integer(actual[names(expected)]),unname(expected)))
stopifnot(nrow(absolute_tiles)==537L,nrow(normalized_tiles)==500L)
arxiv_waffle_plot <- arxiv_waffle_plot & theme(text=element_text(family="DAMADI STIX Two Text",size=9),strip.text=element_text(size=9.5,vjust=0,margin=margin(t=0,r=1,b=4,l=1)),plot.title=element_text(size=11,face="bold"),plot.subtitle=element_text(size=9))
save_damadi_pdf(arxiv_waffle_plot,file.path(root,"reproduced/figures/fig-arxiv-citing-patterns.pdf"),width_mm=156,height_mm=(156 * 5.68 / 7.25))
write.csv(arxiv_plot_data,file.path(root,"reproduced/analysis/waffle-subclass-counts.csv"),row.names=FALSE)
write.csv(data.frame(arxiv_category=names(category_palette),color=unname(category_palette)),file.path(root,"reproduced/analysis/figure-category-palette.csv"),row.names=FALSE)
writeLines(capture.output(sessionInfo()),file.path(root,"verification/r-session-info.txt"))
writeLines(c("displayed_records: 537","classified_records: 551","outside_display: 14","font: bundled STIX Two Text; Medium emphasis","profile_mm: 156 x 122.2179310345 (author aspect 7.25:5.68)"),file.path(root,"verification/waffle-checks.txt"))
dev.off(metrics_device)
cat("STIX waffle PDF exported; four-field counts and tile totals verified.\n")

#!/usr/bin/env Rscript
# E.03: the author's ggplot specification, using verified revision values.
args <- commandArgs(trailingOnly = TRUE)
root <- normalizePath(if (length(args)) args[[1]] else getwd(), mustWork = TRUE)
stopifnot(file.exists(file.path(root, "config/workspace.json")))
suppressPackageStartupMessages({library(ggplot2);library(scales)})
if (!capabilities("aqua")) stop("The requested Quartz PDF export requires macOS R with Quartz support.")
source(file.path(root, "scripts/modules/quartz-fonts.R"))
register_damadi_quartz_fonts(
  file.path(root, "resources/fonts"),
  file.path(root, ".runtime/quartz-fonts"),
  file.path(root, "scripts/modules/register-quartz-fonts.c")
)
# Preserve the revised model predictions; only adapt column names for this plot.
verified <- read.csv(file.path(root, "reproduced/analysis/stratum-differences.csv"), check.names = FALSE)
stratum_plot_data <- verified
stratum_plot_data$diff_adj <- stratum_plot_data$difference
keep_categories <- c("astrophysics", "high-energy physics", "general relativity and quantum cosmology", "physics", "nuclear physics", "quantum physics")
stratum_plot_data$arxiv_category <- factor(stratum_plot_data$arxiv_category, levels = rev(keep_categories))
stratum_plot_data$year_bin <- factor(stratum_plot_data$year_bin, levels = c("2018–2019", "2020–2021", "2022–2023", "2024–2025"))
palette_data <- read.csv(file.path(root, "reproduced/analysis/figure-category-palette.csv"))
category_palette <- setNames(palette_data$color, palette_data$arxiv_category)
stopifnot(nrow(stratum_plot_data) == 12L, all(stratum_plot_data$n_bh >= 10), all(stratum_plot_data$n_nonbh >= 100))
stopifnot(identical(stratum_plot_data$diff_adj, verified$difference), all(as.character(stratum_plot_data$arxiv_category) %in% names(category_palette)))

stratum_plot <- ggplot(
  stratum_plot_data,
  aes(x = diff_adj, y = arxiv_category, size = n_bh, color = arxiv_category)
) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.3) +
  geom_point(alpha = 0.85) +
  facet_grid(~ year_bin) +
  scale_y_discrete(labels = scales::label_wrap(15)) +
  scale_x_continuous(
    labels = percent_format(accuracy = 1),
    breaks = seq(0.0, 0.10, by = 0.02)
  ) +
  coord_cartesian(xlim = c(0.0, 0.10)) +
  scale_color_manual(values = category_palette) +
  scale_size_continuous(range = c(5, 15), name = "# BH-citing papers") +
  guides(
    color = "none",
    size = guide_legend(
      title = "Review-citing papers:",
      position = "bottom",
      title.position = "left",
      nrow = 1
    )
  ) +
  theme_minimal(base_family = "STIX Two Text") +
  labs(
    x = "Adjusted probability difference (percentage points)",
    y = NULL,
    title = "Length-adjusted canon uptake among BH citers"
  ) +
  theme(
    panel.grid.minor.x = element_blank(),
    plot.title = element_text(size = 11, hjust = 0.5),
    axis.text.x = element_text(size = 8.5, hjust = 0.2),
    axis.text.y = element_text(size = 9, lineheight = 0.96),
    strip.text.x = element_text(size = 9.5, hjust = 0.5),
    strip.text.y = element_text(angle = 0),
    axis.title.x = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 10),
    legend.margin = margin(t = 0, r = 0, b = 0, l = 0)
  )

# Outer page padding keeps the last 10% tick inside the PDF; supplied settings above are retained.
stratum_plot <- stratum_plot + theme(plot.margin = margin(t = 5.5, r = 9, b = 5.5, l = 5.5))

quartz(
  type = "pdf",
  file = file.path(root, "reproduced/figures/fig-cocitation-strata.pdf"),
  width = 7,
  height = 4,
  family = "STIX Two Text",
  bg = "white"
)
print(stratum_plot)
dev.off()
# Keep the actual display data, scale values, and session as inspectable records.
built <- ggplot_build(stratum_plot)
write.csv(stratum_plot_data, file.path(root, "reproduced/analysis/stratum-plot-data.csv"), row.names = FALSE)
drawn <- built$data[[2]]
stopifnot(nrow(drawn) == nrow(stratum_plot_data), all(is.finite(drawn$size)))
stopifnot(isTRUE(all.equal(drawn$x, stratum_plot_data$diff_adj)))
expected_sizes <- 5 + 10 * sqrt((stratum_plot_data$n_bh - min(stratum_plot_data$n_bh)) / diff(range(stratum_plot_data$n_bh)))
stopifnot(isTRUE(all.equal(drawn$size, expected_sizes)))
write.csv(data.frame(n_bh = stratum_plot_data$n_bh, diff_adj = stratum_plot_data$diff_adj, diameter_mm = drawn$size), file.path(root, "verification/stratum-size-scale.csv"), row.names = FALSE)
writeLines(capture.output(sessionInfo()), file.path(root, "verification/stratum-r-session-info.txt"))
writeLines(c("backend: Quartz PDF", "width_in: 7", "height_in: 4", "font: bundled STIX Two Text (process-local CoreText registration)", "plot: author-supplied ggplot specification, unchanged", "data: verified revision stratum predictions; no refit"), file.path(root, "verification/stratum-export.txt"))
cat("Author-specified ggplot exported with Quartz at 7 x 4 inches.\n")

def write_run_record(ctx, verification=None):
    record = {
        "entry_code": "E.03--MS-DA-CR",
        "edition": "0.1.0",
        "companion_status": "unreleased normalized update",
        "inputs": ctx["checked_inputs"],
        "software": {
            n: metadata.version(n)
            for n in [
                "numpy",
                "pandas",
                "pyarrow",
                "scipy",
                "statsmodels",
                "matplotlib",
            ]
        },
        "summary_file": "analysis/bibliometric-summary.json",
        "verification": verification,
        "scope": "Reconstruction from retained revision inputs; no live retrieval or new analytical choices",
    }
    path = ctx["output_root"] / "run-record.json"
    path.write_text(json.dumps(record, indent=2) + "\n")
    return record

if os.environ.get("DAMADI_ANALYSIS_ONLY") == "1":
    print("Numerical-only mode: approved PDF export skipped.")
else:
    figures = plotting.export_all(root)
    for figure in figures:
        print(figure["pdf"], figure["width_mm"], figure["height_mm"])
record = analysis.write_run_record(state, verification=report)
assert record["verification"]["passed"]
print("Checked reconstruction complete. Supplied evidence remains unchanged.")