Draw random medial tire decompositions

papers/medial_tire_decompositions_of_plane_triangulations/experiments/draw_random_medial_tire_decompositions.py

@@ -0,0 +1,480 @@
+"""Draw medial tire decompositions of random 5-connected triangulations.
+
+The source graphs come from ``plantri -c5`` in graph6 format.  For each sampled
+30-vertex triangulation, this script chooses a random source vertex, builds the
+BFS depth-component tire tree, recognizes every full medial tire graph in the
+decomposition, and draws both the tire tree and the realized full medial tire
+graphs.
+"""
+
+from __future__ import annotations
+
+import argparse
+import math
+import os
+import random
+import subprocess
+import sys
+from collections import defaultdict
+from dataclasses import dataclass
+from pathlib import Path
+
+PAPER_DIR = Path(__file__).resolve().parents[1]
+REPO_ROOT = PAPER_DIR.parents[1]
+os.environ.setdefault(
+    "MPLCONFIGDIR", str(PAPER_DIR / "experiments" / ".matplotlib-cache")
+)
+os.environ.setdefault("XDG_CACHE_HOME", str(PAPER_DIR / "experiments" / ".cache"))
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
+import networkx as nx
+
+if str(PAPER_DIR) not in sys.path:
+    sys.path.insert(0, str(PAPER_DIR))
+
+from lib.medial_tire_decomposition import ekey, medial_tire_facemodel, recognise
+from lib.full_medial_tire_generator import FullMedialTireGraph
+
+
+@dataclass(frozen=True)
+class TreadNode:
+    idx: int
+    depth: int
+    face_indices: tuple[int, ...]
+    annular: frozenset
+    up: frozenset
+    down: frozenset
+    bites: frozenset
+    tires: tuple[tuple[FullMedialTireGraph, dict], ...]
+
+
+def sample_plantri_graphs(n: int, count: int, seed: int, scan_limit: int) -> list[nx.Graph]:
+    cmd = [str(REPO_ROOT / "plantri" / "plantri"), "-g", "-c5", str(n)]
+    rng = random.Random(seed)
+    sample: list[tuple[int, nx.Graph]] = []
+    seen = 0
+    with subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as proc:
+        assert proc.stdout is not None
+        for raw in proc.stdout:
+            line = raw.strip()
+            if not line or line.startswith(b">>"):
+                continue
+            graph = nx.from_graph6_bytes(line)
+            if nx.node_connectivity(graph) < 5:
+                continue
+            seen += 1
+            if len(sample) < count:
+                sample.append((seen, graph))
+            else:
+                j = rng.randrange(seen)
+                if j < count:
+                    sample[j] = (seen, graph)
+            if seen >= scan_limit:
+                proc.terminate()
+                break
+        proc.wait(timeout=10)
+    if len(sample) < count:
+        raise RuntimeError(f"only found {len(sample)} graphs after scanning {seen}")
+    return [graph for _ordinal, graph in sample]
+
+
+def triangular_faces(g: nx.Graph):
+    ok, emb = nx.check_planarity(g)
+    if not ok:
+        raise ValueError("not planar")
+    seen = set()
+    faces = []
+    for u, v in list(emb.edges()):
+        if (u, v) in seen:
+            continue
+        face = tuple(emb.traverse_face(u, v, mark_half_edges=seen))
+        faces.append(face)
+    return faces
+
+
+def edge_face_data(faces):
+    face_edges = []
+    edge_faces: dict[tuple, list[int]] = defaultdict(list)
+    for i, face in enumerate(faces):
+        edges = {
+            ekey(face[0], face[1]),
+            ekey(face[1], face[2]),
+            ekey(face[2], face[0]),
+        }
+        face_edges.append(edges)
+        for edge in edges:
+            edge_faces[edge].append(i)
+    return face_edges, edge_faces
+
+
+def depth_components(faces, face_edges, edge_faces, levels):
+    depths = [min(levels[v] for v in face) for face in faces]
+    dual_adj: dict[int, set[int]] = defaultdict(set)
+    for incident in edge_faces.values():
+        for a in range(len(incident)):
+            for b in range(a + 1, len(incident)):
+                dual_adj[incident[a]].add(incident[b])
+                dual_adj[incident[b]].add(incident[a])
+
+    comps = []
+    seen = [False] * len(faces)
+    for start in range(len(faces)):
+        if seen[start]:
+            continue
+        depth = depths[start]
+        stack = [start]
+        comp = []
+        seen[start] = True
+        while stack:
+            face = stack.pop()
+            comp.append(face)
+            for other in dual_adj[face]:
+                if not seen[other] and depths[other] == depth:
+                    seen[other] = True
+                    stack.append(other)
+        comps.append((depth, tuple(sorted(comp))))
+    return comps, depths, dual_adj
+
+
+def tread_from_component(faces, levels, face_indices):
+    tread_faces = [faces[i] for i in face_indices]
+    if not tread_faces:
+        return None
+    depth = min(min(levels[v] for v in face) for face in tread_faces)
+    annular, up, down = set(), set(), set()
+    face_of_down = defaultdict(int)
+    for face in tread_faces:
+        for x, y in ((face[0], face[1]), (face[1], face[2]), (face[2], face[0])):
+            e = ekey(x, y)
+            lx, ly = levels[x], levels[y]
+            if {lx, ly} == {depth, depth + 1}:
+                annular.add(e)
+            elif lx == ly == depth:
+                up.add(e)
+            elif lx == ly == depth + 1:
+                down.add(e)
+                face_of_down[e] += 1
+    if len(annular) < 3:
+        return None
+    return {
+        "tread_faces": tread_faces,
+        "annular": annular,
+        "up": up,
+        "down": down,
+        "bites": {e for e in down if face_of_down[e] == 2},
+    }
+
+
+def build_tire_tree(g: nx.Graph, source: int):
+    faces = triangular_faces(g)
+    face_edges, edge_faces = edge_face_data(faces)
+    levels = nx.single_source_shortest_path_length(g, source)
+    comps, depths, dual_adj = depth_components(faces, face_edges, edge_faces, levels)
+    comp_of_face = {}
+    for comp_idx, (_depth, face_indices) in enumerate(comps):
+        for face in face_indices:
+            comp_of_face[face] = comp_idx
+
+    nodes: list[TreadNode] = []
+    comp_to_node = {}
+    for comp_idx, (depth, face_indices) in enumerate(comps):
+        tread = tread_from_component(faces, levels, face_indices)
+        if tread is None or len(tread["up"]) < 3:
+            continue
+        mt = medial_tire_facemodel(tread["tread_faces"])
+        tires = tuple(recognise(mt, tread))
+        if not tires:
+            continue
+        node = TreadNode(
+            idx=len(nodes),
+            depth=depth,
+            face_indices=face_indices,
+            annular=frozenset(tread["annular"]),
+            up=frozenset(tread["up"]),
+            down=frozenset(tread["down"]),
+            bites=frozenset(tread["bites"]),
+            tires=tires,
+        )
+        comp_to_node[comp_idx] = node.idx
+        nodes.append(node)
+
+    tree_edges = set()
+    for comp_idx, (depth, face_indices) in enumerate(comps):
+        if comp_idx not in comp_to_node:
+            continue
+        child = comp_to_node[comp_idx]
+        parent_candidates = set()
+        for face in face_indices:
+            for other in dual_adj[face]:
+                other_comp = comp_of_face[other]
+                if depths[other] == depth - 1 and other_comp in comp_to_node:
+                    parent_candidates.add(comp_to_node[other_comp])
+        for parent in parent_candidates:
+            tree_edges.add((parent, child))
+    return faces, levels, nodes, sorted(tree_edges)
+
+
+def graph_layout(g: nx.Graph):
+    try:
+        return nx.planar_layout(g)
+    except nx.NetworkXException:
+        return nx.spring_layout(g, seed=0)
+
+
+def draw_base_graph(ax, g, levels, source):
+    pos = graph_layout(g)
+    max_level = max(levels.values())
+    cmap = plt.get_cmap("viridis", max_level + 1)
+    node_colors = [cmap(levels[v]) for v in g.nodes()]
+    nx.draw_networkx_edges(g, pos, ax=ax, edge_color="#cbd5e1", width=0.8)
+    nx.draw_networkx_nodes(
+        g,
+        pos,
+        ax=ax,
+        node_color=node_colors,
+        node_size=[150 if v == source else 72 for v in g.nodes()],
+        edgecolors=["#dc2626" if v == source else "#111827" for v in g.nodes()],
+        linewidths=[1.8 if v == source else 0.45 for v in g.nodes()],
+    )
+    labels = {v: str(v) for v in g.nodes()}
+    nx.draw_networkx_labels(g, pos, labels=labels, ax=ax, font_size=5)
+    ax.set_title(f"G, source {source}; vertex levels 0..{max_level}", fontsize=10)
+    ax.set_aspect("equal")
+    ax.axis("off")
+
+
+def tree_positions(nodes: list[TreadNode], tree_edges):
+    children: dict[int, list[int]] = defaultdict(list)
+    has_parent = set()
+    for parent, child in tree_edges:
+        children[parent].append(child)
+        has_parent.add(child)
+    roots = [node.idx for node in nodes if node.idx not in has_parent]
+    for child_list in children.values():
+        child_list.sort(key=lambda idx: (nodes[idx].depth, idx))
+
+    x_counter = 0
+    pos = {}
+
+    def place(idx, depth):
+        nonlocal x_counter
+        if not children[idx]:
+            pos[idx] = (x_counter, -depth)
+            x_counter += 1
+            return pos[idx][0]
+        xs = [place(child, depth + 1) for child in children[idx]]
+        x = sum(xs) / len(xs)
+        pos[idx] = (x, -depth)
+        return x
+
+    for root in sorted(roots, key=lambda idx: (nodes[idx].depth, idx)):
+        place(root, 0)
+        x_counter += 1
+    return pos
+
+
+def draw_tire_tree(ax, nodes: list[TreadNode], tree_edges):
+    pos = tree_positions(nodes, tree_edges)
+    for parent, child in tree_edges:
+        x0, y0 = pos[parent]
+        x1, y1 = pos[child]
+        ax.plot([x0, x1], [y0, y1], color="#374151", lw=1.0, zorder=1)
+    for node in nodes:
+        x, y = pos[node.idx]
+        ax.text(
+            x,
+            y,
+            f"T{node.idx}\nd={node.depth}\n{len(node.tires)} tire(s)",
+            ha="center",
+            va="center",
+            fontsize=8,
+            bbox={
+                "boxstyle": "round,pad=0.32",
+                "facecolor": "#fef3c7",
+                "edgecolor": "#111827",
+                "linewidth": 0.9,
+            },
+            zorder=3,
+        )
+    ax.set_title("Depth-component tire tree", fontsize=10)
+    if pos:
+        xs = [p[0] for p in pos.values()]
+        ys = [p[1] for p in pos.values()]
+        ax.set_xlim(min(xs) - 1.0, max(xs) + 1.0)
+        ax.set_ylim(min(ys) - 0.7, max(ys) + 0.7)
+    ax.axis("off")
+
+
+def vertex_xy(k: int, n: int, radius: float) -> tuple[float, float]:
+    angle = math.pi / 2 - 2 * math.pi * k / n
+    return radius * math.cos(angle), radius * math.sin(angle)
+
+
+def edge_midpoint_angle(i: int, n: int) -> float:
+    return math.pi / 2 - 2 * math.pi * (i + 0.5) / n
+
+
+def draw_full_medial_tire(ax, graph: FullMedialTireGraph, title: str):
+    n = graph.n
+    ann = [vertex_xy(k, n, 1.0) for k in range(n)]
+    matched = graph.bite_edges
+    cyc_x = [p[0] for p in ann] + [ann[0][0]]
+    cyc_y = [p[1] for p in ann] + [ann[0][1]]
+    ax.plot(cyc_x, cyc_y, color="black", lw=1.3, zorder=2)
+    for i, tooth in enumerate(graph.tooth_word):
+        if tooth == "U":
+            r, color = 1.42, "#2563eb"
+        elif i not in matched:
+            r, color = 0.58, "#dc2626"
+        else:
+            continue
+        ang = edge_midpoint_angle(i, n)
+        apex = (r * math.cos(ang), r * math.sin(ang))
+        for corner in (ann[i], ann[(i + 1) % n]):
+            ax.plot([apex[0], corner[0]], [apex[1], corner[1]], color="#9ca3af", lw=0.5)
+        ax.scatter([apex[0]], [apex[1]], s=12, color=color, zorder=3)
+    for i, j in sorted(graph.bites):
+        corners = [ann[i], ann[(i + 1) % n], ann[j], ann[(j + 1) % n]]
+        apex = (0.82 * sum(p[0] for p in corners) / 4, 0.82 * sum(p[1] for p in corners) / 4)
+        for corner in corners:
+            ax.plot([apex[0], corner[0]], [apex[1], corner[1]], color="#9ca3af", lw=0.5)
+        ax.scatter([apex[0]], [apex[1]], s=22, color="#7f1d1d", edgecolors="black", lw=0.4)
+    ax.scatter([p[0] for p in ann], [p[1] for p in ann], s=9, color="black", zorder=4)
+    bites = ",".join(f"{i}{j}" for i, j in sorted(graph.bites)) or "-"
+    ax.set_title(f"{title}\n{graph.tooth_word} b:{bites}", fontsize=5.8, pad=1.5)
+    ax.set_xlim(-1.6, 1.6)
+    ax.set_ylim(-1.6, 1.6)
+    ax.set_aspect("equal")
+    ax.axis("off")
+
+
+def draw_medial_tire_grid(fig, outer_spec, nodes):
+    tires = []
+    for node in nodes:
+        for comp_idx, (graph, _bij) in enumerate(node.tires):
+            tires.append((node.idx, comp_idx, graph))
+    if not tires:
+        ax = fig.add_subplot(outer_spec)
+        ax.text(0.5, 0.5, "No full medial tire graphs recognized", ha="center")
+        ax.axis("off")
+        return
+    cols = min(5, max(1, math.ceil(math.sqrt(len(tires)))))
+    rows = math.ceil(len(tires) / cols)
+    sub = outer_spec.subgridspec(rows, cols, wspace=0.08, hspace=0.35)
+    for i in range(rows * cols):
+        ax = fig.add_subplot(sub[i // cols, i % cols])
+        if i < len(tires):
+            node_idx, comp_idx, graph = tires[i]
+            draw_full_medial_tire(ax, graph, f"T{node_idx}.{comp_idx} n={graph.n}")
+        else:
+            ax.axis("off")
+
+
+def write_index(path: Path, graph_idx: int, source: int, g: nx.Graph, nodes, tree_edges):
+    lines = [
+        f"# Random medial tire decomposition {graph_idx}",
+        "",
+        f"- vertices: {g.number_of_nodes()}",
+        f"- edges: {g.number_of_edges()}",
+        f"- node connectivity: {nx.node_connectivity(g)}",
+        f"- source vertex: {source}",
+        f"- tire-tree nodes: {len(nodes)}",
+        f"- tire-tree edges: {len(tree_edges)}",
+        "",
+        "| node | depth | faces | annular | up | down | bites | full medial tires |",
+        "|--:|--:|--:|--:|--:|--:|--:|--:|",
+    ]
+    for node in nodes:
+        lines.append(
+            f"| T{node.idx} | {node.depth} | {len(node.face_indices)} | "
+            f"{len(node.annular)} | {len(node.up)} | {len(node.down)} | "
+            f"{len(node.bites)} | {len(node.tires)} |"
+        )
+        for comp_idx, (tire, _bij) in enumerate(node.tires):
+            bites = ",".join(f"({i},{j})" for i, j in sorted(tire.bites)) or "-"
+            lines.append(
+                f"| T{node.idx}.{comp_idx} |  |  | | {len(tire.up_edges)} | "
+                f"{len(tire.down_edges)} | {bites} | `{tire.tooth_word}` |"
+            )
+    path.write_text("\n".join(lines) + "\n")
+
+
+def draw_case(out_dir: Path, graph_idx: int, g: nx.Graph, source: int):
+    _faces, levels, nodes, tree_edges = build_tire_tree(g, source)
+    fig = plt.figure(figsize=(17, 10))
+    spec = fig.add_gridspec(2, 2, width_ratios=[1.15, 1.0], height_ratios=[1.0, 1.25])
+    ax_graph = fig.add_subplot(spec[0, 0])
+    ax_tree = fig.add_subplot(spec[1, 0])
+    draw_base_graph(ax_graph, g, levels, source)
+    draw_tire_tree(ax_tree, nodes, tree_edges)
+    draw_medial_tire_grid(fig, spec[:, 1], nodes)
+    fig.suptitle(
+        f"Random 5-connected maximal planar graph {graph_idx}: "
+        f"n={g.number_of_nodes()}, source={source}",
+        fontsize=13,
+    )
+    legend = [
+        Line2D([0], [0], marker="o", color="w", label="source",
+               markerfacecolor="#fde68a", markeredgecolor="#dc2626", markersize=8),
+        Line2D([0], [0], color="black", lw=1.3, label="annular cycle A(T)"),
+        Line2D([0], [0], marker="o", color="w", label="up tooth",
+               markerfacecolor="#2563eb", markersize=6),
+        Line2D([0], [0], marker="o", color="w", label="down tooth",
+               markerfacecolor="#dc2626", markersize=6),
+        Line2D([0], [0], marker="o", color="w", label="bite apex",
+               markerfacecolor="#7f1d1d", markeredgecolor="black", markersize=6),
+    ]
+    fig.legend(handles=legend, loc="lower center", ncol=5, fontsize=9)
+    fig.subplots_adjust(left=0.03, right=0.99, top=0.92, bottom=0.08, wspace=0.08, hspace=0.16)
+
+    png = out_dir / f"random_c5_n30_medial_tire_decomposition_{graph_idx}.png"
+    pdf = out_dir / f"random_c5_n30_medial_tire_decomposition_{graph_idx}.pdf"
+    fig.savefig(png, dpi=180)
+    fig.savefig(pdf)
+    plt.close(fig)
+    write_index(
+        out_dir / f"random_c5_n30_medial_tire_decomposition_{graph_idx}.md",
+        graph_idx,
+        source,
+        g,
+        nodes,
+        tree_edges,
+    )
+    return png, pdf, len(nodes), sum(len(node.tires) for node in nodes)
+
+
+def run(args: argparse.Namespace):
+    out_dir = Path(args.out_dir)
+    out_dir.mkdir(parents=True, exist_ok=True)
+    graphs = sample_plantri_graphs(args.n, args.count, args.seed, args.scan_limit)
+    rng = random.Random(args.seed + 101)
+    for i, graph in enumerate(graphs, start=1):
+        source = rng.choice(list(graph.nodes()))
+        png, pdf, node_count, tire_count = draw_case(out_dir, i, graph, source)
+        print(
+            f"case {i}: source={source}, connectivity={nx.node_connectivity(graph)}, "
+            f"tire nodes={node_count}, full medial tires={tire_count}"
+        )
+        print(f"  wrote {png}")
+        print(f"  wrote {pdf}")
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--n", type=int, default=30)
+    parser.add_argument("--count", type=int, default=2)
+    parser.add_argument("--seed", type=int, default=20260615)
+    parser.add_argument("--scan-limit", type=int, default=500)
+    parser.add_argument(
+        "--out-dir",
+        default=str(PAPER_DIR / "experiments" / "random_medial_tire_decompositions"),
+    )
+    run(parser.parse_args())
+
+
+if __name__ == "__main__":
+    main()

papers/medial_tire_decompositions_of_plane_triangulations/experiments/random_medial_tire_decompositions/random_c5_n30_medial_tire_decomposition_1.md

@@ -0,0 +1,18 @@
+# Random medial tire decomposition 1
+
+- vertices: 30
+- edges: 84
+- node connectivity: 5
+- source vertex: 9
+- tire-tree nodes: 3
+- tire-tree edges: 2
+
+| node | depth | faces | annular | up | down | bites | full medial tires |
+|--:|--:|--:|--:|--:|--:|--:|--:|
+| T0 | 1 | 16 | 16 | 6 | 10 | 0 | 1 |
+| T0.0 |  |  | | 6 | 10 | - | `DUDUDUDDUDDDUDDU` |
+| T1 | 2 | 20 | 20 | 10 | 10 | 0 | 1 |
+| T1.0 |  |  | | 10 | 10 | - | `UUDUUDUDDUDUDUDUUDDD` |
+| T2 | 3 | 14 | 13 | 12 | 1 | 1 | 2 |
+| T2.0 |  |  | | 6 | 2 | (1,5) | `UDUUUDUU` |
+| T2.1 |  |  | | 5 | 0 | - | `UUUUU` |

papers/medial_tire_decompositions_of_plane_triangulations/experiments/random_medial_tire_decompositions/random_c5_n30_medial_tire_decomposition_2.md

@@ -0,0 +1,17 @@
+# Random medial tire decomposition 2
+
+- vertices: 30
+- edges: 84
+- node connectivity: 5
+- source vertex: 4
+- tire-tree nodes: 3
+- tire-tree edges: 2
+
+| node | depth | faces | annular | up | down | bites | full medial tires |
+|--:|--:|--:|--:|--:|--:|--:|--:|
+| T0 | 1 | 16 | 16 | 7 | 9 | 0 | 1 |
+| T0.0 |  |  | | 7 | 9 | - | `DUDUDUDUDUDDUDUD` |
+| T1 | 2 | 17 | 17 | 9 | 8 | 0 | 1 |
+| T1.0 |  |  | | 9 | 8 | - | `UUUDDUDUDUDUDDUUD` |
+| T2 | 3 | 14 | 14 | 8 | 5 | 1 | 1 |
+| T2.0 |  |  | | 8 | 6 | (1,5) | `DDUUUDDUUDUDUU` |