Diverging Paths to
Dark-Matter Discovery

def input_checks(root):
    root = Path(root)
    load_local(
        root / "code/materialize_inputs.py", "diverging_paths_materialize"
    ).materialize(root)
    pins = json.loads((root / "checks/input-hashes.json").read_text())
    for rel, digest in pins.items():
        assert hashlib.sha256((root / rel).read_bytes()).hexdigest() == digest, rel
    release = json.loads(
        (root / "provenance/v3.0.0/release_manifest.json").read_text()
    )
    members = 0
    for archive in release["network_archives"]:
        path = root / archive["path"]
        assert path.stat().st_size == archive["bytes"]
        assert hashlib.sha256(path.read_bytes()).hexdigest() == archive["sha256"]
        with zipfile.ZipFile(path) as z:
            assert set(z.namelist()) == set(archive["members"])
            for name, digest in archive["members"].items():
                assert hashlib.sha256(z.read(name)).hexdigest() == digest, name
            members += len(z.namelist())
    return write_report(
        root,
        "input-integrity.json",
        {
            "status": "PASS",
            "pinned_data_products": len(pins),
            "network_archives": len(release["network_archives"]),
            "archive_members": members,
            "python": platform.python_version(),
            "scope": "Integrity only; no holdout annotations read",
        },
    )

"""Restore exact frozen tables from lossless transport parts."""

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


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


def materialize(root):
    root = Path(root).resolve()
    data = root / "data"
    spec = json.loads((data / "INPUT-STORAGE.json").read_text())
    restored = []
    for name, entry in spec["tables"].items():
        target = data / name
        if target.exists():
            if (
                target.stat().st_size != entry["bytes"]
                or sha256(target) != entry["sha256"]
            ):
                raise ValueError(
                    "Existing input has a different identity; preserve it and use a clean checkout: "
                    + name
                )
            continue
        temporary = target.with_suffix(target.suffix + ".assembling")
        try:
            with temporary.open("wb") as output:
                for part in entry["parts"]:
                    path = data / part["file"]
                    if (
                        path.stat().st_size != part["bytes"]
                        or sha256(path) != part["sha256"]
                    ):
                        raise ValueError(
                            "Transport-part identity differs: " + part["file"]
                        )
                    with path.open("rb") as source:
                        shutil.copyfileobj(source, output)
            if (
                temporary.stat().st_size != entry["bytes"]
                or sha256(temporary) != entry["sha256"]
            ):
                raise ValueError("Restored table identity differs: " + name)
            os.replace(temporary, target)
            restored.append(name)
        finally:
            if temporary.exists():
                temporary.unlink()
    return restored


if __name__ == "__main__":
    print(
        json.dumps(
            {"restored": materialize(Path(__file__).resolve().parents[1])}, indent=2
        )
    )

reports["input_checks"] = input_checks(ROOT)
reports["input_checks"]

from pathlib import Path
import sys, os

ROOT = Path.cwd().resolve()
if ROOT.name == "notebooks":
    ROOT = ROOT.parent
assert (ROOT / "config/workspace.json").is_file()
sys.path.insert(0, str(ROOT / "code/modules"))
from workflow import (
    input_checks,
    corpus_checks,
    candidate_checks,
    network_checks,
    figure_checks,
    write_report,
)

reports = {}

def corpus_checks(root):
    root = Path(root)
    data = root / "data"
    out = root / "reproduced/tables"
    out.mkdir(parents=True, exist_ok=True)
    papers = pd.read_parquet(
        data / "papers.parquet", columns=["bibcode", "year", "abstract"]
    )
    papers = papers[papers.year <= 2025]
    assert len(papers) == 201627 and papers.bibcode.nunique() == len(papers)
    classes = pd.read_parquet(data / "paper_arxiv_classes.parquet").sort_values(
        ["bibcode", "class_pos"]
    )
    primary = classes.drop_duplicates("bibcode").merge(
        papers[["bibcode", "year"]], on="bibcode", validate="one_to_one"
    )
    assert len(primary) == 143390
    rebuilt = (
        primary.groupby(["year", "arxiv_class", "arxiv_category"])
        .size()
        .rename("n_papers")
        .reset_index()
    )
    frozen = pd.read_parquet(
        data / "analysis/classification_primary_yearly.parquet"
    )
    joined = rebuilt.merge(
        frozen,
        on=["year", "arxiv_class", "arxiv_category"],
        how="outer",
        suffixes=("_rebuilt", "_frozen"),
        indicator=True,
    )
    assert (
        joined["_merge"].eq("both").all()
        and joined.n_papers_rebuilt.eq(joined.n_papers_frozen).all()
    )
    volume = (
        papers.groupby("year")
        .size()
        .rename("n_papers_rebuilt")
        .reset_index()
        .merge(
            pd.read_parquet(data / "analysis/paper_volume_yearly.parquet"),
            on="year",
            validate="one_to_one",
        )
    )
    assert volume.n_papers_rebuilt.eq(volume.n_papers).all()
    categories = (
        primary.groupby("arxiv_category")
        .size()
        .sort_values(ascending=False)
        .rename("papers")
        .reset_index()
    )
    counts = (
        primary.groupby("arxiv_class")
        .size()
        .sort_values(ascending=False)
        .rename("papers")
        .reset_index()
    )
    for name, table in [
        ("primary-categories", categories),
        ("primary-classes", counts),
    ]:
        table["percent"] = 100 * table.papers / len(primary)
        table.to_csv(out / (name + ".csv"), index=False)
        lines = [
            r"\begin{tabular}{lrr}",
            r"\toprule",
            "Primary assignment & Papers & Share " + r"\\",
            r"\midrule",
        ]
        for _, row in table.head(
            4 if name.endswith("categories") else 9
        ).iterrows():
            label = str(row.iloc[0]).replace("_", r"\_")
            lines.append(
                f"{label} & {int(row.papers):,} & {row.percent:.1f}" + r"\% \\"
            )
        (out / (name + ".tex")).write_text(
            "\n".join(lines + [r"\bottomrule", r"\end{tabular}"]) + "\n"
        )
    expected = {
        "astrophysics": 78673,
        "high-energy physics": 39191,
        "general relativity and quantum cosmology": 14668,
        "physics": 6427,
    }
    actual = categories.set_index("arxiv_category").papers.to_dict()
    assert all(actual.get(k) == v for k, v in expected.items())
    return write_report(
        root,
        "corpus-checks.json",
        {
            "status": "PASS",
            "records": len(papers),
            "classified_records": len(primary),
            "primary_class_year_cells_rebuilt": len(rebuilt),
            "volume_years_rebuilt": len(volume),
            "main_table_counts_match": True,
            "largest_primary_class": counts.iloc[0, 0],
        },
    )

reports["corpus_checks"] = corpus_checks(ROOT)
reports["corpus_checks"]

def candidate_checks(root):
    root = Path(root)
    d = root / "data"
    raw = pd.read_parquet(d / "dm_model_candidates_long.parquet")
    raw = raw[raw.year.between(1995, 2025)]
    assert len(raw) == 27680 and raw.bibcode.nunique() == 22247
    tracked = (
        (d / "analysis/tracked_candidate_species.txt").read_text().splitlines()
    )
    assert len(tracked) == 13
    pairs = raw[
        raw.SpeciesLabel.isin(tracked)
        & raw.year.between(1995, 2025)
        & raw.arxiv_category.isin(["astrophysics", "high-energy physics"])
    ].drop_duplicates(["year", "arxiv_category", "bibcode", "SpeciesLabel"])
    rebuilt = (
        pairs.groupby(["year", "SpeciesLabel", "arxiv_category"])
        .size()
        .unstack(fill_value=0)
        .reset_index()
    )
    ref = pd.read_parquet(d / "analysis/candidate_field_yearly.parquet")
    merged = ref.merge(rebuilt, on=["year", "SpeciesLabel"], how="left").fillna(
        {"astrophysics": 0, "high-energy physics": 0}
    )
    assert (merged.astro_mentions == merged["astrophysics"]).all() and (
        merged.hep_mentions == merged["high-energy physics"]
    ).all()
    # Reconstruct smoothed proportions using the retained 0.5 per-family pseudocount.
    n = 13
    astro = (merged.astro_mentions + 0.5) / (
        merged.groupby("year").astro_mentions.transform("sum") + 0.5 * n
    )
    hep = (merged.hep_mentions + 0.5) / (
        merged.groupby("year").hep_mentions.transform("sum") + 0.5 * n
    )
    assert np.allclose(astro, merged.astro_share) and np.allclose(
        hep, merged.hep_share
    )
    assert np.allclose(np.log2(hep / astro), merged.log2_hep_astro_ratio)
    assert (merged.meets_minimum_count == (merged.total_mentions >= 20)).all()
    merged[
        [
            "year",
            "SpeciesLabel",
            "astro_mentions",
            "hep_mentions",
            "astro_share",
            "hep_share",
            "log2_hep_astro_ratio",
        ]
    ].to_csv(
        root / "reproduced/tables/candidate-field-reconstructed.csv", index=False
    )
    return write_report(
        root,
        "candidate-checks.json",
        {
            "status": "PASS",
            "candidate_tag_rows_1995_2025": len(raw),
            "candidate_papers": int(raw.bibcode.nunique()),
            "tracked_families": len(tracked),
            "candidate_year_cells": len(merged),
            "paper_family_deduplication": True,
            "pseudocount": 0.5,
            "display_threshold": 20,
            "bootstrap_intervals": "Retained and hash-verified; not freshly resampled",
            "validation": "No independent accuracy estimate; no untouched holdout answers inspected",
        },
    )

reports["candidate_checks"] = candidate_checks(ROOT)
reports["candidate_checks"]

def reconstruct_graph(nodes, edges):
    graph = nx.Graph()
    for row in edges:
        weight = 0.5 * (
            int(row["astrophysics_paper_count"])
            / int(row["astrophysics_all_dm_papers"])
            + int(row["hep_paper_count"]) / int(row["hep_all_dm_papers"])
        )
        check_close(
            weight, row["mean_field_cooc_prevalence"], "Equal-field edge weight"
        )
        distance = float(row[DISTANCE])
        check_close(distance, 1 / weight, "Inverse-prevalence distance")
        source, target = row["Source"], row["Target"]
        if graph.has_edge(source, target) or source == target:
            raise ValueError("Duplicate edge or self-loop")
        graph.add_edge(source, target, **{DISTANCE: distance})
    if set(graph) != {row["Id"] for row in nodes}:
        raise ValueError("Canonical node and edge sets disagree")
    astro, hep = [], []
    for row in nodes:
        candidate = row["node_type"] == "candidate"
        ratio = float(row["log10_astro_hep_prevalence_ratio"])
        for name, condition, endpoints in [
            ("is_astro_anchor", not candidate and ratio >= math.log10(2), astro),
            ("is_hep_anchor", not candidate and ratio <= -math.log10(2), hep),
        ]:
            if (row[name] == "True") != condition:
                raise ValueError(f'Endpoint rule mismatch: {row["Id"]}, {name}')
            if condition:
                endpoints.append(row["Id"])
    if not astro or not hep:
        raise ValueError(
            "No endpoints; archive fallback would require explicit treatment"
        )
    return graph, astro, hep

def efficiency(graph, sources, targets):
    """Mean reciprocal distance; disconnected pairs remain in the denominator."""
    sources, targets = sorted(sources), sorted(targets)
    if not sources or not targets or not set(sources + targets) <= set(graph):
        raise ValueError("Endpoint sets must be nonempty and present in the graph")
    total, reachable = 0.0, 0
    for source in sources:
        lengths = nx.single_source_dijkstra_path_length(
            graph, source, weight=DISTANCE
        )
        for target in targets:
            distance = lengths.get(target)
            if distance is not None and distance > 0:
                total += 1.0 / distance
                reachable += 1
    pairs = len(sources) * len(targets)
    return total / pairs, reachable / pairs

def removal(graph, sources, targets, node, baseline):
    if node in set(sources) | set(targets):
        raise ValueError(
            "Removal requires a non-anchor node, keeping endpoints fixed"
        )
    reduced = graph.copy()
    reduced.remove_node(node)
    after, reachability = efficiency(reduced, sources, targets)
    return after, reachability, max(0.0, (baseline - after) / baseline)

def spectral_scores(rows):
    """Role-specific HITS via the leading left singular vector, max scaled."""
    candidates = sorted({r["Source"] for r in rows})
    concepts = sorted({r["Target"] for r in rows})
    ci, ti = {v: i for i, v in enumerate(candidates)}, {
        v: i for i, v in enumerate(concepts)
    }
    matrix = np.zeros((len(candidates), len(concepts)))
    seen = set()
    for row in rows:
        edge = row["Source"], row["Target"]
        if edge in seen:
            raise ValueError(f"Duplicate candidate-concept edge {edge}")
        seen.add(edge)
        weight = int(row["paper_count"]) / int(row["candidate_bearing_papers"])
        check_close(
            weight, row["cooc_doc_prevalence_candidate_papers"], "HITS weight"
        )
        matrix[ci[edge[0]], ti[edge[1]]] = weight
    left, _, _ = np.linalg.svd(matrix, full_matrices=False)
    values = np.abs(left[:, 0])
    values /= values.max()
    return dict(zip(candidates, map(float, values)))

def verify(release_root, pins_path):
    pins = json.loads(pins_path.read_text())
    manifest_path = release_root / "provenance/v3.0.0/release_manifest.json"
    if not manifest_path.is_file():
        manifest_path = release_root / "release_manifest.json"
    if sha256(manifest_path) not in pins["accepted_release_manifest_sha256"]:
        raise ValueError("Release manifest does not match the pinned v3.0.0 input")
    manifest = json.loads(manifest_path.read_text())
    entry = next(
        x for x in manifest["network_archives"] if x["path"] == pins["archive_path"]
    )
    archive_path = release_root / entry["path"]
    digest = sha256(archive_path)
    if digest != pins["archive_sha256"] or digest != entry["sha256"]:
        raise ValueError("Frozen archive SHA-256 mismatch")
    result = {
        "status": "passed",
        "python": platform.python_version(),
        "networkx": nx.__version__,
        "numpy": np.__version__,
        "archive_sha256": digest,
        "release_manifest_sha256": sha256(manifest_path),
        "scope": "Conditional recalculation from frozen graphs and count tables; not extraction validation",
        "windows": [],
        "common_endpoints": [],
    }
    with zipfile.ZipFile(archive_path) as archive:
        if set(archive.namelist()) != set(entry["members"]):
            raise ValueError("Archive member list differs from release manifest")
        for name, expected in entry["members"].items():
            if hashlib.sha256(archive.read(name)).hexdigest() != expected:
                raise ValueError(f"Archive member SHA-256 mismatch: {name}")
        result["archive_members_verified"] = len(entry["members"])
        nodes = grouped(read_rows(archive, BROKERAGE + "nodes_by_window.csv"))
        edges = grouped(read_rows(archive, BROKERAGE + "edges_by_window.csv"))
        trajectories = grouped(
            read_rows(archive, BROKERAGE + "wimp_brokerage_trajectory.csv")
        )
        removals = grouped(
            read_rows(archive, BROKERAGE + "node_removal_metrics.csv")
        )
        importance_edges = grouped(
            [
                r
                for r in read_rows(archive, IMPORTANCE + "edges_by_window.csv")
                if r["evidence_scope"] == "abstract"
            ]
        )
        importance_nodes = grouped(
            read_rows(archive, IMPORTANCE + "candidate_metrics_by_window.csv")
        )
        importance_trajectory = grouped(
            read_rows(archive, IMPORTANCE + "wimp_importance_trajectory.csv")
        )
        if not (
            set(nodes)
            == set(edges)
            == set(trajectories)
            == set(importance_edges)
            == set(importance_nodes)
            == set(importance_trajectory)
        ):
            raise ValueError("Window/definition groups differ across input tables")
        stable = {}
        removal_checks = 0
        for key in sorted(nodes):
            definition, kind, window = key
            graph, astro, hep = reconstruct_graph(nodes[key], edges[key])
            if kind == "stable":
                suffix = "" if definition == "published_umbrella" else "robustness/"
                name = (
                    BROKERAGE
                    + f'graphs/{suffix}{definition}--{window.replace("-", "_")}.graphml'
                )
                saved = nx.read_graphml(io.BytesIO(archive.read(name)))
                if set(saved) != set(graph) or {
                    frozenset(e) for e in saved.edges
                } != {frozenset(e) for e in graph.edges}:
                    raise ValueError(f"GraphML topology mismatch: {key}")
                for a, b, attrs in graph.edges(data=True):
                    check_close(
                        attrs[DISTANCE], saved[a][b][DISTANCE], "GraphML distance"
                    )
                # Keep the archived GraphML floats/order for exact shortest-path tie handling.
                graph = saved
            scores = nx.betweenness_centrality_subset(
                graph, sorted(astro), sorted(hep), normalized=True, weight=DISTANCE
            )
            candidates = {
                r["Id"]: scores[r["Id"]]
                for r in nodes[key]
                if r["node_type"] == "candidate"
            }
            for row in nodes[key]:
                check_close(
                    scores[row["Id"]],
                    row["cross_field_betweenness_prevalence"],
                    f'Betweenness {key} {row["Id"]}',
                )
            focal = FOCAL[definition]
            archived = trajectories[key][0]
            baseline, reachable = efficiency(graph, astro, hep)
            after, after_reachable, loss = removal(
                graph, astro, hep, focal, baseline
            )
            check_close(
                loss,
                archived["removal_effect_efficiency_prevalence"],
                f"WIMP removal {key}",
            )
            check_close(
                scores[focal],
                archived["cross_field_betweenness_prevalence"],
                f"WIMP trajectory {key}",
            )
            for row in removals[key]:
                if kind != "stable" and row["node_id"] != focal:
                    continue
                after_node, reach_node, loss_node = removal(
                    graph, astro, hep, row["node_id"], baseline
                )
                for value, column in [
                    (baseline, "baseline_efficiency"),
                    (after_node, "removed_efficiency"),
                    (reachable, "baseline_reachable_share"),
                    (reach_node, "removed_reachable_share"),
                    (loss_node, "relative_efficiency_loss"),
                ]:
                    check_close(
                        value, row[column], f'{column} {key} {row["node_id"]}'
                    )
                removal_checks += 1
            hits = spectral_scores(importance_edges[key])
            for row in importance_nodes[key]:
                check_close(
                    hits[row["Id"]],
                    row["hits_score_prevalence"],
                    f'HITS {key} {row["Id"]}',
                )
            hits_pct = percentile(hits, focal)
            check_close(
                hits_pct,
                importance_trajectory[key][0]["hits_candidate_percentile"],
                f"HITS percentile {key}",
            )
            result["windows"].append(
                {
                    "definition": definition,
                    "kind": kind,
                    "window": window,
                    "nodes": len(graph),
                    "edges": graph.number_of_edges(),
                    "astro_endpoints": len(astro),
                    "hep_endpoints": len(hep),
                    "wimp_betweenness": scores[focal],
                    "wimp_candidate_rank": 1
                    + sum(v > scores[focal] for v in candidates.values()),
                    "candidate_count": len(candidates),
                    "wimp_removal_effect": loss,
                    "wimp_hits_score": hits[focal],
                    "wimp_hits_percentile": hits_pct,
                }
            )
            if kind == "stable" and definition == "published_umbrella":
                stable[window] = graph, set(astro), set(hep), candidates
            print(f"Checked {definition} {kind} {window}", flush=True)
        common_astro = set.intersection(*(v[1] for v in stable.values()))
        common_hep = set.intersection(*(v[2] for v in stable.values()))
        for row in read_rows(
            archive, BROKERAGE + "common_endpoint_sensitivity.csv"
        ):
            graph, _, _, candidates = stable[row["window"]]
            scores = nx.betweenness_centrality_subset(
                graph,
                sorted(common_astro),
                sorted(common_hep),
                normalized=True,
                weight=DISTANCE,
            )
            baseline, reachable = efficiency(graph, common_astro, common_hep)
            focal = FOCAL["published_umbrella"]
            after, after_reachable, loss = removal(
                graph, common_astro, common_hep, focal, baseline
            )
            rank = 1 + sum(scores[n] > scores[focal] for n in candidates)
            for value, column in [
                (len(common_astro), "common_astro_endpoint_count"),
                (len(common_hep), "common_hep_endpoint_count"),
                (scores[focal], "common_endpoint_betweenness"),
                (rank, "common_endpoint_candidate_rank"),
                (baseline, "common_endpoint_baseline_efficiency"),
                (after, "common_endpoint_removed_efficiency"),
                (loss, "common_endpoint_removal_effect"),
                (reachable, "common_endpoint_baseline_reachable_share"),
                (after_reachable, "common_endpoint_removed_reachable_share"),
            ]:
                check_close(
                    value, row[column], f'Common endpoints {row["window"]} {column}'
                )
            result["common_endpoints"].append(
                {
                    "window": row["window"],
                    "astro_endpoints": len(common_astro),
                    "hep_endpoints": len(common_hep),
                    "betweenness": scores[focal],
                    "candidate_rank": rank,
                    "removal_effect": loss,
                }
            )
        result["node_removal_rows_recalculated"] = removal_checks
        result["betweenness_nodes_checked"] = sum(map(len, nodes.values()))
        result["hits_candidates_checked"] = sum(map(len, importance_nodes.values()))
    return result

def network_checks(root):
    root = Path(root)
    mod = load_local(
        root / "code/network-reconstruction/recompute_metrics.py",
        "damadi_network_recompute",
    )
    report = mod.verify(root, root / "code/network-reconstruction/inputs.json")
    assert report["status"] == "passed" and len(report["windows"]) == 57
    path = root / "reproduced/network-metrics.json"
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(json.dumps(report, indent=2, allow_nan=False) + "\n")
    return {
        "status": "PASS",
        "windows": len(report["windows"]),
        "nodes_checked": report["betweenness_nodes_checked"],
        "removal_rows": report["node_removal_rows_recalculated"],
        "scope": report["scope"],
    }

reports["network_checks"] = network_checks(ROOT)
reports["network_checks"]

def figure_checks(root, render=True):
    root = Path(root)
    if render:
        subprocess.run(
            ["Rscript", str(root / "code/04-build-figures.R")], cwd=root, check=True
        )
        mod = load_local(
            root / "code/modules/network_graphics.py", "damadi_network_graphics"
        )
        mod.export_all(root)
    from pypdf import PdfReader

    files = sorted((root / "reproduced/stepwise-figures").glob("*.pdf")) + sorted(
        (root / "reproduced/figures").glob("*.pdf")
    )
    assert len(files) == 9, len(files)
    for file in files:
        assert len(PdfReader(file).pages) == 1
    config = json.loads((root / "config/main-figures.json").read_text())
    assert (
        hashlib.sha256(
            (root / "code/02-main-figures-stepwise.R").read_bytes()
        ).hexdigest()
        == config["active_script_sha256"]
    )
    for name, dimensions in config["figures"].items():
        page = PdfReader(root / "reproduced/stepwise-figures" / name).pages[0]
        actual = [
            float(page.mediabox.width) * 25.4 / 72,
            float(page.mediabox.height) * 25.4 / 72,
        ]
        assert np.allclose(actual, dimensions, atol=0.02), (name, actual)
    for name, digest in config["table_sha256"].items():
        assert (
            hashlib.sha256(
                (root / "reproduced/tables" / name).read_bytes()
            ).hexdigest()
            == digest
        ), name
    return write_report(
        root,
        "figure-checks.json",
        {
            "status": "PASS",
            "pdfs": len(files),
            "statistical_renderings": 7,
            "network_presentations": 2,
            "network_geometry": "Unchanged",
            "main_figure_source_sha256": config["source_script_sha256"],
            "main_figure_dimensions_mm": config["figures"],
            "main_tables_byte_identical_to_retained_fragments": True,
            "authority": "Main Figures 1–3 regenerated from the user-designated script at final arXiv dimensions, with the preprint class colors fixed explicitly and bundled STIX typography. Four SI manuscript figures remain retained preprint assets; separate SI R companions remain available.",
        },
    )

composition_plot_wide_guide <- (
  (broad_plot +
     one_row_guides +
     theme(legend.position = "top", legend.justification = "left")) +
    (paper_plot_wide +
       one_row_guides +
       theme(legend.position = "top", legend.justification = "left")) +
    (papers_strip + theme(legend.position = "none")) +
    guide_area() +
    plot_layout(
      design = "
      DD
      AB
      CC
    ",
      widths = c(1, 1.55),
      heights = c(0.5, 9.3, 0.7),
      guides = "keep",
      axes = "collect"
    ) +
    plot_annotation(
      title = "Frequency of primary research classifications",
      subtitle = "arXiv class and category − 2010-2025",
    )
) &
  theme(
    legend.box = "horizontal",
    legend.direction = "horizontal",
    legend.byrow = TRUE,
    #legend.justification = "center",
    legend.key.spacing = grid::unit(1, "pt"),
    legend.spacing.x = grid::unit(3, "pt"),
    legend.text = element_text(size = 9, margin = margin(l = 3, r = 10)),
    legend.margin = margin(0, 0, 0, 0),
    
    plot.title = element_text(
      size = 11, family = REPORT_FONT,
      hjust = 0.05, vjust = 0.0
    ),
    plot.subtitle = element_text(
      size = 9, family = REPORT_FONT,
      hjust = 0.045, vjust = 0.0
    ),
    plot.tag = element_text(size = 7, family = REPORT_FONT),
    plot.tag.position = c(0.01, 0.98),
    axis.title.x = element_blank(),
    plot.margin = margin(0, 7, 1, 0)
  )

composition_plot_wide_guide <- composition_plot_wide_guide &
  theme(
    plot.background = element_blank(),
    panel.background = element_blank(),
    legend.background = element_blank(),
    legend.box.background = element_blank()
  )

if (interactive()) print(composition_plot_wide_guide)
if (SAVE_OUTPUTS) {
  save_stepwise_pdf(composition_plot_wide_guide, "fig-research-classifications.pdf", 165, 90, bg = "transparent")
}

MIN_TOTAL_MENTIONS <- 20L

candidate_field_plot <- candidate_field |>
  mutate(SpeciesLabel = stringr::str_replace(SpeciesLabel, fixed("Sterile nu"), "Sterile ν"))

stopifnot(all(
  candidate_field_plot$meets_minimum_count ==
    (candidate_field_plot$total_mentions >= MIN_TOTAL_MENTIONS)
))

species_order <- candidate_field_plot |>
  filter(year >= 1995, year <= 2025, total_mentions >= MIN_TOTAL_MENTIONS) |>
  group_by(SpeciesLabel) |>
  summarise(
    mean_ratio = weighted.mean(log2_hep_astro_ratio, w = total_mentions, na.rm = TRUE),
    .groups = "drop"
  ) |>
  arrange(mean_ratio) |>
  pull(SpeciesLabel)

log_ratio_data <- candidate_field_plot |>
  filter(year >= 1995, year <= 2025, SpeciesLabel %in% species_order) |>
  mutate(
    SpeciesLabel = factor(SpeciesLabel, levels = species_order),
    shown_ratio = if_else(total_mentions >= MIN_TOTAL_MENTIONS, log2_hep_astro_ratio, NA_real_)
  )

div <- c(
  "#00429d", "#2b57a7", "#426cb0", "#5681b9", "#6997c2",
  "#7daeca", "#93c4d2", "#abdad9", "#caefdf", "#ffffe0",
  "#ffe2ca", "#ffc4b4", "#ffa59e", "#f98689", "#ed6976",
  "#dd4c65", "#ca2f55", "#b11346", "#93003a"
)


log_ratio_plot <- log_ratio_data |>
  ggplot(aes(year, SpeciesLabel, fill = shown_ratio)) +
  geom_tile(color = "white", linewidth = 0.15, alpha = 0.85) +
  scale_x_continuous(
    breaks = seq(1995, 2025, 2),
    expand = expansion(mult = c(0.01, 0))
  ) +
  scale_fill_gradientn(
    colours = div,
    name = "Share ratio",
    limits = c(-log2(10), log2(10)),
    breaks = c(-log2(10), -log2(3), 0, log2(3), log2(10)),
    labels = c("10× Astro", "3× Astro", "1×", "3× HEP", "10× HEP"),
    oob = scales::squish,
    na.value = NA,
    guide = guide_colorbar(
      direction = "horizontal",
      title.position = "top",
      title.hjust = 0,
      title.vjust = -0.3,
      barwidth = grid::unit(70, "mm"),
      barheight = grid::unit(3.5, "mm"),
      ticks.colour = "white"
    )
  ) +
  labs(
    x = NULL,
    y = NULL,
    title = "Relative prominence of candidates",
    subtitle = "Log2 ratio of HEP to astrophysics shares, 1995-2025"
  ) +
  theme_minimal(base_family = REPORT_FONT, base_size = 11) +
  theme(
    plot.title = element_text(size = 11, family = REPORT_FONT, hjust = 0),
    plot.subtitle = element_text(
      size = 8, family = REPORT_FONT, hjust = 0, vjust = 2.0),
    plot.margin = margin(2, 5, 5, 0),
    panel.grid.major = element_line(color = "#282828", linewidth = 0.05),
    panel.grid.minor = element_blank(),
    legend.margin = margin(-7, 0, 0, 0),
    legend.text = element_text(size = 8, family = REPORT_FONT, color = "black"),
    legend.title = element_text(size = 8, family = REPORT_FONT, color = "black"),
    legend.position = "bottom",
    plot.background = element_blank(),
    panel.background = element_blank(),
    legend.background = element_blank(),
    legend.box.background = element_blank()
  )

if (interactive()) print(log_ratio_plot)
if (SAVE_OUTPUTS) {
  save_stepwise_pdf(log_ratio_plot, "fig-candidate-prominence.pdf", 156, 80, bg = "transparent")
}

tracked_species <- stringr::str_replace(tracked_species, fixed("Sterile nu"), "Sterile ν")
candidate_norm <- candidates |>
  mutate(SpeciesLabel = stringr::str_replace(SpeciesLabel, fixed("Sterile nu"), "Sterile ν")) |>
  filter(
    year >= 1995,
    year <= 2025,
    arxiv_category %in% names(FIELD_LABELS),
    SpeciesLabel %in% tracked_species
  ) |>
  distinct(year, arxiv_category, bibcode, SpeciesLabel) |>
  count(year, arxiv_category, SpeciesLabel, name = "n") |>
  group_by(year, arxiv_category) |>
  mutate(share_tracked_candidates = n / sum(n)) |>
  ungroup() |>
  mutate(
    field = recode(arxiv_category, !!!FIELD_LABELS),
    SpeciesLabel = factor(SpeciesLabel, levels = tracked_species)
  )

candidate_colors <- setNames(
  khroma::color("batlow", type = "qualitative", reverse = TRUE)(length(tracked_species)),
  tracked_species
)

p_norm <- candidate_norm |>
  ggplot(aes(year, share_tracked_candidates, fill = SpeciesLabel)) +
  geom_area(position = "fill", alpha = 0.89, linewidth = 0.15, color = "white") +
  facet_wrap(~ field, ncol = 2, scales = "fixed") +
  scale_fill_manual(values = candidate_colors, name = NULL, drop = FALSE) +
  scale_x_continuous(
    breaks = seq(1995, 2025, 5),
    expand = expansion(mult = c(0.01, 0))
  ) +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    expand = expansion(mult = c(0.01, 0))
  ) +
  labs(x = NULL, 
       title = "Share of paper-candidate-family mentions", 
       subtitle = "Normalized to total mentions in each field, 1995-2025",
       y = NULL) +
  theme_minimal(base_size = 11, base_family = REPORT_FONT) +
  theme(
    panel.spacing = grid::unit(1.5, "lines"),
    plot.margin = margin(0, 0, 0, 0),
    strip.background = element_rect(color = NA, fill = NA),
    strip.text = element_text(color = "#141414", size = 10, family = REPORT_FONT),
    strip.text.x.bottom = element_text(hjust = 0.5, vjust = 0.5),
    axis.ticks.x = element_blank(),
    legend.position = "bottom",
    legend.key.spacing.y = grid::unit(0.2, "lines")
  ) +
  guides(fill = guide_legend(
    nrow = 1,
    byrow = TRUE,
    title.position = "bottom",
    title.hjust = 0.5,
    label.position = "bottom",
    override.aes = list(alpha = 1, size = 1.5, shape = 15),
    keywidth = grid::unit(5, "pt"),
    keyheight = grid::unit(10, "pt")
  ))

candidate_robustness <- p_norm +
  patchwork::plot_layout(widths = 1, axes = "collect", guides = "keep") &
  theme(
    plot.title = element_text(size = 10, color = "black", hjust = 0, vjust = -2.8),
    plot.subtitle = element_text(size = 9, color = "black", hjust = 0, vjust = -1.0),
    plot.margin = margin(2, 12, 5, 0),
    legend.margin = margin(-5, 0, 0, 0),
    legend.position = "bottom",
    legend.box = "horizontal",
    #legend.box.just = "left",
    panel.grid.major = element_line(color = "black", linewidth = 0.15),
    panel.grid.minor = element_blank(),
    plot.background = element_blank(),
    panel.background = element_blank(),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.key.spacing = grid::unit(2.0, "pt"),
    legend.text = element_text(size = 7.7, margin = margin(t = 2.5, l = 0.7, r = 0.7, unit = "pt")),
    legend.title = element_text(size = 9, family = BASE_FONT, color = "black", hjust = 0)
  )

if (interactive()) print(candidate_robustness)
if (SAVE_OUTPUTS) {

def normalize_one(source, destination, font_dir):
    source = Path(source)
    destination = Path(destination)
    font_dir = Path(font_dir)
    original = fitz.open(source)
    page = original[0]
    assert len(original) == 1 and page.rotation == 0
    overlays = []
    for block in page.get_text("dict")["blocks"]:
        for line in block.get("lines", []):
            assert abs(line["dir"][1]) < 1e-8 and line["dir"][0] > 0
            for span in line["spans"]:
                font = span["font"]
                face = (
                    "math"
                    if "Math" in font
                    else (
                        "italic"
                        if "Italic" in font
                        else (
                            "medium"
                            if "IBM" in font or "Bold" in font or "Medium" in font
                            else "regular"
                        )
                    )
                )
                overlays.append(
                    {
                        "text": span["text"],
                        "origin": span["origin"],
                        "bbox": span["bbox"],
                        "size": span["size"],
                        "color": span["color"],
                        "face": face,
                        "source_font": font,
                    }
                )
    reader = PdfReader(source)
    writer = PdfWriter()
    writer.clone_document_from_reader(reader)
    fingerprints = []
    visited = set()

    def clean(stream, resources, label):
        key = id(stream)
        if key in visited:
            return
        visited.add(key)
        content = ContentStream(stream, writer)
        before = []
        after = []
        kept = []
        for operands, op in content.operations:
            if op == b"Tr":
                assert (
                    float(operands[0]) < 4
                ), "Text clipping requires explicit treatment"
            if op not in TEXT_OPS:
                before.append((repr(operands), op.decode("latin1")))
                kept.append((operands, op))
                after.append((repr(operands), op.decode("latin1")))
        assert before == after
        fingerprints.append(
            {
                "stream": label,
                "nontext_operations": len(before),
                "sha256": hashlib.sha256(repr(before).encode()).hexdigest(),
            }
        )
        content.operations = kept
        stream.set_data(content.get_data())
        if "/Font" in resources:
            del resources["/Font"]
        for name, value in resources.get("/XObject", {}).items():
            obj = value.get_object()
            if obj.get("/Subtype") == "/Form":
                clean(
                    obj,
                    obj.get("/Resources", resources).get_object(),
                    label + "/" + str(name),
                )

    pdfpage = writer.pages[0]
    # Flatten the page content array into one stream without changing drawing operations.
    contents = ContentStream(pdfpage.get_contents(), writer)
    pdfpage[NameObject("/Contents")] = writer._add_object(contents)
    clean(contents, pdfpage["/Resources"].get_object(), "page")
    overlay = fitz.open()
    newpage = overlay.new_page(width=page.rect.width, height=page.rect.height)
    faces = {
        "regular": "STIXTwoText-Regular.ttf",
        "italic": "STIXTwoText-Italic.ttf",
        "medium": "STIXTwoText-Medium.ttf",
        "math": "STIXTwoMath-Regular.otf",
    }
    fonts = {k: fitz.Font(fontfile=str(font_dir / v)) for k, v in faces.items()}
    for k, v in faces.items():
        newpage.insert_font(fontname="DAMADI" + k, fontfile=str(font_dir / v))
    for span in overlays:
        text = span["text"]
        face = span["face"]
        font = fonts[face]
        assert all(c.isspace() or font.has_glyph(ord(c)) for c in text), (
            text,
            face,
        )
        width = font.text_length(text, fontsize=span["size"])
        source_width = span["bbox"][2] - span["bbox"][0]
        factor = source_width / width if width > 0 else 1
        origin = fitz.Point(span["origin"])
        rgb = fitz.sRGB_to_pdf(span["color"])
        newpage.insert_text(
            origin,
            text,
            fontsize=span["size"],
            fontname="DAMADI" + face,
            color=rgb,
            morph=(origin, fitz.Matrix(factor, 1)),
        )
        span["horizontal_scale"] = factor
    overlay_writer = PdfWriter()
    overlay_writer.clone_document_from_reader(
        PdfReader(io.BytesIO(overlay.tobytes()))
    )
    # MuPDF's CFF insertion needs an explicit UTF-16 mapping for non-BMP math letters.
    math_chars = sorted(
        {c for span in overlays if span["face"] == "math" for c in span["text"]}
    )
    if math_chars:
        pairs = [
            f'<{fonts["math"].has_glyph(ord(c)):04X}> <{c.encode("utf-16-be").hex().upper()}>'
            for c in math_chars
        ]
        cmap = (
            "/CIDInit /ProcSet findresource begin\n12 dict begin\nbegincmap\n/CIDSystemInfo << /Registry (Adobe) /Ordering (UCS) /Supplement 0 >> def\n/CMapName /DAMADI-Math-Unicode def\n/CMapType 2 def\n1 begincodespacerange\n<0000> <FFFF>\nendcodespacerange\n"
            + str(len(pairs))
            + " beginbfchar\n"
            + "\n".join(pairs)
            + "\nendbfchar\nendcmap\nCMapName currentdict /CMap defineresource pop\nend\nend\n"
        )
        stream = DecodedStreamObject()
        stream.set_data(cmap.encode())
        overlay_writer.pages[0]["/Resources"]["/Font"]["/DAMADImath"].get_object()[
            NameObject("/ToUnicode")
        ] = overlay_writer._add_object(stream)
    pdfpage.merge_page(overlay_writer.pages[0])
    writer.add_metadata(
        {
            "/Title": "DAMADI STIX presentation of retained network figure",
            "/Author": "Simon Allzén",
            "/Subject": "Original graph geometry and source label positions retained; STIX typography only",
        }
    )
    destination.parent.mkdir(parents=True, exist_ok=True)
    with destination.open("wb") as output:
        writer.write(output)
    rendered = fitz.open(destination)
    expected = "".join(x["text"] for x in overlays)
    # Verify every source label's characters remain in the presentation text layer.
    actual = rendered[0].get_text()
    from collections import Counter

    assert Counter(c for c in expected if not c.isspace()) == Counter(
        c for c in actual if not c.isspace()
    )
    return {
        "source": source.name,
        "pdf": destination.name,
        "label_spans": overlays,
        "nontext_stream_fingerprints": fingerprints,
        "all_nontext_operations_unchanged": True,
        "source_label_characters_retained": True,
        "font_sizes_and_origins_retained": True,
        "width_mm": page.rect.width * 25.4 / 72,
        "height_mm": page.rect.height * 25.4 / 72,
    }