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

"""Test the CW-order-forced parity relation at each shared vertex.

At each v in V(K_b) cap V(K_c), under the constancy hypothesis h(v) = +1
(CW colour order (a, b, c) at v), Lemma A predicts:

  s_b(v) = i_b(v) mod 2       # side of c-edge relative to K_b walk
  s_c(v) = (i_c(v) + 1) mod 2  # side of b-edge relative to K_c walk

The actual sides are determined by the embedding and the actual h(v).
The formulas combine to (regardless of h(v)):

  s_b(v) XOR s_c(v) = (i_b(v) XOR i_c(v) XOR 1) mod 2

which is a structural identity. We verify this identity and tally the
joint (h, i_b mod 2, i_c mod 2, s_b, s_c) over all shared vertices.

Most informative: under the constancy hypothesis (h(v) = +1 fixed), the
formulas s_b = i_b and s_c = 1 XOR i_c are the prediction. We check at
each shared v whether the actual (s_b, s_c) matches this prediction --
and where the mismatches occur. Mismatches signal Heawood h(v) = -1 at v.

Run with: sage experiments/check_cw_parity_prediction.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_final_scaled import (
    apply_reduction,
    proper_3_edge_colorings,
    matches_chord_apex_kempe,
    trace_kempe_cycle,
    edge_idx,
)
from check_heawood_on_kempe import dual_of, heawood_numbers
from check_heawood_local_side import c_edge_local_side


def walk_positions(walk):
    pos = {}
    for k, (_, leave_v) in enumerate(walk):
        if leave_v not in pos:
            pos[leave_v] = k
    return pos


def test_one(D):
    D.is_planar(set_embedding=True)
    n_col = 0
    # For each shared vertex: record (h_b, i_b mod 2, i_c mod 2, s_b, s_c).
    joint = {}
    # Verify the identity s_b XOR s_c == (i_b XOR i_c XOR 1) for each
    # shared v.
    identity_holds = 0; identity_fails = 0
    # Under constancy h = +1 prediction: predict s_b = i_b, s_c = 1 XOR i_c.
    # Count how often the prediction matches per-vertex.
    constancy_matches = 0
    constancy_total = 0
    # Per-colouring: count #shared v where constancy prediction matches.
    constancy_dist = {}     # # of shared v matching -> #colourings
    constancy_perfect = 0   # #colourings where prediction matches at all v
    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)
            emb = H.get_embedding()
            edges, colorings = proper_3_edge_colorings(H)
            cand = [c for c in colorings
                    if matches_chord_apex_kempe(edges, c, named)]
            for col in cand:
                n_col += 1
                try:
                    h = heawood_numbers(H, edges, col)
                except RuntimeError:
                    continue
                merged_idx = edge_idx(edges, named['merged'])
                a = col[merged_idx]
                bs = [c for c in range(3) if c != a]
                b_color, c_color = bs[0], bs[1]
                # The "third colour" at each cycle.
                walk_b = trace_kempe_cycle(edges, col, merged_idx, (a, b_color))
                walk_c = trace_kempe_cycle(edges, col, merged_idx, (a, c_color))
                pos_b = walk_positions(walk_b)
                pos_c = walk_positions(walk_c)
                V_b = set(pos_b.keys())
                V_c = set(pos_c.keys())
                shared = V_b & V_c
                # Compute s_b for each v on K_b
                Lb = len(walk_b)
                sides_b = {}
                for k in range(Lb):
                    v = walk_b[k][1]
                    in_e = edges[walk_b[k][0]]
                    out_e = edges[walk_b[(k+1) % Lb][0]]
                    u_in = in_e[0] if in_e[1] == v else in_e[1]
                    u_out = out_e[0] if out_e[1] == v else out_e[1]
                    s = c_edge_local_side(v, c_color, col, edges, emb,
                                          u_in, u_out)
                    sides_b[v] = (0 if s == 'L' else 1) if s else None
                # Compute s_c for each v on K_c. "Third colour off K_c" = b_color.
                Lc = len(walk_c)
                sides_c = {}
                for k in range(Lc):
                    v = walk_c[k][1]
                    in_e = edges[walk_c[k][0]]
                    out_e = edges[walk_c[(k+1) % Lc][0]]
                    u_in = in_e[0] if in_e[1] == v else in_e[1]
                    u_out = out_e[0] if out_e[1] == v else out_e[1]
                    s = c_edge_local_side(v, b_color, col, edges, emb,
                                          u_in, u_out)
                    sides_c[v] = (0 if s == 'L' else 1) if s else None
                # For each shared v, record.
                col_matches = 0
                col_shared_count = 0
                for v in shared:
                    pb = pos_b[v] % 2
                    pc = pos_c[v] % 2
                    sb = sides_b.get(v)
                    sc = sides_c.get(v)
                    hv = h[v]
                    h_b = 1 if hv == 1 else 0
                    if sb is None or sc is None: continue
                    key = (h_b, pb, pc, sb, sc)
                    joint[key] = joint.get(key, 0) + 1
                    # Identity check: s_b XOR s_c == (i_b XOR i_c XOR 1)
                    if (sb ^ sc) == (pb ^ pc ^ 1):
                        identity_holds += 1
                    else:
                        identity_fails += 1
                    # Constancy prediction: s_b = i_b, s_c = 1 XOR i_c
                    pred_match = (sb == pb) and (sc == (1 - pc))
                    constancy_total += 1
                    col_shared_count += 1
                    if pred_match:
                        constancy_matches += 1
                        col_matches += 1
                if col_shared_count > 0:
                    constancy_dist[(col_matches, col_shared_count)] = \
                        constancy_dist.get(
                            (col_matches, col_shared_count), 0) + 1
                    if col_matches == col_shared_count:
                        constancy_perfect += 1
    return (n_col, joint, identity_holds, identity_fails,
            constancy_matches, constancy_total,
            constancy_dist, constancy_perfect)


def main(max_n=18, time_budget_per_n=1800):
    print(f"CW-parity prediction at shared K_b cap K_c, n in [12, {max_n}]\n")
    grand_col = 0
    grand_joint = {}
    grand_id_ok = 0; grand_id_fail = 0
    grand_cm = 0; grand_ct = 0
    grand_cd = {}
    grand_cp = 0
    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_col_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)
            (ni, ji, ido, idf, cm, ct,
             cd, cp) = test_one(D)
            n_col_n += ni
            for k, v in ji.items(): grand_joint[k] = grand_joint.get(k, 0) + v
            grand_id_ok += ido; grand_id_fail += idf
            grand_cm += cm; grand_ct += ct
            for k, v in cd.items(): grand_cd[k] = grand_cd.get(k, 0) + v
            grand_cp += cp
        elapsed = time.time() - start
        print(f"n={n}: {n_col_n} col., [{elapsed:.0f}s]")
        sys.stdout.flush()
        grand_col += n_col_n

    print()
    print("=" * 78)
    print(f"Grand totals (n in [12, {max_n}], {grand_col} colourings)")
    print(f"\n  Structural identity s_b XOR s_c == i_b XOR i_c XOR 1:")
    print(f"    holds: {grand_id_ok}, fails: {grand_id_fail}")
    print(f"\n  Constancy prediction (s_b == i_b mod 2 AND s_c == 1 - i_c mod 2)")
    print(f"  match rate per shared vertex: "
          f"{grand_cm}/{grand_ct} ({100*grand_cm/max(1,grand_ct):.2f}%)")
    print(f"  perfect (all shared v match) colourings: "
          f"{grand_cp}/{grand_col} ({100*grand_cp/max(1,grand_col):.2f}%)")
    print(f"\n  Joint (h_b, i_b mod 2, i_c mod 2, s_b, s_c) distribution:")
    keys = sorted(grand_joint.keys())
    total = sum(grand_joint.values())
    # Group by (h_b, i_b mod 2, i_c mod 2) and show s_b, s_c spread:
    for k in keys:
        v = grand_joint[k]
        print(f"    {k}: {v} ({100*v/max(1,total):.2f}%)")


if __name__ == '__main__':
    main()