# Independent exact verification of antinorm certificates: recompute (not reuse) the
# convex combinations by LP, then certify each exactly from an optimal basis in Fractions.
# Membership p in t*R  (R = conv(Z) + K', K' = cumulative order) <=> exists mu>=0, sum mu >= t,
# cum(p) - sum mu_i cum(z_i) >= 0.  We maximise sum mu (=ell(p)).
import json, sys, itertools
from fractions import Fraction as Fr
import numpy as np
from scipy.optimize import linprog

def merge(P,Q):
    A,B,C=P; a,b,c=Q
    return (A*(2*a+b+c), B*(a+b+c/2)+C*a, C*(b+c/2))
def cum(P): return (P[0], P[0]+P[1], P[0]+P[1]+P[2])

def ell_exact_lower(p, CZ):
    """return an exact rational lower bound on ell(p) from an LP basis (or None)."""
    cp = cum(p)
    A = np.array([[float(cz[d]) for cz in CZ] for d in range(3)])
    b = np.array([float(x) for x in cp])
    res = linprog(-np.ones(len(CZ)), A_ub=A, b_ub=b, bounds=[(0,None)]*len(CZ), method='highs')
    assert res.status == 0
    supp = [i for i in range(len(CZ)) if res.x[i] > 1e-12]
    best = None
    # try all subsets of size <=3 of support (and, if needed, of near-support) with tight rows
    cand = sorted(range(len(CZ)), key=lambda i: -res.x[i])[:6]
    for k in (1,2,3):
        for S in itertools.combinations(cand, k):
            for rows in itertools.combinations(range(3), k):
                M = [[CZ[i][d] for i in S] for d in rows]
                r = [cp[d] for d in rows]
                sol = solve(M, r)
                if sol is None or any(x < 0 for x in sol): continue
                ok = all(cp[d] - sum(sol[j]*CZ[S[j]][d] for j in range(k)) >= 0 for d in range(3))
                if ok:
                    v = sum(sol)
                    if best is None or v > best: best = v
    return best, -res.fun

def solve(M, r):
    k = len(M)
    M = [list(row)+[rv] for row, rv in zip(M, r)]
    for col in range(k):
        piv = next((i for i in range(col,k) if M[i][col] != 0), None)
        if piv is None: return None
        M[col], M[piv] = M[piv], M[col]
        for i in range(k):
            if i != col and M[i][col] != 0:
                f = M[i][col]/M[col][col]
                M[i] = [x - f*yv for x, yv in zip(M[i], M[col])]
    return [M[i][k]/M[i][i] for i in range(k)]

for f in sys.argv[1:]:
    c = json.load(open(f))
    y, e1, e2, m0 = Fr(*c['y']), Fr(*c['e1']), Fr(*c['e2']), c['m0']
    Z = [tuple(Fr(*v) for v in z) for z in c['Z']]
    assert all(min(z) >= 0 for z in Z) and 0 < e2 < 1 and e1 > 0
    CZ = [cum(z) for z in Z]
    worst = None
    for zi in Z:
        for zj in Z:
            lb, fl = ell_exact_lower(merge(zi, zj), CZ)
            assert lb is not None and lb >= 1, (lb, fl)
            worst = lb if worst is None else min(worst, lb)
    winit = None
    for m in list(range(m0+1)) + ['inf']:
        p = (y/2, Fr(0), Fr(0) if m == 'inf' else Fr(1, 2**m))
        need = e1 * e2**(m0+1 if m == 'inf' else m)
        lb, fl = ell_exact_lower(p, CZ)
        assert lb is not None and lb >= need, (m, float(lb), float(need))
        r = lb/need; winit = r if winit is None else min(winit, r)
    # also check that Init(m) for m in (m0, 200] is directly fine (sanity)
    cmin = min(z[0]+z[1]+z[2]/2 for z in Z)
    print(f"{f}: PRIMAL-EXACT VERIFIED; min ell(z_i o z_j) >= {float(worst):.8f}; min ell(Init)/needed = {float(winit):.8f}; c = {float(cmin):.6f}; |Z| = {len(Z)}")