"""Supporting computation (floating point): the obstruction family of prop:fa-dst-lower.

The family: a spine r -> p_1 -> p_2 -> ... -> p_{2L}.  An odd spine vertex has one terminal leaf and
one Steiner child (a "tooth") carrying two terminal leaves; an even spine vertex has two terminal
leaves; the last spine vertex has no spine child.  Per period of two spine vertices there are s = 3
Steiner vertices and p = 5 terminals; rho padding non-terminals per period are added (each with arcs
to and from two sibling terminal leaves), so n = 1 + L (8 + rho) and the universe has m = L (3 + rho)
vertices.  The two-state recursion of lem:fa-dst-rec gives, for the message (a', b') of p_{i+1},
the message of p_i as M (a', b') with M = [[A, A], [B/2, B]], where (A, B) is the message of p_i after
its leaves and teeth.  With T(y) = M_odd(y) M_even(y) and lambda(y) its spectral radius, the
generating function grows like lambda(y)^L, and the rate per vertex of the family is
    R(rho) = min_y [ (3 + rho) log2(1 + y) - log2 lambda(y) ] / (8 + rho),
maximised over rho >= 0.  This is a numerical evaluation, not a proof.  The maximum
sits where the interior minimum in y meets the limit y -> 0, whose value 5/(8 + rho) is the cost of
the trial t = 0 (inclusion-exclusion over all 5L terminals).
Usage: python3 dst_obstruction.py
"""
from math import log2, sqrt, exp


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


def parts(y):
    term = (1.0, 0.0)
    tooth = comp(comp((y, 1.0), term), term)          # Steiner with two terminal leaves
    odd = comp(comp((y, 1.0), term), tooth)           # one leaf and one tooth
    even = comp(comp((y, 1.0), term), term)           # two leaves
    return odd, even


def lam(y):
    (A1, B1), (A2, B2) = parts(y)
    M1 = [[A1, A1], [B1 / 2, B1]]
    M2 = [[A2, A2], [B2 / 2, B2]]
    T = [[sum(M1[i][k] * M2[k][j] for k in range(2)) for j in range(2)] for i in range(2)]
    tr, det = T[0][0] + T[1][1], T[0][0] * T[1][1] - T[0][1] * T[1][0]
    return (tr + sqrt(tr * tr - 4 * det)) / 2


def inner(rho, step=0.01):
    """min over y > 0, by a grid in u = ln y on [-30, 6] and golden-section refinement.
    As y -> 0 the bracket tends to 5 / (8 + rho) (the trial t = 0), which the grid reaches."""
    f = lambda u: ((3 + rho) * log2(1 + exp(u)) - log2(lam(exp(u)))) / (8 + rho)
    N = int(36 / step)
    vals = [(f(-30 + step * i), -30 + step * i) for i in range(N + 1)]
    v, u = min(vals)
    lo, hi = u - step, u + step
    g = (sqrt(5) - 1) / 2
    for _ in range(80):
        a, b = hi - g * (hi - lo), lo + g * (hi - lo)
        if f(a) < f(b):
            hi = b
        else:
            lo = a
    u2 = (lo + hi) / 2
    return min(v, f(u2)), exp(u2 if f(u2) < v else u)


if __name__ == "__main__":
    grid = [i * 0.05 for i in range(0, 161)]
    v, r = max((inner(r, 0.02)[0], r) for r in grid)
    lo, hi = max(r - 0.05, 0), r + 0.05
    for _ in range(60):
        a, b = lo + (hi - lo) / 3, hi - (hi - lo) / 3
        if inner(a, 0.002)[0] < inner(b, 0.002)[0]:
            lo = a
        else:
            hi = b
    r = (lo + hi) / 2
    v, ystar = inner(r, 0.002)
    print(f"max_rho R(rho) = {v:.9f} bits = {2 ** v:.6f}^n at rho = {r:.4f}, y = {ystar:.5f}")
    print(f"value of the trial t = 0 there: 5/(8 + rho) = {5 / (8 + r):.9f}")