"""Gaspers-Lee compound analysis for maximum degree 4 and general graphs, rerun from the subcubic weights of
Theorem F.

Port of the authors' program deg5plus.py (2016; the degree-3 piece is injected from the subcubic analysis,
then five degree-4 sections with pivots 3, 3.2, 3.5, 3.75, 4, then degrees 5 and 6), with the constraints of
Algorithmica 85 (2023), Table 1 (Fig. 4) and Lemmas 7-8. Steps:
  1. reproduce the published chain from the published subcubic weights (w2', w3') = (0.02279113575, 0.18752030809);
  2. rerun from the certified weights of Theorem F, (w2', w3') = (0.004171, 0.176113), and certify explicit
     rational weights: linear conditions exactly, branching vectors in mpmath interval arithmetic, section running
     times by exact vertex enumeration of the degree-distribution polytope;
  3. the same for the exponential-space variant (their Theorem 2: subcubic measure 2^(n3/6) from the pathwidth
     dynamic program, reproducing 1.2330), and with 2^(c3 n3) from Lemma F1;
  4. exponential space with the pathwidth dynamic program for maximum degree 4 (Theorem 2(c)) and one measure for
     degrees 5-7; degree 7 is added to their analysis with out(v) >= 3 of the right parity (as in their Lemma 7).
Run from the repository root:  .venv/bin/python experiments/algorithmic-consequences/gaspers_lee_general.py
"""
from fractions import Fraction as Fr
from itertools import combinations_with_replacement as cwr
import numpy as np
from mpmath import iv
from scipy.optimize import minimize

iv.dps = 40
fails = []

# ------------------------------------------------------------------ constraint families (generic arithmetic)


def deg4_vectors(w2, w3, w4, psi):
    """Table 1 of the journal (Fig. 4), by section. psi is the potential of the (2,2,2,2) cases."""
    return [
        [(5 * w4 - 4 * w3 + 4 * w2 + psi,) * 2,            # (2,2,2,2)(a); the potential is removed
         (4 * w4 - 2 * w3 + 3 * w2 - psi,) * 2,            # (2,2,2,2)(b); the potential is added
         (4 * w4 - 2 * w3 + 2 * w2, 4 * w4 - 2 * w3 + 3 * w2),   # (2,2,2,3); journal prints "2 w1"
         (5 * w4 - 4 * w3 + 3 * w2, 6 * w4 - 4 * w3 + 3 * w2),   # (2,2,2,4)
         (3 * w4, 5 * w4 - 2 * w3 + 2 * w2),                     # (2,2,3,3)
         (4 * w4 - 2 * w3 + w2, 5 * w4 - 2 * w3 + 2 * w2),       # (2,2,3,4)
         (5 * w4 - 4 * w3 + 2 * w2, 7 * w4 - 4 * w3 + 2 * w2)],  # (2,2,4,4)
        [(2 * w4 + 2 * w3 - 2 * w2, 4 * w4 + w2)],               # (2,3,3,3)
        [(3 * w4 - w2, 6 * w4 - 2 * w3 + w2),                    # (2,3,3,4)
         (4 * w4 - 2 * w3, 6 * w4 - 2 * w3 + w2),                # (2,3,4,4)
         (5 * w4 - 4 * w3 + w2, 8 * w4 - 4 * w3 + w2),           # (2,4,4,4)
         (w4 + 4 * w3 - 4 * w2, 5 * w4),                         # (3,3,3,3)
         (2 * w4 + 2 * w3 - 3 * w2, 5 * w4)],                    # (3,3,3,4)
        [(3 * w4 - 2 * w2, 7 * w4 - 2 * w3),                     # (3,3,4,4)
         (4 * w4 - 2 * w3 - w2, 7 * w4 - 2 * w3)],               # (3,4,4,4)
        [(5 * w4 - 4 * w3, 9 * w4 - 4 * w3)]]                    # (4,4,4,4)


SECTIONS = [(2, 3), (3, Fr(16, 5)), (Fr(16, 5), Fr(7, 2)), (Fr(7, 2), Fr(15, 4)), (Fr(15, 4), 4)]


def out_edges(i):
    """Lemma 7: number of edges leaving N[v] for neighbor degrees i of v (degree 5 or 6)."""
    if len(i) == 5:
        o = 3 if sum(i) % 2 == 0 else 4
        return 5 if i == (2, 2, 2, 2, 2) else o
    o = 3 if sum(i) % 2 == 1 else 4
    if i == (2, 2, 2, 2, 2, 2):
        return 6
    return 5 if i == (2, 2, 2, 2, 2, 3) else o


def deg56_vectors(w):
    """Lemma 8 (as in deg5plus.py): w = [0, 0, w2, ..., w6]."""
    res = []
    for d, rng in ((5, [2, 3, 4, 5]), (6, [2, 3, 4, 5, 6])):
        for i in cwr(rng, d):
            if i == (d,) * d:
                continue                       # regular graph: happens once per component
            deg2 = sum(1 for k in i if k == 2)
            f1 = w[d] + sum(w[k] for k in i) + out_edges(i) * (w[d] - w[d - 1])
            f2 = w[d] + sum(w[k] - w[k - 1] for k in i) + deg2 * (w[d] - w[d - 1])
            res.append((f1, f2))
    return res


def deg56_linear(w, w4final):
    """Handover from degree 4 (weights at least the last degree-4 section), monotonicity and concavity."""
    return [w[2] - w4final[0], w[3] - w4final[1], w[4] - w4final[2],
            w[6] - w[5], w[5] - w[4], w[4] - w[3], w[3] - w[2],
            (w[5] - w[4]) - (w[6] - w[5]), (w[4] - w[3]) - (w[5] - w[4]), (w[3] - w[2]) - (w[4] - w[3]),
            w[2] - (w[3] - w[2])]


# ------------------------------------------------------------------ floating-point optimization


def solve(obj, ineqs, eqs, x0, starts=25):
    cons = [{"type": "ineq", "fun": lambda x: np.array(ineqs(x))}]
    if eqs:
        cons.append({"type": "eq", "fun": lambda x: np.array(eqs(x))})
    best = None
    rng = np.random.default_rng(0)
    for k in range(starts):
        xs = np.array(x0) * (1 + (0.1 * rng.standard_normal(len(x0)) if k else 0))
        r = minimize(obj, xs, constraints=cons, method="SLSQP", options={"ftol": 1e-14, "maxiter": 5000})
        if r.success and min(ineqs(r.x)) > -1e-9 and (not eqs or max(map(abs, eqs(r.x))) < 1e-8):
            if best is None or r.fun < best.fun:
                best = r
    return best.x


def br(v):
    return [1 - 2.0 ** -a - 2.0 ** -b for a, b in v]


def float_chain(W2, W3, margin=0.0):
    rows = []
    # section 1: x = (w4, psi)
    x = solve(lambda x: x[0], lambda x: br(deg4_vectors(W2, W3, x[0], x[1])[0]), None, [0.33, 0.0])
    w = [W2, W3, x[0]]
    rows.append((x[0], x[1], 0.0))
    for k in range(1, 5):
        p0 = float(SECTIONS[k][0])
        f = lambda a, w=w, p0=p0: [w[0] - (2 - p0) * a[0], w[1] - (3 - p0) * a[0], w[2] - (4 - p0) * a[0]]
        a = solve(lambda a: -a[0], lambda a, f=f, k=k: br(deg4_vectors(*f(a), 0)[k]), None, [0.01])
        w = f(a)
        rows.append((w[2], 0.0, a[0]))
    wf = w
    x = solve(lambda x: x[4], lambda x: [g - margin for g in br(deg56_vectors([0, 0, *x]))]
              + [g - margin for g in deg56_linear([0, 0, *x], wf)], None, [wf[0], wf[1], wf[2], 0.3, 0.31])
    return rows, wf, x


def section_time(w, p0, p1):
    """max sum a_i w_i over distributions on degrees {2,3,4} with average degree in [p0, p1] (exact)."""
    best = None
    degs = [2, 3, 4]
    cands = [{d: Fr(1)} for d in degs if p0 <= d <= p1]
    for i in degs:
        for j in degs:
            if i < j:
                for p in (Fr(p0), Fr(p1)):
                    if i <= p <= j:
                        t = (p - i) / (j - i)
                        cands.append({i: 1 - t, j: t})
    for c in cands:
        v = sum(c[d] * w[d - 2] for d in c)
        best = v if best is None or v > best else best
    return best


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


def check_vectors(vs, label):
    worst = 0
    for a, b in vs:
        s = iv.mpf(2) ** (-ivq(a)) + iv.mpf(2) ** (-ivq(b))
        worst = max(worst, float(s.b))
        if not s.b <= 1:
            fails.append((label, float(a), float(b)))
    return worst


def q(x, digits=7, up=None):
    s = 10 ** digits
    v = Fr(round(x * s), s)
    if up is True:
        v += Fr(1, s)
    if up is False:
        v -= Fr(1, s)
    return v


# ------------------------------------------------------------------ 1. reproduce the published chain
print("1. Published subcubic weights (0.02279113575, 0.18752030809), floating point:")
rows, wf, x = float_chain(0.02279113575, 0.18752030809)
w = [0.02279113575, 0.18752030809]
for (p0, p1), (w4, psi, a) in zip(SECTIONS, rows):
    print("   section [%.2f, %.2f]: w4 = %.7f, alpha = %.5f" % (float(p0), float(p1), w4, a))
print("   degree 4: final w2, w3, w4 = %s (published Fig. 5: 0.1057321, 0.1998925, 0.2713302)" % np.round(wf, 7))
print("   degrees 5-6: w2..w6 = %s; 2^w6 = %.5f (published 1.2356)" % (np.round(x, 7), 2 ** x[4]))

# ------------------------------------------------------------------ 2.-4. rerun and certify


def certify(W2, W3, label, chain_rows=None):
    rows, wf, x = float_chain(float(W2), float(W3))
    print("%s, floating point: degree-4 final %s, w6 = %.7f, 2^w6 = %.5f" % (label, np.round(wf, 7), x[4], 2 ** x[4]))
    _, _, x = float_chain(float(W2), float(W3), margin=1e-6)       # degrees 5-6 with slack, for rounding
    w4, psi = q(rows[0][0], up=True), q(rows[0][1])
    if psi < 0:
        fails.append((label, "psi < 0"))
    w = [W2, W3, w4]
    secw = [w]
    worst = [check_vectors(deg4_vectors(W2, W3, w4, psi)[0], label + " section 1")]
    times = [section_time(w, *SECTIONS[0])]
    for k in range(1, 5):
        p0 = Fr(SECTIONS[k][0])
        alpha = q(rows[k][2], up=False)
        if alpha < 0:
            fails.append((label, "alpha < 0", k))
        w = [w[0] - (2 - p0) * alpha, w[1] - (3 - p0) * alpha, w[2] - (4 - p0) * alpha]   # pivot shift, exact
        secw.append(w)
        worst.append(check_vectors(deg4_vectors(*w, 0)[k], "%s section %d" % (label, k + 1)))
        times.append(section_time(w, *SECTIONS[k]))
    # A case with highest average degree h can occur whenever d(G) <= h, so the weights of every section must also
    # satisfy the vectors of all later sections (the measure is the minimum over sections).
    cross = 0
    for k in range(5):
        for j in range(k + 1, 5):
            cross = max(cross, check_vectors(deg4_vectors(*secw[k], 0)[j], "%s cross %d-%d" % (label, k + 1, j + 1)))
    t4 = max(times)
    W = [0, 0] + [q(v, up=True) for v in x[:4]] + [q(x[4], up=True)]
    W[2], W[3], W[4] = max(W[2], w[0]), max(W[3], w[1]), max(W[4], w[2])
    for i, g in enumerate(deg56_linear(W, w)):
        if g < 0 or (i >= 3 and g == 0):
            fails.append((label, "degree 5-6 linear", i, float(g)))
    worst.append(check_vectors(deg56_vectors(W), label + " degrees 5-6"))
    base4 = float((iv.mpf(2) ** ivq(t4)).b)
    base = (iv.mpf(2) ** ivq(W[6]))
    if not TAU18 < float(base.a):
        fails.append((label, "tau(1,8) larger"))
    print("   certified: section weights (w2, w3, w4) ending at %s; largest branching sum %.7f; later-section cases"
          " under earlier-section weights: largest sum %.4f" % ([round(float(v), 7) for v in w], max(worst), cross))
    print("   maximum degree 4: 2^(%.7f n) < %.5f^n; general graphs: w2..w6 = %s, 2^(%s n) < %.5f^n"
          % (float(t4), base4, [round(float(v), 7) for v in W[2:]], W[6], float(base.b)))
    return base4, float(base.b)


TAU18 = 1.2320547            # upper bound for the root 1.23205463 of x^8 = x^7 + 1 (degree >= 7, their Theorem 1)
print("2. Polynomial space, from the weights (w2', w3') = (0.004171, 0.176113) of Theorem F")
b4, b = certify(Fr("0.004171"), Fr("0.176113"), "   poly. space")
print("   polynomial-space coloring via Bjorklund-Husfeldt-Koivisto: O((1 + %.5f)^n)" % b)
print("3. Exponential space (their Theorem 2): subcubic measure (0, 1/6), reproducing 1.2330, then (0, c3)")
certify(Fr(0), Fr(1, 6), "   exp. space, 1/6")
certify(Fr(0), Fr(1398222, 10 ** 7), "   exp. space, c3 <= 0.1398222")


def deg7_vectors(w):
    """Degree-7 branching (extension of Lemmas 7-8): out(v) >= 3 with the parity of sum_u d(u) - 7."""
    res = []
    for i in cwr([2, 3, 4, 5, 6, 7], 7):
        if i == (7,) * 7:
            continue
        deg2 = sum(1 for k in i if k == 2)
        o = 3 if (sum(i) - 7) % 2 == 1 else 4
        res.append((w[7] + sum(w[k] for k in i) + o * (w[7] - w[6]),
                    w[7] + sum(w[k] - w[k - 1] for k in i) + deg2 * (w[7] - w[6])))
    return res


print("4. Exponential space, maximum degree 4 by the pathwidth dynamic program with Theorem 2(c),"
      " (w2, w3, w4) >= (0, 5/32, 1/4), and degrees 5-7 in one measure")
A = [0, Fr(5, 32), Fr(1, 4)]


def lin7(w):
    return deg56_linear(w[:7], A) + [w[7] - w[6], (w[6] - w[5]) - (w[7] - w[6])]


x = solve(lambda x: x[5], lambda x: [g - 1e-6 for g in br(deg56_vectors([0, 0, *x[:5]])) + br(deg7_vectors([0, 0, *x]))]
          + [g - 1e-6 for g in lin7([0, 0, *x])], None, [0.11, 0.19, 0.26, 0.29, 0.30, 0.31])
W = [0, 0] + [q(v, up=True) for v in x]
for i, g in enumerate(lin7(W)):
    if g < 0 or (3 <= i and g == 0 and i != 11):
        fails.append(("exp. space degree 7 linear", i, float(g)))
worst = max(check_vectors(deg56_vectors(W[:7]), "exp. space degrees 5-6"), check_vectors(deg7_vectors(W), "exp. space degree 7"))
base = iv.mpf(2) ** ivq(W[7])
TAU19 = 1.2132                                        # root of x^9 = x^8 + 1 is 1.21315, for degree >= 8
print("   w2..w7 = %s; largest branching sum %.7f; general graphs 2^(%s n) < %.5f^n; degree >= 8: tau(1,9) < %.4f"
      % ([round(float(v), 7) for v in W[2:]], worst, W[7], float(base.b), TAU19))
if not TAU19 < float(base.a):
    fails.append("tau(1,9)")
print("ALL CHECKS PASSED" if not fails else "FAILURES: %s" % fails)