pd: wire cover traffic through the sweep as a fourth table

engine gains a cover_rates axis: each rate plays a timeline through the same
graph and pairs it with the epoch emission budget, which needs no graph and so is
computed alongside rather than inside the window. Seeds are separate streams
(traffic_seedseq for the timeline and clocks, stake_seedseq for the stake draw
and budget), so the stake distribution is independent of the topology and of the
message schedule.

sweep writes traffic.parquet only when a cover-traffic study actually ran, so
every existing config keeps producing exactly three tables. quota_summary reports
the measured ceiling beside the predicted one in the same row, so a run can be
checked against the closed form instead of asked to be believed.

Two figures: blending against cover rate and release delay with the
rate*(2M+1)/3 law overlaid, and the quota ceiling with the measured transition
band against the prediction.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

tools/simulators/blend/pd/configs/cover-traffic.yaml

@@ -0,0 +1,36 @@
+# Cover traffic: the anonymity set it buys, and the stake concentration it tolerates.
+#
+# Every node emits on its own uniformly-chosen slots at rate cover_rate_mult/N, so the default of
+# 1.0 puts one emission per second on the whole network. Winning the block lottery consumes the
+# next scheduled cover, which is what keeps every node's emission COUNT identical whether or not it
+# produces blocks -- the property cover traffic exists to buy.
+#
+# Two things are measured. On a timeline: blending, the broadcasts a relay saw between consecutive
+# releases (the anonymity set, since an observer cannot tell which one it forwarded), against
+# mixing, the messages it held at once. At the baseline rate mixing is essentially nil, so the rate
+# is swept over three decades to find where per-relay mixing actually begins. Over an epoch: which
+# nodes win more proposals than their emission quota and so cannot stay uniform.
+#
+# max_blend_delay is swept alongside the rate because delay is the cheaper lever -- blending is
+# rate*(2M+1)/3, so a longer release interval multiplies the anonymity set without adding a single
+# message to the wire. The cost of that delay is priced in section 3.1.
+n_nodes: [20000]
+degree: [8]
+blend_hops: [3]
+max_blend_delay: [3, 10, 30]
+unresponsive_frac: [0.0]
+cover_rate_mult: [1.0, 4.0, 16.0, 64.0, 256.0]
+f_adv: [0.2]
+adversary_mode: [random]
+seeds: 4
+base:
+  # heavy-tailed stake is what makes the quota ceiling bite: the head sits orders of magnitude
+  # above it while the tail sits far below, so the breakpoint is measurable rather than assumed.
+  stake_dist: zipf
+  stake_zipf_a: 1.0
+  stake_inference_ratio: 1.0    # D_hat/D from the consensus study; 1.0 = an accurate estimator
+  block_interval_slots: 30
+  slots_per_epoch: 648000
+  traffic_window_slots: 900
+  n_rounds: 20                  # propagation is not the subject here; keep it cheap
+  n_placements: 1

tools/simulators/blend/pd/src/pd/config.py

@@ -164,6 +164,7 @@ class SweepConfig:
     unresponsive_frac: list[float] = field(default_factory=lambda: [0.0])
     churn_mode: list[str] = field(default_factory=lambda: ["uniform"])
     redundancy: list[int] = field(default_factory=lambda: [1])
+    cover_rate_mult: list[float] = field(default_factory=list)   # empty = no cover-traffic study
     f_adv: list[float] = field(default_factory=lambda: [0.1, 0.2, 0.33, 0.5])
     adversary_mode: list[str] = field(default_factory=lambda: ["random"])
     seeds: int = 8                       # number of graph_seed values (topology ensemble)
@@ -201,7 +202,8 @@ class SweepConfig:
         d = dict(d)
         base = d.pop("base", {})
         known = {"n_nodes", "degree", "blend_hops", "max_blend_delay", "unresponsive_frac",
-                 "churn_mode", "redundancy", "f_adv", "adversary_mode", "seeds"}
+                 "churn_mode", "redundancy", "cover_rate_mult", "f_adv", "adversary_mode",
+                 "seeds"}
         unknown = set(d) - known
         if unknown:
             raise ValueError(f"unknown sweep keys: {sorted(unknown)}")

tools/simulators/blend/pd/src/pd/engine.py

@@ -7,29 +7,42 @@ propagation and adversary sub-grids.
 
 from __future__ import annotations
 
+import dataclasses
+
 import numpy as np
 
 from .adversary import adversary_metrics, deanon_metrics, place_adversary
 from .config import WORSTCASE_MODES, SimConfig
 from .graph import build_graph
-from .metrics import adversary_row, deanon_row, propagation_row
+from .metrics import adversary_row, deanon_row, propagation_row, traffic_row
 from .propagation import assign_responsive, propagation_metrics
-from .rng import placement_seedseq, responsive_seedseq, round_seedseq
+from .quota import assign_stake, quota_summary
+from .rng import (
+    placement_seedseq,
+    responsive_seedseq,
+    round_seedseq,
+    stake_seedseq,
+    traffic_seedseq,
+)
+from .traffic import simulate_window, traffic_metrics
 
 
 def run_graph_cell(base: SimConfig, prop_grid: list[tuple[int, int]],
                    unresponsive_fracs: list[float], redundancies: list[int],
                    adv_grid: list[tuple[float, str]],
                    churn_modes: list[str] | None = None,
-                   ) -> tuple[list[dict], list[dict], list[dict]]:
-    """Build ``base``'s topology once; return (propagation, adversary, deanonymization rows).
+                   cover_rates: list[float] | None = None,
+                   ) -> tuple[list[dict], list[dict], list[dict], list[dict]]:
+    """Build ``base``'s topology once; return (propagation, adversary, deanon, traffic rows).
 
     ``base`` carries the topology (n_nodes, degree, graph_seed) and all shared knobs;
     ``prop_grid`` = [(blend_hops, max_blend_delay)], ``unresponsive_fracs`` = the relay-dropout
     axis, ``redundancies`` = the messaging-redundancy axis (R independent cascades per emission),
     ``adv_grid`` = [(f_adv, mode)]. Deanonymization crosses each adversary placement with the
     propagation grid's blend-path lengths and redundancies, so it is emitted alongside the
-    adversary rows.
+    adversary rows. ``cover_rates`` (empty by default) turns on the cover-traffic study: each rate
+    plays a timeline through the same graph and pairs it with the epoch emission budget, which is
+    graph-free and therefore computed separately.
     """
     graph = build_graph(base)
     blend_hops_set = sorted({bh for bh, _ in prop_grid})
@@ -66,7 +79,24 @@ def run_graph_cell(base: SimConfig, prop_grid: list[tuple[int, int]],
                     dz = deanon_metrics(graph.n, adv["n_adv"], adv["observed_frac"], bh, R)
                     deanon_rows.append(deanon_row(base, bh, f_adv, mode, rep, R, adv, dz))
 
-    return prop_rows, adv_rows, deanon_rows
+    traffic_rows: list[dict] = []
+    for rate in (cover_rates or []):
+        cfg = dataclasses.replace(base, cover_rate_mult=rate)
+        f = 1.0 / base.block_interval_slots
+        srng = np.random.default_rng(stake_seedseq(base, rate))
+        stake = assign_stake(base.n_nodes, base.stake_dist, srng, base.stake_zipf_a)
+        quota = quota_summary(stake, f, base.n_nodes, base.slots_per_epoch, srng,
+                              base.stake_inference_ratio, rate)
+        for blend_hops, max_blend_delay in prop_grid:
+            trng = np.random.default_rng(
+                traffic_seedseq(base, blend_hops, max_blend_delay, rate))
+            win = simulate_window(graph, cfg, trng, base.traffic_window_slots,
+                                  max_blend_delay, blend_hops)
+            tm = traffic_metrics(win, cfg, max_blend_delay)
+            traffic_rows.append(
+                traffic_row(base, blend_hops, max_blend_delay, rate, tm, quota))
+
+    return prop_rows, adv_rows, deanon_rows, traffic_rows
 
 
 def run_trajectory(config: SimConfig) -> dict:

tools/simulators/blend/pd/src/pd/plotting/figures.py

@@ -638,3 +638,74 @@ def redundancy_tradeoff(prop, adv, deanon):
     ax.set_title(f"Redundancy: reliability vs anonymity (N={n:,}, degree={deg}, blend_hops={bh})")
     ax.legend()
     return fig
+
+
+# --- cover traffic (traffic table only) -----------------------------------------------------------
+
+def blending_vs_rate_and_delay(traffic: pd.DataFrame):
+    """The anonymity set against the two knobs that buy it: cover rate and release delay.
+
+    Blending is the number of broadcasts a relay saw between consecutive releases -- an observer
+    cannot tell which of them it forwarded. Intervals sampled at a release are size-biased, so the
+    law is ``rate * (2M+1)/3``, exactly twice the mean hold; the dashed lines are that prediction.
+    Delay is the cheaper lever: it multiplies the set without adding a single message to the wire.
+    """
+    if traffic is None or not len(traffic) or "blending_mean" not in traffic:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    d = traffic[traffic.n_nodes == traffic.n_nodes.max()]
+    if d["cover_rate_mult"].nunique() < 2 and d["max_blend_delay"].nunique() < 2:
+        return None
+    fig, ax = plt.subplots()
+    for i, M in enumerate(sorted(d.max_blend_delay.unique())):
+        s = d[d.max_blend_delay == M].groupby("cover_rate_mult").blending_mean.mean().reset_index()
+        c = style.color_for(i)
+        ax.plot(s.cover_rate_mult, s.blending_mean, "-o", ms=4, color=c,
+                label=f"max_blend_delay={int(M)}s")
+        ax.plot(s.cover_rate_mult, s.cover_rate_mult * (2 * M + 1) / 3, "--", lw=0.8,
+                color=c, alpha=0.6)
+    ax.set_xscale("log")
+    ax.set_yscale("log")
+    ax.set_xlabel("cover-traffic rate  (messages/second, network-wide)")
+    ax.set_ylabel("blending: broadcasts seen between releases")
+    ax.set_title("Anonymity set vs cover rate and release delay\n"
+                 "dashed = $rate\\cdot(2M{+}1)/3$")
+    ax.legend()
+    return fig
+
+
+def quota_stake_ceiling(traffic: pd.DataFrame):
+    """Where the uniform-emission guarantee breaks: measured ceiling against the closed form.
+
+    Cover traffic keeps every node's emission count identical only while its block proposals fit
+    inside the quota. Nodes above the ceiling must emit more often than everyone else, which is the
+    very signal the scheme exists to hide. Raising the cover rate raises the ceiling in proportion.
+    """
+    if traffic is None or not len(traffic) or "s_max_predicted" not in traffic:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    d = traffic[traffic.n_nodes == traffic.n_nodes.max()]
+    g = d.groupby("cover_rate_mult").agg(
+        pred=("s_max_predicted", "mean"), safe=("alpha_max_99", "mean"),
+        meas_hi=("max_compliant_stake", "mean"), meas_lo=("min_overrun_stake", "mean"),
+        top=("top_stake", "mean")).reset_index()
+    if not len(g):
+        return None
+    fig, ax = plt.subplots()
+    ax.plot(g.cover_rate_mult, g.pred * 100, "-", lw=1.4, color=style.color_for(0),
+            label="predicted ceiling $s_{max}$")
+    ax.plot(g.cover_rate_mult, g.safe * 100, ":", lw=1.0, color=style.color_for(0),
+            label="99%-safe ceiling")
+    ax.fill_between(g.cover_rate_mult, g.meas_lo * 100, g.meas_hi * 100, alpha=0.25,
+                    color=style.color_for(1), label="measured transition band")
+    ax.plot(g.cover_rate_mult, g.top * 100, "--", lw=1.0, color=style.color_for(2),
+            label="largest staker present")
+    ax.set_xscale("log")
+    ax.set_yscale("log")
+    ax.set_xlabel("cover-traffic rate  (messages/second, network-wide)")
+    ax.set_ylabel("stake (%)")
+    ax.set_title("Emission quota: the most stake a node can hold and stay uniform")
+    ax.legend(fontsize=7)
+    return fig

tools/simulators/blend/pd/src/pd/plotting/make_figures.py

@@ -40,13 +40,20 @@ _DEANON_BUILDERS = [
     ("21_redundancy_time_to_link", figures.redundancy_time_to_link),
 ]
 
+# (traffic,) builders — the cover-traffic study, rendered only when that table exists.
+_TRAFFIC_BUILDERS = [
+    ("23_blending_vs_rate_and_delay", figures.blending_vs_rate_and_delay),
+    ("24_quota_stake_ceiling", figures.quota_stake_ceiling),
+]
+
 
 def render(prop_df: pd.DataFrame, adv_df: pd.DataFrame, deanon_df: pd.DataFrame,
-           out_dir: Path) -> list[Path]:
+           out_dir: Path, traffic_df: pd.DataFrame | None = None) -> list[Path]:
     out_dir = Path(out_dir)
     prov = figures._prov(prop_df, adv_df, deanon_df)
     jobs = ([(name, fn, (prop_df, adv_df)) for name, fn in _BUILDERS]
-            + [(name, fn, (prop_df, adv_df, deanon_df)) for name, fn in _DEANON_BUILDERS])
+            + [(name, fn, (prop_df, adv_df, deanon_df)) for name, fn in _DEANON_BUILDERS]
+            + [(name, fn, (traffic_df,)) for name, fn in _TRAFFIC_BUILDERS])
     written: list[Path] = []
     for name, fn, fn_args in jobs:
         try:

tools/simulators/blend/pd/src/pd/quota.py

@@ -61,6 +61,30 @@ def assign_stake(n_nodes: int, dist: str, rng: np.random.Generator,
     return s
 
 
+def quota_summary(stake: np.ndarray, f: float, n_nodes: int, slots_per_epoch: int,
+                  rng: np.random.Generator, stake_inference_ratio: float = 1.0,
+                  cover_rate_mult: float = 1.0) -> dict:
+    """Scalar view of :func:`simulate_epoch_emissions`, for a result row.
+
+    Reports the measured ceiling (where compliance actually breaks) beside the predicted one, so a
+    run can be checked against the closed form rather than asked to be believed.
+    """
+    r = simulate_epoch_emissions(stake, f, n_nodes, slots_per_epoch, rng,
+                                 stake_inference_ratio, cover_rate_mult)
+    return {
+        "quota_per_epoch": float(r["quota"]),
+        "compliant_frac": r["compliant_frac"],
+        "max_compliant_stake": r["max_compliant_stake"],
+        "min_overrun_stake": r["min_overrun_stake"],
+        "total_overrun": int(r["overrun"].sum()),
+        "top_stake": float(np.max(r["stake"])),
+        "alpha_max_predicted": alpha_max(n_nodes, f, cover_rate_mult),
+        "s_max_predicted": s_max_true(n_nodes, f, stake_inference_ratio, cover_rate_mult),
+        "alpha_max_99": max_alpha_for_confidence(f, n_nodes, slots_per_epoch, 0.99,
+                                                 cover_rate_mult),
+    }
+
+
 def inferred_alpha(stake: np.ndarray, stake_inference_ratio: float = 1.0) -> np.ndarray:
     """Convert true relative stake to the **inferred** relative stake the lottery actually uses.
 

tools/simulators/blend/pd/src/pd/rng.py

@@ -61,6 +61,25 @@ def round_seedseq(config: SimConfig, blend_hops: int, max_blend_delay: int,
     ))
 
 
+def traffic_seedseq(config: SimConfig, blend_hops: int, max_blend_delay: int,
+                    cover_rate_mult: float) -> np.random.SeedSequence:
+    """Cover-traffic timeline seed: emissions, relay paths, and every node's release clock."""
+    return np.random.SeedSequence(_digest(
+        config.root_seed, "traffic", config.n_nodes, config.degree, config.graph_seed,
+        blend_hops, max_blend_delay, cover_rate_mult, config.traffic_window_slots,
+        config.block_interval_slots, config.transport_jitter_mean_ms,
+    ))
+
+
+def stake_seedseq(config: SimConfig, cover_rate_mult: float) -> np.random.SeedSequence:
+    """Stake draw + epoch emission budget; independent of the timeline and of the graph."""
+    return np.random.SeedSequence(_digest(
+        config.root_seed, "stake", config.n_nodes, config.graph_seed, config.stake_dist,
+        config.stake_zipf_a, config.stake_inference_ratio, config.slots_per_epoch,
+        config.block_interval_slots, cover_rate_mult,
+    ))
+
+
 def placement_seedseq(config: SimConfig, f_adv: float, mode: str,
                       placement_rep: int) -> np.random.SeedSequence:
     """Adversary-placement seed, independent of the graph draw and the rounds."""

tools/simulators/blend/pd/src/pd/sweep.py

@@ -38,36 +38,40 @@ def new_run_dir(outdir: Path, label: str) -> Path:
 
 
 def _cell_worker(base: SimConfig, prop_grid, unresponsive_fracs, redundancies, adv_grid,
-                 churn_modes):
+                 churn_modes, cover_rates):
     return run_graph_cell(base, prop_grid, unresponsive_fracs, redundancies, adv_grid,
-                          churn_modes)
+                          churn_modes, cover_rates)
 
 
-def run_sweep(sweep: SweepConfig,
-              n_jobs: int = -1) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]:
+def run_sweep(sweep: SweepConfig, n_jobs: int = -1) -> tuple[pd.DataFrame, ...]:
     cells = sweep.graph_cells()
     prop_grid = sweep.prop_grid()
     unresponsive_fracs = list(sweep.unresponsive_frac)
     redundancies = list(sweep.redundancy)
     churn_modes = list(sweep.churn_mode)
+    cover_rates = list(sweep.cover_rate_mult)
     adv_grid = sweep.adv_grid()
     bases = [sweep.base_config(n, d, g) for (n, d, g) in cells]
     results = Parallel(n_jobs=n_jobs, prefer="processes")(
         delayed(_cell_worker)(base, prop_grid, unresponsive_fracs, redundancies, adv_grid,
-                              churn_modes)
+                              churn_modes, cover_rates)
         for base in tqdm(bases, desc="topologies")
     )
-    prop_rows = [r for pr, _, _ in results for r in pr]
-    adv_rows = [r for _, ar, _ in results for r in ar]
-    deanon_rows = [r for _, _, dr in results for r in dr]
-    return pd.DataFrame(prop_rows), pd.DataFrame(adv_rows), pd.DataFrame(deanon_rows)
+    prop_rows = [r for pr, _, _, _ in results for r in pr]
+    adv_rows = [r for _, ar, _, _ in results for r in ar]
+    deanon_rows = [r for _, _, dr, _ in results for r in dr]
+    traffic_rows = [r for _, _, _, tr in results for r in tr]
+    return (pd.DataFrame(prop_rows), pd.DataFrame(adv_rows), pd.DataFrame(deanon_rows),
+            pd.DataFrame(traffic_rows))
 
 
 def persist(prop_df: pd.DataFrame, adv_df: pd.DataFrame, deanon_df: pd.DataFrame,
-            run_dir: Path) -> None:
+            traffic_df: pd.DataFrame, run_dir: Path) -> None:
     prop_df.to_parquet(run_dir / "propagation.parquet", index=False)
     adv_df.to_parquet(run_dir / "adversary.parquet", index=False)
     deanon_df.to_parquet(run_dir / "deanon.parquet", index=False)
+    if len(traffic_df):                      # only written when a cover-traffic study ran
+        traffic_df.to_parquet(run_dir / "traffic.parquet", index=False)
 
 
 def main(argv: list[str] | None = None) -> int:
@@ -82,14 +86,15 @@ def main(argv: list[str] | None = None) -> int:
     sweep = load_sweep_yaml(args.config)
     label = args.label or Path(args.config).stem
     run_dir = new_run_dir(Path(args.outdir), label)
-    prop_df, adv_df, deanon_df = run_sweep(sweep, n_jobs=args.n_jobs)
-    persist(prop_df, adv_df, deanon_df, run_dir)
+    prop_df, adv_df, deanon_df, traffic_df = run_sweep(sweep, n_jobs=args.n_jobs)
+    persist(prop_df, adv_df, deanon_df, traffic_df, run_dir)
+    extra = f" + {len(traffic_df)} traffic" if len(traffic_df) else ""
     print(f"wrote {len(prop_df)} propagation + {len(adv_df)} adversary + "
-          f"{len(deanon_df)} deanon rows -> {run_dir}")
+          f"{len(deanon_df)} deanon{extra} rows -> {run_dir}")
 
     if not args.no_figures:
         from .plotting.make_figures import render
-        figs = render(prop_df, adv_df, deanon_df, run_dir / "figures")
+        figs = render(prop_df, adv_df, deanon_df, run_dir / "figures", traffic_df)
         print(f"wrote {len(figs)} figures -> {run_dir / 'figures'}")
     return 0
 

tools/simulators/blend/pd/tests/test_deanon.py

@@ -85,7 +85,7 @@ def test_engine_emits_deanon_rows():
     base = SimConfig(n_nodes=1000, degree=8, graph_seed=0, n_placements=2)
     prop_grid = [(2, 0), (3, 0)]                 # distinct blend_hops = {2, 3}
     adv_grid = [(0.2, "random"), (0.0, "random")]
-    prop_rows, adv_rows, deanon_rows = run_graph_cell(base, prop_grid, [0.0], [1], adv_grid)
+    prop_rows, adv_rows, deanon_rows, _ = run_graph_cell(base, prop_grid, [0.0], [1], adv_grid)
 
     # one deanon row per (placement, distinct blend_hops, redundancy)
     assert len(deanon_rows) == len(adv_rows) * 2