Small-n bridge-derivability probe: classification + invariant search

Findings at n=9 (50 triangulations, orbits fully exhaustible):
- 36 bridge-derived, 14 NOT bridge-derived. So bridge-derived is a PROPER
  subclass of derived (49 derived at n=9). All 14 non-bridge graphs are
  intertwining trees -- as are all 50, necessarily: intertwining tree
  <=> dual Hamiltonian, and the smallest non-Hamiltonian 3-connected cubic
  planar graph has 38 vertices, i.e. dual on 2n-4=38 => n=21. Hence every
  triangulation with n<=20 is an intertwining tree, and the disjunction
  "bridge-derived OR intertwining" is trivially true below n=21. The 4
  Holton-McKay duals are the first non-intertwining triangulations.
- Static parity-subgraph invariants (Betti numbers, component counts,
  cross-edge count, existence of an all-forest partition) do NOT separate
  bridge-derived from non-bridge-derived -- both classes realize beta=0
  partitions and identical ranges. Bridge-derivability is dynamical, not a
  simple static invariant; no easy obstruction.
- Side lemma: every valid parity partition of an n-vertex triangulation has
  exactly 2n-4 cross edges (intra-edges = n-2). Holds for all n=9 graphs.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

papers/plane_depth_sequencing/experiments/plane_depth_sequencing.py

@@ -0,0 +1,492 @@
+"""Plane depth sequencing on maximal planar graphs."""
+import random
+from typing import Any, TypedDict
+from sage.all import Graph, graphs  # type: ignore[attr-defined]  # pylint: disable=no-name-in-module
+from lib.colored_graphs import canonize_and_save_graph
+
+
+class DeeplyEmbeddedGraph(TypedDict):
+    graph: Graph
+    outer_cycle: list[Any]
+    plane_depth_labelling: dict[Any, int]
+    deep_embedding: Graph
+
+
+class QuadrilateralSequence(TypedDict):
+    deep_embedding: Graph
+    triangular_faces: list[frozenset[Any]]
+    depth_labelling: dict[Any, int]
+    outer_cap_vertex: Any
+    quadrilaterals: list[frozenset[frozenset[Any]]]
+    sequence: list[frozenset[frozenset[Any]]]
+    move_codes: list[int]
+
+
+def get_plane_depth_labelling(g: Graph, outer_cycle: list[Any]) -> dict[Any, int]:
+    """Return the plane depth of each vertex relative to the given outer cycle."""
+    # equivalent to the commented out naive implementation:
+    # return {v: min(g.distance(v, u) for u in outer_cycle) for v in g.vertices()}
+    source = max(g.vertices()) + 1
+    g_prime = g.copy()
+    g_prime.add_vertex(source)
+    g_prime.add_edges([(source, v) for v in outer_cycle])
+    distances = g_prime.breadth_first_search(source, report_distance=True)
+    return {v: d - 1 for v, d in distances if v != source}
+
+
+def deep_embedding(g: Graph, outer_cycle: list[Any], plane_depth_labelling: dict[Any, int] | None = None) -> Graph:
+    """
+    Return the deep embedding of g relative to outer_cycle.
+
+    For every interior triangular face whose three vertices all have the same
+    plane depth, a new vertex is added adjacent to each vertex of that face.
+    """
+    if plane_depth_labelling is None:
+        plane_depth_labelling = get_plane_depth_labelling(g, outer_cycle)
+    outer_vertices = set(outer_cycle)
+
+    embedding = g.get_embedding()
+    if embedding is None:
+        g.is_planar(set_embedding=True)
+        embedding = g.get_embedding()
+
+    g_prime = g.copy()
+    new_vertex = max(g.vertices()) + 1
+
+    for face in g.faces(embedding):
+        face_vertices = [u for u, v in face]
+        if set(face_vertices) == outer_vertices:
+            continue
+        if plane_depth_labelling[face_vertices[0]] == plane_depth_labelling[face_vertices[1]] == plane_depth_labelling[face_vertices[2]]:
+            g_prime.add_vertex(new_vertex)
+            g_prime.add_edges([(new_vertex, v) for v in face_vertices])
+            new_vertex += 1
+
+    return g_prime
+
+
+def _triangle_type(face: frozenset[Any], depth_labelling: dict[Any, int]) -> str:
+    """Return 'up' or 'down' for a triangular face (up: {d, d+1, d+1}; down: {d, d, d+1})."""
+    depths = sorted(depth_labelling[v] for v in face)
+    if depths[0] == depths[1]:
+        return 'down'
+    return 'up'
+
+
+def _level_edge_of_face(face: frozenset[Any], depth_labelling: dict[Any, int]) -> frozenset[Any]:
+    """Return the unique level edge of an up/down triangular face."""
+    vs = list(face)
+    for i in range(3):
+        for j in range(i + 1, 3):
+            if depth_labelling[vs[i]] == depth_labelling[vs[j]]:
+                return frozenset([vs[i], vs[j]])
+    raise ValueError(f"Face {set(face)} has no level edge (not up or down)")
+
+
+def _quad_vertices(quad: frozenset[frozenset[Any]]) -> frozenset[Any]:
+    """Return the 4 vertices of a quadrilateral."""
+    f1, f2 = list(quad)
+    return f1 | f2
+
+
+def _quad_perimeter_edges(
+    quad: frozenset[frozenset[Any]],
+    depth_labelling: dict[Any, int],
+) -> list[frozenset[Any]]:
+    """Return the 4 perimeter edges (non-level) of a quadrilateral."""
+    f1 = next(iter(quad))
+    level_edge = _level_edge_of_face(f1, depth_labelling)
+    edges: list[frozenset[Any]] = []
+    for f in quad:
+        vs = list(f)
+        for i in range(3):
+            edge = frozenset([vs[i], vs[(i + 1) % 3]])
+            if edge != level_edge:
+                edges.append(edge)
+    return edges
+
+
+def _quad_type(quad: frozenset[frozenset[Any]], depth_labelling: dict[Any, int]) -> str:
+    """Return 'shallow_diamond', 'deep_diamond', or 's_quad' for a quadrilateral."""
+    types = tuple(sorted(_triangle_type(f, depth_labelling) for f in quad))
+    if types == ('up', 'up'):
+        return 'shallow_diamond'
+    if types == ('down', 'down'):
+        return 'deep_diamond'
+    return 's_quad'
+
+
+def extended_deep_embedding(
+    g: Graph,
+    outer_cycle: list[Any],
+    plane_depth_labelling: dict[Any, int] | None = None,
+) -> tuple[Graph, list[frozenset[Any]], dict[Any, int], Any, dict[Any, list[Any]]]:
+    """
+    Return the extended deep embedding of g (including the outer face).
+
+    Subdivides every neutral triangular face --- including the outer face,
+    which is always neutral for a maximal planar graph --- by adding a new
+    interior vertex. The vertex added inside the outer face is the
+    outer-cap vertex x*.
+
+    Returns (g_prime, faces, depth_labelling, outer_cap_vertex, embedding) where:
+      - g_prime: the deep embedding graph G'
+      - faces: list of triangular faces of G' (each a frozenset of 3 vertices),
+        derived from Sage's planar embedding of g_prime
+      - depth_labelling: depth of every vertex in G' (extends the original)
+      - outer_cap_vertex: the new vertex placed inside the outer face
+      - embedding: the CCW rotation system of g_prime (vertex -> ordered
+        neighbors) as computed by Sage's planarity routine
+    """
+    if plane_depth_labelling is None:
+        plane_depth_labelling = get_plane_depth_labelling(g, outer_cycle)
+    outer_vertices = frozenset(outer_cycle)
+
+    embedding_g = g.get_embedding()
+    if embedding_g is None:
+        g.is_planar(set_embedding=True)
+        embedding_g = g.get_embedding()
+
+    g_prime = g.copy()
+    next_vertex = max(g.vertices()) + 1
+    depth_labelling = dict(plane_depth_labelling)
+    outer_cap_vertex: Any = None
+
+    for face in g.faces(embedding_g):
+        face_vertices = [u for u, v in face]
+        a, b, c = face_vertices
+        if depth_labelling[a] == depth_labelling[b] == depth_labelling[c]:
+            x = next_vertex
+            next_vertex += 1
+            g_prime.add_vertex(x)
+            g_prime.add_edges([(x, a), (x, b), (x, c)])
+            depth_labelling[x] = depth_labelling[a] + 1
+            if frozenset(face_vertices) == outer_vertices:
+                outer_cap_vertex = x
+
+    g_prime.is_planar(set_embedding=True)
+    embedding = g_prime.get_embedding()
+    assert embedding is not None, "g_prime must be planar after construction"
+    faces = [frozenset([u for u, v in face]) for face in g_prime.faces(embedding)]
+
+    return g_prime, faces, depth_labelling, outer_cap_vertex, embedding
+
+
+def quadrilateral_decomposition(
+    faces: list[frozenset[Any]],
+    depth_labelling: dict[Any, int],
+) -> tuple[list[frozenset[frozenset[Any]]], dict[frozenset[Any], frozenset[frozenset[Any]]]]:
+    """
+    Pair each face with the face on the other side of its level edge.
+
+    Returns (quads, quad_of_face) where:
+      - quads: list of {F1, F2} pairs (each a quadrilateral of the decomposition)
+      - quad_of_face: dict mapping each face to its containing quad
+    """
+    faces_by_edge: dict[frozenset[Any], list[frozenset[Any]]] = {}
+    for face in faces:
+        vs = list(face)
+        for i in range(3):
+            edge = frozenset([vs[i], vs[(i + 1) % 3]])
+            faces_by_edge.setdefault(edge, []).append(face)
+
+    quads: list[frozenset[frozenset[Any]]] = []
+    quad_of_face: dict[frozenset[Any], frozenset[frozenset[Any]]] = {}
+    seen: set[frozenset[Any]] = set()
+    for face in faces:
+        if face in seen:
+            continue
+        level_edge = _level_edge_of_face(face, depth_labelling)
+        adjacent = faces_by_edge[level_edge]
+        assert len(adjacent) == 2, f"Level edge {set(level_edge)} has {len(adjacent)} adjacent faces"
+        other = adjacent[0] if adjacent[1] == face else adjacent[1]
+        quad = frozenset([face, other])
+        quads.append(quad)
+        quad_of_face[face] = quad
+        quad_of_face[other] = quad
+        seen.add(face)
+        seen.add(other)
+
+    return quads, quad_of_face
+
+
+def _boundary_edges(
+    slice_faces: set[frozenset[Any]],
+    faces_by_edge: dict[frozenset[Any], list[frozenset[Any]]],
+) -> set[frozenset[Any]]:
+    """Return the set of edges that lie on the boundary of the slice."""
+    boundary: set[frozenset[Any]] = set()
+    for edge, fs in faces_by_edge.items():
+        inside = sum(1 for f in fs if f in slice_faces)
+        if inside == 1:
+            boundary.add(edge)
+    return boundary
+
+
+def _face_of_wedge(
+    embedding: dict[Any, list[Any]],
+) -> dict[tuple[Any, Any, Any], frozenset[Any]]:
+    """For each (v, n_a, n_b) where n_a, n_b are CCW-consecutive in v's rotation,
+    return the triangular face at that wedge."""
+    fow: dict[tuple[Any, Any, Any], frozenset[Any]] = {}
+    for v, rotation in embedding.items():
+        k = len(rotation)
+        for i in range(k):
+            n_a = rotation[i]
+            n_b = rotation[(i + 1) % k]
+            fow[(v, n_a, n_b)] = frozenset({v, n_a, n_b})
+    return fow
+
+
+def _find_starting_boundary_edge(
+    slice_faces: set[frozenset[Any]],
+    embedding: dict[Any, list[Any]],
+    face_of_wedge: dict[tuple[Any, Any, Any], frozenset[Any]],
+) -> tuple[Any, Any]:
+    """Find a directed boundary edge (v, n_b) with the slice on its LEFT."""
+    for v, rotation in embedding.items():
+        k = len(rotation)
+        for i in range(k):
+            n_a = rotation[i]
+            n_b = rotation[(i + 1) % k]
+            n_c = rotation[(i + 2) % k]
+            left_wedge = face_of_wedge[(v, n_a, n_b)]
+            right_wedge = face_of_wedge[(v, n_b, n_c)]
+            if left_wedge in slice_faces and right_wedge not in slice_faces:
+                return (v, n_b)
+    raise RuntimeError("No boundary edge found for slice")
+
+
+def _boundary_walk(
+    slice_faces: set[frozenset[Any]],
+    embedding: dict[Any, list[Any]],
+    face_of_wedge: dict[tuple[Any, Any, Any], frozenset[Any]],
+) -> list[tuple[Any, Any]]:
+    """Trace the slice's boundary CCW with the slice on the LEFT.
+
+    Returns an ordered list of directed edges forming the closed boundary walk.
+    """
+    start = _find_starting_boundary_edge(slice_faces, embedding, face_of_wedge)
+    walk: list[tuple[Any, Any]] = [start]
+    curr_u, curr_v = start
+    while True:
+        rotation = embedding[curr_v]
+        k = len(rotation)
+        i = rotation.index(curr_u)
+        next_v = None
+        for j in range(1, k + 1):
+            n_jm1 = rotation[(i + j - 1) % k]
+            n_j = rotation[(i + j) % k]
+            if face_of_wedge[(curr_v, n_jm1, n_j)] not in slice_faces:
+                next_v = n_jm1
+                break
+        if next_v is None:
+            raise RuntimeError(f"All wedges at {curr_v} are in slice; no boundary edge")
+        next_edge = (curr_v, next_v)
+        if next_edge == start:
+            break
+        walk.append(next_edge)
+        curr_u, curr_v = next_edge
+    return walk
+
+
+def _canonicalize_walk(
+    walk: list[tuple[Any, Any]],
+    depth_labelling: dict[Any, int],
+) -> list[tuple[Any, Any]]:
+    """Rotate walk to start at the (smallest depth, smallest vertex id) position.
+
+    This fixes the canonical 'top' of the slice for the top-to-bottom scan.
+    """
+    starts = [(depth_labelling[u], u) for u, _ in walk]
+    canon_idx = min(range(len(walk)), key=lambda i: starts[i])
+    return walk[canon_idx:] + walk[:canon_idx]
+
+
+def _attachment_position(
+    quad: frozenset[frozenset[Any]],
+    walk: list[tuple[Any, Any]],
+    depth_labelling: dict[Any, int],
+) -> int:
+    """Return the latest index in walk at which a perimeter edge of quad appears.
+
+    Realizes 'bottommost attachment on the right boundary scanned top-to-bottom'.
+    Returns -1 if the quad has no perimeter edge on the boundary walk.
+    """
+    perimeter = {frozenset(e) for e in _quad_perimeter_edges(quad, depth_labelling)}
+    latest = -1
+    for i, (u, v) in enumerate(walk):
+        if frozenset([u, v]) in perimeter:
+            latest = i
+    return latest
+
+
+def _boundary_deep_diamonds(
+    quads: list[frozenset[frozenset[Any]]],
+    outer_cycle: list[Any],
+    outer_cap_vertex: Any,
+) -> list[frozenset[frozenset[Any]]]:
+    """Return the three boundary deep diamonds (each spans an edge of C)."""
+    outer_set = set(outer_cycle)
+    diamonds: list[frozenset[frozenset[Any]]] = []
+    for q in quads:
+        vs = _quad_vertices(q)
+        if outer_cap_vertex in vs and len(vs & outer_set) == 2:
+            diamonds.append(q)
+    assert len(diamonds) == 3, f"Expected 3 boundary deep diamonds, got {len(diamonds)}"
+    return diamonds
+
+
+def _run_sequence(
+    initial_quad: frozenset[frozenset[Any]],
+    quads: list[frozenset[frozenset[Any]]],
+    depth_labelling: dict[Any, int],
+    faces_by_edge: dict[frozenset[Any], list[frozenset[Any]]],
+    embedding: dict[Any, list[Any]],
+    face_of_wedge: dict[tuple[Any, Any, Any], frozenset[Any]],
+) -> tuple[list[frozenset[frozenset[Any]]], list[int]]:
+    """Run the sequencing loop from initial_quad and return (sequence, move_codes)."""
+    sequence: list[frozenset[frozenset[Any]]] = [initial_quad]
+    move_codes: list[int] = []
+    slice_quads: set[frozenset[frozenset[Any]]] = {initial_quad}
+    slice_faces: set[frozenset[Any]] = set(initial_quad)
+
+    while len(slice_quads) < len(quads):
+        slice_v: set[Any] = set().union(*slice_faces)
+        boundary = _boundary_edges(slice_faces, faces_by_edge)
+        walk = _canonicalize_walk(
+            _boundary_walk(slice_faces, embedding, face_of_wedge),
+            depth_labelling,
+        )
+
+        anchor_drop: list[frozenset[frozenset[Any]]] = []
+        level_add: list[frozenset[frozenset[Any]]] = []
+        join: list[frozenset[frozenset[Any]]] = []
+        ring_completion: list[frozenset[frozenset[Any]]] = []
+
+        for q in quads:
+            if q in slice_quads:
+                continue
+            vs = _quad_vertices(q)
+            k = len(vs & slice_v)
+            edges = _quad_perimeter_edges(q, depth_labelling)
+            j = sum(1 for e in edges if e in boundary)
+            qt = _quad_type(q, depth_labelling)
+
+            if qt == 's_quad' and k == 2 and j == 1:
+                anchor_drop.append(q)
+            if k == 3 and j == 2:
+                level_add.append(q)
+            if qt == 'deep_diamond' and k == 2 and j == 1:
+                join.append(q)
+            if k == 4:
+                ring_completion.append(q)
+
+        def pick(cands: list[frozenset[frozenset[Any]]]) -> frozenset[frozenset[Any]]:
+            return max(cands, key=lambda q: _attachment_position(q, walk, depth_labelling))
+
+        if anchor_drop:
+            next_quad = pick(anchor_drop)
+            next_code = 0
+        elif level_add:
+            next_quad = pick(level_add)
+            next_code = 1
+        elif join:
+            next_quad = pick(join)
+            next_code = 2
+        elif ring_completion:
+            next_quad = pick(ring_completion)
+            next_code = 3
+        else:
+            raise RuntimeError(
+                f"No applicable move at step {len(sequence)}; slice has {len(slice_quads)}/{len(quads)} quads"
+            )
+
+        sequence.append(next_quad)
+        move_codes.append(next_code)
+        slice_quads.add(next_quad)
+        slice_faces.update(next_quad)
+
+    return sequence, move_codes
+
+
+def quadrilateral_sequencing(
+    g: Graph,
+    outer_cycle: list[Any],
+    plane_depth_labelling: dict[Any, int] | None = None,
+) -> QuadrilateralSequence:
+    """
+    Build the quadrilateral sequence of g relative to outer_cycle.
+
+    Constructs the extended deep embedding G' (including the outer-cap vertex
+    x* in the outer face), decomposes G' into quadrilaterals via level-edge
+    pairing, and produces the deterministic sequence Q_1, ..., Q_N by
+    repeatedly applying the move-selection rule:
+
+        anchor drop (0) > level add (1) > join (2) > ring completion (3).
+
+    Within each move's candidate list, the bottommost attachment on the
+    canonical boundary walk is selected (the largest index in the walk where
+    a perimeter edge of the candidate appears). Among the three boundary deep
+    diamonds, Q_1 is the start that produces the lexicographically smallest
+    move-code string.
+
+    Simplification: anchor drop / join detection uses (k, j) = (vertices in
+    slice, perimeter edges on slice boundary) only; the orientation-specific
+    'left edge = right edge' clause is not separately enforced.
+    """
+    g_prime, faces, depth_labelling, outer_cap_vertex, embedding = extended_deep_embedding(
+        g, outer_cycle, plane_depth_labelling
+    )
+    quads, _ = quadrilateral_decomposition(faces, depth_labelling)
+
+    faces_by_edge: dict[frozenset[Any], list[frozenset[Any]]] = {}
+    for face in faces:
+        vs = list(face)
+        for i in range(3):
+            edge = frozenset([vs[i], vs[(i + 1) % 3]])
+            faces_by_edge.setdefault(edge, []).append(face)
+
+    face_of_wedge = _face_of_wedge(embedding)
+    candidates = _boundary_deep_diamonds(quads, outer_cycle, outer_cap_vertex)
+
+    best_sequence: list[frozenset[frozenset[Any]]] | None = None
+    best_codes: list[int] | None = None
+    for q1 in candidates:
+        seq, codes = _run_sequence(q1, quads, depth_labelling, faces_by_edge, embedding, face_of_wedge)
+        if best_codes is None or codes < best_codes:
+            best_sequence = seq
+            best_codes = codes
+    assert best_sequence is not None and best_codes is not None
+
+    return QuadrilateralSequence(
+        deep_embedding=g_prime,
+        triangular_faces=faces,
+        depth_labelling=depth_labelling,
+        outer_cap_vertex=outer_cap_vertex,
+        quadrilaterals=quads,
+        sequence=best_sequence,
+        move_codes=best_codes,
+    )
+
+
+def generate_example(n: int) -> DeeplyEmbeddedGraph:
+    """Generate a random maximal planar graph of size n and return the triangulation, outer cycle, and deep embedding."""
+    g = graphs.RandomTriangulation(n)
+    g.is_planar(set_embedding=True)
+    embedding = g.get_embedding()
+    faces = g.faces(embedding)
+    outer_cycle = [u for u, v in random.choice(faces)]
+    plane_depth_labelling = get_plane_depth_labelling(g, outer_cycle)
+    return DeeplyEmbeddedGraph(graph=g, outer_cycle=outer_cycle, plane_depth_labelling=plane_depth_labelling, deep_embedding=deep_embedding(g, outer_cycle, plane_depth_labelling))
+
+
+if __name__ == "__main__":
+    example = generate_example(10)
+    result = quadrilateral_sequencing(example['graph'], example['outer_cycle'])
+    canonical, graph_dir = canonize_and_save_graph(example['graph'])
+    (graph_dir / "plane_depth_sequence").mkdir(parents=True, exist_ok=True)
+    print(canonical)
+    print(f"Number of quadrilaterals: {len(result['quadrilaterals'])}")
+    print(f"Sequence length: {len(result['sequence'])}")
+    print(f"Move codes: {result['move_codes']}")

papers/plane_depth_sequencing/experiments/plane_depth_sequencing_figure.py

@@ -0,0 +1,221 @@
+"""Generate the labelled quadrilateral sequencing example figure for the paper.
+
+Renders the deep embedding G' of a small maximal planar graph, fills each
+quadrilateral with a color encoding its type, and overlays the index it
+occupies in the canonical sequence Q_1, ..., Q_N. The output is a PDF placed
+next to paper.tex.
+"""
+import argparse
+from pathlib import Path
+from typing import Any
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.patches as patches  # noqa: E402  pylint: disable=wrong-import-position
+import matplotlib.pyplot as plt  # noqa: E402  pylint: disable=wrong-import-position
+
+from sage.all import Graph, graphs  # type: ignore[attr-defined]  # pylint: disable=no-name-in-module,wrong-import-position  # noqa: E402
+from sage.misc.randstate import set_random_seed  # type: ignore[import-not-found]  # noqa: E402
+
+from lib.tutte_embedding import tutte_embedding  # noqa: E402
+from plane_depth_sequencing import (  # noqa: E402
+    _level_edge_of_face,
+    _quad_type,
+    _quad_vertices,
+    quadrilateral_sequencing,
+)
+
+
+def _planar_pos(
+    g_prime: Graph,
+    outer_face: list[Any],
+    layout: str,
+) -> dict[Any, tuple[float, float]]:
+    """Compute positions for g_prime under either 'tutte' or 'sage_planar' layout."""
+    if layout == "tutte":
+        return tutte_embedding(g_prime, outer_face)
+    if layout == "sage_planar":
+        g_prime.is_planar(set_embedding=True, set_pos=True)
+        pos = g_prime.get_pos()
+        assert pos is not None, "Sage failed to set planar positions"
+        return {v: (float(p[0]), float(p[1])) for v, p in pos.items()}
+    raise ValueError(f"Unknown layout: {layout}")
+
+
+MOVE_NAMES = {0: "AD", 1: "LA", 2: "J", 3: "RC"}
+TYPE_COLORS = {
+    "shallow_diamond": "#FFE0B2",
+    "deep_diamond": "#B2DFDB",
+    "s_quad": "#F8BBD0",
+}
+
+
+def _ordered_quad_vertices(quad: frozenset, depth_labelling: dict[Any, int]) -> list[Any]:
+    """Return the 4 quad vertices in cyclic order: e1, a, e2, b."""
+    f1, f2 = list(quad)
+    level_edge = _level_edge_of_face(f1, depth_labelling)
+    e1, e2 = list(level_edge)
+    a = next(v for v in f1 if v not in level_edge)
+    b = next(v for v in f2 if v not in level_edge)
+    return [e1, a, e2, b]
+
+
+def _pick_outer_cap_face(
+    g_prime: Graph,
+    outer_cap_vertex: Any,
+    outer_cycle: list[Any],
+) -> list[Any]:
+    """Return the vertex list of an outer-cap face (the face used as outer in the plane drawing)."""
+    g_prime.is_planar(set_embedding=True)
+    embedding = g_prime.get_embedding()
+    outer_set = set(outer_cycle)
+    for face in g_prime.faces(embedding):
+        verts = [u for u, _ in face]
+        if outer_cap_vertex in verts and len(set(verts) & outer_set) == 2:
+            return verts
+    raise RuntimeError("No outer-cap face found in G' embedding")
+
+
+def _draw_figure(
+    g: Graph,
+    outer_cycle: list[Any],
+    out_path: Path,
+    layout: str = "sage_planar",
+) -> dict[str, Any]:
+    """Render and save the labelled sequencing figure. Returns summary stats."""
+    result = quadrilateral_sequencing(g, outer_cycle)
+    g_prime: Graph = result["deep_embedding"]
+    depth_labelling: dict[Any, int] = result["depth_labelling"]
+    sequence = result["sequence"]
+    move_codes = result["move_codes"]
+    outer_cap_vertex = result["outer_cap_vertex"]
+
+    outer_face = _pick_outer_cap_face(g_prime, outer_cap_vertex, outer_cycle)
+    pos = _planar_pos(g_prime, outer_face, layout)
+
+    fig, ax = plt.subplots(figsize=(9, 9))
+
+    for i, quad in enumerate(sequence):
+        ordered = _ordered_quad_vertices(quad, depth_labelling)
+        poly_pts = [pos[v] for v in ordered]
+        qt = _quad_type(quad, depth_labelling)
+        polygon = patches.Polygon(
+            poly_pts,
+            closed=True,
+            facecolor=TYPE_COLORS[qt],
+            edgecolor="none",
+            alpha=0.65,
+            zorder=1,
+        )
+        ax.add_patch(polygon)
+
+    for u, v in g_prime.edges(labels=False):
+        x1, y1 = pos[u]
+        x2, y2 = pos[v]
+        if depth_labelling[u] == depth_labelling[v]:
+            ax.plot([x1, x2], [y1, y2], linestyle=(0, (3, 2)), color="#666666", linewidth=1.0, zorder=2)
+        else:
+            ax.plot([x1, x2], [y1, y2], "-", color="black", linewidth=1.0, zorder=2)
+
+    outer_set = set(outer_cycle)
+    for v, (x, y) in pos.items():
+        if v == outer_cap_vertex:
+            color = "#D32F2F"
+        elif v in outer_set:
+            color = "#1976D2"
+        else:
+            color = "black"
+        ax.scatter([x], [y], s=28, color=color, zorder=5, edgecolors="white", linewidths=0.8)
+        name = "$x^{*}$" if v == outer_cap_vertex else str(v)
+        ax.annotate(
+            name,
+            (x, y),
+            xytext=(7, 6),
+            textcoords="offset points",
+            fontsize=9,
+            zorder=8,
+            color="#222222",
+            bbox={"boxstyle": "round,pad=0.05", "facecolor": "white", "edgecolor": "none", "alpha": 0.85},
+        )
+
+    for i, quad in enumerate(sequence):
+        ordered = _ordered_quad_vertices(quad, depth_labelling)
+        poly_pts = [pos[v] for v in ordered]
+        cx = sum(p[0] for p in poly_pts) / 4
+        cy = sum(p[1] for p in poly_pts) / 4
+        if i == 0:
+            label = "$Q_{1}$"
+        else:
+            label = f"$Q_{{{i + 1}}}^{{\\mathrm{{{MOVE_NAMES[move_codes[i - 1]]}}}}}$"
+        ax.text(
+            cx,
+            cy,
+            label,
+            ha="center",
+            va="center",
+            fontsize=12,
+            zorder=7,
+            bbox={"boxstyle": "round,pad=0.22", "facecolor": "white", "edgecolor": "#444444", "alpha": 0.95},
+        )
+
+    legend_handles = [
+        patches.Patch(color=TYPE_COLORS["shallow_diamond"], label="shallow diamond"),
+        patches.Patch(color=TYPE_COLORS["deep_diamond"], label="deep diamond"),
+        patches.Patch(color=TYPE_COLORS["s_quad"], label="S quad"),
+    ]
+    ax.legend(handles=legend_handles, loc="lower right", fontsize=10, framealpha=0.95)
+
+    ax.set_aspect("equal")
+    ax.axis("off")
+    plt.tight_layout()
+    out_path.parent.mkdir(parents=True, exist_ok=True)
+    plt.savefig(out_path, format="pdf", bbox_inches="tight")
+    plt.close(fig)
+
+    return {
+        "n_vertices_g": g.order(),
+        "n_vertices_gprime": g_prime.order(),
+        "n_quads": len(sequence),
+        "move_codes": move_codes,
+        "outer_cycle": outer_cycle,
+        "outer_cap_vertex": outer_cap_vertex,
+        "depth_labelling": depth_labelling,
+    }
+
+
+def _build_example(seed: int, n: int) -> tuple[Graph, list[Any]]:
+    """Build a random triangulation and pick a deterministic outer cycle."""
+    set_random_seed(seed)
+    g = graphs.RandomTriangulation(n)
+    g.is_planar(set_embedding=True)
+    embedding = g.get_embedding()
+    faces = g.faces(embedding)
+    outer_cycle = [u for u, _ in faces[0]]
+    return g, outer_cycle
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--seed", type=int, default=141)
+    parser.add_argument("--n", type=int, default=9)
+    parser.add_argument(
+        "--out",
+        type=Path,
+        default=Path("papers/plane_depth_sequencing/example_figure.pdf"),
+    )
+    parser.add_argument("--layout", choices=("tutte", "sage_planar"), default="sage_planar")
+    args = parser.parse_args()
+
+    g, outer_cycle = _build_example(args.seed, args.n)
+    summary = _draw_figure(g, outer_cycle, args.out, layout=args.layout)
+
+    print(f"Wrote {args.out}")
+    print(f"  |V(G)|={summary['n_vertices_g']}, |V(G')|={summary['n_vertices_gprime']}, "
+          f"N={summary['n_quads']}")
+    print(f"  outer cycle: {summary['outer_cycle']}, x* = {summary['outer_cap_vertex']}")
+    code_str = "".join(str(c) for c in summary["move_codes"])
+    print(f"  move-code string: {code_str}")
+
+
+if __name__ == "__main__":
+    main()

papers/plane_depth_sequencing/experiments/quad_sequence_coloring_check.py

@@ -0,0 +1,155 @@
+"""Empirical check: does pure greedy online 4-coloring succeed under the canonical quadrilateral sequencing?
+
+The user's conjecture-of-the-day: for the canonical sequence, every non-ring-completion
+move (anchor drop, level add, join) admits a local 4-coloring choice, and ring completions
+introduce no new vertices, so the only question is whether the existing coloring agrees
+on the ring-completion quad's edges.
+
+Under proper greedy online coloring this latter check is automatic: every G'-edge
+constraint is enforced as the second endpoint is colored. So the empirical question
+becomes: does pure greedy (smallest color avoiding already-colored G'-neighbors) ever
+get stuck at a non-ring-completion move with all four colors forbidden?
+
+We run greedy in a canonical vertex order (depth, then vertex id) and classify every
+failure as one of:
+  - 'no_color_available': a new vertex saw all 4 colors among colored G'-neighbors
+  - 'ring_completion_monochromatic': a ring completion's quad has a monochromatic edge
+  - 'success': global proper 4-coloring of G'.
+"""
+from typing import Any
+from sage.all import Graph, graphs  # type: ignore[attr-defined]  # pylint: disable=no-name-in-module
+from lib.colored_graphs import canonize_and_save_graph
+from plane_depth_sequencing import (
+    quadrilateral_sequencing,
+    _quad_vertices,
+)
+
+
+def greedy_4_color_under_sequence(
+    g_prime: Graph,
+    sequence: list[frozenset[frozenset[Any]]],
+    depth_labelling: dict[Any, int],
+) -> tuple[dict[Any, int], str, int | None]:
+    """Run greedy online 4-coloring under the canonical sequence.
+
+    Returns (partial_coloring, status, failing_step).
+    status is one of:
+      - 'success'
+      - 'no_color_available' (failing_step is the sequence index)
+      - 'ring_completion_monochromatic' (failing_step is the sequence index)
+    """
+    coloring: dict[Any, int] = {}
+    slice_vertices: set[Any] = set()
+    adj: dict[Any, set[Any]] = {v: set(g_prime.neighbors(v)) for v in g_prime.vertices()}
+
+    for step, quad in enumerate(sequence):
+        quad_v = _quad_vertices(quad)
+        new_vertices = quad_v - slice_vertices
+
+        if not new_vertices:
+            # Ring completion: verify the existing coloring is proper on the quad's edges.
+            quad_v_list = list(quad_v)
+            for i in range(len(quad_v_list)):
+                for j in range(i + 1, len(quad_v_list)):
+                    u, v = quad_v_list[i], quad_v_list[j]
+                    if v in adj[u] and coloring[u] == coloring[v]:
+                        return coloring, 'ring_completion_monochromatic', step
+            continue
+
+        # Color new vertices in a canonical order (depth, then vertex id).
+        for v in sorted(new_vertices, key=lambda x: (depth_labelling[x], x)):
+            used = {coloring[u] for u in adj[v] if u in coloring}
+            avail = [c for c in (0, 1, 2, 3) if c not in used]
+            if not avail:
+                return coloring, 'no_color_available', step
+            coloring[v] = min(avail)
+            slice_vertices.add(v)
+
+    return coloring, 'success', None
+
+
+def _outer_cycle_from_embedding(g: Graph) -> list[Any]:
+    """Pick a deterministic outer cycle: the first face Sage returns from the embedding."""
+    g.is_planar(set_embedding=True)
+    embedding = g.get_embedding()
+    faces = g.faces(embedding)
+    return [u for u, v in faces[0]]
+
+
+def check_graph(g: Graph) -> dict[str, Any]:
+    """Run sequencing + greedy coloring on a single graph. Returns a result dict."""
+    outer_cycle = _outer_cycle_from_embedding(g)
+    try:
+        seq_result = quadrilateral_sequencing(g, outer_cycle)
+    except Exception as exc:  # pylint: disable=broad-except
+        return {
+            'graph6': g.graph6_string(),
+            'status': 'sequencing_error',
+            'error': str(exc),
+        }
+    coloring, status, failing_step = greedy_4_color_under_sequence(
+        seq_result['deep_embedding'],
+        seq_result['sequence'],
+        seq_result['depth_labelling'],
+    )
+    return {
+        'graph6': g.graph6_string(),
+        'status': status,
+        'failing_step': failing_step,
+        'coloring_size': len(coloring),
+        'num_quads': len(seq_result['quadrilaterals']),
+        'move_codes': seq_result['move_codes'],
+    }
+
+
+def run_enumeration(min_order: int, max_order: int) -> dict[int, dict[str, Any]]:
+    """Iterate over all 3-connected triangulations of order in [min_order, max_order]."""
+    summary: dict[int, dict[str, Any]] = {}
+    for n in range(min_order, max_order + 1):
+        total = 0
+        by_status: dict[str, int] = {}
+        failures: list[dict[str, Any]] = []
+        for g in graphs.planar_graphs(n, minimum_connectivity=3, maximum_face_size=3):
+            total += 1
+            result = check_graph(g)
+            status = result['status']
+            by_status[status] = by_status.get(status, 0) + 1
+            if status != 'success':
+                failures.append(result)
+            if total % 200 == 0:
+                print(f"  n={n}: checked {total}, statuses so far: {by_status}")
+        summary[n] = {'total': total, 'by_status': by_status, 'failures': failures}
+        print(f"n={n}: total={total}, by_status={by_status}")
+        for f in failures[:5]:
+            print(f"  failure {f['graph6']}: {f['status']} at step {f['failing_step']}")
+        if len(failures) > 5:
+            print(f"  ... ({len(failures) - 5} more failures)")
+    return summary
+
+
+if __name__ == "__main__":
+    import sys
+
+    min_order = int(sys.argv[1]) if len(sys.argv) > 1 else 4
+    max_order = int(sys.argv[2]) if len(sys.argv) > 2 else 9
+
+    summary = run_enumeration(min_order, max_order)
+
+    print()
+    print("=" * 60)
+    print("Final summary:")
+    for n, s in summary.items():
+        print(f"  n={n}: total={s['total']}, by_status={s['by_status']}")
+
+    # Save up to first 10 failure graphs per order for inspection.
+    saved = 0
+    for n, s in summary.items():
+        for f in s['failures'][:10]:
+            try:
+                g = Graph(f['graph6'])
+                canonical, graph_dir = canonize_and_save_graph(g)
+                print(f"  saved n={n} {f['graph6']} ({f['status']}) -> {graph_dir}")
+                saved += 1
+            except Exception as exc:  # pylint: disable=broad-except
+                print(f"  failed to save {f['graph6']}: {exc}")
+    print(f"Saved {saved} failure graphs total.")