"""Constants of Part II of the paper.

1. The Gaspers-Sorkin (TALG 2017) linear program for cubic Max 2-CSP, constraints (7)-(24) of Section 3.2 of
   arXiv 1404.0753v2, with the separator constant c in constraint (12); solved numerically, and the closed-form
   weights w_s=(1+c)/D, w_s2=1/D, w_r=w_b=2c/D, w_c=c/D (D=1+4c) checked exactly for rational c.
2. The lifts of Scott-Sorkin Theorem 22: beta(alpha)=7/50+3alpha/10, beta_4(alpha)=1/8+3alpha/8 for 1/9<=alpha<=1/5.
3. Exponential-space exponents from the project's width bounds, and the bases r^e for r=2,3.
Run from the repository root:  .venv/bin/python experiments/algorithmic-consequences/constants.py
"""
import math
from fractions import Fraction as Fr
from scipy.optimize import linprog

C3 = 3 / (4 * math.pi) * math.acos(5 / 6)
F = 0.14015825233145217
F4 = 0.24577197475           # Arb upper bound
fails = []


def up(x, d=4):
    """x rounded up in the d-th decimal, as the manuscript prints bases."""
    return "%.*f" % (d, math.ceil(x * 10 ** d - 1e-9) / 10 ** d)


def constraints(c):
    """Rows (coefficients on ws, ws2, wr, wb, wc; constant) meaning coef.w + const <= 0."""
    return {
        8: ([0, 0, 0, -1, 1], 0), 9: ([-1, 1, 0, 0, 0], 0), 10: ([0, 0, -1, 0, 1], 0), 11: ([0, 0, -1, 1, -1], 0),
        12: ([c, 0, (1 - c) / 2 - 1, 0, 0], 0), 13: ([1, -1, 0, -1, 1], 0), 14: ([1, -1, -1, 0, 1], 0),
        15: ([1, -1, -1, 1, -1], 0), 16: ([-1, 0, 1, 0, 0], 0), 17: ([-1, 0, 0, 0, 1], 0), 18: ([-1, 0, 0, 1, -1], 0),
        19: ([-2, 1, -1, 0, 0], 1), 20: ([-1, 0, -1, -1, 1], 1), 21: ([-1, 0, -1, 0, -1], 1), 22: ([0, 0, -1, 1, 0], 0),
        23: ([-1, 0, -2, 1, -1], 1), 24: ([-1, 0, -2, 0, 1], 1)}


def lp(c):
    rows = constraints(c)
    A = [r[0] for r in rows.values()]; b = [-r[1] for r in rows.values()]
    res = linprog([0, 0, 1, 0, 0], A_ub=A, b_ub=b, bounds=[(0, None)] * 5, method="highs")
    return res.fun, res.x


def closed_form(c):
    D = 1 + 4 * c
    return [(1 + c) / D, 1 / D, 2 * c / D, 2 * c / D, c / D]


# 1. exact feasibility of the closed form for rational c in (0, 1]; LP optimum equals 2c/(1+4c)
for c in [Fr(1, 6), Fr(1398221, 10 ** 7), Fr(1, 100), Fr(1, 2), Fr(1)]:
    w = closed_form(c)
    for k, (coef, const) in constraints(c).items():
        v = sum(a * x for a, x in zip(coef, w)) + const
        if v > 0:
            fails.append((c, k, v))
    if min(w) < 0:
        fails.append((c, "nonneg"))
for c in (1 / 6, C3, 0.1, 0.3):
    v, x = lp(c)
    if abs(v - 2 * c / (1 + 4 * c)) > 1e-9:
        fails.append(("lp", c, v))
print("Gaspers-Sorkin LP: closed-form weights feasible for c in {1/6, 0.1398221, 1/100, 1/2, 1}; "
      "LP optimum = 2c/(1+4c) at c = 1/6, c3, 0.1, 0.3")
v, x = lp(C3)
print("  c = c3: w_r = %.6f (2c/(1+4c) = %.6f); weights ws, ws2, wr, wb, wc = %s" % (
    v, 2 * C3 / (1 + 4 * C3), [round(t, 6) for t in x]))

# 2. polynomial-space exponents
for name, c in (("published c = 1/6", 1 / 6), ("c3", C3)):
    wn = 2 * c / (1 + 4 * c); alpha = 2 * wn / 3
    beta = 7 / 50 + 3 * alpha / 10; beta4 = 1 / 8 + 3 * alpha / 8
    assert 1 / 9 <= alpha <= 1 / 5
    print("%-18s cubic r^(%.6f n) = r^(%.6f m) = r^(m/%.4f); Delta<=4 r^(%.6f m) = r^(m/%.4f); general r^(%.6f m) = r^(m/%.4f)"
          % (name, wn, alpha, 1 / alpha, beta4, 1 / beta4, beta, 1 / beta))
    print("%-18s bases (rounded up): Max-Cut cubic 2^(%.6f n) < %s^n; #DS subcubic 3^(%.6f n) < %s^n" % (
        "", wn, up(2 ** wn), wn, up(3 ** wn)))

# 3. exponential space
a = {3: 5 / 32, 4: 1 / 4, 5: 11 / 32, 6: 7 / 16}
print("bases below are rounded up in the last printed digit, as in the manuscript")
print("exp space cubic: 2^(c3 n) < %s^n (was 2^(n/6) < %s^n); 3^(c3 n) < %s^n (was %s^n)" % (
    up(2 ** C3), up(2 ** (1 / 6)), up(3 ** C3), up(3 ** (1 / 6))))
print("exp space Delta<=4: 2^(F4 n) < %s^n (was 2^(n/3) < %s^n); 3^(F4 n) < %s^n" % (up(2 ** F4), up(2 ** (1 / 3)), up(3 ** F4)))
print("exp space general: 2^(7m/48) < %s^m = 2^(m/%.4f) (was 2^(13m/75) < %s^m = 2^(m/%.4f))" % (
    up(2 ** (7 / 48)), 48 / 7, up(2 ** (13 / 75)), 75 / 13))
# Delta <= 4 in terms of m: max over x = n4/n of min(F4 n, (5 n3 + 8 n4)/32) / m
best = max((min(F4, (5 * (1 - x) + 8 * x) / 32) * 2 / (3 + x), x) for x in [i / 10 ** 5 for i in range(10 ** 5 + 1)])
xs = (32 * F4 - 5) / 3
print("exp space Delta<=4, delta>=3, in m: max_x min(F4 n, 5n3/32 + n4/4)/m = %.5f at n4/n = %.4f (closed form x = %.4f, value %.5f)"
      % (best[0], best[1], xs, 2 * F4 / (3 + xs)))
print("ALL CHECKS PASSED" if not fails else "FAILURES: %s" % fails)