math-research/papers/face_monochromatic_pairs/experiments/search_kempe_property.py

"""Search for a proper 3-edge-coloring of a reduced dual satisfying:

  (i)   color(spike) == color(merged)
  (ii)  the {c_spike, c_side_0}-Kempe cycle through the spike contains
        side_0 AND merged
  (iii) the {c_spike, c_side_1}-Kempe cycle through the spike contains
        side_1 AND merged

Iterates over all min-degree-5 simple planar triangulations on n vertices
(via plantri through Sage), then every pentagonal face of the dual G', then
every index i in {0,...,4} for the reduced-dual construction, then every
proper 3-edge-coloring of the resulting reduced dual. Stops on the first hit.

The icosahedron (n = 12) is the only n = 12 triangulation with min degree 5
and was already known to fail (see check_dodecahedron_kempe.py).

Run with:
    sage experiments/search_kempe_property.py
"""
from sage.all import Graph
from sage.graphs.graph_generators import graphs
import sys
import time


def dual_of(G):
    """Combinatorial dual of a planar G with the embedding set on G."""
    faces = G.faces()
    edge_to_faces = {}
    for fi, face in enumerate(faces):
        for u, v in face:
            e = frozenset((u, v))
            edge_to_faces.setdefault(e, []).append(fi)
    dual_edges = []
    for e, fs in edge_to_faces.items():
        if len(fs) == 2:
            dual_edges.append((fs[0], fs[1]))
    return Graph(dual_edges, multiedges=False, loops=False)


def apply_reduction(G, face, i):
    """Apply Definition 2.1 at the pentagonal `face` (list of directed edges)
    with index i, returning (H, named) or (None, None) if the construction
    breaks."""
    boundary = [u for (u, v) in face]
    if len(set(boundary)) != 5:
        return None, None
    A = []
    for B_k in boundary:
        outer = [w for w in G.neighbor_iterator(B_k) if w not in boundary]
        if len(outer) != 1:
            return None, None
        A.append(outer[0])
    if len(set(A)) != 5:
        return None, None
    H = G.copy()
    for v in boundary:
        H.delete_vertex(v)
    v_n = '__v_n__'
    H.add_vertex(v_n)
    side_0 = (v_n, A[i % 5])
    spike  = (v_n, A[(i + 1) % 5])
    side_1 = (v_n, A[(i + 2) % 5])
    merged = (A[(i + 3) % 5], A[(i + 4) % 5])
    # If merged would be a self-loop or duplicate of an existing edge, skip
    if A[(i + 3) % 5] == A[(i + 4) % 5]:
        return None, None
    H.add_edges([side_0, spike, side_1, merged])
    if H.has_multiple_edges() or H.has_loops():
        return None, None
    if not H.is_planar(set_embedding=True):
        return None, None
    # Cubic check
    if not all(H.degree(v) == 3 for v in H.vertex_iterator()):
        return None, None
    named = {
        'spike':  frozenset(spike),
        'side_0': frozenset(side_0),
        'side_1': frozenset(side_1),
        'merged': frozenset(merged),
    }
    return H, named


def proper_3_edge_colorings(G):
    edges = list(G.edges(labels=False))
    n_edges = len(edges)
    adj = [[] for _ in range(n_edges)]
    for i in range(n_edges):
        u, v = edges[i][0], edges[i][1]
        for j in range(i):
            x, y = edges[j][0], edges[j][1]
            if u in (x, y) or v in (x, y):
                adj[i].append(j)
                adj[j].append(i)
    coloring = [-1] * n_edges

    def back(k):
        if k == n_edges:
            yield tuple(coloring)
            return
        for c in range(3):
            if all(coloring[j] != c for j in adj[k]):
                coloring[k] = c
                yield from back(k + 1)
        coloring[k] = -1

    return edges, back(0)


def kempe_cycle(edges, coloring, start_idx, color_pair):
    a, b = color_pair
    if coloring[start_idx] not in (a, b):
        return None
    in_sub = [i for i in range(len(edges)) if coloring[i] in (a, b)]
    visited = {start_idx}
    stack = [start_idx]
    while stack:
        cur = stack.pop()
        u, v = edges[cur][0], edges[cur][1]
        for j in in_sub:
            if j in visited:
                continue
            x, y = edges[j][0], edges[j][1]
            if u in (x, y) or v in (x, y):
                visited.add(j)
                stack.append(j)
    return visited


def matches(edges, col, named_idx):
    c_spike  = col[named_idx['spike']]
    c_merged = col[named_idx['merged']]
    if c_spike != c_merged:
        return False
    c_s0 = col[named_idx['side_0']]
    c_s1 = col[named_idx['side_1']]
    kc0 = kempe_cycle(edges, col, named_idx['spike'], (c_spike, c_s0))
    if named_idx['side_0'] not in kc0 or named_idx['merged'] not in kc0:
        return False
    kc1 = kempe_cycle(edges, col, named_idx['spike'], (c_spike, c_s1))
    if named_idx['side_1'] not in kc1 or named_idx['merged'] not in kc1:
        return False
    return True


def report(edges, col, named_idx, named, info):
    print()
    print("*** MATCH FOUND ***")
    for k, v in info.items():
        print(f"  {k} = {v}")
    print("Named edges and their colours:")
    for role in ('spike', 'side_0', 'side_1', 'merged'):
        e = edges[named_idx[role]]
        print(f"  {role:7s}: edge {tuple(e)}, colour {col[named_idx[role]]}")
    print(f"Full coloring (indexed by edge order):")
    for i, e in enumerate(edges):
        print(f"  [{col[i]}] {tuple(e)}")


def search(max_n=16, time_budget_sec=600):
    start = time.time()
    n = 12
    while n <= max_n:
        elapsed = time.time() - start
        if elapsed > time_budget_sec:
            print(f"\nTime budget {time_budget_sec}s exhausted at n = {n}.")
            return
        print(f"\n=== n = {n} ===  (elapsed {elapsed:.1f}s)")
        tcount = 0
        total_reductions = 0
        total_colorings = 0
        try:
            it = graphs.triangulations(n, minimum_degree=5)
        except Exception as ex:
            print(f"  cannot enumerate triangulations: {ex}")
            n += 1
            continue
        for G in it:
            tcount += 1
            if not G.is_planar(set_embedding=True):
                continue
            D = dual_of(G)
            if not D.is_planar(set_embedding=True):
                continue
            faces_D = D.faces()
            pentagonal = [f for f in faces_D if len(f) == 5]
            for face in pentagonal:
                for i_red in range(5):
                    H, named = apply_reduction(D, face, i_red)
                    if H is None:
                        continue
                    total_reductions += 1
                    edges, gen = proper_3_edge_colorings(H)
                    # Find indices of named edges
                    named_idx = {}
                    for ii, e in enumerate(edges):
                        es = frozenset((e[0], e[1]))
                        for role, ns in named.items():
                            if es == ns:
                                named_idx[role] = ii
                    if len(named_idx) != 4:
                        continue
                    for col in gen:
                        total_colorings += 1
                        if matches(edges, col, named_idx):
                            report(edges, col, named_idx, named, {
                                'n': n,
                                'triangulation_index': tcount,
                                'face_size': len(face),
                                'i_red': i_red,
                                'reduced_dual_V': H.order(),
                                'reduced_dual_E': H.size(),
                            })
                            return
            sys.stdout.flush()
        print(f"  n = {n}: {tcount} triangulation(s), "
              f"{total_reductions} reductions, "
              f"{total_colorings} colorings; no hit.")
        n += 1
    print(f"\nSearch exhausted up to n = {max_n}.")


if __name__ == '__main__':
    search()