taylorn.py
Cited on: FLATTEN.
To run it: save taylorn.py, with band.py, bandend.py, bandfamily.py, circlemean.py, compbound.py, exacttest.py, flatten.py, prunerr.py, takeph.py beside it; needs mpmath, numpy, scipy, sympy.
"""taylorn.py -- does one fixed Taylor order, on nu's exact moments, carry
every shape the squared test's exact E Q_S clears at e = 3..5?
python taylorn.py # r = 21..33: controls, KILL 1
python taylorn.py --to 61 > run.txt # every rank, about 4,240 s
python taylorn.py --kill run.txt 60,8 # both kills off the saved rows
python taylorn.py --reread 60 "k=0 2^37 9 7 5 5" 1000000
The bare run reads r = 21..33 in about 25 s: C2 and C1's E sin^2,
E tan^2 part then cover r <= 33 alone, and it prints KILL 1 over those
ranks and no KILL 2, which needs r = 50..60. The full run is the arm
that reads both: no rows are kept, so --kill reads the rank rows of a
run saved to a file, a cut failure counted only once named as re-read
(60,8: r = 60 at n = 8), and --reread re-reads one shape's prefixes
with the dead-scan close rule (scan) at the DEADCAP given.
QUESTION. The squared test reads E_nu Q_S = sum_v P(X = v) psi(2v)
exact in integers (psi(t) = E_nu cos tu; Q_S the product over the
parts g in S of cos^2(gu), odd g, or (sin(gu)/((g/2) sin 2u))^2, even
g; nu the density cos^A sin^Sq on [0, pi/2]) and clears every shape of
compbound's failing set (exacttest.failing_shapes) at e >= 3 but the
part-5 family.
A fourth-order Taylor of the low frequencies about ubar on nu's exact
moments, the arm below at n = 4, fails shapes past that family at
r = 29..33. At every r = 21..61 whose e (takeph.theorem_e) is 3, 4 or
5, over the shapes the exact test clears, does one order n in 4, 6, 8,
10, 12 clear them all?
DESIGN (by hand, frozen before this code). The arm at order n
(n even): E Q_low >= sum_(j < n, j != 1) Q^(j)(ubar) mo_j / j!
- M_n mo_n / n!, Q^(j)(u) = sum_L P(v) (2v)^j cos(2vu + j pi/2), M_n =
sum_L P(v) (2v)^n, mo_j nu's exact central moments (Gauss-Legendre, 256
nodes), L the frequencies with 4v^2/kappa <= tau; the high part bounded
by the total variation of nu's density and of its derivative,
|E cos 2vu| <= min(1, TV/2v, TV2/4v^2); the best over tau in TAUS. The
domain lemma: the
arm is a lower bound on E Q_S and LB is increasing in E, so the arm can
clear only where the exact test does; it is read on D(S), the p of the
exact scan (patience 40, DEADCAP 1000) whose LB clears the
bar, the scan extended while it does. A shape fails at n when no prefix
S and no p in D(S) clears; a failing shape with a prefix whose exact
scan was cut is counted apart (its D unread past the cut).
PREDICTIONS (frozen). P1 no arm value above the exact E Q_S.
P2 at n = 4 the failures rise with r at e = 3 (more at even r = 50..60
than odd r = 29..33). P3 e = 4, 5 clear at n = 6. P4 N_odd = 6, N_even
<= 8. P5 every rank with e >= 3 clears at some n <= 10. P6 e = 2 fails
more, some rank at no n <= 12. P7 the exact pipeline about 1290 s, the
arm half that again.
KILL (frozen, printed as KILL 1 and KILL 2). (1) some rank with e >= 3
fails, at every n, some shape the exact test clears, with no cut prefix.
(2) N_odd, N_even = the least n failing no shape at any odd r = 21..33,
even r = 50..60 (inf if none); N_even - N_odd >= 4 fires.
CONTROLS, read first. C1 E sin^2 u, E tan^2 u against al/(al + be),
al/(be - 1) at every (A, Sq) met (relative 1e-10); the moments 2..12 at
256 nodes against 512 (relative 1e-9) at three (A, Sq). C2 the arm at
or below the exact E Q_S at every evaluated point; PLANT: the remainder
dropped puts it above somewhere. C3 at n = 4, none failing at r =
21..27, 2, 4, 6 at r = 29, 31, 33 with the shapes C3 names (the
order-4 arm's known failures). C4 the witnesses at r = 33, p = 0,
computed by an earlier, unsaved arm (S = (9, 7), A = 338, Sq = 18: -0.0155,
0.0226, 0.0271 at n = 4, 6, 8, exact 0.0288; S = (11, 7), Sq = 16:
-0.0351, 0.0087 at n = 4, 6, exact 0.0276); the exact Horner against
mpmath. C5 the population 14 (r = 21), 442 (r = 33); the exact test's
failures as want_exact lists them: 5 to 22 at e = 2 (a count, the
shapes unchecked), none at the other ranks r <= 27 and at e >= 4, the
part-5 family k=0 2^(r-5) 5 at the e = 3 ranks past 27, with
k=0 2^51 5 5 at r = 60.
FINDINGS (copied from the printed output; rows "n:genuine+cut").
(1) The controls hold. C1 the even moments at 256 against 512 nodes
within 7.8e-15 at (338, 18), (338, 16), (1200, 40), not the
frozen r = 61's largest A nor the odd moments; E sin^2, E tan^2
within 1.1e-14 (42445 (A, Sq), r <= 33) and 6.3e-14 (r = 61). C4
the witnesses reproduced to four
digits: (9, 7) -0.0155 0.0226 0.0271, then 0.0286, 0.0288 at
n = 10, 12 (exact 0.0288); (11, 7) -0.0351 0.0087, then 0.0245,
0.0267, 0.0275 (exact 0.0276); Horner = mpmath to 12 digits. C3 at
n = 4: none at r = 21..27, 2, 4, 6 at r = 29, 31, 33 with every
named shape. C5 14 and 442; the exact test's failures match
want_exact at every rank. C2 the arm never above exact at 5151490 points
(r <= 33) and at r = 61; the plant above at 2618473 (r <= 33).
C2 over r = 34..60 is unread (the run's print at its last rank
was lost, below); the kill reads failures, which a valid bound
does not need.
(2) KILL 2 FIRES: N_odd = 6, N_even = 10. KILL 1 misses: every rank
with e >= 3 clears at some n <= 10. At e = 3 the least n climbs
along even r: 6 at r = 38..44 (n = 4 failing 23 .. 84), 8 at
r = 46..58 (n = 6 failing 1 .. 7 genuine; the first k=0 2^27
9 7 5), 10 at r = 60, where n = 8 fails one shape, k=0 2^37 9 7 5
5, by a cut its re-read made genuine: the three shorter prefixes'
scans close with D empty, and on (9, 7, 5, 5) (|D| 187, exact
best 5.35 above the bar in log) the arm's best f_n/need is 0.979
at n = 8, 32.5 at n = 10 (-29.8, -2.72 at n = 4, 6). Odd r at
e = 3: 4 at r = 21..27, 6 at 29..33. e = 4 (odd r = 35..51): 4 at
r = 35..39, 6 at 41..51. e = 5 (odd r = 53..61): 4 at 53, 55; 6 at
57, 59, 61 by one cut alone at n = 4 (k=0 2^(r-15) 9 7), unread.
e = 2 (beside): 8, 6, 10, 10, 12, 12, 12, 12 at r = 22..36.
(3) Predictions: P1 holds where read (above); P2 holds (n = 4: 2 .. 6
genuine at odd r = 29..33, 236 .. 334 at even r = 50..60); P3
holds (the three e = 5 cuts unread); P4 fails (N_even 10); P5
holds; P6 fails (every e = 2 rank clears at n = 12); P7 fails
(the ranks 3989.4 + 153.6 = 4143 s against about 1935).
READS (argued). Order 10 carries every rank with e >= 3 to r = 61,
order 12 every rank (C2 unread over r = 34..60, so there the clearing
rests on the arm's validity): what fired is the rise, the order needed
at e = 3 climbing 6, 8, 10 along even r = 38..60, so whether any fixed order
holds past r = 61 is unread and the trend says it may not. The
printed witnesses, at most three a rank, carry a 7 beside 5s, 3s or a
9 at e = 3 and 11 7 at e = 4; the failing sets are counted, not
listed.
RUN RECORD. r = 21..60: the ranks 3989.4 s, the run 4082.9 s, refused
at the 512 MB ceiling inside r = 61 (peak 514.3 MB: one quadrature call
on 7464 A at 256 nodes); r = 60's rank 1302.5 s. nu is now read in
chunks of 1000 A: r = 61 alone 153.6 s, peak 202.3 MB. The re-read
0.4 s; --kill reads the saved rows. The dead scans were most of the
cost, by argument (3173 of r = 60's n = 4 failures carry a cut prefix,
each cut 1000 exact steps); the
re-read's close rule, LB <= head(p) < bar and falling, ends a dead
scan rigorously and was not in the main run. The bare run, r = 21..33:
23.5 s, peak 133.6 MB, the rows and controls as above.
"""
import math
import os
import sys
import time
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import mpmath as mp # noqa: E402
import numpy as np # noqa: E402
from scipy.special import gammaln # noqa: E402
sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "code"))
import circlemean # noqa: E402
import compbound # noqa: E402
import exacttest # noqa: E402
import takeph # noqa: E402
LN2 = math.log(2)
LNPI = math.log(math.pi)
ORDERS = (4, 6, 8, 10, 12)
TOP = max(ORDERS)
TAUS = (0.0, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0, 16.0, math.inf)
NODES = 256
DEADCAP = 1000
FACT = [math.factorial(j) for j in range(TOP + 1)]
CHECK = {"sin2": 0.0, "tan2": 0.0, "n": 0}
C2 = {"pts": 0, "above": 0, "plant": 0, "worst": -math.inf}
# --- the exact test ---
def share6(g):
"""6 times a part's remainder share: odd 3g^2, even g^2 - 4."""
return 3 * g * g if g % 2 else g * g - 4
def prefixes(parts):
"""(S, 6T') over the prefixes of the parts above 2 ranked by share."""
rk = sorted((g for g in parts if g > 2), key=lambda g: (-share6(g), -g))
t6 = sum(share6(g) for g in parts)
for n in range(1, len(rk) + 1):
S = tuple(rk[:n])
yield S, t6 - sum(share6(g) for g in S)
_Q = {}
_E = {}
def qpoly(S):
"""(even-index coefficients, denominator): Q_S = poly(c)/den."""
if S not in _Q:
out, den = [1], 1
for g in S:
if g % 2:
f = exacttest._cheb(exacttest._T, g)
else:
f = exacttest._cheb(exacttest._U, g - 1)[1:]
den *= g * g
out = exacttest.mul(out, exacttest.mul(f, f))
_Q[S] = (out[0::2], den)
return _Q[S]
def eq(S, A, Sq):
"""E_nu Q_S exact in integers, rounded once; nu ~ cos^A sin^Sq."""
key = (S, A, Sq)
if key not in _E:
if len(_E) > 200000:
_E.clear()
co, den = qpoly(S)
H = len(co) - 1
N, D = co[H], 1
for t in range(H - 1, -1, -1):
y = A + Sq + 2 + 2 * t
N = co[t] * D * y + N * (A + 1 + 2 * t)
D *= y
_E[key] = N / (D * den)
return _E[key]
def eq_quad(S, A, Sq):
mp.mp.dps = 30
def q(u):
v = mp.mpf(1)
for g in S:
if g % 2:
v *= mp.cos(g * u) ** 2
else:
v *= (mp.sin(g * u) / ((g // 2) * mp.sin(2 * u))) ** 2
return v
w = lambda u: mp.cos(u) ** A * mp.sin(u) ** Sq # noqa: E731
return mp.quad(lambda u: q(u) * w(u), [0, mp.pi / 4, mp.pi / 2]) \
/ mp.quad(w, [0, mp.pi / 4, mp.pi / 2])
def scan(J, k, s, o, S, t6r, thr, close=False):
"""The scan of the exact LB_S over p, recording D: (best, the
p with LB > thr, cut). close (the re-read): a dead scan also stops,
uncut, once head(p) < thr where head falls for every larger p:
f <= E Q_S <= 1. At q = 0 head falls always (lgamma(be) -
lgamma(al + be) falls in be); at q = 1 its step to p + 2 is
ln[be (p + 2)/((al + be)(p + 1))] at even p and ln[be (p + 3)/
((al + be)(p + 2))] at odd p, both negative when be < al (p + 1),
which then holds at every larger p, be rising by 1 a step and
al (p + 1) by 2 al >= 1."""
q = s % 2
al = (s + q + 1) / 2
os_ = sum(1 for g in S if g % 2)
pmin = math.floor(2 * (t6r / 6 * al + 1) - (J + k - o + 1)) + 1
par = (J + k + os_) % 2
p = max(par, pmin)
p += (p - par) % 2
best, dom = -math.inf, []
since = dead = 0
while True:
be = (J + k - o + p + 1) / 2
f = float(eq(S, J + k - o + p, s + q)) - t6r / 6 * al / (be - 1)
v = head(J, k, s, o, p) + math.log(f) if f > 0 else -math.inf
if v > thr:
dom.append(p)
if v > best:
best, since = v, 0
elif best > -math.inf:
since += 1
else:
dead += 1
if close and head(J, k, s, o, p) < thr and (
q == 0 or (J + k - o + p + 1) / 2 < al * (p + 1)):
return best, dom, False
if dead >= DEADCAP:
return best, dom, True
if since >= 40 and not v > thr:
return best, dom, False
p += 2
def head(J, k, s, o, p):
"""log of LB's factor before f at pi = 1."""
q = s % 2
al = (s + q + 1) / 2
be = (J + k - o + p + 1) / 2
return ((J + k + s + p + q) * LN2 - LNPI + math.lgamma(al)
+ math.lgamma(be) - math.lgamma(al + be)
- compbound.log_w1(p, q))
# --- the arm at order n ---
_X = {}
_NU = {}
def xdist(S):
"""(v >= 0, weight of |X| = v)."""
if S in _X:
return _X[S]
d = {0: 1.0}
for g in S:
if g % 2:
comp = {0: 0.5, g: 0.25, -g: 0.25}
else:
j = g // 2
comp = {2 * m: (j - abs(m)) / (j * j) for m in range(1 - j, j)}
out = {}
for a, pa in d.items():
for b, pb in comp.items():
out[a + b] = out.get(a + b, 0.0) + pa * pb
d = out
top = max(d)
v = np.arange(top + 1, dtype=float)
w = np.array([d.get(0, 0.0)] + [d.get(i, 0.0) + d.get(-i, 0.0)
for i in range(1, top + 1)])
_X[S] = (v, w)
return _X[S]
def nu(A, Sq, nodes=NODES):
"""nu's numbers at each A: rows ubar, TV, TV2, kappa, mo_0..mo_TOP."""
A = np.asarray(A, dtype=float)
gx, gw = np.polynomial.legendre.leggauss(nodes)
al = (Sq + 1) / 2
be = (A + 1) / 2
us = np.arctan(np.sqrt(Sq / A))
kap = (np.sqrt(A) + math.sqrt(Sq)) ** 2
v2 = 2 * (A + Sq) if Sq > 0 else A
half = 16 / np.sqrt(v2)
lo = np.maximum(0.0, us - half)
hi = np.minimum(math.pi / 2, us + half)
x = lo[:, None] + (hi - lo)[:, None] * (gx[None, :] + 1) / 2
ld = A[:, None] * np.log(np.cos(x))
if Sq:
ld = ld + Sq * np.log(np.sin(x))
d = np.exp(ld - ld.max(axis=1, keepdims=True)) * gw[None, :] \
* ((hi - lo) / 2)[:, None]
W = d / d.sum(axis=1, keepdims=True)
ubar = (W * x).sum(axis=1)
dx = x - ubar[:, None]
mos = [np.ones_like(A), np.zeros_like(A)]
pw = dx * dx
for j in range(2, TOP + 1):
mos.append((W * pw).sum(axis=1))
pw = pw * dx
s2 = (W * np.sin(x) ** 2).sum(axis=1)
t2 = (W * np.tan(x) ** 2).sum(axis=1)
CHECK["sin2"] = max(CHECK["sin2"], float(np.max(np.abs(
s2 / (al / (al + be)) - 1))))
CHECK["tan2"] = max(CHECK["tan2"], float(np.max(np.abs(
t2 / (al / (be - 1)) - 1))))
CHECK["n"] += len(A)
lz = gammaln(al) + gammaln(be) - gammaln(al + be) - LN2
def lrho(u):
out = A * np.log(np.cos(u)) - lz
return out + Sq * np.log(np.sin(u)) if Sq else out
rmax = np.exp(lrho(us)) if Sq else np.exp(-lz)
tv = 2 * rmax - (rmax if Sq == 0 else 0)
if Sq >= 2:
a2, b2, c2 = A * A - A, 2 * A * Sq + A + Sq, Sq * Sq - Sq
disc = np.sqrt(b2 * b2 - 4 * a2 * c2)
i1 = np.arctan(np.sqrt((b2 - disc) / (2 * a2)))
i2 = np.arctan(np.sqrt((b2 + disc) / (2 * a2)))
def drho(u):
return np.exp(lrho(u)) * (Sq / np.tan(u) - A * np.tan(u))
tv2 = 2 * (drho(i1) - drho(i2))
else:
tv2 = np.full_like(A, np.inf)
return np.vstack([ubar, tv, tv2, kap] + mos)
def nu_cached(A, Sq):
A = [int(a) for a in A]
miss = [a for a in A if (a, Sq) not in _NU]
if len(_NU) > 300000:
_NU.clear()
miss = A
for c in range(0, len(miss), 1000):
rows = nu(miss[c:c + 1000], Sq)
for i, a in enumerate(miss[c:c + 1000]):
_NU[(a, Sq)] = rows[:, i]
return np.array([_NU[(a, Sq)] for a in A]).T
def arm(S, A, Sq):
"""(bound, bound with the remainder dropped: the PLANT), each of shape
(len ORDERS, len A), the best over TAUS."""
v, w = xdist(S)
rows = nu_cached(A, Sq)
ubar, tv, tv2, kap = rows[:4]
mo = rows[4:]
vs = np.maximum(v, 1.0)[None, :]
bnd = np.minimum(1.0, np.minimum(tv[:, None] / (2 * vs),
tv2[:, None] / (4 * vs * vs)))
bnd[:, 0] = 0.0
ang = 2 * v[None, :] * ubar[:, None]
trig = (np.cos(ang), -np.sin(ang), -np.cos(ang), np.sin(ang))
f = 2 * v[None, :]
x = 4 * v[None, :] ** 2 / kap[:, None]
best = np.full((len(ORDERS), len(ubar)), -np.inf)
bestp = best.copy()
for tau in TAUS:
L = x <= tau
wl = np.where(L, w[None, :], 0.0)
high = (np.where(L, 0.0, w[None, :]) * bnd).sum(axis=1)
terms, mts = [], []
fj = np.ones_like(f)
for j in range(TOP + 1):
if j < TOP:
terms.append((wl * fj * trig[j % 4]).sum(axis=1)
* mo[j] / FACT[j] if j != 1 else 0 * ubar)
mts.append((wl * fj).sum(axis=1) * mo[j] / FACT[j])
fj = fj * f
cum = np.cumsum(np.array(terms), axis=0)
for i, n in enumerate(ORDERS):
trunc = cum[n - 1]
best[i] = np.maximum(best[i], trunc - mts[n] - high)
bestp[i] = np.maximum(bestp[i], trunc - high)
return best, bestp
# --- one rank ---
def judge(J, k, s, o, S, t6r, thr):
"""(exact best, cut, the orders clearing on D(S), arm's best margin
per order as log(f_n / f_needed))."""
best, dom, cut = scan(J, k, s, o, S, t6r, thr)
ok = set()
if not dom:
return best, cut, ok
q = s % 2
al = (s + q + 1) / 2
ps = np.array(dom)
A = J + k - o + ps
be = (A + 1) / 2
val, valp = arm(S, A, s + q)
ex = np.array([float(eq(S, int(a), s + q)) for a in A])
C2["pts"] += val.size
C2["above"] += int(np.sum(val > ex[None, :] * (1 + 1e-9) + 1e-12))
C2["plant"] += int(np.sum(valp > ex[None, :] * (1 + 1e-9) + 1e-12))
C2["worst"] = max(C2["worst"], float(np.max(val - ex[None, :])))
pen = t6r / 6 * al / (be - 1)
need = np.array([math.exp(thr - head(J, k, s, o, int(p)))
for p in dom])
for i, n in enumerate(ORDERS):
if np.any(val[i] - pen > need):
ok.add(n)
return best, cut, ok
def rank(r):
"""(e, |fs|, exact failures [(name, cut)], population, failures per
order [(name, cut)])."""
J, e, lmean, nk, fs = exacttest.failing_shapes(r)
bar = lmean + math.log(1 + 1e-6)
_E.clear()
_Q.clear()
_X.clear()
_NU.clear()
memo = {}
xfail, pop = [], 0
tfail = {n: [] for n in ORDERS}
full = set(ORDERS)
for k, parts, s, o, pi, t6, plb in fs:
thr = bar - math.log(pi)
xok, cut, ok = False, False, set()
for S, t6r in prefixes(parts):
key = (k, s, o, S, t6r, pi)
if key not in memo:
memo[key] = judge(J, k, s, o, S, t6r, thr)
b, c, o2 = memo[key]
xok |= b > thr
cut |= c
ok |= o2
if xok and ok == full:
break
nm = takeph.name(k, parts)
if not xok:
xfail.append((nm, cut))
continue
pop += 1
for n in ORDERS:
if n not in ok:
tfail[n].append((nm, cut))
return e, len(fs), xfail, pop, tfail
# --- a cut re-read ---
def reread(r, name, cap):
"""The shape `name` at r, every prefix's exact scan with DEADCAP =
cap: (best - thr, |D|, cut, the orders clearing) per prefix."""
global DEADCAP
DEADCAP = cap
J, e = takeph.theorem_J(r), takeph.theorem_e(r)
bar = float(circlemean.log_mean(r, J, e)) + math.log(1 + 1e-6)
w = name.split()
k, twos = int(w[0][2:]), int(w[1][2:])
parts = tuple(sorted([int(x) for x in w[2:]] + [2] * twos,
reverse=True))
s, o, pi, t6 = compbound.stats(k, parts)
thr = bar - math.log(pi)
_E.clear()
out, ok = [], set()
for S, t6r in prefixes(parts):
best, dom, cut = scan(J, k, s, o, S, t6r, thr, close=cap > 1000)
o2, marg = set(), []
if dom:
q = s % 2
al = (s + q + 1) / 2
A = J + k - o + np.array(dom)
val, _ = arm(S, A, s + q)
pen = t6r / 6 * al / ((A + 1) / 2 - 1)
need = np.array([math.exp(thr - head(J, k, s, o, int(p)))
for p in dom])
for i, n in enumerate(ORDERS):
m = float(np.max((val[i] - pen) / need))
marg.append("%d:%.3f" % (n, m))
if m > 1:
o2.add(n)
ok |= o2
out.append((S, best - thr, len(dom), cut, " ".join(marg)))
return out, sorted(ok)
# --- the controls ---
def controls():
t0 = time.time()
ok = True
for A, Sq in ((338, 18), (338, 16), (1200, 40)):
a = nu([A], Sq)[4:, 0]
b = nu([A], Sq, 512)[4:, 0]
off = max(abs(a[j] / b[j] - 1) for j in range(2, TOP + 1, 2))
print("C1 A = %d Sq = %d: even moments 2..12 at 256 against 512"
" nodes, worst relative %.1e; mo_12 %.3e" % (A, Sq, off, b[TOP]),
flush=True)
ok &= off < 1e-9
for S, Sq, want in (((9, 7), 18, (-0.0155, 0.0226, 0.0271)),
((11, 7), 16, (-0.0351, 0.0087))):
val, _ = arm(S, np.array([338]), Sq)
ex = float(eq(S, 338, Sq))
got = [val[i, 0] for i in range(len(ORDERS))]
print("C4 S = %s A = 338 Sq = %d: arm by n %s; exact %.4f; the"
" witness %s" % (S, Sq, " ".join("%d:%.4f" % (n, g) for n, g
in zip(ORDERS, got)), ex, want), flush=True)
ok &= all(abs(g - w) < 6e-5 for g, w in zip(got, want))
for S, A, Sq in (((5,), 40, 10), ((27,), 100, 30), ((9, 5), 61, 0),
((8,), 50, 4)):
a = float(eq(S, A, Sq))
b = float(eq_quad(S, A, Sq))
print("C4 Horner %s A = %d Sq = %d: %.12g quadrature %.12g"
% (S, A, Sq, a, b), flush=True)
ok &= abs(a - b) < 1e-10 * max(b, 1e-30)
print("controls %s [%.1f s]" % ("hold" if ok else "FAIL",
time.time() - t0), flush=True)
return ok
def want_exact(r, e, xfail):
"""C5: the exact test's failures, as recorded."""
names = sorted(n for n, _ in xfail)
if e == 2:
return 5 <= len(names) <= 22
if e >= 4 or r <= 27:
return not names
if r == 60:
return names == sorted(["k=0 2^55 5", "k=0 2^51 5 5"])
return names == ["k=0 2^%d 5" % (r - 5)]
C3 = {21: (), 23: (), 25: (), 27: (),
29: ("k=0 2^18 7 5", "k=0 2^16 9 5"),
31: ("k=0 2^18 7 7", "k=0 2^16 9 7"),
33: ("k=0 2^20 7 7", "k=0 2^18 9 7", "k=0 2^16 11 7")}
C3N = {21: 0, 23: 0, 25: 0, 27: 0, 29: 2, 31: 4, 33: 6}
def main():
t0 = time.time()
if "--kill" in sys.argv:
i = sys.argv.index("--kill")
rr = [tuple(int(x) for x in t.split(",")) for t in sys.argv[i + 2:]]
rows = {}
for path in sys.argv[i + 1].split(";"):
rows.update(read_rows(path, rr))
print("rows read: r = %s; cut failures re-read as genuine: %s"
% (sorted(rows), rr))
kills(rows)
return True
if "--reread" in sys.argv:
i = sys.argv.index("--reread")
r, name = int(sys.argv[i + 1]), sys.argv[i + 2]
for cap in (1000, int(sys.argv[i + 3])):
out, ok = reread(r, name, cap)
print("REREAD r = %d %s, DEADCAP %d: clears at n = %s"
% (r, name, cap, ok), flush=True)
for S, gap, nd, cut, o2 in out:
print(" S = %s: exact best - bar %.4f, |D| %d, cut %s;"
" the arm's best f_n/need by n: %s"
% (S, gap, nd, cut, o2 or "-"), flush=True)
print("%.1f s" % (time.time() - t0))
return True
if not controls():
return False
lo = int(sys.argv[sys.argv.index("--from") + 1]) \
if "--from" in sys.argv else 21
hi = int(sys.argv[sys.argv.index("--to") + 1]) \
if "--to" in sys.argv else 33
print()
print(" r e fs xfail(cut) pop | failing by n, genuine+cut"
" | least n | at the largest failing n")
least, res = {}, {}
ok = True
for r in range(lo, hi + 1):
t = time.time()
e, nfs, xfail, pop, tf = rank(r)
res[r] = (e, tf)
good = [n for n in ORDERS if not tf[n]]
least[r] = good[0] if good else math.inf
worst = [n for n in ORDERS if tf[n]]
wit = ", ".join(nm + ("*" if c else "") for nm, c in
tf[worst[-1]][:3]) if worst else "-"
print("%3d %2d %6d %5d(%d) %6d | %s | %s | %s%s [%.1f]"
% (r, e, nfs, len(xfail), sum(c for _, c in xfail), pop,
" ".join("%d:%d+%d" % (n, sum(1 for _, c in tf[n] if not c),
sum(1 for _, c in tf[n] if c))
for n in ORDERS),
least[r], ("n=%d " % worst[-1]) if worst else "", wit,
time.time() - t), flush=True)
if len(xfail) <= 4 and xfail:
print(" exact fails: %s" % ", ".join(n for n, _ in xfail))
if not want_exact(r, e, xfail):
print(" C5 FAIL: the exact test's failures differ from"
" the record's", flush=True)
ok = False
if r in (21, 33):
want = {21: 14, 33: 442}[r]
print(" C5 population %d (want %d)" % (nfs, want))
ok &= nfs == want
if r in C3:
got = [nm for nm, _ in tf[4]]
good3 = len(got) == C3N[r] and all(x in got for x in C3[r])
print(" C3 n = 4: %d failing (want %d) %s: %s"
% (len(got), C3N[r], ", ".join(got),
"holds" if good3 else "FAIL"), flush=True)
ok &= good3
if r == 33 or r == hi:
print(" C1 quadrature at %d (A, Sq): E sin^2 off %.1e,"
" E tan^2 off %.1e; C2 %d points, arm above exact %d,"
" worst arm - exact %.2e; PLANT remainder dropped: above"
" at %d" % (CHECK["n"], CHECK["sin2"], CHECK["tan2"],
C2["pts"], C2["above"], C2["worst"],
C2["plant"]), flush=True)
ok &= CHECK["sin2"] < 1e-10 and CHECK["tan2"] < 1e-10 \
and C2["above"] == 0 and C2["plant"] > 0
if not ok:
print("controls FAIL: stopped before the expensive stage")
return False
print()
kills({r: (res[r][0], {n: (sum(1 for _, c in res[r][1][n] if not c),
sum(1 for _, c in res[r][1][n] if c))
for n in ORDERS}) for r in res})
print("%.1f s" % (time.time() - t0))
return ok
def kills(rows):
"""KILL 1 and KILL 2 from rows r -> (e, {n: (genuine, cut)}): a cut
failure counts only once re-read (passed in as genuine), and is
printed apart."""
k1 = [r for r in sorted(rows) if rows[r][0] >= 3
and all(rows[r][1][n][0] for n in ORDERS)]
cuts = [(r, n) for r in sorted(rows) for n in ORDERS
if rows[r][0] >= 3 and rows[r][1][n][1] and not rows[r][1][n][0]]
print("KILL 1: ranks with e >= 3 failing a genuine shape at every n:"
" %s -> %s" % (k1 or "none", "FIRES" if k1 else "misses"))
def window(rs, genuine=True):
for n in ORDERS:
if all(not rows[r][1][n][0] and (genuine or not rows[r][1][n][1])
for r in rs):
return n
return math.inf
lo, hi = range(21, 34, 2), range(50, 61, 2)
if all(r in rows for r in lo) and all(r in rows for r in hi):
no, ne = window(lo), window(hi)
print("KILL 2: N_odd (odd r = 21..33) = %s, N_even (even r ="
" 50..60) = %s -> %s" % (no, ne, "FIRES" if ne - no >= 4
else "misses"))
print(" with every cut counted: N_odd %s, N_even %s; e >= 3 (r, n)"
" failing by a cut alone: %s"
% (window(lo, False), window(hi, False), cuts))
def read_rows(path, reread=()):
"""The rank rows of a saved run; reread: the (r, n) whose cut
failures were re-read as genuine."""
rows = {}
for line in open(path, encoding="utf-8"):
w = line.split()
if len(w) > 6 and w[0].isdigit() and w[1].isdigit() and "|" in w:
r, e = int(w[0]), int(w[1])
i = w.index("|")
cell = {}
for t in w[i + 1:i + 1 + len(ORDERS)]:
n, gc = t.split(":")
g, c = (int(x) for x in gc.split("+"))
if (r, int(n)) in reread:
g, c = g + c, 0
cell[int(n)] = (g, c)
rows[r] = (e, cell)
return rows
if __name__ == "__main__":
sys.exit(0 if main() else 1)