Add pd: peering-degree Blend Monte-Carlo graph simulator

Static-graph simulator quantifying how a node's peering degree trades off
propagation speed, adversary exposure, deanonymization, and reliability in the
Blend network. Scales to 1e6 nodes (sparse CSR + sampled Dijkstra); the
adversary and deanonymization metrics are exact at every N.

Model (ms): seeded d-regular peer graph (matching-union), Blend cascade
(sender -> blend_hops timed-release mix relays -> final flood), geographic link
base + exponential transport jitter, per-node processing lag, free-running
release-clock mixing.

Metrics:
- propagation: full-delay mean/p50/p90/p99, path/broadcast split, coverage times
- reliability: message success-delivery-rate ~ (1-unresponsive_frac)^blend_hops
  and flood coverage, with unresponsive nodes modelled as routing holes
- adversary (exact): observed/eclipsed fractions, random + worst-case placement
- deanonymization (exact): P(whole blend path adversarial) ~ f_adv^blend_hops,
  and full deanonymization (path adversarial AND honest sender peered with an
  adversary) = deanon_rate * observed_frac

Deterministic blake2b seed streams, three parquet tables, joblib parallelism,
memguard, an analytic verify harness, 50 unit tests, and an auto-installing
Makefile.

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

tools/simulators/blend/pd/.gitignore

@@ -0,0 +1,17 @@
+# Generated artifacts
+runs/
+results/
+figures/
+!results/.gitkeep
+
+# Python
+__pycache__/
+*.py[cod]
+*.egg-info/
+.venv/
+venv/
+.pytest_cache/
+.mypy_cache/
+.ruff_cache/
+*.parquet
+*.csv

tools/simulators/blend/pd/Makefile

@@ -0,0 +1,45 @@
+VENV ?= .venv
+PY := $(VENV)/bin/python
+STAMP := $(VENV)/.installed
+
+# Keep numpy/scipy BLAS single-threaded so joblib process parallelism doesn't oversubscribe.
+export OMP_NUM_THREADS := 1
+export OPENBLAS_NUM_THREADS := 1
+export MKL_NUM_THREADS := 1
+export NUMEXPR_NUM_THREADS := 1
+
+.PHONY: install smoke sweep sweep-fullscale figures verify test lint clean
+
+# The stamp is the real install; targets below depend on it so `make sweep` (etc.) auto-installs
+# on a fresh checkout and re-installs whenever pyproject.toml changes.
+$(STAMP): pyproject.toml
+	python3 -m venv $(VENV)
+	$(PY) -m pip install -U pip
+	$(PY) -m pip install -e ".[dev]"
+	@touch $(STAMP)
+
+install: $(STAMP)
+
+smoke: $(STAMP) ## fast end-to-end (seconds): tiny N, few rounds/seeds
+	$(PY) -m pd.sweep --config configs/smoke.yaml
+
+sweep: $(STAMP)
+	$(PY) -m pd.sweep --config configs/default.yaml
+
+sweep-fullscale: $(STAMP)
+	$(PY) -m pd.sweep --config configs/fullscale.yaml
+
+figures: $(STAMP) ## make figures RUN=runs/<dir>
+	$(PY) -m pd.plotting.make_figures --run $(RUN)
+
+verify: $(STAMP)
+	$(PY) -m pd.verify
+
+test: $(STAMP)
+	$(PY) -m pytest
+
+lint: $(STAMP)
+	$(VENV)/bin/ruff check src scripts tests
+
+clean:
+	rm -rf runs/* figures/* .pytest_cache .ruff_cache .mypy_cache

tools/simulators/blend/pd/README.md

@@ -0,0 +1,68 @@
+# pd — peering-degree Monte-Carlo graph simulator
+
+Quantifies how a node's **peering degree** trades off, in the Blend network:
+
+- **propagation speed** — the full delay (ms) of a message: a random sender routes it along a
+  `blend_hops`-relay Blend path (each relay a free-running timed-release mix node) and the last
+  relay floods the whole network;
+- **adversary exposure** — with a fraction `f_adv` of adversarial nodes, how many honest nodes are
+  peered with ≥1 adversary (**observed**) and how many are fully surrounded (**eclipsed**);
+- **deanonymization** — tying propagation to the adversary: how often a message's *whole* blend path
+  is adversarial (**deanonymization** — the adversary owns the cascade end-to-end) and how often the
+  honest sender is *additionally* peered with an adversary (**full deanonymization** — the message is
+  tied back to its originator); and
+- **reliability under churn** — with a fraction `unresponsive_frac` of nodes that relay nothing, the
+  **message success-delivery-rate** (fraction of messages that survive the whole blend cascade to a
+  responsive final relay) and the flood **coverage** of those that do.
+
+The peer graph is a seeded random **d-regular** graph (exactly `degree` symmetric peers, identical
+for everyone from one global seed). This is static-graph analysis — no consensus — so it is far
+lighter than the TSI simulators and scales to **10⁶ nodes** (sparse CSR + sampled Dijkstra; the
+adversary metrics are exact at every N).
+
+## Model (all delays in ms)
+- **Link delay:** geographic base (metro 15 → antipodal 200 ms) + exponential transport jitter.
+- **Processing lag:** each node draws a fixed lag from a categorical distribution (default
+  {10, 50, 100} ms at {0.5, 0.4, 0.1}), incurred every time it relays.
+- **Blend mixing:** each relay releases on a free-running clock whose successive intervals are
+  Uniform{0…`max_blend_delay`} whole seconds; a held message waits for the relay's next release
+  (the renewal residual). Mixing happens only at the `blend_hops` relays; the final flood is plain.
+- **Unresponsive nodes:** a random `unresponsive_frac` of the population relays nothing (its outgoing
+  edges are removed). Relays are drawn from the whole node list *blind to responsiveness*, so a
+  message dies if any relay on its path is unresponsive — the delivery-rate then tracks
+  `(1−unresponsive_frac)^blend_hops`. Unresponsive nodes still *receive*, but they are routing holes,
+  so a delivered flood can strand pockets; a higher peering degree supplies redundant paths that keep
+  coverage high. This axis affects propagation only, not the adversary metrics.
+- **Deanonymization:** relays are picked *blind to who is adversarial*, so P(the whole blend path is
+  adversarial) is the exact hypergeometric `C(n_adv, blend_hops) / C(N−1, blend_hops)` ≈
+  `f_adv^blend_hops` (**deanon_rate**) — placement-independent, driven by path length, not degree.
+  Multiplying by the fraction of honest nodes with ≥1 adversary peer (`observed_frac`, which the
+  worst-case-coverage placement maximizes) gives **full_deanon_rate** — the honest sender is *also*
+  directly exposed, so the message is tied to its originator. Lengthening the blend path is the
+  dominant defence; a higher degree speeds propagation but *raises* the chance a sender directly
+  touches the adversary. Both are exact at every N (no Monte-Carlo), like the other adversary metrics.
+
+## Quick start
+```
+make install                 # or reuse a sibling venv: PYTHONPATH=src <python> -m pd.sweep ...
+make smoke                   # fast end-to-end -> runs/<ts>_smoke/{propagation,adversary}.parquet + figures/
+make verify                  # analytic checks (closed forms + graph invariants)
+make test                    # unit tests
+make sweep                   # configs/default.yaml (N up to 1e5, both adversary modes)
+make sweep-fullscale         # configs/fullscale.yaml (N up to 1e6, random-mode exact)
+make figures RUN=runs/<dir>
+```
+
+## Outputs
+Three parquets per run: `propagation.parquet` (`full_delay_ms_*`, `delivery_rate`, `frac_reached`,
+`coverN_ms` vs degree / blend_hops / N / unresponsive_frac), `adversary.parquet` (`observed_frac` /
+`eclipsed_frac` vs degree / f_adv / mode, random + worst-case envelope), and `deanon.parquet`
+(`deanon_rate` / `full_deanon_rate` vs degree / blend_hops / f_adv / mode — propagation paths crossed
+with the adversary set). Figures render all three, including delivery-rate and flood-coverage vs the
+unresponsive fraction and the deanonymization rates vs blend-path length, f_adv, and degree.
+
+## Layout
+`src/pd/`: `graph` (matching-union CSR d-regular), `propagation` (Blend cascade), `mixclock`
+(release-clock residual), `adversary` (exact observation/eclipse + deanonymization + placement),
+`config`/`engine`/`sweep`/`metrics`, `plotting`. `configs/` sweeps, `tests/`, `scripts/` shims.
+Reports of record live outside the sim at `reports/blend/pd/`.

tools/simulators/blend/pd/configs/default.yaml

@@ -0,0 +1,13 @@
+# Main study grid (N up to 1e5; both adversary modes; degree is the primary axis).
+n_nodes: [1000, 10000, 100000]
+degree: [3, 4, 6, 8, 12, 16]
+blend_hops: [1, 2, 3, 5]
+max_blend_delay: [3]
+unresponsive_frac: [0.0, 0.05, 0.1, 0.2, 0.3, 0.5]
+f_adv: [0.05, 0.1, 0.2, 0.33, 0.5]
+adversary_mode: [random, worstcase_coverage, worstcase_eclipse]
+seeds: 8
+base:
+  n_rounds: 200
+  n_placements: 8
+  worstcase_max_n: 100000

tools/simulators/blend/pd/configs/fullscale.yaml

@@ -0,0 +1,14 @@
+# Internet-scale: N up to 1e6. Propagation is sampled (n_rounds); adversary metrics are exact.
+# Worst-case placement is capped (worstcase_max_n) so only `random` runs at 1e6.
+n_nodes: [100000, 1000000]
+degree: [4, 6, 8, 12, 16]
+blend_hops: [3]
+max_blend_delay: [3]
+unresponsive_frac: [0.0, 0.1, 0.2, 0.5]
+f_adv: [0.1, 0.2, 0.33]
+adversary_mode: [random]
+seeds: 3
+base:
+  n_rounds: 64
+  n_placements: 3
+  worstcase_max_n: 100000

tools/simulators/blend/pd/configs/smoke.yaml

@@ -0,0 +1,12 @@
+# Fast end-to-end smoke run (seconds): tiny N, few rounds/seeds, both adversary modes.
+n_nodes: [1000]
+degree: [4, 8]
+blend_hops: [2, 3]
+max_blend_delay: [0, 3]
+unresponsive_frac: [0.0, 0.2, 0.5]
+f_adv: [0.1, 0.33]
+adversary_mode: [random, worstcase_coverage]
+seeds: 2
+base:
+  n_rounds: 20
+  n_placements: 2

tools/simulators/blend/pd/pyproject.toml

@@ -0,0 +1,43 @@
+[build-system]
+requires = ["hatchling"]
+build-backend = "hatchling.build"
+
+[project]
+name = "pd"
+version = "0.1.0"
+description = "Peering-degree Monte-Carlo graph simulator for the Blend network"
+requires-python = ">=3.11"
+dependencies = [
+    "numpy",
+    "pandas",
+    "pyarrow",
+    "matplotlib",
+    "scipy",
+    "pyyaml",
+    "tqdm",
+    "joblib",
+]
+
+[project.optional-dependencies]
+dev = ["pytest", "pytest-xdist", "ruff", "mypy"]
+
+[project.scripts]
+pd-sweep = "pd.sweep:main"
+pd-verify = "pd.verify:main"
+pd-figures = "pd.plotting.make_figures:main"
+
+[tool.hatch.build.targets.wheel]
+packages = ["src/pd"]
+
+[tool.ruff]
+line-length = 100
+target-version = "py311"
+
+[tool.ruff.lint]
+select = ["E", "F", "I", "UP", "B", "NPY"]
+ignore = ["E741"]
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
+addopts = "-q"
+markers = ["slow: marks slow tests (deselect with -m 'not slow')"]

tools/simulators/blend/pd/requirements-dev.txt

@@ -0,0 +1 @@
+-e .[dev]

tools/simulators/blend/pd/requirements.txt

@@ -0,0 +1 @@
+-e .

tools/simulators/blend/pd/scripts/make_figures.py

@@ -0,0 +1,11 @@
+#!/usr/bin/env python3
+"""Shim: add src/ to sys.path, then run pd.plotting.make_figures:main."""
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
+
+from pd.plotting.make_figures import main  # noqa: E402
+
+if __name__ == "__main__":
+    raise SystemExit(main())

tools/simulators/blend/pd/scripts/run_sweep.py

@@ -0,0 +1,11 @@
+#!/usr/bin/env python3
+"""Shim: add src/ to sys.path, then run pd.sweep:main."""
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
+
+from pd.sweep import main  # noqa: E402
+
+if __name__ == "__main__":
+    raise SystemExit(main())

tools/simulators/blend/pd/scripts/verify.py

@@ -0,0 +1,11 @@
+#!/usr/bin/env python3
+"""Shim: add src/ to sys.path, then run pd.verify:main."""
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
+
+from pd.verify import main  # noqa: E402
+
+if __name__ == "__main__":
+    raise SystemExit(main())

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

@@ -0,0 +1,3 @@
+"""pd — peering-degree Monte-Carlo graph simulator for the Blend network."""
+
+__version__ = "0.1.0"

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

@@ -0,0 +1,162 @@
+"""Adversary observation + eclipse metrics (exact, O(N*degree)) and placement strategies.
+
+- ``adversary_metrics``: given a boolean adversary mask, count the honest nodes peered with >=1
+  adversary (**observed**) and the honest nodes whose EVERY peer is adversarial (**eclipsed**),
+  via one sparse reduction over the CSR. Exact at every N (incl. 1e6).
+- ``place_adversary``: pick the adversary set. ``random`` (average case) scales to 1e6; the
+  worst-case greedy strategies (``worstcase_coverage``/``worstcase_eclipse``) characterize the
+  security *envelope* and are meant for N <= ``worstcase_max_n``. ``worstcase_degree`` is
+  degenerate on a strict d-regular graph (all degrees equal) and coincides with ``random``.
+"""
+
+from __future__ import annotations
+
+import heapq
+
+import numpy as np
+
+from .config import WORSTCASE_MODES
+from .graph import Graph
+
+
+def adversary_metrics(graph: Graph, adv_mask: np.ndarray) -> dict:
+    """Exact observation/eclipse counts for the given adversary mask."""
+    adv = adv_mask
+    honest = ~adv
+    counts = np.add.reduceat(adv[graph.indices].astype(np.int32), graph.indptr[:-1])
+    observed = honest & (counts >= 1)
+    eclipsed = honest & (counts == graph.degree)
+    n_adv = int(adv.sum())
+    n_honest = int(honest.sum())
+    obs = int(observed.sum())
+    ecl = int(eclipsed.sum())
+    return {
+        "n_adv": n_adv,
+        "n_honest": n_honest,
+        "observed_count": obs,
+        "observed_frac": float(obs / n_honest) if n_honest else 0.0,
+        "eclipsed_count": ecl,
+        "eclipsed_frac": float(ecl / n_honest) if n_honest else 0.0,
+        "mean_adv_peers_honest": float(counts[honest].mean()) if n_honest else 0.0,
+        "max_adv_peers": int(counts.max()) if graph.n else 0,
+    }
+
+
+def deanon_metrics(n: int, n_adv: int, observed_frac: float, blend_hops: int) -> dict:
+    """Exact deanonymization rates for an honest sender whose message traverses ``blend_hops``
+    relays chosen uniformly *blind to who is adversarial* (the sender cannot know who is honest).
+
+    * ``deanon_rate`` — P(**every** relay on the path is adversarial): the adversary then controls
+      the whole blend cascade and links the message from entry to exit. Because the relays are a
+      uniform draw, this is the hypergeometric ``C(n_adv, k) / C(n-1, k)`` (an honest sender leaves
+      all ``n_adv`` adversaries in the ``n-1``-node relay pool). It therefore depends only on the
+      adversary *count*, not the placement, and ~ ``f_adv**blend_hops``.
+    * ``full_deanon_rate`` — additionally the honest sender is directly peered with >=1 adversary,
+      so the adversary also ties the message to its originator. The sender is uniform over honest
+      nodes, so this is ``deanon_rate * observed_frac``; ``observed_frac`` (the honest-node fraction
+      with an adversary peer) is placement-dependent, so full deanonymization carries the
+      placement's fingerprint (worst-case coverage drives it up).
+
+    Both are exact at every N -- no Monte-Carlo -- matching the other adversary metrics.
+    """
+    k = int(blend_hops)
+    deanon = 0.0
+    if 1 <= k <= n_adv and k <= n - 1:
+        deanon = 1.0
+        for i in range(k):
+            deanon *= (n_adv - i) / (n - 1 - i)   # C(n_adv,k)/C(n-1,k), stable for small k
+    return {"deanon_rate": float(deanon), "full_deanon_rate": float(deanon * observed_frac)}
+
+
+def _peers(graph: Graph, v: int) -> np.ndarray:
+    return graph.indices[graph.indptr[v]:graph.indptr[v + 1]]
+
+
+def _greedy_coverage(graph: Graph, n_adv: int, rng: np.random.Generator) -> np.ndarray:
+    """Lazy-greedy (CELF) maximum-coverage: pick adversaries covering the most honest neighbours."""
+    n, degree = graph.n, graph.degree
+    adv = np.zeros(n, dtype=bool)
+    covered = np.zeros(n, dtype=bool)
+
+    def gain(v: int) -> int:
+        pu = _peers(graph, v)
+        return int(np.count_nonzero(~adv[pu] & ~covered[pu]))
+
+    heap = [(-degree, int(v)) for v in range(n)]
+    heapq.heapify(heap)
+    chosen: list[int] = []
+    while len(chosen) < n_adv and heap:
+        neg, v = heapq.heappop(heap)
+        if adv[v]:
+            continue
+        g = gain(v)
+        if -neg != g:
+            heapq.heappush(heap, (-g, v))
+            continue
+        adv[v] = True
+        chosen.append(v)
+        pu = _peers(graph, v)
+        covered[pu[~adv[pu] & ~covered[pu]]] = True
+    return np.asarray(chosen, dtype=np.int64)
+
+
+def _greedy_eclipse(graph: Graph, n_adv: int, rng: np.random.Generator) -> np.ndarray:
+    """Greedy fully-surround: repeatedly finish the honest node closest to eclipse.
+
+    Heuristic (eclipse maximization is NP-hard) — an upper-envelope, not exact.
+    """
+    n, degree = graph.n, graph.degree
+    adv = np.zeros(n, dtype=bool)
+    remaining = np.full(n, degree, dtype=np.int64)   # honest node's peers not yet adversary
+    budget = int(n_adv)
+    heap = [(degree, int(u)) for u in range(n)]
+    heapq.heapify(heap)
+    while budget > 0 and heap:
+        rem_u, u = heapq.heappop(heap)
+        if adv[u]:
+            continue
+        cur = int(remaining[u])
+        if cur != rem_u:
+            heapq.heappush(heap, (cur, u))
+            continue
+        if cur == 0:
+            continue
+        pu = _peers(graph, u)
+        tozap = pu[~adv[pu]]
+        if tozap.shape[0] > budget:
+            tozap = tozap[:budget]
+        for z in tozap.tolist():
+            if adv[z]:
+                continue
+            adv[z] = True
+            budget -= 1
+            remaining[_peers(graph, z)] -= 1
+    return np.where(adv)[0]
+
+
+def place_adversary(graph: Graph, f_adv: float, mode: str, rng: np.random.Generator,
+                    worstcase_max_n: int) -> np.ndarray:
+    """Return a boolean adversary mask of size ``graph.n``."""
+    n = graph.n
+    n_adv = int(round(f_adv * n))
+    mask = np.zeros(n, dtype=bool)
+    if n_adv <= 0:
+        return mask
+    if n_adv >= n:
+        mask[:] = True
+        return mask
+    if mode in ("random", "worstcase_degree"):
+        # d-regular: all degrees equal, so degree-ranked == uniform-random.
+        mask[rng.choice(n, size=n_adv, replace=False)] = True
+        return mask
+    if mode not in WORSTCASE_MODES:
+        raise ValueError(f"unknown adversary_mode {mode!r}")
+    if n > worstcase_max_n:
+        raise ValueError(
+            f"worst-case mode {mode!r} at N={n} exceeds worstcase_max_n={worstcase_max_n}; "
+            f"run worst-case only at smaller N (the engine skips it above the cap).")
+    if mode == "worstcase_coverage":
+        mask[_greedy_coverage(graph, n_adv, rng)] = True
+    else:  # worstcase_eclipse
+        mask[_greedy_eclipse(graph, n_adv, rng)] = True
+    return mask

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

@@ -0,0 +1,151 @@
+"""Configuration dataclasses for single runs and parameter sweeps."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+from typing import Any, Literal
+
+from . import constants
+
+LatencyDist = Literal["geo", "fixed", "uniform", "exp"]
+AdversaryMode = Literal[
+    "random", "worstcase_coverage", "worstcase_eclipse", "worstcase_degree"
+]
+_DISTS = ("geo", "fixed", "uniform", "exp")
+_MODES = ("random", "worstcase_coverage", "worstcase_eclipse", "worstcase_degree")
+WORSTCASE_MODES = ("worstcase_coverage", "worstcase_eclipse", "worstcase_degree")
+
+
+@dataclass(frozen=True)
+class SimConfig:
+    """A single fully-specified graph cell (one topology + one propagation/adversary setting)."""
+
+    # --- network ---
+    n_nodes: int = 1000                 # must be even (matching-union construction)
+    degree: int = 8                     # peering degree — the primary study axis
+
+    # --- propagation (Blend cascade, delays in ms) ---
+    blend_hops: int = 3                 # relay-path length (swept)
+    max_blend_delay: int = 3            # free-running release-clock max interval, whole SECONDS
+    unresponsive_frac: float = 0.0      # ratio of nodes that do NOT relay any messages (swept)
+    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)
+    processing_lag_probs: tuple[float, ...] = (0.5, 0.4, 0.1)
+    link_latency_dist: LatencyDist = "geo"
+    link_latency_mean_ms: float = constants.GEO_LATENCY_MEAN_MS   # only used by non-geo dists
+    coverage_pcts: tuple[float, ...] = (50.0, 90.0, 99.0)
+
+    # --- adversary ---
+    f_adv: float = 0.0
+    adversary_mode: AdversaryMode = "random"
+    n_placements: int = 4               # random-placement sub-replicates (worstcase -> 1)
+    worstcase_max_n: int = 100_000      # cap greedy worst-case strategies above this N
+
+    # --- bookkeeping ---
+    graph_seed: int = 0                 # the "global seed" / topology-ensemble index
+    replicate: int = 0
+    root_seed: int = 12345
+
+    def __post_init__(self) -> None:
+        if self.n_nodes % 2 != 0:
+            raise ValueError(f"n_nodes must be even, got {self.n_nodes}")
+        if not (1 <= self.degree < self.n_nodes):
+            raise ValueError(f"need 1 <= degree < n_nodes ({self.degree} vs {self.n_nodes})")
+        if not (1 <= self.blend_hops < self.n_nodes):
+            raise ValueError(f"need 1 <= blend_hops < n_nodes, got {self.blend_hops}")
+        if self.max_blend_delay < 0:
+            raise ValueError("max_blend_delay must be >= 0 (whole seconds)")
+        if not (0.0 <= self.unresponsive_frac < 1.0):
+            raise ValueError(f"need 0 <= unresponsive_frac < 1, got {self.unresponsive_frac}")
+        if not (0.0 <= self.f_adv < 1.0):
+            raise ValueError(f"need 0 <= f_adv < 1, got {self.f_adv}")
+        if self.n_rounds < 1 or self.n_placements < 1:
+            raise ValueError("n_rounds and n_placements must be >= 1")
+        if len(self.processing_lags_ms) != len(self.processing_lag_probs):
+            raise ValueError("processing_lags_ms and processing_lag_probs must have equal length")
+        if abs(sum(self.processing_lag_probs) - 1.0) > 1e-9:
+            raise ValueError(
+                f"processing_lag_probs must sum to 1, got {sum(self.processing_lag_probs)}")
+        if any(p < 0 for p in self.processing_lag_probs):
+            raise ValueError("processing_lag_probs must be non-negative")
+        if self.link_latency_dist not in _DISTS:
+            raise ValueError(f"link_latency_dist must be one of {_DISTS}")
+        if self.adversary_mode not in _MODES:
+            raise ValueError(f"adversary_mode must be one of {_MODES}")
+
+    @property
+    def n_adv(self) -> int:
+        return int(round(self.f_adv * self.n_nodes))
+
+    @property
+    def n_honest(self) -> int:
+        return self.n_nodes - self.n_adv
+
+    def key(self) -> tuple:
+        """Hashable identity used to seed RNGs deterministically."""
+        return (
+            self.n_nodes, self.degree, self.blend_hops, self.max_blend_delay,
+            self.unresponsive_frac, 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,
+            self.graph_seed, self.replicate, self.root_seed,
+        )
+
+
+_TUPLE_FIELDS = ("processing_lags_ms", "processing_lag_probs", "coverage_pcts")
+
+
+@dataclass
+class SweepConfig:
+    """Grids for a sweep: topologies (n x degree x graph_seed) plus the propagation and
+    adversary sub-grids that every topology is measured over."""
+
+    n_nodes: list[int] = field(default_factory=lambda: [1000])
+    degree: list[int] = field(default_factory=lambda: [3, 4, 6, 8, 12, 16])
+    blend_hops: list[int] = field(default_factory=lambda: [3])
+    max_blend_delay: list[int] = field(default_factory=lambda: [3])
+    unresponsive_frac: list[float] = field(default_factory=lambda: [0.0])
+    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)
+    base: dict[str, Any] = field(default_factory=dict)
+
+    def base_config(self, n_nodes: int, degree: int, graph_seed: int) -> SimConfig:
+        b = dict(self.base)
+        for k in _TUPLE_FIELDS:
+            if k in b and b[k] is not None:
+                b[k] = tuple(b[k])
+        b.pop("n_nodes", None)
+        b.pop("degree", None)
+        b.pop("graph_seed", None)
+        return SimConfig(n_nodes=n_nodes, degree=degree, graph_seed=graph_seed, **b)
+
+    def graph_cells(self) -> list[tuple[int, int, int]]:
+        """Distinct topologies to build: (n_nodes, degree, graph_seed)."""
+        return [(n, d, g) for n in self.n_nodes for d in self.degree for g in range(self.seeds)]
+
+    def prop_grid(self) -> list[tuple[int, int]]:
+        """(blend_hops, max_blend_delay) settings each topology is measured over."""
+        return [(bh, md) for bh in self.blend_hops for md in self.max_blend_delay]
+
+    def adv_grid(self) -> list[tuple[float, str]]:
+        """(f_adv, adversary_mode) settings; f_adv==0 keeps only one (mode-irrelevant) row."""
+        out: list[tuple[float, str]] = []
+        for f in self.f_adv:
+            modes = self.adversary_mode if f > 0 else self.adversary_mode[:1]
+            for m in modes:
+                out.append((f, m))
+        return out
+
+    @classmethod
+    def from_dict(cls, d: dict[str, Any]) -> SweepConfig:
+        d = dict(d)
+        base = d.pop("base", {})
+        known = {"n_nodes", "degree", "blend_hops", "max_blend_delay", "unresponsive_frac",
+                 "f_adv", "adversary_mode", "seeds"}
+        unknown = set(d) - known
+        if unknown:
+            raise ValueError(f"unknown sweep keys: {sorted(unknown)}")
+        return cls(base=base, **{k: v for k, v in d.items() if k in known})

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

@@ -0,0 +1,24 @@
+"""Network-latency constants for pd. All delays are in **milliseconds**.
+
+One-way, application-level latency between two directly-peered nodes, bucketed by the
+geographic relationship of the peers (~ RTT/2 from public latency measurements plus a little
+gossip processing overhead). In a globally distributed node set a random peer is usually on
+another continent, so most peer links fall in the long-latency bands. A message gossip-floods
+over the peering graph, so its delay to a far node is the sum of a few such per-link latencies
+along the fastest (Dijkstra) path — see graph.py / propagation.py.
+"""
+
+from __future__ import annotations
+
+GEO_LATENCY_BANDS_MS = (
+    15.0,    # metro / same country            (~15 ms one-way)
+    40.0,    # same continent, e.g. EU<->EU    (~40 ms)
+    90.0,    # transatlantic, e.g. EU<->US-East (~90 ms)
+    200.0,   # antipodal, e.g. EU<->AU / EU<->JP (~200 ms)
+)
+# Share of random peer links falling in each band for a globally distributed node set.
+GEO_LATENCY_WEIGHTS = (0.15, 0.35, 0.35, 0.15)
+# Mean one-way latency of a random global peer link under the mixture above (~78.75 ms).
+GEO_LATENCY_MEAN_MS = sum(
+    b * w for b, w in zip(GEO_LATENCY_BANDS_MS, GEO_LATENCY_WEIGHTS, strict=True)
+)

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

@@ -0,0 +1,81 @@
+"""Build-once engine: one topology -> propagation rows (per blend setting) + adversary rows.
+
+The graph build and propagation depend only on ``(n_nodes, degree, graph_seed)``; the adversary
+metrics are cheap and exact. So a topology is built once and measured across the whole
+propagation and adversary sub-grids.
+"""
+
+from __future__ import annotations
+
+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 .propagation import assign_responsive, propagation_metrics
+from .rng import placement_seedseq, responsive_seedseq, round_seedseq
+
+
+def run_graph_cell(base: SimConfig, prop_grid: list[tuple[int, int]],
+                   unresponsive_fracs: list[float],
+                   adv_grid: list[tuple[float, str]],
+                   ) -> tuple[list[dict], list[dict], list[dict]]:
+    """Build ``base``'s topology once; return (propagation, adversary, deanonymization 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 (propagation-only), ``adv_grid`` = [(f_adv, mode)]. Deanonymization crosses each adversary
+    placement with the propagation grid's blend-path lengths, so it is emitted alongside the
+    adversary rows.
+    """
+    graph = build_graph(base)
+    blend_hops_set = sorted({bh for bh, _ in prop_grid})
+
+    prop_rows: list[dict] = []
+    for uf in unresponsive_fracs:
+        responsive = assign_responsive(
+            base.n_nodes, uf, np.random.default_rng(responsive_seedseq(base, uf)))
+        for blend_hops, max_blend_delay in prop_grid:
+            rng = np.random.default_rng(round_seedseq(base, blend_hops, max_blend_delay, uf))
+            prop = propagation_metrics(
+                graph, blend_hops, max_blend_delay, uf, responsive, base, rng)
+            prop_rows.append(propagation_row(base, blend_hops, max_blend_delay, uf, prop))
+
+    adv_rows: list[dict] = []
+    deanon_rows: list[dict] = []
+    for f_adv, mode in adv_grid:
+        if mode in WORSTCASE_MODES and base.n_nodes > base.worstcase_max_n:
+            continue  # worst-case is an envelope characterized at N <= worstcase_max_n
+        n_placements = base.n_placements if mode == "random" else 1
+        for rep in range(n_placements):
+            rng = np.random.default_rng(placement_seedseq(base, f_adv, mode, rep))
+            adv_mask = place_adversary(graph, f_adv, mode, rng, base.worstcase_max_n)
+            adv = adversary_metrics(graph, adv_mask)
+            adv_rows.append(adversary_row(base, f_adv, mode, rep, adv))
+            for bh in blend_hops_set:
+                dz = deanon_metrics(graph.n, adv["n_adv"], adv["observed_frac"], bh)
+                deanon_rows.append(deanon_row(base, bh, f_adv, mode, rep, adv, dz))
+
+    return prop_rows, adv_rows, deanon_rows
+
+
+def run_trajectory(config: SimConfig) -> dict:
+    """Single-cell convenience for tests/verify: build the graph, run one propagation cell
+    (``config.blend_hops``/``config.max_blend_delay``) and one adversary cell."""
+    graph = build_graph(config)
+    uf = config.unresponsive_frac
+    responsive = assign_responsive(
+        config.n_nodes, uf, np.random.default_rng(responsive_seedseq(config, uf)))
+    prng = np.random.default_rng(
+        round_seedseq(config, config.blend_hops, config.max_blend_delay, uf))
+    prop = propagation_metrics(
+        graph, config.blend_hops, config.max_blend_delay, uf, responsive, config, prng)
+    arng = np.random.default_rng(
+        placement_seedseq(config, config.f_adv, config.adversary_mode, config.replicate))
+    adv_mask = place_adversary(
+        graph, config.f_adv, config.adversary_mode, arng, config.worstcase_max_n)
+    adv = adversary_metrics(graph, adv_mask)
+    deanon = deanon_metrics(graph.n, adv["n_adv"], adv["observed_frac"], config.blend_hops)
+    return {"graph": graph, "propagation": prop, "adversary": adv, "adv_mask": adv_mask,
+            "deanon": deanon}

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

@@ -0,0 +1,179 @@
+"""Deterministic random d-regular peer graph as a sparse CSR, scalable to 1e6 nodes.
+
+Construction: a **union of `degree` random perfect matchings** (fully vectorized), repaired to
+drop the few parallel edges (expected count ~ C(degree,2), independent of N) so the result is
+EXACTLY d-regular, simple, undirected — and reconstructible from one global seed. Preferred over
+the tsi ring-lattice + Maslov-Sneppen generator: O(N*degree) (~seconds at 1e6 vs minutes of
+Python swaps) and genuinely low-diameter immediately.
+
+The ``Graph`` holds the undirected CSR adjacency (``indptr``/``indices`` — drives the exact
+adversary reductions) plus the directed per-edge base latency (``base``) and source node
+(``src``), and the per-node processing lags (``p``), which drive the per-round directed Dijkstra
+in ``propagation.py``.
+"""
+
+from __future__ import annotations
+
+from collections import defaultdict
+from dataclasses import dataclass
+
+import numpy as np
+from scipy.sparse import csr_matrix
+
+from . import latency
+from .config import SimConfig
+from .memguard import check_alloc
+from .rng import graph_seedseq
+
+
+@dataclass(frozen=True)
+class Graph:
+    n: int
+    degree: int
+    indptr: np.ndarray      # (n+1,)  CSR row pointers (contiguous degree-length blocks)
+    indices: np.ndarray     # (2E,)   peers of each node, grouped by source; 2E = n*degree
+    base: np.ndarray        # (2E,)   geo base latency (ms) per directed edge, CSR-aligned
+    src: np.ndarray         # (2E,)   source node per directed edge, CSR-aligned
+    p: np.ndarray           # (n,)    per-node processing lag (ms)
+
+    @property
+    def n_edges(self) -> int:
+        return int(self.indices.shape[0] // 2)
+
+    def weighted_csr(self, data: np.ndarray) -> csr_matrix:
+        """A directed CSR sharing this graph's sparsity, with the given per-edge ``data``."""
+        return csr_matrix((data, self.indices, self.indptr), shape=(self.n, self.n))
+
+
+def build_regular_edges(n: int, degree: int, rng: np.random.Generator) -> np.ndarray:
+    """Exactly d-regular simple undirected edge list, shape (E, 2) with u < v.
+
+    Requires ``n`` even and ``1 <= degree < n``. Deterministic in ``rng``.
+    """
+    if n % 2 != 0:
+        raise ValueError("n must be even")
+    if not (1 <= degree < n):
+        raise ValueError("need 1 <= degree < n")
+    h = n // 2
+    lo = np.empty(degree * h, dtype=np.int64)
+    hi = np.empty(degree * h, dtype=np.int64)
+    for m in range(degree):
+        perm = rng.permutation(n)
+        a, b = perm[0::2], perm[1::2]
+        lo[m * h:(m + 1) * h] = np.minimum(a, b)
+        hi[m * h:(m + 1) * h] = np.maximum(a, b)
+    # drop parallel edges (keep first occurrence of each undirected pair)
+    key = lo * n + hi
+    order = np.argsort(key, kind="stable")
+    ks = key[order]
+    first = np.ones(ks.shape[0], dtype=bool)
+    first[1:] = ks[1:] != ks[:-1]
+    edges = np.stack([lo[order][first], hi[order][first]], axis=1)
+    deg = np.bincount(edges.ravel(), minlength=n)
+    if np.all(deg == degree):
+        return edges
+    return _repair_to_regular(edges, deg, n, degree, rng)
+
+
+def _repair_to_regular(edges_arr: np.ndarray, deg: np.ndarray, n: int, degree: int,
+                       rng: np.random.Generator, max_passes: int = 100_000) -> np.ndarray:
+    """Restore exact d-regularity after dedup dropped a few parallel edges.
+
+    Re-pairs the freed stubs (a tiny set, size ~ 2*C(degree,2)) into new simple edges, using a
+    single Maslov-Sneppen swap against a full-degree edge whenever a pass makes no progress.
+    """
+    deficient = np.where(deg < degree)[0]
+    defset = set(int(x) for x in deficient)
+    nbr: dict[int, set[int]] = defaultdict(set)
+    mask = np.isin(edges_arr[:, 0], deficient) | np.isin(edges_arr[:, 1], deficient)
+    for a, b in edges_arr[mask].tolist():
+        if a in defset:
+            nbr[a].add(b)
+        if b in defset:
+            nbr[b].add(a)
+    edges: list[list[int]] = edges_arr.tolist()
+    deg = deg.astype(np.int64).copy()
+
+    def _swap_connect(x: int, y: int) -> bool:
+        for _ in range(500):
+            idx = int(rng.integers(0, len(edges)))
+            a, b = edges[idx]
+            if deg[a] != degree or deg[b] != degree:   # only rewire full-degree edges
+                continue
+            if a in (x, y) or b in (x, y):
+                continue
+            if a in nbr[x] or b in nbr[y]:
+                continue
+            # remove (a,b); add (x,a) and (y,b) -> x,y gain a peer; a,b unchanged
+            if a in defset:
+                nbr[a].discard(b)
+                nbr[a].add(x)
+            if b in defset:
+                nbr[b].discard(a)
+                nbr[b].add(y)
+            edges[idx] = [min(x, a), max(x, a)]
+            edges.append([min(y, b), max(y, b)])
+            nbr[x].add(a)
+            nbr[y].add(b)
+            deg[x] += 1
+            deg[y] += 1
+            return True
+        return False
+
+    for _ in range(max_passes):
+        stubs: list[int] = []
+        for node in defset:
+            stubs.extend([node] * int(degree - deg[node]))
+        if not stubs:
+            break
+        rng.shuffle(stubs)
+        progressed = False
+        leftover: list[int] = []
+        i = 0
+        while i < len(stubs) - 1:
+            x, y = int(stubs[i]), int(stubs[i + 1])
+            i += 2
+            if x != y and y not in nbr[x]:
+                edges.append([min(x, y), max(x, y)])
+                nbr[x].add(y)
+                nbr[y].add(x)
+                deg[x] += 1
+                deg[y] += 1
+                progressed = True
+            else:
+                leftover.extend([x, y])
+        if i == len(stubs) - 1:
+            leftover.append(int(stubs[-1]))
+        if not progressed and leftover:
+            x = leftover[0]
+            y = leftover[1] if len(leftover) > 1 else leftover[0]
+            if not _swap_connect(x, y):
+                raise RuntimeError("graph repair swap failed")
+
+    out = np.asarray(edges, dtype=np.int64)
+    if not np.all(np.bincount(out.ravel(), minlength=n) == degree):
+        raise RuntimeError("graph repair failed to reach exact d-regularity")
+    return out
+
+
+def build_graph(config: SimConfig) -> Graph:
+    """Build the exact d-regular peer graph + directed base/src arrays + per-node lags.
+
+    Everything is a pure function of ``graph_seedseq(config)`` (topology fields only).
+    """
+    n, degree = config.n_nodes, config.degree
+    check_alloc(int(n * degree * 24), "d-regular CSR (indices+base+src)",
+                f"N={n}, degree={degree}")
+    rng = np.random.default_rng(graph_seedseq(config))
+    edges = build_regular_edges(n, degree, rng)
+    base_u = latency.sample_link_latencies(edges.shape[0], config, rng)
+    p = latency.assign_processing_lags(n, config, rng)
+    rows = np.concatenate([edges[:, 0], edges[:, 1]])
+    cols = np.concatenate([edges[:, 1], edges[:, 0]])
+    data = np.concatenate([base_u, base_u])
+    csr = csr_matrix((data, (rows, cols)), shape=(n, n))
+    csr.sort_indices()
+    indptr = csr.indptr.astype(np.int64)
+    src = np.repeat(np.arange(n, dtype=np.int64), np.diff(indptr))
+    return Graph(n=n, degree=degree, indptr=indptr, indices=csr.indices,
+                 base=csr.data, src=src, p=p)

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

@@ -0,0 +1,39 @@
+"""Per-link base latencies (ms) and per-node processing lags (ms)."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from . import constants
+from .config import SimConfig
+
+
+def sample_link_latencies(n_edges: int, config: SimConfig, rng: np.random.Generator) -> np.ndarray:
+    """Per-(undirected-)edge one-way base latency in ms.
+
+    ``geo`` draws each link from the real-world geographic band mixture (used as-is, in ms);
+    ``fixed``/``uniform``/``exp`` have mean ``link_latency_mean_ms``.
+    """
+    dist = config.link_latency_dist
+    if dist == "geo":
+        bands = np.asarray(constants.GEO_LATENCY_BANDS_MS, dtype=float)
+        weights = np.asarray(constants.GEO_LATENCY_WEIGHTS, dtype=float)
+        idx = rng.choice(bands.shape[0], size=n_edges, p=weights)
+        return bands[idx]
+    mean = config.link_latency_mean_ms
+    if dist == "fixed":
+        return np.full(n_edges, mean, dtype=float)
+    if dist == "uniform":
+        return rng.uniform(0.0, 2.0 * mean, size=n_edges)
+    if dist == "exp":
+        return rng.exponential(mean, size=n_edges)
+    raise ValueError(f"unknown link_latency_dist {dist!r}")
+
+
+def assign_processing_lags(n_nodes: int, config: SimConfig, rng: np.random.Generator) -> np.ndarray:
+    """Per-node fixed processing lag (ms): a categorical draw over ``processing_lags_ms`` with
+    weights ``processing_lag_probs`` (e.g. {10,50,100} ms at {0.5,0.4,0.1})."""
+    lags = np.asarray(config.processing_lags_ms, dtype=float)
+    probs = np.asarray(config.processing_lag_probs, dtype=float)
+    idx = rng.choice(lags.shape[0], size=n_nodes, p=probs)
+    return lags[idx]

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

@@ -0,0 +1,55 @@
+"""Fail-loud memory guard for large allocations (copied from tsi-sim-pernode).
+
+The dominant arrays here are the sparse CSR adjacency (``2E = N*degree`` entries) and the
+sampled single-source distance matrices (``S x N`` or ``(blend_hops+1) x N``), which grow with
+``N``. Every worker checks the size *before* allocating and raises ``AllocationTooLarge`` if it
+would exceed its budget, so an under-sized config fails with a clear message instead of freezing
+the machine.
+
+Budget (``budget_bytes``): ``PD_BYTES_BUDGET`` > 0 -> that many bytes (the sweep sets this to
+each worker's RAM share); otherwise ``DEFAULT_BUDGET_FRAC`` of physical RAM.
+"""
+
+from __future__ import annotations
+
+import os
+import subprocess
+
+DEFAULT_BUDGET_FRAC = 0.9
+
+
+class AllocationTooLarge(MemoryError):
+    """A pd array would exceed the memory budget; raised before allocating."""
+
+
+def total_ram_bytes() -> int:
+    """Best-effort physical RAM in bytes (POSIX sysconf, then Darwin sysctl, then 8 GB)."""
+    try:
+        return os.sysconf("SC_PAGE_SIZE") * os.sysconf("SC_PHYS_PAGES")
+    except (ValueError, OSError, AttributeError):
+        pass
+    try:  # macOS lacks SC_PHYS_PAGES
+        out = subprocess.run(["sysctl", "-n", "hw.memsize"], capture_output=True, text=True)
+        return int(out.stdout.strip())
+    except (OSError, ValueError):
+        return 8 * 1024**3
+
+
+def budget_bytes() -> int:
+    """Per-process byte budget for a single big array (see module docstring)."""
+    try:
+        explicit = int(os.environ.get("PD_BYTES_BUDGET", "0"))
+    except ValueError:
+        explicit = 0
+    if explicit > 0:
+        return explicit
+    return int(DEFAULT_BUDGET_FRAC * total_ram_bytes())
+
+
+def check_alloc(nbytes: int, label: str, detail: str = "") -> None:
+    """Raise ``AllocationTooLarge`` if allocating ``nbytes`` would exceed the budget."""
+    budget = budget_bytes()
+    if nbytes > budget:
+        raise AllocationTooLarge(
+            f"{label} needs {nbytes / 1024**3:.1f} GB > per-process budget "
+            f"{budget / 1024**3:.1f} GB.{(' ' + detail) if detail else ''}")

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

@@ -0,0 +1,57 @@
+"""Flat parquet-row builders for the two result tables."""
+
+from __future__ import annotations
+
+from .config import SimConfig
+
+
+def propagation_row(config: SimConfig, blend_hops: int, max_blend_delay: int,
+                    unresponsive_frac: float, prop: dict) -> dict:
+    return {
+        "n_nodes": config.n_nodes,
+        "degree": config.degree,
+        "blend_hops": blend_hops,
+        "max_blend_delay": max_blend_delay,
+        "unresponsive_frac": unresponsive_frac,
+        "graph_seed": config.graph_seed,
+        "n_rounds": config.n_rounds,
+        "transport_jitter_mean_ms": config.transport_jitter_mean_ms,
+        "processing_lags_ms": str(tuple(config.processing_lags_ms)),
+        "processing_lag_probs": str(tuple(config.processing_lag_probs)),
+        **prop,
+    }
+
+
+def adversary_row(config: SimConfig, f_adv: float, mode: str, placement_rep: int,
+                  adv: dict) -> dict:
+    return {
+        "n_nodes": config.n_nodes,
+        "degree": config.degree,
+        "f_adv": f_adv,
+        "adversary_mode": mode,
+        "graph_seed": config.graph_seed,
+        "placement_rep": placement_rep,
+        **adv,
+    }
+
+
+def deanon_row(config: SimConfig, blend_hops: int, f_adv: float, mode: str,
+               placement_rep: int, adv: dict, deanon: dict) -> dict:
+    """One row of the deanonymization table: a (placement x blend-path-length) cell.
+
+    ``blend_hops`` comes from the propagation grid, the rest from the adversary placement; the two
+    are crossed here because deanonymization is where propagation paths meet the adversary set.
+    """
+    return {
+        "n_nodes": config.n_nodes,
+        "degree": config.degree,
+        "blend_hops": blend_hops,
+        "f_adv": f_adv,
+        "adversary_mode": mode,
+        "graph_seed": config.graph_seed,
+        "placement_rep": placement_rep,
+        "n_adv": adv["n_adv"],
+        "n_honest": adv["n_honest"],
+        "observed_frac": adv["observed_frac"],
+        **deanon,
+    }

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

@@ -0,0 +1,36 @@
+"""Free-running release-clock mixing delay for a blend relay.
+
+Each relay releases on its own ongoing schedule whose successive intervals are
+``S ~ Uniform{0,1,...,max_blend_delay}`` whole seconds (0 = release now), re-drawn after each
+release. A message arrives at a stationary random phase relative to this clock, so the mixing
+delay it experiences is the renewal-process **residual life** to the next release:
+
+    covering interval  S*  is size-biased:  P(S* = s) proportional to s   (s in 1..M);
+    phase within it is uniform, so residual  R = Uniform(0, S*)  seconds.
+
+Mean residual = (2M+1)/6 seconds (= E[S^2]/(2 E[S]) for S ~ Uniform{0..M}). Returned in ms.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+
+def mix_wait(rng: np.random.Generator, max_blend_delay: int, size: int) -> np.ndarray:
+    """``size`` i.i.d. mixing-delay residuals (ms) for a Uniform{0..max_blend_delay}-sec clock."""
+    m = int(max_blend_delay)
+    if m <= 0 or size <= 0:
+        return np.zeros(max(size, 0), dtype=float)
+    s = np.arange(1, m + 1, dtype=float)
+    probs = s / s.sum()                              # size-biased over positive intervals
+    covering = rng.choice(s, size=size, p=probs)
+    residual_seconds = rng.uniform(0.0, covering)    # uniform phase within the covering interval
+    return residual_seconds * 1000.0                 # -> milliseconds
+
+
+def mean_residual_ms(max_blend_delay: int) -> float:
+    """Analytic mean mixing delay (ms): (2M+1)/6 seconds."""
+    m = int(max_blend_delay)
+    if m <= 0:
+        return 0.0
+    return (2.0 * m + 1.0) / 6.0 * 1000.0

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

@@ -0,0 +1 @@
+"""Plotting: shared style + figure builders for pd."""

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

@@ -0,0 +1,382 @@
+"""Figure builders for pd. Each takes (prop_df, adv_df) and returns a Figure or None.
+
+Propagation full delay (ms) vs peering degree / blend-path length / network size, and adversary
+observation + eclipse fractions vs adversary fraction / degree (with the worst-case envelope) plus
+heatmaps. Slices default to the largest N and a representative setting.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+
+from . import style
+
+_MS = "full propagation delay (ms)"
+
+
+def _prov(*dfs: pd.DataFrame) -> str:
+    ns, seeds = set(), set()
+    for df in dfs:
+        if df is not None and len(df):
+            ns |= set(df["n_nodes"].unique())
+            seeds |= set(df["graph_seed"].unique())
+    return f"pd  |  N={sorted(int(x) for x in ns)}  seeds={len(seeds)}"
+
+
+def _largest_n(df: pd.DataFrame) -> int:
+    return int(sorted(df["n_nodes"].unique())[-1])
+
+
+def delay_vs_degree(prop: pd.DataFrame, adv: pd.DataFrame):
+    if prop is None or not len(prop):
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(prop)
+    mbd = int(sorted(prop["max_blend_delay"].unique())[0])
+    d = prop[(prop.n_nodes == n) & (prop.max_blend_delay == mbd)]
+    fig, ax = plt.subplots()
+    for i, bh in enumerate(sorted(d.blend_hops.unique())):
+        s = d[d.blend_hops == bh].groupby("degree").full_delay_ms_mean.mean().reset_index()
+        ax.plot(s.degree, s.full_delay_ms_mean, "-o", ms=4, color=style.color_for(i),
+                label=f"blend_hops={bh}")
+    ax.set_xlabel("peering degree")
+    ax.set_ylabel(_MS)
+    ax.set_title(f"Blend full delay vs peering degree (N={n:,}, max_blend_delay={mbd}s)")
+    ax.legend()
+    return fig
+
+
+def _median_degree(d: pd.DataFrame) -> int:
+    degs = sorted(d.degree.unique())
+    return int(degs[len(degs) // 2])
+
+
+def delivery_vs_unresponsive(prop: pd.DataFrame, adv: pd.DataFrame):
+    """Message success-delivery-rate vs the unresponsive fraction, one line per blend-path length.
+
+    A message is delivered only if every relay on its (responsiveness-blind) path forwards, so the
+    rate tracks the analytic ``(1-u)^blend_hops`` cascade-survival law (dashed) and is dominated by
+    path length, not degree.
+    """
+    if prop is None or not len(prop) or prop["unresponsive_frac"].nunique() < 2:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(prop)
+    mbd = int(sorted(prop["max_blend_delay"].unique())[0])
+    d = prop[(prop.n_nodes == n) & (prop.max_blend_delay == mbd)]
+    deg = _median_degree(d)
+    d = d[d.degree == deg]
+    fig, ax = plt.subplots()
+    for i, bh in enumerate(sorted(d.blend_hops.unique())):
+        s = d[d.blend_hops == bh].groupby("unresponsive_frac").delivery_rate.mean().reset_index()
+        c = style.color_for(i)
+        ax.plot(s.unresponsive_frac, s.delivery_rate, "-o", ms=4, color=c,
+                label=f"blend_hops={bh}")
+        u = np.linspace(0.0, float(s.unresponsive_frac.max()), 50)
+        ax.plot(u, (1.0 - u) ** bh, "--", lw=0.8, color=c, alpha=0.6)
+    ax.set_xlabel("unresponsive fraction  u")
+    ax.set_ylabel("message delivery rate")
+    ax.set_ylim(0.0, 1.02)
+    ax.set_title(f"Cascade delivery vs unresponsive nodes "
+                 f"(N={n:,}, degree={deg}); dashed = $(1-u)^{{hops}}$")
+    ax.legend()
+    return fig
+
+
+def coverage_vs_unresponsive(prop: pd.DataFrame, adv: pd.DataFrame):
+    """Flood coverage of *delivered* messages vs the unresponsive fraction, one line per degree.
+
+    Unresponsive nodes still receive but do not forward, so they strand pockets of the network; a
+    higher peering degree supplies redundant paths that keep coverage high as ``u`` rises.
+    """
+    if prop is None or not len(prop) or prop["unresponsive_frac"].nunique() < 2:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(prop)
+    mbd = int(sorted(prop["max_blend_delay"].unique())[0])
+    bh = int(sorted(prop["blend_hops"].unique())[0])
+    d = prop[(prop.n_nodes == n) & (prop.max_blend_delay == mbd) & (prop.blend_hops == bh)]
+    fig, ax = plt.subplots()
+    for i, deg in enumerate(sorted(d.degree.unique())):
+        s = d[d.degree == deg].groupby("unresponsive_frac").frac_reached.mean().reset_index()
+        ax.plot(s.unresponsive_frac, s.frac_reached, "-o", ms=4, color=style.color_for(i),
+                label=f"degree={deg}")
+    ax.set_xlabel("unresponsive fraction  u")
+    ax.set_ylabel("flood coverage  (fraction reached | delivered)")
+    ax.set_title(f"Flood coverage vs unresponsive nodes (N={n:,}, blend_hops={bh})")
+    ax.legend()
+    return fig
+
+
+def delay_vs_blendhops(prop: pd.DataFrame, adv: pd.DataFrame):
+    if prop is None or not len(prop) or prop["blend_hops"].nunique() < 2:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(prop)
+    mbd = int(sorted(prop["max_blend_delay"].unique())[0])
+    d = prop[(prop.n_nodes == n) & (prop.max_blend_delay == mbd)]
+    fig, ax = plt.subplots()
+    for i, deg in enumerate(sorted(d.degree.unique())):
+        s = d[d.degree == deg].groupby("blend_hops").full_delay_ms_mean.mean().reset_index()
+        ax.plot(s.blend_hops, s.full_delay_ms_mean, "-o", ms=4, color=style.color_for(i),
+                label=f"degree={deg}")
+    ax.set_xlabel("blend-path hops")
+    ax.set_ylabel(_MS)
+    ax.set_title(f"Blend full delay vs path length (N={n:,}, max_blend_delay={mbd}s)")
+    ax.legend()
+    return fig
+
+
+def delay_vs_N(prop: pd.DataFrame, adv: pd.DataFrame):
+    if prop is None or not len(prop) or prop["n_nodes"].nunique() < 2:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    mbd = int(sorted(prop["max_blend_delay"].unique())[0])
+    bh = int(sorted(prop["blend_hops"].unique())[0])
+    d = prop[(prop.max_blend_delay == mbd) & (prop.blend_hops == bh)]
+    fig, ax = plt.subplots()
+    for i, deg in enumerate(sorted(d.degree.unique())):
+        s = d[d.degree == deg].groupby("n_nodes").full_delay_ms_mean.mean().reset_index()
+        ax.plot(s.n_nodes, s.full_delay_ms_mean, "-o", ms=4, color=style.color_for(i),
+                label=f"degree={deg}")
+    ax.set_xscale("log")
+    ax.set_xlabel("network size N")
+    ax.set_ylabel(_MS)
+    ax.set_title(f"Blend full delay vs network size (blend_hops={bh}, max_blend_delay={mbd}s)")
+    ax.legend()
+    return fig
+
+
+def _adv_vs_fadv(adv: pd.DataFrame, col: str, ylabel: str, title: str):
+    if adv is None or not len(adv):
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(adv)
+    degs = sorted(adv[adv.n_nodes == n].degree.unique())
+    deg = int(degs[len(degs) // 2])
+    d = adv[(adv.n_nodes == n) & (adv.degree == deg)]
+    fig, ax = plt.subplots()
+    for i, mode in enumerate(sorted(d.adversary_mode.unique())):
+        s = d[d.adversary_mode == mode].groupby("f_adv")[col].mean().reset_index()
+        ax.plot(s.f_adv, s[col], "-o", ms=4, color=style.color_for(i), label=mode)
+    ax.set_xlabel("adversary fraction  f_adv")
+    ax.set_ylabel(ylabel)
+    ax.set_title(f"{title} (N={n:,}, degree={deg})")
+    ax.legend()
+    return fig
+
+
+def observed_vs_fadv(prop, adv):
+    return _adv_vs_fadv(adv, "observed_frac", "honest observed (fraction)",
+                        "Adversary observation vs f_adv")
+
+
+def eclipse_vs_fadv(prop, adv):
+    return _adv_vs_fadv(adv, "eclipsed_frac", "honest eclipsed (fraction)",
+                        "Honest eclipse vs f_adv")
+
+
+def _adv_vs_degree(adv: pd.DataFrame, col: str, ylabel: str, title: str):
+    if adv is None or not len(adv) or adv["degree"].nunique() < 2:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(adv)
+    d = adv[adv.n_nodes == n]
+    favs = sorted(d.f_adv.unique())
+    f = favs[len(favs) // 2]
+    d = d[d.f_adv == f]
+    fig, ax = plt.subplots()
+    for i, mode in enumerate(sorted(d.adversary_mode.unique())):
+        s = d[d.adversary_mode == mode].groupby("degree")[col].mean().reset_index()
+        ax.plot(s.degree, s[col], "-o", ms=4, color=style.color_for(i), label=mode)
+    ax.set_xlabel("peering degree")
+    ax.set_ylabel(ylabel)
+    ax.set_title(f"{title} (N={n:,}, f_adv={f})")
+    ax.legend()
+    return fig
+
+
+def observed_vs_degree(prop, adv):
+    return _adv_vs_degree(adv, "observed_frac", "honest observed (fraction)",
+                          "Adversary observation vs degree")
+
+
+def eclipse_vs_degree(prop, adv):
+    return _adv_vs_degree(adv, "eclipsed_frac", "honest eclipsed (fraction)",
+                          "Honest eclipse vs degree")
+
+
+def _heatmap(adv: pd.DataFrame, col: str, title: str):
+    if adv is None or not len(adv) or adv["degree"].nunique() < 2 or adv["f_adv"].nunique() < 2:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(adv)
+    d = adv[(adv.n_nodes == n) & (adv.adversary_mode == "random")]
+    if not len(d):
+        d = adv[adv.n_nodes == n]
+    piv = d.groupby(["degree", "f_adv"])[col].mean().unstack("f_adv")
+    fig, ax = plt.subplots()
+    im = ax.imshow(piv.values, origin="lower", aspect="auto", cmap=style.SEQUENTIAL_CMAP)
+    ax.set_xticks(range(len(piv.columns)), [f"{c:g}" for c in piv.columns])
+    ax.set_yticks(range(len(piv.index)), [str(int(i)) for i in piv.index])
+    ax.set_xlabel("adversary fraction  f_adv")
+    ax.set_ylabel("peering degree")
+    ax.set_title(f"{title} (N={n:,})")
+    fig.colorbar(im, ax=ax, shrink=0.85)
+    return fig
+
+
+def heatmap_observed(prop, adv):
+    return _heatmap(adv, "observed_frac", "Honest observed fraction")
+
+
+def heatmap_eclipse(prop, adv):
+    return _heatmap(adv, "eclipsed_frac", "Honest eclipsed fraction")
+
+
+# --- deanonymization: propagation paths x the adversary set (3-arg builders) -------------------
+
+def _pos(s: pd.Series) -> pd.Series:
+    """Blank out non-positive rates so they vanish on a log axis instead of erroring."""
+    return s.where(s > 0)
+
+
+def deanon_vs_blendhops(prop, adv, deanon):
+    """P(whole blend path adversarial) vs path length, one line per f_adv (log-y).
+
+    Relays are drawn blind to who is adversarial, so the rate is the exact hypergeometric ~
+    ``f_adv**blend_hops`` (dashed) -- lengthening the blend path is the dominant defence, and it is
+    independent of peering degree.
+    """
+    if deanon is None or not len(deanon) or deanon["blend_hops"].nunique() < 2:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(deanon)
+    d = deanon[(deanon.n_nodes == n) & (deanon.adversary_mode == "random") & (deanon.f_adv > 0)]
+    if not len(d):
+        return None
+    deg = _median_degree(d)
+    d = d[d.degree == deg]
+    fig, ax = plt.subplots()
+    for i, f in enumerate(sorted(d.f_adv.unique())):
+        s = d[d.f_adv == f].groupby("blend_hops").deanon_rate.mean().reset_index()
+        c = style.color_for(i)
+        ax.plot(s.blend_hops, _pos(s.deanon_rate), "-o", ms=4, color=c, label=f"f_adv={f:g}")
+        ax.plot(s.blend_hops, f ** s.blend_hops, "--", lw=0.8, color=c, alpha=0.6)
+    ax.set_yscale("log")
+    ax.set_xlabel("blend-path hops")
+    ax.set_ylabel("deanonymization rate\nP(whole path adversarial)")
+    ax.set_title(f"Deanonymization vs path length (N={n:,}, degree={deg}); "
+                 f"dashed = $f_{{adv}}^{{hops}}$")
+    ax.legend()
+    return fig
+
+
+def full_deanon_vs_blendhops(prop, adv, deanon):
+    """P(whole path adversarial AND sender peered with an adversary) vs path length, per f_adv.
+
+    Dashed = the closed form ``f_adv**hops * (1-(1-f_adv)**degree)`` at the plotted degree.
+    """
+    if deanon is None or not len(deanon) or deanon["blend_hops"].nunique() < 2:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(deanon)
+    d = deanon[(deanon.n_nodes == n) & (deanon.adversary_mode == "random") & (deanon.f_adv > 0)]
+    if not len(d):
+        return None
+    deg = _median_degree(d)
+    d = d[d.degree == deg]
+    fig, ax = plt.subplots()
+    for i, f in enumerate(sorted(d.f_adv.unique())):
+        s = d[d.f_adv == f].groupby("blend_hops").full_deanon_rate.mean().reset_index()
+        c = style.color_for(i)
+        ax.plot(s.blend_hops, _pos(s.full_deanon_rate), "-o", ms=4, color=c, label=f"f_adv={f:g}")
+        ax.plot(s.blend_hops, f ** s.blend_hops * (1.0 - (1.0 - f) ** deg), "--",
+                lw=0.8, color=c, alpha=0.6)
+    ax.set_yscale("log")
+    ax.set_xlabel("blend-path hops")
+    ax.set_ylabel("full-deanonymization rate\nP(path adversarial & sender exposed)")
+    ax.set_title(f"Full deanonymization vs path length (N={n:,}, degree={deg})")
+    ax.legend()
+    return fig
+
+
+def full_deanon_vs_fadv(prop, adv, deanon):
+    """Full deanonymization vs adversary fraction, random vs worst-case placement (log-y).
+
+    The whole-path-adversarial rate is placement-independent (thin grey ceiling); full
+    deanonymization adds the sender-peer factor ``observed_frac``, which the worst-case-coverage
+    adversary maximizes -- so the placements fan out below that ceiling.
+    """
+    if deanon is None or not len(deanon) or deanon["f_adv"].nunique() < 2:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(deanon)
+    bhs = sorted(deanon.blend_hops.unique())
+    bh = int(bhs[len(bhs) // 2])
+    d = deanon[(deanon.n_nodes == n) & (deanon.blend_hops == bh) & (deanon.f_adv > 0)]
+    if not len(d):
+        return None
+    deg = _median_degree(d)
+    d = d[d.degree == deg]
+    fig, ax = plt.subplots()
+    for i, mode in enumerate(sorted(d.adversary_mode.unique())):
+        s = d[d.adversary_mode == mode].groupby("f_adv").full_deanon_rate.mean().reset_index()
+        ax.plot(s.f_adv, _pos(s.full_deanon_rate), "-o", ms=4, color=style.color_for(i),
+                label=f"full ({mode})")
+    ceil = d[d.adversary_mode == "random"].groupby("f_adv").deanon_rate.mean().reset_index()
+    if len(ceil):
+        ax.plot(ceil.f_adv, _pos(ceil.deanon_rate), ":", lw=1.0, color="0.5",
+                label="whole path (any placement)")
+    ax.set_yscale("log")
+    ax.set_xlabel("adversary fraction  f_adv")
+    ax.set_ylabel("deanonymization rate")
+    ax.set_title(f"Full deanonymization vs f_adv (N={n:,}, degree={deg}, blend_hops={bh})")
+    ax.legend()
+    return fig
+
+
+def full_deanon_vs_degree(prop, adv, deanon):
+    """Full deanonymization vs peering degree, per f_adv (fixed path length, log-y).
+
+    ``deanon_rate`` (dashed) is degree-flat -- relays are uniform -- but the sender-peer factor
+    ``1-(1-f_adv)**degree`` rises with degree, so *full* deanonymization worsens as degree grows
+    even though a higher degree speeds propagation. The peering-degree tension, in one plot.
+    """
+    if deanon is None or not len(deanon) or deanon["degree"].nunique() < 2:
+        return None
+    import matplotlib.pyplot as plt
+    style.apply_style()
+    n = _largest_n(deanon)
+    bhs = sorted(deanon.blend_hops.unique())
+    bh = int(bhs[len(bhs) // 2])
+    d = deanon[(deanon.n_nodes == n) & (deanon.blend_hops == bh)
+               & (deanon.adversary_mode == "random") & (deanon.f_adv > 0)]
+    if not len(d):
+        return None
+    fig, ax = plt.subplots()
+    for i, f in enumerate(sorted(d.f_adv.unique())):
+        c = style.color_for(i)
+        s = d[d.f_adv == f].groupby("degree").agg(
+            full=("full_deanon_rate", "mean"), whole=("deanon_rate", "mean")).reset_index()
+        ax.plot(s.degree, _pos(s.full), "-o", ms=4, color=c, label=f"full, f_adv={f:g}")
+        ax.plot(s.degree, _pos(s.whole), "--", lw=0.8, color=c, alpha=0.6)
+    ax.set_yscale("log")
+    ax.set_xlabel("peering degree")
+    ax.set_ylabel("deanonymization rate")
+    ax.set_title(f"Full deanonymization vs degree (N={n:,}, blend_hops={bh}); "
+                 f"dashed = whole-path (degree-flat)")
+    ax.legend()
+    return fig

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

@@ -0,0 +1,71 @@
+"""Render pd figures from the propagation + adversary parquets (pd-figures)."""
+
+from __future__ import annotations
+
+import argparse
+from pathlib import Path
+
+import pandas as pd
+
+from . import figures, style
+
+# (prop, adv) builders.
+_BUILDERS = [
+    ("01_delay_vs_degree", figures.delay_vs_degree),
+    ("02_delay_vs_blendhops", figures.delay_vs_blendhops),
+    ("03_delay_vs_N", figures.delay_vs_N),
+    ("04_observed_vs_fadv", figures.observed_vs_fadv),
+    ("05_eclipse_vs_fadv", figures.eclipse_vs_fadv),
+    ("06_observed_vs_degree", figures.observed_vs_degree),
+    ("07_eclipse_vs_degree", figures.eclipse_vs_degree),
+    ("08_heatmap_observed", figures.heatmap_observed),
+    ("09_heatmap_eclipse", figures.heatmap_eclipse),
+    ("10_delivery_vs_unresponsive", figures.delivery_vs_unresponsive),
+    ("11_coverage_vs_unresponsive", figures.coverage_vs_unresponsive),
+]
+
+# (prop, adv, deanon) builders — deanonymization crosses propagation paths with the adversary set.
+_DEANON_BUILDERS = [
+    ("12_deanon_vs_blendhops", figures.deanon_vs_blendhops),
+    ("13_full_deanon_vs_blendhops", figures.full_deanon_vs_blendhops),
+    ("14_full_deanon_vs_fadv", figures.full_deanon_vs_fadv),
+    ("15_full_deanon_vs_degree", figures.full_deanon_vs_degree),
+]
+
+
+def render(prop_df: pd.DataFrame, adv_df: pd.DataFrame, deanon_df: pd.DataFrame,
+           out_dir: Path) -> 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])
+    written: list[Path] = []
+    for name, fn, fn_args in jobs:
+        try:
+            fig = fn(*fn_args)
+        except Exception as e:  # noqa: BLE001 — a bad slice shouldn't kill the whole render
+            print(f"skip {name}: {e}")
+            continue
+        if fig is not None:
+            written += style.save(fig, out_dir / name, provenance=prov)
+    return written
+
+
+def main(argv: list[str] | None = None) -> int:
+    ap = argparse.ArgumentParser(description="Render pd figures from a run directory")
+    ap.add_argument("--run", required=True, help="runs/<ts>_<label>/ containing the parquets")
+    ap.add_argument("--out", default=None, help="output dir (default <run>/figures)")
+    args = ap.parse_args(argv)
+    run = Path(args.run)
+    prop_df = pd.read_parquet(run / "propagation.parquet")
+    adv_df = pd.read_parquet(run / "adversary.parquet")
+    dz_path = run / "deanon.parquet"
+    deanon_df = pd.read_parquet(dz_path) if dz_path.exists() else pd.DataFrame()
+    out = Path(args.out) if args.out else run / "figures"
+    written = render(prop_df, adv_df, deanon_df, out)
+    print(f"wrote {len(written)} figures -> {out}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

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

@@ -0,0 +1,90 @@
+"""Shared academic matplotlib theme, palette, and helpers (copied from tsi-sim-mc).
+
+Palette: Okabe-Ito — the standard colorblind-safe qualitative set — for categorical
+series; ``cividis`` (perceptually uniform, CVD-safe) for heatmaps. Figures are saved as
+300-dpi PNG.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Sequence
+from pathlib import Path
+
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+
+# Okabe-Ito colorblind-safe qualitative palette
+OKABE_ITO = [
+    "#0072B2",  # blue
+    "#D55E00",  # vermillion
+    "#009E73",  # bluish green
+    "#CC79A7",  # reddish purple
+    "#E69F00",  # orange
+    "#56B4E9",  # sky blue
+    "#F0E442",  # yellow
+    "#000000",  # black
+]
+SEQUENTIAL_CMAP = "cividis"
+DIVERGING_CMAP = "RdBu_r"
+
+
+def apply_style() -> None:
+    """Install the shared rcParams theme (idempotent)."""
+    mpl.rcParams.update(
+        {
+            "figure.dpi": 150,
+            "savefig.dpi": 300,
+            "savefig.bbox": "tight",
+            "figure.figsize": (6.4, 4.0),
+            "font.size": 10,
+            "font.family": "sans-serif",
+            "axes.titlesize": 11,
+            "axes.labelsize": 10,
+            "axes.spines.top": False,
+            "axes.spines.right": False,
+            "axes.grid": True,
+            "grid.alpha": 0.25,
+            "grid.linewidth": 0.6,
+            "legend.frameon": False,
+            "legend.fontsize": 8.5,
+            "lines.linewidth": 1.8,
+            "axes.prop_cycle": mpl.cycler(color=OKABE_ITO),
+            "mathtext.default": "regular",
+        }
+    )
+
+
+def color_for(index: int) -> str:
+    return OKABE_ITO[index % len(OKABE_ITO)]
+
+
+def band_plot(
+    ax: plt.Axes,
+    x: Sequence[float],
+    series: np.ndarray,
+    *,
+    color: str,
+    label: str | None = None,
+    percentiles: tuple[float, float] = (10, 90),
+) -> None:
+    """Plot the mean of ``series`` (shape ``(n_replicates, len(x))``) with a percentile band."""
+    x = np.asarray(x, dtype=float)
+    mean = np.nanmean(series, axis=0)
+    lo = np.nanpercentile(series, percentiles[0], axis=0)
+    hi = np.nanpercentile(series, percentiles[1], axis=0)
+    ax.plot(x, mean, color=color, label=label)
+    ax.fill_between(x, lo, hi, color=color, alpha=0.18, linewidth=0)
+
+
+def save(fig: plt.Figure, out_stem: str | Path, provenance: str | None = None) -> list[Path]:
+    """Save ``fig`` as a 300-dpi PNG. ``out_stem`` has no suffix. Returns written paths."""
+    out_stem = Path(out_stem)
+    out_stem.parent.mkdir(parents=True, exist_ok=True)
+    if provenance:
+        fig.subplots_adjust(bottom=0.16)
+        fig.text(0.99, 0.005, provenance, fontsize=6, alpha=0.5, va="bottom", ha="right")
+    p = out_stem.with_suffix(".png")
+    fig.savefig(p)
+    plt.close(fig)
+    return [p]

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

@@ -0,0 +1,163 @@
+"""Blend-cascade propagation: per-round delivery success and full delay (ms), aggregated.
+
+One round: a (responsive) sender injects a message that must traverse a **blend path** of
+``blend_hops`` relays chosen *blind to responsiveness* (the path is drawn from the consensus node
+list, not from who happens to be up) -> the last relay floods the whole network. Transport legs are
+directed shortest paths with edge weights ``base + Exp(jitter) + p(relaying node)`` (ms); each relay
+adds a free-running-clock mixing wait (``mixclock.mix_wait``). Full delay = sum of transport legs +
+sum of mixing waits + the final broadcast delay. Only the ``blend_hops`` relays mix; the final flood
+is plain.
+
+**Unresponsive nodes** (a configurable ``unresponsive_frac`` of the population) relay nothing: their
+outgoing edges are removed (weight -> inf), so nothing routes *through* them, though they can still
+*receive* (be reached as a leaf). Two consequences, both measured:
+
+* **message success-delivery-rate** (``delivery_rate``) — the fraction of rounds whose message is
+  delivered through the whole cascade to a **responsive** final relay that then floods. If any relay
+  on the drawn path is unresponsive the message dies there, so the rate falls ~
+  ``(1-f)^blend_hops``; a shorter path or a higher degree that keeps legs routable raises it.
+* **coverage** (``frac_reached``, ``coverN_ms``) of a delivered flood, which degrades as routing
+  holes strand pockets of the network.
+
+Delay statistics are conditioned on delivery; ``delivery_rate`` is over all rounds.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+from scipy.sparse.csgraph import dijkstra
+
+from .config import SimConfig
+from .graph import Graph
+from .memguard import check_alloc
+from .mixclock import mix_wait
+
+
+def assign_responsive(n: int, unresponsive_frac: float, rng: np.random.Generator) -> np.ndarray:
+    """Boolean mask (True == responsive/relaying). A random ``unresponsive_frac`` are set False."""
+    responsive = np.ones(n, dtype=bool)
+    n_unresp = int(round(unresponsive_frac * n))
+    if n_unresp > 0:
+        responsive[rng.choice(n, size=n_unresp, replace=False)] = False
+    return responsive
+
+
+def blend_round(graph: Graph, sender: int, relays: np.ndarray, jitter_mean_ms: float,
+                max_blend_delay: int, rng: np.random.Generator,
+                coverage_pcts: tuple[float, ...], responsive: np.ndarray | None = None) -> dict:
+    """One Blend cascade. ``delivered`` is True iff every relay forwards and the final relay (which
+    must be responsive) floods; delay fields are NaN on a dropped message."""
+    data = (graph.base
+            + rng.exponential(jitter_mean_ms, size=graph.base.shape[0])
+            + graph.p[graph.src])
+    if responsive is not None:
+        data[~responsive[graph.src]] = np.inf   # unresponsive nodes relay nothing (no out-edges)
+    csr = graph.weighted_csr(data)
+    k = int(relays.shape[0])
+    sources = np.empty(k + 1, dtype=np.int64)
+    sources[0] = sender
+    sources[1:] = relays
+    dist = dijkstra(csr, directed=True, indices=sources)   # (k+1, n)
+
+    # legs s->r1->...->rk. A leg from an unresponsive relay is inf (no outgoing edges), so legs_ok
+    # already encodes "every intermediate relay forwarded"; the final relay is checked below.
+    legs = 0.0
+    legs_ok = True
+    for i in range(k):
+        d = dist[i, relays[i]]
+        if not np.isfinite(d):
+            legs_ok = False
+            break
+        legs += float(d)
+    mix_total = float(mix_wait(rng, max_blend_delay, k).sum())
+
+    final_relay = int(sources[k])
+    final_ok = responsive is None or bool(responsive[final_relay])
+
+    flood = dist[k]
+    finite = np.isfinite(flood)
+    delivered = legs_ok and final_ok and bool(finite.any())
+    if delivered:
+        reached = flood[finite]
+        broadcast = float(reached.max())
+        covers = [float(np.percentile(reached, pc)) for pc in coverage_pcts]
+        frac_reached = float(finite.mean())
+        path = legs + mix_total
+        full = path + broadcast
+    else:
+        broadcast = float("nan")
+        covers = [float("nan")] * len(coverage_pcts)
+        frac_reached = 0.0
+        path = float("nan")
+        full = float("nan")
+    return {"full": full, "path": path, "broadcast": broadcast, "covers": covers,
+            "frac_reached": frac_reached, "delivered": delivered}
+
+
+def propagation_metrics(graph: Graph, blend_hops: int, max_blend_delay: int,
+                        unresponsive_frac: float, responsive: np.ndarray,
+                        config: SimConfig, rng: np.random.Generator) -> dict:
+    """Aggregate the Blend cascade over ``config.n_rounds`` rounds.
+
+    Each round a responsive node injects a message and the ``blend_hops`` relays are drawn from the
+    whole node list (blind to responsiveness). ``delivery_rate`` is the fraction that complete the
+    cascade; delay/coverage statistics are conditioned on those delivered rounds.
+    """
+    n = graph.n
+    check_alloc(int((blend_hops + 1) * n * 8), "sampled (blend_hops+1) x N distance matrix",
+                f"N={n}, blend_hops={blend_hops}")
+    pcts = tuple(config.coverage_pcts)
+    resp_ids = np.where(responsive)[0]
+
+    def _empty() -> dict:
+        out = {name: float("nan") for name in
+               ("full_delay_ms_mean", "full_delay_ms_p50", "full_delay_ms_p90",
+                "full_delay_ms_p99", "path_delay_ms_mean", "broadcast_delay_ms_mean")}
+        out["frac_reached"] = 0.0
+        out["delivery_rate"] = 0.0
+        for pc in pcts:
+            out[f"cover{int(pc)}_ms"] = float("nan")
+        return out
+
+    if resp_ids.shape[0] < 1 or n < blend_hops + 1:
+        return _empty()   # no responsive sender, or too few nodes to draw a distinct path
+
+    fulls, paths, bcasts, fracs = [], [], [], []
+    covers = [[] for _ in pcts]
+    delivered = 0
+    for _ in range(config.n_rounds):
+        sender = int(rng.choice(resp_ids))
+        relays = rng.choice(n - 1, size=blend_hops, replace=False)
+        relays[relays >= sender] += 1          # blend_hops distinct nodes, all != sender
+        r = blend_round(graph, sender, relays, config.transport_jitter_mean_ms,
+                        max_blend_delay, rng, pcts, responsive)
+        if not r["delivered"]:
+            continue
+        delivered += 1
+        fulls.append(r["full"])
+        paths.append(r["path"])
+        bcasts.append(r["broadcast"])
+        fracs.append(r["frac_reached"])
+        for j, c in enumerate(r["covers"]):
+            covers[j].append(c)
+
+    delivery_rate = delivered / config.n_rounds
+    if delivered == 0:
+        out = _empty()
+        out["delivery_rate"] = delivery_rate
+        return out
+
+    fulls = np.asarray(fulls)
+    out = {
+        "full_delay_ms_mean": float(np.mean(fulls)),
+        "full_delay_ms_p50": float(np.percentile(fulls, 50)),
+        "full_delay_ms_p90": float(np.percentile(fulls, 90)),
+        "full_delay_ms_p99": float(np.percentile(fulls, 99)),
+        "path_delay_ms_mean": float(np.mean(paths)),
+        "broadcast_delay_ms_mean": float(np.mean(bcasts)),
+        "frac_reached": float(np.mean(fracs)),
+        "delivery_rate": delivery_rate,
+    }
+    for pc, col in zip(pcts, covers, strict=True):
+        out[f"cover{int(pc)}_ms"] = float(np.mean(col))
+    return out

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

@@ -0,0 +1,68 @@
+"""Deterministic, order-independent RNG derivation (blake2b -> SeedSequence).
+
+Three independent seed streams so that graph, propagation and adversary randomness are
+reproducible and separable:
+- ``graph_seedseq`` depends ONLY on the topology fields, so the peer graph + per-node
+  processing lags are a pure function of the global seed (the "consensus topology" property)
+  and are IDENTICAL across every ``f_adv``/``adversary_mode`` cell that shares the topology;
+- ``round_seedseq`` seeds the per-round sender/relay/jitter/mix draws;
+- ``placement_seedseq`` seeds the adversary placement, independent of the graph/rounds.
+"""
+
+from __future__ import annotations
+
+import hashlib
+
+import numpy as np
+
+from .config import SimConfig
+
+
+def _digest(*parts: object) -> int:
+    payload = repr(parts).encode()
+    return int.from_bytes(hashlib.blake2b(payload, digest_size=16).digest(), "big")
+
+
+def seedseq_for(config: SimConfig) -> np.random.SeedSequence:
+    """Full-key root SeedSequence for this exact config (parity with the sibling sims)."""
+    return np.random.SeedSequence(_digest(config.root_seed, config.key()))
+
+
+def rng_for(config: SimConfig) -> np.random.Generator:
+    return np.random.default_rng(seedseq_for(config))
+
+
+def graph_seedseq(config: SimConfig) -> np.random.SeedSequence:
+    """Topology-only seed: peer graph + processing lags depend on these fields alone."""
+    return np.random.SeedSequence(_digest(
+        config.root_seed, "graph", config.n_nodes, config.degree, config.graph_seed,
+        config.link_latency_dist, config.link_latency_mean_ms,
+        config.processing_lags_ms, config.processing_lag_probs,
+    ))
+
+
+def responsive_seedseq(config: SimConfig, unresponsive_frac: float) -> np.random.SeedSequence:
+    """Which nodes are responsive: fixed per (topology, unresponsive_frac), not per round."""
+    return np.random.SeedSequence(_digest(
+        config.root_seed, "responsive", config.n_nodes, config.degree, config.graph_seed,
+        unresponsive_frac,
+    ))
+
+
+def round_seedseq(config: SimConfig, blend_hops: int, max_blend_delay: int,
+                  unresponsive_frac: float) -> np.random.SeedSequence:
+    """Per-cell propagation seed (senders/relays/jitter/mix over n_rounds)."""
+    return np.random.SeedSequence(_digest(
+        config.root_seed, "rounds", config.n_nodes, config.degree, config.graph_seed,
+        blend_hops, max_blend_delay, unresponsive_frac, config.n_rounds,
+        config.transport_jitter_mean_ms,
+    ))
+
+
+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."""
+    return np.random.SeedSequence(_digest(
+        config.root_seed, "place", config.n_nodes, config.degree, config.graph_seed,
+        f_adv, mode, placement_rep,
+    ))

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

@@ -0,0 +1,93 @@
+"""Run a pd sweep from a YAML config -> propagation.parquet + adversary.parquet + figures.
+
+Each topology ``(n_nodes, degree, graph_seed)`` is one embarrassingly-parallel work item; the
+engine builds it once and measures the propagation and adversary sub-grids on it.
+"""
+
+from __future__ import annotations
+
+import argparse
+import datetime as _dt
+from pathlib import Path
+
+import pandas as pd
+import yaml
+from joblib import Parallel, delayed
+from tqdm import tqdm
+
+from .config import SimConfig, SweepConfig
+from .engine import run_graph_cell
+
+HERE = Path(__file__).resolve().parents[2]   # pd/
+
+
+def load_sweep_yaml(path: str | Path) -> SweepConfig:
+    with open(path) as f:
+        return SweepConfig.from_dict(yaml.safe_load(f) or {})
+
+
+def new_run_dir(outdir: Path, label: str) -> Path:
+    ts = _dt.datetime.now().strftime("%Y-%m-%d_%H%M%S")
+    run_dir = outdir / f"{ts}_{label}"
+    suffix = 2
+    while run_dir.exists():
+        run_dir = outdir / f"{ts}_{label}_{suffix}"
+        suffix += 1
+    run_dir.mkdir(parents=True)
+    return run_dir
+
+
+def _cell_worker(base: SimConfig, prop_grid, unresponsive_fracs, adv_grid):
+    return run_graph_cell(base, prop_grid, unresponsive_fracs, adv_grid)
+
+
+def run_sweep(sweep: SweepConfig,
+              n_jobs: int = -1) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]:
+    cells = sweep.graph_cells()
+    prop_grid = sweep.prop_grid()
+    unresponsive_fracs = list(sweep.unresponsive_frac)
+    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, adv_grid)
+        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)
+
+
+def persist(prop_df: pd.DataFrame, adv_df: pd.DataFrame, deanon_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)
+
+
+def main(argv: list[str] | None = None) -> int:
+    ap = argparse.ArgumentParser(description="Run a pd peering-degree sweep")
+    ap.add_argument("--config", required=True, help="sweep YAML (see configs/)")
+    ap.add_argument("--outdir", default=str(HERE / "runs"))
+    ap.add_argument("--label", default=None)
+    ap.add_argument("--n-jobs", type=int, default=-1)
+    ap.add_argument("--no-figures", action="store_true")
+    args = ap.parse_args(argv)
+
+    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)
+    print(f"wrote {len(prop_df)} propagation + {len(adv_df)} adversary + "
+          f"{len(deanon_df)} deanon 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")
+        print(f"wrote {len(figs)} figures -> {run_dir / 'figures'}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

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

@@ -0,0 +1,119 @@
+"""Analytic sanity checks for pd (pd-verify): closed forms + graph invariants.
+
+Random-placement observation/eclipse have exact closed forms on a random d-regular graph:
+    observed_frac ~ 1 - (1 - f_adv)^degree      (an honest node has >=1 adversarial peer)
+    eclipsed_frac ~ f_adv^degree                (all `degree` peers adversarial)
+and the propagation full delay (mixing off) drops as the peering degree rises. Each check prints
+PASS/FAIL; a non-zero exit signals failure.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .adversary import adversary_metrics, deanon_metrics, place_adversary
+from .config import SimConfig
+from .graph import build_graph
+from .rng import placement_seedseq
+
+
+def _check(name: str, ok: bool, detail: str) -> bool:
+    print(f"[{'PASS' if ok else 'FAIL'}] {name}: {detail}")
+    return ok
+
+
+def main(argv: list[str] | None = None) -> int:
+    ok = True
+
+    # 1. graph exactly d-regular + symmetric
+    cfg = SimConfig(n_nodes=2000, degree=8, f_adv=0.0)
+    g = build_graph(cfg)
+    deg = np.diff(g.indptr)
+    regular = bool(np.all(deg == g.degree))
+    csr = g.weighted_csr(np.ones_like(g.base))
+    symmetric = int((csr != csr.T).nnz) == 0
+    ok &= _check("graph d-regular", regular, f"all degrees == {g.degree}: {regular}")
+    ok &= _check("graph symmetric", symmetric, f"adj == adj.T: {symmetric}")
+
+    # 2. random observation / eclipse vs closed form (avg over placements + a couple of seeds)
+    for degree in (4, 8):
+        for f in (0.1, 0.3):
+            obs, ecl = [], []
+            for seed in range(4):
+                gg = build_graph(SimConfig(n_nodes=4000, degree=degree, graph_seed=seed))
+                pc = SimConfig(n_nodes=4000, degree=degree, graph_seed=seed, f_adv=f)
+                rng = np.random.default_rng(placement_seedseq(pc, f, "random", 0))
+                m = adversary_metrics(gg, place_adversary(gg, f, "random", rng, 100_000))
+                obs.append(m["observed_frac"])
+                ecl.append(m["eclipsed_frac"])
+            obs_m, ecl_m = float(np.mean(obs)), float(np.mean(ecl))
+            obs_th = 1 - (1 - f) ** degree
+            ecl_th = f ** degree
+            ok &= _check(f"observed_frac d={degree} f={f}", abs(obs_m - obs_th) < 0.02,
+                         f"sim {obs_m:.3f} vs theory {obs_th:.3f}")
+            ok &= _check(f"eclipsed_frac d={degree} f={f}",
+                         abs(ecl_m - ecl_th) < max(0.01, 0.2 * ecl_th),
+                         f"sim {ecl_m:.4f} vs theory {ecl_th:.4f}")
+
+    # 3. propagation full delay decreases as degree rises (mixing off, small jitter)
+    from .engine import run_trajectory
+    delays = {}
+    for degree in (4, 16):
+        c = SimConfig(n_nodes=4000, degree=degree, blend_hops=3, max_blend_delay=0,
+                      transport_jitter_mean_ms=0.0, n_rounds=40)
+        delays[degree] = run_trajectory(c)["propagation"]["full_delay_ms_mean"]
+    ok &= _check("delay decreases with degree", delays[16] < delays[4],
+                 f"d=4 {delays[4]:.0f} ms > d=16 {delays[16]:.0f} ms")
+
+    # 4. message delivery-rate ~ (1 - u)^blend_hops (relays drawn blind to responsiveness; a high
+    #    degree keeps legs routable so the only loss is a relay landing on an unresponsive node)
+    for u in (0.1, 0.3):
+        rates = []
+        for seed in range(3):
+            c = SimConfig(n_nodes=5000, degree=16, blend_hops=3, max_blend_delay=0,
+                          transport_jitter_mean_ms=0.0, unresponsive_frac=u,
+                          n_rounds=400, graph_seed=seed)
+            rates.append(run_trajectory(c)["propagation"]["delivery_rate"])
+        rate_m = float(np.mean(rates))
+        rate_th = (1 - u) ** 3
+        ok &= _check(f"delivery_rate u={u}", abs(rate_m - rate_th) < 0.03,
+                     f"sim {rate_m:.3f} vs theory {rate_th:.3f}")
+
+    # 5. deanonymization: the exact closed forms reproduce a direct Monte-Carlo of the SAME draw
+    #    (honest sender + blend_hops relays picked blind to who is adversarial). deanon_rate ~
+    #    f_adv^blend_hops (whole cascade adversarial); full adds the sender-has-an-adversary-peer
+    #    factor. degree lifts the full rate (more peers -> sender more exposed), not deanon_rate.
+    for degree, f, k in [(8, 0.33, 2), (16, 0.33, 2), (8, 0.33, 3), (16, 0.2, 3)]:
+        cfg = SimConfig(n_nodes=3000, degree=degree, graph_seed=2, f_adv=f, blend_hops=k)
+        g = build_graph(cfg)
+        prng = np.random.default_rng(placement_seedseq(cfg, f, "random", 0))
+        mask = place_adversary(g, f, "random", prng, cfg.worstcase_max_n)
+        adv = adversary_metrics(g, mask)
+        dz = deanon_metrics(g.n, adv["n_adv"], adv["observed_frac"], k)
+        counts = np.add.reduceat(mask[g.indices].astype(np.int32), g.indptr[:-1])
+        observed_node = counts >= 1
+        honest = np.where(~mask)[0]
+        n = g.n
+        srng = np.random.default_rng(777 + degree + k)
+        trials, d_hit, fd_hit = 60_000, 0, 0
+        for _ in range(trials):
+            s = int(srng.choice(honest))
+            r = srng.choice(n - 1, size=k, replace=False)
+            r[r >= s] += 1                       # blend_hops distinct nodes, all != sender
+            if mask[r].all():
+                d_hit += 1
+                fd_hit += int(observed_node[s])
+        d_emp, fd_emp = d_hit / trials, fd_hit / trials
+        ok &= _check(f"deanon_rate d={degree} f={f} k={k}",
+                     abs(dz["deanon_rate"] - d_emp) < max(0.006, 0.12 * d_emp),
+                     f"closed {dz['deanon_rate']:.4f} vs MC {d_emp:.4f} (~f^k={f ** k:.4f})")
+        ok &= _check(f"full_deanon d={degree} f={f} k={k}",
+                     abs(dz["full_deanon_rate"] - fd_emp) < max(0.006, 0.15 * fd_emp),
+                     f"closed {dz['full_deanon_rate']:.4f} vs MC {fd_emp:.4f}")
+
+    print("OK" if ok else "FAILURES PRESENT")
+    return 0 if ok else 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

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

@@ -0,0 +1,62 @@
+from itertools import combinations
+
+import numpy as np
+
+from pd.adversary import _greedy_coverage, adversary_metrics, place_adversary
+from pd.config import SimConfig
+from pd.graph import Graph, build_graph
+
+
+def _cycle4():
+    indptr = np.array([0, 2, 4, 6, 8], dtype=np.int64)
+    indices = np.array([1, 3, 0, 2, 1, 3, 0, 2], dtype=np.int64)
+    return Graph(n=4, degree=2, indptr=indptr, indices=indices,
+                 base=np.ones(8), src=np.array([0, 0, 1, 1, 2, 2, 3, 3]), p=np.zeros(4))
+
+
+def test_coverage_eclipse_hand_checked():
+    g = _cycle4()  # 0-1-2-3-0
+    m = adversary_metrics(g, np.array([False, True, False, True]))  # adv {1,3}
+    assert m["observed_count"] == 2 and m["eclipsed_count"] == 2   # honest {0,2} fully surrounded
+    m = adversary_metrics(g, np.array([False, True, False, False]))  # adv {1}
+    assert m["observed_count"] == 2 and m["eclipsed_count"] == 0    # {0,2} observed, none eclipsed
+
+
+def test_random_closed_form():
+    g = build_graph(SimConfig(n_nodes=5000, degree=6, graph_seed=0))
+    rng = np.random.default_rng(0)
+    def _obs():
+        return adversary_metrics(g, place_adversary(g, 0.2, "random", rng, 10**9))["observed_frac"]
+    obs = np.mean([_obs() for _ in range(5)])
+    assert abs(obs - (1 - 0.8 ** 6)) < 0.02
+
+
+def test_worstcase_coverage_is_an_envelope():
+    g = build_graph(SimConfig(n_nodes=400, degree=4, graph_seed=0))
+    rng = np.random.default_rng(0)
+    rand = adversary_metrics(g, place_adversary(g, 0.2, "random", rng, 10**9))["observed_frac"]
+    wc = adversary_metrics(
+        g, place_adversary(g, 0.2, "worstcase_coverage", rng, 10**9))["observed_frac"]
+    assert wc >= rand - 1e-9
+
+
+def test_greedy_coverage_near_optimal():
+    g = build_graph(SimConfig(n_nodes=10, degree=3, graph_seed=0))
+    best = max(adversary_metrics(g, _mask(10, c))["observed_count"]
+               for c in combinations(range(10), 2))
+    idx = _greedy_coverage(g, 2, np.random.default_rng(0))
+    got = adversary_metrics(g, _mask(10, idx))["observed_count"]
+    assert got >= (1 - 1 / np.e) * best - 1e-9
+
+
+def _mask(n, idx):
+    m = np.zeros(n, dtype=bool)
+    m[list(idx)] = True
+    return m
+
+
+def test_worstcase_cap_raises():
+    import pytest
+    g = build_graph(SimConfig(n_nodes=200, degree=4))
+    with pytest.raises(ValueError):
+        place_adversary(g, 0.2, "worstcase_coverage", np.random.default_rng(0), worstcase_max_n=100)

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

@@ -0,0 +1,59 @@
+import dataclasses
+
+import pytest
+
+from pd.config import SimConfig, SweepConfig
+
+
+def test_key_covers_every_field():
+    fields = [f.name for f in dataclasses.fields(SimConfig)]
+    assert len(SimConfig().key()) == len(fields)
+    base = SimConfig()
+    for name in fields:
+        cur = getattr(base, name)
+        alt = {"n_nodes": 2000, "degree": 4, "blend_hops": 2, "max_blend_delay": 5,
+               "unresponsive_frac": 0.2, "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,
+               "coverage_pcts": (25.0,), "f_adv": 0.1, "adversary_mode": "worstcase_coverage",
+               "n_placements": 1, "worstcase_max_n": 5, "graph_seed": 99, "replicate": 1,
+               "root_seed": 7}[name]
+        assert alt != cur
+        assert dataclasses.replace(base, **{name: alt}).key() != base.key(), name
+
+
+@pytest.mark.parametrize("kw", [
+    {"n_nodes": 999},                       # odd
+    {"degree": 1000},                       # >= n
+    {"blend_hops": 0},                      # < 1
+    {"f_adv": 1.0},                         # >= 1
+    {"max_blend_delay": -1},
+    {"processing_lag_probs": (0.5, 0.4)},   # doesn't sum to 1 (with default 3 lags -> len mismatch)
+    {"link_latency_dist": "bogus"},
+    {"adversary_mode": "bogus"},
+])
+def test_validation_rejects(kw):
+    with pytest.raises(ValueError):
+        SimConfig(**kw)
+
+
+def test_sweep_grids_and_collapse():
+    sw = SweepConfig(n_nodes=[1000, 10000], degree=[4, 8], blend_hops=[2, 3],
+                     max_blend_delay=[0, 3], f_adv=[0.0, 0.2],
+                     adversary_mode=["random", "worstcase_coverage"], seeds=3)
+    assert len(sw.graph_cells()) == 2 * 2 * 3
+    assert len(sw.prop_grid()) == 2 * 2
+    # f_adv=0 collapses to a single (mode-irrelevant) row; f_adv=0.2 keeps both modes
+    assert sw.adv_grid() == [(0.0, "random"), (0.2, "random"), (0.2, "worstcase_coverage")]
+
+
+def test_from_dict_rejects_unknown():
+    with pytest.raises(ValueError):
+        SweepConfig.from_dict({"nonsense": [1]})
+
+
+def test_base_config_coerces_tuples():
+    sw = SweepConfig(base={"processing_lags_ms": [10.0, 90.0], "processing_lag_probs": [0.3, 0.7]})
+    cfg = sw.base_config(1000, 8, 0)
+    assert cfg.processing_lags_ms == (10.0, 90.0)
+    assert isinstance(cfg.key(), tuple)

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

@@ -0,0 +1,100 @@
+"""Deanonymization metrics: exact closed forms + a Monte-Carlo tie to the actual draw.
+
+A *deanonymization* event is a round whose whole blend path is adversarial; *full* deanonymization
+additionally requires the honest sender to be directly peered with an adversary. Relays are drawn
+uniformly blind to who is adversarial, so both rates are exact (no sampling in production)."""
+
+import numpy as np
+
+from pd.adversary import adversary_metrics, deanon_metrics, place_adversary
+from pd.config import SimConfig
+from pd.engine import run_graph_cell
+from pd.graph import build_graph
+
+
+def test_deanon_rate_hand_computed():
+    # n=4, 2 adversaries, honest sender leaves 3 nodes (2 adversarial) in the relay pool;
+    # k=2 distinct relays both adversarial: C(2,2)/C(3,2) = 1/3.
+    dz = deanon_metrics(n=4, n_adv=2, observed_frac=0.5, blend_hops=2)
+    assert abs(dz["deanon_rate"] - 1.0 / 3.0) < 1e-12
+    assert abs(dz["full_deanon_rate"] - (1.0 / 3.0) * 0.5) < 1e-12
+
+
+def test_deanon_rate_zero_when_too_few_adversaries():
+    assert deanon_metrics(n=100, n_adv=1, observed_frac=0.9, blend_hops=2)["deanon_rate"] == 0.0
+    assert deanon_metrics(n=100, n_adv=0, observed_frac=0.0, blend_hops=1)["deanon_rate"] == 0.0
+    # too few adversaries -> no full deanonymization either
+    too_few = deanon_metrics(n=100, n_adv=1, observed_frac=0.9, blend_hops=2)
+    assert too_few["full_deanon_rate"] == 0.0
+
+
+def test_full_deanon_is_deanon_times_observed():
+    dz = deanon_metrics(n=5000, n_adv=1000, observed_frac=0.73, blend_hops=3)
+    assert abs(dz["full_deanon_rate"] - dz["deanon_rate"] * 0.73) < 1e-12
+    assert dz["full_deanon_rate"] <= dz["deanon_rate"] + 1e-12
+
+
+def test_deanon_rate_is_placement_independent_but_full_is_not():
+    """The whole-path-adversarial rate depends only on the adversary COUNT; the full rate also
+    tracks how many honest nodes are peered with an adversary, which the worst case maximizes."""
+    g = build_graph(SimConfig(n_nodes=2000, degree=6, graph_seed=0))
+    rng = np.random.default_rng(0)
+    rand = adversary_metrics(g, place_adversary(g, 0.2, "random", rng, 10**9))
+    wc = adversary_metrics(g, place_adversary(g, 0.2, "worstcase_coverage", rng, 10**9))
+    assert rand["n_adv"] == wc["n_adv"]                       # same budget
+    dz_rand = deanon_metrics(g.n, rand["n_adv"], rand["observed_frac"], 3)
+    dz_wc = deanon_metrics(g.n, wc["n_adv"], wc["observed_frac"], 3)
+    assert abs(dz_rand["deanon_rate"] - dz_wc["deanon_rate"]) < 1e-12   # placement-independent
+    assert dz_wc["full_deanon_rate"] >= dz_rand["full_deanon_rate"] - 1e-12  # worst case >= random
+
+
+def test_deanon_asymptotic_fadv_power():
+    # C(A,k)/C(n-1,k) -> f_adv^k for large n.
+    f, k, n = 0.3, 3, 20000
+    dz = deanon_metrics(n=n, n_adv=int(round(f * n)), observed_frac=0.5, blend_hops=k)
+    assert abs(dz["deanon_rate"] - f ** k) < 0.002
+
+
+def test_deanon_matches_direct_sampling():
+    """Closed form == empirical rate of the exact honest-sender/blind-relay draw the sim uses."""
+    f, k = 0.33, 2
+    cfg = SimConfig(n_nodes=1500, degree=8, graph_seed=3, f_adv=f, blend_hops=k)
+    g = build_graph(cfg)
+    mask = place_adversary(g, f, "random", np.random.default_rng(1), cfg.worstcase_max_n)
+    adv = adversary_metrics(g, mask)
+    dz = deanon_metrics(g.n, adv["n_adv"], adv["observed_frac"], k)
+
+    counts = np.add.reduceat(mask[g.indices].astype(np.int32), g.indptr[:-1])
+    observed_node = counts >= 1
+    honest = np.where(~mask)[0]
+    n = g.n
+    rng = np.random.default_rng(42)
+    trials, d_hit, fd_hit = 40_000, 0, 0
+    for _ in range(trials):
+        s = int(rng.choice(honest))
+        r = rng.choice(n - 1, size=k, replace=False)
+        r[r >= s] += 1
+        if mask[r].all():
+            d_hit += 1
+            fd_hit += int(observed_node[s])
+    assert abs(dz["deanon_rate"] - d_hit / trials) < max(0.006, 0.1 * dz["deanon_rate"])
+    assert abs(dz["full_deanon_rate"] - fd_hit / trials) < max(0.006, 0.12 * dz["full_deanon_rate"])
+
+
+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], adv_grid)
+
+    # one deanon row per (placement, distinct blend_hops)
+    assert len(deanon_rows) == len(adv_rows) * 2
+    cols = {"n_nodes", "degree", "blend_hops", "f_adv", "adversary_mode", "graph_seed",
+            "placement_rep", "n_adv", "n_honest", "observed_frac",
+            "deanon_rate", "full_deanon_rate"}
+    assert cols <= set(deanon_rows[0])
+    assert {row["blend_hops"] for row in deanon_rows} == {2, 3}
+    for row in deanon_rows:
+        assert 0.0 <= row["full_deanon_rate"] <= row["deanon_rate"] + 1e-12
+        if row["f_adv"] == 0.0:
+            assert row["deanon_rate"] == 0.0            # no adversary -> no deanonymization

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

@@ -0,0 +1,57 @@
+import numpy as np
+import pytest
+
+from pd.config import SimConfig
+from pd.graph import build_graph, build_regular_edges
+
+
+@pytest.mark.parametrize("n,degree", [(10, 1), (10, 2), (10, 3), (100, 4), (100, 7),
+                                      (1000, 8), (500, 16), (256, 15)])
+def test_exactly_d_regular_simple_undirected(n, degree):
+    rng = np.random.default_rng(0)
+    edges = build_regular_edges(n, degree, rng)
+    deg = np.bincount(edges.ravel(), minlength=n)
+    assert np.all(deg == degree), "every node must have exactly `degree` peers"
+    assert np.all(edges[:, 0] < edges[:, 1]), "no self-loops; canonical u<v"
+    # no duplicate undirected edges
+    keys = edges[:, 0] * n + edges[:, 1]
+    assert len(np.unique(keys)) == len(keys), "graph must be simple"
+    assert edges.shape[0] == n * degree // 2
+
+
+def test_graph_symmetric_and_connected():
+    g = build_graph(SimConfig(n_nodes=2000, degree=6))
+    assert np.all(np.diff(g.indptr) == g.degree)
+    csr = g.weighted_csr(np.ones_like(g.base))
+    assert (csr != csr.T).nnz == 0, "adjacency must be symmetric"
+    from scipy.sparse.csgraph import connected_components
+    ncomp, _ = connected_components(csr, directed=False)
+    assert ncomp == 1, "degree>=3 random d-regular should be connected"
+
+
+def test_reconstructible_from_seed():
+    a = build_graph(SimConfig(n_nodes=1000, degree=8, graph_seed=7))
+    b = build_graph(SimConfig(n_nodes=1000, degree=8, graph_seed=7))
+    assert np.array_equal(a.indptr, b.indptr)
+    assert np.array_equal(a.indices, b.indices)
+    assert np.array_equal(a.base, b.base)
+    assert np.array_equal(a.p, b.p)
+    c = build_graph(SimConfig(n_nodes=1000, degree=8, graph_seed=8))
+    assert not np.array_equal(a.indices, c.indices), "different seed -> different topology"
+
+
+def test_topology_invariant_to_adversary_and_blend_fields():
+    a = build_graph(SimConfig(n_nodes=800, degree=6, graph_seed=1, f_adv=0.0, blend_hops=2))
+    b = build_graph(SimConfig(n_nodes=800, degree=6, graph_seed=1, f_adv=0.4,
+                              blend_hops=5, adversary_mode="worstcase_coverage"))
+    assert np.array_equal(a.indices, b.indices)
+    assert np.array_equal(a.p, b.p)
+
+
+def test_processing_lags_follow_distribution():
+    g = build_graph(SimConfig(n_nodes=20000, degree=6,
+                              processing_lags_ms=(10.0, 50.0, 100.0),
+                              processing_lag_probs=(0.5, 0.4, 0.1)))
+    for lag, prob in zip((10.0, 50.0, 100.0), (0.5, 0.4, 0.1), strict=True):
+        frac = np.mean(g.p == lag)
+        assert abs(frac - prob) < 0.03, f"lag {lag}: {frac:.3f} vs {prob}"

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

@@ -0,0 +1,24 @@
+import numpy as np
+
+from pd.mixclock import mean_residual_ms, mix_wait
+
+
+def test_zero_max_delay_is_zero():
+    rng = np.random.default_rng(0)
+    w = mix_wait(rng, 0, 1000)
+    assert np.all(w == 0.0)
+
+
+def test_residual_within_bounds():
+    rng = np.random.default_rng(1)
+    m = 5
+    w = mix_wait(rng, m, 100000)
+    assert w.min() >= 0.0
+    assert w.max() <= m * 1000.0 + 1e-6  # residual within a covering interval (<= M seconds)
+
+
+def test_mean_matches_analytic():
+    rng = np.random.default_rng(2)
+    for m in (1, 3, 8):
+        w = mix_wait(rng, m, 400000)
+        assert abs(w.mean() - mean_residual_ms(m)) < 0.03 * mean_residual_ms(m)

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

@@ -0,0 +1,95 @@
+import numpy as np
+
+from pd.graph import Graph
+from pd.propagation import assign_responsive, blend_round
+
+
+def _k4(p):
+    """Complete graph on 4 nodes (degree 3), base latency 10 ms on every link, node lags `p`."""
+    indptr = np.array([0, 3, 6, 9, 12], dtype=np.int64)
+    indices = np.array([1, 2, 3, 0, 2, 3, 0, 1, 3, 0, 1, 2], dtype=np.int64)
+    base = np.full(12, 10.0)
+    src = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3], dtype=np.int64)
+    return Graph(n=4, degree=3, indptr=indptr, indices=indices, base=base, src=src,
+                 p=np.asarray(p, dtype=float))
+
+
+def test_single_relay_delay_with_node_lags():
+    # jitter=0, max_blend_delay=0. Directed edge (u->v) = base(10) + p(u).
+    g = _k4([1.0, 2.0, 3.0, 4.0])
+    rng = np.random.default_rng(0)
+    r = blend_round(g, sender=0, relays=np.array([1]), jitter_mean_ms=0.0,
+                    max_blend_delay=0, rng=rng, coverage_pcts=(50.0, 90.0, 99.0))
+    # leg 0->1 = 10 + p(0) = 11 ; broadcast from 1 to farthest = 10 + p(1) = 12
+    assert r["path"] == 11.0
+    assert r["broadcast"] == 12.0
+    assert r["full"] == 23.0
+    assert r["frac_reached"] == 1.0
+
+
+def test_two_relay_path_sums_legs():
+    g = _k4([1.0, 2.0, 3.0, 4.0])
+    rng = np.random.default_rng(0)
+    r = blend_round(g, sender=0, relays=np.array([1, 2]), jitter_mean_ms=0.0,
+                    max_blend_delay=0, rng=rng, coverage_pcts=(50.0,))
+    # legs: 0->1 = 11, 1->2 = 10 + p(1) = 12  => path 23 ; broadcast from 2 = 10 + p(2) = 13
+    assert r["path"] == 23.0
+    assert r["broadcast"] == 13.0
+    assert r["full"] == 36.0
+
+
+def test_mixing_adds_positive_delay():
+    g = _k4([0.0, 0.0, 0.0, 0.0])
+    rng = np.random.default_rng(1)
+    no_mix = blend_round(g, 0, np.array([1]), 0.0, 0, rng, (50.0,))["full"]
+    mixed = np.mean([blend_round(g, 0, np.array([1]), 0.0, 5, rng, (50.0,))["full"]
+                     for _ in range(500)])
+    assert mixed > no_mix  # the free-running clock adds a positive mixing residual
+
+
+def _path4():
+    """Line graph 0-1-2-3 (base 10 ms each way, no node lags)."""
+    indptr = np.array([0, 1, 3, 5, 6], dtype=np.int64)
+    indices = np.array([1, 0, 2, 1, 3, 2], dtype=np.int64)
+    base = np.full(6, 10.0)
+    src = np.array([0, 1, 1, 2, 2, 3], dtype=np.int64)
+    return Graph(n=4, degree=2, indptr=indptr, indices=indices, base=base, src=src,
+                 p=np.zeros(4))
+
+
+def test_assign_responsive_count_and_edges():
+    rng = np.random.default_rng(0)
+    mask = assign_responsive(1000, 0.3, rng)
+    assert mask.dtype == bool
+    assert int(mask.sum()) == 700                      # exactly 30% dropped
+    assert assign_responsive(1000, 0.0, rng).all()     # frac 0 -> everyone responsive
+
+
+def test_unresponsive_final_relay_drops_message():
+    # final relay (node 1) unresponsive -> it receives but cannot flood: not delivered.
+    g = _k4([0.0, 0.0, 0.0, 0.0])
+    responsive = np.array([True, False, True, True])
+    r = blend_round(g, sender=0, relays=np.array([1]), jitter_mean_ms=0.0, max_blend_delay=0,
+                    rng=np.random.default_rng(0), coverage_pcts=(50.0,), responsive=responsive)
+    assert r["delivered"] is False
+    assert np.isnan(r["full"])
+
+
+def test_unresponsive_intermediate_relay_drops_message():
+    # first relay (node 1) unresponsive -> the second leg 1->2 is inf: not delivered.
+    g = _k4([0.0, 0.0, 0.0, 0.0])
+    responsive = np.array([True, False, True, True])
+    r = blend_round(g, sender=0, relays=np.array([1, 2]), jitter_mean_ms=0.0, max_blend_delay=0,
+                    rng=np.random.default_rng(0), coverage_pcts=(50.0,), responsive=responsive)
+    assert r["delivered"] is False
+
+
+def test_unresponsive_node_strands_flood_pocket():
+    # path 0-1-2-3; relay 1 is responsive so the message is delivered, but node 2 is a routing hole
+    # so node 3 (only reachable through 2) never receives the flood.
+    g = _path4()
+    responsive = np.array([True, True, False, True])
+    r = blend_round(g, sender=0, relays=np.array([1]), jitter_mean_ms=0.0, max_blend_delay=0,
+                    rng=np.random.default_rng(0), coverage_pcts=(50.0,), responsive=responsive)
+    assert r["delivered"] is True
+    assert r["frac_reached"] == 0.75                   # node 3 stranded behind unresponsive node 2

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

@@ -0,0 +1,39 @@
+from pd.config import SimConfig
+from pd.rng import (
+    graph_seedseq,
+    responsive_seedseq,
+    rng_for,
+    round_seedseq,
+    seedseq_for,
+)
+
+
+def test_deterministic():
+    c = SimConfig(n_nodes=1000, degree=8, graph_seed=3)
+    assert seedseq_for(c).entropy == seedseq_for(c).entropy
+    a = rng_for(c).integers(0, 10**9, size=5)
+    b = rng_for(c).integers(0, 10**9, size=5)
+    assert list(a) == list(b)
+
+
+def test_graph_seed_is_topology_only():
+    base = SimConfig(n_nodes=1000, degree=8, graph_seed=3, f_adv=0.0, blend_hops=2)
+    same_topo = SimConfig(n_nodes=1000, degree=8, graph_seed=3, f_adv=0.4, blend_hops=5,
+                          adversary_mode="worstcase_coverage", n_rounds=999)
+    assert graph_seedseq(base).entropy == graph_seedseq(same_topo).entropy
+    diff_topo = SimConfig(n_nodes=1000, degree=8, graph_seed=4)
+    assert graph_seedseq(base).entropy != graph_seedseq(diff_topo).entropy
+
+
+def test_responsive_seed_depends_on_frac_only():
+    c = SimConfig(n_nodes=1000, degree=8, graph_seed=3)
+    # fixed per (topology, unresponsive_frac); different frac -> different responsive draw
+    assert responsive_seedseq(c, 0.1).entropy == responsive_seedseq(c, 0.1).entropy
+    assert responsive_seedseq(c, 0.1).entropy != responsive_seedseq(c, 0.2).entropy
+
+
+def test_round_seed_includes_unresponsive_frac():
+    c = SimConfig(n_nodes=1000, degree=8, graph_seed=3)
+    a = round_seedseq(c, blend_hops=3, max_blend_delay=3, unresponsive_frac=0.0)
+    b = round_seedseq(c, blend_hops=3, max_blend_delay=3, unresponsive_frac=0.2)
+    assert a.entropy != b.entropy