multipart.py

"""multipart.py -- past r = 51, at the takeover theorem's own e(r), does
the exact centred test keep every shape the closed-form bound misses
above F_e's circle mean, and does a product of cosines bind below the
single parts?

QUESTION. J* = (4e - 1) l - 1, l = r - 2, e = takeph.theorem_e. Every
shape compbound's bound misses, read by exacttest's best X over p: at
r = 52..61 is any below F_e's mean, and does a multi-part shape bind below
the single-part least (danger.py), as k=0 2^26 7 3 did at r = 35 (7.0171
under the part 5's 7.8844)?

BY HAND (re-derived before the code). A shape: t = k extra (x - 1)'s, d
even parts, o odd parts, Pi the even parts' halves, q = d mod 2,
s0 = r - 3. Its pairing is (2^(J+t+d+p+q)/pi) Pi B(al, be) E Phi,
al = (d + q + 1)/2, be = (J + t - o + p + 1)/2, the mean (2^(J+s0)/pi)
B(a, b) A, a = l/2, b = (J + 1)/2, A = E|1 - 4eU|, so
  X/mean = 2^(t+d+p+q-s0) Pi B(al, be)/(|w|_1 B(a, b)) |E Phi| / A,
which at k=0 2^(r-g) g is danger's B1 |E cos gu|/A. U = sin^2 u ~
Beta(al, be). Phi is a cosine polynomial: cos(gu) at odd g, D_j(2u) =
(1/j) sum_(m<j) cos(2(j-1-2m)u) at even g = 2j, products by cos a cos b =
(cos(a+b) + cos(a-b))/2. So Phi = sum_n c_n cos(nu), and the hand form
H = B1 |sum_n c_n cos(n ubar) exp(-n^2 Var u/2)| / A, ubar and Var u
danger.hand's: exact for a Gaussian u, danger.hand itself at one part.
At k=0 2^26 7 3, u0 = arcsin(1/4), Phi = (cos 10u + cos 4u)/2 is about
0.14 against |cos 7u0| = 0.20: the product dips through a
near-cancellation of two frequencies.

PREDICTIONS (frozen before the code). P1 H/exact in [0.95, 1.05] at
k=0 2^26 7 3 at its exact best p (r = 35). P2 at r = 54..61 no failing
shape is below the mean. P3 at r = 54..61 every rank binds at a single
part, the least multi-part shape at least 1.3 times it. P4 the least
multi-part shape per rank carries a part 3 beside its big odd part.
KILL. Printed per rank: the least X/mean over every failing shape, its
shape and p; the least over multi-part shapes, its shape; H/exact at
both. (a) Any X/mean <= 1 + 1e-6: the centred test does not carry
J_hi >= J* there, and the rank owes evenarm's height(F_e) route. (b)
H/exact outside [0.9, 1.1] at a rank's least multi-part shape: the
product hand form cannot carry a uniform margin.
CONTROLS, read first. C1 the hand form equals danger.hand at every
single-part shape of r = 53 and 54 (log, to 1e-12). C2 the pipeline
reproduces exacttest's leasts 7.0171 (r = 35, k=0 2^26 7 3) and 14.0320
(r = 51, k=0 2^44 7). C3 H at k=0 2^26 7 3, r = 35. PLANT: the hand form
with the part 3's cos 3u dropped leaves [0.9, 1.1] of exact there. (The
frozen plant missed, 0.9286, at the first run: cos 7u sits near its zero
there, so the 3 moves H by 7%. Replaced after that run by the product
rule's difference term dropped, Phi read as cos 10u, a plant on the
mechanism the hand form adds, a gross plant and so a weak control; the
frozen one still prints, and its miss says KILL (b)'s band cannot see a
whole factor dropped at that shape.)
A shape is multi-part when it has two parts above 2 or k > 0 (set with
the code, after the freeze).

FINDINGS (copied from the printed output).
  (1) The controls hold. C1: max |log difference| 1.4e-14. C2: 7.0171
      at k=0 2^26 7 3 (p = 44), 14.0320 at k=0 2^44 7 (p = 0). C3:
      H/exact 0.9958 at k=0 2^26 7 3 (P1 holds); the frozen plant read
      0.9286, missed; the difference-term plant 13.1836, caught.
  (2) Both KILLs miss. Every failing shape is above the mean at every
      r = 52..61 (P2 holds); H/exact at the binding and least
      multi-part shapes lies in 0.9990 .. 0.9997. Failing shapes per rank:
      26543 .. 110949 at e = 3 (even r), 2062 .. 9282 at e = 5 (odd r).
  (3) Every rank binds at a single part, danger's least: k=0 2^(r-5) 5 at
      even r, 2.3144 (r = 52) rising to 2.4014 (r = 60); k=0 2^(r-7) 7 at
      odd r, 11.1125 (r = 53) falling to 9.9328 (r = 61). The least
      multi-part shape is k=0 2^(r-7) 5 3 at even r, 6.6473 .. 6.6643
      (2.76 .. 2.88 times the single), and k=0 2^44 7 3 (84.5056) at
      r = 53, then k=1 2^(r-8) 7 at odd r = 55..61, 50.8999 .. 46.8574
      (4.7 times). P3 holds (>= 1.3). P4 fails at odd r = 55..61, where
      the least multi-part shape is an extra (x - 1) beside one part 7:
      a single cosine, multi-part only by the definition set after the
      freeze; the least product of cosines there is unprinted.
READS. BY HAND's near-cancellation was read at u0, the TRANSPLANT it
forbade: at the shape's own ubar (0.2187 at r = 35, p = 44, against
u0 = 0.2527) the hand terms are cos 7u 0.039 and Phi 0.042 (printed
off hand's terms beside the run, not by main), so k=0 2^26 7 3 is
small because 7 ubar sits just under pi/2, and its 3 scales it, as the frozen plant's miss shows.
Whether that dip recurs is where g ubar of the big part sits; unread.
With the bound clearing the rest above the mean, every non-champion pure
product is above F_e's mean at J*, so J_hi(r) >= J* at every r =
52..61, computed. The enumeration grows about 1.38x per two ranks at
e = 3's even r (20.8 s at r = 52 to 75.3 s at r = 60), and the failing
counts move with e (26543 .. 110949 at e = 3, 2062 .. 9282 at e = 5),
so a uniform statement needs the shapes bounded without listing them.

RUN RECORD. 358.6 s, peak commit 132.4 MB.
"""
import math
import os
import sys
import time

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

LN2 = math.log(2)


def phi(parts, drop=None, plant=False):
    """Phi's cosine coefficients {n: c_n}, n >= 0; drop: a part whose
    factor is left out (the frozen PLANT); plant: the product rule's
    difference term dropped (the PLANT)."""
    out = {0: 1.0}
    for g in parts:
        if g == 2 or g == drop:
            if g == drop:
                drop = None
            continue
        if g % 2:
            f = {g: 1.0}
        else:
            j = g // 2
            f = {}
            for m in range(j):
                n = abs(2 * (j - 1 - 2 * m))
                f[n] = f.get(n, 0.0) + 1.0 / j
        new = {}
        for a, x in out.items():
            for b, y in f.items():
                for n in (a + b,) if plant else (a + b, abs(a - b)):
                    new[n] = new.get(n, 0.0) + x * y / (1 if plant else 2)
        out = new
    return out


def hand(r, J, logA, k, parts, p, drop=None, plant=False):
    """log H at (shape, p)."""
    d = takeph.order(parts)
    o = len(parts) - d
    q = d % 2
    pi = 1
    for g in parts:
        if g % 2 == 0:
            pi *= g // 2
    al, be = (d + q + 1) / 2, (J + k - o + p + 1) / 2
    a, b = (r - 2) / 2, (J + 1) / 2
    lB1 = ((k + d + p + q - (r - 3)) * LN2 + math.log(pi)
           - compbound.log_w1(p, q) + danger.lbeta(al, be)
           - danger.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))
    E = sum(c * math.cos(n * ub) * math.exp(-n * n * Vu / 2)
            for n, c in phi(parts, drop, plant).items())
    if E == 0:
        return -math.inf
    return lB1 + math.log(abs(E)) - logA


def multi(k, parts):
    return k > 0 or sum(1 for g in parts if g > 2) > 1


def rank(r):
    """(e, J, failing keys, failing shapes, least, least multi, least
    single), each least (log X/mean, name, p, H/exact)."""
    e, l, J, lm, logA = danger.setup(r)
    _, _, lmean, nkeys, fs = exacttest.failing_shapes(r)
    best = {"all": (math.inf,), "multi": (math.inf,),
            "single": (math.inf,)}
    for k, parts, s, o, pi, t6, plb in fs:
        v, p = exacttest.best_x(J, k, parts)
        v -= lmean
        for key in ("all", "multi" if multi(k, parts) else "single"):
            if v < best[key][0]:
                best[key] = (v, k, parts, p)
    out = {}
    for key, b in best.items():
        if b[0] == math.inf:
            out[key] = None
            continue
        v, k, parts, p = b
        h = hand(r, J, logA, k, parts, p)
        out[key] = (v, takeph.name(k, parts), p, math.exp(h - v))
    return e, J, nkeys, len(fs), out


def controls():
    ok = True
    worst = 0.0
    for r in (53, 54):
        e, l, J, lm, logA = danger.setup(r)
        for g in range(5, r + 1, 2):
            for p in range(J % 2, J % 2 + 60, 2):
                a = hand(r, J, logA, 0, (g,) + (2,) * (r - g), p)
                b = danger.hand(r, J, logA, g, p)[0]
                if a == b == -math.inf:
                    continue
                worst = max(worst, abs(a - b))
    good = worst < 1e-12
    ok &= good
    print("C1 r = 53, 54: hand against danger.hand at every single part,"
          " p < J mod 2 + 60: max |log difference| %.1e  %s"
          % (worst, "holds" if good else "FAILS"), flush=True)
    for r, want, name in ((35, 7.0171, "k=0 2^26 7 3"),
                          (51, 14.0320, "k=0 2^44 7")):
        e, J, nk, nf, out = rank(r)
        v, nm, p, ratio = out["all"]
        good = abs(math.exp(v) - want) < 5e-5 and nm == name
        ok &= good
        print("C2 r = %d: least %.4f at %s p = %d (want %.4f at %s)  %s"
              % (r, math.exp(v), nm, p, want, name,
                 "holds" if good else "FAILS"), flush=True)
    r = 35
    e, l, J, lm, logA = danger.setup(r)
    parts = (7, 3) + (2,) * 26
    v, p = exacttest.best_x(J, 0, parts)
    v -= lm
    h = math.exp(hand(r, J, logA, 0, parts, p) - v)
    fr = math.exp(hand(r, J, logA, 0, parts, p, drop=3) - v)
    pl = math.exp(hand(r, J, logA, 0, parts, p, plant=True) - v)
    caught = not 0.9 <= pl <= 1.1
    ok &= caught
    print("C3 r = 35 k=0 2^26 7 3 p = %d: X/mean %.4f, H/exact %.4f (P1 %s),"
          " frozen plant (cos 3u dropped) %.4f %s, plant (difference term"
          " dropped) %.4f %s"
          % (p, math.exp(v), h, "holds" if 0.95 <= h <= 1.05 else "fails",
             fr, "caught" if not 0.9 <= fr <= 1.1 else "missed", pl,
             "caught" if caught else "MISSED"), flush=True)
    return ok


def main():
    t0 = time.time()
    ok = controls()
    print("controls %s  [%.1f s]" % ("hold" if ok else "FAIL",
                                     time.time() - t0), flush=True)
    if not ok:
        return False
    print()
    print("  r  e  keys shapes   least  shape (p, H/exact)  |  least multi"
          "  shape (p, H/exact)  |  single  [s]")
    ka, kb = [], []
    for r in range(52, 62):
        t = time.time()
        e, J, nk, nf, out = rank(r)
        a, m, s = out["all"], out["multi"], out["single"]
        if a[0] <= math.log(1 + 1e-6):
            ka.append(r)
        if m and not 0.9 <= m[3] <= 1.1:
            kb.append(r)
        print("%3d %2d %5d %6d  %.4f %s (%d, %.4f)  |  %s  |  %s  [%.1f]"
              % (r, e, nk, nf,
                 math.exp(a[0]), a[1], a[2], a[3],
                 "%.4f %s (%d, %.4f)" % (math.exp(m[0]), m[1], m[2], m[3])
                 if m else "-",
                 "%.4f %s" % (math.exp(s[0]), s[1]) if s else "-",
                 time.time() - t), flush=True)
    print("KILL (a) %s" % ("fires at r = %s" % ka if ka else "misses"))
    print("KILL (b) %s" % ("fires at r = %s" % kb if kb else "misses"))
    print("%.1f s" % (time.time() - t0))
    return ok


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