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

"""For the 1,314 bad chord-apex+Kempe colourings (where
Lemma flank-covering-hex empirically fails), characterize the set
S = V \\ (V(K_b) ∪ V(K_c)) -- the "uncovered" vertices.

For each bad colouring, report:
  - |S| (size of S);
  - distribution of S-vertices' graph-distance from v_n (the spike
    vertex), and from the merged edge;
  - on which faces of the reduced dual the S-vertices sit (G'-pentagon,
    G'-hexagon, F^♭_flank, F^♭_outer, F^♭_merged, etc.);
  - per-colouring: how MANY G'-pentagons are hit by S (= how close to
    "covering all G'-pentagons" we are; the pigeonhole would need
    |S|*3 ≥ #pentagons).

The goal is to find structural patterns we could exploit in a
discharging argument.

Run with: sage experiments/characterize_S_vertices.py
"""
import os
import sys
import time
from collections import defaultdict

from sage.all import Graph
from sage.graphs.graph_generators import graphs

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, HERE)

from check_conj_final_scaled import (
    apply_reduction,
    proper_3_edge_colorings,
    matches_chord_apex_kempe,
    kempe_cycle_set,
    edge_idx,
)
from check_heawood_on_kempe import dual_of, vertices_of_kempe


def classify_face(f, named):
    """Classify a face of the reduced dual."""
    edges_set = {frozenset(e) for e in f}
    has_s0 = named['side_0'] in edges_set
    has_s1 = named['side_1'] in edges_set
    has_sp = named['spike'] in edges_set
    has_m = named['merged'] in edges_set
    if has_s0 and has_sp and not has_s1 and not has_m: return 'flank-lower'
    if has_sp and has_s1 and not has_s0 and not has_m: return 'flank-upper'
    if has_s0 and has_s1 and has_m and not has_sp: return 'outer'
    if has_m and not (has_s0 or has_s1 or has_sp): return 'merged'
    if not (has_s0 or has_s1 or has_sp or has_m): return 'G-prime'
    return 'mixed'


def test_one(D):
    D.is_planar(set_embedding=True)
    bad_count = 0
    S_size_dist = {}
    S_dist_from_vn = {}    # graph dist of S-vertex from v_n
    S_face_types = {}      # face types containing S-vertices
    pentagons_hit_dist = {} # how many G'-pentagons hit by S, per colouring
    for face in D.faces():
        if len(face) != 5: continue
        for i_red in range(5):
            res = apply_reduction(D, face, i_red, 9999)
            if res is None: continue
            H = res['H']; named = res['named']
            H.is_planar(set_embedding=True)
            edges, colorings = proper_3_edge_colorings(H)
            cand = [c for c in colorings
                    if matches_chord_apex_kempe(edges, c, named)]
            v_n = 9999
            for col in cand:
                # Identify P_1, P_2 on lower flank, check sub-case (ii.B)
                target = {named['side_0'], named['spike']}
                lower_flank = None
                for f in H.faces():
                    if target.issubset({frozenset(e) for e in f}):
                        lower_flank = f; break
                if lower_flank is None or len(lower_flank) != 5: continue
                arc_verts = [e[0] for e in lower_flank]
                if v_n not in arc_verts: continue
                k = arc_verts.index(v_n)
                cyc = arc_verts[k:] + arc_verts[:k]
                A_i = next(iter(named['side_0'] - {v_n}))
                A_ip1 = next(iter(named['spike'] - {v_n}))
                if cyc[1] == A_i and cyc[4] == A_ip1:
                    P_1, P_2 = cyc[2], cyc[3]
                elif cyc[1] == A_ip1 and cyc[4] == A_i:
                    P_2, P_1 = cyc[2], cyc[3]
                else: continue
                merged_idx = edge_idx(edges, named['merged'])
                c_col = col[merged_idx]
                c_0_col = col[edge_idx(edges, named['side_0'])]
                c_1_col = col[edge_idx(edges, named['side_1'])]
                e_AiP1 = edge_idx(edges, frozenset((A_i, P_1)))
                e_P1P2 = edge_idx(edges, frozenset((P_1, P_2)))
                if e_AiP1 is None or e_P1P2 is None: continue
                if col[e_AiP1] != c_1_col or col[e_P1P2] != c_0_col:
                    continue
                # Compute K_b, K_c, S
                a = c_col
                other = [x for x in range(3) if x != a]
                kc_b = kempe_cycle_set(edges, col, merged_idx, (a, other[0]))
                kc_c = kempe_cycle_set(edges, col, merged_idx, (a, other[1]))
                V_b = vertices_of_kempe(edges, kc_b)
                V_c = vertices_of_kempe(edges, kc_c)
                V_union = V_b | V_c
                S = set(H.vertices()) - V_union
                # P_1 ∉ V_union check
                if P_1 in V_union: continue
                bad_count += 1
                S_size = len(S)
                S_size_dist[S_size] = S_size_dist.get(S_size, 0) + 1
                # Distances from v_n
                for s in S:
                    try:
                        d = H.distance(s, v_n)
                    except:
                        d = -1
                    S_dist_from_vn[d] = S_dist_from_vn.get(d, 0) + 1
                # Face types containing S-vertices
                for s in S:
                    for f in H.faces():
                        verts = {u for (u, v) in f} | {v for (u, v) in f}
                        if s in verts:
                            t = classify_face(f, named)
                            S_face_types[t] = S_face_types.get(t, 0) + 1
                # # G'-pentagons hit by S
                n_pent_hit = 0
                n_pent_total = 0
                for f in H.faces():
                    if classify_face(f, named) != 'G-prime': continue
                    if len(f) != 5: continue
                    n_pent_total += 1
                    verts = {u for (u, v) in f} | {v for (u, v) in f}
                    if verts & S: n_pent_hit += 1
                pent_pair = (n_pent_hit, n_pent_total)
                pentagons_hit_dist[pent_pair] = (
                    pentagons_hit_dist.get(pent_pair, 0) + 1)
    return (bad_count, S_size_dist, S_dist_from_vn, S_face_types,
            pentagons_hit_dist)


def main(max_n=20, time_budget_per_n=1800):
    print("Structural characterization of S = V \\ (V(K_b) ∪ V(K_c))\n"
          "on the 1,314 bad chord-apex+Kempe colourings.\n")
    grand_bad = 0
    grand_S_size = {}
    grand_S_dist = {}
    grand_S_face = {}
    grand_pent_hit = {}
    for n in range(12, max_n + 1):
        start = time.time()
        try:
            triangulations = list(graphs.triangulations(n, minimum_degree=5))
        except Exception as ex:
            print(f"n={n}: cannot enumerate ({ex})")
            continue
        n_bad_n = 0
        for tri_idx, G in enumerate(triangulations):
            if time.time() - start > time_budget_per_n:
                print(f"  n={n}: timeout at tri {tri_idx}/{len(triangulations)}")
                break
            G.is_planar(set_embedding=True)
            D = dual_of(G)
            nb, ssize, sdist, sface, phit = test_one(D)
            n_bad_n += nb
            for k, v in ssize.items(): grand_S_size[k] = grand_S_size.get(k, 0) + v
            for k, v in sdist.items(): grand_S_dist[k] = grand_S_dist.get(k, 0) + v
            for k, v in sface.items(): grand_S_face[k] = grand_S_face.get(k, 0) + v
            for k, v in phit.items(): grand_pent_hit[k] = grand_pent_hit.get(k, 0) + v
        elapsed = time.time() - start
        print(f"n={n}: {n_bad_n} bad colourings [{elapsed:.0f}s]")
        sys.stdout.flush()
        grand_bad += n_bad_n
    print()
    print("=" * 70)
    print(f"Total bad colourings: {grand_bad}")
    print("\n|S| distribution:")
    for s in sorted(grand_S_size):
        c = grand_S_size[s]
        print(f"  |S| = {s}: {c} ({100*c/max(grand_bad, 1):.2f}%)")
    print("\nGraph distance of S-vertices from v_n:")
    for d in sorted(grand_S_dist):
        c = grand_S_dist[d]
        total_S_verts = sum(grand_S_dist.values())
        print(f"  dist = {d}: {c} ({100*c/max(total_S_verts, 1):.2f}%)")
    print("\nFace types containing S-vertices:")
    total = sum(grand_S_face.values())
    for t in sorted(grand_S_face, key=lambda k: -grand_S_face[k]):
        c = grand_S_face[t]
        print(f"  {t}: {c} ({100*c/max(total, 1):.2f}%)")
    print("\n(# G'-pentagons hit by S, # G'-pentagons total) per colouring:")
    for pair in sorted(grand_pent_hit, key=lambda k: -grand_pent_hit[k])[:20]:
        c = grand_pent_hit[pair]
        h, t = pair
        ratio = h / t if t > 0 else 0
        print(f"  hit/total = {h}/{t} ({100*ratio:.0f}%): {c}")


if __name__ == '__main__':
    main()