"""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)