hull.py

Cited on: FLUSH.

To run it: save hull.py, with flush.py, induct.py, ostrowski.py, tie.py beside it; standard library only.

"""
hull.py -- the pinned closure against the carry's contracting
coordinate: the tops are not the bound's, but the window is inductive on
the bound's whole set at c_inf.

QUESTION. saturate.py prints that at c_inf the pinned closure S_PIN
(induct.py) is, per phase t and flag f, every valid key with G in
[G_lo, G_top](t, f), the ends read off the closure. Along a play Z = G +
kappa_t y, the limit of |theta_(t-1)| Q, contracts: Z' = r_t Z + m d
A_t - e B_t, e >= 0, B_t > 0. Dropping e bounds Z by the fixed point
Z_hi(t, f) of a
max-affine contraction. Is G_top(t, f) = T(t, f), the greatest G with a
valid key (G, H) and G + kappa_t y <= Z_hi(t, f), at every row?

DESIGN (frozen before the engine).
Z(G, H; t) = lambda_t G + kappa_t H, kappa_t = s_t / (1 + r_t s_t), s_t =
[0; a_t, a_(t-1), ...] (the limit of q_(n-1)/q_n, t = n mod p), lambda_t
= 1 - kappa_t r_t, all exact in the quadratic field. A_t = lambda_(t+1)
U_t + kappa_(t+1) V_t (the reply's step m d (U, V)), B_t = lambda_(t+1)
(the move's step -e on G). THE HULL: Z_hi(t+1, R) = max over f with
d_max(f) >= 1 of r_t Z_hi(t, f) + m d_max(f) A_t, d_max(full) =
a_(t+c+2), d_max(reduced) = a_(t+c+2) - 1; Z_hi(t+1, F) = max(r_t max_f
Z_hi(t, f), the pinned starts' Z at t+1); solved exactly on the policy a
float iteration picks, the policy then checked exactly (Bellman).
c_inf is induct.py's c_pin, equal to the pinned game's at all 48 cells
(induct.py), so no game is solved. PREDICTION: T = G_top at every (t,
f) row of the 48 cells. KILL: "T != G_top at n of N (t, R) rows" or
"... (t, F) rows", n > 0 on either. A row with T < G_top is an
instrument fault (the bound is proved), printed apart. THE LOW SIDE (no kill; predicted loose): the hull with d_min and
e_max = a_(t+2) gives Z_lo, and G_lo's hull value is the least G with a
valid key and Z >= Z_lo. CONTROLS: (a) Z's linear part, kappa_(t+1) =
r_t lambda_t and lambda_(t+1) - a_(t+1) kappa_(t+1) = r_t kappa_t, exact at
every phase; A_t > 0; kappa_t, lambda_t, A_t against q_(n-1)
|theta_(n-1)|, q_n |theta_(n-1)|, q_(n+c+1) |theta_n| at n near 40 p;
(b) THE BOUND: every key of S has Z <= Z_hi(t, f), exactly; (c) THE
PLANT: d_max(full) lowered by 1 must put keys of S above the planted
Z_hi at some cell and its T below G_top at some row.

SECOND QUESTION, HULL-WI (frozen after the first verdict, before its
engine). A move's validity reads the offset alone: y' = -y / r_t + m d
(V_t - U_t r_(t+1)) + e r_(t+1), G's coefficient 1 - (a_(t+1) + r_(t+1))
r_t being 0. THE HULL SET H(c): every valid key with Z_lo <= Z <= Z_hi,
closed under every reply and valid move by the fixed points' Bellman
inequalities, holding the pinned starts, so S_PIN(c) inside H(c) and
c_hull (least c with the pinned moves and WI on H(c)) >= c_inf.
PREDICTION: c_hull = c_inf at all 48 cells; every pre-read's valid
landing in H(c_inf). KILL: "c_hull != c_inf at n of 48 cells" or
"landings outside H(c_inf)" > 0. CONTROLS: (a) the WI routine on S_PIN
gives c_inf; (b) S_PIN(c_inf) inside H(c_inf); (c) THE PLANT, in the
kill's direction: Z_lo lowered by 3 at golden x2, c = 3, admits the key
(-5, -3), whose offset has no valid move after the zero reply, and must
fail.

FINDINGS (the --wide run, copied from its prints).
  KILL FIRED. T != G_top at 104 of 104 (t, R) rows and at 100 of 104
      (t, F) rows; T - G_top is 1 at 84, 2 at 15, 3 at 5 reduced rows,
      and 0 at 4, 1 at 72, 2 at 17, 3 at 11 full rows (golden x2: T 10
      against G_top 9 at t0R, Z_hi 9.9193; T 6 against 5 at t0F). T <
      G_top at 0 rows. The bound dropping e is loose by 0 to 3 carries
      at the top (1 to 3 at 204 of 208 rows): e >= 0 alone is not
      enough; what the bound misses (moves the window forces, a sup not
      attained, the lattice's grain) is not separated.
  THE HULL IS THE CARRY'S LIMIT. Z_hi(t, R) equals the carry bound's max
      over the last 3 p of t <= 40 p, m |theta_(t-1)| (q_(t+c+1) - 1), at
      104 of 104 (t, R) rows: the bound with e dropped is the limit the
      carry bound tends to. On the full flag Z_hi is far below it
      (golden x2: 6.13 against 9.92).
  THE LOW SIDE. G_lo - the low hull's value is 1 to 4 on the reduced
      flag and 2 to 4 on the full flag over the 104 rows each, 0 at none:
      the crude floor (d_min, e_max) is loose on both flags.
  CONTROLS. (a) 0 exact failures of Z's linear part, A > 0 and B > 0;
      kappa, lambda and A against q |theta| at n near 40 p, worst 2.2e-12;
      0 Bellman failures of the exact fixed points. (b) 0 of 17983 keys of
      S with Z > Z_hi. (c) the plant puts keys of S above it at 48 of 48
      cells and its T below G_top at 208 of 208 rows: (b) can fire.
  HULL-WI, KILL MISSED. c_hull != c_inf at 0 of 48 cells; valid landings
      outside H(c_inf): 0; pre-reads with no valid landing at c_inf: 0.
      H(c_inf) holds 3 to 20 keys beyond S_PIN (470 over the 48 cells), and
      the window is inductive on them all. So at each cell c_pin = c_inf
      is the least c at which the window is inductive on a set cut out by
      the window and two exact inequalities on Z, no closure grown and no
      game solved, and with COVER, the landings inside, the fixed points'
      Bellman certificates and the completion reader's theorem, c_saf <=
      c_inf rests on that check.
  HULL-WI CONTROLS. (a) the WI routine on S_PIN gives c_inf at 48 of 48.
      (b) H(c_inf) holds S_PIN(c_inf) at 48 of 48. (c) Z_lo - 4 at golden
      x2 adds 4 keys, (-5, -3) among them, and WI fails first at (-5, -3,
      full, t 0) under the zero reply: the kill can fire. (d) successors
      of keys of H(c_inf) outside it: 0 (the closure, proved, read
      directly; no plant shows that print can fire).

RUN RECORD. Bare 0.5 s, --wide 1.0 s. The first bare run broke on the
plant at golden (d_max(full) - 1 = 0 leaves the reduced node with no
predecessor); such a node now takes a constant far below every key. The
OVERSHOOT print was added after the bare run, before the --wide one.
HULL-WI: bare 0.7 s, --wide 1.9 s. Its frozen plant, Z_lo - 3, missed
(-5, -3) (Z_lo(0, full) = -1.894, the key's Z = -4.960); it was resized
to Z_lo - 4 after the bare run, before the verdict was read. COVER and
(d) were added after the first --wide run, which printed the kill and
(a) to (c) as now. The landing kill as first frozen counted cells; this
DESIGN and the print count landings and list their cells. Added to the
frozen DESIGN after the verdict, with no change to the run: "c_inf is
induct.py's c_pin" and the plant's frozen amount "by 3". The Bellman
count at every c stage_wi solves, and the label "max over the last 3 p"
(it read "sup over t <= 40 p"), were added after the verdict: 0
failures, the rows unchanged.

Usage: python hull.py [--wide]   (--wide adds m = 6 .. 9)
"""
import math
import os
import sys
import time
from fractions import Fraction as Fr

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
from flush import PERIODIC  # noqa: E402
from induct import Field, Window, least, QF, le, carry_bounds  # noqa: E402

WIDE = "--wide" in sys.argv


def inv(x):
    n = x.a * x.a - x.b * x.b * QF.D
    return QF(x.a / n, -x.b / n)


def eq(x, y):
    return (x - y).sign() == 0


def isqrt_fr(q):
    """sqrt of a non-negative rational square, else None."""
    a, b = q.numerator, q.denominator
    ra, rb = math.isqrt(a), math.isqrt(b)
    return Fr(ra, rb) if ra * ra == a and rb * rb == b else None


def back_ratio(F, t):
    """s_t = [0; a_t, a_(t-1), ...], exact."""
    M = [[1, 0], [0, 1]]
    for j in range(t, t - F.p, -1):
        b = F.a(j)
        M = [[M[0][1], M[0][0] + b * M[0][1]],
             [M[1][1], M[1][0] + b * M[1][1]]]
    (m00, m01), (m10, m11) = M
    d = Fr((m11 - m00) ** 2 + 4 * m10 * m01)
    k = isqrt_fr(d / QF.D)
    assert k is not None, "backward period outside the field"
    return QF(Fr(-(m11 - m00), 2 * m10), k / (2 * m10))


class Coords:
    """kappa_t, lambda_t, A_t, B_t of one cell."""

    def __init__(self, F, m, c):
        self.F, self.m, self.c, p = F, m, c, F.p
        self.s = [back_ratio(F, t) for t in range(p)]
        self.kap = [self.s[t] * inv(QF(1) + F.r[t] * self.s[t])
                    for t in range(p)]
        self.lam = [QF(1) - self.kap[t] * F.r[t] for t in range(p)]
        self.A, self.B = [], []
        for t in range(p):
            t1 = (t + 1) % p
            U, V = F.uv(t, c)
            self.A.append(self.lam[t1] * U + self.kap[t1] * V)
            self.B.append(self.lam[t1])

    def Z(self, G, H, t):
        return self.lam[t] * G + self.kap[t] * H

    def control_a(self):
        F, p = self.F, self.F.p
        bad = 0
        for t in range(p):
            t1 = (t + 1) % p
            bad += not eq(self.kap[t1], F.r[t] * self.lam[t])
            bad += not eq(self.lam[t1] - self.kap[t1] * F.a(t + 1),
                          F.r[t] * self.kap[t])
            bad += self.A[t].sign() <= 0 or self.B[t].sign() <= 0
        T = 40 * p + p
        q = {-1: 0, 0: 1}
        for k in range(0, T + self.c + 3):
            q[k + 1] = F.a(k + 1) * q[k] + q[k - 1]
        th = {-1: 1.0}
        for k in range(0, T + 1):
            th[k] = th[k - 1] * F.r[k % p].f
        worst = 0.0
        for n in range(40 * p, 41 * p):
            t = n % p
            worst = max(worst, abs(q[n - 1] * th[n - 1] - self.kap[t].f),
                        abs(q[n] * th[n - 1] - self.lam[t].f),
                        abs(q[n + self.c + 1] * th[n] - self.A[t].f))
        return bad, worst


def dmax(F, c, t, f):
    cap = F.a(t + c + 2)
    return cap if f else cap - 1


def hull(win, co, starts, plant=False, low=False):
    """{(t, f): Z bound}: the top (or, low, the floor) hull, exact."""
    F, m, c, p = win.F, win.m, win.c, win.F.p
    sZ = {}
    for G, H, f, t in starts:
        z = co.Z(G, H, t)
        if (t, f) not in sZ or (z.f < sZ[(t, f)].f if low
                                else z.f > sZ[(t, f)].f):
            sZ[(t, f)] = z
    sgn = -1 if low else 1

    def options(t1, f1):
        """[(pred node or None, coefficient, constant)]: the node's value
        is the best of coef * value(pred) + const."""
        t = (t1 - 1) % p
        out = []
        if (t1, f1) in sZ:
            out.append((None, QF(0), sZ[(t1, f1)]))
        emax = F.a(t1 + 1) if low else 0
        for f in (True, False):
            if f1:
                d = 0
            else:
                dm = dmax(F, c, t, f) - (1 if plant and f else 0)
                if dm < 1:
                    continue
                d = 1 if low else dm
            out.append(((t, f), F.r[t], co.A[t] * (m * d)
                        - co.B[t] * emax))
        if not out:                    # unreachable node (the plant)
            out.append((None, QF(0), QF(-sgn * 10 ** 6)))
        return out

    nodes = [(t, f) for t in range(p) for f in (True, False)]
    opts = {n: options(*n) for n in nodes}
    val = {n: 0.0 for n in nodes}
    for _ in range(3000):
        val = {n: sgn * max(sgn * (k.f * (val[pr] if pr else 0) + b.f)
                            for pr, k, b in opts[n] if pr is None
                            or pr in val) for n in nodes}
    pol = {}
    for n in nodes:
        pol[n] = max(opts[n], key=lambda o: sgn * (o[1].f * (
            val[o[0]] if o[0] else 0) + o[2].f))
    exact = {}

    def solve(n):
        path, seen = [], {}
        cur = n
        while cur is not None and cur not in exact and cur not in seen:
            seen[cur] = len(path)
            path.append(cur)
            cur = pol[cur][0]
        if cur is None:
            base = None
        elif cur in exact:
            base = exact[cur]
        else:
            cyc = path[seen[cur]:]
            S, P = QF(0), QF(1)
            for x in cyc:
                _pr, k, b = pol[x]
                S = S + P * b
                P = P * k
            exact[cur] = S * inv(QF(1) - P)
            path = path[:seen[cur]]
            base = exact[cur]
        for x in reversed(path):
            pr, k, b = pol[x]
            exact[x] = b + (k * base if pr is not None else QF(0))
            base = exact[x]
    for n in nodes:
        if n not in exact:
            solve(n)
    bell = 0
    for n in nodes:
        for pr, k, b in opts[n]:
            v = b + (k * exact[pr] if pr is not None else QF(0))
            bell += not (le(exact[n], v) if low else le(v, exact[n]))
    return exact, bell


def top(win, co, t, f, zhi):
    """T(t, f): the greatest G with a valid key and Z <= zhi, exact."""
    lo, hi = win.V[(t, f)]
    r, k = win.F.r[t], co.kap[t]
    G = math.floor(zhi.f - k.f * lo.f) + 2
    for _ in range(100000):
        base = r * G
        for H in range(math.floor((lo + base).f) - 1,
                       math.ceil((hi + base).f) + 2):
            if win.ok(G, H, f, t):
                if le(co.Z(G, H, t), zhi):
                    return G
                break
        G -= 1
    raise RuntimeError("no top")


def bottom(win, co, t, f, zlo):
    """The least G with a valid key and Z >= zlo (floats: a side print)."""
    lo, hi = win.V[(t, f)]
    r, k = win.F.r[t], co.kap[t]
    G = math.ceil(zlo.f - k.f * hi.f) - 2
    for _ in range(100000):
        base = r * G
        for H in range(math.ceil((hi + base).f) + 1,
                       math.floor((lo + base).f) - 2, -1):
            if win.ok(G, H, f, t):
                if co.Z(G, H, t).f >= zlo.f - 1e-9:
                    return G
                break
        G += 1
    raise RuntimeError("no bottom")


def hull_set(win, co, zlo, zhi):
    """Every valid key with Z_lo <= Z <= Z_hi, exact."""
    out = set()
    for (t, f), z1 in zhi.items():
        z0 = zlo[(t, f)]
        lo, hi = win.V[(t, f)]
        r, k = win.F.r[t], co.kap[t]
        for G in range(math.floor(z0.f - k.f * hi.f) - 2,
                       math.floor(z1.f - k.f * lo.f) + 3):
            base = r * G
            for H in range(math.floor((lo + base).f) - 1,
                           math.ceil((hi + base).f) + 2):
                if win.ok(G, H, f, t):
                    z = co.Z(G, H, t)
                    if le(z0, z) and le(z, z1):
                        out.add((G, H, f, t))
    return out


def wi(win, keys):
    """The first key with a reply and no valid move, or None."""
    p = win.F.p
    for G, H, f, t in sorted(keys):
        for d in win.replies(f, t):
            G1, H1 = win.pre(G, H, t, d)
            if not win.moves(G1, H1, d == 0, (t + 1) % p):
                return (G, H, f, t, d)
    return None


def stage_wi(maps):
    """HULL-WI: the least c with WI on the hull set, against c_inf."""
    t0 = time.time()
    kill, n, ctl_a, ctl_b, land_out, land_cells = 0, 0, 0, 0, 0, []
    plant = None
    sizes = []
    cover = esc = bell = 0
    for name, per in PERIODIC:
        F = Field(per)
        p = F.p
        for m in maps:
            c_inf = least(F, m, "pin")
            c_hull = c_mine = None
            for c in range(14):
                win = Window(F, m, c)
                st, miss = win.pin_starts()
                if miss:
                    continue
                co = Coords(F, m, c)
                zhi, b1 = hull(win, co, st)
                zlo, b2 = hull(win, co, st, low=True)
                bell += b1 + b2
                _ok, S, _bad = win.closure(st)
                if c_mine is None and wi(win, S) is None:
                    c_mine = c
                Hs = hull_set(win, co, zlo, zhi)
                if c == c_inf:
                    ctl_b += not S <= Hs
                    cover += win.pre_starts()[1]
                    for G, H, f, t in Hs:
                        t1 = (t + 1) % p
                        for d in win.replies(f, t):
                            G1, H1 = win.pre(G, H, t, d)
                            for G2, H2 in win.moves(G1, H1, d == 0, t1):
                                esc += (G2, H2, d == 0, t1) not in Hs
                    sizes.append(len(Hs) - len(S))
                    for digs, N, _P in win.prereads():
                        f = digs[-1] == 0
                        for G, H in win.landings(N, f):
                            if (G, H, f, 0) not in Hs:
                                land_out += 1
                                cell = f"{name.split()[0]} x{m}"
                                if cell not in land_cells:
                                    land_cells.append(cell)
                    if name.startswith("golden") and m == 2:
                        z3 = {k: v - 4 for k, v in zlo.items()}
                        H3 = hull_set(win, co, z3, zhi)
                        plant = (len(H3 - Hs), wi(win, H3),
                                 any(k[0] == -5 and k[1] == -3 for k in H3))
                if c_hull is None and wi(win, Hs) is None:
                    c_hull = c
                if c_hull is not None and c_mine is not None:
                    break
            n += 1
            kill += c_hull != c_inf
            ctl_a += c_mine != c_inf
            print(f"    {name:20s} x{m} c_inf {c_inf} c_hull {c_hull} "
                  f"(WI routine on S_PIN: {c_mine}); |H| - |S| at c_inf "
                  f"{sizes[-1]}", flush=True)
    print("== HULL-WI: the window on the hull set")
    print(f"  (a) the WI routine on S_PIN gives c_inf at {n - ctl_a} of {n}")
    print(f"  (b) H(c_inf) holds S_PIN(c_inf) at {n - ctl_b} of {n}")
    print(f"  (c) THE PLANT (Z_lo - 4, golden x2): {plant[0]} keys added, "
          f"G = -5, H = -3 among them {plant[2]}, first WI failure "
          f"{plant[1]}")
    print(f"  (d) THE CLOSURE OF H, read: successors of keys of H(c_inf) "
          f"outside it {esc}")
    print(f"  COVER: pre-reads with no valid landing at c_inf {cover}")
    print(f"  Bellman failures of the fixed points at every c solved: "
          f"{bell}")
    print(f"  keys of H beyond S at c_inf: {min(sizes)} to {max(sizes)}, "
          f"sum {sum(sizes)}")
    print(f"  LANDING: valid landings outside H(c_inf): {land_out} (cells "
          f"{', '.join(land_cells) or '-'})")
    print(f"  KILL: c_hull != c_inf at {kill} of {n} cells")
    print(f"  {time.time() - t0:.1f} s")


def ranges(S):
    out = {}
    for G, _H, f, t in S:
        lo, hi = out.get((t, f), (G, G))
        out[(t, f)] = (min(lo, G), max(hi, G))
    return out


def main():
    t0 = time.time()
    maps = (2, 3, 4, 5) + ((6, 7, 8, 9) if WIDE else ())
    a_bad, a_worst, bell_bad = 0, 0.0, 0
    b_bad, b_keys = 0, 0
    pl_above_cells, pl_low_rows = 0, 0
    kill = {True: 0, False: 0}
    nrows = {True: 0, False: 0}
    under = 0
    hunch_eq, hunch_n = 0, 0
    low_rows = {True: [], False: []}
    over = {True: [], False: []}
    fails = []
    for name, per in PERIODIC:
        F = Field(per)
        for m in maps:
            c = least(F, m, "pin")
            win = Window(F, m, c)
            co = Coords(F, m, c)
            ab, aw = co.control_a()
            a_bad += ab
            a_worst = max(a_worst, aw)
            st, _ = win.pin_starts()
            _ok, S, _bad = win.closure(st)
            rng = ranges(S)
            zhi, bb = hull(win, co, st)
            bell_bad += bb
            zlo, bb = hull(win, co, st, low=True)
            bell_bad += bb
            bad_here = sum(not le(co.Z(G, H, t), zhi[(t, f)])
                           for G, H, f, t in S)
            b_bad += bad_here
            b_keys += len(S)
            zpl, _ = hull(win, co, st, plant=True)
            pl_above_cells += any(not le(co.Z(G, H, t), zpl[(t, f)])
                                  for G, H, f, t in S)
            bounds = carry_bounds(F, m, c, 40 * F.p)
            row = []
            for (t, f), (g0, g1) in sorted(rng.items()):
                T = top(win, co, t, f, zhi[(t, f)])
                over[f].append(T - g1)
                nrows[f] += 1
                if T != g1:
                    kill[f] += 1
                    under += T < g1
                    if len(fails) < 12:
                        fails.append(f"    {name:20s} x{m} t{t}"
                                     f"{'F' if f else 'R'}: T {T}, "
                                     f"G_top {g1}, Z_hi {zhi[(t, f)].f:.4f}")
                pl_low_rows += top(win, co, t, f, zpl[(t, f)]) < g1
                B = bottom(win, co, t, f, zlo[(t, f)])
                low_rows[f].append(g0 - B)
                if not f:
                    L = max(hi for ph, _lo, hi, _k in bounds[-3 * F.p:]
                            if ph == t)
                    hunch_n += 1
                    hunch_eq += abs(L - zhi[(t, f)].f) < 1e-6
                row.append(f"t{t}{'F' if f else 'R'} T {T}/{g1} "
                           f"Z {zhi[(t, f)].f:.2f}")
            print(f"    {name:20s} x{m} c {c}  " + "  ".join(row),
                  flush=True)
    print("== HULL: the closure's tops against the contracting bound")
    print(f"  (a) Z's linear part, A > 0, B > 0: {a_bad} exact failures; "
          f"kappa, lambda, A against q |theta| at n near 40 p: worst "
          f"{a_worst:.1e}; Bellman failures of the exact fixed points "
          f"{bell_bad}")
    print(f"  (b) THE BOUND: keys of S with Z > Z_hi: {b_bad} of {b_keys}")
    print(f"  (c) THE PLANT (d_max(full) - 1): cells with keys of S above "
          f"it {pl_above_cells}; rows with its T below G_top {pl_low_rows}")
    print(f"  the carry bound's max over the last 3 p of t <= 40 p "
          f"(m |theta| (q - 1)) equals Z_hi(t, R) at {hunch_eq} of {hunch_n} (t, R) rows")
    for f in (False, True):
        g = low_rows[f]
        print(f"  LOW SIDE, {'full' if f else 'reduced'}: G_lo - hull "
              f"value over {len(g)} rows: min {min(g)}, max {max(g)}, "
              f"0 at {g.count(0)}")
    for f in (False, True):
        g = over[f]
        print(f"  OVERSHOOT, {'full' if f else 'reduced'}: T - G_top over "
              f"{len(g)} rows: " + ", ".join(
                  f"{v} at {g.count(v)}" for v in sorted(set(g))))
    print(f"  T < G_top (instrument fault): {under} rows")
    if fails:
        print("  first mismatches:")
        print("\n".join(fails))
    print(f"  KILL: T != G_top at {kill[False]} of {nrows[False]} (t, R) "
          f"rows; at {kill[True]} of {nrows[True]} (t, F) rows")
    print(f"  {time.time() - t0:.1f} s")
    stage_wi(maps)


if __name__ == "__main__":
    main()