"""Repaired analysis of polynomial-space #IS for maximum degree 3 (uniform separator weight).

Checks for the revised subsection "Counting independent sets in graphs of maximum degree 3" of
the paper (in Appendix C.6 (`app:is3`); 2026-10-10 repair of the open points of the published analysis).

The repaired subroutine #3IS gives every separator vertex the same weight s3 (no spider weights). Its analysis
needs only the following constraints, where c is the separator constant and kappa = c s3 + (1-c) r3/2:
  (U1) s3 + 3 r3 / 2 >= 1        every branching of #3IS has the vector (1, 1) [lem:is3-credits]
  (U2) s3 > r3, 2 s3 >= 1        moves of a separator vertex to a side; separator credits in the zone
  (U3) kappa < r3                new separations [eq:gl-newsep with the single weight s3]
  (U4) w2' >= 0, 2 w3' + w2' >= 2 kappa, w3' >= w2'   hand-over at average degree 2 and 8/3
  (U5) (B11) and (B12): the (3,3,3) branching vectors outside #3IS, unchanged.
Parts:
  1. certification of explicit rational weights at c = 0.13985 (exact rationals, mpmath interval arithmetic);
  2. the optimum of the program in closed form (kappa -> 2c/(1+4c)), floating point, for c3, 0.13985, 1/6;
  3. numbers for the remark on the published analysis (weights of the earlier version of the manuscript);
  4. a randomized test (not a proof) of the counting claims of lem:is3-credits on random subcubic graphs.
Run (from the repository root):  .venv/bin/python experiments/algorithmic-consequences/is3_uniform_repair.py
"""
import math
import random
from collections import defaultdict
from fractions import Fraction as Fr
from mpmath import iv, mp, mpf, findroot

iv.dps = 40
mp.dps = 30
fails = []
C3 = 3 / (4 * math.pi) * math.acos(5 / 6)


def ivq(q):
    return iv.mpf(q.numerator) / iv.mpf(q.denominator)


def bsum(a, b):
    """Interval enclosure of 2^-a + 2^-b for rationals a, b."""
    return iv.mpf(2) ** (-ivq(Fr(a))) + iv.mpf(2) ** (-ivq(Fr(b)))


def need(label, cond, value=None):
    print("  %-62s %s%s" % (label, "ok" if cond else "FAIL", "" if value is None else "   (%s)" % value))
    if not cond:
        fails.append(label)


# ---------------------------------------------------------------- 1. certified weights
print("1. Explicit weights at c = 13985/100000 (exact rational arithmetic; vectors in interval arithmetic)")
c = Fr(13985, 100000)
need("c3 < c", Fr(C3) < c, "c3 = %.10f" % C3)
r3 = Fr("0.179364")
s3 = 1 - Fr(3, 2) * r3                        # = 0.730954, makes (U1) an equality
kappa = c * s3 + (1 - c) * r3 / 2
w3p = Fr("0.176755")
w2p = Fr("0.0052182")
print("  r3 = %s, s3 = %s, kappa = c s3 + (1-c) r3/2 = %.12f, w2' = %s, w3' = %s"
      % (r3, s3, float(kappa), w2p, w3p))
need("(U1) s3 + 3 r3/2 >= 1", s3 + Fr(3, 2) * r3 >= 1, float(s3 + Fr(3, 2) * r3))
s = bsum(s3 + Fr(3, 2) * r3, s3 + Fr(3, 2) * r3)
need("(U1) vector (1,1): 2^-a + 2^-b <= 1", s.b <= 1, "upper end %.10f" % float(s.b))
need("(U2) s3 > r3", s3 > r3)
need("(U2) 2 s3 >= 1 (hence s3 >= 3 r3/2 and s3 >= r3/2)", 2 * s3 >= 1)
need("(U3) kappa < r3 (strict new-separator constraint)", kappa < r3, "r3 - kappa = %.3e" % float(r3 - kappa))
need("(U4) w2' >= 0", w2p >= 0)
need("(U4) 2 w3' + w2' >= 2 kappa", 2 * w3p + w2p >= 2 * kappa, "slack %.3e" % float(2 * w3p + w2p - 2 * kappa))
need("(U4) w3' >= w2'", w3p >= w2p)
b11 = bsum(4 * w3p - 3 * w2p, 8 * w3p - 4 * w2p)
need("(B11) (4w3'-3w2', 8w3'-4w2')", b11.b <= 1, "sum <= %.9f" % float(b11.b))
b12 = bsum(6 * w3p - w2p, 6 * w3p - w2p)
need("(B12) (6w3'-w2', 6w3'-w2')", b12.b <= 1, "sum <= %.9f" % float(b12.b))
base = iv.mpf(2) ** ivq(w3p)
need("2^w3' < 1.1304", base.b < iv.mpf("1.1304"), "2^w3' <= %.8f" % float(base.b))
D, B = 40 * r3, 42 * r3
need("constants D = 40 r3, B = 42 r3: B > D + r3", B > D + r3, "B = %.4f" % float(B))
print("  N0 must exceed 16 (lem:is3-credits), 1/w3' = %.3f, 1/kappa = %.3f, 2/r3 = %.3f"
      % (1 / float(w3p), 1 / float(kappa), 2 / float(r3)))

# ---------------------------------------------------------------- 2. optimum
print("2. Optimum of the program (floating point; kappa can be taken arbitrarily close to 2c/(1+4c))")


def optimum(cc):
    k = mpf(2) * cc / (1 + 4 * mpf(cc))
    w = findroot(lambda w: mpf(2) ** (-(10 * w - 6 * k)) + mpf(2) ** (-(16 * w - 8 * k)) - 1, mpf("0.17"))
    return k, w, 2 * k - 2 * w


for name, cc in (("c3", mpf(C3)), ("0.13985", mpf("0.13985")), ("1/6", mpf(1) / 6)):
    k, w, w2 = optimum(cc)
    print("  c = %-8s r3 -> 2c/(1+4c) = %.9f, w3' = %.9f, w2' = %.9f, base 2^w3' = %.7f"
          % (name, float(k), float(w), float(w2), float(mpf(2) ** w)))
print("  (with B11 tight and 2w3' + w2' = 2 kappa; the binding constraints are (U1), (U3), (U4), (B11))")

# ---------------------------------------------------------------- 3. published analysis
print("3. The published local analysis, with the spider weights of the earlier version of the manuscript")
R3, S2, S3, SB = Fr("0.179364"), Fr("0.641272"), Fr("0.716333"), Fr("0.730954")
BB = 6 * S3
dlt = SB - S3
print("  weights r3 = %s, s2 = %s, s3 = %s, s3bar = %s, B = 6 s3, delta = s3bar - s3 = %s" % (R3, S2, S3, SB, dlt))
# (i) left spider at imbalance 47 r3 (crossing)
dec = SB + R3 + (BB - 47 * R3 / 2)
x = bsum(dec, dec)
need("(i) left-spider branching at imbalance 47 r3: sum > 1", x.a > 1, "sum >= %.5f" % float(x.a))
# K1 at the threshold: branching on an unweighted (2,2,2) separator vertex with two right neighbours, or the pull
x = bsum(S3 + Fr(3, 2) * R3, S3 + Fr(3, 2) * R3)
need("(vi) unweighted (2,2,2) vertex at imbalance 2B, branching: sum > 1", x.a > 1, "sum >= %.5f" % float(x.a))
print("      the pull of spider at imbalance 2B raises the measure by r3 = %s (and the logarithmic term)" % R3)
# (iv) drags of the center case
cost = R3 / 2 - (S3 - S2)
print("  (iv) one unweighted drag raises the measure by r3/2 - Delta s3 = %.6f" % float(cost))
for k in (1, 2):
    x = bsum(SB + Fr(3, 2) * R3 - k * cost, SB + Fr(3, 2) * R3 - k * cost)
    need("(iv) (B1) charged with %d drag(s): sum > 1" % k, x.a > 1, "sum >= %.5f" % float(x.a))
# (v) the other vertex of a center pair becomes a right spider
x = bsum(S3 + Fr(3, 2) * R3, S3 + Fr(3, 2) * R3)
need("(v) (B1) when the partner of the pair becomes a spider: sum > 1", x.a > 1, "sum >= %.5f" % float(x.a))

# ---------------------------------------------------------------- 4. randomized test of the credit lemma
print("4. Randomized test of the counting claims of lem:is3-credits (not a proof)")


def random_instance(n3, rng, pdirect):
    while True:
        stubs = [v for v in range(n3) for _ in range(3)]
        rng.shuffle(stubs)
        pairs = [(stubs[2 * i], stubs[2 * i + 1]) for i in range(len(stubs) // 2)]
        if all(a != b for a, b in pairs):
            break
    adj = defaultdict(set)
    nxt, seen = n3, set()
    for a, b in pairs:
        key = (min(a, b), max(a, b))
        k = 0 if (key not in seen and rng.random() < pdirect) else rng.choice([1, 1, 2])
        seen.add(key)
        prev = a
        for _ in range(k):
            adj[prev].add(nxt)
            adj[nxt].add(prev)
            prev, nxt = nxt, nxt + 1
        adj[prev].add(b)
        adj[b].add(prev)
    for v in range(n3):                       # no (3,3,3) vertex
        if all(len(adj[u]) == 3 for u in adj[v]):
            u = next(iter(adj[v]))
            adj[v].discard(u)
            adj[u].discard(v)
            for x_, y_ in ((v, nxt), (nxt, u)):
                adj[x_].add(y_)
                adj[y_].add(x_)
            nxt += 1
    return {v: set(t) for v, t in adj.items()}


def skel(adj, v):
    """(far end, inner vertices) of the skeleton edges at the degree-3 vertex v."""
    res = []
    for u in adj[v]:
        path, prev, cur = [], v, u
        while len(adj[cur]) == 2:
            path.append(cur)
            prev, cur = cur, [w for w in adj[cur] if w != prev][0]
        res.append((cur, path))
    return res


def random_separation(adj, rng):
    side = {v: rng.choice("LR") for v in adj if len(adj[v]) == 3}
    for v in list(side):
        for f, _ in skel(adj, v):
            if f != v and side[v] in "LR" and side[f] in "LR" and side[v] != side[f]:
                side[rng.choice([v, f])] = "S"
    for v in list(side):
        for f, path in skel(adj, v):
            ends = (side[v], side[f])
            sd = "L" if "L" in ends else ("R" if "R" in ends else rng.choice("LR"))
            for p in path:
                side.setdefault(p, sd)
    for v in adj:
        side.setdefault(v, rng.choice("LR"))
    return side


def valid(adj, side):
    return all({side[v], side[u]} != {"L", "R"} for v in adj for u in adj[v])


def counts(adj, side):
    return (sum(1 for v in adj if side[v] == "S"),
            sum(1 for v in adj if side[v] == "R" and len(adj[v]) == 3),
            sum(1 for v in adj if side[v] == "L" and len(adj[v]) == 3))


def delete(adj, side, X):
    for x_ in X:
        for y_ in adj[x_]:
            if y_ not in X:
                adj[y_].discard(x_)
    for x_ in X:
        del adj[x_], side[x_]


def fold(adj, side):
    changed = True
    while changed:
        changed = False
        for v in list(adj):
            if v in adj and len(adj[v]) <= 1:
                delete(adj, side, {v})
                changed = True


def normalize(adj, side):
    """Rule (N1): move the maximal path of a separator vertex of degree 2 to a side."""
    while True:
        bad = [v for v in adj if side[v] == "S" and len(adj[v]) <= 2]
        if not bad:
            return
        v = bad[0]
        inner, ends = [v], []
        for start in list(adj[v]):
            prev, cur = v, start
            while len(adj[cur]) == 2 and cur != v:
                inner.append(cur)
                prev, cur = cur, [w for w in adj[cur] if w != prev][0]
            if cur == v:
                ends = None
                break
            ends.append(cur)
        es = [side[e] for e in ends] if ends is not None else []
        to = "R" if ("R" in es and "L" not in es) else "L"
        for x_ in inner:
            side[x_] = to
        if "R" in es and "L" in es:
            side[ends[es.index("R")]] = "S"


def comp3(adj, v):
    seen, st = {v}, [v]
    while st:
        x_ = st.pop()
        for y_ in adj[x_]:
            if y_ not in seen:
                seen.add(y_)
                st.append(y_)
    return sum(1 for x_ in seen if len(adj[x_]) == 3)


tested, bad, patterns = 0, [], defaultdict(int)
for pd in (0.0, 0.3, 0.6, 0.85):
    for seed in range(250):
        rng = random.Random(20261010 + int(1000 * pd) * 1000 + seed)
        adj = random_instance(rng.choice([22, 26, 34, 50]), rng, pd)
        side = random_separation(adj, rng)
        assert valid(adj, side)
        if any(f == v for v in adj if len(adj[v]) == 3 for f, _ in skel(adj, v)):
            continue                          # a loop: the multiplier reduction applies
        for v in [x_ for x_ in adj if side[x_] == "S"]:
            if comp3(adj, v) <= 20:
                continue                      # small component: rule (S3)
            e = defaultdict(int)
            for f, _ in skel(adj, v):
                e[side[f]] += 1
            if e["L"] < 1 or e["R"] < 1:
                continue                      # rule (N2) applies
            patterns[(e["L"], e["R"], e["S"])] += 1
            for branch in (1, 2):
                a2, s2 = {x_: set(t) for x_, t in adj.items()}, dict(side)
                S0, R0, L0 = counts(a2, s2)
                delete(a2, s2, {v} if branch == 1 else {v} | set(a2[v]))
                fold(a2, s2)
                normalize(a2, s2)
                S1, R1, L1 = counts(a2, s2)
                dS, dR, dL = S0 - S1, R0 - R1, L0 - L1
                tested += 1
                ok = valid(a2, s2) and all(not (s2[x_] == "S" and len(a2[x_]) != 3) for x_ in a2)
                ok = ok and dS >= 1 and dR >= 0 and dL >= 0 and (dS - 1) + dR + dL >= 3
                if e["R"] >= e["L"]:
                    ok = ok and (dS >= 2 or (dR >= 2 and dR + dL >= 3))
                if not ok:
                    bad.append((pd, seed, v, branch, dict(e), dS, dR, dL))
print("  branches tested: %d (patterns e_L, e_R, e_S: %s); violations: %d"
      % (tested, dict(sorted(patterns.items())), len(bad)))
if bad:
    fails.append(("credit lemma, randomized", bad[:3]))

print("ALL CHECKS PASSED" if not fails else "FAILURES: %s" % fails)