danger.py

"""danger.py -- at the takeover theorem's own e(r), does the exact centred
test put every single-odd-part shape above F_e's circle mean?

QUESTION. J* = (4e - 1) l - 1, l = r - 2, e = floor(sqrt(l/2)) at odd r,
floor(sqrt(l)/2) at even r, e = 3..11 (odd r = 21..289, even r =
38..576). The shapes k=0 2^(r-g) g, g odd >= 5: is the exact centred
test, best over p, above F_e's mean at every one?

BY HAND (re-derived before the code). al = (r - g + 1 + q)/2, be =
(J + p)/2, q = (r - g) mod 2, a = l/2, b = (J + 1)/2:
  X/mean = B1 |E cos gu| / A,  B1 = 2^(p+q+3-g)/|w|_1 B(al, be)/B(a, b),
u drawn as cos^(2be-1) sin^(2al-1) on [0, pi/2], A = E|1 - 4eU|,
U ~ Beta(a, b). The hand form H replaces E cos gu by cos(g ubar)
exp(-g^2 Var u/2), ubar and Var u read off the shape's own Beta(al, be)
with Jensen's correction (partfive.py's C5h at general g). The class
read off u0 = arcsin(1/(2 sqrt e)) dies past the first crossing: the
shape's own mean moves g u by O(1) once g ~ sqrt(l) (e = 5, r = 53,
g = 21: g u0 = 4.736, g ubar = 3.834). Scale: X/mean ~ 0.63 l
|cos g ubar| near the first crossing, the per-rank least O(1) with no
e-trend; p rescues at most ~0.33 g/sqrt(8e).

PREDICTIONS (frozen before the code). P1 H/exact in [0.95, 1.05] at
every shape with g <= r/2 and |cos g ubar| >= 0.2, at its exact best p.
P2 per rank the least exact X/mean sits at a part of the first or second
crossing, in [1, 4] at most ranks, with no rising trend in e. P3 at
e = 5 the binding shape is the part-7 shape, least in [1.2, 2.2]. P4
somewhere in e = 3..11 a shape falls below the mean.
KILL. Printed per rank: the least exact X/mean over the family, its g
and p, and H at the same (g, p). H/exact outside [0.9, 1.1] at a rank's
binding shape fires: the hand form cannot carry a uniform margin.
CONTROLS, read first. C1 at e = 3, odd r = 21..33: the part-5 shape's
exact X/mean equals partfive's printed 2.6738 .. 2.7100 (4 digits). C2
at e = 4, odd r = 39..51: the part-7 shape reads exacttest's 9.8599 ..
14.0320. PLANT: H with cos((g - 2) ubar) leaves [0.9, 1.1] of exact at
the part-5 shape at every r of C1. C3 psi1 equals
exacttest.psi at r = 53, every g; at r = 53, 54, 101, 102, 289 the
pipeline's least and below-list equal the full read's (best_x at every
shape; the below-lists are empty at all five).
The pipeline (changed after the first run, which read exact only the
shapes whose hand best fell below 6 and each rank's three least by
hand): every shape is read exact at its hand-best p; exacttest.best_x,
a scan over p that stops 40 steps past its best (not the maximum: at
r = 258, g = 251 the hand-best p reads 6.3 times its value), is read
where that read is below the running least or the mean. Every shape
then has some p whose exact X/mean is at least the printed least, and
none skipped is below the mean; that the least equals the full read's
is C3's, at its five ranks. "Least" is the least best_x read, a lower
bound on the least maximum over p.

FINDINGS (copied from the printed output).
  (1) The controls hold. C1: the part-5 shape reads 2.6738 .. 2.7100 at
      r = 21..33, H/exact 0.9981 .. 0.9991, the cos 3u plant 3.4751 ..
      4.1417 (caught at every r). C2: 9.8599 (r = 39), 14.0320 (r = 51).
      C3: psi1 differs at 0 of 25 shapes; the pipeline's least equals the
      full read at r = 53 (11.1125, g = 7), 54 (2.3361, g = 5), 101
      (38.5160, g = 5), 102 (5.5378, g = 7), 289 (124.0961, g = 5).
  (2) The KILL misses: at all 405 ranks H/exact at the binding shape
      lies in 0.9981 .. 1.0002. P1 at its frozen p, best_x's, is
      read only at the binding shapes (the KILL's H/exact); read at each
      shape's hand-best p it fails at 33 shapes, all at even r =
      192..572 with g >= 0.46 r, H/exact 0.942 .. 1.067.
  (3) All 51165 shapes are read exact at their hand-best p, 405 (one a
      rank) by best_x: no shape falls below the mean at any rank (P4
      fails). The least per (e, parity), at r and g:
        e  odd r                 even r
        3  2.6617 (25, 5)        2.1747 (38, 5)
        4  7.8844 (35, 5)        7.8243 (66, 5)
        5  8.5955 (73, 7)        4.1169 (144, 7)
        6  25.8052 (75, 5)       27.3201 (146, 5)
        7  38.5160 (101, 5)      41.1192 (198, 5)
        8  32.1444 (163, 9)      57.5718 (258, 9)
        9  71.4459 (165, 5)      76.8094 (326, 5)
        10 91.6668 (203, 5)      98.6990 (402, 5)
        11 114.3926 (245, 5)     123.2842 (486, 5)
      The binding part is 5 at 358 ranks, 7 at 33 (all of e = 5), 9 at
      14 (of e = 8: odd r = 139..163 and r = 258). P2 fails: the
      least rises with e, with dips at e = 5 and 8. P3 fails: e = 5 binds at the part-7 shape as
      predicted, but at 4.1169 (even) and 8.5955 (odd), above 2.2.
  (4) Within e = 5 the part-7 shape falls in r: 11.1125 (r = 53) to
      8.5955 (r = 73), 5.5378 (r = 102) to 4.1169 (r = 144).
READS. The dips sit where the first-crossing part's g u0 is just above
pi/2: 7u0 = 1.579 at e = 5 (0.008 above), 9u0 = 1.599 at e = 8. The
shape's own downward shift, about g(g - 3)/(4 sqrt(e) l) (its mean
alone, odd r, p = 0; Jensen adds a further g/(8 sqrt(e) l) downward,
and p moves it further at even r), holds g ubar below the zero, and it
shrinks as l grows, so the margin falls in r at fixed e (at e = 5 the
mean alone would close it near l ~ 400, with Jensen near l ~ 450,
argued, not printed; e = 5's window ends at l = 142). Away from a dip the binding shape is
the part 5, where 5 ubar falls toward 0 as e grows: X/mean over l reads
0.47 (r = 245) to 0.43 (r = 289) at e = 11's odd ranks, 0.25 to 0.23
at its even ones. BY HAND's scale 0.63 l |cos g ubar| held the part
5's prefactor; B1 grows with g, and the part 7 stands above that
scale at e = 5. A dip is an accident
of e: one window of l moves the shift by O(e^-2.5), so the margin at e
is set by how near g u0, read mod pi, sits to pi/2 + shift: a
Diophantine condition on pi/(2 arcsin(1/(2 sqrt e))) against the odd
integers, which the centred test alone cannot bound uniformly. The
design variable is e itself: J_hi >= c r^(3/2) needs only some e of
order sqrt(l) per rank, not the theorem's.

RUN RECORD. 15.7 s, peak commit 29.0 MB. The full read
(best_x at every shape) costs 20.4 s at r = 576 alone, about 25 min
over the range (a separate timing run).
"""
import math
import os
import sys
import time

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import circlemean  # noqa: E402
import compbound  # noqa: E402
import exacttest  # noqa: E402
import takeph  # noqa: E402

LN2 = math.log(2)


def setup(r):
    e, J = takeph.theorem_e(r), takeph.theorem_J(r)
    l = r - 2
    lm = circlemean.log_mean(r, J, e)
    logA = float(circlemean.mp.log(circlemean.mad(
        circlemean.mp.mpf(l) / 2, circlemean.mp.mpf(J + 1) / 2, 4 * e)))
    return e, l, J, lm, logA


def lbeta(x, y):
    return math.lgamma(x) + math.lgamma(y) - math.lgamma(x + y)


def hand(r, J, logA, g, p, shift=0):
    """log H and cos(g ubar) at (g, p); shift: the plant's g - 2."""
    l = r - 2
    q = (r - g) % 2
    a, b = l / 2, (J + 1) / 2
    al, be = (r - g + 1 + q) / 2, (J + p) / 2
    lB1 = ((p + q + 3 - g) * LN2 - compbound.log_w1(p, q)
           + lbeta(al, be) - lbeta(a, b))
    U = al / (al + be)
    VU = al * be / ((al + be) ** 2 * (al + be + 1))
    upp = -(1 - 2 * U) / (4 * (U * (1 - U)) ** 1.5)
    ub = math.asin(math.sqrt(U)) + upp * VU / 2
    Vu = VU / (4 * U * (1 - U))
    c = math.cos((g - shift) * ub)
    if c == 0:
        return -math.inf, c
    return lB1 + math.log(abs(c)) - g * g * Vu / 2 - logA, c


def hand_best(r, J, logA, g, patience=40):
    p = J % 2
    best, since = (-math.inf, p), 0
    while since < patience:
        v = hand(r, J, logA, g, p)[0]
        if v > best[0]:
            best, since = (v, p), 0
        else:
            since += 1
        p += 2
    return best


def exact(r, J, lm, g):
    v, p = exacttest.best_x(J, 0, (g,) + (2,) * (r - g))
    return v - lm, p


def psi1(r, g):
    """exacttest.psi of the shape k=0 2^(r-g) g: 2^(r-g) T_g, built
    without r - g list products (C3 reads it against psi)."""
    return [v << (r - g) for v in exacttest._cheb(exacttest._T, g)]


def controls():
    ok = True
    want1 = {21: 2.6738, 23: 2.6630, 25: 2.6617, 27: 2.6673, 29: 2.6780,
             31: 2.6925, 33: 2.7100}
    for r, w in want1.items():
        e, l, J, lm, logA = setup(r)
        v, p = exact(r, J, lm, 5)
        x = math.exp(v)
        h = math.exp(hand(r, J, logA, 5, p)[0] - v)
        pl = math.exp(hand(r, J, logA, 5, p, shift=2)[0] - v)
        good = e == 3 and abs(x - w) < 5e-5 and not 0.9 <= pl <= 1.1
        ok &= good
        print("C1 r = %d e = %d: part 5 X/mean %.4f (want %.4f) p = %d,"
              " H/exact %.4f, plant %.4f %s  %s"
              % (r, e, x, w, p, h, pl, "caught" if not 0.9 <= pl <= 1.1
                 else "MISSED", "holds" if good else "FAILS"), flush=True)
    want2 = {39: 9.8599, 51: 14.0320}
    for r, w in want2.items():
        e, l, J, lm, logA = setup(r)
        x = math.exp(exact(r, J, lm, 7)[0])
        good = e == 4 and abs(x - w) < 5e-5
        ok &= good
        print("C2 r = %d e = %d: part 7 X/mean %.4f (want %.4f)  %s"
              % (r, e, x, w, "holds" if good else "FAILS"), flush=True)
    return ok


def rank(r, full=False):
    """(e, (least exact X/mean, g, p, H/exact at it), shapes read at every
    p, shapes below the mean, P1 misses, shapes). Every shape is read
    exact at its hand-best p; best_x is read where that read is below the
    running least or the mean (every shape when full)."""
    e, l, J, lm, logA = setup(r)
    rows = []
    for g in range(5, r + 1, 2):
        hp = hand_best(r, J, logA, g)[1]
        lo = exacttest.log_x(J, 0, r - g, 1, psi1(r, g), hp) - lm
        rows.append((lo, g, hp))
    rows.sort()
    least, below, p1, n = (math.inf, 0, 0, 0.0), [], [], 0
    for lo, g, hp in rows:
        hv, c = hand(r, J, logA, g, hp)
        ratio = math.exp(hv - lo)
        if g <= r / 2 and abs(c) >= 0.2 and not 0.95 <= ratio <= 1.05:
            p1.append((g, round(ratio, 3)))
        if not (full or lo < max(least[0], 0.0)):
            continue
        v, p = exact(r, J, lm, g)
        n += 1
        if v < 0:
            below.append((g, p, round(math.exp(v), 4)))
        if v < least[0]:
            least = (v, g, p, math.exp(hand(r, J, logA, g, p)[0] - v))
    return e, least, n, below, p1, len(rows)


def main():
    t0 = time.time()
    ok = controls()
    bad = sum(psi1(53, g) != exacttest.psi((g,) + (2,) * (53 - g))
              for g in range(5, 54, 2))
    ok &= bad == 0
    print("C3 r = 53: psi1 against exacttest.psi, %d of 25 shapes differ  %s"
          % (bad, "holds" if bad == 0 else "FAILS"), flush=True)
    for r in (53, 54, 101, 102, 289):
        a, b = rank(r), rank(r, full=True)
        good = (a[1][1] == b[1][1] and abs(a[1][0] - b[1][0]) < 1e-12
                and a[3] == b[3])
        ok &= good
        print("C3 r = %d: pipeline least %.4f at g = %d (%d of %d read at"
              " by best_x), full %.4f at g = %d  %s"
              % (r, math.exp(a[1][0]), a[1][1], a[2], a[5],
                 math.exp(b[1][0]), b[1][1],
                 "holds" if good else "FAILS"), flush=True)
    print("controls %s  [%.1f s]" % ("hold" if ok else "FAIL",
                                     time.time() - t0))
    if not ok:
        return False
    print()
    print("  r  e  least   g    p  H/exact  read of  below / P1 misses")
    kill, per_e, tot = [], {}, [0, 0]
    rs = [r for r in range(21, 290, 2)] + [r for r in range(38, 577, 2)]
    for r in rs:
        e, (v, g, p, ratio), n, below, p1, m = rank(r)
        tot[0] += n
        tot[1] += m
        x = math.exp(v)
        if not 0.9 <= ratio <= 1.1:
            kill.append(r)
        cur = per_e.setdefault((e, r % 2), [math.inf, 0, 0, 0, 0])
        if x < cur[0]:
            cur[:3] = [x, r, g]
        cur[3] += len(below)
        cur[4] += len(p1)
        print("%3d %2d  %.4f %3d %4d  %.4f  %4d %3d  %s %s"
              % (r, e, x, g, p, ratio, n, m, below or "", p1 or ""),
              flush=True)
    print()
    print(" e parity  least   at r   g  shapes below  P1 misses")
    for (e, par), (x, r, g, nb, np1) in sorted(per_e.items()):
        print("%2d %s  %.4f  %4d %3d  %d  %d"
              % (e, "odd " if par else "even", x, r, g, nb, np1))
    print("%d shapes over %d ranks, each read exact at its hand-best p;"
          " %d read by best_x" % (tot[1], len(rs), tot[0]))
    print("KILL %s" % ("fires at r = %s" % kill if kill else "misses"))
    print("%.1f s" % (time.time() - t0))
    return ok


if __name__ == "__main__":
    sys.exit(0 if main() else 1)