diff --git a/tools/simulators/blend/pd/.gitignore b/tools/simulators/blend/pd/.gitignore new file mode 100644 index 0000000..8a18509 --- /dev/null +++ b/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 diff --git a/tools/simulators/blend/pd/Makefile b/tools/simulators/blend/pd/Makefile new file mode 100644 index 0000000..4e8ee44 --- /dev/null +++ b/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/ + $(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 diff --git a/tools/simulators/blend/pd/README.md b/tools/simulators/blend/pd/README.md new file mode 100644 index 0000000..282ab24 --- /dev/null +++ b/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 -m pd.sweep ... +make smoke # fast end-to-end -> runs/_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/ +``` + +## 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/`. diff --git a/tools/simulators/blend/pd/configs/default.yaml b/tools/simulators/blend/pd/configs/default.yaml new file mode 100644 index 0000000..850ef76 --- /dev/null +++ b/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 diff --git a/tools/simulators/blend/pd/configs/fullscale.yaml b/tools/simulators/blend/pd/configs/fullscale.yaml new file mode 100644 index 0000000..92e2f23 --- /dev/null +++ b/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 diff --git a/tools/simulators/blend/pd/configs/smoke.yaml b/tools/simulators/blend/pd/configs/smoke.yaml new file mode 100644 index 0000000..0a3ebca --- /dev/null +++ b/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 diff --git a/tools/simulators/blend/pd/pyproject.toml b/tools/simulators/blend/pd/pyproject.toml new file mode 100644 index 0000000..458a292 --- /dev/null +++ b/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')"] diff --git a/tools/simulators/blend/pd/requirements-dev.txt b/tools/simulators/blend/pd/requirements-dev.txt new file mode 100644 index 0000000..aefbcb6 --- /dev/null +++ b/tools/simulators/blend/pd/requirements-dev.txt @@ -0,0 +1 @@ +-e .[dev] diff --git a/tools/simulators/blend/pd/requirements.txt b/tools/simulators/blend/pd/requirements.txt new file mode 100644 index 0000000..d6e1198 --- /dev/null +++ b/tools/simulators/blend/pd/requirements.txt @@ -0,0 +1 @@ +-e . diff --git a/tools/simulators/blend/pd/scripts/make_figures.py b/tools/simulators/blend/pd/scripts/make_figures.py new file mode 100644 index 0000000..be5fbe5 --- /dev/null +++ b/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()) diff --git a/tools/simulators/blend/pd/scripts/run_sweep.py b/tools/simulators/blend/pd/scripts/run_sweep.py new file mode 100644 index 0000000..9e18f86 --- /dev/null +++ b/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()) diff --git a/tools/simulators/blend/pd/scripts/verify.py b/tools/simulators/blend/pd/scripts/verify.py new file mode 100644 index 0000000..7863933 --- /dev/null +++ b/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()) diff --git a/tools/simulators/blend/pd/src/pd/__init__.py b/tools/simulators/blend/pd/src/pd/__init__.py new file mode 100644 index 0000000..8214a75 --- /dev/null +++ b/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" diff --git a/tools/simulators/blend/pd/src/pd/adversary.py b/tools/simulators/blend/pd/src/pd/adversary.py new file mode 100644 index 0000000..8a37766 --- /dev/null +++ b/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 diff --git a/tools/simulators/blend/pd/src/pd/config.py b/tools/simulators/blend/pd/src/pd/config.py new file mode 100644 index 0000000..7f75376 --- /dev/null +++ b/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}) diff --git a/tools/simulators/blend/pd/src/pd/constants.py b/tools/simulators/blend/pd/src/pd/constants.py new file mode 100644 index 0000000..acb0034 --- /dev/null +++ b/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) +) diff --git a/tools/simulators/blend/pd/src/pd/engine.py b/tools/simulators/blend/pd/src/pd/engine.py new file mode 100644 index 0000000..26417b6 --- /dev/null +++ b/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} diff --git a/tools/simulators/blend/pd/src/pd/graph.py b/tools/simulators/blend/pd/src/pd/graph.py new file mode 100644 index 0000000..9f0b2bb --- /dev/null +++ b/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) diff --git a/tools/simulators/blend/pd/src/pd/latency.py b/tools/simulators/blend/pd/src/pd/latency.py new file mode 100644 index 0000000..eeb3846 --- /dev/null +++ b/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] diff --git a/tools/simulators/blend/pd/src/pd/memguard.py b/tools/simulators/blend/pd/src/pd/memguard.py new file mode 100644 index 0000000..5a4bddf --- /dev/null +++ b/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 ''}") diff --git a/tools/simulators/blend/pd/src/pd/metrics.py b/tools/simulators/blend/pd/src/pd/metrics.py new file mode 100644 index 0000000..65b3a20 --- /dev/null +++ b/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, + } diff --git a/tools/simulators/blend/pd/src/pd/mixclock.py b/tools/simulators/blend/pd/src/pd/mixclock.py new file mode 100644 index 0000000..d480af1 --- /dev/null +++ b/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 diff --git a/tools/simulators/blend/pd/src/pd/plotting/__init__.py b/tools/simulators/blend/pd/src/pd/plotting/__init__.py new file mode 100644 index 0000000..2d0f432 --- /dev/null +++ b/tools/simulators/blend/pd/src/pd/plotting/__init__.py @@ -0,0 +1 @@ +"""Plotting: shared style + figure builders for pd.""" diff --git a/tools/simulators/blend/pd/src/pd/plotting/figures.py b/tools/simulators/blend/pd/src/pd/plotting/figures.py new file mode 100644 index 0000000..0d176e9 --- /dev/null +++ b/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 diff --git a/tools/simulators/blend/pd/src/pd/plotting/make_figures.py b/tools/simulators/blend/pd/src/pd/plotting/make_figures.py new file mode 100644 index 0000000..9b4506a --- /dev/null +++ b/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/_