"""Build the rational antinorm certificates for Directed Steiner Tree (lem:fa-dst-cert, prop:fa-dst-certs).

This program only *produces* candidate data; nothing in the paper depends on its correctness.
The certificates are checked in exact rational arithmetic by dst_certificate_verify.py.

A certificate is (y, e1, e2, Z) with y >= 0, e1, e2 > 0 rational and Z a finite set of nonnegative
rational points of the plane.  With R = conv(Z) + R_{>=0}^2 and the bilinear product
    (a, b) o (a', b') = (a (a' + b'), b (a'/2 + b')),
the hypotheses of the certificate lemma are
    (C1) z o z' in R for all z, z' in Z,
    (C2) (1, 0) in e2 R,
    (C3) (y, 1) o z in e1 R for all z in Z.
Construction (floating point, following an earlier search over certificates, not included
here): start from the lower hull of {(1/e2, 0), (y,1)/e1 o (1/e2, 0)},
add every product that violates (C1) or (C3) and iterate; bisect for the largest e1 for which the
iteration closes.  The points are then rounded to rationals, Z is rescaled so that (C1) holds exactly,
and e1, e2 are replaced by rationals below the exact antinorms required by (C3) and (C2).
Finally the weights lambda of the rate bound (see the verifier) are found by linear programming
(scipy) and rounded to rationals.

Usage: python3 dst_certificate_build.py dst_certificate_params.json
       (writes dst_certificates.json; about 15 seconds)
"""
import json
from fractions import Fraction as F

EPS = -1e-10


def comp(P, Q):
    return (P[0] * (Q[0] + Q[1]), P[1] * (Q[0] / 2 + Q[1]))


def lower_hull(pts):
    pts = sorted(set(pts))
    pareto = []
    for p in pts:
        if not pareto or p[1] < pareto[-1][1]:
            pareto.append(p)
    hull = []
    for p in pareto:
        while len(hull) >= 2:
            (x1, y1), (x2, y2) = hull[-2], hull[-1]
            if (x2 - x1) * (p[1] - y1) - (y2 - y1) * (p[0] - x1) <= 0:
                hull.pop()
            else:
                break
        hull.append(p)
    return hull


def facets(hull):
    Fs = [((1, 0), hull[0][0]), ((0, 1), hull[-1][1])]
    for (x1, y1), (x2, y2) in zip(hull, hull[1:]):
        w = (y1 - y2, x2 - x1)
        Fs.append((w, w[0] * x1 + w[1] * y1))
    return Fs


def antinorm(P, Fs):
    return min((w[0] * P[0] + w[1] * P[1]) / h for (w, h) in Fs if h > 0)


def certify(y, e1, e2, maxit=80, maxpts=200):
    iS = (y / e1, 1.0 / e1)
    iT = (1.0 / e2, 0.0)
    Z = lower_hull([iT, comp(iS, iT)])
    for _ in range(maxit):
        Fs = facets(Z)
        new = [comp(z, w) for z in Z for w in Z] + [comp(iS, w) for w in Z]
        bad = [q for q in new if antinorm(q, Fs) < 1 + EPS]
        if not bad:
            return Z
        Z = lower_hull(Z + bad)
        if len(Z) > maxpts or min(z[0] + z[1] for z in Z) < 1e-6 * iT[0]:
            return None
    return None


def best_e1(y, e2, it=60):
    lo, hi = 1e-3, (1 + y) * 2
    if certify(y, lo, e2) is None:
        return None
    for _ in range(it):
        mid = (lo * hi) ** 0.5
        if certify(y, mid, e2) is not None:
            lo = mid
        else:
            hi = mid
    return lo


def floor_rat(x, den=10 ** 12):
    return F((x.numerator * den) // x.denominator, den)


def exact(y, e2f):
    """Return the exact rational certificate built from the floating-point search."""
    e1f = best_e1(float(y), float(e2f))
    Zf = certify(float(y), e1f, float(e2f))
    Z = lower_hull([(F(a).limit_denominator(10 ** 12), F(b).limit_denominator(10 ** 12)) for a, b in Zf])
    Fs = facets(Z)
    alpha = min(antinorm(comp(z, w), Fs) for z in Z for w in Z)
    alpha = min(F(1), floor_rat(alpha, 10 ** 15))
    Z = [(a / alpha, b / alpha) for a, b in Z]          # now (C1) holds exactly
    Fs = facets(Z)
    beta = min(antinorm(comp((y, F(1)), z), Fs) for z in Z)
    gamma = antinorm((F(1), F(0)), Fs)
    return {"y": y, "e1": floor_rat(beta), "e2": floor_rat(gamma), "Z": Z}


def rates(c):
    from math import log2
    L = log2(1 + float(c["y"]))
    return (L - log2(float(c["e1"])), -log2(float(c["e2"])), L)


if __name__ == "__main__":
    import sys
    from scipy.optimize import linprog
    # (y, e2) pairs: the three certificates that carry positive weight in the linear programme over the
    # 1,155 floating-point certificates of an earlier search, which is not included here.  This program does not repeat that search; it only rebuilds these certificates.
    PARAMS = json.load(open(sys.argv[1]))
    certs = [exact(F(y), F(e2)) for (y, e2) in PARAMS]
    P = [rates(c) for c in certs]
    n = len(P)
    res = linprog([0] * n + [1], A_ub=[[P[j][i] for j in range(n)] + [-1] for i in range(3)],
                  b_ub=[0, 0, 0], A_eq=[[1] * n + [0]], b_eq=[1],
                  bounds=[(0, None)] * n + [(None, None)], method="highs")
    lam = [max(F(0), F(x).limit_denominator(10 ** 6)) for x in res.x[:n]]
    s = sum(lam)
    lam = [l / s for l in lam]
    out = {"certificates": [{"y": str(c["y"]), "e1": str(c["e1"]), "e2": str(c["e2"]),
                             "Z": [[str(a), str(b)] for a, b in c["Z"]]} for c in certs],
           "lambda": [str(l) for l in lam]}
    json.dump(out, open("dst_certificates.json", "w"), indent=1)
    print(f"LP value {res.fun:.9f}; sizes {[len(c['Z']) for c in certs]}")