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

"""Simpler version: enumerate the residual (|S|=8, hit=p_total=8)
colourings without requiring v_parent identification. For each,
determine:
  - the n_k sequence of the reduction (from the reduced dual's
    F_k structure),
  - whether any G'-face (length ≢ 0 mod 3) is uncovered.

The earlier check_S_face_structure.py showed at most 30 |S|=8 cases
with hit = 8, but didn't constrain p_total. Of those, ≤30 have
p_total = 8 (= the residual).

Run with: sage experiments/check_30_residual_v2.py
"""
import os
import sys
import time

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_3_8_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 is_g_prime_pentagon(f, named):
    if len(f) != 5: return False
    fset = {frozenset(e) for e in f}
    return not (named['side_0'] in fset or named['side_1'] in fset
                or named['spike'] in fset or named['merged'] in fset)


def is_g_prime_face(f, named):
    fset = {frozenset(e) for e in f}
    return not (named['side_0'] in fset or named['side_1'] in fset
                or named['spike'] in fset or named['merged'] in fset)


def test_one(D):
    D.is_planar(set_embedding=True)
    residual_colourings = []
    other_bad = []
    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
            # Compute the n_k sequence for this reduction
            # n_k = degree of the face F_k of original G' = length of
            # the corresponding face in the reduced dual after subtracting
            # the reduction effect.
            # Simpler: the flank face F^♭_{i, i+1} has length n_i - 1.
            # So n_i = length of F^♭_{i,i+1} + 1.
            # Get all faces:
            named_face_lengths = {}
            for f in H.faces():
                fset = {frozenset(e) for e in f}
                if named['side_0'] in fset and named['spike'] in fset:
                    named_face_lengths['flank_lower'] = len(f)
                if named['spike'] in fset and named['side_1'] in fset:
                    named_face_lengths['flank_upper'] = len(f)
                if (named['side_0'] in fset and named['side_1'] in fset
                    and named['merged'] in fset):
                    named_face_lengths['outer'] = len(f)
                if (named['merged'] in fset and named['side_0'] not in fset
                    and named['side_1'] not in fset
                    and named['spike'] not in fset):
                    named_face_lengths['merged'] = len(f)
            n_i = named_face_lengths.get('flank_lower', 0) + 1   # = n_i, where the flank covers
            n_ip1 = named_face_lengths.get('flank_upper', 0) + 1
            # n_{i+3} = F_merged length + 2
            n_ip3 = named_face_lengths.get('merged', 0) + 2
            # n_{i+2} + n_{i+4} from outer
            outer_len = named_face_lengths.get('outer', 0)
            # outer_len = n_{i+2} + n_{i+4} - 3

            for col in cand:
                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
                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
                if P_1 in V_union: continue
                # bad colouring
                # Count G'-pentagons total and hit
                p_total = 0
                p_hit = 0
                non_pent_uncovered = []
                for f in H.faces():
                    if not is_g_prime_face(f, named): continue
                    L = len(f)
                    verts = {u for (u, v) in f} | {v for (u, v) in f}
                    if L == 5:
                        p_total += 1
                        if verts & S: p_hit += 1
                    else:
                        if L % 3 != 0 and verts.issubset(V_union):
                            non_pent_uncovered.append(L)
                # Is this a "residual" case?
                S_size = len(S)
                if S_size == 8 and p_total == p_hit:
                    residual_colourings.append({
                        'n_i': n_i, 'n_ip1': n_ip1, 'n_ip3': n_ip3,
                        'outer_len': outer_len,
                        'p_total': p_total, 'p_hit': p_hit,
                        'S_size': S_size,
                        'non_pent_uncovered': non_pent_uncovered,
                    })
                elif S_size == 8 and p_hit == p_total - 1 and p_total >= 7:
                    # Border case: only 1 pentagon uncovered
                    pass

    return residual_colourings


def main(max_n=20, time_budget_per_n=1800):
    print("Detailed analysis of |S|=8, p_hit = p_total residual "
          "chord-apex+Kempe colourings.\n")
    grand_residuals = []
    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_count = 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}")
                break
            G.is_planar(set_embedding=True)
            D = dual_of(G)
            resids = test_one(D)
            for r in resids:
                r['n_G'] = n
                r['tri_idx'] = tri_idx
            n_count += len(resids)
            grand_residuals.extend(resids)
        elapsed = time.time() - start
        print(f"n={n}: {n_count} residual colourings [{elapsed:.0f}s]")
        sys.stdout.flush()
    print()
    print("=" * 70)
    print(f"Total residual colourings: {len(grand_residuals)}")
    if grand_residuals:
        # n_k sequence distribution
        seq_dist = {}
        for r in grand_residuals:
            key = (r['n_i'], r['n_ip1'], r['n_ip3'], r['outer_len'])
            seq_dist[key] = seq_dist.get(key, 0) + 1
        print("\n(n_i, n_{i+1}, n_{i+3}, F_outer length) distribution:")
        for k, c in sorted(seq_dist.items(), key=lambda x: -x[1]):
            print(f"  {k}: {c}")
        # For each, does a non-pentagon G'-face provide a deciding face?
        has_non_pent = sum(1 for r in grand_residuals if r['non_pent_uncovered'])
        print(f"\nResidual colourings with at least one length-≢0-mod-3 "
              f"G'-face uncovered: {has_non_pent} / {len(grand_residuals)} "
              f"({100*has_non_pent/len(grand_residuals):.2f}%)")
        # Lengths of those non-pent G'-faces
        len_dist = {}
        for r in grand_residuals:
            for L in r['non_pent_uncovered']:
                len_dist[L] = len_dist.get(L, 0) + 1
        print("Lengths of uncovered non-pentagon G'-faces:")
        for L, c in sorted(len_dist.items()):
            print(f"  |f| = {L}: {c}")


if __name__ == '__main__':
    main()