diff --git a/tools/simulators/tsi/tsi-sim/.gitignore b/tools/simulators/tsi/tsi-sim/.gitignore
new file mode 100644
index 0000000..c4ee506
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/.gitignore
@@ -0,0 +1,18 @@
+# Generated artifacts
+runs/
+results/
+figures/
+!results/.gitkeep
+!figures/.gitkeep
+
+# Python
+__pycache__/
+*.py[cod]
+*.egg-info/
+.venv/
+venv/
+.pytest_cache/
+.mypy_cache/
+.ruff_cache/
+*.parquet
+*.csv
diff --git a/tools/simulators/tsi/tsi-sim/Makefile b/tools/simulators/tsi/tsi-sim/Makefile
new file mode 100644
index 0000000..a1b8ae6
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/Makefile
@@ -0,0 +1,41 @@
+.PHONY: install smoke sweep figures verify test lint clean
+
+VENV ?= .venv
+PY   := $(VENV)/bin/python
+PIP  := $(VENV)/bin/pip
+
+$(VENV):
+	python3 -m venv $(VENV)
+	$(PIP) install --upgrade pip
+
+install: $(VENV)
+	$(PIP) install -e ".[dev]"
+
+# Every sweep writes results + figures into a fresh dated folder under runs/ so new runs
+# never overwrite old ones (runs/<YYYY-MM-DD_HHMMSS>_<label>/).
+
+# Fast end-to-end check (tiny scaled-k grid); auto-generates figures.
+smoke: install
+	$(PY) scripts/run_sweep.py --config configs/smoke.yaml --label smoke
+
+# Full scaled-k parameter sweep (auto-generates figures).
+sweep: install
+	$(PY) scripts/run_sweep.py --config configs/default.yaml --label default
+
+# Re-render figures from an existing run's parquet into a fresh dated figures/ folder.
+figures: install
+	$(PY) scripts/make_figures.py --results $(RUN)/results.parquet
+
+# Analytic sanity checks (simulator vs closed-form theory)
+verify: install
+	$(PY) scripts/verify.py
+
+test: install
+	$(PY) -m pytest
+
+lint: install
+	$(VENV)/bin/ruff check src scripts tests
+
+clean:
+	rm -rf runs/* results/*.parquet figures/* .pytest_cache .ruff_cache .mypy_cache
+	find . -name __pycache__ -type d -prune -exec rm -rf {} +
diff --git a/tools/simulators/tsi/tsi-sim/README.md b/tools/simulators/tsi/tsi-sim/README.md
new file mode 100644
index 0000000..4542a71
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/README.md
@@ -0,0 +1,58 @@
+# tsi-sim — Cryptarchia Total Stake Inference simulator (uncle references)
+
+Monte-Carlo simulation framework used to choose safe values for the **uncle-reference**
+parameters of Cryptarchia's Total Stake Inference (TSI):
+
+- `U` — max uncles referenced per block (`MAX_UNCLES`); `U=0` is the no-uncle baseline.
+- `W` — uncle reference window in slots (spec default 300).
+- swept against network size `N`, stake distribution `S` (uniform / Pareto), and
+  network **latency** `L` (in slots — deliberately *not* `D`, which denotes the stake estimate).
+
+It measures how well the inferred total active stake `D` tracks the true total stake, and
+whether uncle references recover the active slots that network latency loses to forks.
+
+> This lives under a `raw/` docs path, so it is invisible to the repository's
+> markdown-lint CI. It is a standalone Python package with its own tooling.
+
+## Model
+
+Reduced **canonical-chain-with-orphans** model: we simulate the global winning-slot
+sequence (stake-weighted φ lottery), build a real block tree with latency- and
+multi-winner-induced forks, resolve the canonical chain (honest longest-chain), let
+canonical blocks reference uncles per the spec's selection rules, and count TSI density
+`m = honest-chain blocks + deduplicated referenced uncles` in the measurement window.
+All honest nodes converge to the same deep chain (k-finality), so a single per-epoch
+`D` is faithful. A full per-node model is the planned next phase (`per_node_dest` flag
+scaffolds it).
+
+See the sibling spec `../` and `../../cryptarchia-total-stake-inference.md` for the math.
+
+## Quick start
+
+```bash
+make install         # create .venv and install (editable) with dev deps
+make test            # unit tests + fast theory checks
+make verify          # simulator vs closed-form analytic checks
+make smoke           # tiny scaled-k sweep + figures (end-to-end smoke test)
+make sweep figures   # full scaled-k parameter sweep + academic figures
+```
+
+Outputs: `results/*.parquet` (one row per config×epoch) and `figures/*.{pdf,png}`
+(both git-ignored).
+
+## Scale
+
+True constants (`k=2160`, `f=1/30`) give 648,000-slot epochs — too large to sweep.
+Mean accuracy is provably `k`-invariant (only variance scales `~1/T`), so sweeps use a
+**scaled `k`** (`configs/default.yaml`); the final accuracy/variance figures re-run at
+true `k` (`configs/fullscale.yaml`). `configs/smoke.yaml` is a tiny dev grid.
+
+## Layout
+
+```
+src/tsi_sim/   constants config rng stake lottery latency blocktree uncles
+               tsi epoch engine metrics theory sweep  plotting/{style,figures}
+scripts/       run_sweep.py  make_figures.py  verify.py
+configs/       smoke.yaml  default.yaml  fullscale.yaml
+tests/         test_{lottery,uncles,tsi_counting,blocktree,theory_convergence}.py
+```
diff --git a/tools/simulators/tsi/tsi-sim/configs/default.yaml b/tools/simulators/tsi/tsi-sim/configs/default.yaml
new file mode 100644
index 0000000..051dca9
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/configs/default.yaml
@@ -0,0 +1,16 @@
+# Full scaled-k parameter sweep. Mean accuracy / q_eff / convergence are k-invariant,
+# so k is scaled down for tractability; re-run headline + variance figures at full scale
+# with fullscale.yaml.
+# Sweep axes (cartesian product x replicates); every value is a list.
+n_nodes: [1000, 2000, 4000]      # number of nodes / stake holders
+stake_dist: [uniform, pareto]    # stake distribution (uniform = equal, pareto = heavy-tailed)
+latency: [0, 1, 2, 4, 8, 16]     # L: network latency in slots; block visible to others at t+L
+max_uncles: [0, 1, 2, 3, 4]      # U: max uncle references per block (0 = baseline, no uncles)
+uncle_strategy: [oldest, random] # uncle selection: oldest-first fill vs random coin-flip
+replicates: 12                   # independent RNG replicates per grid cell
+base:                            # per-run settings shared by every cell (not swept)
+  k: 64                          # security parameter (scaled; T = 6*floor(64/f) = 11520 slots)
+  epochs: 45                     # epochs simulated per trajectory
+  f: 0.03333333333333333        # slot activation coefficient (default 1/30); configurable
+  genesis_d_factor: 0.5          # genesis D_est = factor x true total stake
+  pareto_shape: 1.16             # Pareto (Lomax) tail index for the pareto stake distribution
diff --git a/tools/simulators/tsi/tsi-sim/configs/fullscale.yaml b/tools/simulators/tsi/tsi-sim/configs/fullscale.yaml
new file mode 100644
index 0000000..f680d3a
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/configs/fullscale.yaml
@@ -0,0 +1,16 @@
+# Full-scale (true k=2160) runs for the final headline + variance figures only.
+# ~100x slower per config than default.yaml, so the grid is deliberately small.
+# Note: T = 6*floor(2160/f) = 388,800 slots per epoch — expect minutes per config.
+# Sweep axes (cartesian product x replicates); every value is a list.
+n_nodes: [1000]                  # number of nodes / stake holders
+stake_dist: [uniform, pareto]    # stake distribution (uniform = equal, pareto = heavy-tailed)
+latency: [0, 2, 8]               # L: network latency in slots; block visible to others at t+L
+max_uncles: [0, 1, 2, 4]         # U: max uncle references per block (0 = baseline, no uncles)
+uncle_strategy: [oldest]         # uncle selection: oldest-first fill
+replicates: 4                    # independent RNG replicates per grid cell
+base:                            # per-run settings shared by every cell (not swept)
+  k: 2160                        # true security parameter
+  epochs: 12                     # epochs simulated per trajectory
+  f: 0.03333333333333333        # slot activation coefficient (default 1/30); configurable
+  genesis_d_factor: 0.5          # genesis D_est = factor x true total stake
+  pareto_shape: 1.16             # Pareto (Lomax) tail index for the pareto stake distribution
diff --git a/tools/simulators/tsi/tsi-sim/configs/smoke.yaml b/tools/simulators/tsi/tsi-sim/configs/smoke.yaml
new file mode 100644
index 0000000..15947fd
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/configs/smoke.yaml
@@ -0,0 +1,14 @@
+# Tiny end-to-end smoke grid (fast; for dev + CI-style checks).
+# Sweep axes (cartesian product x replicates); every value is a list.
+n_nodes: [1000]                  # number of nodes / stake holders
+stake_dist: [uniform, pareto]    # stake distribution (uniform = equal, pareto = heavy-tailed)
+latency: [0, 4, 8]               # L: network latency in slots; block visible to others at t+L
+max_uncles: [0, 1, 4]            # U: max uncle references per block (0 = baseline, no uncles)
+uncle_strategy: [oldest]         # uncle selection: oldest-first fill
+replicates: 4                    # independent RNG replicates per grid cell
+base:                            # per-run settings shared by every cell (not swept)
+  k: 16                          # security parameter (scaled; T = 6*floor(16/f) = 2880 slots)
+  epochs: 25                     # epochs simulated per trajectory
+  f: 0.03333333333333333        # slot activation coefficient (default 1/30); configurable
+  genesis_d_factor: 0.5          # genesis D_est = factor x true total stake
+  pareto_shape: 1.16             # Pareto (Lomax) tail index for the pareto stake distribution
diff --git a/tools/simulators/tsi/tsi-sim/pyproject.toml b/tools/simulators/tsi/tsi-sim/pyproject.toml
new file mode 100644
index 0000000..380065b
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/pyproject.toml
@@ -0,0 +1,44 @@
+[build-system]
+requires = ["hatchling"]
+build-backend = "hatchling.build"
+
+[project]
+name = "tsi-sim"
+version = "0.1.0"
+description = "Monte-Carlo simulation of Cryptarchia Total Stake Inference with uncle references"
+readme = "README.md"
+requires-python = ">=3.11"
+license = { text = "MIT" }
+authors = [{ name = "Logos" }]
+dependencies = [
+    "numpy>=1.26",
+    "pandas>=2.1",
+    "pyarrow>=14",
+    "matplotlib>=3.8",
+    "scipy>=1.11",
+    "pyyaml>=6.0",
+    "tqdm>=4.66",
+    "joblib>=1.3",
+]
+
+[project.optional-dependencies]
+dev = ["pytest>=8.0", "pytest-xdist>=3.5", "ruff>=0.5", "mypy>=1.8"]
+
+[project.scripts]
+tsi-sweep = "tsi_sim.sweep:main"
+
+[tool.hatch.build.targets.wheel]
+packages = ["src/tsi_sim"]
+
+[tool.ruff]
+line-length = 100
+target-version = "py311"
+
+[tool.ruff.lint]
+select = ["E", "F", "I", "UP", "B", "NPY"]
+ignore = ["E741"]  # allow single-char names like L (latency), q, m — domain symbols
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
+markers = ["slow: long-running (full-scale k) checks"]
+addopts = "-q"
diff --git a/tools/simulators/tsi/tsi-sim/requirements-dev.txt b/tools/simulators/tsi/tsi-sim/requirements-dev.txt
new file mode 100644
index 0000000..2a39dbe
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/requirements-dev.txt
@@ -0,0 +1,5 @@
+-r requirements.txt
+pytest>=8.0
+pytest-xdist>=3.5
+ruff>=0.5
+mypy>=1.8
diff --git a/tools/simulators/tsi/tsi-sim/requirements.txt b/tools/simulators/tsi/tsi-sim/requirements.txt
new file mode 100644
index 0000000..3212e73
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/requirements.txt
@@ -0,0 +1,8 @@
+numpy>=1.26
+pandas>=2.1
+pyarrow>=14
+matplotlib>=3.8
+scipy>=1.11
+pyyaml>=6.0
+tqdm>=4.66
+joblib>=1.3
diff --git a/tools/simulators/tsi/tsi-sim/scripts/make_figures.py b/tools/simulators/tsi/tsi-sim/scripts/make_figures.py
new file mode 100644
index 0000000..314a219
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/scripts/make_figures.py
@@ -0,0 +1,14 @@
+#!/usr/bin/env python
+"""Thin shim so `python scripts/make_figures.py` works without installing."""
+
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
+
+from tsi_sim.plotting.make_figures import main  # noqa: E402
+
+if __name__ == "__main__":
+    main()
diff --git a/tools/simulators/tsi/tsi-sim/scripts/run_sweep.py b/tools/simulators/tsi/tsi-sim/scripts/run_sweep.py
new file mode 100644
index 0000000..77d98a0
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/scripts/run_sweep.py
@@ -0,0 +1,14 @@
+#!/usr/bin/env python
+"""Run a TSI parameter sweep from a YAML config -> results parquet."""
+
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
+
+from tsi_sim.sweep import main  # noqa: E402
+
+if __name__ == "__main__":
+    main()
diff --git a/tools/simulators/tsi/tsi-sim/scripts/verify.py b/tools/simulators/tsi/tsi-sim/scripts/verify.py
new file mode 100644
index 0000000..fa75f1f
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/scripts/verify.py
@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+"""Analytic sanity checks: simulator vs closed-form theory.
+
+Run with the project venv:  python scripts/verify.py
+Exits non-zero if any check fails.
+"""
+
+from __future__ import annotations
+
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
+
+import numpy as np  # noqa: E402
+import pandas as pd  # noqa: E402
+
+from tsi_sim.config import SimConfig  # noqa: E402
+from tsi_sim.engine import run_trajectory  # noqa: E402
+from tsi_sim.epoch import simulate_epoch  # noqa: E402
+from tsi_sim.rng import rng_for  # noqa: E402
+from tsi_sim.stake import make_stake  # noqa: E402
+from tsi_sim.theory import expected_ratio  # noqa: E402
+
+F = 1.0 / 30.0
+K = 128            # scaled: T = 6*floor(128/f) = 23040 slots
+EPOCHS = 45
+REPS = 12
+BURN = 25
+
+
+def tail_mean(cfg: SimConfig, col: str, reps: int = REPS) -> tuple[float, float]:
+    vals = []
+    for r in range(reps):
+        df = pd.DataFrame(run_trajectory(cfg.__class__(**{**cfg.__dict__, "replicate": r})))
+        vals.append(df[col].iloc[BURN:].mean())
+    return float(np.mean(vals)), float(np.std(vals) / np.sqrt(reps))
+
+
+def check(name: str, ok: bool, detail: str) -> bool:
+    print(f"[{'PASS' if ok else 'FAIL'}] {name}: {detail}")
+    return ok
+
+
+def main() -> int:
+    results = []
+
+    # 1. Active-slot rate ~= f when D_est = D_true, L=0, U=0.
+    cfg = SimConfig(n_nodes=2000, stake_dist="uniform", latency=0, max_uncles=0,
+                    k=K, epochs=1, genesis_d_factor=1.0)
+    rng = rng_for(cfg)
+    stake = make_stake(cfg, rng)
+    er = simulate_epoch(cfg, stake, float(stake.sum()), rng)
+    active_rate = er.n_active / cfg.period_T
+    results.append(check("active-slot rate ~ f (L=0,U=0,D=D_true)",
+                         abs(active_rate - F) / F < 0.05,
+                         f"active_rate={active_rate:.5f} f={F:.5f}"))
+
+    # 2. U=0 equilibrium ratio ~= expected_ratio(f, measured q).
+    cfg = SimConfig(n_nodes=1000, stake_dist="uniform", latency=4, max_uncles=0,
+                    k=K, epochs=EPOCHS, genesis_d_factor=0.5)
+    ratio, se = tail_mean(cfg, "ratio")
+    q, _ = tail_mean(cfg, "q")
+    pred = float(expected_ratio(F, q))
+    results.append(check("U=0 ratio ~ theory(q)",
+                         abs(ratio - pred) < 0.02 + 2 * se,
+                         f"sim={ratio:.4f}±{se:.4f} theory(q={q:.3f})={pred:.4f}"))
+
+    # 3. Underestimate at higher latency (q < 1 => ratio < 1), U=0.
+    cfg = SimConfig(n_nodes=1000, stake_dist="uniform", latency=10, max_uncles=0,
+                    k=K, epochs=EPOCHS, genesis_d_factor=0.5)
+    ratio, se = tail_mean(cfg, "ratio")
+    q, _ = tail_mean(cfg, "q")
+    results.append(check("U=0 underestimates true stake at latency",
+                         ratio < 0.98 and q < 0.98,
+                         f"ratio={ratio:.4f} q={q:.3f}"))
+
+    # 4. q_eff -> 1 and ratio -> ~1 as U grows (uncles recover forks).
+    base = dict(n_nodes=1000, stake_dist="uniform", latency=8, uncle_strategy="oldest",
+                k=K, epochs=EPOCHS, genesis_d_factor=0.5)
+    r0, _ = tail_mean(SimConfig(max_uncles=0, **base), "ratio")
+    r4, _ = tail_mean(SimConfig(max_uncles=4, **base), "ratio")
+    qe4, _ = tail_mean(SimConfig(max_uncles=4, **base), "q_eff")
+    results.append(check("uncles recover accuracy (q_eff->1, |ratio-1| shrinks)",
+                         qe4 > 0.99 and abs(r4 - 1) < abs(r0 - 1) and abs(r4 - 1) < 0.03,
+                         f"ratio U0={r0:.4f} -> U4={r4:.4f}; q_eff(U4)={qe4:.4f}"))
+
+    print()
+    n_pass = sum(results)
+    print(f"{n_pass}/{len(results)} checks passed")
+    return 0 if n_pass == len(results) else 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/__init__.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/__init__.py
new file mode 100644
index 0000000..9816edf
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/__init__.py
@@ -0,0 +1,3 @@
+"""Cryptarchia Total Stake Inference simulator (uncle references)."""
+
+__version__ = "0.1.0"
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/blocktree.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/blocktree.py
new file mode 100644
index 0000000..1f807d2
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/blocktree.py
@@ -0,0 +1,127 @@
+"""Block tree, latency-driven forks, and honest longest-chain fork choice.
+
+Blocks are stored in parallel arrays (id == index). A virtual genesis is block 0 at
+slot -1, height 0. Every real block is produced at an active slot by one winning node
+and points at the best tip *visible to that node at production time*, which is what makes
+network latency (and same-slot multi-winners) produce forks.
+
+Fork choice is honest longest-chain with a first-seen tie-break (prefer higher height,
+then earlier slot, then lower id) — no adversary is modelled, so the spec's density /
+deep-fork rules never engage.
+"""
+
+from __future__ import annotations
+
+import heapq
+from dataclasses import dataclass
+
+import numpy as np
+
+from .latency import LatencyModel
+
+GENESIS = 0
+
+
+@dataclass
+class BlockTree:
+    slot: np.ndarray          # int64, slot of each block (genesis = -1)
+    parent: np.ndarray        # int64, parent id (genesis = -1)
+    height: np.ndarray        # int64, chain height (genesis = 0)
+    leader: np.ndarray        # int64, producing node id (genesis = -1)
+    uncles: list[tuple[int, ...]]   # referenced uncle ids per block (filled later)
+
+    @property
+    def n_blocks(self) -> int:
+        return self.slot.shape[0]
+
+    def ancestors(self, block_id: int) -> list[int]:
+        """Ancestor chain of ``block_id`` from itself down to (excluding) genesis."""
+        out: list[int] = []
+        b = block_id
+        while b > GENESIS:
+            out.append(b)
+            b = int(self.parent[b])
+        return out
+
+    def canonical_chain(self) -> list[int]:
+        """Honest longest-chain: ancestors of the best tip over the whole tree.
+
+        Returns real block ids (genesis excluded), tip-first.
+        """
+        tip = self._best_over_all()
+        return self.ancestors(tip)
+
+    def _rank(self, bid: int) -> tuple[int, int, int]:
+        # Preference order for "better tip": higher height, earlier slot, lower id.
+        return (int(self.height[bid]), -int(self.slot[bid]), -bid)
+
+    def _best_over_all(self) -> int:
+        best = GENESIS
+        best_rank = self._rank(GENESIS)
+        for bid in range(1, self.n_blocks):
+            r = self._rank(bid)
+            if r > best_rank:
+                best_rank, best = r, bid
+        return best
+
+
+def build_tree(
+    active_slots: np.ndarray,
+    winners_per_slot: list[np.ndarray],
+    latency: LatencyModel,
+    rng: np.random.Generator,
+) -> BlockTree:
+    """Construct the block tree from grouped lottery winners under a latency model."""
+    # Preallocate with genesis in slot 0.
+    slot = [-1]
+    parent = [-1]
+    height = [0]
+    leader = [-1]
+
+    # global_best = best publicly-visible tip so far, as (height, slot, id).
+    def better(a: tuple[int, int, int], b: tuple[int, int, int]) -> tuple[int, int, int]:
+        # higher height, then earlier slot, then lower id
+        ah, as_, ai = a
+        bh, bs, bi = b
+        if ah != bh:
+            return a if ah > bh else b
+        if as_ != bs:
+            return a if as_ < bs else b
+        return a if ai < bi else b
+
+    global_best = (0, -1, GENESIS)
+    own_best: dict[int, tuple[int, int, int]] = {}
+    # min-heap of (visible_at, block_id) awaiting public visibility
+    pending: list[tuple[int, int]] = []
+
+    next_id = 1
+    for si in range(active_slots.shape[0]):
+        t = int(active_slots[si])
+        # advance visibility frontier to slot t
+        while pending and pending[0][0] <= t:
+            _, bid = heapq.heappop(pending)
+            cand = (height[bid], slot[bid], bid)
+            global_best = better(global_best, cand)
+        for v in winners_per_slot[si].tolist():
+            gb = global_best
+            ob = own_best.get(v, (0, -1, GENESIS))
+            chosen = better(gb, ob)
+            p_id = chosen[2]
+            h = chosen[0] + 1
+            bid = next_id
+            next_id += 1
+            slot.append(t)
+            parent.append(p_id)
+            height.append(h)
+            leader.append(v)
+            own_best[v] = (h, t, bid)
+            va = latency.visible_at(t, rng)
+            heapq.heappush(pending, (va, bid))
+
+    return BlockTree(
+        slot=np.asarray(slot, np.int64),
+        parent=np.asarray(parent, np.int64),
+        height=np.asarray(height, np.int64),
+        leader=np.asarray(leader, np.int64),
+        uncles=[() for _ in range(next_id)],
+    )
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/config.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/config.py
new file mode 100644
index 0000000..d6642d9
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/config.py
@@ -0,0 +1,116 @@
+"""Configuration dataclasses for single runs and parameter sweeps."""
+
+from __future__ import annotations
+
+import itertools
+from dataclasses import dataclass, field, replace
+from typing import Any, Literal
+
+from . import constants
+
+StakeDist = Literal["uniform", "pareto"]
+UncleStrategy = Literal["oldest", "random"]
+
+
+@dataclass(frozen=True)
+class SimConfig:
+    """A single fully-specified simulation run (one grid cell, one replicate)."""
+
+    # --- network / stake ---
+    n_nodes: int = 1000
+    stake_dist: StakeDist = "uniform"
+    pareto_shape: float = 1.16          # Pareto (Lomax) tail index; ~80/20 by default
+    uniform_random: bool = False        # if True, draw i.i.d. uniform stakes; else equal
+    total_stake: float = 1.0e9          # FIXED across distributions for comparability
+
+    # --- network latency (slots) ---
+    latency: int = 0                    # L: block visible to others at t + L
+    latency_stochastic: bool = False    # if True, L is the mean of a stochastic model
+
+    # --- uncle references ---
+    uncle_window: int = constants.W_DEFAULT   # W
+    max_uncles: int = 0                       # U (0 = baseline, no uncles)
+    uncle_strategy: UncleStrategy = "oldest"
+    uncle_random_p: float = 0.5               # coin-flip inclusion prob (random strategy)
+
+    # --- consensus / TSI ---
+    f: float = constants.F
+    beta: float = constants.BETA_DEFAULT
+    k: int = 64                          # scaled by default; full scale = 2160
+    genesis_d_factor: float = 0.5        # genesis D = factor * true total stake
+    epochs: int = 40
+    per_node_dest: bool = False          # Phase-2 hook: per-node D_est (unused in reduced model)
+
+    # --- bookkeeping ---
+    replicate: int = 0
+    root_seed: int = 12345
+
+    # derived geometry -------------------------------------------------------
+    @property
+    def epoch_len(self) -> int:
+        return constants.epoch_len(self.k, self.f)
+
+    @property
+    def period_T(self) -> int:
+        return constants.period_T(self.k, self.f)
+
+    def key(self) -> tuple:
+        """Hashable identity used to seed the RNG deterministically."""
+        return (
+            self.n_nodes, self.stake_dist, self.pareto_shape, self.uniform_random,
+            self.total_stake, self.latency, self.latency_stochastic, self.uncle_window,
+            self.max_uncles, self.uncle_strategy, self.uncle_random_p, self.f, self.beta,
+            self.k, self.genesis_d_factor, self.epochs, self.per_node_dest, self.replicate,
+        )
+
+
+@dataclass
+class SweepConfig:
+    """A cartesian grid of runs plus replicates, all sharing ``base`` settings."""
+
+    n_nodes: list[int] = field(default_factory=lambda: [1000])
+    stake_dist: list[StakeDist] = field(default_factory=lambda: ["uniform", "pareto"])
+    latency: list[int] = field(default_factory=lambda: [0, 1, 2, 4, 8])
+    max_uncles: list[int] = field(default_factory=lambda: [0, 1, 2, 3, 4])
+    uncle_strategy: list[UncleStrategy] = field(default_factory=lambda: ["oldest", "random"])
+    f: list[float] = field(default_factory=lambda: [constants.F])
+    replicates: int = 8
+    base: dict[str, Any] = field(default_factory=dict)
+
+    def expand(self) -> list[SimConfig]:
+        """Materialise every ``SimConfig`` in the grid × replicates."""
+        base = SimConfig(**self.base)
+        cells: list[SimConfig] = []
+        axes = itertools.product(
+            self.n_nodes, self.stake_dist, self.latency, self.max_uncles,
+            self.uncle_strategy, self.f,
+        )
+        for n, dist, lat, u, strat, fval in axes:
+            # U=0 is strategy-independent; keep only one strategy to avoid duplicate work.
+            if u == 0 and strat != self.uncle_strategy[0]:
+                continue
+            for rep in range(self.replicates):
+                cells.append(
+                    replace(
+                        base,
+                        n_nodes=n,
+                        stake_dist=dist,
+                        latency=lat,
+                        max_uncles=u,
+                        uncle_strategy=strat,
+                        f=fval,
+                        replicate=rep,
+                    )
+                )
+        return cells
+
+    @classmethod
+    def from_dict(cls, d: dict[str, Any]) -> SweepConfig:
+        d = dict(d)
+        base = d.pop("base", {})
+        known = {
+            "n_nodes", "stake_dist", "latency", "max_uncles",
+            "uncle_strategy", "f", "replicates",
+        }
+        kwargs = {k: v for k, v in d.items() if k in known}
+        return cls(base=base, **kwargs)
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/constants.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/constants.py
new file mode 100644
index 0000000..baba55b
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/constants.py
@@ -0,0 +1,39 @@
+"""Protocol constants and epoch/window geometry.
+
+All slot geometry derives from the pair ``(k, f)`` so a scaled-down ``k`` (used for
+parameter sweeps) automatically shrinks the epoch and measurement window. See
+``cryptarchia-v1-protocol.md`` and ``cryptarchia-total-stake-inference.md``.
+"""
+
+from __future__ import annotations
+
+# --- True protocol values (full scale) -------------------------------------
+K_TRUE = 2160          # security parameter (blocks)
+F = 1.0 / 30.0         # slot activation coefficient (fixed; never tuned)
+W_DEFAULT = 300        # uncle reference window w_u (slots)
+BETA_DEFAULT = 1.0     # TSI learning rate
+SLOT_SECONDS = 1       # slot length
+
+
+def floor_k_over_f(k: int, f: float = F) -> int:
+    """``floor(k / f)`` — the base quantum of the epoch schedule."""
+    return int(k / f)
+
+
+def epoch_len(k: int, f: float = F) -> int:
+    """Epoch length in slots: ``10 * floor(k/f)``."""
+    return 10 * floor_k_over_f(k, f)
+
+
+def period_T(k: int, f: float = F) -> int:
+    """TSI measurement window length ``T`` in slots: ``6 * floor(k/f)``.
+
+    This is the first ``6*floor(k/f)`` slots of the (previous) epoch over which the
+    block density is measured.
+    """
+    return 6 * floor_k_over_f(k, f)
+
+
+def expected_blocks_in_window(k: int, f: float = F) -> float:
+    """Expected honest-chain block count in the measurement window at equilibrium."""
+    return period_T(k, f) * f
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/engine.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/engine.py
new file mode 100644
index 0000000..54531db
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/engine.py
@@ -0,0 +1,34 @@
+"""Multi-epoch trajectory driver for a single config."""
+
+from __future__ import annotations
+
+from typing import Any
+
+from . import tsi
+from .config import SimConfig
+from .epoch import simulate_epoch
+from .metrics import metric_row
+from .rng import rng_for
+from .stake import make_stake
+
+
+def run_trajectory(config: SimConfig) -> list[dict[str, Any]]:
+    """Run ``config.epochs`` epochs of TSI, returning one metric row per epoch.
+
+    Stake is drawn once (it is fixed; only ``D_est`` evolves). ``D_est`` starts at the
+    hardcoded genesis value ``genesis_d_factor * D_true`` and is updated each epoch from
+    the measured density.
+    """
+    rng = rng_for(config)
+    stake = make_stake(config, rng)
+    d_true = float(stake.sum())
+    d_est = config.genesis_d_factor * d_true
+    T = config.period_T
+
+    rows: list[dict[str, Any]] = []
+    for epoch in range(config.epochs):
+        er = simulate_epoch(config, stake, d_est, rng)
+        d_next = tsi.update_D(d_est, er.m, T, config.f, config.beta)
+        rows.append(metric_row(config, epoch, d_est, d_next, d_true, er))
+        d_est = d_next
+    return rows
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/epoch.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/epoch.py
new file mode 100644
index 0000000..344d416
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/epoch.py
@@ -0,0 +1,62 @@
+"""Single-epoch simulation: lottery -> block tree -> uncles -> density counting."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from . import lottery, tsi
+from .blocktree import build_tree
+from .config import SimConfig
+from .latency import make_latency
+
+
+@dataclass
+class EpochResult:
+    m: int                 # TSI block count in window
+    q: float               # honest active-slot fraction (window)
+    q_eff: float           # uncle-recovered active-slot fraction (window)
+    n_active: int          # active slots in window
+    n_honest: int          # honest slots in window
+    n_recovered: int       # orphan slots recovered by uncles in window
+    total_winners_window: int   # total lottery wins in window (incl. multi-winner)
+    n_blocks: int          # real blocks produced this epoch
+    n_canonical: int       # canonical chain length
+    n_orphans: int         # orphaned blocks
+
+
+def simulate_epoch(
+    config: SimConfig, stake: np.ndarray, d_est: float, rng: np.random.Generator
+) -> EpochResult:
+    f = config.f
+    T = config.period_T
+    p_win = lottery.win_probs(stake, d_est, f)
+    winner_slots, winner_nodes = lottery.sample_wins(p_win, config.epoch_len, rng)
+    active_slots, groups = lottery.group_by_slot(winner_slots, winner_nodes)
+
+    latency = make_latency(config)
+    tree = build_tree(active_slots, groups, latency, rng)
+    canonical = tree.canonical_chain()
+    from .uncles import annotate_uncles
+
+    annotate_uncles(tree, canonical, config, rng)
+
+    m = tsi.density_m(tree, canonical, T)
+    ref = tsi.referenced_uncle_ids(tree, canonical)
+    ss = tsi.slot_stats(tree, canonical, ref, active_slots, T)
+
+    total_winners_window = int((winner_slots < T).sum())
+    n_real = tree.n_blocks - 1
+    return EpochResult(
+        m=m,
+        q=ss.q,
+        q_eff=ss.q_eff,
+        n_active=ss.n_active,
+        n_honest=ss.n_honest,
+        n_recovered=ss.n_recovered,
+        total_winners_window=total_winners_window,
+        n_blocks=n_real,
+        n_canonical=len(canonical),
+        n_orphans=n_real - len(canonical),
+    )
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/latency.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/latency.py
new file mode 100644
index 0000000..213c7c1
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/latency.py
@@ -0,0 +1,54 @@
+"""Network latency models.
+
+Latency ``L`` is the number of slots between a block being produced and it becoming
+visible to the rest of the network. ``L`` is deliberately named to avoid clashing with
+``D`` (the stake estimate). A leader at slot ``t`` can only build on blocks whose
+``visible_at <= t`` (its own block is visible to itself immediately), which is what
+produces latency-induced forks.
+"""
+
+from __future__ import annotations
+
+from typing import Protocol
+
+import numpy as np
+
+
+class LatencyModel(Protocol):
+    def visible_at(self, produced_slot: int, rng: np.random.Generator) -> int:
+        """Slot at which a block produced at ``produced_slot`` becomes visible to others."""
+        ...
+
+
+class FixedSlotLatency:
+    """Deterministic integer-slot latency: visible to all others at ``t + L``."""
+
+    def __init__(self, latency: int) -> None:
+        self.latency = int(latency)
+
+    def visible_at(self, produced_slot: int, rng: np.random.Generator) -> int:
+        return produced_slot + self.latency
+
+
+class RealisticLatency:
+    """Stochastic latency with mean ``L`` slots (optional sensitivity model).
+
+    Rounds an exponential draw (mean ``L``) up to whole slots. A stand-in for the
+    reference notebook's blend/broadcast delay model; not used by the primary sweep.
+    """
+
+    def __init__(self, mean_latency: float) -> None:
+        self.mean_latency = float(mean_latency)
+
+    def visible_at(self, produced_slot: int, rng: np.random.Generator) -> int:
+        if self.mean_latency <= 0:
+            return produced_slot
+        draw = rng.exponential(self.mean_latency)
+        return produced_slot + int(np.ceil(draw))
+
+
+def make_latency(config) -> LatencyModel:  # noqa: ANN001 - avoid import cycle with config
+    """Build the latency model for a config."""
+    if config.latency_stochastic:
+        return RealisticLatency(config.latency)
+    return FixedSlotLatency(config.latency)
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/lottery.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/lottery.py
new file mode 100644
index 0000000..da7ebe0
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/lottery.py
@@ -0,0 +1,63 @@
+"""Stake-weighted slot lottery.
+
+Per node ``i`` and slot, an independent Bernoulli win with probability
+``phi_f(alpha_i) = 1 - (1 - f)^alpha_i`` where ``alpha_i = w_i / D_est`` is the node's
+relative stake against its inferred total active stake. Multiple winners in a slot are
+possible (a guaranteed fork). ``phi`` is taken verbatim from the reference notebook.
+
+Winners are returned as *sparse coordinates* — flat ``(winner_slots, winner_nodes)``
+arrays sorted by slot — so downstream code iterates only over the ``O(f * n_slots)``
+winners, never over every slot.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+
+def phi(f: float, alpha: np.ndarray | float) -> np.ndarray | float:
+    """Leader-lottery win probability ``1 - (1 - f)^alpha``."""
+    return 1.0 - (1.0 - f) ** alpha
+
+
+def win_probs(stake: np.ndarray, d_est: float, f: float) -> np.ndarray:
+    """Per-node win probability ``phi_f(w_i / D_est)``."""
+    return phi(f, stake / d_est)
+
+
+def sample_wins(
+    p_win: np.ndarray, n_slots: int, rng: np.random.Generator, chunk: int = 8192
+) -> tuple[np.ndarray, np.ndarray]:
+    """Sample lottery wins over ``n_slots`` slots.
+
+    Returns ``(winner_slots, winner_nodes)``: parallel int arrays with one entry per
+    (slot, winning node), sorted by slot ascending. Draws are chunked over slots to bound
+    peak memory to ``n_nodes * chunk`` bools.
+    """
+    n = p_win.shape[0]
+    p_col = p_win[:, None]
+    slot_parts: list[np.ndarray] = []
+    node_parts: list[np.ndarray] = []
+    for start in range(0, n_slots, chunk):
+        width = min(chunk, n_slots - start)
+        hits = rng.random((n, width)) < p_col
+        node_idx, slot_idx = np.nonzero(hits)
+        slot_parts.append(slot_idx.astype(np.int64) + start)
+        node_parts.append(node_idx.astype(np.int64))
+    if not slot_parts:  # pragma: no cover - n_slots == 0
+        return np.empty(0, np.int64), np.empty(0, np.int64)
+    winner_slots = np.concatenate(slot_parts)
+    winner_nodes = np.concatenate(node_parts)
+    order = np.argsort(winner_slots, kind="stable")
+    return winner_slots[order], winner_nodes[order]
+
+
+def group_by_slot(
+    winner_slots: np.ndarray, winner_nodes: np.ndarray
+) -> tuple[np.ndarray, list[np.ndarray]]:
+    """Group sorted winner coordinates into ``(active_slots, winners_per_active_slot)``."""
+    if winner_slots.size == 0:
+        return np.empty(0, np.int64), []
+    active_slots, starts = np.unique(winner_slots, return_index=True)
+    groups = np.split(winner_nodes, starts[1:])
+    return active_slots, groups
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/metrics.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/metrics.py
new file mode 100644
index 0000000..0b84112
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/metrics.py
@@ -0,0 +1,70 @@
+"""Per-epoch metric rows and equilibrium summaries."""
+
+from __future__ import annotations
+
+from typing import Any
+
+import numpy as np
+
+from .config import SimConfig
+from .epoch import EpochResult
+
+# Config fields recorded on every row for grouping/plotting.
+_CONFIG_FIELDS = (
+    "n_nodes", "stake_dist", "pareto_shape", "latency", "uncle_window",
+    "max_uncles", "uncle_strategy", "uncle_random_p", "f", "beta", "k",
+    "genesis_d_factor", "epochs", "replicate",
+)
+
+
+def metric_row(
+    config: SimConfig,
+    epoch: int,
+    d_in: float,
+    d_out: float,
+    d_true: float,
+    er: EpochResult,
+) -> dict[str, Any]:
+    row: dict[str, Any] = {field: getattr(config, field) for field in _CONFIG_FIELDS}
+    row.update(
+        epoch=epoch,
+        d_in=d_in,
+        d_out=d_out,
+        d_true=d_true,
+        ratio=d_out / d_true,
+        m=er.m,
+        measured_density=er.m / config.period_T,
+        q=er.q,
+        q_eff=er.q_eff,
+        n_active=er.n_active,
+        n_honest=er.n_honest,
+        n_recovered=er.n_recovered,
+        total_winners_window=er.total_winners_window,
+        n_blocks=er.n_blocks,
+        n_canonical=er.n_canonical,
+        n_orphans=er.n_orphans,
+        orphan_rate=er.n_orphans / er.n_blocks if er.n_blocks else float("nan"),
+    )
+    return row
+
+
+def equilibrium_stats(ratios: np.ndarray, burn_in: int) -> dict[str, float]:
+    """Mean/variance of the stake ratio after ``burn_in`` epochs."""
+    tail = ratios[burn_in:]
+    if tail.size == 0:
+        tail = ratios[-1:]
+    return {
+        "mean_ratio": float(np.mean(tail)),
+        "var_ratio": float(np.var(tail)),
+        "std_ratio": float(np.std(tail)),
+    }
+
+
+def epochs_to_within(ratios: np.ndarray, target: float, eps: float) -> int:
+    """First epoch index after which ``|ratio - target| <= eps`` holds for the rest."""
+    within = np.abs(ratios - target) <= eps
+    n = within.size
+    for i in range(n):
+        if within[i:].all():
+            return i
+    return n  # never converged within the run
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/__init__.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/__init__.py
new file mode 100644
index 0000000..527b68e
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/__init__.py
@@ -0,0 +1 @@
+"""Academic-quality figure generation."""
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/figures.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/figures.py
new file mode 100644
index 0000000..99bf1aa
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/figures.py
@@ -0,0 +1,241 @@
+"""Figure builders. Each takes the results DataFrame and returns a Matplotlib Figure.
+
+Rows are per (config, epoch). We summarise each config to its *equilibrium* by averaging
+the tail epochs (after a burn-in), keeping replicates so we can draw percentile bands.
+"""
+
+from __future__ import annotations
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+from . import style
+
+CONFIG_COLS = ["n_nodes", "stake_dist", "latency", "max_uncles", "uncle_strategy", "k"]
+
+
+def equilibrium(df: pd.DataFrame, burn_frac: float = 0.5) -> pd.DataFrame:
+    """Per-(config, replicate) equilibrium means over the tail epochs."""
+    cutoff = df["epochs"] * burn_frac
+    tail = df[df["epoch"] >= cutoff]
+    keys = [*CONFIG_COLS, "replicate"]
+    return (
+        tail.groupby(keys, as_index=False)
+        .agg(
+            ratio=("ratio", "mean"),
+            q=("q", "mean"),
+            q_eff=("q_eff", "mean"),
+            var_ratio=("ratio", "var"),
+            orphan_rate=("orphan_rate", "mean"),
+            n_active=("n_active", "mean"),
+        )
+    )
+
+
+def _series_over_replicates(eq: pd.DataFrame, xcol: str, ycol: str, xvals) -> np.ndarray:
+    """Build a ``(n_replicates, len(xvals))`` matrix of ``ycol`` for band plots."""
+    reps = sorted(eq["replicate"].unique())
+    mat = np.full((len(reps), len(xvals)), np.nan)
+    for i, r in enumerate(reps):
+        sub = eq[eq["replicate"] == r].set_index(xcol)[ycol]
+        for j, x in enumerate(xvals):
+            if x in sub.index:
+                mat[i, j] = sub.loc[x]
+    return mat
+
+
+def _provenance(df: pd.DataFrame) -> str:
+    k = [int(x) for x in sorted(df["k"].unique())]
+    n = [int(x) for x in sorted(df["n_nodes"].unique())]
+    reps = int(df["replicate"].nunique())
+    return f"tsi-sim  |  k={k}  N={n}  reps={reps}  f={float(df['f'].iloc[0]):.4g}"
+
+
+# --- Figure 1: accuracy vs U per latency ------------------------------------
+def accuracy_vs_u(df: pd.DataFrame, stake_dist: str, strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["n_nodes"] == df["n_nodes"].max())]
+    eq = eq[(eq["uncle_strategy"] == strategy) | (eq["max_uncles"] == 0)]
+    latencies = sorted(eq["latency"].unique())
+    uvals = sorted(eq["max_uncles"].unique())
+
+    fig, ax = plt.subplots()
+    for i, lat in enumerate(latencies):
+        sub = eq[eq["latency"] == lat]
+        mat = _series_over_replicates(sub, "max_uncles", "ratio", uvals)
+        style.band_plot(ax, uvals, mat, color=style.color_for(i), label=f"L={lat}")
+    ax.axhline(1.0, color="0.4", lw=1.0, ls="--", zorder=0)
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel(r"inferred / true stake  $\langle \hat D / D_{\mathrm{true}} \rangle$")
+    ax.set_title(f"TSI accuracy vs uncle cap ({stake_dist} stake)")
+    ax.set_xticks(uvals)
+    ax.legend(title="latency (slots)", ncol=2)
+    return fig
+
+
+# --- Figure 2: q_eff vs U per latency ---------------------------------------
+def qeff_vs_u(df: pd.DataFrame, stake_dist: str, strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["n_nodes"] == df["n_nodes"].max())]
+    eq = eq[(eq["uncle_strategy"] == strategy) | (eq["max_uncles"] == 0)]
+    latencies = sorted(eq["latency"].unique())
+    uvals = sorted(eq["max_uncles"].unique())
+
+    fig, ax = plt.subplots()
+    for i, lat in enumerate(latencies):
+        sub = eq[eq["latency"] == lat]
+        mat = _series_over_replicates(sub, "max_uncles", "q_eff", uvals)
+        style.band_plot(ax, uvals, mat, color=style.color_for(i), label=f"L={lat}")
+        base_q = sub[sub["max_uncles"] == 0]["q"].mean()
+        if not np.isnan(base_q):
+            ax.axhline(base_q, color=style.color_for(i), lw=0.8, ls=":", alpha=0.6)
+    ax.axhline(1.0, color="0.4", lw=1.0, ls="--", zorder=0)
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel(r"effective utilisation  $q_{\mathrm{eff}}$")
+    ax.set_title(f"Active-slot recovery vs uncle cap ({stake_dist} stake)")
+    ax.set_xticks(uvals)
+    ax.legend(title="latency (slots)", ncol=2)
+    return fig
+
+
+# --- Figure 3: convergence (ratio vs epoch) per U ---------------------------
+def convergence(df: pd.DataFrame, stake_dist: str, latency: int,
+                strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    sub = df[(df["stake_dist"] == stake_dist) & (df["latency"] == latency)
+             & (df["n_nodes"] == df["n_nodes"].max())]
+    sub = sub[(sub["uncle_strategy"] == strategy) | (sub["max_uncles"] == 0)]
+    uvals = sorted(sub["max_uncles"].unique())
+    epochs = sorted(sub["epoch"].unique())
+
+    fig, ax = plt.subplots()
+    for i, u in enumerate(uvals):
+        su = sub[sub["max_uncles"] == u]
+        mat = np.full((su["replicate"].nunique(), len(epochs)), np.nan)
+        for ri, r in enumerate(sorted(su["replicate"].unique())):
+            s = su[su["replicate"] == r].set_index("epoch")["ratio"]
+            for j, e in enumerate(epochs):
+                if e in s.index:
+                    mat[ri, j] = s.loc[e]
+        style.band_plot(ax, epochs, mat, color=style.color_for(i), label=f"U={u}")
+    ax.axhline(1.0, color="0.4", lw=1.0, ls="--", zorder=0)
+    ax.set_xlabel("epoch")
+    ax.set_ylabel(r"$\hat D / D_{\mathrm{true}}$")
+    ax.set_title(f"Convergence ({stake_dist}, latency L={latency})")
+    ax.legend(title="uncle cap", ncol=2)
+    return fig
+
+
+# --- Figure 4: accuracy heatmap over latency x U ----------------------------
+def heatmap_accuracy(df: pd.DataFrame, stake_dist: str, n_nodes: int,
+                     strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["n_nodes"] == n_nodes)]
+    eq = eq[(eq["uncle_strategy"] == strategy) | (eq["max_uncles"] == 0)]
+    piv = eq.groupby(["latency", "max_uncles"])["ratio"].mean().unstack("max_uncles")
+    latencies = piv.index.to_numpy()
+    uvals = piv.columns.to_numpy()
+    data = piv.to_numpy()
+
+    fig, ax = plt.subplots()
+    vmax = np.nanmax(np.abs(data - 1.0))
+    im = ax.imshow(data, origin="lower", aspect="auto", cmap=style.DIVERGING_CMAP,
+                   vmin=1 - vmax, vmax=1 + vmax)
+    ax.set_xticks(range(len(uvals)), uvals)
+    ax.set_yticks(range(len(latencies)), latencies)
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel("network latency  $L$ (slots)")
+    ax.set_title(f"Accuracy $\\hat D/D_{{\\mathrm{{true}}}}$ ({stake_dist}, N={n_nodes})")
+    for yi in range(len(latencies)):
+        for xi in range(len(uvals)):
+            v = data[yi, xi]
+            if not np.isnan(v):
+                safe = 0.98 <= v <= 1.02
+                ax.text(xi, yi, f"{v:.2f}", ha="center", va="center", fontsize=7,
+                        color="black", fontweight="bold" if safe else "normal")
+    fig.colorbar(im, ax=ax, label=r"$\hat D / D_{\mathrm{true}}$")
+    ax.grid(False)
+    return fig
+
+
+# --- Figure 5: distribution comparison --------------------------------------
+def dist_comparison(df: pd.DataFrame, latency: int, strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["latency"] == latency) & (eq["n_nodes"] == df["n_nodes"].max())]
+    eq = eq[(eq["uncle_strategy"] == strategy) | (eq["max_uncles"] == 0)]
+    uvals = sorted(eq["max_uncles"].unique())
+    fig, ax = plt.subplots()
+    for i, dist in enumerate(sorted(eq["stake_dist"].unique())):
+        sub = eq[eq["stake_dist"] == dist]
+        mat = _series_over_replicates(sub, "max_uncles", "ratio", uvals)
+        style.band_plot(ax, uvals, mat, color=style.color_for(i), label=dist)
+    ax.axhline(1.0, color="0.4", lw=1.0, ls="--", zorder=0)
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel(r"$\hat D / D_{\mathrm{true}}$")
+    ax.set_title(f"Stake-distribution comparison (latency L={latency})")
+    ax.set_xticks(uvals)
+    ax.legend(title="distribution")
+    return fig
+
+
+# --- Figure 6: variance vs U ------------------------------------------------
+def variance_vs_u(df: pd.DataFrame, stake_dist: str, latency: int,
+                  strategy: str = "oldest") -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["latency"] == latency)
+            & (eq["n_nodes"] == df["n_nodes"].max())]
+    eq = eq[(eq["uncle_strategy"] == strategy) | (eq["max_uncles"] == 0)]
+    uvals = sorted(eq["max_uncles"].unique())
+    var_by_u = eq.groupby("max_uncles")["ratio"].var()
+    fig, ax = plt.subplots()
+    ax.plot(uvals, [var_by_u.get(u, np.nan) for u in uvals], "o-",
+            color=style.color_for(0), label="empirical (across replicates)")
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel(r"$\mathrm{Var}[\hat D / D_{\mathrm{true}}]$")
+    ax.set_title(f"Estimate variance vs uncle cap ({stake_dist}, L={latency})")
+    ax.set_xticks(uvals)
+    ax.legend()
+    return fig
+
+
+# --- Figure 7: strategy comparison ------------------------------------------
+def strategy_comparison(df: pd.DataFrame, stake_dist: str, latency: int) -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["latency"] == latency)
+            & (eq["n_nodes"] == df["n_nodes"].max()) & (eq["max_uncles"] > 0)]
+    uvals = sorted(eq["max_uncles"].unique())
+    fig, ax = plt.subplots()
+    for i, strat in enumerate(sorted(eq["uncle_strategy"].unique())):
+        sub = eq[eq["uncle_strategy"] == strat]
+        mat = _series_over_replicates(sub, "max_uncles", "ratio", uvals)
+        style.band_plot(ax, uvals, mat, color=style.color_for(i), label=strat)
+    ax.axhline(1.0, color="0.4", lw=1.0, ls="--", zorder=0)
+    ax.set_xlabel("max uncles per block  $U$")
+    ax.set_ylabel(r"$\hat D / D_{\mathrm{true}}$")
+    ax.set_title(f"Uncle-selection strategy ({stake_dist}, L={latency})")
+    ax.set_xticks(uvals)
+    ax.legend(title="strategy")
+    return fig
+
+
+# --- Figure 8: diagnostics (orphan rate vs latency) -------------------------
+def orphan_diagnostics(df: pd.DataFrame, stake_dist: str) -> plt.Figure:
+    style.apply_style()
+    eq = equilibrium(df)
+    eq = eq[(eq["stake_dist"] == stake_dist) & (eq["max_uncles"] == 0)
+            & (eq["n_nodes"] == df["n_nodes"].max())]
+    g = eq.groupby("latency")["orphan_rate"].agg(["mean", "std"])
+    fig, ax = plt.subplots()
+    ax.errorbar(g.index, g["mean"], yerr=g["std"], marker="o", color=style.color_for(1),
+                capsize=3)
+    ax.set_xlabel("network latency  $L$ (slots)")
+    ax.set_ylabel("orphan rate (orphans / blocks)")
+    ax.set_title(f"Fork/orphan rate vs latency ({stake_dist} stake)")
+    return fig
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/make_figures.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/make_figures.py
new file mode 100644
index 0000000..b37a193
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/make_figures.py
@@ -0,0 +1,69 @@
+"""Render academic figures from a results frame (importable + CLI)."""
+
+from __future__ import annotations
+
+import argparse
+import datetime as _dt
+from pathlib import Path
+
+import pandas as pd
+
+from . import figures as F
+from . import style
+
+
+def timestamped_figdir(outdir: str | Path, label: str) -> Path:
+    """A fresh ``<outdir>/<YYYY-MM-DD_HHMMSS>_<label>`` folder (never overwrites)."""
+    ts = _dt.datetime.now().strftime("%Y-%m-%d_%H%M%S")
+    base = Path(outdir) / f"{ts}_{label}"
+    d, i = base, 2
+    while d.exists():
+        d = base.with_name(f"{base.name}_{i}")
+        i += 1
+    return d
+
+
+def render(df: pd.DataFrame, out: str | Path) -> int:
+    """Render every applicable figure from ``df`` into directory ``out``; return the count."""
+    out = Path(out)
+    prov = F._provenance(df)
+    dists = sorted(df["stake_dist"].unique())
+    latencies = sorted(df["latency"].unique())
+    lat_focus = latencies[len(latencies) // 2] if latencies else 0
+    written: list[Path] = []
+    for dist in dists:
+        written += style.save(F.accuracy_vs_u(df, dist), out / f"01_accuracy_vs_u_{dist}", prov)
+        written += style.save(F.qeff_vs_u(df, dist), out / f"02_qeff_vs_u_{dist}", prov)
+        written += style.save(F.convergence(df, dist, lat_focus),
+                              out / f"03_convergence_{dist}_L{lat_focus}", prov)
+        written += style.save(F.orphan_diagnostics(df, dist), out / f"08_orphans_{dist}", prov)
+        for n in sorted(df["n_nodes"].unique()):
+            written += style.save(F.heatmap_accuracy(df, dist, int(n)),
+                                  out / f"04_heatmap_{dist}_N{int(n)}", prov)
+        written += style.save(F.variance_vs_u(df, dist, lat_focus),
+                              out / f"06_variance_vs_u_{dist}_L{lat_focus}", prov)
+        if df[(df["stake_dist"] == dist) & (df["max_uncles"] > 0)]["uncle_strategy"].nunique() > 1:
+            written += style.save(F.strategy_comparison(df, dist, lat_focus),
+                                  out / f"07_strategy_{dist}_L{lat_focus}", prov)
+    written += style.save(
+        F.dist_comparison(df, lat_focus), out / f"05_dist_comparison_L{lat_focus}", prov
+    )
+    return len(written)
+
+
+def main(argv: list[str] | None = None) -> None:
+    ap = argparse.ArgumentParser(description="Render TSI figures")
+    ap.add_argument("--results", required=True)
+    ap.add_argument("--outdir", default="figures",
+                    help="parent dir; a dated sub-folder is created (never overwrites)")
+    ap.add_argument("--label", default=None, help="run label (default: results file stem)")
+    ap.add_argument("--out", default=None, help="explicit output dir (skips the dated folder)")
+    args = ap.parse_args(argv)
+    label = args.label or Path(args.results).stem
+    out = Path(args.out) if args.out else timestamped_figdir(args.outdir, label)
+    n = render(pd.read_parquet(args.results), out)
+    print(f"wrote {n} figures -> {out}")
+
+
+if __name__ == "__main__":  # pragma: no cover
+    main()
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/style.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/style.py
new file mode 100644
index 0000000..f012417
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/plotting/style.py
@@ -0,0 +1,95 @@
+"""Shared academic matplotlib theme, palette, and helpers.
+
+Palette: Okabe-Ito — the standard colorblind-safe qualitative set — for categorical
+series (the ``U`` lines); ``cividis`` (perceptually uniform, CVD-safe) for heatmaps.
+Figures are saved as vector PDF (primary) plus 300-dpi PNG preview.
+"""
+
+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 map centered at ratio = 1 (accuracy heatmaps)
+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 PDF + 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:
+        # leave headroom so the footnote never collides with the x-axis label
+        fig.subplots_adjust(bottom=0.16)
+        fig.text(0.99, 0.005, provenance, fontsize=6, alpha=0.5, va="bottom", ha="right")
+    paths = []
+    for ext in ("pdf", "png"):
+        p = out_stem.with_suffix(f".{ext}")
+        fig.savefig(p)
+        paths.append(p)
+    plt.close(fig)
+    return paths
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/rng.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/rng.py
new file mode 100644
index 0000000..d106dfb
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/rng.py
@@ -0,0 +1,26 @@
+"""Deterministic, order-independent RNG derivation.
+
+Each ``SimConfig`` maps to an independent NumPy ``Generator`` seeded from a hash of its
+identity plus the global root seed. This guarantees a config yields the same random
+stream regardless of the order in which parallel workers run it.
+"""
+
+from __future__ import annotations
+
+import hashlib
+
+import numpy as np
+
+from .config import SimConfig
+
+
+def _entropy(config: SimConfig) -> int:
+    payload = repr((config.root_seed, config.key())).encode()
+    digest = hashlib.blake2b(payload, digest_size=16).digest()
+    return int.from_bytes(digest, "big")
+
+
+def rng_for(config: SimConfig) -> np.random.Generator:
+    """Return the reproducible ``Generator`` for this exact config+replicate."""
+    seed = np.random.SeedSequence(_entropy(config))
+    return np.random.default_rng(seed)
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/stake.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/stake.py
new file mode 100644
index 0000000..a9e3a32
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/stake.py
@@ -0,0 +1,32 @@
+"""Stake distribution generation.
+
+Total stake is held FIXED across distributions (via renormalisation) so that accuracy
+comparisons between ``uniform`` and ``pareto`` isolate the *shape* effect on the lottery
+(winner multiplicity / forking), not a difference in aggregate stake.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .config import SimConfig
+
+
+def make_stake(config: SimConfig, rng: np.random.Generator) -> np.ndarray:
+    """Return an ``(n_nodes,)`` non-negative stake vector summing to ``total_stake``."""
+    n = config.n_nodes
+    if config.stake_dist == "uniform":
+        if config.uniform_random:
+            w = rng.random(n)
+        else:
+            w = np.ones(n)
+    elif config.stake_dist == "pareto":
+        # numpy.pareto draws Lomax = Pareto(shape) - 1, heavy-tailed for small shape.
+        w = rng.pareto(config.pareto_shape, n) + 1.0
+    else:  # pragma: no cover - guarded by Literal typing
+        raise ValueError(f"unknown stake_dist: {config.stake_dist}")
+
+    total = w.sum()
+    if total <= 0:  # pragma: no cover - degenerate
+        raise ValueError("stake vector summed to zero")
+    return w * (config.total_stake / total)
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/sweep.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/sweep.py
new file mode 100644
index 0000000..c48eae0
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/sweep.py
@@ -0,0 +1,85 @@
+"""Parameter-sweep expansion, parallel execution, and result persistence."""
+
+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 SweepConfig
+from .engine import run_trajectory
+
+
+def run_sweep(sweep: SweepConfig, n_jobs: int = -1, progress: bool = True) -> pd.DataFrame:
+    """Expand the grid, run every config (embarrassingly parallel), return one big frame."""
+    configs = sweep.expand()
+    runner = Parallel(n_jobs=n_jobs, return_as="generator")(
+        delayed(run_trajectory)(c) for c in configs
+    )
+    if progress:
+        runner = tqdm(runner, total=len(configs), desc="configs")
+    rows: list[dict] = []
+    for traj in runner:
+        rows.extend(traj)
+    return pd.DataFrame(rows)
+
+
+def persist(df: pd.DataFrame, path: str | Path) -> None:
+    path = Path(path)
+    path.parent.mkdir(parents=True, exist_ok=True)
+    if path.suffix == ".parquet":
+        df.to_parquet(path, index=False)
+    else:
+        df.to_csv(path, index=False)
+
+
+def load_sweep_yaml(path: str | Path) -> SweepConfig:
+    with open(path) as fh:
+        data = yaml.safe_load(fh)
+    return SweepConfig.from_dict(data)
+
+
+def new_run_dir(outdir: str | Path, label: str) -> Path:
+    """Create a fresh timestamped run directory so runs never overwrite each other."""
+    ts = _dt.datetime.now().strftime("%Y-%m-%d_%H%M%S")
+    base = Path(outdir) / f"{ts}_{label}"
+    run_dir, i = base, 2
+    while run_dir.exists():
+        run_dir = base.with_name(f"{base.name}_{i}")
+        i += 1
+    run_dir.mkdir(parents=True)
+    return run_dir
+
+
+def main(argv: list[str] | None = None) -> None:
+    parser = argparse.ArgumentParser(description="Run a TSI parameter sweep + figures")
+    parser.add_argument("--config", required=True, help="sweep YAML")
+    parser.add_argument("--outdir", default="runs",
+                        help="parent dir; a dated sub-folder is created per run")
+    parser.add_argument("--label", default=None, help="run label (default: config name)")
+    parser.add_argument("--n-jobs", type=int, default=-1)
+    parser.add_argument("--no-figures", action="store_true", help="skip auto figure generation")
+    args = parser.parse_args(argv)
+
+    label = args.label or Path(args.config).stem
+    run_dir = new_run_dir(args.outdir, label)
+    sweep = load_sweep_yaml(args.config)
+    df = run_sweep(sweep, n_jobs=args.n_jobs)
+    results_path = run_dir / "results.parquet"
+    persist(df, results_path)
+    n_cfg = df[["n_nodes", "stake_dist", "latency", "max_uncles", "uncle_strategy", "replicate"]]
+    print(f"wrote {len(df)} rows ({len(n_cfg.drop_duplicates())} configs) -> {results_path}")
+    if not args.no_figures:
+        from .plotting.make_figures import render
+        n_fig = render(df, run_dir / "figures")
+        print(f"wrote {n_fig} figures -> {run_dir / 'figures'}")
+    print(f"run dir: {run_dir}")
+
+
+if __name__ == "__main__":  # pragma: no cover
+    main()
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/theory.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/theory.py
new file mode 100644
index 0000000..5b05e74
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/theory.py
@@ -0,0 +1,42 @@
+"""Closed-form TSI results from ``analysis-total-stake-inference.md``.
+
+Used both as figure overlays and as ground truth for the verification checks. ``q`` is
+the honest active-slot utilisation; with uncle references, substitute the effective
+``q_eff`` to predict the improved accuracy.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+ArrayLike = np.ndarray | float
+
+
+def expected_ratio(f: float, q: ArrayLike) -> ArrayLike:
+    """Equilibrium ``E[D_inf] / D_true = log(1-f) / log(1 - f/q)`` for ``q in (f, 1]``."""
+    q = np.asarray(q, dtype=float)
+    return np.log(1.0 - f) / np.log(1.0 - f / q)
+
+
+def variance_ratio(f: float, q: ArrayLike, T: int, beta: float = 1.0) -> ArrayLike:
+    """Equilibrium ``Var[D_inf / D_true]``."""
+    q = np.asarray(q, dtype=float)
+    er = expected_ratio(f, q)
+    return (beta / f) ** 2 * (q / T) * er**2 * (1.0 - f) * f
+
+
+def variance_bound(f: float, T: int, beta: float = 1.0) -> float:
+    """Upper bound on ``Var[D_inf / D_true]`` at ``q = 1`` (perfect network)."""
+    return (beta / f) ** 2 / T * (1.0 - f) * f
+
+
+def beta_stability_bound(f: float, q: ArrayLike) -> ArrayLike:
+    """Stability threshold: convergence requires ``beta < 2f/((q-f) log(1/(1-f/q)))``."""
+    q = np.asarray(q, dtype=float)
+    return 2.0 * f / ((q - f) * np.log(1.0 / (1.0 - f / q)))
+
+
+def optimal_beta(f: float, q: ArrayLike) -> ArrayLike:
+    """Convergence-optimal learning rate ``f/((q-f) log(1/(1-f/q)))``."""
+    q = np.asarray(q, dtype=float)
+    return f / ((q - f) * np.log(1.0 / (1.0 - f / q)))
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/tsi.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/tsi.py
new file mode 100644
index 0000000..6442430
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/tsi.py
@@ -0,0 +1,77 @@
+"""Total Stake Inference: density counting and the per-epoch estimate update.
+
+The estimate update counts *blocks* exactly as the spec's ``density_over_slots`` does:
+``m = honest-chain blocks in window + deduplicated referenced uncle blocks (by their own
+slot) in window``. We additionally report slot-based ``q`` / ``q_eff`` (honest and
+uncle-recovered active-slot fractions) for comparison against the closed-form theory,
+which is written in terms of active-slot utilisation. The two differ only at rare
+multi-winner slots.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from .blocktree import BlockTree
+
+
+def referenced_uncle_ids(tree: BlockTree, canonical_ids: list[int]) -> set[int]:
+    """Deduplicated set of uncle ids referenced by the canonical chain."""
+    ref: set[int] = set()
+    for b in canonical_ids:
+        ref.update(tree.uncles[b])
+    return ref
+
+
+def _in_window(slot: int, T: int) -> bool:
+    return 0 <= slot < T
+
+
+def density_m(tree: BlockTree, canonical_ids: list[int], T: int) -> int:
+    """Block count ``m`` for the TSI update (honest blocks + deduped uncles, in window)."""
+    s = tree.slot[canonical_ids]
+    honest = int(((s >= 0) & (s < T)).sum())
+    ref = referenced_uncle_ids(tree, canonical_ids)
+    uncle = sum(1 for u in ref if _in_window(int(tree.slot[u]), T))
+    return honest + uncle
+
+
+def update_D(d_prev: float, m: int, T: int, f: float, beta: float) -> float:
+    """Spec TSI recursion: ``max(1, D_prev * (1 - beta*(f - m/T)/f))``."""
+    measured_density = m / T
+    d_new = d_prev * (1.0 - beta * (f - measured_density) / f)
+    return max(d_new, 1.0)
+
+
+@dataclass
+class SlotStats:
+    n_active: int          # active slots (>=1 winner) in window
+    n_honest: int          # honest-chain-occupied slots in window
+    n_recovered: int       # orphan-only slots recovered via referenced uncles
+    q: float               # n_honest / n_active
+    q_eff: float           # (n_honest + n_recovered) / n_active
+
+
+def slot_stats(
+    tree: BlockTree,
+    canonical_ids: list[int],
+    ref_uncle_ids: set[int],
+    active_slots: np.ndarray,
+    T: int,
+) -> SlotStats:
+    """Slot-based utilisation stats used for theory overlays."""
+    active_in = active_slots[(active_slots >= 0) & (active_slots < T)]
+    n_active = int(active_in.size)
+    honest_slots = {int(tree.slot[b]) for b in canonical_ids if _in_window(int(tree.slot[b]), T)}
+    recovered: set[int] = set()
+    for u in ref_uncle_ids:
+        su = int(tree.slot[u])
+        if _in_window(su, T) and su not in honest_slots:
+            recovered.add(su)
+    n_honest = len(honest_slots)
+    n_rec = len(recovered)
+    q = n_honest / n_active if n_active else float("nan")
+    q_eff = (n_honest + n_rec) / n_active if n_active else float("nan")
+    return SlotStats(n_active=n_active, n_honest=n_honest, n_recovered=n_rec, q=q, q_eff=q_eff)
diff --git a/tools/simulators/tsi/tsi-sim/src/tsi_sim/uncles.py b/tools/simulators/tsi/tsi-sim/src/tsi_sim/uncles.py
new file mode 100644
index 0000000..1daa6f7
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/src/tsi_sim/uncles.py
@@ -0,0 +1,80 @@
+"""Proposer-local uncle selection.
+
+For each canonical block ``B`` (processed oldest-first so ancestors' references are
+known), candidates are orphan (non-canonical) blocks ``U`` with
+``0 < slot_B - slot_U <= W`` that have not already been referenced by an ancestor of
+``B``. Two strategies match the spec: deterministic oldest-first, and random (oldest-first
+order, unbiased coin per candidate, capped at ``U``). Dedup across ancestors is enforced
+by threading a ``referenced`` set down the canonical chain.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .blocktree import BlockTree
+from .config import SimConfig
+
+
+def _orphans_sorted(tree: BlockTree, canonical_ids: list[int]) -> tuple[np.ndarray, np.ndarray]:
+    """Return orphan block ids sorted by (slot, id) and their slots."""
+    canonical = np.zeros(tree.n_blocks, dtype=bool)
+    canonical[canonical_ids] = True
+    all_real = np.arange(1, tree.n_blocks)
+    orphan_ids = all_real[~canonical[1:]]
+    orphan_slots = tree.slot[orphan_ids]
+    order = np.lexsort((orphan_ids, orphan_slots))  # by slot, then id
+    return orphan_ids[order], orphan_slots[order]
+
+
+def annotate_uncles(
+    tree: BlockTree, canonical_ids: list[int], config: SimConfig, rng: np.random.Generator
+) -> None:
+    """Fill ``tree.uncles[B]`` for every canonical block ``B`` per the selection rule."""
+    u_max = config.max_uncles
+    if u_max <= 0:
+        return
+    w = config.uncle_window
+    orphan_ids, orphan_slots = _orphans_sorted(tree, canonical_ids)
+    if orphan_ids.size == 0:
+        return
+
+    referenced: set[int] = set()
+    # oldest canonical block first
+    for b in reversed(canonical_ids):
+        sb = int(tree.slot[b])
+        lo = int(np.searchsorted(orphan_slots, sb - w, side="left"))   # slot_U >= sb - W
+        hi = int(np.searchsorted(orphan_slots, sb, side="left"))       # slot_U <  sb
+        if hi <= lo:
+            continue
+        window_ids = orphan_ids[lo:hi]  # already oldest-first
+        selected = _select(window_ids, referenced, config, rng)
+        if selected:
+            tree.uncles[b] = tuple(selected)
+            referenced.update(selected)
+
+
+def _select(
+    window_ids: np.ndarray, referenced: set[int], config: SimConfig, rng: np.random.Generator
+) -> list[int]:
+    u_max = config.max_uncles
+    out: list[int] = []
+    if config.uncle_strategy == "oldest":
+        for bid in window_ids.tolist():
+            if bid in referenced:
+                continue
+            out.append(bid)
+            if len(out) >= u_max:
+                break
+    elif config.uncle_strategy == "random":
+        p = config.uncle_random_p
+        for bid in window_ids.tolist():
+            if bid in referenced:
+                continue
+            if rng.random() < p:
+                out.append(bid)
+                if len(out) >= u_max:
+                    break
+    else:  # pragma: no cover - guarded by Literal
+        raise ValueError(config.uncle_strategy)
+    return out
diff --git a/tools/simulators/tsi/tsi-sim/tests/test_blocktree.py b/tools/simulators/tsi/tsi-sim/tests/test_blocktree.py
new file mode 100644
index 0000000..367581a
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/tests/test_blocktree.py
@@ -0,0 +1,54 @@
+import numpy as np
+
+from tsi_sim.blocktree import build_tree
+from tsi_sim.latency import FixedSlotLatency
+
+
+def _winners(*groups):
+    return [np.array(g, dtype=np.int64) for g in groups]
+
+
+def test_latency_induces_fork_and_longest_chain():
+    # L=2. slot0: node0; slot1: node1 (can't see block1 yet -> forks on genesis);
+    # slot3: node2 (sees both, builds on the earlier-slot tip -> block1).
+    active = np.array([0, 1, 3], dtype=np.int64)
+    winners = _winners([0], [1], [2])
+    tree = build_tree(active, winners, FixedSlotLatency(2), np.random.default_rng(0))
+
+    assert tree.n_blocks == 4                      # genesis + 3
+    assert tree.height.tolist() == [0, 1, 1, 2]
+    # block 3 built on block 1 (earlier slot wins the height-1 tie)
+    assert tree.parent[3] == 1
+    assert tree.canonical_chain() == [3, 1]        # tip-first
+    # block 2 is the orphan
+    canon = set(tree.canonical_chain())
+    orphans = [b for b in range(1, tree.n_blocks) if b not in canon]
+    assert orphans == [2]
+
+
+def test_same_slot_multiwinner_forks_even_at_zero_latency():
+    active = np.array([0, 1], dtype=np.int64)
+    winners = _winners([0, 1], [2])   # two winners in slot 0 -> guaranteed fork
+    tree = build_tree(active, winners, FixedSlotLatency(0), np.random.default_rng(0))
+    # blocks 1 and 2 are siblings at height 1 on genesis
+    assert tree.parent[1] == 0 and tree.parent[2] == 0
+    assert tree.height[1] == 1 and tree.height[2] == 1
+    # block 3 at slot 1 extends one of them (height 2)
+    assert tree.height[3] == 2
+    assert len(tree.canonical_chain()) == 2
+
+
+def test_self_extension_within_latency():
+    # A single node winning consecutive slots builds on its own block despite latency.
+    active = np.array([0, 1], dtype=np.int64)
+    winners = _winners([5], [5])
+    tree = build_tree(active, winners, FixedSlotLatency(10), np.random.default_rng(0))
+    assert tree.parent[2] == 1          # node 5 self-extends
+    assert tree.height.tolist() == [0, 1, 2]
+
+
+def test_ancestors():
+    active = np.array([0, 1, 2], dtype=np.int64)
+    winners = _winners([0], [0], [0])   # one node, clean chain
+    tree = build_tree(active, winners, FixedSlotLatency(0), np.random.default_rng(0))
+    assert tree.ancestors(3) == [3, 2, 1]
diff --git a/tools/simulators/tsi/tsi-sim/tests/test_lottery.py b/tools/simulators/tsi/tsi-sim/tests/test_lottery.py
new file mode 100644
index 0000000..a9a3a1b
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/tests/test_lottery.py
@@ -0,0 +1,41 @@
+import numpy as np
+
+from tsi_sim import lottery
+
+
+def test_phi_bounds():
+    f = 1 / 30
+    assert lottery.phi(f, 0.0) == 0.0
+    # a single all-stake node (alpha=1) wins at exactly rate f
+    assert abs(lottery.phi(f, 1.0) - f) < 1e-12
+
+
+def test_win_probs_monotone_in_stake():
+    stake = np.array([1.0, 2.0, 3.0])
+    p = lottery.win_probs(stake, d_est=6.0, f=1 / 30)
+    assert np.all(np.diff(p) > 0)
+
+
+def test_sample_wins_sorted_and_rate():
+    rng = np.random.default_rng(0)
+    n, slots = 500, 4000
+    p = np.full(n, 0.001)
+    ws, wn = lottery.sample_wins(p, slots, rng, chunk=512)
+    assert np.all(np.diff(ws) >= 0)                     # sorted by slot
+    assert ws.shape == wn.shape
+    assert np.all((wn >= 0) & (wn < n))
+    # expected wins ~ n * slots * p
+    assert abs(ws.size - n * slots * 0.001) < 4 * np.sqrt(n * slots * 0.001)
+
+
+def test_group_by_slot():
+    ws = np.array([0, 0, 2, 5, 5, 5])
+    wn = np.array([3, 7, 1, 2, 4, 9])
+    active, groups = lottery.group_by_slot(ws, wn)
+    assert list(active) == [0, 2, 5]
+    assert [g.tolist() for g in groups] == [[3, 7], [1], [2, 4, 9]]
+
+
+def test_group_by_slot_empty():
+    active, groups = lottery.group_by_slot(np.empty(0, int), np.empty(0, int))
+    assert active.size == 0 and groups == []
diff --git a/tools/simulators/tsi/tsi-sim/tests/test_theory_convergence.py b/tools/simulators/tsi/tsi-sim/tests/test_theory_convergence.py
new file mode 100644
index 0000000..5490d96
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/tests/test_theory_convergence.py
@@ -0,0 +1,49 @@
+"""End-to-end statistical checks against closed-form theory (scaled k)."""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from tsi_sim.config import SimConfig
+from tsi_sim.engine import run_trajectory
+from tsi_sim.theory import expected_ratio
+
+F = 1 / 30
+
+
+def _tail_mean(col, **cfg):
+    reps = cfg.pop("reps", 6)
+    burn = cfg["epochs"] // 2
+    vals = []
+    for r in range(reps):
+        df = pd.DataFrame(run_trajectory(SimConfig(replicate=r, **cfg)))
+        vals.append(df[col].iloc[burn:].mean())
+    return float(np.mean(vals))
+
+
+def test_baseline_exact_without_forks():
+    # No latency, no uncles: active-slot rate == f, so the estimate is unbiased.
+    ratio = _tail_mean("ratio", n_nodes=1000, stake_dist="uniform", latency=0,
+                       max_uncles=0, k=64, epochs=35, genesis_d_factor=0.5, reps=6)
+    assert abs(ratio - 1.0) < 0.02
+
+
+@pytest.mark.slow
+def test_u0_matches_expected_ratio():
+    cfg = dict(n_nodes=1000, stake_dist="uniform", latency=6, max_uncles=0,
+               k=96, epochs=45, genesis_d_factor=0.5)
+    ratio = _tail_mean("ratio", reps=8, **cfg)
+    q = _tail_mean("q", reps=8, **cfg)
+    assert abs(ratio - float(expected_ratio(F, q))) < 0.02
+
+
+@pytest.mark.slow
+def test_uncles_recover_accuracy():
+    common = dict(n_nodes=1000, stake_dist="uniform", latency=8, uncle_strategy="oldest",
+                  k=96, epochs=45, genesis_d_factor=0.5)
+    r0 = _tail_mean("ratio", max_uncles=0, reps=8, **common)
+    r4 = _tail_mean("ratio", max_uncles=4, reps=8, **common)
+    qe4 = _tail_mean("q_eff", max_uncles=4, reps=8, **common)
+    assert qe4 > 0.99
+    assert abs(r4 - 1) < abs(r0 - 1)     # uncles reduce the error
+    assert abs(r4 - 1) < 0.03            # residual is the small block-count overshoot
diff --git a/tools/simulators/tsi/tsi-sim/tests/test_tsi_counting.py b/tools/simulators/tsi/tsi-sim/tests/test_tsi_counting.py
new file mode 100644
index 0000000..97fee9e
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/tests/test_tsi_counting.py
@@ -0,0 +1,74 @@
+import numpy as np
+
+from tsi_sim.blocktree import BlockTree
+from tsi_sim.tsi import density_m, referenced_uncle_ids, slot_stats, update_D
+
+
+def make_tree(slots, parents, heights, uncles):
+    n = len(slots)
+    return BlockTree(
+        slot=np.array(slots, np.int64),
+        parent=np.array(parents, np.int64),
+        height=np.array(heights, np.int64),
+        leader=np.zeros(n, np.int64),
+        uncles=uncles,
+    )
+
+
+def test_density_counts_honest_plus_deduped_uncles_in_window():
+    # canonical 1(slot0),4(slot3); orphans 2(slot1),3(slot2). block4 refs uncles 2 and 3.
+    tree = make_tree(
+        slots=[-1, 0, 1, 2, 3],
+        parents=[-1, 0, 0, 0, 1],
+        heights=[0, 1, 1, 1, 2],
+        uncles=[(), (), (), (), (2, 3)],
+    )
+    canonical = [4, 1]  # tip-first
+    # window T=10 includes all slots
+    assert density_m(tree, canonical, T=10) == 4      # 2 honest (slots 0,3) + 2 uncles (1,2)
+    # window T=2 excludes slots 2,3 -> honest {slot0}=1, uncle slot1=1 (slot2 excluded)
+    assert density_m(tree, canonical, T=2) == 2
+
+
+def test_uncle_counted_by_own_slot_and_deduped():
+    tree = make_tree(
+        slots=[-1, 0, 5, 1],
+        parents=[-1, 0, 1, 0],
+        heights=[0, 1, 2, 1],
+        uncles=[(), (), (3,), ()],   # block2 (slot5) references orphan 3 (slot1)
+    )
+    canonical = [2, 1]
+    assert referenced_uncle_ids(tree, canonical) == {3}
+    # uncle counted by its OWN slot (1), so window T=2 includes it
+    assert density_m(tree, canonical, T=2) == 2       # honest slot0 + uncle slot1
+    # window that excludes the uncle's own slot
+    assert density_m(tree, canonical, T=1) == 1       # only honest slot0
+
+
+def test_slot_stats_q_and_qeff():
+    # active slots 0,1,2 in window; honest occupies 0,2; orphan at slot1 recovered by uncle
+    tree = make_tree(
+        slots=[-1, 0, 1, 2],
+        parents=[-1, 0, 0, 1],
+        heights=[0, 1, 1, 2],
+        uncles=[(), (), (), (2,)],   # block3 refs orphan 2 (slot1)
+    )
+    canonical = [3, 1]
+    active = np.array([0, 1, 2], np.int64)
+    ref = referenced_uncle_ids(tree, canonical)
+    ss = slot_stats(tree, canonical, ref, active, T=10)
+    assert ss.n_active == 3
+    assert ss.n_honest == 2            # slots 0 and 2
+    assert ss.n_recovered == 1         # slot 1 recovered via uncle
+    assert abs(ss.q - 2 / 3) < 1e-9
+    assert abs(ss.q_eff - 1.0) < 1e-9
+
+
+def test_update_D_fixed_point():
+    f, T = 1 / 30, 3000
+    m = int(round(T * f))              # measured density == f -> D unchanged
+    assert abs(update_D(1000.0, m, T, f, beta=1.0) - 1000.0) < 1e-6
+    # measured below f -> estimate drops
+    assert update_D(1000.0, m - 20, T, f, 1.0) < 1000.0
+    # clamp at 1
+    assert update_D(1.0, 0, T, f, 1.0) == 1.0
diff --git a/tools/simulators/tsi/tsi-sim/tests/test_uncles.py b/tools/simulators/tsi/tsi-sim/tests/test_uncles.py
new file mode 100644
index 0000000..e32b382
--- /dev/null
+++ b/tools/simulators/tsi/tsi-sim/tests/test_uncles.py
@@ -0,0 +1,68 @@
+import numpy as np
+
+from tsi_sim.blocktree import BlockTree
+from tsi_sim.config import SimConfig
+from tsi_sim.uncles import annotate_uncles
+
+
+def make_tree(slots, parents, heights, leaders):
+    n = len(slots)
+    return BlockTree(
+        slot=np.array(slots, np.int64),
+        parent=np.array(parents, np.int64),
+        height=np.array(heights, np.int64),
+        leader=np.array(leaders, np.int64),
+        uncles=[() for _ in range(n)],
+    )
+
+
+def _canonical_and_orphan_tree():
+    # genesis(0); canonical chain 1(slot0)->3(slot3)->4(slot5); orphan 2(slot1)
+    tree = make_tree(
+        slots=[-1, 0, 1, 3, 5],
+        parents=[-1, 0, 0, 1, 3],
+        heights=[0, 1, 1, 2, 3],
+        leaders=[-1, 0, 1, 2, 3],
+    )
+    canonical = [4, 3, 1]  # tip-first
+    return tree, canonical
+
+
+def test_oldest_selection_and_window():
+    tree, canonical = _canonical_and_orphan_tree()
+    cfg = SimConfig(max_uncles=1, uncle_window=300, uncle_strategy="oldest")
+    annotate_uncles(tree, canonical, cfg, np.random.default_rng(0))
+    # orphan 2 (slot1) is within window of block 3 (slot3) -> referenced there
+    referenced = {u for b in canonical for u in tree.uncles[b]}
+    assert referenced == {2}
+
+
+def test_no_uncles_when_u_zero():
+    tree, canonical = _canonical_and_orphan_tree()
+    cfg = SimConfig(max_uncles=0)
+    annotate_uncles(tree, canonical, cfg, np.random.default_rng(0))
+    assert all(tree.uncles[b] == () for b in canonical)
+
+
+def test_window_excludes_out_of_range_orphan():
+    tree, canonical = _canonical_and_orphan_tree()
+    cfg = SimConfig(max_uncles=1, uncle_window=1, uncle_strategy="oldest")
+    annotate_uncles(tree, canonical, cfg, np.random.default_rng(0))
+    # orphan 2 at slot1; nearest canonical after it is block3 at slot3 -> gap 2 > W=1
+    referenced = {u for b in canonical for u in tree.uncles[b]}
+    assert referenced == set()
+
+
+def test_dedup_across_ancestors():
+    # Two canonical blocks both within window of the single orphan: only one references it.
+    tree = make_tree(
+        slots=[-1, 0, 1, 2, 3],
+        parents=[-1, 0, 0, 1, 3],
+        heights=[0, 1, 1, 2, 3],
+        leaders=[-1, 0, 9, 2, 3],
+    )
+    canonical = [4, 3, 1]
+    cfg = SimConfig(max_uncles=4, uncle_window=300, uncle_strategy="oldest")
+    annotate_uncles(tree, canonical, cfg, np.random.default_rng(0))
+    counts = sum(len(tree.uncles[b]) for b in canonical)
+    assert counts == 1  # orphan 2 referenced exactly once despite two eligible blocks