#!/usr/bin/env python3
"""Reproduce the local nonexistence obstruction implying C(19,5,3) >= 104.

Python >= 3.10; standard library only. No input database, floating point,
optimization package, time cutoff, node cutoff, or symmetry assumption.

IMPORTANT PRIORITY NOTE: the local obstruction was independently computed
here, then found already proved more generally in Kovar and Zhang,
arXiv:2507.06745v2 (30 August 2026), Theorem 2.5. See proof.md.

Usage: python prove_c19.py --out rerun
"""
from __future__ import annotations
import argparse
import json
from collections import Counter, defaultdict
from itertools import combinations
from pathlib import Path

Adj = tuple[int, ...]
Pair = tuple[int, int]
MATCHINGS = (
    ((13, 14), (15, 16)),
    ((13, 15), (14, 16)),
    ((13, 16), (14, 15)),
)


def validate_graph(adj: Adj) -> None:
    n = len(adj)
    assert all(0 <= a < 1 << n for a in adj)
    assert all(not (adj[i] >> i & 1) for i in range(n))
    assert all(adj[i].bit_count() <= 3 for i in range(n))
    assert all((adj[i] >> j & 1) == (adj[j] >> i & 1)
               for i in range(n) for j in range(n))


def canonical(adj: Adj) -> tuple[int, Adj]:
    """Exact individualization/refinement, retaining an actual adjacency code.

    Equal returned codes at a fixed order always imply isomorphism, regardless
    of search-order details: each code encodes every edge of a relabelled graph.
    Twin pruning is only a transposition fixing every other vertex.
    """
    n = len(adj)

    def refine(partition):
        while True:
            new = []
            masks = [sum(1 << x for x in cell) for cell in partition]
            for cell in partition:
                buckets = defaultdict(list)
                for x in cell:
                    sig = tuple((adj[x] & mask).bit_count() for mask in masks)
                    buckets[sig].append(x)
                new.extend(tuple(buckets[s]) for s in sorted(buckets))
            if len(new) == len(partition):
                return new
            partition = new

    def visit(partition):
        partition = refine(partition)
        pos = next((i for i, c in enumerate(partition) if len(c) > 1), None)
        if pos is None:
            order = [c[0] for c in partition]
            code = 0
            for j in range(n):
                for k in range(j):
                    code = 2 * code + ((adj[order[j]] >> order[k]) & 1)
            return code, order
        cell = partition[pos]
        best = None
        tried = []
        for x in cell:
            if any((adj[x] & ~(1 << z)) == (adj[z] & ~(1 << x))
                   for z in tried):
                continue
            tried.append(x)
            answer = visit(partition[:pos] + [(x,), tuple(z for z in cell if z != x)]
                           + partition[pos + 1:])
            if best is None or answer[0] < best[0]:
                best = answer
        assert best is not None
        return best

    degrees = defaultdict(list)
    for x in range(n):
        degrees[adj[x].bit_count()].append(x)
    code, order = visit([tuple(degrees[d]) for d in sorted(degrees)])
    relabelled = tuple(sum(((adj[x] >> order[k]) & 1) << k for k in range(n))
                       for x in order)
    return code, relabelled


def generate_graphs() -> tuple[list[Adj], list[int]]:
    """All simple graphs on 10 vertices with 13 edges and maximum degree <=3.

    Every graph has a vertex deletion predecessor. All 0..3 neighbor sets
    among vertices of degree <3 are tried. An extension is pruned only if its
    edges exceed 13, or even 3 edges per remaining vertex cannot reach 13.
    """
    level = {0: ()}
    counts = []
    for n in range(1, 11):
        new = {}
        for adj in level.values():
            edges = sum(a.bit_count() for a in adj) // 2
            available = [i for i, a in enumerate(adj) if a.bit_count() < 3]
            for d in range(4):
                updated = edges + d
                if updated > 13 or updated + 3 * (10 - n) < 13:
                    continue
                for neighbors in combinations(available, d):
                    extended = list(adj) + [sum(1 << i for i in neighbors)]
                    for i in neighbors:
                        extended[i] |= 1 << (n - 1)
                    code, graph = canonical(tuple(extended))
                    new[code] = graph
        level = new
        counts.append(len(level))
        print(f"graph_level={n} representatives={len(level)}", flush=True)
    graphs = list(level.values())
    for adj in graphs:
        validate_graph(adj)
        assert sum(a.bit_count() for a in adj) == 26
    return graphs, counts


def reconstruct(adj: Adj, matching: int) -> dict:
    """A dual edge is an ordinary point of degree 2; a stub has degree 1."""
    edges = [(i, j) for i in range(10) for j in range(i + 1, 10)
             if adj[i] >> j & 1]
    assert len(edges) == 13
    stubs = [i for i in range(10) for _ in range(3 - adj[i].bit_count())]
    assert len(stubs) == 4
    triangles = [tuple([j for j, edge in enumerate(edges) if i in edge]
                       + [13 + j for j, owner in enumerate(stubs) if owner == i])
                 for i in range(10)]
    assert all(len(t) == 3 for t in triangles)
    shadow = Counter(p for t in triangles for p in combinations(sorted(t), 2))
    assert len(shadow) == 30 and all(d == 1 for d in shadow.values())
    matching_edges = set(MATCHINGS[matching])
    target = {p: 1 - shadow[p] + int(p in matching_edges)
              for p in combinations(range(17), 2)}
    assert all(0 <= d <= 2 for d in target.values())
    assert sum(target.values()) == 108
    return make_instance(target)


def make_instance(target: dict[Pair, int]) -> dict:
    """Construct every admissible 4-subset, not a restricted block family."""
    pairs = sorted(p for p, demand in target.items() if demand)
    index = {p: i for i, p in enumerate(pairs)}
    blocks, masks = [], []
    for quad in combinations(range(17), 4):
        six = tuple(combinations(quad, 2))
        if all(target.get(p, 0) > 0 for p in six):
            blocks.append(quad)
            masks.append(sum(1 << index[p] for p in six))
    edge_candidates = [0] * len(pairs)
    for q, mask in enumerate(masks):
        while mask:
            bit = mask & -mask
            mask -= bit
            edge_candidates[bit.bit_length() - 1] |= 1 << q
    demands = [target[p] for p in pairs]
    assert all(d in (1, 2) for d in demands)
    assert demands.count(2) <= 2
    return dict(blocks=blocks, masks=masks, edge_candidates=edge_candidates,
                demands=demands, target=target)


def exact_cover(inst: dict) -> tuple[list[int] | None, int]:
    """Complete finite search; None is returned only after exhausting branches.

    R records positive residual pairs and D the pairs with residual demand 2.
    Each chosen quad has six different pairs and D contains at most two pairs,
    so every selected quad immediately becomes inadmissible. Consequently the
    active candidates depend only on R, and memoization on (R,D) is sound.
    """
    masks, ec = inst['masks'], inst['edge_candidates']
    R = (1 << len(ec)) - 1
    D = sum(1 << i for i, d in enumerate(inst['demands']) if d == 2)
    failed = set()
    nodes = 0

    def visit(R: int, D: int, active: int):
        nonlocal nodes
        nodes += 1
        if not R:
            return []
        key = (R, D)
        if key in failed:
            return None
        remaining, best, best_count = R, 0, len(masks) + 1
        while remaining:
            bit = remaining & -remaining
            remaining -= bit
            choices = ec[bit.bit_length() - 1] & active
            count = choices.bit_count()
            if count < 1 + bool(D & bit):
                failed.add(key)
                return None
            if count < best_count:
                best_count, best = count, choices
        while best:
            choice = best & -best
            best -= choice
            q = choice.bit_length() - 1
            mask = masks[q]
            saturated, banned = mask & ~D, 0
            while saturated:
                bit = saturated & -saturated
                saturated -= bit
                banned |= ec[bit.bit_length() - 1]
            answer = visit((R & ~mask) | (D & mask), D & ~mask, active & ~banned)
            if answer is not None:
                return [q] + answer
        failed.add(key)
        return None

    return visit(R, D, (1 << len(masks)) - 1), nodes


def check_solution(inst: dict, solution: list[int]) -> None:
    assert len(set(solution)) == len(solution)
    counts = Counter(p for q in solution for p in combinations(inst['blocks'][q], 2))
    assert all(counts[p] == d for p, d in inst['target'].items())
    assert all(p in inst['target'] for p in counts)


def self_test() -> None:
    """Positive controls include all possible sizes 0,1,2 of the doubled set."""
    fixtures = [[], [(0, 1, 2, 3)],
                [(0, 1, 2, 3), (0, 1, 4, 5)],
                [(0, 1, 2, 3), (0, 1, 4, 5), (2, 3, 6, 7)]]
    # An affine plane of order 4 supplies 20 pair-disjoint quads.
    def multiply(a, b):
        answer = 0
        while b:
            if b & 1:
                answer ^= a
            b >>= 1
            a <<= 1
            if a & 4:
                a ^= 7
        return answer
    lines = [tuple(sorted(4*x + (multiply(m, x) ^ b) for x in range(4)))
             for m in range(4) for b in range(4)]
    lines += [tuple(4*x+y for y in range(4)) for x in range(4)]
    for shift in range(20):
        fixtures.append([line for i, line in enumerate(lines)
                         if i not in (shift, (shift+1) % 20)])
    for quads in fixtures:
        target = dict(Counter(p for q in quads for p in combinations(q, 2)))
        inst = make_instance(target)
        solution, _ = exact_cover(inst)
        assert solution is not None
        check_solution(inst, solution)
    # Two impossible targets, including a doubled demand with one candidate.
    for target in [{(0, 1): 1},
                   {p: (2 if p == (0, 1) else 1) for p in combinations(range(4), 2)}]:
        solution, _ = exact_cover(make_instance(target))
        assert solution is None
    # Every possible extension is generated; canonicalization only discards
    # equal actual edge encodings. Relabelling tests provide an extra guard.
    from random import Random
    rng = Random(20260906)
    for _ in range(50):
        n = 8
        adj = [0]*n
        possible = list(combinations(range(n), 2))
        rng.shuffle(possible)
        for i, j in possible:
            if adj[i].bit_count() < 3 and adj[j].bit_count() < 3 and rng.randrange(2):
                adj[i] |= 1 << j
                adj[j] |= 1 << i
        order = list(range(n))
        rng.shuffle(order)
        other = tuple(sum(((adj[order[i]] >> order[j]) & 1) << j for j in range(n))
                      for i in range(n))
        assert canonical(tuple(adj))[0] == canonical(other)[0]
    print('self_tests=PASS (24 positive, 2 negative, 50 graph relabellings)', flush=True)


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--out', type=Path, default=Path('results'))
    args = parser.parse_args()
    args.out.mkdir(parents=True, exist_ok=True)
    self_test()
    graphs, level_counts = generate_graphs()
    # Save actual representatives for the separately implemented C++ checker.
    (args.out / 'graphs.json').write_text(json.dumps(graphs, indent=2)+'\n')
    (args.out / 'graphs.txt').write_text(str(len(graphs))+'\n'+''.join(
        ' '.join(map(str, adj))+'\n' for adj in graphs))
    results = []
    for g, adj in enumerate(graphs):
        for matching in range(3):
            inst = reconstruct(adj, matching)
            solution, nodes = exact_cover(inst)
            row = dict(graph=g, matching=matching,
                       status='SAT' if solution is not None else 'UNSAT',
                       nodes=nodes, candidates=len(inst['blocks']))
            results.append(row)
            print(f"{g} {matching} {row['status']} {nodes}", flush=True)
            if solution is not None:
                check_solution(inst, solution)
                row['selected'] = [inst['blocks'][q] for q in solution]
                (args.out/'counterexample.json').write_text(json.dumps(row, indent=2))
                raise RuntimeError('A feasible local instance exists; no proof obtained.')
    degrees = Counter(tuple(sorted(a.bit_count() for a in adj)) for adj in graphs)
    report = dict(graph_level_counts=level_counts, graphs=len(graphs),
                  instances=len(results), satisfiable=0,
                  total_nodes=sum(row['nodes'] for row in results),
                  maximum_nodes=max(row['nodes'] for row in results),
                  degree_sequences={','.join(map(str, k)): v for k, v in sorted(degrees.items())},
                  results=results)
    (args.out / 'python_results.json').write_text(json.dumps(report, indent=2)+'\n')
    print('COMPLETE '+json.dumps({k:v for k,v in report.items() if k!='results'}), flush=True)
    assert level_counts == [1, 2, 4, 11, 23, 61, 141, 344, 598, 441]
    assert len(results) == 1323
    print('PROVED local nonexistence; see proof.md for C(19,5,3) >= 104.', flush=True)


if __name__ == '__main__':
    main()