diff --git a/papers/medial_tire_decompositions_of_plane_triangulations/experiments/draw_evened_leaf.py b/papers/medial_tire_decompositions_of_plane_triangulations/experiments/draw_evened_leaf.py
new file mode 100644
index 0000000..c5a3bc2
--- /dev/null
+++ b/papers/medial_tire_decompositions_of_plane_triangulations/experiments/draw_evened_leaf.py
@@ -0,0 +1,144 @@
+"""Picture: evening a terminal (leaf) triangle by the two-vertex operation:
+add y at the midpoint of uv and z at the centroid of uvt, delete uv, add edges
+xy, uy, vy, zy, zu, zv, zt.  The leaf becomes a 4-wheel tread with hub z.
+
+Panels:
+  A  before: terminal face uvt, level cycle C_k = the 3-cycle (odd seam)
+  B  after: seam u-y-v-t (length 4, even); leaf = 4-wheel with hub z (level k+1)
+  C  medial overlay with the canonical colouring: seam apexes mono-3,
+     leaf annular 4-cycle alternating 1,2 -- proper, no ears, no defect.
+"""
+
+import os
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+PAL = {0: "#e6550d", 1: "#3182bd", 2: "#31a354"}  # colours "1","2","3"
+GRAY = "#999999"
+
+u, v, t, x = (-1.0, 0.0), (1.0, 0.0), (0.0, -1.6), (0.0, 1.0)
+y = (0.0, 0.0)                       # midpoint of uv
+z = (0.0, -1.6 / 3)                  # centroid of uvt
+
+
+def mid(a, b):
+    return ((a[0] + b[0]) / 2, (a[1] + b[1]) / 2)
+
+
+def vertex(ax, p, name, dx, dy, color="black"):
+    ax.plot(*p, "o", color=color, ms=5.5, zorder=6)
+    ax.annotate(name, p, textcoords="offset points", xytext=(dx, dy),
+                fontsize=10, fontweight="bold", zorder=6)
+
+
+def panel_before(ax):
+    ax.fill([u[0], v[0], t[0]], [u[1], v[1], t[1]], color="#fde9d9")
+    for a, b in [(u, v), (v, t), (t, u)]:
+        ax.plot([a[0], b[0]], [a[1], b[1]], color="black", lw=2.4)
+    for p in (u, v):
+        ax.plot([x[0], p[0]], [x[1], p[1]], color=GRAY, lw=0.9)
+    ax.annotate("terminal face\n(leaf of tire tree)", (0, -0.62), ha="center",
+                fontsize=8, color="#b06030")
+    vertex(ax, u, "u", -12, -2); vertex(ax, v, "v", 8, -2)
+    vertex(ax, t, "t", 0, -14); vertex(ax, x, "x", 8, 2)
+    ax.annotate("(apex in tread above)", x, textcoords="offset points",
+                xytext=(16, -2), fontsize=7, color=GRAY)
+
+
+def panel_after(ax, faint=1.0):
+    # wheel faces
+    for tri, c in [((u, y, z), "#fde9d9"), ((y, v, z), "#fdf3d9"),
+                   ((v, t, z), "#fde9d9"), ((t, u, z), "#fdf3d9")]:
+        ax.fill([p[0] for p in tri], [p[1] for p in tri], color=c, alpha=faint)
+    # seam (level cycle) bold: u-y, y-v, v-t, t-u
+    for a, b in [(u, y), (y, v), (v, t), (t, u)]:
+        ax.plot([a[0], b[0]], [a[1], b[1]], color="black", lw=2.4, alpha=faint)
+    # parent spokes xu, xy, xv
+    for p in (u, y, v):
+        ax.plot([x[0], p[0]], [x[1], p[1]], color=GRAY, lw=0.9, alpha=faint)
+    # leaf spokes zu, zy, zv, zt
+    for p in (u, y, v, t):
+        ax.plot([z[0], p[0]], [z[1], p[1]], color="black", lw=1.0, alpha=faint)
+    vertex(ax, u, "u", -12, -2); vertex(ax, v, "v", 8, -2)
+    vertex(ax, t, "t", 0, -14); vertex(ax, x, "x", 8, 2)
+    vertex(ax, y, "y", 6, 6); vertex(ax, z, "z", 7, -4)
+
+
+def panel_medial(ax):
+    panel_after(ax, faint=0.35)
+    apexes = {"uy": mid(u, y), "yv": mid(y, v), "vt": mid(v, t), "tu": mid(t, u)}
+    leaf_ann = {"zu": mid(z, u), "zy": mid(z, y), "zv": mid(z, v), "zt": mid(z, t)}
+    par_ann = {"ux": mid(u, x), "xy": mid(x, y), "xv": mid(x, v)}
+    # leaf annular 4-cycle (faces uyz, yvz, vtz, tuz)
+    ring = ["zu", "zy", "zv", "zt"]
+    for i in range(4):
+        a, b = leaf_ann[ring[i]], leaf_ann[ring[(i + 1) % 4]]
+        ax.plot([a[0], b[0]], [a[1], b[1]], color="#555555", lw=1.6, zorder=4)
+    # parent annular path m_ux - m_xy - m_xv (faces xuy, xyv)
+    for a, b in [("ux", "xy"), ("xy", "xv")]:
+        ax.plot([par_ann[a][0], par_ann[b][0]], [par_ann[a][1], par_ann[b][1]],
+                color="#555555", lw=1.6, zorder=4)
+    # apex spokes: each seam apex to its two leaf-annular and two parent nbrs
+    spokes = [("uy", "zu"), ("uy", "zy"), ("yv", "zy"), ("yv", "zv"),
+              ("vt", "zv"), ("vt", "zt"), ("tu", "zt"), ("tu", "zu")]
+    for a, b in spokes:
+        pa, pb = apexes[a], leaf_ann[b]
+        ax.plot([pa[0], pb[0]], [pa[1], pb[1]], color="#aaaaaa", lw=1.0, zorder=3)
+    for a, b in [("uy", "ux"), ("uy", "xy"), ("yv", "xy"), ("yv", "xv")]:
+        pa, pb = apexes[a], par_ann[b]
+        ax.plot([pa[0], pb[0]], [pa[1], pb[1]], color="#aaaaaa", lw=1.0, zorder=3)
+    # off-picture parent stubs for m_vt, m_tu
+    for k, d in [("vt", (0.25, -0.18)), ("tu", (-0.25, -0.18))]:
+        p = apexes[k]
+        ax.plot([p[0], p[0] + d[0]], [p[1], p[1] + d[1]], color="#cccccc",
+                lw=0.9, linestyle=":", zorder=2)
+    # colours: apexes mono-3 (green); leaf ring alternating 1,2; parent 1,2,1
+    col = {}
+    for k in apexes: col[("a", k)] = 2
+    for k, c in zip(ring, (0, 1, 0, 1)): col[("l", k)] = c
+    for k, c in zip(("ux", "xy", "xv"), (0, 1, 0)): col[("p", k)] = c
+    for k, p in apexes.items():
+        ax.plot(*p, "o", color=PAL[2], ms=11, markeredgecolor="black", zorder=5)
+    for k, p in leaf_ann.items():
+        ax.plot(*p, "o", color=PAL[col[("l", k)]], ms=9,
+                markeredgecolor="black", zorder=5)
+    for k, p in par_ann.items():
+        ax.plot(*p, "o", color=PAL[col[("p", k)]], ms=9,
+                markeredgecolor="black", zorder=5)
+    lbl = {"uy": (-26, 4), "yv": (12, 4), "vt": (12, -2), "tu": (-30, -2)}
+    for k, p in apexes.items():
+        ax.annotate(f"m_{k}", p, textcoords="offset points", xytext=lbl[k],
+                    fontsize=7, zorder=6)
+
+
+fig, axes = plt.subplots(1, 3, figsize=(14, 4.8))
+for ax in axes:
+    ax.set_xlim(-1.7, 1.7)
+    ax.set_ylim(-2.1, 1.45)
+    ax.set_aspect("equal")
+    ax.axis("off")
+
+panel_before(axes[0])
+axes[0].set_title("A. before: leaf = terminal face uvt\nseam C_k = 3-cycle (odd)",
+                  fontsize=9)
+panel_after(axes[1])
+axes[1].set_title("B. add y = mid(uv), z = centroid; delete uv;\n"
+                  "add xy, uy, vy, zy, zu, zv, zt\n"
+                  "seam u-y-v-t (even);  leaf = 4-wheel, hub z (level k+1)",
+                  fontsize=9)
+panel_medial(axes[2])
+axes[2].set_title("C. medial + canonical colouring:\nseam apexes all 3 (green), "
+                  "leaf ring alternates 1,2 — proper, no ears", fontsize=9)
+
+fig.suptitle(
+    "Evening a terminal leaf with the two-vertex operation (y splits uv under x; "
+    "z stellates the leaf as a wheel hub).\n"
+    "Both new vertices have degree 4; the leaf tread is a 4-wheel with an even "
+    "annular cycle, so the monochromatic-3 seam is VALID — no leaf defect.",
+    fontsize=10)
+fig.tight_layout(rect=(0, 0, 1, 0.86))
+out = os.path.join(HERE, "evened_leaf.png")
+fig.savefig(out, dpi=170)
+print("wrote", out)
diff --git a/papers/medial_tire_decompositions_of_plane_triangulations/experiments/even_program_findings.md b/papers/medial_tire_decompositions_of_plane_triangulations/experiments/even_program_findings.md
new file mode 100644
index 0000000..32dc052
--- /dev/null
+++ b/papers/medial_tire_decompositions_of_plane_triangulations/experiments/even_program_findings.md
@@ -0,0 +1,76 @@
+# The even-level-cycle colouring program
+
+A constructive route distinct from the `R_T` composition line. Idea: surger a
+triangulation `G` so that **every level cycle is even**, take the resulting
+*canonical even colouring* of `M(G')` (no 4CT used), then **remove the planted
+vertices** by Kempe switches, landing on a proper 3-colouring of `M(G)` — i.e. a
+Tait/4CT colouring of `G`.
+
+Scripts: `kempe_even_program_harness.py`, `draw_evened_leaf.py` (`evened_leaf.png`).
+
+## The two surgeries
+
+- **Leaf gadget (two vertices).** On a terminal triangle `uvt` with outer apex
+  `x`: add `y = mid(uv)` and hub `z`; delete `uv`; add `xy, uy, vy, zy, zu, zv,
+  zt`. Both new vertices have degree 4; the seam becomes `u-y-v-t` (even) and the
+  leaf becomes a **4-wheel** with hub `z`. No ears, no chord — the monochromatic-3
+  seam stays valid, so **leaves create no colouring defect**. (Earlier one-vertex
+  chord version forced a `{0011,0101}` defect; this is strictly better.)
+- **Diamond.** On an odd internal seam edge `uv` with apexes `x,t`: delete `uv`,
+  add `w ~ u,v,x,t` (degree 4). Flips that seam's parity.
+
+By `n_T = p + Σq_i + 2b`, evening every internal seam makes every annular cycle
+even **except the root** (the outer triangle's odd charge `Σ_T n_T ≡ 3 (mod 2)` is
+invariant — confirmed; the root is handled as the one unavoidable defect region,
+solved by local backtracking).
+
+## Canonical even colouring (constructive, no 4CT)
+
+Every level-edge medial vertex → colour 3; every non-root annular cycle alternates
+1,2; the root region solved by DFS. Proper because each apex is forced 3 between
+two `{1,2}` pairs and (in the non-degenerate tread model) no two level edges are
+consecutive around a vertex or face.
+
+## Removal conditions (degree-4 Kempe reduction — the historically *safe* case)
+
+- **Diamond** `w` (quad `u-x-v-t`, restore diagonal `uv`): removable iff the pair
+  `(m_ux,m_xv)` is distinct, `(m_vt,m_tu)` is distinct, ≤2 colours total; then
+  `m_uv` takes the third.
+- **Gadget**: collapse `z` then `y` (or `y` then `z`), ending in a degree-3
+  unstellation needing a rainbow triangle. Two orders = free choice.
+
+## Status (synthetic ring triangulations, the clean-level-structure domain)
+
+Pipeline runs end to end. Surgery, canonical colouring, and gadget removal all
+work. The program now lands squarely on the **cycle layer**:
+
+```
+60 random ring triangulations:  39 ok,  21 fail:diamond-switch
+```
+
+**Crucial diagnostic:** for a failing case, a simultaneously-removable proper
+3-colouring of `M(G')` was shown to **exist** (it must — `M(G)` is 3-colourable).
+So `fail:diamond-switch` is **not** non-existence; it is **Kempe-reachability** —
+whether switches carry the *canonical even* colouring to a descendable one. That is
+exactly the conjecture's core, and the harness has localised the entire program
+difficulty to it, with everything upstream constructive.
+
+**Why greedy fails (and what's next).** Diamonds on different odd seams share
+*vertical* `{1,3}`-Kempe cycles, so per-diamond local switching cannot satisfy them
+simultaneously. The principled solve is joint: vertices = `{1,3}`-Kempe cycles,
+one edge per diamond joining its two side cycles; removability for all diamonds at
+once = a consistent XOR assignment = **bipartiteness** of that graph (no self-loop =
+the side cycles differ; no odd cycle = no three diamonds whose side cycles form a
+triangle). Insertion-site choice (which seam edge) and tread phase are the control
+knobs. Building this joint solver — and finding the smallest configuration, if any,
+forcing a self-loop or odd cycle — is the next step and the exact thing a proof
+would need to rule out.
+
+## Caveats / domain
+
+- Real plantri triangulations mostly `skip:chord-level-edge` under BFS-from-outer
+  level structure — a reflection of how restrictive the clean nested-tire level
+  structure is, not a harness bug. The synthetic concentric-ring generator produces
+  the clean domain the program is stated for.
+- Root defect and the (deferred) outer-face handling are localised; the user has a
+  separate idea for the outer face.
diff --git a/papers/medial_tire_decompositions_of_plane_triangulations/experiments/evened_leaf.png b/papers/medial_tire_decompositions_of_plane_triangulations/experiments/evened_leaf.png
new file mode 100644
index 0000000..54023a9
Binary files /dev/null and b/papers/medial_tire_decompositions_of_plane_triangulations/experiments/evened_leaf.png differ
diff --git a/papers/medial_tire_decompositions_of_plane_triangulations/experiments/kempe_even_program_harness.py b/papers/medial_tire_decompositions_of_plane_triangulations/experiments/kempe_even_program_harness.py
new file mode 100644
index 0000000..1622d2f
--- /dev/null
+++ b/papers/medial_tire_decompositions_of_plane_triangulations/experiments/kempe_even_program_harness.py
@@ -0,0 +1,847 @@
+"""Harness for the even-level-cycle colouring program.
+
+Pipeline, per plane triangulation G (from plantri, sphere embedding, first
+face chosen as outer):
+
+  1. level structure from the outer triangle (BFS); seams (level cycles) per
+     depth; classify faces; detect unsupported degeneracies (chords, flat
+     faces that are not terminal triangles, malformed seams) and skip those.
+  2. surgeries: the two-vertex LEAF GADGET on every terminal triangle face
+     (add y on a seam edge under apex x, hub z at the face centre; delete uv;
+     leaf becomes a 4-wheel) and a DIAMOND insertion on every remaining odd
+     internal seam.  All internal seams become even; by n_T = p + sum(q) + 2b
+     every annular cycle is then even except the root tread's (the outer
+     triangle's odd charge is invariant).
+  3. canonical colouring of M(G'): every level-edge medial vertex (seam
+     apexes and bite fins) gets colour 3; every non-root annular cycle
+     alternates {1,2}; the root region (root annulus + the outer triangle's
+     three mutually-adjacent medials) is solved by backtracking -- the one
+     unavoidable defect.
+  4. descent: remove the planted vertices one at a time.  Diamonds collapse
+     by the degree-4 condition (corner pairs distinct, <=2 colours, restored
+     edge takes the third); the gadget collapses in two steps (either order),
+     ending with a degree-3 unstellation needing a rainbow triangle.  Before
+     each step a bounded Kempe-switch search (components through the support,
+     all three colour pairs, subsets of size <= 2) establishes the condition.
+  5. verify the final colouring is a proper 3-colouring of M(G) -- which is a
+     Tait/4CT colouring of the original triangulation.
+
+Reports, per graph: ok / skip:<degeneracy> / fail:<stage>.
+
+Run:  python3 kempe_even_program_harness.py --min-n 6 --max-n 10
+"""
+
+from __future__ import annotations
+
+import argparse
+import os
+import subprocess
+import sys
+from collections import defaultdict, deque
+from itertools import combinations
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+PLANTRI = os.path.join(HERE, "..", "..", "..", "plantri", "plantri")
+
+Edge = frozenset
+
+
+def e(a, b):
+    return Edge((a, b))
+
+
+# ---------------------------------------------------------------------------
+# Plane triangulation with rotation system.
+# ---------------------------------------------------------------------------
+
+class Tri:
+    def __init__(self, rot: dict[int, list[int]]):
+        self.rot = rot  # vertex -> neighbours in cyclic (embedding) order
+
+    def copy(self):
+        return Tri({v: list(ns) for v, ns in self.rot.items()})
+
+    def vertices(self):
+        return list(self.rot)
+
+    def edges(self):
+        return {e(u, v) for u in self.rot for v in self.rot[u]}
+
+    def faces(self):
+        """All faces as vertex triples (cyclic), via next-in-rotation walk."""
+        seen = set()
+        out = []
+        for u in self.rot:
+            for v in self.rot[u]:
+                if (u, v) in seen:
+                    continue
+                face = []
+                cu, cv = u, v
+                while (cu, cv) not in seen:
+                    seen.add((cu, cv))
+                    face.append(cu)
+                    ns = self.rot[cv]
+                    cw = ns[(ns.index(cu) + 1) % len(ns)]
+                    cu, cv = cv, cw
+                out.append(tuple(face))
+        return out
+
+    def check(self):
+        nv = len(self.rot)
+        ne = len(self.edges())
+        fs = self.faces()
+        assert ne == 3 * nv - 6, f"edge count {ne} != {3*nv-6}"
+        assert all(len(f) == 3 for f in fs), "non-triangular face"
+        assert len(fs) == 2 * nv - 4, f"face count {len(fs)}"
+
+    def apexes_of(self, u, v):
+        """The two third-vertices of the faces on edge uv."""
+        out = []
+        for f in (self._face_from(u, v), self._face_from(v, u)):
+            out.append(next(w for w in f if w not in (u, v)))
+        return out  # [apex left of u->v, apex left of v->u]
+
+    def _face_from(self, u, v):
+        face = []
+        cu, cv = u, v
+        for _ in range(3):
+            face.append(cu)
+            ns = self.rot[cv]
+            cw = ns[(ns.index(cu) + 1) % len(ns)]
+            cu, cv = cv, cw
+        return face
+
+    # -- surgeries ----------------------------------------------------------
+
+    def _new_vertex(self):
+        return max(self.rot) + 1
+
+    def insert_diamond(self, u, v):
+        """Delete uv; add w adjacent to u, v and both apexes.  Returns
+        (w, u, v, x, t) with x, t the two apexes."""
+        x, t = self.apexes_of(u, v)
+        w = self._new_vertex()
+        for ordering in ([u, x, v, t], [u, t, v, x]):
+            trial = self.copy()
+            trial.rot[u][trial.rot[u].index(v)] = w
+            trial.rot[v][trial.rot[v].index(u)] = w
+            for a in (x, t):
+                i, j = trial.rot[a].index(u), trial.rot[a].index(v)
+                # insert w between u and v (they are cyclically adjacent at a)
+                k = i if (i + 1) % len(trial.rot[a]) == j else j
+                trial.rot[a].insert(k + 1, w)
+            trial.rot[w] = ordering
+            try:
+                trial.check()
+            except AssertionError:
+                continue
+            self.rot = trial.rot
+            return (w, u, v, x, t)
+        raise RuntimeError("diamond insertion failed both orientations")
+
+    def insert_leaf_gadget(self, u, v, t):
+        """Terminal face (u,v,t): add y on uv (under outer apex x) and hub z;
+        delete uv; add xy,uy,vy,zy,zu,zv,zt.  Returns (y, z, u, v, x, t)."""
+        aps = self.apexes_of(u, v)
+        x = next(a for a in aps if a != t)
+        y, z = self._new_vertex(), self._new_vertex() + 1
+        base = self.copy()
+        for uord in ([y, z], [z, y]):
+            for yord in ([x, v, z, u], [x, u, z, v]):
+                for zord in ([y, v, t, u], [y, u, t, v]):
+                    trial = base.copy()
+                    iu = trial.rot[u].index(v)
+                    trial.rot[u][iu:iu + 1] = uord
+                    iv = trial.rot[v].index(u)
+                    trial.rot[v][iv:iv + 1] = list(reversed(uord))
+                    ix, jx = trial.rot[x].index(u), trial.rot[x].index(v)
+                    k = ix if (ix + 1) % len(trial.rot[x]) == jx else jx
+                    trial.rot[x].insert(k + 1, y)
+                    it_, jt = trial.rot[t].index(u), trial.rot[t].index(v)
+                    k = it_ if (it_ + 1) % len(trial.rot[t]) == jt else jt
+                    trial.rot[t].insert(k + 1, z)
+                    trial.rot[y] = yord
+                    trial.rot[z] = zord
+                    try:
+                        trial.check()
+                        # the gadget must create the four wheel faces
+                        fs = {frozenset(f) for f in trial.faces()}
+                        need = [{u, y, z}, {y, v, z}, {v, t, z}, {t, u, z},
+                                {x, u, y}, {x, y, v}]
+                        if not all(frozenset(s) in fs for s in need):
+                            continue
+                    except AssertionError:
+                        continue
+                    self.rot = trial.rot
+                    return (y, z, u, v, x, t)
+        raise RuntimeError("leaf gadget insertion failed")
+
+    def remove_degree4(self, w, a, c):
+        """Remove degree-4 w, restoring diagonal ac (a, c opposite in rot[w])."""
+        ns = self.rot[w]
+        assert len(ns) == 4 and a in ns and c in ns
+        assert (ns.index(a) - ns.index(c)) % 4 == 2, "diagonal not opposite"
+        for p in ns:
+            i = self.rot[p].index(w)
+            if p == a:
+                self.rot[p][i] = c
+            elif p == c:
+                self.rot[p][i] = a
+            else:
+                del self.rot[p][i]
+        del self.rot[w]
+
+    def remove_degree3(self, w):
+        for p in self.rot[w]:
+            self.rot[p].remove(w)
+        del self.rot[w]
+
+
+# ---------------------------------------------------------------------------
+# plantri reader.
+# ---------------------------------------------------------------------------
+
+def plantri_triangulations(n, limit=None):
+    raw = subprocess.run([PLANTRI, str(n)], capture_output=True).stdout
+    assert raw.startswith(b">>planar_code<<")
+    data = raw[len(b">>planar_code<<"):]
+    i = 0
+    out = []
+    while i < len(data) and (limit is None or len(out) < limit):
+        nv = data[i]
+        i += 1
+        rot = {}
+        for v in range(1, nv + 1):
+            ns = []
+            while data[i] != 0:
+                ns.append(data[i])
+                i += 1
+            i += 1
+            rot[v] = ns
+        out.append(Tri(rot))
+    return out
+
+
+# ---------------------------------------------------------------------------
+# Level / seam analysis.
+# ---------------------------------------------------------------------------
+
+class Analysis:
+    def __init__(self, g: Tri, outer: tuple[int, int, int]):
+        self.g = g
+        self.outer = outer
+        self.level = self._levels()
+        self.faces = g.faces()
+        self.face_min = {frozenset(f): min(self.level[v] for v in f)
+                         for f in self.faces}
+        self.degenerate = None
+        self._classify()
+
+    def _levels(self):
+        lev = {v: None for v in self.g.rot}
+        q = deque()
+        for v in self.outer:
+            lev[v] = 0
+            q.append(v)
+        while q:
+            u = q.popleft()
+            for w in self.g.rot[u]:
+                if lev[w] is None:
+                    lev[w] = lev[u] + 1
+                    q.append(w)
+        return lev
+
+    def _classify(self):
+        g, lev = self.g, self.level
+        outer_face = frozenset(self.outer)
+        self.flat_faces = []      # all-equal-level faces (excluding outer)
+        for f in self.faces:
+            ls = [lev[v] for v in f]
+            if ls[0] == ls[1] == ls[2] and frozenset(f) != outer_face:
+                self.flat_faces.append(tuple(f))
+        # seam edges per depth k>=1: level-k edges with exactly one face of
+        # min-level k-1
+        self.seam_edges = defaultdict(set)   # k -> set of edges
+        self.fins = set()
+        for edge in g.edges():
+            a, b = tuple(edge)
+            if lev[a] != lev[b]:
+                continue
+            k = lev[a]
+            if k == 0:
+                continue  # outer triangle
+            x, t = g.apexes_of(a, b)
+            mins = sorted([min(lev[x], k), min(lev[t], k)])
+            above = sum(1 for ap in (x, t) if lev[ap] == k - 1)
+            if above == 1:
+                self.seam_edges[k].add(edge)
+            elif above == 2:
+                self.fins.add(edge)
+            else:
+                self.degenerate = "chord-level-edge"
+                return
+        # seams must decompose into disjoint simple cycles
+        self.seams = []  # list of (k, [edges], [vertices in cyclic order])
+        for k, es in self.seam_edges.items():
+            deg = defaultdict(list)
+            for edge in es:
+                a, b = tuple(edge)
+                deg[a].append(b)
+                deg[b].append(a)
+            if any(len(ns) != 2 for ns in deg.values()):
+                self.degenerate = "non-simple-seam"
+                return
+            left = set(es)
+            while left:
+                a0, b0 = tuple(next(iter(left)))
+                cyc = [a0, b0]
+                left.discard(e(a0, b0))
+                while True:
+                    nxt = next(c for c in deg[cyc[-1]] if c != cyc[-2])
+                    if nxt == cyc[0]:
+                        break
+                    cyc.append(nxt)
+                    left.discard(e(cyc[-2], cyc[-1]))
+                left.discard(e(cyc[-1], cyc[0]))
+                self.seams.append((k, cyc))
+        # terminal triangles: flat faces whose 3 edges are all seam edges
+        self.terminal = []
+        for f in self.flat_faces:
+            a, b, c = f
+            k = self.level[a]
+            es = [e(a, b), e(b, c), e(c, a)]
+            if all(x in self.seam_edges.get(k, ()) for x in es):
+                self.terminal.append(f)
+            else:
+                self.degenerate = "non-terminal-flat-face"
+                return
+
+
+# ---------------------------------------------------------------------------
+# Medial graph + canonical colouring.
+# ---------------------------------------------------------------------------
+
+def medial_adj(g: Tri):
+    adj = defaultdict(set)
+    for v, ns in g.rot.items():
+        d = len(ns)
+        for i in range(d):
+            a, b = e(v, ns[i]), e(v, ns[(i + 1) % d])
+            adj[a].add(b)
+            adj[b].add(a)
+    return adj
+
+
+def annular_cycles(g: Tri, level):
+    """Spoke medial vertices grouped into annular cycles (per tread)."""
+    spokes = [ed for ed in g.edges() if abs(level[min(ed)] - level[max(ed)]) == 1
+              or True]
+    spokes = [ed for ed in g.edges()
+              if level[tuple(ed)[0]] != level[tuple(ed)[1]]]
+    adj = medial_adj(g)
+    sset = set(spokes)
+    comp = []
+    seen = set()
+    for s in spokes:
+        if s in seen:
+            continue
+        stack, cur = [s], set()
+        seen.add(s)
+        while stack:
+            a = stack.pop()
+            cur.add(a)
+            for b in adj[a]:
+                if b in sset and b not in seen:
+                    seen.add(b)
+                    stack.append(b)
+        comp.append(cur)
+    return comp, adj
+
+
+def order_cycle(verts, adj):
+    """Order a vertex set known to induce a cycle."""
+    v0 = next(iter(verts))
+    cyc = [v0]
+    prev = None
+    while True:
+        nbrs = [w for w in adj[cyc[-1]] if w in verts and w != prev]
+        if not nbrs:
+            return None
+        prev = cyc[-1]
+        nxt = nbrs[0]
+        if nxt == cyc[0]:
+            return cyc if len(cyc) == len(verts) else None
+        cyc.append(nxt)
+        if len(cyc) > len(verts):
+            return None
+
+
+def canonical_coloring(g: Tri, level, outer, rng=None):
+    """Canonical colouring of M(g): level-edge medials -> 2 ("colour 3"),
+    non-root annuli alternate 0,1 (random phases), root region solved by DFS.
+    Returns (coloring dict, adj) or (None, reason)."""
+    adj = medial_adj(g)
+    col = {}
+    level_edges = [ed for ed in g.edges()
+                   if level[tuple(ed)[0]] == level[tuple(ed)[1]]]
+    outer_es = {e(outer[0], outer[1]), e(outer[1], outer[2]),
+                e(outer[2], outer[0])}
+    for ed in level_edges:
+        if ed not in outer_es:
+            col[ed] = 2
+    comps, _ = annular_cycles(g, level)
+    root_comp = None
+    free = set(outer_es)
+    for comp in comps:
+        cyc = order_cycle(comp, adj)
+        if cyc is None:
+            return None, "annulus-not-cycle"
+        depth = min(level[v] for ed in comp for v in ed)
+        if depth == 0:  # root annulus (spokes touching level 0)
+            if root_comp is not None:
+                return None, "multiple-root-annuli"
+            root_comp = cyc
+            free |= set(cyc)
+            continue
+        if len(cyc) % 2 != 0:
+            return None, f"odd-non-root-annulus(len {len(cyc)})"
+        phase = rng.randint(0, 1) if rng else 0
+        for i, m in enumerate(cyc):
+            col[m] = (i + phase) % 2
+    if root_comp is None:
+        return None, "no-root-annulus"
+    # solve the free region (root annulus + outer trio) by DFS; if that
+    # fails, grow the defect region one level at a time.
+    max_level = max(level.values())
+    for grow in range(0, max_level + 1):
+        free_now = set(free)
+        if grow:
+            for ed in g.edges():
+                a, b = tuple(ed)
+                if min(level[a], level[b]) < grow:
+                    free_now.add(ed)
+        trial = {m: c for m, c in col.items() if m not in free_now}
+        order = sorted(free_now,
+                       key=lambda m: -sum(1 for x in adj[m] if x in trial))
+        colorder = [0, 1, 2]
+        if rng:
+            rng.shuffle(colorder)
+
+        def ok(m, c):
+            return all(trial.get(x) != c for x in adj[m])
+
+        def dfs(i):
+            if i == len(order):
+                return True
+            m = order[i]
+            for c in colorder:
+                if ok(m, c):
+                    trial[m] = c
+                    if dfs(i + 1):
+                        return True
+                    del trial[m]
+            return False
+
+        if dfs(0):
+            for m, c in trial.items():
+                for x in adj[m]:
+                    assert trial.get(x) != c, "canonical colouring improper"
+            return trial, None
+        if len(free_now) > 60:  # defect region growing unmanageable
+            break
+    return None, "root-unsolvable"
+
+
+# ---------------------------------------------------------------------------
+# Kempe switching + removal conditions.
+# ---------------------------------------------------------------------------
+
+def kempe_component(col, adj, start, pair):
+    if col[start] not in pair:
+        return None
+    comp = {start}
+    stack = [start]
+    while stack:
+        a = stack.pop()
+        for b in adj[a]:
+            if b in col and col[b] in pair and b not in comp:
+                comp.add(b)
+                stack.append(b)
+    return frozenset(comp)
+
+
+def switch(col, comp, pair):
+    a, b = pair
+    for m in comp:
+        col[m] = b if col[m] == a else a
+
+
+def diamond_condition(col, w_quad):
+    """w_quad = (a, b, c, d) cyclic rot of planted vertex, diagonal (a, c).
+    Quad edges ab, bc, cd, da; new faces a-b-c and a-c-d."""
+    a, b, c, d = w_quad
+    p, q, r, s = col[e(a, b)], col[e(b, c)], col[e(c, d)], col[e(d, a)]
+    if p == q or r == s:
+        return None
+    used = {p, q, r, s}
+    if len(used) > 2:
+        return None
+    return ({0, 1, 2} - used).pop()
+
+
+def rainbow_condition(col, tri):
+    a, b, c = tri
+    cols = {col[e(a, b)], col[e(b, c)], col[e(c, a)]}
+    return len(cols) == 3
+
+
+def try_establish(col, adj, support, test, max_switch=3):
+    """Bounded search: switch <= max_switch Kempe components (any pair)
+    anchored at the support medials (and their coloured neighbours) to make
+    test(col) true.  Mutates col on success; restores on failure."""
+    if test(col):
+        return True
+    anchors = set(support)
+    for m in support:
+        anchors.update(adj[m])
+    cands = []
+    seenc = set()
+    for m in anchors:
+        if m not in col:
+            continue
+        for pair in ((0, 1), (0, 2), (1, 2)):
+            comp = kempe_component(col, adj, m, pair)
+            if comp and (comp, pair) not in seenc:
+                seenc.add((comp, pair))
+                cands.append((comp, pair))
+    for k in range(1, max_switch + 1):
+        for subset in combinations(range(len(cands)), k):
+            # only co-switch components that are pairwise vertex-disjoint
+            verts = set()
+            clash = False
+            for idx in subset:
+                if verts & cands[idx][0]:
+                    clash = True
+                    break
+                verts |= cands[idx][0]
+            if clash:
+                continue
+            for idx in subset:
+                switch(col, *cands[idx])
+            if test(col):
+                return True
+            for idx in subset:
+                switch(col, *cands[idx])
+    return False
+
+
+# ---------------------------------------------------------------------------
+# The per-graph pipeline.
+# ---------------------------------------------------------------------------
+
+def quad_of(g, w, diag_a, diag_c):
+    ns = g.rot[w]
+    i = ns.index(diag_a)
+    if ns[(i + 2) % 4] != diag_c:
+        return None
+    return (ns[i], ns[(i + 1) % 4], ns[(i + 2) % 4], ns[(i + 3) % 4])
+
+
+def collapse_degree4(g, col, w, a, c):
+    """Remove planted degree-4 w restoring diagonal ac; col gets col[ac]."""
+    quad = quad_of(g, w, a, c)
+    third = diamond_condition(col, quad)
+    assert third is not None
+    for p in g.rot[w]:
+        col.pop(e(w, p), None)
+    g.remove_degree4(w, a, c)
+    col[e(a, c)] = third
+
+
+def collapse_degree3(g, col, w):
+    for p in g.rot[w]:
+        col.pop(e(w, p), None)
+    g.remove_degree3(w)
+
+
+def verify_proper(g, col):
+    adj = medial_adj(g)
+    for m in adj:
+        assert m in col, f"uncoloured medial {m}"
+        for b in adj[m]:
+            if col[m] == col[b]:
+                return False
+    return True
+
+
+def run_graph(g0: Tri, outer=None, verbose=False, attempts=4):
+    import random
+    g0.check()
+    if outer is None:
+        outer = tuple(g0.faces()[0])
+    an0 = Analysis(g0.copy(), outer)
+    if an0.degenerate:
+        return f"skip:{an0.degenerate}"
+    last = "fail:unknown"
+    for att in range(attempts):
+        rng = random.Random(1000 + att)
+        last = _attempt(g0, outer, rng, verbose)
+        if last == "ok" or last.startswith("skip"):
+            return last
+    return last
+
+
+def _attempt(g0: Tri, outer, rng, verbose=False):
+    g = g0.copy()
+    an = Analysis(g, outer)
+    if an.degenerate:
+        return f"skip:{an.degenerate}"
+
+    # --- surgeries -------------------------------------------------------
+    diamonds = []   # (w, quad(a,b,c,d) with diagonal (u,v)) for removal
+    gadgets = []    # (y, z, u, v, x, t)
+    for f in an.terminal:
+        a, b, c = f
+        y, z, u, v, x, t = g.insert_leaf_gadget(a, b, c)
+        gadgets.append((y, z, u, v, x, t))
+    # re-analyse (gadgets change seams)
+    an = Analysis(g, outer)
+    if an.degenerate:
+        return f"skip:post-gadget-{an.degenerate}"
+    for k, cyc in an.seams:
+        if len(cyc) % 2 == 0:
+            continue
+        # choose a seam edge whose below-apex is strictly deeper
+        choice = None
+        for i in range(len(cyc)):
+            a, b = cyc[i], cyc[(i + 1) % len(cyc)]
+            x, t = g.apexes_of(a, b)
+            lx, lt = an.level[x], an.level[t]
+            if {lx, lt} == {k - 1, k + 1}:
+                choice = (a, b)
+                break
+        if choice is None:
+            return "skip:no-diamond-site"
+        w, u, v, x, t = g.insert_diamond(*choice)
+        diamonds.append((w, u, v, x, t))
+    an = Analysis(g, outer)
+    if an.degenerate:
+        return f"skip:post-diamond-{an.degenerate}"
+    if any(len(cyc) % 2 for _, cyc in an.seams):
+        return "fail:seam-evening"
+
+    # --- canonical colouring ----------------------------------------------
+    col, reason = canonical_coloring(g, an.level, outer, rng=rng)
+    if col is None:
+        return f"fail:canonical-{reason}"
+
+    # --- descent -----------------------------------------------------------
+    level = an.level
+    # diamonds: remove each (diagonal uv)
+    for w, u, v, x, t in diamonds:
+        adj = medial_adj(g)
+        quad = quad_of(g, w, u, v)
+        if quad is None:
+            return "fail:diamond-quad"
+        support = [e(quad[i], quad[(i + 1) % 4]) for i in range(4)]
+        if not try_establish(col, adj, support,
+                             lambda c: diamond_condition(c, quad) is not None):
+            return "fail:diamond-switch"
+        collapse_degree4(g, col, w, u, v)
+        if not verify_proper(g, col):
+            return "fail:improper-after-diamond"
+    # gadgets: try order A (z then y) else order B (y then z)
+    for y, z, u, v, x, t in gadgets:
+        done = False
+        for first, second, tri in ((z, y, (x, u, v)), (y, z, (u, v, t))):
+            gtrial, ctrial = g.copy(), dict(col)
+            adj = medial_adj(gtrial)
+            quad = quad_of(gtrial, first, u, v)
+            if quad is None:
+                continue
+            support = [e(quad[i], quad[(i + 1) % 4]) for i in range(4)]
+            if not try_establish(ctrial, adj, support,
+                                 lambda c: diamond_condition(c, quad) is not None):
+                continue
+            collapse_degree4(gtrial, ctrial, first, u, v)
+            if not verify_proper(gtrial, ctrial):
+                continue
+            adj2 = medial_adj(gtrial)
+            a_, b_, c_ = tri
+            support2 = [e(a_, b_), e(b_, c_), e(c_, a_)]
+            if not try_establish(ctrial, adj2, support2,
+                                 lambda c: rainbow_condition(c, tri)):
+                continue
+            collapse_degree3(gtrial, ctrial, second)
+            if not verify_proper(gtrial, ctrial):
+                continue
+            g, col = gtrial, ctrial
+            done = True
+            break
+        if not done:
+            return "fail:gadget-removal"
+
+    # --- final check --------------------------------------------------------
+    if set(g.rot) != set(g0.rot):
+        return "fail:vertex-mismatch"
+    if not verify_proper(g, col):
+        return "fail:final-improper"
+    return "ok"
+
+
+# ---------------------------------------------------------------------------
+# Driver.
+# ---------------------------------------------------------------------------
+
+# ---------------------------------------------------------------------------
+# Synthetic clean ring triangulations (the programme's natural domain).
+# ---------------------------------------------------------------------------
+
+def tri_from_faces(faces):
+    """Build a Tri from a consistently oriented face list covering the sphere
+    (every directed edge used exactly once)."""
+    succ = defaultdict(dict)
+    for a, b, c in faces:
+        succ[a][b] = c
+        succ[b][c] = a
+        succ[c][a] = b
+    rot = {}
+    for v, nx in succ.items():
+        start = next(iter(nx))
+        cyc = [start]
+        for _ in range(len(nx)):
+            w = nx.get(cyc[-1])
+            if w is None or w == start:
+                break
+            cyc.append(w)
+        if len(cyc) != len(nx) or nx.get(cyc[-1]) != start:
+            return None
+        rot[v] = cyc
+    return Tri(rot)
+
+
+def ring_triangulation(ring_sizes, leaf, rng):
+    """Concentric ring triangulation: ring_sizes[0] must be 3 (outer
+    triangle); annuli between consecutive rings get random tooth words;
+    leaf in {'hub','face'} ('face' needs innermost size 3)."""
+    assert ring_sizes[0] == 3
+    rings = []
+    nid = 0
+    for r in ring_sizes:
+        rings.append(list(range(nid, nid + r)))
+        nid += r
+    faces = []
+    for k in range(len(rings) - 1):
+        A, B = rings[k], rings[k + 1]
+        a, b = len(A), len(B)
+        for _ in range(50):  # reject words that revisit a spoke
+            word = ["A"] * a + ["B"] * b
+            rng.shuffle(word)
+            spokes = {(0, 0)}
+            i = j = 0
+            good = True
+            for mv in word[:-1]:
+                if mv == "A":
+                    i += 1
+                else:
+                    j += 1
+                if (i % a, j % b) in spokes:
+                    good = False
+                    break
+                spokes.add((i % a, j % b))
+            if good:
+                break
+        else:
+            return None, None
+        i = j = 0
+        for mv in word:
+            if mv == "A":
+                faces.append((A[i % a], A[(i + 1) % a], B[j % b]))
+                i += 1
+            else:
+                faces.append((B[(j + 1) % b], B[j % b], A[i % a]))
+                j += 1
+    inner = rings[-1]
+    if leaf == "hub":
+        h = nid
+        for j in range(len(inner)):
+            faces.append((inner[j], inner[(j + 1) % len(inner)], h))
+    else:
+        assert len(inner) == 3
+        faces.append((inner[0], inner[1], inner[2]))
+    outer = rings[0]
+    faces.append((outer[2], outer[1], outer[0]))  # the unbounded face
+    for variant in (faces, [tuple(reversed(f)) for f in faces]):
+        g = tri_from_faces(variant)
+        if g is None:
+            continue
+        try:
+            g.check()
+            return g, tuple(outer)
+        except AssertionError:
+            continue
+    return None, None
+
+
+def random_profile(rng):
+    depth = rng.randint(2, 4)
+    sizes = [3] + [rng.randint(3, 8) for _ in range(depth)]
+    leaf = "face" if (sizes[-1] == 3 and rng.random() < 0.5) else "hub"
+    if leaf == "face" and sizes[-1] != 3:
+        leaf = "hub"
+    return sizes, leaf
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("--min-n", type=int, default=6)
+    ap.add_argument("--max-n", type=int, default=10)
+    ap.add_argument("--limit", type=int, default=None)
+    ap.add_argument("--synthetic", type=int, default=0,
+                    help="number of random ring triangulations to test")
+    ap.add_argument("--seed", type=int, default=1)
+    ap.add_argument("--verbose", action="store_true")
+    args = ap.parse_args()
+
+    import random
+    grand = defaultdict(int)
+    if args.synthetic:
+        rng = random.Random(args.seed)
+        tally = defaultdict(int)
+        for idx in range(args.synthetic):
+            sizes, leaf = random_profile(rng)
+            g, outer = ring_triangulation(sizes, leaf, rng)
+            if g is None:
+                tally["error:construction"] += 1
+                continue
+            try:
+                res = run_graph(g, outer=outer, verbose=args.verbose)
+            except Exception as ex:  # noqa: BLE001
+                res = f"error:{type(ex).__name__}"
+            tally[res] += 1
+            grand[res] += 1
+            if args.verbose and not res.startswith("ok"):
+                print(f"  synth #{idx} sizes={sizes} leaf={leaf}: {res}")
+        line = "  ".join(f"{k}={v}" for k, v in sorted(tally.items()))
+        print(f"synthetic ({args.synthetic}): {line}")
+    else:
+        for n in range(args.min_n, args.max_n + 1):
+            tally = defaultdict(int)
+            gs = plantri_triangulations(n, args.limit)
+            for idx, g in enumerate(gs):
+                try:
+                    res = run_graph(g, verbose=args.verbose)
+                except Exception as ex:  # noqa: BLE001
+                    res = f"error:{type(ex).__name__}"
+                tally[res] += 1
+                grand[res] += 1
+                if args.verbose and not res.startswith(("ok", "skip")):
+                    print(f"  n={n} #{idx}: {res}")
+            line = "  ".join(f"{k}={v}" for k, v in sorted(tally.items()))
+            print(f"n={n} ({len(gs)} graphs): {line}")
+            sys.stdout.flush()
+    print("\nTOTAL: " + "  ".join(f"{k}={v}" for k, v in sorted(grand.items())))
+
+
+if __name__ == "__main__":
+    main()