pd: stake distribution, so the emission ceiling is measured not asserted

The quota ceiling was closed-form only. This adds per-node stake so a run can
show nodes actually breaking it.

- assign_stake: uniform, or heavy-tailed zipf (s ~ 1/rank^a), which is what makes
  the ceiling bite -- the head sits orders of magnitude above it, the tail far below;
- inferred_alpha: converts true relative stake to the sigma/D_hat the lottery
  actually weighs, so a low estimate inflates every node alpha;
- simulate_epoch_emissions: measures the budget over a full epoch. Overrun happens
  at epoch scale and needs no graph, so this is cheap: proposals are Binomial over
  the epoch slots, a proposal cancels the next cover, and a node stays at exactly
  its quota until its wins no longer fit -- at which point it emits more often than
  everyone else, which is the signal cover traffic exists to suppress.

Measured against the closed form at N=20,000, zipf stake, over an epoch: the
predicted ceiling falls inside the transition band every time, and at D_hat/D = 1
the smallest overrunning node sits at 0.1468% against a predicted 0.1475%. The
D_hat/D normalisation is confirmed empirically -- deflating the estimate to 0.64
pulls the measured ceiling down with it, as the (D_hat/D)*alpha_max form requires.
With heavy-tailed stake 99.7% of nodes comply and only the head breaks; the
largest holder at 9.5% stake is some 65x over its allowance.

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

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

@@ -12,6 +12,8 @@ AdversaryMode = Literal[
     "random", "worstcase_coverage", "worstcase_eclipse", "worstcase_degree"
 ]
 ChurnMode = Literal["uniform", "regional"]
+StakeDist = Literal["uniform", "zipf"]
+_STAKE_DISTS = ("uniform", "zipf")
 _DISTS = ("geo", "fixed", "uniform", "exp")
 _MODES = ("random", "worstcase_coverage", "worstcase_eclipse", "worstcase_degree")
 WORSTCASE_MODES = ("worstcase_coverage", "worstcase_eclipse", "worstcase_degree")
@@ -41,6 +43,8 @@ class SimConfig:
     slots_per_epoch: int = 648_000      # epoch length, for the per-node emission quota
     stake_inference_ratio: float = 1.0  # D_hat/D from the consensus study; scales the stake ceiling
     traffic_window_slots: int = 600     # simulated timeline length (seconds) per cell
+    stake_dist: StakeDist = "uniform"   # per-node stake: "uniform" or heavy-tailed "zipf"
+    stake_zipf_a: float = 1.0           # zipf exponent; larger = more concentrated
     n_rounds: int = 200                 # random-sender rounds per topology
     transport_jitter_mean_ms: float = 5.0
     processing_lags_ms: tuple[float, ...] = (10.0, 50.0, 100.0)
@@ -82,6 +86,10 @@ class SimConfig:
         if self.stake_inference_ratio <= 0.0:
             raise ValueError(
                 f"stake_inference_ratio (D_hat/D) must be > 0, got {self.stake_inference_ratio}")
+        if self.stake_dist not in _STAKE_DISTS:
+            raise ValueError(f"stake_dist must be one of {_STAKE_DISTS}")
+        if self.stake_zipf_a <= 0.0:
+            raise ValueError(f"stake_zipf_a must be > 0, got {self.stake_zipf_a}")
         if self.n_regions < 1:
             raise ValueError(f"n_regions must be >= 1, got {self.n_regions}")
         if self.n_regions > 1:
@@ -132,7 +140,8 @@ class SimConfig:
             self.blend_hops, self.max_blend_delay,
             self.unresponsive_frac, self.churn_mode, self.redundancy,
             self.cover_rate_mult, self.block_interval_slots, self.slots_per_epoch,
-            self.stake_inference_ratio, self.traffic_window_slots, self.n_rounds,
+            self.stake_inference_ratio, self.traffic_window_slots,
+            self.stake_dist, self.stake_zipf_a, self.n_rounds,
             self.transport_jitter_mean_ms, self.processing_lags_ms, self.processing_lag_probs,
             self.link_latency_dist, self.link_latency_mean_ms, self.coverage_pcts,
             self.f_adv, self.adversary_mode, self.n_placements, self.worstcase_max_n,

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

@@ -1,4 +1,4 @@
-"""Flat parquet-row builders for the two result tables."""
+"""Flat parquet-row builders for the result tables."""
 
 from __future__ import annotations
 
@@ -27,6 +27,30 @@ def propagation_row(config: SimConfig, blend_hops: int, max_blend_delay: int,
     }
 
 
+def traffic_row(config: SimConfig, blend_hops: int, max_blend_delay: int,
+                cover_rate_mult: float, traffic: dict, quota: dict) -> dict:
+    """One cover-traffic cell: what the timeline measured, plus the epoch emission budget.
+
+    ``traffic`` comes from the windowed simulation (blending, mixing, counts) and ``quota`` from
+    the epoch-scale emission budget, which needs no graph and so is computed separately.
+    """
+    return {
+        "n_nodes": config.n_nodes,
+        "degree": config.degree,
+        "blend_hops": blend_hops,
+        "max_blend_delay": max_blend_delay,
+        "cover_rate_mult": cover_rate_mult,
+        "block_interval_slots": config.block_interval_slots,
+        "slots_per_epoch": config.slots_per_epoch,
+        "stake_dist": config.stake_dist,
+        "stake_inference_ratio": config.stake_inference_ratio,
+        "graph_seed": config.graph_seed,
+        "traffic_window_slots": config.traffic_window_slots,
+        **traffic,
+        **quota,
+    }
+
+
 def adversary_row(config: SimConfig, f_adv: float, mode: str, placement_rep: int,
                   adv: dict) -> dict:
     return {

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

@@ -33,6 +33,76 @@ from __future__ import annotations
 
 import math
 
+import numpy as np
+
+
+def assign_stake(n_nodes: int, dist: str, rng: np.random.Generator,
+                 zipf_a: float = 1.0) -> np.ndarray:
+    """Per-node **true** relative stake ``s = sigma/D``, summing to 1.
+
+    ``uniform`` gives every node ``1/n_nodes`` -- the case where the quota binds on nobody until
+    the network is small. ``zipf`` makes stake heavy-tailed (``s ~ 1/rank**zipf_a``), which is what
+    real stake looks like and what makes the ceiling bite: the head of the distribution sits orders
+    of magnitude above it while the tail sits far below.
+    """
+    if n_nodes < 1:
+        raise ValueError("n_nodes must be >= 1")
+    if dist == "uniform":
+        s = np.full(n_nodes, 1.0 / n_nodes)
+    elif dist == "zipf":
+        if zipf_a <= 0:
+            raise ValueError("zipf_a must be > 0")
+        ranks = np.arange(1, n_nodes + 1, dtype=float)
+        s = ranks ** (-zipf_a)
+        s = s / s.sum()
+        rng.shuffle(s)                       # stake is not correlated with node id
+    else:
+        raise ValueError(f"unknown stake distribution {dist!r}")
+    return s
+
+
+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.
+
+    The threshold is derived from ``D_hat``, so a node's lottery weight is ``sigma/D_hat``, i.e.
+    ``s * D/D_hat = s / stake_inference_ratio``. An estimator that runs low makes every node win
+    more often, which is why it tightens the true-stake ceiling.
+    """
+    if stake_inference_ratio <= 0.0:
+        raise ValueError("stake_inference_ratio must be > 0")
+    return np.asarray(stake, dtype=float) / stake_inference_ratio
+
+
+def simulate_epoch_emissions(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:
+    """Measure the quota over one epoch: who wins more proposals than their emission budget.
+
+    Needs no graph -- only the counts matter. Each node's proposals are Binomial over the epoch's
+    slots at its lottery probability; a proposal cancels the next scheduled cover, so a node stays
+    at exactly its quota while its wins fit inside it and **overruns** once they do not. An
+    overrunning node emits more often than everybody else, which is precisely the signal cover
+    traffic exists to suppress.
+    """
+    alpha = inferred_alpha(stake, stake_inference_ratio)
+    phi = 1.0 - (1.0 - f) ** alpha
+    blocks = rng.binomial(slots_per_epoch, phi)
+    quota = quota_per_epoch(n_nodes, slots_per_epoch, cover_rate_mult)
+    overrun = np.maximum(0, blocks - math.floor(quota))
+    compliant = overrun == 0
+    return {
+        "stake": np.asarray(stake, dtype=float),
+        "alpha": alpha,
+        "blocks": blocks,
+        "quota": quota,
+        "overrun": overrun,
+        "compliant": compliant,
+        "emissions": np.where(compliant, math.floor(quota), blocks),
+        "compliant_frac": float(compliant.mean()),
+        "max_compliant_stake": float(stake[compliant].max()) if compliant.any() else 0.0,
+        "min_overrun_stake": float(stake[~compliant].min()) if (~compliant).any() else float("nan"),
+    }
+
 
 def emission_quota_per_slot(n_nodes: int, cover_rate_mult: float = 1.0) -> float:
     """Emissions allowed per node per slot. The default rate puts one emission per slot on the

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

@@ -18,6 +18,7 @@ def test_key_covers_every_field():
                "unresponsive_frac": 0.2, "churn_mode": "regional", "redundancy": 2,
                "cover_rate_mult": 2.0, "block_interval_slots": 60, "slots_per_epoch": 1000,
                "stake_inference_ratio": 0.7, "traffic_window_slots": 100,
+               "stake_dist": "zipf", "stake_zipf_a": 1.5,
                "n_rounds": 10, "transport_jitter_mean_ms": 1.0,
                "processing_lags_ms": (11.0, 51.0, 101.0), "processing_lag_probs": (0.6, 0.3, 0.1),
                "link_latency_dist": "fixed", "link_latency_mean_ms": 1.0,

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

@@ -6,12 +6,15 @@ import numpy as np
 
 from pd.quota import (
     alpha_max,
+    assign_stake,
     emission_quota_per_slot,
     expected_blocks_per_epoch,
+    inferred_alpha,
     max_alpha_for_confidence,
     quota_exceedance_prob,
     quota_per_epoch,
     s_max_true,
+    simulate_epoch_emissions,
     win_prob,
 )
 
@@ -83,3 +86,74 @@ def test_exceedance_matches_a_direct_simulation():
     wins = rng.binomial(S, win_prob(a, F), size=20_000)
     emp = float(np.mean(wins > math.floor(quota)))
     assert abs(closed - emp) < max(0.01, 0.1 * closed)
+
+
+# --- stake distribution and the measured ceiling --------------------------------------------------
+
+def test_stake_distributions_normalise_and_zipf_is_heavy_tailed():
+    n = 5_000
+    rng = np.random.default_rng(0)
+    uni = assign_stake(n, "uniform", rng)
+    zipf = assign_stake(n, "zipf", rng, zipf_a=1.0)
+    for s in (uni, zipf):
+        assert abs(s.sum() - 1.0) < 1e-12
+        assert (s > 0).all()
+    assert np.allclose(uni, 1.0 / n)
+    assert zipf.max() > 50 * uni.max()          # a real head, unlike the flat case
+
+
+def test_inferred_alpha_divides_by_the_estimator_ratio():
+    """The lottery weighs sigma/D_hat, so a low estimate inflates every node's alpha."""
+    s = np.array([0.001, 0.01])
+    assert np.allclose(inferred_alpha(s, 1.0), s)
+    assert np.allclose(inferred_alpha(s, 0.5), s * 2.0)
+
+
+def test_uniform_stake_stays_inside_the_quota_at_scale():
+    """At 1/N each, every node's block rate is f/N -- far under a 1/N emission budget."""
+    n, S = 20_000, 648_000
+    s = assign_stake(n, "uniform", np.random.default_rng(1))
+    r = simulate_epoch_emissions(s, F, n, S, np.random.default_rng(2))
+    assert r["compliant_frac"] == 1.0
+    assert r["overrun"].sum() == 0
+
+
+def test_heavy_tailed_stake_makes_the_head_overrun_its_quota():
+    n, S = 20_000, 648_000
+    s = assign_stake(n, "zipf", np.random.default_rng(3), zipf_a=1.0)
+    r = simulate_epoch_emissions(s, F, n, S, np.random.default_rng(4))
+    assert 0.0 < r["compliant_frac"] < 1.0            # the head breaks, the tail does not
+    assert r["min_overrun_stake"] > r["stake"].min()  # it is the large holders that break
+    assert r["overrun"][np.argmax(s)] > 0             # the biggest staker certainly does
+
+
+def test_the_measured_ceiling_matches_the_closed_form():
+    """Where compliance actually breaks must bracket the analytic alpha_max."""
+    n, S = 20_000, 648_000
+    s = assign_stake(n, "zipf", np.random.default_rng(5), zipf_a=0.8)
+    r = simulate_epoch_emissions(s, F, n, S, np.random.default_rng(6))
+    predicted = alpha_max(n, F)
+    assert r["max_compliant_stake"] < 3.0 * predicted
+    assert r["min_overrun_stake"] > 0.3 * predicted
+
+
+def test_a_low_stake_estimate_tightens_the_measured_ceiling():
+    """D_hat/D is an input, and lowering it must push more nodes over their quota."""
+    n, S = 20_000, 648_000
+    s = assign_stake(n, "zipf", np.random.default_rng(7), zipf_a=1.0)
+    accurate = simulate_epoch_emissions(s, F, n, S, np.random.default_rng(8),
+                                        stake_inference_ratio=1.0)
+    deflated = simulate_epoch_emissions(s, F, n, S, np.random.default_rng(8),
+                                        stake_inference_ratio=0.64)
+    assert deflated["compliant_frac"] < accurate["compliant_frac"]
+    assert deflated["max_compliant_stake"] <= accurate["max_compliant_stake"]
+
+
+def test_more_cover_traffic_raises_the_ceiling():
+    """The quota is the budget, so paying more cover traffic admits more concentrated stake."""
+    n, S = 20_000, 648_000
+    s = assign_stake(n, "zipf", np.random.default_rng(9), zipf_a=1.0)
+    lean = simulate_epoch_emissions(s, F, n, S, np.random.default_rng(10), cover_rate_mult=1.0)
+    rich = simulate_epoch_emissions(s, F, n, S, np.random.default_rng(10), cover_rate_mult=32.0)
+    assert rich["compliant_frac"] > lean["compliant_frac"]
+    assert rich["max_compliant_stake"] > lean["max_compliant_stake"]