"""Check what fraction of V(Ĝ'_{v,i}) is covered by V(K_b) ∪ V(K_c)
across all chord-apex+Kempe colourings of all reduced duals up to
|V(G)| ≤ 18.

This is the key prerequisite for the Heawood-face-sum proof attempt:
if V(K_b) ∪ V(K_c) = V(reduced dual) in many cases, then under
constancy on V(K_b) ∪ V(K_c) every face of the reduced dual must
satisfy ε * |f| ≡ 0 (mod 3); since ε = ±1, this forces |f| ≡ 0 mod 3
for every face -- and any reduced dual with a face of length 4, 5, or
7 (etc., not multiple of 3) would yield an immediate contradiction.

Run with: sage experiments/check_kb_kc_coverage.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, heawood_numbers, vertices_of_kempe,
)


def test_one(D):
    D.is_planar(set_embedding=True)
    n_col = 0
    coverage_dist = {}  # |V(K_b)∪V(K_c)| / |V| as ratio (rounded) -> count
    full_coverage = 0
    face_lengths_at_full = {}  # face length distribution when full coverage
    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)]
            n_v = H.order()
            face_lens = sorted([len(f) for f in H.faces()])
            for col in cand:
                n_col += 1
                merged_idx = edge_idx(edges, named['merged'])
                a = col[merged_idx]
                bs = [c for c in range(3) if c != a]
                kc_b = kempe_cycle_set(edges, col, merged_idx, (a, bs[0]))
                kc_c = kempe_cycle_set(edges, col, merged_idx, (a, bs[1]))
                V_b = vertices_of_kempe(edges, kc_b)
                V_c = vertices_of_kempe(edges, kc_c)
                cov = len(V_b | V_c)
                cov_ratio = cov / n_v
                bucket = round(cov_ratio * 20) / 20   # bucket in 5% increments
                coverage_dist[bucket] = coverage_dist.get(bucket, 0) + 1
                if cov == n_v:
                    full_coverage += 1
                    key = tuple(face_lens)
                    face_lengths_at_full[key] = (
                        face_lengths_at_full.get(key, 0) + 1)
    return n_col, coverage_dist, full_coverage, face_lengths_at_full


def main(max_n=18, time_budget_per_n=300):
    print(f"V(K_b) ∪ V(K_c) coverage / V(Ĝ'_{{v,i}}) "
          f"on chord-apex+Kempe colourings, n_G ∈ [12, {max_n}]\n")
    grand_col = 0
    grand_dist = {}
    grand_full = 0
    grand_full_faces = {}
    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
        n_full_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, dist, full, full_faces = test_one(D)
            n_col_n += ni
            n_full_n += full
            for k, v in dist.items():
                grand_dist[k] = grand_dist.get(k, 0) + v
            for k, v in full_faces.items():
                grand_full_faces[k] = grand_full_faces.get(k, 0) + v
        elapsed = time.time() - start
        print(f"n={n}: {n_col_n} col., {n_full_n} full-coverage "
              f"({100 * n_full_n / max(n_col_n, 1):.1f}%) [{elapsed:.0f}s]")
        sys.stdout.flush()
        grand_col += n_col_n
        grand_full += n_full_n

    print()
    print("=" * 70)
    print(f"Grand totals: {grand_col} chord-apex+Kempe colourings")
    print(f"  full coverage (V(K_b)∪V(K_c) = V):  {grand_full} "
          f"({100 * grand_full / max(grand_col, 1):.2f}%)")
    print()
    print("  Coverage ratio distribution "
          "(|V(K_b)∪V(K_c)| / |V(Ĝ')|):")
    for r in sorted(grand_dist):
        c = grand_dist[r]
        pct = 100 * c / max(grand_col, 1)
        bar = '#' * max(1, int(pct))
        print(f"    {r:5.2f}: {c:>7} ({pct:5.2f}%) {bar[:50]}")
    if grand_full > 0:
        print()
        print("  Face-length distributions among full-coverage cases:")
        for key in sorted(grand_full_faces, key=lambda k: -grand_full_faces[k]):
            print(f"    {key}: {grand_full_faces[key]}")


if __name__ == '__main__':
    main()