# Search periodic "spine + attached subtrees" families for the largest RS-theta rate
# min_y [log2(1+y) - lim log2 W(y)/n]  (exact asymptotics via transfer matrix).
import math, itertools
import numpy as np
from scipy.optimize import minimize_scalar
def merge(P,Q):
    A,B,C=P; a,b,c=Q
    return np.array([A*(2*a+b+c), B*(a+b+c/2)+C*a, C*(b+c/2)])
def init(y,m): return np.array([y/2,0.0,2.0**-m])
def Mmat(P):
    A,B,C=P
    return np.array([[2*A,A,A],[B+C,B,B/2],[0,C,C/2]])
# library of small rooted subtrees: (m_root, [children subtrees])
def build(y, t):
    m, ch = t
    P = init(y, m)
    for c in ch: P = merge(P, build(y, c))
    return P
def size(t): return 1 + sum(size(c) for c in t[1])
def outs(t): return t[0] + sum(outs(c) for c in t[1])
lib = []
for m in range(0,5): lib.append((m, ()))
for m in range(0,3):
    for k in (1,2):
        for gm in itertools.combinations_with_replacement(range(0,4), k):
            lib.append((m, tuple((g, ()) for g in gm)))
def rate_unit(y, mspine, teeth):
    P = init(y, mspine)
    for t in teeth: P = merge(P, build(y, t))
    rho = max(abs(np.linalg.eigvals(Mmat(P))))
    n = 1 + mspine + sum(size(t)+outs(t) for t in teeth)
    return math.log2(1+y) - math.log2(rho)/n
def best_y(f):
    r = minimize_scalar(f, bounds=(0.02, 3), method='bounded', options={'xatol':1e-10})
    return r.fun, r.x
res = []
for mspine in range(0,3):
    for k in range(0,4):
        for teeth in itertools.combinations_with_replacement(range(len(lib)), k):
            T = [lib[i] for i in teeth]
            v, y = best_y(lambda y: rate_unit(y, mspine, T))
            n = 1 + mspine + sum(size(t)+outs(t) for t in T); s = 1 + sum(size(t) for t in T)
            res.append((v, y, mspine, T, s/n))
res.sort(key=lambda r: -r[0])
for r in res[:8]:
    print(f"rate {r[0]:.7f} base {2**r[0]:.6f} y*={r[1]:.4f} sigma={r[4]:.4f} spine m={r[2]} teeth={r[3]}")
print("searched", len(res), "periodic families; max base", 2**res[0][0], "(claimed envelope of method 1.72028, certified 1.72035)")