espg · espg · Jul 22, 2026 · Jul 22, 2026 · Jul 22, 2026 · Jul 22, 2026
diff --git a/README.md b/README.md
@@ -27,13 +27,16 @@ series — are documented in
 
 ## Performance
 
-Mortie uses **Rust-accelerated** morton indexing functions for high performance. The Rust implementation provides dramatic speedups:
-
-| Dataset Size | Rust | Python (reference) | Speedup |
-|--------------|------|--------------------|---------|
-| 1,000 values | 1.93 ms | 4.14 ms | **2.1x** |
-| 100,000 values | 1.85 ms | 410.59 ms | **222x** |
-| 1.2M coordinates | 102.51 ms | 5.1 sec | **50x** |
+Mortie's morton core is a Rust extension and the sole runtime path — there is no
+Python implementation to fall back on — so performance is reported as **absolute
+throughput** rather than a speedup ratio. Encoding (`geo2mort`) and decoding
+(`mort2geo`) run at **tens of millions of morton indices per second** on one
+core, staying within roughly 2× across orders 4–29.
+
+See **[docs/benchmarks.md](docs/benchmarks.md)** for the full cross-order table
+(raw encode / decode throughput and coverage timing at orders 4 / 12 / 18 / 29),
+regenerated in place by a committed script. Cell counts there are deterministic;
+timings are machine/run dependent.
 
 Pre-built wheels are available for Linux, macOS, and Windows. The Rust extension is required and is included in all pip-installed wheels.
 

diff --git a/bench_cross_order.py b/bench_cross_order.py
@@ -0,0 +1,107 @@
+"""Cross-order throughput benchmark — writes the table straight into
+docs/benchmarks.md between the BENCH_CROSS_ORDER markers (and prints it).
+
+    python bench_cross_order.py
+
+Self-contained: fixed-seed coordinate arrays, no external data. Reports raw
+throughput (millions of morton indices per second) for encode (`geo2mort`) and
+decode (`mort2geo`), plus mixed-order coverage (`morton_coverage_moc`) timing,
+at representative orders. Throughput/timings are machine/run dependent; cell
+counts are deterministic (fixed seed / fixed polygon).
+
+Note on the coverage column: coverage cost scales with the polygon's *boundary
+length measured in cells* (~2**order per boundary edge for the MOC), so a fixed
+polygon covered at order 29 explodes. The coverage input is therefore a small
+fixed ~0.01 degree (~1 km) box, chosen so order 29 stays tractable; see
+docs/coverage_methods.md for the flat-vs-MOC / precision / budget trade-offs on
+real-world polygons. Flat `morton_coverage` is not benchmarked cross-order
+because it scales as ~4**order along the boundary and exhausts memory well
+before order 29.
+"""
+import time
+from pathlib import Path
+
+import numpy as np
+
+import mortie
+
+DOC = Path("docs/benchmarks.md")
+START, END = "<!-- BENCH_CROSS_ORDER:START -->", "<!-- BENCH_CROSS_ORDER:END -->"
+
+SEED = 20260722
+N = 1_000_000          # coordinates encoded/decoded per order
+ORDERS = [4, 12, 18, 29]
+
+# Small fixed box (~0.01 degree ~ 1 km) for the coverage column; deliberately
+# small so order-29 MOC coverage stays tractable (see module docstring).
+BOX_LAT = np.array([10.0, 10.0, 10.01, 10.01])
+BOX_LON = np.array([20.0, 20.01, 20.01, 20.0])
+
+
+def coords(n=N, seed=SEED):
+    """Fixed-seed lat/lon arrays covering the whole sphere."""
+    rng = np.random.default_rng(seed)
+    lat = rng.uniform(-89.0, 89.0, n)
+    lon = rng.uniform(-180.0, 180.0, n)
+    return lat, lon
+
+
+def timed(fn, rep):
+    ts, out = [], None
+    for _ in range(rep):
+        t0 = time.perf_counter()
+        out = fn()
+        ts.append(time.perf_counter() - t0)
+    return float(np.median(ts)), out
+
+
+def build_table():
+    lat, lon = coords()
+    # Warm the extension once before timing anything: the first parallel call
+    # spins up the rayon threadpool and touches cold caches -- a one-time cost
+    # that is ~5x the warm time on a small cover (and negligible on a large
+    # batch). We report steady-state, so a throwaway call absorbs it here.
+    mortie.morton_coverage_moc(BOX_LAT, BOX_LON, order=6)
+    rows = ["| order | encode (M idx/s) | decode (M idx/s) | coverage (cells / ms) |",
+            "|--:|--:|--:|--:|"]
+    for order in ORDERS:
+        # encode: geo2mort(lat, lon, order) -> uint64 morton array
+        morton = mortie.geo2mort(lat, lon, order=order)  # warmup + decode input
+        t_enc, _ = timed(lambda: mortie.geo2mort(lat, lon, order=order), 5)
+        enc_mps = N / t_enc / 1e6
+
+        # decode: mort2geo(morton) -> (lat, lon)
+        mortie.mort2geo(morton)  # warmup
+        t_dec, _ = timed(lambda: mortie.mort2geo(morton), 5)
+        dec_mps = N / t_dec / 1e6
+
+        # coverage: warmed above; median of a few reps at cheap orders, a
+        # single timed call at order 29 (seconds-scale, one sample is plenty).
+        reps = 1 if order >= 24 else 5
+        t_cov, cov = timed(
+            lambda: mortie.morton_coverage_moc(BOX_LAT, BOX_LON, order=order),
+            reps)
+
+        rows.append(f"| {order} | {enc_mps:,.1f} | {dec_mps:,.1f} | "
+                    f"{len(cov):,}c / {t_cov * 1e3:.0f}ms |")
+    return "\n".join(rows)
+
+
+def main():
+    table = build_table()
+    print(table)
+    doc = DOC.read_text()
+    if START not in doc or END not in doc:
+        raise SystemExit(f"markers {START} / {END} not found in {DOC}")
+    pre = doc[: doc.index(START) + len(START)]
+    post = doc[doc.index(END):]
+    note = (f"Encode/decode throughput measured over {N:,} fixed-seed coordinates; "
+            "coverage over a fixed ~0.01 degree box. `M idx/s` = millions of morton "
+            "indices per second, `c` = cell count, `ms` = milliseconds. Throughput "
+            "and timings are machine/run dependent; cell counts are deterministic.")
+    DOC.write_text(f"{pre}\n\n{table}\n\n{note}\n\n{post}")
+    print(f"\nwrote table into {DOC}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/bench_matrix.py b/bench_matrix.py
@@ -9,14 +9,20 @@
 flood-fill baseline is not reproduced here (that code was removed) — it lives as
 static prose in the doc.
 """
+import json
+import subprocess
+import sys
 import time
 from pathlib import Path
+
 import numpy as np
+
 import mortie
 
 DATA = Path("mortie/tests/Ant_Grounded_DrainageSystem_Polygons.txt")
 DOC = Path("docs/coverage_methods.md")
 START, END = "<!-- BENCH_MATRIX:START -->", "<!-- BENCH_MATRIX:END -->"
+WSTART, WEND = "<!-- BENCH_WARMUP:START -->", "<!-- BENCH_WARMUP:END -->"
 
 
 def basin(n_verts=None):
@@ -70,18 +76,62 @@ def build_table():
     return "\n".join(rows)
 
 
+# A genuine "cold" (first-call) measurement can only be the *first* parallel call
+# in a process — the rayon threadpool and caches warm globally after that. So each
+# warm-up row is measured in its own fresh subprocess. Contrasts a tiny cover
+# (fixed startup dominates) with the full basin (startup amortized away).
+WARMUP_CASES = [("~1 km box", "box", 11), ("81.6k-vert basin", "basin", 10)]
+
+_PROBE = """
+import time, json, sys, numpy as np, mortie
+kind, order = sys.argv[1], int(sys.argv[2])
+if kind == "box":
+    la = np.array([10.0, 10.0, 10.01, 10.01]); lo = np.array([20.0, 20.01, 20.01, 20.0])
+else:
+    d = np.loadtxt("mortie/tests/Ant_Grounded_DrainageSystem_Polygons.txt")
+    b = d[d[:, 2] == 1]; la, lo = b[:, 0].copy(), b[:, 1].copy()
+t = time.perf_counter(); mortie.morton_coverage_moc(la, lo, order=order)
+cold = (time.perf_counter() - t) * 1e3
+ts = []
+for _ in range(7):
+    t = time.perf_counter(); mortie.morton_coverage_moc(la, lo, order=order)
+    ts.append((time.perf_counter() - t) * 1e3)
+print(json.dumps({"cold": cold, "warm": float(np.median(ts))}))
+"""
+
+
+def warmup_table():
+    rows = ["| MOC cover | cold (first call) | warm (steady state) | ratio |",
+            "|--|--:|--:|--:|"]
+    for label, kind, order in WARMUP_CASES:
+        out = subprocess.run([sys.executable, "-c", _PROBE, kind, str(order)],
+                             capture_output=True, text=True, check=True)
+        d = json.loads(out.stdout)
+        rows.append(f"| {label}, order {order} | {d['cold']:.1f} ms | "
+                    f"{d['warm']:.1f} ms | {d['cold'] / d['warm']:.1f}x |")
+    return "\n".join(rows)
+
+
+def replace_block(doc, start, end, body):
+    if start not in doc or end not in doc:
+        raise SystemExit(f"markers {start} / {end} not found in {DOC}")
+    return f"{doc[: doc.index(start) + len(start)]}\n\n{body}\n\n{doc[doc.index(end):]}"
+
+
 def main():
     table = build_table()
     print(table)
+    warmup = warmup_table()
+    print("\n" + warmup)
+    note = ("`c` = cell count, `ms` = milliseconds. Matrix timings are the warm "
+            "median (each method is called once to warm up, then timed); see the "
+            "first-call warm-up table for the one-time cold cost. Timings are "
+            "machine/run dependent; cell counts are deterministic.")
     doc = DOC.read_text()
-    if START not in doc or END not in doc:
-        raise SystemExit(f"markers {START} / {END} not found in {DOC}")
-    pre = doc[: doc.index(START) + len(START)]
-    post = doc[doc.index(END):]
-    note = ("`c` = cell count, `ms` = milliseconds (machine/run dependent; cell "
-            "counts are deterministic).")
-    DOC.write_text(f"{pre}\n\n{table}\n\n{note}\n\n{post}")
-    print(f"\nwrote table into {DOC}")
+    doc = replace_block(doc, START, END, f"{table}\n\n{note}")
+    doc = replace_block(doc, WSTART, WEND, warmup)
+    DOC.write_text(doc)
+    print(f"\nwrote tables into {DOC}")
 
 
 if __name__ == "__main__":

diff --git a/docs/benchmarks.md b/docs/benchmarks.md
@@ -0,0 +1,101 @@
+# Benchmarks
+
+Absolute throughput for the core morton operations across representative HEALPix
+orders. `mortie`'s morton core is Rust-only — there is no Python fallback to
+compare against (the Rust extension is the sole runtime path), so these are raw
+numbers, not a speedup ratio.
+
+The two headline operations are the coordinate <-> morton conversions:
+
+- **encode** — `geo2mort(lats, lons, order)`: lat/lon -> `uint64` morton indices.
+- **decode** — `mort2geo(morton)`: morton indices -> lat/lon of the cell centre.
+
+Both are reported as **millions of morton indices per second** (`M idx/s`),
+measured over a fixed-seed array of one million coordinates spanning the sphere.
+Higher is better; the order-encoded morton means throughput is nearly flat
+across orders (the work is the same per element).
+
+The **coverage** column times `morton_coverage_moc` — the compact mixed-order
+polygon cover — on a small fixed box. It scales with the polygon's *boundary
+length in cells* rather than element count, so unlike encode/decode it grows
+steeply with order; see [coverage_methods.md](coverage_methods.md) for the
+flat-vs-MOC / precision / budget trade-offs and a per-method matrix.
+
+## Cross-order throughput
+
+The table below is regenerated by `bench_cross_order.py` (run from the repo
+root) — it writes itself in place between the markers:
+
+<!-- BENCH_CROSS_ORDER:START -->
+
+| order | encode (M idx/s) | decode (M idx/s) | coverage (cells / ms) |
+|--:|--:|--:|--:|
+| 4 | 47.9 | 13.1 | 1c / 0ms |
+| 12 | 38.9 | 10.5 | 4c / 0ms |
+| 18 | 29.9 | 9.3 | 186c / 1ms |
+| 29 | 26.9 | 7.6 | 419,292c / 4306ms |
+
+Encode/decode throughput measured over 1,000,000 fixed-seed coordinates; coverage over a fixed ~0.01 degree box. `M idx/s` = millions of morton indices per second, `c` = cell count, `ms` = milliseconds. Throughput and timings are machine/run dependent; cell counts are deterministic.
+
+<!-- BENCH_CROSS_ORDER:END -->
+
+### Notes on the numbers
+
+- **Cell/coverage counts are deterministic** (a fixed seed and a fixed polygon);
+  **throughput and timings are machine/run dependent** — treat them as
+  order-of-magnitude, not reproducible constants.
+- **encode/decode are element-wise**, so `M idx/s` is roughly order-independent;
+  the coordinate array is one million points and each order re-encodes the same
+  points. Each throughput figure is the warm median of five timed runs.
+- **The coverage `ms` is a warm median** (five reps at cheap orders; a single
+  timed call at order 29, which is seconds-scale) taken after a throwaway warm-up
+  call — see *First-call warm-up* below. The cell count beside it is exact and
+  reproducible.
+- **Coverage is deliberately scoped small.** The coverage input is a fixed
+  ~0.01 degree (~1 km) box, chosen so order 29 stays tractable. Coverage cost
+  grows with boundary length measured in cells (~2\*\*order per boundary edge for
+  the MOC), so a large real-world polygon covered at order 29 would be far more
+  expensive — cover it at a working order or with a `max_cells=` budget instead
+  (see [coverage_methods.md](coverage_methods.md)). Flat `morton_coverage` is not
+  benchmarked cross-order because it scales as ~4\*\*order along the boundary and
+  exhausts memory well before order 29.
+
+### First-call warm-up
+
+The **first** call into the extension pays a one-time cost — the `rayon`
+threadpool spins up and caches are cold — that later calls do not. The size of
+the penalty is inversely proportional to the per-call work, so it matters most
+for small, latency-sensitive calls:
+
+| operation (measured) | cold (first call) | warm (steady state) | ratio |
+|---|--:|--:|--:|
+| `morton_coverage_moc`, ~1 km box, order 11 | ~1.1 ms | ~0.2 ms | ~5x |
+| `mort2geo`, 1M indices, order 12 | ~196 ms | ~132 ms | ~1.5x |
+| `geo2mort`, 1M indices, order 12 | ~47 ms | ~44 ms | ~1.05x |
+
+So a big batch encode/decode barely notices it, but a small cover (or a single
+point, or a tight per-call loop) can run several times slower on its first
+invocation. Two practical consequences:
+
+- **When benchmarking**, make one throwaway call before timing (this page's
+  generator warms the extension once up front, then reports the steady state).
+- **In production**, if first-call latency matters (a request path, an
+  interactive tool), warm the functions you use once at startup with a small
+  throwaway call; steady-state throughput is what the table above reports.
+
+## Regenerating benchmarks
+
+Both benchmark tables regenerate themselves in place, between HTML-comment
+markers, on the current build. Rebuild the Rust extension first
+(`maturin develop --release`, see [../BUILDING.md](../BUILDING.md)), then run
+from the repo root:
+
+```bash
+python bench_cross_order.py   # this page: cross-order encode/decode/coverage
+python bench_matrix.py        # docs/coverage_methods.md: coverage-method matrix
+```
+
+Each script prints its table and rewrites its target doc between the markers, so
+committing the doc after a run keeps the published numbers current. Cell counts
+are deterministic; timings are machine dependent — regenerate on the machine you
+want to quote.
diff --git a/docs/coverage_methods.md b/docs/coverage_methods.md
@@ -10,6 +10,16 @@ one-time `O(V)` setup and per-boundary-cell work grows with local edge density
 but far more gently than the old `O(cells × vertices)` flood-fill (a 1M-vertex
 polygon covers ~40× faster); see the benchmark matrix below.
 
+> **First-call warm-up.** The first coverage call in a process spins up the
+> `rayon` threadpool and runs on cold caches — a one-time cost that is a large
+> fraction of the runtime for a *small* cover (several times the warm time),
+> though negligible for a large one. If first-call latency matters (a request path or
+> interactive tool), warm it once at startup with a throwaway call —
+> e.g. `morton_coverage_moc(box_lats, box_lons, order=6)` — before the calls you
+> care about. The *First-call warm-up cost* table under the benchmark matrix
+> below measures it on a real MOC cover; steady-state timings are what the matrix
+> and [benchmarks.md](benchmarks.md) report.
+
 Two output shapes and two adaptive stop criteria are available.
 
 ## Output shapes
@@ -128,15 +138,36 @@ it writes itself in place between the markers:
 
 | verts | order | flat | MOC | MOC tol 0.5° | MOC tol 0.05° | MOC budget 2k | MOC budget 500 |
 |--:|--:|--|--|--|--|--|--|
-| 81,595 | 8 | 883c / 49ms | 196c / 55ms | 79c / 58ms | 196c / 60ms | 196c / 49ms | 196c / 49ms |
-| 81,595 | 10 | 12,461c / 45ms | 1,058c / 53ms | 79c / 38ms | 1,058c / 62ms | 867c / 58ms | 200c / 59ms |
-| 81,595 | 12 | 191,710c / 73ms | 5,146c / 77ms | 79c / 49ms | 2,039c / 61ms | 867c / 72ms | 200c / 61ms |
-| 1,000,000 | 10 | 12,461c / 1128ms | 1,058c / 1258ms | 79c / 946ms | 1,058c / 1127ms | 867c / 1505ms | 200c / 1362ms |
+| 81,595 | 8 | 883c / 148ms | 196c / 148ms | 79c / 132ms | 196c / 136ms | 196c / 182ms | 196c / 178ms |
+| 81,595 | 10 | 12,461c / 168ms | 1,058c / 162ms | 79c / 139ms | 1,058c / 156ms | 867c / 200ms | 200c / 179ms |
+| 81,595 | 12 | 191,710c / 195ms | 5,146c / 199ms | 79c / 137ms | 2,039c / 172ms | 867c / 188ms | 200c / 175ms |
+| 1,000,000 | 10 | 12,461c / 2893ms | 1,058c / 2820ms | 79c / 2302ms | 1,058c / 2759ms | 867c / 3208ms | 200c / 2997ms |
 
-`c` = cell count, `ms` = milliseconds (machine/run dependent; cell counts are deterministic).
+`c` = cell count, `ms` = milliseconds. Matrix timings are the warm median (each method is called once to warm up, then timed); see the first-call warm-up table for the one-time cold cost. Timings are machine/run dependent; cell counts are deterministic.
 
 <!-- BENCH_MATRIX:END -->
 
+### First-call warm-up cost
+
+The matrix timings above are the **warm** (steady-state) median. The **first**
+`morton_coverage_moc` call in a process additionally pays a one-time cost (the
+`rayon` threadpool spins up, caches are cold). Because that cost is fixed, its
+weight is inversely proportional to the cover's size — it dominates a tiny cover
+and vanishes on a large one. Measured as a genuine first call (each row in its
+own fresh process):
+
+<!-- BENCH_WARMUP:START -->
+
+| MOC cover | cold (first call) | warm (steady state) | ratio |
+|--|--:|--:|--:|
+| ~1 km box, order 11 | 1.1 ms | 0.2 ms | 6.2x |
+| 81.6k-vert basin, order 10 | 178.6 ms | 159.7 ms | 1.1x |
+
+<!-- BENCH_WARMUP:END -->
+
+So a realistic basin cover barely notices it, but a small cover runs several
+times slower on its first call — warm once at startup (above) if that matters.
+
 ### Reading the matrix
 
 - **MOC vs. flat:** identical coverage, far fewer cells — at order 12 the flat