Add level resolutions of maximal planar graphs paper

Migrate the paper content into the amsart template and include the
supporting experiments scripts.

papers/level_resolutions_of_maximal_planar_graphs/experiments/README.md

@@ -0,0 +1,58 @@
+# Level Resolution Experiments
+
+Computational investigation of a structural proof strategy for the four
+color theorem via *level resolutions* of maximal planar graphs.
+
+See `paper.tex` for full definitions, conjectures, and findings.
+
+## Files
+
+### Core library
+- `level_cycles.py` — levels, level subgraphs, level cycles, resolution
+  enumeration (used by old-definition coverage).
+- `triangulation_gen.py` — vertex-insertion + flip closure (good to n=10).
+- `triangulation_gen_fast.py` — WL-hash pre-filter for n ≥ 11.
+- `balanced_layout.py` — Tutte-init random-search planar layout.
+- `four_color.py` — level 4-coloring via parity 2-coloring of L_k.
+
+### Experiments
+- `coverage_new_def.py` — **coverage under the cleaner definition**:
+  G' is a level resolution of G via S iff its parity subgraphs are
+  bipartite. Reachability reduces to "G' admits a bipartite 2-partition
+  with cardinality matching some BFS-realizable parity split."
+- `coverage.py`, `coverage_fast.py`, `coverage_chunked.py` — coverage
+  under the OLD (stricter) definition involving specific edge flips on
+  level cycles.
+- `face_counting.py` — per-target preimage counts (N_iso, N_paths) under
+  the old definition.
+- `orbit_check.py` — orbit-counting with k-flip reverse-preimages (used
+  for old-definition icosahedron analysis).
+
+### Visualizations
+- `plot_oct.py`, `n7_examples.py`, `four_color_viz.py`.
+
+## Summary under the new definition
+
+| n  | iso-classes | reachable | md4 reachable |
+|----|-------------|-----------|---------------|
+| 6  | 2           | 2         | 1/1           |
+| 7  | 5           | 5         | 1/1           |
+| 8  | 14          | 14        | 2/2           |
+| 9  | 50          | 50        | 5/5           |
+| 10 | 233         | 233       | 12/12         |
+| 11 | 1249        | 1249      | 34/34         |
+| 12 | icosahedron | reachable | yes           |
+
+**Every iso-class is reachable** at every tested size. The previously
+"uncovered" classes T1 (n=7) and T6 (n=8) under the old definition are
+both reachable under the cleaner definition.
+
+The new definition makes coverage equivalent to 4CT plus a BFS-realizable
+partition cardinality constraint, raising the question of what additional
+structure on the preimage G would make the framework non-circular.
+
+## Dependencies
+
+```
+pip install networkx matplotlib numpy scipy
+```

papers/level_resolutions_of_maximal_planar_graphs/experiments/balanced_layout.py

@@ -0,0 +1,212 @@
+"""
+Face-area-balanced planar layout for maximal planar graphs.
+
+Starts from a Tutte embedding (outer face on a circle) and uses random-search
+optimization to equalize interior face areas while maintaining planarity.
+"""
+import networkx as nx
+import numpy as np
+import scipy.linalg
+
+
+def _get_all_faces(emb):
+    """Enumerate all faces of a PlanarEmbedding."""
+    seen, faces = set(), []
+    for v in emb.nodes():
+        for w in emb[v]:
+            if (v, w) not in seen:
+                face = emb.traverse_face(v, w)
+                for i in range(len(face)):
+                    seen.add((face[i], face[(i+1) % len(face)]))
+                faces.append(tuple(face))
+    return faces
+
+
+def _tutte_layout(G, outer_face):
+    """Standard uniform Tutte embedding with outer face on unit circle."""
+    outer = list(outer_face)
+    k = len(outer)
+    fixed_pos = {v: np.array([np.cos(2*np.pi*i/k - np.pi/2),
+                              np.sin(2*np.pi*i/k - np.pi/2)])
+                 for i, v in enumerate(outer)}
+    interior = [v for v in G.nodes() if v not in fixed_pos]
+    if not interior:
+        return fixed_pos
+    int_idx = {v: i for i, v in enumerate(interior)}
+    m = len(interior)
+    A = np.zeros((m, m)); bx = np.zeros(m); by = np.zeros(m)
+    for v in interior:
+        i = int_idx[v]; nbrs = list(G.neighbors(v)); deg = len(nbrs)
+        A[i, i] = -1.0
+        for u in nbrs:
+            if u in int_idx:
+                A[i, int_idx[u]] += 1.0 / deg
+            else:
+                bx[i] -= fixed_pos[u][0] / deg
+                by[i] -= fixed_pos[u][1] / deg
+    px = scipy.linalg.solve(A, bx)
+    py = scipy.linalg.solve(A, by)
+    pos = dict(fixed_pos)
+    for v in interior:
+        pos[v] = np.array([px[int_idx[v]], py[int_idx[v]]])
+    return pos
+
+
+def _segments_cross(p1, p2, p3, p4):
+    """True iff open segments p1-p2 and p3-p4 cross (not at endpoints)."""
+    def cross(o, a, b):
+        return (a[0]-o[0])*(b[1]-o[1]) - (a[1]-o[1])*(b[0]-o[0])
+    d1 = cross(p3, p4, p1); d2 = cross(p3, p4, p2)
+    d3 = cross(p1, p2, p3); d4 = cross(p1, p2, p4)
+    return ((d1 > 0 and d2 < 0) or (d1 < 0 and d2 > 0)) and \
+           ((d3 > 0 and d4 < 0) or (d3 < 0 and d4 > 0))
+
+
+def _count_crossings(G, pos):
+    """Number of pairs of edges whose straight-line segments cross."""
+    edges = list(G.edges())
+    n = 0
+    for i in range(len(edges)):
+        a, b = edges[i]
+        for j in range(i+1, len(edges)):
+            c, d = edges[j]
+            if len({a, b, c, d}) < 4:
+                continue
+            if _segments_cross(pos[a], pos[b], pos[c], pos[d]):
+                n += 1
+    return n
+
+
+def _face_area(pos, face):
+    """Unsigned area of a triangular (or polygonal) face."""
+    pts = [pos[v] for v in face]
+    return 0.5 * abs((pts[1][0]-pts[0][0])*(pts[2][1]-pts[0][1])
+                   - (pts[1][1]-pts[0][1])*(pts[2][0]-pts[0][0]))
+
+
+def _edge_lengths(G, pos):
+    return [np.linalg.norm(pos[u] - pos[v]) for u, v in G.edges()]
+
+
+def _score(G, pos, interior_faces, w_area=1.0, w_edge=0.3):
+    """Combined penalty: normalized variance of face areas + edge lengths."""
+    areas = np.array([_face_area(pos, f) for f in interior_faces])
+    a_mean = areas.mean()
+    area_score = np.sum((areas - a_mean)**2) / (a_mean**2)
+    lens = np.array(_edge_lengths(G, pos))
+    l_mean = lens.mean()
+    edge_score = np.sum((lens - l_mean)**2) / (l_mean**2)
+    return w_area * area_score + w_edge * edge_score
+
+
+def balanced_planar_layout(
+    G,
+    outer_face,
+    n_explore=8000,
+    n_refine=4000,
+    explore_step=0.30,
+    refine_step=0.05,
+    w_area=1.0,
+    w_edge=0.3,
+    seed=42,
+    verbose=False,
+):
+    """
+    Compute a planar layout for a maximal planar graph G whose interior faces
+    have roughly equal area and whose edges have roughly equal length.
+
+    Starts from a uniform Tutte embedding (outer face on a unit circle in CCW
+    order starting at the south pole), then runs random-search optimization,
+    accepting only moves that keep the layout planar.
+
+    Parameters
+    ----------
+    G : networkx.Graph
+        A simple connected planar graph; intended for maximal planar graphs.
+    outer_face : tuple
+        Cyclic ordering of vertices on the outer face. Pinned during
+        optimization.
+    n_explore, n_refine : int
+        Number of iterations in the exploration and refinement phases.
+    explore_step, refine_step : float
+        Standard deviation of Gaussian moves in each phase.
+    w_area, w_edge : float
+        Weights for the face-area and edge-length penalties in the score.
+    seed : int
+        RNG seed for reproducibility.
+    verbose : bool
+        If True, print progress.
+
+    Returns
+    -------
+    dict[node, np.ndarray]
+        Position of each vertex.
+    """
+    is_planar, emb = nx.check_planarity(G)
+    if not is_planar:
+        raise ValueError("G is not planar")
+
+    pos = _tutte_layout(G, outer_face)
+
+    all_faces = _get_all_faces(emb)
+    outer_set = set(outer_face)
+    interior_faces = [f for f in all_faces if set(f) != outer_set]
+
+    if _count_crossings(G, pos) != 0:
+        raise RuntimeError("Initial Tutte layout is not planar; check inputs")
+
+    interior_nodes = [v for v in G.nodes() if v not in outer_set]
+    if not interior_nodes:
+        return pos
+
+    rng = np.random.default_rng(seed)
+    best_score = _score(G, pos, interior_faces, w_area, w_edge)
+    best_pos = {k: v.copy() for k, v in pos.items()}
+    initial_score = best_score
+
+    for phase_name, n_iter, step in [("explore", n_explore, explore_step),
+                                     ("refine",  n_refine,  refine_step)]:
+        pos = {k: v.copy() for k, v in best_pos.items()}
+        accepted = 0
+        for it in range(n_iter):
+            v = interior_nodes[rng.integers(len(interior_nodes))]
+            delta = rng.normal(0, step, 2)
+            old = pos[v].copy()
+            pos[v] = old + delta
+            if _count_crossings(G, pos) == 0:
+                new_score = _score(G, pos, interior_faces, w_area, w_edge)
+                if new_score < best_score:
+                    best_score = new_score
+                    best_pos = {k: vv.copy() for k, vv in pos.items()}
+                    accepted += 1
+                else:
+                    pos[v] = old
+            else:
+                pos[v] = old
+        if verbose:
+            print(f"  {phase_name}: score {best_score:.4f}, "
+                  f"accepted {accepted}/{n_iter}")
+
+    if verbose:
+        print(f"Initial score: {initial_score:.4f}")
+        print(f"Final score:   {best_score:.4f}")
+
+    return best_pos
+
+
+# ── Demo ─────────────────────────────────────────────────────────────────────
+
+if __name__ == "__main__":
+    edges = [(0,1),(0,2),(0,3),(0,4),(0,5),(0,6),
+             (1,2),(1,3),(1,4),(1,5),(1,6),
+             (2,3),(2,4),(3,5),(4,6)]
+    G = nx.Graph()
+    G.add_nodes_from(range(7))
+    G.add_edges_from(edges)
+
+    pos = balanced_planar_layout(G, outer_face=(0, 3, 5), verbose=True)
+    print("\nFinal positions:")
+    for v in sorted(pos):
+        p = pos[v]
+        print(f"  {v}: ({p[0]:+.3f}, {p[1]:+.3f})")
+    print(f"Crossings: {_count_crossings(G, pos)}")

papers/level_resolutions_of_maximal_planar_graphs/experiments/coverage.py

@@ -0,0 +1,85 @@
+"""
+Coverage analysis: for each pair (n, target-class restriction), check whether
+every iso-class in the restriction is reachable as a level resolution of
+some triangulation on n vertices.
+"""
+import networkx as nx
+import time
+from level_cycles import all_level_resolutions
+from triangulation_gen import enumerate_all_triangulations
+
+
+def iso_class(G, reps):
+    for i, r in enumerate(reps):
+        if nx.is_isomorphic(G, r):
+            return i
+    return -1
+
+
+def resolution_classes(G, reps):
+    return {iso_class(Gp, reps) for Gp, _, _, _ in all_level_resolutions(G)}
+
+
+def min_degree(G):
+    return min(d for _, d in G.degree())
+
+
+def coverage_report(n, target_filter=None, source_filter=None):
+    """
+    target_filter / source_filter: callables G -> bool, or None for "any".
+    """
+    print(f"\n{'='*60}\nn = {n}\n{'='*60}")
+    t0 = time.time()
+    reps = enumerate_all_triangulations(n)
+    print(f"Total iso-classes: {len(reps)}")
+
+    if source_filter is None:
+        sources = list(range(len(reps)))
+        print(f"Sources: all {len(sources)}")
+    else:
+        sources = [i for i, G in enumerate(reps) if source_filter(G)]
+        print(f"Sources (filtered): {len(sources)}")
+
+    if target_filter is None:
+        targets = set(range(len(reps)))
+    else:
+        targets = {i for i, G in enumerate(reps) if target_filter(G)}
+    print(f"Targets: {len(targets)} iso-classes")
+    for i in sorted(targets):
+        deg = sorted((d for _, d in reps[i].degree()), reverse=True)
+        print(f"  T{i}: degree {deg}")
+
+    reached, sources_per_target = set(), {i: [] for i in targets}
+    for src_i in sources:
+        prod = resolution_classes(reps[src_i], reps) & targets
+        for p in prod:
+            sources_per_target[p].append(src_i)
+        reached |= prod
+
+    print("\nCoverage:")
+    for i in sorted(targets):
+        sources_list = sources_per_target[i]
+        status = "REACHED" if sources_list else "UNREACHABLE"
+        s = ", ".join(f"T{j}" for j in sources_list[:6])
+        if len(sources_list) > 6:
+            s += f", ... ({len(sources_list)} total)"
+        elif not sources_list:
+            s = "(none)"
+        print(f"  T{i}: {status} via {s}")
+
+    uncov = targets - reached
+    print(f"\nUncovered: {sorted(uncov)}")
+    print(f"Time: {time.time()-t0:.1f}s")
+
+
+if __name__ == "__main__":
+    # General coverage (any source, any target)
+    for n in [6, 7]:
+        coverage_report(n)
+
+    # md3 sources -> md4 targets
+    print("\n\n" + "#"*60)
+    print("md3 sources -> md4 targets")
+    print("#"*60)
+    for n in [6, 7, 8]:
+        coverage_report(n, target_filter=lambda G: min_degree(G) >= 4)

papers/level_resolutions_of_maximal_planar_graphs/experiments/coverage_chunked.py

@@ -0,0 +1,99 @@
+"""Chunked coverage analysis. Processes a slice of sources per invocation.
+Saves and resumes progress via a pickle file. Run repeatedly until done."""
+import sys; sys.path.insert(0, '/home/claude/build')
+import networkx as nx
+import pickle
+import time
+import os
+from collections import defaultdict
+from level_cycles import all_level_resolutions
+from triangulation_gen import enumerate_all_triangulations
+
+STATE_FILE = '/tmp/n11_state.pkl'
+N = 11
+CHUNK_SECONDS = 200  # how long this invocation should work
+
+def degree_seq(G):
+    return tuple(sorted((d for _, d in G.degree()), reverse=True))
+
+def min_degree(G):
+    return min(d for _, d in G.degree())
+
+# Load or initialize state
+if os.path.exists(STATE_FILE):
+    with open(STATE_FILE, 'rb') as f:
+        state = pickle.load(f)
+    print(f"Resumed. Done: {state['next_idx']}/{state['total']}")
+else:
+    t0 = time.time()
+    print("Initializing...")
+    reps = enumerate_all_triangulations(N)
+    print(f"  {len(reps)} iso-classes in {time.time()-t0:.1f}s")
+    deg_b = defaultdict(list); hash_b = defaultdict(list)
+    for i, G in enumerate(reps):
+        ds = degree_seq(G)
+        deg_b[ds].append(i)
+        hash_b[(ds, nx.weisfeiler_lehman_graph_hash(G))].append((i, G))
+    state = {
+        'reps': reps,
+        'deg_b': dict(deg_b),
+        'hash_b': dict(hash_b),
+        'produced_by': {},
+        'next_idx': 0,
+        'total': len(reps),
+    }
+    with open(STATE_FILE, 'wb') as f:
+        pickle.dump(state, f)
+    print(f"Init done in {time.time()-t0:.1f}s")
+
+reps = state['reps']
+deg_b = state['deg_b']
+hash_b = state['hash_b']
+produced_by = state['produced_by']
+
+def fast_iso(G):
+    ds = degree_seq(G)
+    if ds not in deg_b: return -1
+    h = nx.weisfeiler_lehman_graph_hash(G)
+    for idx, rep in hash_b.get((ds, h), []):
+        if nx.is_isomorphic(G, rep): return idx
+    return -1
+
+start = time.time()
+i = state['next_idx']
+processed = 0
+while i < state['total'] and time.time() - start < CHUNK_SECONDS:
+    G = reps[i]
+    s = set()
+    for Gp, _, _, _ in all_level_resolutions(G):
+        c = fast_iso(Gp)
+        if c >= 0: s.add(c)
+    produced_by[i] = s
+    i += 1
+    processed += 1
+
+state['next_idx'] = i
+state['produced_by'] = produced_by
+with open(STATE_FILE, 'wb') as f:
+    pickle.dump(state, f)
+
+elapsed = time.time() - start
+print(f"Processed {processed} sources in {elapsed:.0f}s. "
+      f"Now at {i}/{state['total']}.")
+
+if i >= state['total']:
+    print("\n=== ALL DONE ===")
+    all_prod = set().union(*produced_by.values())
+    uncov = set(range(state['total'])) - all_prod
+    md4_idx = set(j for j, G in enumerate(reps) if min_degree(G) >= 4)
+    md4_reached = set().union(*(produced_by[k] & md4_idx
+                                for k in range(state['total'])))
+    md4_uncov = md4_idx - md4_reached
+    print(f"General: {state['total'] - len(uncov)}/{state['total']} reached")
+    if uncov:
+        for j in sorted(uncov)[:20]:
+            print(f"  uncovered T{j}: degree {degree_seq(reps[j])}")
+    print(f"md4: {len(md4_idx) - len(md4_uncov)}/{len(md4_idx)} reached")
+    if md4_uncov:
+        for j in sorted(md4_uncov):
+            print(f"  uncovered md4 T{j}: degree {degree_seq(reps[j])}")

papers/level_resolutions_of_maximal_planar_graphs/experiments/coverage_fast.py

@@ -0,0 +1,78 @@
+"""
+Fast coverage analysis using WL hash for iso-class lookup.
+"""
+import sys; sys.path.insert(0, '/home/claude/build')
+import networkx as nx
+import time
+from collections import defaultdict
+from level_cycles import all_level_resolutions
+from triangulation_gen_fast import enumerate_all_triangulations_fast, wl_hash
+
+
+def build_iso_index(reps):
+    """Map WL hash -> list of (idx, G) so iso-checks only run within hash bucket."""
+    idx = defaultdict(list)
+    for i, G in enumerate(reps):
+        idx[wl_hash(G)].append((i, G))
+    return idx
+
+
+def iso_class_fast(G, iso_index):
+    h = wl_hash(G)
+    for i, R in iso_index.get(h, []):
+        if nx.is_isomorphic(G, R):
+            return i
+    return -1
+
+
+def resolution_classes_fast(G, iso_index):
+    return {iso_class_fast(Gp, iso_index)
+            for Gp, _, _, _ in all_level_resolutions(G)
+            if iso_class_fast(Gp, iso_index) >= 0}
+
+
+def min_degree(G):
+    return min(d for _, d in G.degree())
+
+
+def coverage_at(n, report_every=100):
+    t0 = time.time()
+    reps = enumerate_all_triangulations_fast(n)
+    iso_index = build_iso_index(reps)
+    print(f"n={n}: {len(reps)} iso-classes enumerated and indexed "
+          f"({time.time()-t0:.1f}s)")
+    md4_idx = set(i for i, G in enumerate(reps) if min_degree(G) >= 4)
+
+    t1 = time.time()
+    produced_by = {}
+    for i, G in enumerate(reps):
+        produced_by[i] = resolution_classes_fast(G, iso_index)
+        if (i + 1) % report_every == 0:
+            elapsed = time.time() - t1
+            eta = elapsed * (len(reps) - i - 1) / (i + 1)
+            print(f"  {i+1}/{len(reps)}  ({elapsed:.0f}s, ETA {eta:.0f}s)")
+    print(f"Resolution computation: {time.time()-t1:.1f}s")
+
+    all_prod = set().union(*produced_by.values())
+    uncovered = set(range(len(reps))) - all_prod
+    md4_reached = set().union(*(produced_by[i] & md4_idx
+                                for i in range(len(reps))))
+    md4_uncov = md4_idx - md4_reached
+
+    print(f"\nGeneral coverage: {len(reps) - len(uncovered)}/{len(reps)}")
+    if uncovered:
+        for i in sorted(uncovered):
+            deg = sorted((d for _, d in reps[i].degree()), reverse=True)
+            print(f"  uncovered T{i}: degree {deg}")
+    print(f"md4 coverage:     {len(md4_idx) - len(md4_uncov)}/{len(md4_idx)}")
+    if md4_uncov:
+        for i in sorted(md4_uncov):
+            deg = sorted((d for _, d in reps[i].degree()), reverse=True)
+            print(f"  uncovered md4 T{i}: degree {deg}")
+    print(f"Total: {time.time()-t0:.1f}s")
+
+
+if __name__ == "__main__":
+    for n in [11]:
+        print("=" * 60)
+        coverage_at(n, report_every=100)

papers/level_resolutions_of_maximal_planar_graphs/experiments/coverage_new_def.py

@@ -0,0 +1,104 @@
+"""
+Coverage analysis under the NEW definition of level resolution:
+G' is a level resolution of G via S iff both parity subgraphs of G' (using
+G's BFS-from-S parities) are bipartite.
+
+Since the iso-class of G' doesn't depend on the specific labeling (only the
+cardinalities |V_e|, |V_o| matter via permutation), G' is reachable iff:
+- there exist achievable parity cardinalities (s, n-s) from some (G, S),
+- G' admits a 2-partition into bipartite-induced subgraphs of those sizes.
+"""
+import sys; sys.path.insert(0, '/home/claude/build')
+import networkx as nx
+import time
+from itertools import combinations
+from collections import defaultdict
+from level_cycles import (
+    compute_levels, get_all_faces, level_sources,
+)
+from triangulation_gen import enumerate_all_triangulations
+from triangulation_gen_fast import enumerate_all_triangulations_fast
+
+
+def min_degree(G):
+    return min(d for _, d in G.degree())
+
+
+def achievable_parity_splits(reps):
+    """Return set of (|V_e|, |V_o|) cardinality tuples achievable across all
+    (G, S) pairs from the given triangulation reps."""
+    splits = set()
+    for G in reps:
+        ip, emb = nx.check_planarity(G)
+        if not ip: continue
+        for kind, label, source_set in level_sources(G, emb):
+            levels = compute_levels(G, source_set)
+            even = sum(1 for v in levels.values() if v % 2 == 0)
+            odd  = sum(1 for v in levels.values() if v % 2 == 1)
+            splits.add((even, odd))
+    return splits
+
+
+def has_bipartite_partition(G, sizes):
+    """Does G admit a 2-partition (S, V\\S) with |S| in `sizes` (or n-|S| in
+    sizes) such that both G[S] and G[V\\S] are bipartite?"""
+    V = list(G.nodes())
+    n = len(V)
+    accept_sizes = set(s for s in sizes if 1 <= s <= n - 1)
+    accept_sizes |= set(n - s for s in sizes if 1 <= n - s <= n - 1)
+    seen_sizes = sorted(accept_sizes)
+    # Avoid duplicate work: check size s and size n-s give same partition
+    seen = set()
+    for size in seen_sizes:
+        if size > n - size:
+            continue
+        for S_tuple in combinations(V, size):
+            S = set(S_tuple)
+            key = frozenset(S)
+            if key in seen: continue
+            seen.add(key)
+            if nx.is_bipartite(G.subgraph(S)) and \
+               nx.is_bipartite(G.subgraph(set(V) - S)):
+                return True
+    return False
+
+
+def analyze(n, use_fast=False, verbose=True):
+    print(f"\n{'='*60}\nn = {n}\n{'='*60}")
+    t0 = time.time()
+    reps = (enumerate_all_triangulations_fast(n) if use_fast
+            else enumerate_all_triangulations(n))
+    md4_idx = set(i for i, G in enumerate(reps) if min_degree(G) >= 4)
+    print(f"{len(reps)} iso-classes, {len(md4_idx)} md4")
+
+    splits = achievable_parity_splits(reps)
+    achievable_sizes = set(s for s, _ in splits)
+    print(f"Achievable parity splits: {sorted(splits)}")
+    print(f"Achievable sizes: {sorted(achievable_sizes)}")
+
+    print(f"\nReachability under NEW definition:")
+    reachable, unreachable_md4 = [], []
+    for i, G in enumerate(reps):
+        if has_bipartite_partition(G, achievable_sizes):
+            reachable.append(i)
+        elif i in md4_idx:
+            unreachable_md4.append(i)
+        if verbose and i not in reachable and i in md4_idx:
+            deg = sorted((d for _, d in G.degree()), reverse=True)
+            print(f"  T{i} (md4, deg {deg}): UNREACHABLE")
+    not_reachable = set(range(len(reps))) - set(reachable)
+    print(f"Reachable: {len(reachable)}/{len(reps)}")
+    print(f"md4 reachable: {len(md4_idx & set(reachable))}/{len(md4_idx)}")
+    if not_reachable:
+        print(f"Unreachable iso-classes: {sorted(not_reachable)}")
+        for i in sorted(not_reachable):
+            deg = sorted((d for _, d in reps[i].degree()), reverse=True)
+            marker = " (md4)" if i in md4_idx else ""
+            print(f"  T{i}: degree {deg}{marker}")
+    print(f"Time: {time.time()-t0:.1f}s")
+    return reachable, unreachable_md4
+
+
+if __name__ == "__main__":
+    for n in [6, 7, 8, 9, 10, 11]:
+        analyze(n, use_fast=(n >= 11))

papers/level_resolutions_of_maximal_planar_graphs/experiments/face_counting.py

@@ -0,0 +1,125 @@
+"""
+Double counting at the face/source level.
+For each (G_iso, source-up-to-aut, flip choice), apply the resolution and
+track which G'_iso it lands on. Then aggregate.
+
+For each md4 target G' we report:
+- N_iso(G')   = # distinct preimage iso-classes G that resolve to G'
+- N_paths(G') = # distinct (G_iso, source orbit, flip choice) triples
+
+For each source G we report:
+- f(G) = # (source, flip choice) pairs producing any md4 target
+"""
+import sys; sys.path.insert(0, '/home/claude/build')
+import networkx as nx
+from collections import defaultdict
+from itertools import product as iproduct
+
+from triangulation_gen import enumerate_all_triangulations
+from triangulation_gen_fast import enumerate_all_triangulations_fast
+from level_cycles import (
+    compute_levels, get_all_faces, get_level_cycles_by_parity,
+    get_flip_candidates, level_sources,
+)
+
+
+def min_degree(G):
+    return min(d for _, d in G.degree())
+
+
+def iso_class(G, reps):
+    for i, r in enumerate(reps):
+        if nx.is_isomorphic(G, r):
+            return i
+    return -1
+
+
+def all_resolutions(G):
+    """Yield (G_resolved, source_kind, source_label) for every resolution.
+    Equivalent paths are NOT deduplicated here — caller can do so via iso."""
+    ip, emb = nx.check_planarity(G)
+    if not ip: return
+    for kind, label, source_set in level_sources(G, emb):
+        levels = compute_levels(G, source_set)
+        odd_cycles, even_cycles = get_level_cycles_by_parity(G, emb, levels)
+        cand_lists = []
+        ok = True
+        for cycle in odd_cycles:
+            cands = get_flip_candidates(G, emb, cycle)
+            if not cands: ok = False; break
+            cand_lists.append(cands)
+        if not ok: continue
+        for cycle in even_cycles:
+            cands = get_flip_candidates(G, emb, cycle)
+            cand_lists.append(cands + [None])
+        if not cand_lists:
+            yield (G.copy(), kind, label); continue
+        for choices in iproduct(*cand_lists):
+            Gp = G.copy()
+            applied = set(); ok2 = True
+            for c in choices:
+                if c is None: continue
+                u, v, w, x = c
+                if frozenset([u,v]) in applied: ok2 = False; break
+                if not Gp.has_edge(u, v) or Gp.has_edge(w, x): ok2 = False; break
+                Gp.remove_edge(u, v); Gp.add_edge(w, x)
+                applied.add(frozenset([u, v]))
+            if not ok2: continue
+            ipp, _ = nx.check_planarity(Gp)
+            if not ipp: continue
+            yield (Gp, kind, label)
+
+
+def analyze(n, use_fast=False):
+    print(f"\n{'='*70}\nn = {n}\n{'='*70}")
+    if use_fast:
+        reps = enumerate_all_triangulations_fast(n)
+    else:
+        reps = enumerate_all_triangulations(n)
+    md4_idx = [i for i, G in enumerate(reps) if min_degree(G) >= 4]
+    md4_set = set(md4_idx)
+
+    # For each source iso-class, compute # paths to each target iso-class
+    pre_count = defaultdict(lambda: defaultdict(int))   # target -> source -> count
+    src_paths_md4 = defaultdict(int)  # source -> total paths landing on md4
+
+    print(f"Computing resolutions for {len(reps)} sources...", flush=True)
+    for src_i, G in enumerate(reps):
+        # Compute all distinct (target iso, source orbit) pairs from this G
+        for Gp, kind, label in all_resolutions(G):
+            tgt = iso_class(Gp, reps)
+            if tgt < 0: continue
+            pre_count[tgt][src_i] += 1
+            if tgt in md4_set:
+                src_paths_md4[src_i] += 1
+
+    # Report by md4 target
+    print(f"\n{'target':28s} {'N_iso':>6} {'N_paths':>9} {'min path/iso':>14}")
+    print("-" * 70)
+    for tgt in md4_idx:
+        deg = sorted((d for _, d in reps[tgt].degree()), reverse=True)
+        sources_hitting = list(pre_count[tgt].keys())
+        N_iso = len(sources_hitting)
+        N_paths = sum(pre_count[tgt].values())
+        min_paths_per_src = (min(pre_count[tgt].values())
+                              if sources_hitting else 0)
+        print(f"{str(deg):28s} {N_iso:>6} {N_paths:>9} {min_paths_per_src:>14}")
+    # Double-counting check:
+    total_md4_paths = sum(sum(d.values()) for d in pre_count.values()
+                           if any(k in md4_set for k in [...]))
+    # Simpler total
+    total_paths = sum(src_paths_md4.values())
+    avg_per_target = total_paths / len(md4_idx) if md4_idx else 0
+    print("-" * 70)
+    print(f"Total md4 paths: {total_paths}")
+    print(f"Average per md4 target: {avg_per_target:.1f}")
+    if md4_idx:
+        min_N_iso = min(len(pre_count[t]) for t in md4_idx)
+        print(f"Minimum N_iso over md4 targets: {min_N_iso}")
+
+    return pre_count, md4_idx
+
+
+if __name__ == "__main__":
+    for n in [6, 7, 8, 9, 10]:
+        analyze(n, use_fast=(n >= 10))

papers/level_resolutions_of_maximal_planar_graphs/experiments/four_color.py

@@ -0,0 +1,98 @@
+"""
+4-color G' using the level structure from G:
+  - Even-level subgraph: 2-color with RED/BLUE via BFS
+  - Odd-level subgraph:  2-color with YELLOW/GREEN via BFS
+
+This succeeds iff each parity subgraph is bipartite — which is the goal of
+the level resolution. If a parity subgraph contains an odd cycle, BFS will
+find a conflict and we report which edge violates the 2-coloring.
+"""
+import networkx as nx
+import numpy as np
+from collections import deque
+
+
+def two_color_subgraph(G_sub, color_a, color_b):
+    """
+    Two-color the (possibly disconnected) subgraph via BFS, alternating
+    color_a / color_b by BFS distance parity from a root in each component.
+
+    Returns:
+        coloring: dict[node, color]
+        bad_edges: list of edges where adjacent vertices got the same colour
+                   (empty iff subgraph is bipartite)
+    """
+    coloring = {}
+    bad_edges = []
+    for start in G_sub.nodes():
+        if start in coloring:
+            continue
+        coloring[start] = color_a
+        queue = deque([start])
+        while queue:
+            v = queue.popleft()
+            for w in G_sub.neighbors(v):
+                if w not in coloring:
+                    coloring[w] = color_b if coloring[v] == color_a else color_a
+                    queue.append(w)
+                elif coloring[w] == coloring[v]:
+                    e = tuple(sorted([v, w]))
+                    if e not in bad_edges:
+                        bad_edges.append(e)
+    return coloring, bad_edges
+
+
+def four_color_via_levels(G_prime, levels):
+    """
+    4-color G' using level labels from G.
+    Even-level vertices get RED/BLUE; odd-level get YELLOW/GREEN.
+
+    Returns:
+        coloring: dict[node, str]
+        bad_edges: dict with keys 'even', 'odd', 'cross' for violations
+                   within each parity subgraph and between them
+                   (the cross list should always be empty by construction)
+    """
+    even_nodes = [v for v in G_prime.nodes() if levels[v] % 2 == 0]
+    odd_nodes  = [v for v in G_prime.nodes() if levels[v] % 2 == 1]
+
+    even_sub = G_prime.subgraph(even_nodes).copy()
+    odd_sub  = G_prime.subgraph(odd_nodes).copy()
+
+    coloring_even, bad_even = two_color_subgraph(even_sub, 'red',    'blue')
+    coloring_odd,  bad_odd  = two_color_subgraph(odd_sub,  'yellow', 'green')
+
+    coloring = {**coloring_even, **coloring_odd}
+
+    cross_bad = []
+    for u, v in G_prime.edges():
+        if coloring[u] == coloring[v]:
+            cross_bad.append(tuple(sorted([u, v])))
+
+    return coloring, {
+        'even':  bad_even,
+        'odd':   bad_odd,
+        'cross': [e for e in cross_bad
+                  if tuple(sorted(e)) not in bad_even
+                  and tuple(sorted(e)) not in bad_odd],
+    }
+
+
+if __name__ == "__main__":
+    # Quick demo on the n=7 graph used earlier
+    edges = [(0,1),(0,2),(0,3),(0,4),(0,5),(0,6),
+             (1,2),(1,3),(1,4),(1,5),(1,6),
+             (2,3),(2,4),(3,5),(4,6)]
+    G = nx.Graph(); G.add_nodes_from(range(7)); G.add_edges_from(edges)
+    Gp = G.copy(); Gp.remove_edge(1, 2); Gp.add_edge(4, 3)
+
+    levels = {0:0, 3:0, 5:0, 1:1, 2:1, 4:1, 6:1}
+    coloring, bad = four_color_via_levels(Gp, levels)
+    print("Coloring:")
+    for v in sorted(coloring):
+        print(f"  vertex {v} (level {levels[v]}): {coloring[v]}")
+    print(f"\nBad edges in even subgraph:  {bad['even']}")
+    print(f"Bad edges in odd subgraph:   {bad['odd']}")
+    print(f"Cross-parity bad edges:      {bad['cross']}")
+    print(f"\nValid 4-coloring: "
+          f"{not any(bad.values())}")

papers/level_resolutions_of_maximal_planar_graphs/experiments/four_color_viz.py

@@ -0,0 +1,256 @@
+"""
+For each n=7 example (G → G' covering all 5 iso-classes), apply the level
+4-coloring to G' using levels from G, and visualize.
+"""
+import sys; sys.path.insert(0, '/home/claude')
+import networkx as nx
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+import numpy as np
+from collections import deque
+from itertools import combinations, product as iproduct
+
+from balanced_layout import balanced_planar_layout, _get_all_faces
+from four_color import four_color_via_levels
+
+def compute_levels(G, outer_face):
+    levels, queue = {}, deque()
+    for v in outer_face:
+        levels[v] = 0; queue.append(v)
+    while queue:
+        v = queue.popleft()
+        for w in G.neighbors(v):
+            if w not in levels:
+                levels[w] = levels[v] + 1; queue.append(w)
+    return levels
+
+def get_odd_level_cycles(emb, levels):
+    odd = []
+    for face in _get_all_faces(emb):
+        lv = {levels.get(v) for v in face}
+        if len(lv) == 1 and None not in lv and len(face) % 2 == 1:
+            odd.append(face)
+    return odd
+
+def get_flip_candidates(G, emb, cycle):
+    cands = []
+    n = len(cycle)
+    for i in range(n):
+        u, v = cycle[i], cycle[(i+1) % n]
+        if not G.has_edge(u, v): continue
+        f1 = emb.traverse_face(u, v); f2 = emb.traverse_face(v, u)
+        if len(f1) != 3 or len(f2) != 3: continue
+        w = next(x for x in f1 if x != u and x != v)
+        x = next(y for y in f2 if y != u and y != v)
+        if w == x or G.has_edge(w, x): continue
+        cands.append((u, v, w, x))
+    return cands
+
+def is_triangulation(G):
+    n = G.number_of_nodes()
+    if G.number_of_edges() != 3*n - 6: return False
+    if not nx.is_connected(G): return False
+    ip, emb = nx.check_planarity(G)
+    if not ip: return False
+    seen = set()
+    for v in emb.nodes():
+        for w in emb[v]:
+            if (v, w) not in seen:
+                f = emb.traverse_face(v, w)
+                for i in range(len(f)):
+                    seen.add((f[i], f[(i+1) % len(f)]))
+                if len(f) != 3: return False
+    return True
+
+# ── Reps ─────────────────────────────────────────────────────────────────────
+nodes = list(range(7))
+all_edges = list(combinations(nodes, 2))
+reps = []
+for sub in combinations(all_edges, 15):
+    G = nx.Graph(); G.add_nodes_from(nodes); G.add_edges_from(sub)
+    if is_triangulation(G):
+        if all(not nx.is_isomorphic(G, r) for r in reps):
+            reps.append(G)
+        if len(reps) == 5: break
+
+def iso_class(G):
+    for i, r in enumerate(reps):
+        if nx.is_isomorphic(G, r): return i
+    return -1
+
+# Same pairs as before
+desired_pairs = [(0, 2), (1, 4), (2, 0), (3, 1), (4, 3)]
+
+def find_example(src_idx, tgt_idx):
+    G = reps[src_idx]
+    ip, emb = nx.check_planarity(G)
+    for outer_face in _get_all_faces(emb):
+        levels = compute_levels(G, outer_face)
+        odd_cycles = get_odd_level_cycles(emb, levels)
+        if not odd_cycles: continue
+        candidate_lists = []
+        ok = True
+        for cycle in odd_cycles:
+            cands = get_flip_candidates(G, emb, cycle)
+            if not cands: ok = False; break
+            candidate_lists.append(cands)
+        if not ok: continue
+        for choices in iproduct(*candidate_lists):
+            Gp = G.copy()
+            valid = True
+            for u, v, w, x in choices:
+                if not Gp.has_edge(u, v) or Gp.has_edge(w, x):
+                    valid = False; break
+                Gp.remove_edge(u, v); Gp.add_edge(w, x)
+            if valid and is_triangulation(Gp) and iso_class(Gp) == tgt_idx:
+                return G, outer_face, choices, Gp
+    return None
+
+examples = [find_example(s, t) for s, t in desired_pairs]
+print("Examples found:", sum(1 for e in examples if e is not None))
+
+# ── 4-color and visualize ────────────────────────────────────────────────────
+
+color_hex = {'red':'#e53935','blue':'#1e88e5',
+             'yellow':'#fdd835','green':'#43a047'}
+
+fig, axes = plt.subplots(len(examples), 2, figsize=(14, 5.5 * len(examples)))
+if len(examples) == 1:
+    axes = np.array([axes])
+
+print("\nColoring each G':")
+for row, ex in enumerate(examples):
+    if ex is None: continue
+    G, outer_face, choices, Gp = ex
+    src, tgt = desired_pairs[row]
+
+    levels_G = compute_levels(G, outer_face)
+    coloring, bad = four_color_via_levels(Gp, levels_G)
+    valid = not any(bad.values())
+    print(f"  T{src}→T{tgt}: valid 4-coloring = {valid}")
+    if not valid:
+        print(f"    bad even: {bad['even']}, bad odd: {bad['odd']}")
+
+    # Layout for G and G'
+    _, emb_G = nx.check_planarity(G)
+    _, emb_p = nx.check_planarity(Gp)
+    Gp_faces = _get_all_faces(emb_p)
+    outer_set = set(outer_face)
+    Gp_outer = max(Gp_faces,
+                   key=lambda f: (len(set(f) & outer_set), -len(f)))
+    pos_G  = balanced_planar_layout(G,  outer_face, n_explore=4000,
+                                    n_refine=2000, seed=42 + src)
+    pos_Gp = balanced_planar_layout(Gp, Gp_outer,  n_explore=4000,
+                                    n_refine=2000, seed=42 + tgt)
+
+    # ── Draw G on the left (no 4-coloring; just structure) ────────────────
+    ax_G = axes[row, 0]
+    odd_cycles_G = get_odd_level_cycles(emb_G, levels_G)
+    ax_G.add_patch(plt.Polygon([pos_G[v] for v in outer_face], closed=True,
+                               facecolor='#E3F2FD', alpha=0.5,
+                               edgecolor='none', zorder=0))
+    for cycle in odd_cycles_G:
+        if set(cycle) == outer_set: continue
+        ax_G.add_patch(plt.Polygon([pos_G[v] for v in cycle], closed=True,
+                                   facecolor='#FFF9C4', alpha=0.6,
+                                   edgecolor='none', zorder=0))
+    flip_edges_G = set(frozenset([u, v]) for u, v, _, _ in choices)
+    for u, v in G.edges():
+        if frozenset([u, v]) in flip_edges_G:
+            nx.draw_networkx_edges(G, pos_G, edgelist=[(u, v)], ax=ax_G,
+                                   edge_color='#e53935', width=3.5,
+                                   style='dashed')
+        elif levels_G[u] == levels_G[v]:
+            nx.draw_networkx_edges(G, pos_G, edgelist=[(u, v)], ax=ax_G,
+                                   edge_color='#2196F3', width=2.5)
+        else:
+            nx.draw_networkx_edges(G, pos_G, edgelist=[(u, v)], ax=ax_G,
+                                   edge_color='#bdbdbd', width=1.5)
+    # Nodes by level (un-coloured)
+    level_palette = ['#1565C0', '#E65100', '#2E7D32']
+    node_lv_cols = [level_palette[min(levels_G[v], len(level_palette)-1)]
+                    for v in G.nodes()]
+    nx.draw_networkx_nodes(G, pos_G, ax=ax_G, node_color=node_lv_cols,
+                           node_size=550)
+    nx.draw_networkx_labels(G, pos_G, ax=ax_G, font_size=11,
+                            font_color='white', font_weight='bold')
+    for v, p in pos_G.items():
+        ax_G.annotate(f'L{levels_G[v]}', xy=p + np.array([0.07, 0.07]),
+                      fontsize=8, color='#333', zorder=5)
+    flip_str = ", ".join(f"({u},{v})→({w},{x})" for u, v, w, x in choices)
+    ax_G.set_title(f"G  (T{src})  |  outer face {outer_face}\n"
+                   f"flips: {flip_str}",
+                   fontsize=9, fontweight='bold', pad=8)
+    ax_G.set_xlim(-1.5, 1.5); ax_G.set_ylim(-1.5, 1.5)
+    ax_G.set_aspect('equal'); ax_G.axis('off')
+
+    # ── Draw G' on the right (4-coloured) ─────────────────────────────────
+    ax = axes[row, 1]
+    ax.add_patch(plt.Polygon([pos_Gp[v] for v in Gp_outer], closed=True,
+                             facecolor='#ECEFF1', alpha=0.6,
+                             edgecolor='none', zorder=0))
+
+    bad_edge_set = set()
+    for e in bad['even'] + bad['odd'] + bad['cross']:
+        bad_edge_set.add(frozenset(e))
+
+    for u, v in Gp.edges():
+        key = frozenset([u, v])
+        same_level = levels_G[u] == levels_G[v]
+        if key in bad_edge_set:
+            nx.draw_networkx_edges(Gp, pos_Gp, edgelist=[(u, v)], ax=ax,
+                                   edge_color='#000', width=4.0)
+        elif same_level:
+            nx.draw_networkx_edges(Gp, pos_Gp, edgelist=[(u, v)], ax=ax,
+                                   edge_color='#666', width=2.5)
+        else:
+            nx.draw_networkx_edges(Gp, pos_Gp, edgelist=[(u, v)], ax=ax,
+                                   edge_color='#bdbdbd', width=1.5)
+
+    node_colors_4 = [color_hex[coloring[v]] for v in Gp.nodes()]
+    nx.draw_networkx_nodes(Gp, pos_Gp, ax=ax,
+                           node_color=node_colors_4,
+                           edgecolors='black', linewidths=2,
+                           node_size=600)
+    nx.draw_networkx_labels(Gp, pos_Gp, ax=ax, font_size=11,
+                            font_color='white', font_weight='bold')
+    for v, p in pos_Gp.items():
+        ax.annotate(f'L{levels_G[v]}', xy=p + np.array([0.07, 0.07]),
+                    fontsize=8, color='#333', zorder=5)
+
+    title = (f"G'  (T{tgt})  |  layout outer {Gp_outer}\n"
+             f"valid 4-coloring: {valid}")
+    if not valid:
+        title += f"\nconflicts: even={bad['even']} odd={bad['odd']}"
+    ax.set_title(title, fontsize=9, fontweight='bold', pad=8)
+    ax.set_xlim(-1.5, 1.5); ax.set_ylim(-1.5, 1.5)
+    ax.set_aspect('equal'); ax.axis('off')
+
+legend_elements = [
+    mpatches.Patch(facecolor='#E3F2FD', edgecolor='none',
+                   label="G outer face"),
+    mpatches.Patch(facecolor='#FFF9C4', edgecolor='none',
+                   label='G odd level cycle'),
+    plt.Line2D([0],[0], color='#e53935', lw=3.0, ls='dashed',
+               label='Edge flipped out (in G)'),
+    plt.Line2D([0],[0], color='#2196F3', lw=2.5, label='Level edge'),
+    mpatches.Patch(facecolor=color_hex['red'],    edgecolor='black',
+                   label="G' red (even-level)"),
+    mpatches.Patch(facecolor=color_hex['blue'],   edgecolor='black',
+                   label="G' blue (even-level)"),
+    mpatches.Patch(facecolor=color_hex['yellow'], edgecolor='black',
+                   label="G' yellow (odd-level)"),
+    mpatches.Patch(facecolor=color_hex['green'],  edgecolor='black',
+                   label="G' green (odd-level)"),
+    plt.Line2D([0],[0], color='#000',    lw=4.0,
+               label='Coloring conflict edge'),
+]
+fig.legend(handles=legend_elements, loc='lower center', ncol=4,
+           fontsize=10, bbox_to_anchor=(0.5, -0.002))
+plt.suptitle("4-coloring of G' using levels of G  "
+             "(even=red/blue, odd=yellow/green)",
+             fontsize=13, fontweight='bold', y=1.0)
+plt.tight_layout()
+plt.savefig('/mnt/user-data/outputs/four_color.png',
+            dpi=130, bbox_inches='tight')
+print("\nSaved.")

papers/level_resolutions_of_maximal_planar_graphs/experiments/level_cycles.py

@@ -0,0 +1,173 @@
+"""
+Level cycle computation for maximal planar graphs.
+
+Definitions
+-----------
+- Level source: either (a) a face of G [all face vertices at level 0]
+                       or (b) a degree-3 vertex of G [singleton at level 0]
+- Levels: BFS distance from the level source.
+- Level cycle: simple face of the level subgraph L_k (subgraph induced by
+  level-k vertices) in the embedding inherited from G's planar embedding.
+- Level resolution: G' obtained from G by flipping
+    - exactly one edge per ODD level cycle (mandatory)
+    - at most one edge per EVEN level cycle (optional)
+"""
+import networkx as nx
+from collections import deque, defaultdict
+from itertools import product as iproduct
+
+
+def compute_levels(G, source_set):
+    """BFS levels from any iterable of source vertices."""
+    levels, queue = {}, deque()
+    for v in source_set:
+        levels[v] = 0
+        queue.append(v)
+    while queue:
+        v = queue.popleft()
+        for w in G.neighbors(v):
+            if w not in levels:
+                levels[w] = levels[v] + 1
+                queue.append(w)
+    return levels
+
+
+def get_all_faces(emb):
+    seen, faces = set(), []
+    for v in emb.nodes():
+        for w in emb[v]:
+            if (v, w) not in seen:
+                face = emb.traverse_face(v, w)
+                for i in range(len(face)):
+                    seen.add((face[i], face[(i+1) % len(face)]))
+                faces.append(tuple(face))
+    return faces
+
+
+def inherited_embedding(emb_G, sub_nodes):
+    """PlanarEmbedding of the induced subgraph with inherited cyclic order."""
+    sub_set = set(sub_nodes)
+    emb = nx.PlanarEmbedding()
+    for v in sub_nodes:
+        emb.add_node(v)
+    for v in sub_nodes:
+        cw = [w for w in emb_G.neighbors_cw_order(v) if w in sub_set]
+        prev = None
+        for w in cw:
+            if prev is None:
+                emb.add_half_edge_first(v, w)
+            else:
+                emb.add_half_edge_cw(v, w, prev)
+            prev = w
+    return emb
+
+
+def _is_simple(face):
+    return len(set(face)) == len(face)
+
+
+def get_level_cycles_by_parity(G, emb_G, levels):
+    """Simple faces of each L_k, split by length parity."""
+    by_level = defaultdict(list)
+    for v, lv in levels.items():
+        by_level[lv].append(v)
+    odd, even, seen = [], [], set()
+    for k, nodes_k in by_level.items():
+        if len(nodes_k) < 3:
+            continue
+        L_k = G.subgraph(nodes_k)
+        if L_k.number_of_edges() < 3:
+            continue
+        try:
+            emb_L = inherited_embedding(emb_G, nodes_k)
+            faces = get_all_faces(emb_L)
+        except Exception:
+            ip, emb_L = nx.check_planarity(L_k)
+            if not ip:
+                continue
+            faces = get_all_faces(emb_L)
+        for face in faces:
+            if not _is_simple(face):
+                continue
+            rots = [tuple(face[i:] + face[:i]) for i in range(len(face))] \
+                 + [tuple(face[::-1][i:] + face[::-1][:i])
+                    for i in range(len(face))]
+            canon = min(rots)
+            if canon in seen:
+                continue
+            seen.add(canon)
+            (odd if len(face) % 2 else even).append(face)
+    return odd, even
+
+
+def get_flip_candidates(G, emb_G, cycle):
+    """Valid edge flips along a level cycle."""
+    cands = []
+    n = len(cycle)
+    for i in range(n):
+        u, v = cycle[i], cycle[(i+1) % n]
+        if not G.has_edge(u, v):
+            continue
+        f1 = emb_G.traverse_face(u, v)
+        f2 = emb_G.traverse_face(v, u)
+        if len(f1) != 3 or len(f2) != 3:
+            continue
+        w = next(x for x in f1 if x != u and x != v)
+        x = next(y for y in f2 if y != u and y != v)
+        if w == x or G.has_edge(w, x):
+            continue
+        cands.append((u, v, w, x))
+    return cands
+
+
+def level_sources(G, emb_G):
+    """Yield every valid level source as (kind, label, vertex_set)."""
+    for face in get_all_faces(emb_G):
+        yield ('face', tuple(face), set(face))
+    for v in G.nodes():
+        if G.degree(v) == 3:
+            yield ('vertex', v, {v})
+
+
+def all_level_resolutions(G):
+    """Yield (G', source_kind, source_label, choices) for every resolution."""
+    is_planar, emb = nx.check_planarity(G)
+    if not is_planar:
+        return
+    for kind, label, source_set in level_sources(G, emb):
+        levels = compute_levels(G, source_set)
+        odd_cycles, even_cycles = get_level_cycles_by_parity(G, emb, levels)
+        if not odd_cycles and not even_cycles:
+            yield (G.copy(), kind, label, [])
+            continue
+        cand_lists, ok = [], True
+        for cycle in odd_cycles:
+            cands = get_flip_candidates(G, emb, cycle)
+            if not cands:
+                ok = False
+                break
+            cand_lists.append(cands)
+        if not ok:
+            continue
+        for cycle in even_cycles:
+            cands = get_flip_candidates(G, emb, cycle)
+            cand_lists.append(cands + [None])
+        for choices in iproduct(*cand_lists):
+            Gp = G.copy()
+            applied, ok2 = set(), True
+            for choice in choices:
+                if choice is None:
+                    continue
+                u, v, w, x = choice
+                if frozenset([u, v]) in applied:
+                    ok2 = False; break
+                if not Gp.has_edge(u, v) or Gp.has_edge(w, x):
+                    ok2 = False; break
+                Gp.remove_edge(u, v)
+                Gp.add_edge(w, x)
+                applied.add(frozenset([u, v]))
+            if not ok2:
+                continue
+            ip2, _ = nx.check_planarity(Gp)
+            if ip2:
+                yield (Gp, kind, label, [c for c in choices if c is not None])

papers/level_resolutions_of_maximal_planar_graphs/experiments/n11_chunk.py

@@ -0,0 +1,57 @@
+"""
+Process a chunk of source indices at n=11. Save accumulating results to JSON.
+Usage: python3 n11_chunk.py <start> <end>
+"""
+import sys, os, json, time
+sys.path.insert(0, '/home/claude/build')
+from coverage_fast import (
+    enumerate_all_triangulations_fast, build_iso_index,
+    resolution_classes_fast, min_degree,
+)
+import pickle
+
+STATE_PKL = '/tmp/n11_state.pkl'
+RESULTS_JSON = '/tmp/n11_results.json'
+
+# Load or build state
+if os.path.exists(STATE_PKL):
+    with open(STATE_PKL, 'rb') as f:
+        state = pickle.load(f)
+    reps = state['reps']; iso_index = state['iso_index']
+    print(f"Loaded {len(reps)} reps from cache")
+else:
+    t0 = time.time()
+    reps = enumerate_all_triangulations_fast(11)
+    iso_index = build_iso_index(reps)
+    print(f"Built {len(reps)} reps in {time.time()-t0:.1f}s")
+    with open(STATE_PKL, 'wb') as f:
+        pickle.dump({'reps': reps, 'iso_index': iso_index}, f)
+
+# Load partial results
+if os.path.exists(RESULTS_JSON):
+    with open(RESULTS_JSON) as f:
+        results = {int(k): v for k, v in json.load(f).items()}
+else:
+    results = {}
+
+start = int(sys.argv[1]); end = int(sys.argv[2])
+end = min(end, len(reps))
+t1 = time.time()
+for i in range(start, end):
+    if i in results:
+        continue
+    res = resolution_classes_fast(reps[i], iso_index)
+    results[i] = sorted(res)
+    if (i + 1) % 25 == 0:
+        elapsed = time.time() - t1
+        done = i + 1 - start
+        eta = elapsed * (end - start - done) / max(done, 1)
+        print(f"  {i+1}/{end} ({elapsed:.0f}s, ETA {eta:.0f}s)", flush=True)
+    if (i + 1) % 50 == 0:
+        with open(RESULTS_JSON, 'w') as f:
+            json.dump({str(k): v for k, v in results.items()}, f)
+
+with open(RESULTS_JSON, 'w') as f:
+    json.dump({str(k): v for k, v in results.items()}, f)
+print(f"Saved {len(results)} source results. Range [{start},{end}). "
+      f"Chunk time: {time.time()-t1:.1f}s")

papers/level_resolutions_of_maximal_planar_graphs/experiments/n7_examples.py

@@ -0,0 +1,251 @@
+"""
+Generate G → G' examples covering all 5 iso-classes on n=7,
+each as source and as target. Uses balanced planar layout.
+"""
+import sys
+sys.path.insert(0, '/home/claude')
+
+import networkx as nx
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+import numpy as np
+from collections import deque
+from itertools import combinations, product as iproduct
+
+from balanced_layout import balanced_planar_layout, _get_all_faces
+
+# ── Helpers ──────────────────────────────────────────────────────────────────
+
+def compute_levels(G, outer_face):
+    levels, queue = {}, deque()
+    for v in outer_face:
+        levels[v] = 0; queue.append(v)
+    while queue:
+        v = queue.popleft()
+        for w in G.neighbors(v):
+            if w not in levels:
+                levels[w] = levels[v] + 1; queue.append(w)
+    return levels
+
+def get_odd_level_cycles(emb, levels):
+    odd = []
+    for face in _get_all_faces(emb):
+        lv = {levels.get(v) for v in face}
+        if len(lv) == 1 and None not in lv and len(face) % 2 == 1:
+            odd.append(face)
+    return odd
+
+def get_flip_candidates(G, emb, cycle):
+    cands = []
+    n = len(cycle)
+    for i in range(n):
+        u, v = cycle[i], cycle[(i+1) % n]
+        if not G.has_edge(u, v): continue
+        f1 = emb.traverse_face(u, v); f2 = emb.traverse_face(v, u)
+        if len(f1) != 3 or len(f2) != 3: continue
+        w = next(x for x in f1 if x != u and x != v)
+        x = next(y for y in f2 if y != u and y != v)
+        if w == x or G.has_edge(w, x): continue
+        cands.append((u, v, w, x))
+    return cands
+
+def is_triangulation(G):
+    n = G.number_of_nodes()
+    if G.number_of_edges() != 3*n - 6: return False
+    if not nx.is_connected(G): return False
+    ip, emb = nx.check_planarity(G)
+    if not ip: return False
+    seen = set()
+    for v in emb.nodes():
+        for w in emb[v]:
+            if (v, w) not in seen:
+                f = emb.traverse_face(v, w)
+                for i in range(len(f)):
+                    seen.add((f[i], f[(i+1) % len(f)]))
+                if len(f) != 3: return False
+    return True
+
+# ── Enumerate iso-class representatives for n=7 ──────────────────────────────
+
+print("Enumerating triangulations on 7 vertices...")
+nodes = list(range(7))
+all_edges = list(combinations(nodes, 2))
+reps = []
+for sub in combinations(all_edges, 15):
+    G = nx.Graph(); G.add_nodes_from(nodes); G.add_edges_from(sub)
+    if is_triangulation(G):
+        if all(not nx.is_isomorphic(G, r) for r in reps):
+            reps.append(G)
+        if len(reps) == 5: break
+print(f"Found {len(reps)} reps.")
+for i, G in enumerate(reps):
+    deg_seq = sorted((d for _, d in G.degree()), reverse=True)
+    print(f"  T{i}: degree {deg_seq}")
+
+def iso_class(G):
+    for i, r in enumerate(reps):
+        if nx.is_isomorphic(G, r): return i
+    return -1
+
+# ── Find one (G_source, outer_face, flip, G_target) for each desired pair ───
+
+desired_pairs = [
+    (0, 2),  # T0 -> T2
+    (1, 4),  # T1 -> T4
+    (2, 0),  # T2 -> T0
+    (3, 1),  # T3 -> T1
+    (4, 3),  # T4 -> T3
+]
+
+def find_example(src_idx, tgt_idx):
+    G = reps[src_idx]
+    ip, emb = nx.check_planarity(G)
+    for outer_face in _get_all_faces(emb):
+        levels = compute_levels(G, outer_face)
+        odd_cycles = get_odd_level_cycles(emb, levels)
+        if not odd_cycles: continue
+        candidate_lists = []
+        ok = True
+        for cycle in odd_cycles:
+            cands = get_flip_candidates(G, emb, cycle)
+            if not cands: ok = False; break
+            candidate_lists.append(cands)
+        if not ok: continue
+        for choices in iproduct(*candidate_lists):
+            Gp = G.copy()
+            valid = True
+            for u, v, w, x in choices:
+                if not Gp.has_edge(u, v) or Gp.has_edge(w, x):
+                    valid = False; break
+                Gp.remove_edge(u, v); Gp.add_edge(w, x)
+            if valid and is_triangulation(Gp) and iso_class(Gp) == tgt_idx:
+                return G, outer_face, choices, Gp
+    return None
+
+print("\nFinding example flips for each pair:")
+examples = []
+for src, tgt in desired_pairs:
+    ex = find_example(src, tgt)
+    if ex is None:
+        print(f"  T{src} -> T{tgt}: NO EXAMPLE FOUND")
+        continue
+    G, outer_face, choices, Gp = ex
+    print(f"  T{src} -> T{tgt}: outer={outer_face}, flips={list(choices)}")
+    examples.append((src, tgt, G, outer_face, choices, Gp))
+
+# ── Layout and plot ──────────────────────────────────────────────────────────
+
+print("\nComputing balanced layouts (this takes a minute)...")
+fig, axes = plt.subplots(len(examples), 2, figsize=(13, 5.5 * len(examples)))
+if len(examples) == 1:
+    axes = np.array([axes])
+
+palette = ['#1565C0', '#E65100', '#2E7D32']
+
+for row, (src, tgt, G, outer_face, choices, Gp) in enumerate(examples):
+    levels = compute_levels(G, outer_face)
+    _, emb = nx.check_planarity(G)
+    _, emb_p = nx.check_planarity(Gp)
+    odd_cycles_G = get_odd_level_cycles(emb, levels)
+    # Re-compute levels for Gp using the same outer face (same level-0 vertices)
+    levels_p = compute_levels(Gp, outer_face)
+    odd_cycles_Gp = get_odd_level_cycles(emb_p, levels_p)
+
+    pos_G  = balanced_planar_layout(G,  outer_face, n_explore=4000,
+                                    n_refine=2000, seed=42 + src)
+
+    # G' may have a different planar embedding; the original outer face
+    # might no longer be a face. Pick a face of G' that maximally overlaps
+    # the original outer face vertices.
+    Gp_faces = _get_all_faces(emb_p)
+    outer_set_orig = set(outer_face)
+    Gp_outer = max(Gp_faces,
+                   key=lambda f: (len(set(f) & outer_set_orig), -len(f)))
+    pos_Gp = balanced_planar_layout(Gp, Gp_outer, n_explore=4000,
+                                    n_refine=2000, seed=42 + tgt)
+    print(f"  T{src}→T{tgt}: layouts done")
+
+    flip_edges = set(frozenset([u, v]) for u, v, _, _ in choices)
+    new_edges  = set(frozenset([w, x]) for _, _, w, x in choices)
+
+    def classify_G(e):
+        if frozenset(e) in flip_edges: return 'flip'
+        eu, ev = e
+        if levels.get(eu) == levels.get(ev): return 'level'
+        return 'normal'
+
+    def classify_Gp(e):
+        if frozenset(e) in new_edges: return 'new'
+        eu, ev = e
+        if levels_p.get(eu) == levels_p.get(ev): return 'level'
+        return 'normal'
+
+    outer_set = set(outer_face)
+
+    def draw(ax, Gd, pos, levels_d, odd_cycles_d, classifier, of, title):
+        # Outer face
+        ax.add_patch(plt.Polygon([pos[nd] for nd in of],
+                                 closed=True, facecolor='#E3F2FD',
+                                 alpha=0.5, edgecolor='none', zorder=0))
+        # Interior odd cycles
+        of_set = set(of)
+        for cycle in odd_cycles_d:
+            if set(cycle) == of_set: continue
+            ax.add_patch(plt.Polygon([pos[nd] for nd in cycle],
+                                     closed=True, facecolor='#FFF9C4',
+                                     alpha=0.6, edgecolor='none', zorder=0))
+        # Edges
+        cats = {'normal': [], 'level': [], 'flip': [], 'new': []}
+        for e in Gd.edges():
+            cats[classifier(e)].append(e)
+        style = {
+            'normal': dict(edge_color='#bdbdbd', width=1.5),
+            'level':  dict(edge_color='#2196F3', width=2.5),
+            'flip':   dict(edge_color='#e53935', width=3.5, style='dashed'),
+            'new':    dict(edge_color='#43a047', width=3.5),
+        }
+        for cat, eds in cats.items():
+            if eds:
+                nx.draw_networkx_edges(Gd, pos, edgelist=eds, ax=ax,
+                                       **style[cat])
+        # Nodes coloured by level
+        node_colors = [palette[min(levels_d.get(nd, 0), len(palette)-1)]
+                       for nd in Gd.nodes()]
+        nx.draw_networkx_nodes(Gd, pos, ax=ax, node_color=node_colors,
+                               node_size=450)
+        nx.draw_networkx_labels(Gd, pos, ax=ax, font_size=11,
+                                font_color='white', font_weight='bold')
+        for nd, p in pos.items():
+            lv = levels_d.get(nd, '?')
+            ax.annotate(f'L{lv}', xy=p + np.array([0.06, 0.06]),
+                        fontsize=7, color='#444', zorder=5)
+        ax.set_title(title, fontsize=10, fontweight='bold', pad=8)
+        ax.set_xlim(-1.5, 1.5); ax.set_ylim(-1.5, 1.5)
+        ax.set_aspect('equal'); ax.axis('off')
+
+    flip_str = ", ".join(f"({u},{v})→({w},{x})" for u,v,w,x in choices)
+    draw(axes[row, 0], G,  pos_G,  levels,   odd_cycles_G,  classify_G,
+         outer_face,
+         f"G  (T{src})  |  outer face {outer_face}\nflip: {flip_str}")
+    draw(axes[row, 1], Gp, pos_Gp, levels_p, odd_cycles_Gp, classify_Gp,
+         Gp_outer,
+         f"G'  (T{tgt})  |  outer face {Gp_outer}")
+
+legend_elements = [
+    mpatches.Patch(facecolor='#E3F2FD', edgecolor='none', label='Outer face'),
+    mpatches.Patch(facecolor='#FFF9C4', edgecolor='none', label='Odd level cycle'),
+    plt.Line2D([0],[0], color='#bdbdbd', lw=1.5, label='Normal edge'),
+    plt.Line2D([0],[0], color='#2196F3', lw=2.5, label='Level edge'),
+    plt.Line2D([0],[0], color='#e53935', lw=3.0, ls='dashed',
+               label='Removed edge'),
+    plt.Line2D([0],[0], color='#43a047', lw=3.0, label='Added edge'),
+]
+fig.legend(handles=legend_elements, loc='lower center', ncol=6,
+           fontsize=10, bbox_to_anchor=(0.5, -0.005))
+plt.suptitle(
+    "Level resolutions on n=7  |  every iso-class appears as source and target",
+    fontsize=14, fontweight='bold', y=1.0)
+plt.tight_layout()
+plt.savefig('/mnt/user-data/outputs/n7_all_classes.png',
+            dpi=130, bbox_inches='tight')
+print("\nSaved.")

papers/level_resolutions_of_maximal_planar_graphs/experiments/orbit_check.py

@@ -0,0 +1,123 @@
+"""
+Extended orbit-counting framework: REVERSE-FLIPS ON ANY SET OF EDGES.
+
+For each md4 target G' and each k ∈ {1, 2, ...}:
+- Enumerate distinct (up to iso) k-flip reverse-preimages G of G'
+- For each, check if G level-resolves to G' via any source
+
+Surjectivity at G' = some k and some preimage works.
+
+Includes a critical test: the icosahedron at n=12 (no deg-4 vertices).
+"""
+import sys; sys.path.insert(0, '/home/claude/build')
+import networkx as nx
+import time
+from collections import defaultdict
+from itertools import combinations, product as iproduct
+
+from level_cycles import (
+    compute_levels, get_all_faces, get_level_cycles_by_parity,
+    get_flip_candidates, level_sources,
+)
+
+
+def reverse_flip_edge(G_prime, e):
+    u, v = tuple(e)
+    if not G_prime.has_edge(u, v): return None
+    _, emb = nx.check_planarity(G_prime)
+    f1 = emb.traverse_face(u, v); f2 = emb.traverse_face(v, u)
+    if len(f1) != 3 or len(f2) != 3: return None
+    a = next(x for x in f1 if x != u and x != v)
+    b = next(x for x in f2 if x != u and x != v)
+    if a == b or G_prime.has_edge(a, b): return None
+    G = G_prime.copy()
+    G.remove_edge(u, v); G.add_edge(a, b)
+    ip, _ = nx.check_planarity(G)
+    if not ip: return None
+    return G
+
+
+def k_flip_preimages_iso(G_prime, k):
+    """Yield distinct (up to iso) triangulations obtainable from G_prime by
+    k successive reverse-flips."""
+    seen = [G_prime]
+    frontier = [G_prime]
+    for step in range(k):
+        new_frontier = []
+        new_seen = []
+        for H in frontier:
+            for e in list(H.edges()):
+                Hp = reverse_flip_edge(H, frozenset(e))
+                if Hp is None: continue
+                if any(nx.is_isomorphic(Hp, S) for S in seen + new_seen):
+                    continue
+                new_seen.append(Hp)
+                new_frontier.append(Hp)
+        seen.extend(new_seen)
+        frontier = new_frontier
+        if not frontier: break
+    return frontier  # candidates at depth EXACTLY k
+
+
+def resolves_to_target(G, target):
+    """Does G level-resolve to iso(target) via any source?"""
+    ip, emb = nx.check_planarity(G)
+    if not ip: return False
+    for kind, label, source_set in level_sources(G, emb):
+        levels = compute_levels(G, source_set)
+        odd_cycles, even_cycles = get_level_cycles_by_parity(G, emb, levels)
+        cand_lists, ok = [], True
+        for cycle in odd_cycles:
+            cands = get_flip_candidates(G, emb, cycle)
+            if not cands: ok = False; break
+            cand_lists.append(cands)
+        if not ok: continue
+        for cycle in even_cycles:
+            cands = get_flip_candidates(G, emb, cycle)
+            cand_lists.append(cands + [None])
+        if not cand_lists:
+            if nx.is_isomorphic(G, target): return True
+            continue
+        for choices in iproduct(*cand_lists):
+            Gp = G.copy(); applied = set(); ok2 = True
+            for c in choices:
+                if c is None: continue
+                u, v, w, x = c
+                if frozenset([u,v]) in applied: ok2 = False; break
+                if not Gp.has_edge(u, v) or Gp.has_edge(w, x): ok2 = False; break
+                Gp.remove_edge(u, v); Gp.add_edge(w, x)
+                applied.add(frozenset([u, v]))
+            if not ok2: continue
+            ipp, _ = nx.check_planarity(Gp)
+            if not ipp: continue
+            if nx.is_isomorphic(Gp, target): return True
+    return False
+
+
+def check_target(G_prime, k_max=3, label=""):
+    """For target G_prime, scan k=1, 2, ..., k_max for a preimage that works.
+    Return (success_k, preimage_idx) or (None, None) if all fail."""
+    for k in range(1, k_max + 1):
+        t0 = time.time()
+        cands = k_flip_preimages_iso(G_prime, k)
+        n_cand = len(cands)
+        for i, G in enumerate(cands):
+            if resolves_to_target(G, G_prime):
+                elapsed = time.time() - t0
+                print(f"  {label}: k={k}, hit at {i+1}/{n_cand} "
+                      f"({elapsed:.1f}s)")
+                return k, i
+        elapsed = time.time() - t0
+        print(f"  {label}: k={k}, no hits among {n_cand} preimages "
+              f"({elapsed:.1f}s)")
+    return None, None
+
+
+if __name__ == "__main__":
+    # Test on icosahedron
+    G_ico = nx.icosahedral_graph()
+    print("=" * 60)
+    print("Icosahedron (n=12, 5-regular)")
+    print("=" * 60)
+    k, idx = check_target(G_ico, k_max=3, label="icosahedron")
+    print(f"Result: reachable at k={k}" if k else "NOT reachable up to k=3")

papers/level_resolutions_of_maximal_planar_graphs/experiments/outerplanar_check.py

@@ -0,0 +1,55 @@
+"""Verify that every level subgraph L_k is outerplanar across n=6..10."""
+import sys; sys.path.insert(0, '/home/claude/build')
+import networkx as nx
+from collections import defaultdict
+from triangulation_gen import enumerate_all_triangulations
+from level_cycles import compute_levels, level_sources
+
+
+def is_outerplanar(G):
+    """Outerplanar iff no K_4 minor and no K_{2,3} minor.
+    Equivalently: planar AND can add a new vertex connected to all existing
+    vertices while remaining planar."""
+    if G.number_of_nodes() <= 3:
+        return True
+    H = G.copy()
+    apex = max(H.nodes()) + 1
+    for v in G.nodes():
+        H.add_edge(apex, v)
+    return nx.check_planarity(H)[0]
+
+
+# Check on a sample of triangulations
+results = defaultdict(int)  # (n, outerplanar) -> count
+non_outerplanar_examples = []
+
+for n in [6, 7, 8, 9, 10]:
+    reps = enumerate_all_triangulations(n)
+    for gi, G in enumerate(reps):
+        ip, emb = nx.check_planarity(G)
+        if not ip: continue
+        for kind, label, source_set in level_sources(G, emb):
+            levels = compute_levels(G, source_set)
+            level_groups = defaultdict(list)
+            for v, lv in levels.items():
+                level_groups[lv].append(v)
+            for k, verts in level_groups.items():
+                Lk = G.subgraph(verts)
+                op = is_outerplanar(Lk)
+                results[(n, op)] += 1
+                if not op and len(non_outerplanar_examples) < 5:
+                    non_outerplanar_examples.append({
+                        'n': n, 'G_idx': gi,
+                        'source': (kind, label),
+                        'level': k,
+                        'L_k_edges': list(Lk.edges()),
+                        'L_k_nodes': verts,
+                    })
+                    print(f"  Non-outerplanar L_{k} at n={n}, T{gi}, "
+                          f"source={kind} {label}: {verts}, edges={list(Lk.edges())}")
+
+print(f"\nSummary:")
+for n in [6, 7, 8, 9, 10]:
+    op = results[(n, True)]
+    nop = results[(n, False)]
+    print(f"  n={n}: {op} outerplanar, {nop} non-outerplanar level subgraphs")

papers/level_resolutions_of_maximal_planar_graphs/experiments/plot_oct.py

@@ -0,0 +1,104 @@
+import networkx as nx
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+import numpy as np
+from collections import deque
+
+def compute_levels(G, outer_face):
+    levels = {}
+    queue = deque()
+    for v in outer_face:
+        levels[v] = 0
+        queue.append(v)
+    while queue:
+        v = queue.popleft()
+        for w in G.neighbors(v):
+            if w not in levels:
+                levels[w] = levels[v] + 1
+                queue.append(w)
+    return levels
+
+def get_all_faces(emb):
+    seen = set()
+    faces = []
+    for v in emb.nodes():
+        for w in emb[v]:
+            if (v, w) not in seen:
+                face = emb.traverse_face(v, w)
+                for i in range(len(face)):
+                    seen.add((face[i], face[(i+1) % len(face)]))
+                faces.append(tuple(face))
+    return faces
+
+G = nx.octahedral_graph()
+_, emb = nx.check_planarity(G)
+all_faces = get_all_faces(emb)
+
+pos = nx.spring_layout(G, seed=7, iterations=200)
+
+fig, axes = plt.subplots(2, 4, figsize=(16, 8))
+axes = axes.flatten()
+
+for idx, outer_face in enumerate(all_faces):
+    ax = axes[idx]
+    outer_face_set = set(outer_face)
+    levels = compute_levels(G, outer_face)
+
+    odd_cycle = None
+    for face in all_faces:
+        if set(face) == outer_face_set:
+            continue
+        face_levels = set(levels.get(v) for v in face)
+        if len(face_levels) == 1 and None not in face_levels and len(face) % 2 == 1:
+            odd_cycle = face
+            break
+
+    odd_edges = set()
+    if odd_cycle:
+        n = len(odd_cycle)
+        for i in range(n):
+            odd_edges.add(frozenset([odd_cycle[i], odd_cycle[(i+1) % n]]))
+
+    outer_edges = set()
+    n_of = len(outer_face)
+    for i in range(n_of):
+        outer_edges.add(frozenset([outer_face[i], outer_face[(i+1) % n_of]]))
+
+    red_edges  = [e for e in G.edges() if frozenset(e) in odd_edges]
+    blue_edges = [e for e in G.edges() if frozenset(e) in outer_edges]
+    grey_edges = [e for e in G.edges()
+                  if frozenset(e) not in odd_edges and frozenset(e) not in outer_edges]
+
+    node_colors = ['#4a90d9' if levels[v] == 0 else '#f0a500' for v in G.nodes()]
+
+    if odd_cycle:
+        tri = np.array([pos[v] for v in odd_cycle])
+        ax.add_patch(plt.Polygon(tri, closed=True, facecolor='#e03030', alpha=0.15, edgecolor='none'))
+
+    out_arr = np.array([pos[v] for v in outer_face])
+    ax.add_patch(plt.Polygon(out_arr, closed=True, facecolor='#4a90d9', alpha=0.15, edgecolor='none'))
+
+    nx.draw_networkx_edges(G, pos, edgelist=grey_edges, ax=ax, edge_color='#aaaaaa', width=1.5)
+    nx.draw_networkx_edges(G, pos, edgelist=blue_edges, ax=ax, edge_color='#4a90d9', width=2.5, style='dashed')
+    nx.draw_networkx_edges(G, pos, edgelist=red_edges,  ax=ax, edge_color='#e03030', width=3.5)
+    nx.draw_networkx_nodes(G, pos, ax=ax, node_color=node_colors, node_size=350)
+    nx.draw_networkx_labels(G, pos, ax=ax, font_size=9, font_color='white', font_weight='bold')
+
+    ax.set_title(f"outer={outer_face}\nodd cycle={odd_cycle}", fontsize=8, pad=4)
+    ax.axis('off')
+
+legend_elements = [
+    mpatches.Patch(facecolor='#4a90d9', alpha=0.5, label='Level 0 (outer face)'),
+    mpatches.Patch(facecolor='#f0a500', alpha=0.8, label='Level 1 node'),
+    mpatches.Patch(facecolor='#e03030', alpha=0.5, label='Odd level cycle (level 1)'),
+    plt.Line2D([0],[0], color='#4a90d9', linewidth=2.5, linestyle='dashed', label='Outer face edges'),
+    plt.Line2D([0],[0], color='#e03030', linewidth=3.5, label='Odd cycle edges'),
+    plt.Line2D([0],[0], color='#aaaaaa', linewidth=1.5, label='Other edges'),
+]
+fig.legend(handles=legend_elements, loc='lower center', ncol=3,
+           fontsize=9, bbox_to_anchor=(0.5, -0.01))
+fig.suptitle("Octahedron: odd level cycle (red) for each outer face",
+             fontsize=13, fontweight='bold')
+plt.tight_layout()
+plt.savefig('/mnt/user-data/outputs/octahedron_level_cycles.png', dpi=150, bbox_inches='tight')
+print("Saved.")

papers/level_resolutions_of_maximal_planar_graphs/experiments/triangulation_gen.py

@@ -0,0 +1,78 @@
+"""
+Complete enumeration of triangulations via flip closure.
+Start from stacked triangulations (vertex insertion), then apply all
+possible edge flips until closed under isomorphism. The flip graph of
+triangulations on n labeled vertices is known to be connected, so this
+covers all iso-classes.
+"""
+import networkx as nx
+from collections import deque
+import time
+
+def get_all_faces(emb):
+    seen, faces = set(), []
+    for v in emb.nodes():
+        for w in emb[v]:
+            if (v, w) not in seen:
+                face = emb.traverse_face(v, w)
+                for i in range(len(face)):
+                    seen.add((face[i], face[(i+1) % len(face)]))
+                faces.append(tuple(face))
+    return faces
+
+def insert_vertex_in_face(G, face, new_v):
+    Gp = G.copy(); Gp.add_node(new_v)
+    for v in face: Gp.add_edge(new_v, v)
+    return Gp
+
+def edge_flips(G):
+    """Yield all triangulations obtainable by a single edge flip from G."""
+    ip, emb = nx.check_planarity(G)
+    if not ip: return
+    yielded_signatures = set()
+    for u, v in G.edges():
+        f1 = emb.traverse_face(u, v); f2 = emb.traverse_face(v, u)
+        if len(f1) != 3 or len(f2) != 3: continue
+        w = next(x for x in f1 if x != u and x != v)
+        x = next(y for y in f2 if y != u and y != v)
+        if w == x or G.has_edge(w, x): continue
+        Gp = G.copy()
+        Gp.remove_edge(u, v); Gp.add_edge(w, x)
+        # signature for dedup of yielded flips
+        sig = frozenset(frozenset(e) for e in Gp.edges())
+        if sig in yielded_signatures: continue
+        yielded_signatures.add(sig)
+        yield Gp
+
+def enumerate_all_triangulations(n):
+    """All non-isomorphic triangulations on n vertices."""
+    if n < 4: return []
+    # Seed with stacked triangulations
+    G = nx.complete_graph(4)
+    current = [G]
+    for k in range(4, n):
+        next_set = []
+        for T in current:
+            _, emb = nx.check_planarity(T)
+            for face in get_all_faces(emb):
+                Tp = insert_vertex_in_face(T, face, k)
+                if all(not nx.is_isomorphic(Tp, q) for q in next_set):
+                    next_set.append(Tp)
+        current = next_set
+    # Flip closure
+    reps = list(current)
+    frontier = list(current)
+    while frontier:
+        new_frontier = []
+        for T in frontier:
+            for Tp in edge_flips(T):
+                if all(not nx.is_isomorphic(Tp, r) for r in reps):
+                    reps.append(Tp); new_frontier.append(Tp)
+        frontier = new_frontier
+    return reps
+
+if __name__ == "__main__":
+    for n in [4, 5, 6, 7, 8]:
+        t0 = time.time()
+        tris = enumerate_all_triangulations(n)
+        print(f"n={n}: {len(tris)} triangulations ({time.time()-t0:.1f}s)")

papers/level_resolutions_of_maximal_planar_graphs/experiments/triangulation_gen_fast.py

@@ -0,0 +1,100 @@
+"""
+Faster triangulation enumeration using Weisfeiler-Lehman graph hash as an
+isomorphism pre-filter. Iso-check is only run when WL hashes match.
+"""
+import networkx as nx
+from collections import defaultdict
+
+
+def get_all_faces(emb):
+    seen, faces = set(), []
+    for v in emb.nodes():
+        for w in emb[v]:
+            if (v, w) not in seen:
+                face = emb.traverse_face(v, w)
+                for i in range(len(face)):
+                    seen.add((face[i], face[(i+1) % len(face)]))
+                faces.append(tuple(face))
+    return faces
+
+
+def insert_vertex_in_face(G, face, new_v):
+    Gp = G.copy(); Gp.add_node(new_v)
+    for v in face: Gp.add_edge(new_v, v)
+    return Gp
+
+
+def edge_flips(G):
+    ip, emb = nx.check_planarity(G)
+    if not ip: return
+    yielded = set()
+    for u, v in G.edges():
+        f1 = emb.traverse_face(u, v); f2 = emb.traverse_face(v, u)
+        if len(f1) != 3 or len(f2) != 3: continue
+        w = next(x for x in f1 if x != u and x != v)
+        x = next(y for y in f2 if y != u and y != v)
+        if w == x or G.has_edge(w, x): continue
+        Gp = G.copy()
+        Gp.remove_edge(u, v); Gp.add_edge(w, x)
+        sig = frozenset(frozenset(e) for e in Gp.edges())
+        if sig in yielded: continue
+        yielded.add(sig)
+        yield Gp
+
+
+def wl_hash(G):
+    return nx.weisfeiler_lehman_graph_hash(G)
+
+
+class IsoBucket:
+    """Maintain a set of graphs deduplicated by isomorphism, using WL hash
+    as pre-filter."""
+    def __init__(self):
+        self.by_hash = defaultdict(list)
+        self.reps = []
+
+    def add(self, G):
+        h = wl_hash(G)
+        for R in self.by_hash[h]:
+            if nx.is_isomorphic(G, R):
+                return False
+        self.by_hash[h].append(G)
+        self.reps.append(G)
+        return True
+
+
+def enumerate_all_triangulations_fast(n, verbose=False):
+    if n < 4: return []
+    G = nx.complete_graph(4)
+    current = [G]
+    for k in range(4, n):
+        bucket = IsoBucket()
+        for T in current:
+            _, emb = nx.check_planarity(T)
+            for face in get_all_faces(emb):
+                bucket.add(insert_vertex_in_face(T, face, k))
+        current = list(bucket.reps)
+        if verbose: print(f"  after vertex-insertion to {k+1}: {len(current)}")
+    bucket = IsoBucket()
+    for T in current: bucket.add(T)
+    frontier = list(bucket.reps)
+    rounds = 0
+    while frontier:
+        rounds += 1
+        new_frontier = []
+        for T in frontier:
+            for Tp in edge_flips(T):
+                if bucket.add(Tp):
+                    new_frontier.append(Tp)
+        if verbose:
+            print(f"  flip-round {rounds}: total {len(bucket.reps)}")
+        frontier = new_frontier
+    return bucket.reps
+
+
+if __name__ == "__main__":
+    import time
+    for n in [9, 10, 11]:
+        t0 = time.time()
+        tris = enumerate_all_triangulations_fast(n, verbose=(n >= 11))
+        print(f"n={n}: {len(tris)} triangulations ({time.time()-t0:.1f}s)")