excess.py
Cited on: READING.
To run it: save excess.py, with heat.py, real_support.py, route.py beside it; needs mpmath, numpy.
"""
excess.py -- the larger sup's excess over L past 400, bounded through
D_1 = sum_m |h_m - h_(m+1)|, h's total variation, by the walk's bumps
and the roots.
QUESTION Z. Does an explicit bound D_1^+ >= D_1 leave the excess inside
the carry bound's budget at d = 400, and how does it move with P?
BY HAND.
(Z1) The excess. H_s - L <= (L_s - L) + sup_t eps_t (D_(t+1) + 1 -
L_s)^+ (the real support). D_s <= s D_1, so the term is 0 while
(t + 1) D_1^+ <= L_s - 1, i.e. t < t1 = floor(t0), t0 = (L_s -
1)/D_1^+. At t = t1 + k, D_(t+1) + 1 - L_s <= (k + 1) D_1^+ and
<= 2L + 1 - L_s, so the excess is at most X = (L_s - L) + eps_t1
min(2L + 1 - L_s, (2/beta) D_1^+), max_k (k + 1) beta^-k = 2/beta
for 3/2 <= beta < 2. The coarser X_R = eps_t1 (2L + 1) is printed
beside. H_x's excess vanishes up to L/D_1^+ >= t0, the smaller.
(Z2) The bumps. On m >= 0, h_m = [m = 0]/2 + sum_(u>=1) T_u(m), T_u =
b(u, n) j/(2n), j = m - uP, n = 2u - j, on j <= u, 0 beyond (the
walk). As a function of j, T_u does not involve P: its total
variation V_u is one sequence. By the local limit theorem T_u ~
-phi_n'(j), phi_n the N(0, n) density, n ~ 2u, so u V_u -> kappa =
2/sqrt(2 pi e) = 0.48394, a TRANSPLANT the script reads.
(Z3) The bound. With U = floor(tau_c P^2), M = U P:
D_1 <= |h_-1| + |h_0 - h_-1| + 1/2 + sum_(u<=U) V_u
+ 2 sum_(u>U) sum_(m<=M) |T_u(m)|
+ sum_k wt_k |c_k| |1 - x_k| |x_k|^M/(1 - |x_k|),
h on m <= 0 being -c_0 beta^m, monotone; the 1/2 the jump at m =
0; the bumps by the triangle inequality, those past U by their L1
at m <= M (the pairs with m < M); the roots past M, about 2
sum_k e^(-4 pi^2 k^2 tau_c), 4.64 at 1/400, 0.17 at 1/16. Each
tau_c gives a bound; D_1^+ is the least over the grid.
(Z4) Sizes: D_1 ~ 2 kappa ln P + b, 2 kappa = 0.968, about 6.5 at d =
400; the root tail's triangle inequality loses 1 to 3; t0 ~
98.75/D_1^+, near 12 to 15.
PREDICTIONS, CONTROL AND KILL (frozen before code).
Z0 CONTROL, read first. (a) h from the bumps equals the FFT's
(real_support.inverse_fft) on [0, M] at d = 40, 100, tau_c = 1/16,
to 1e-12; (b) D_1 rebuilt from those pieces (h on m <= 0 by c_0
and beta, the walk to M, the FFT past M) equals the FFT's D_1 to
1e-10, and route.py's streamed D_1 equals the FFT's to 1e-9; (c)
PLANT: j shifted by one in T_u fails (a).
Z1 u V_u within 2% of kappa at u = 1000, |u V_u - kappa| falling over
u = 10, 100, 1000, 10^4.
Z2 D_1 - 0.968 ln P within a band of width 0.6 over d = 101 and
110..400 by 10; D_1 in [5.5, 7.5] at d = 400.
Z3 D_1^+ - D_1 in [0.5, 3.5] at d = 400; the best tau_c printed.
Z4 at d = 400: X_R in [3e-3, 0.1] (the planned bar, X_R below 1e-3,
predicted to miss: a transplant of route.py's X1a, itself missed;
the budget is what the theorem needs), X in [1e-4, 5e-3].
Z5 X falls with d over the d = 0 mod 4 read.
KILL: (K1) D_1^+ below a printed D_1 at any d read (those of Z2 and
800); or (K2) R_maj + F + X at d = 400 at or above 1/4 + 1/(4P) -
1/P^2, R_maj + F = 0.063511 (uniform.py). Past 400, X is read at P
= 799, 1599, 3199, 9999 with V_u past its computed range as
kappa'/u, kappa' the largest u V_u on its last decade: printed.
FINDINGS (copied from the second run's prints).
- Z0 holds: the bumps' h equals the FFT's on [0, M] to 4.6e-16 (d =
40) and 1.1e-15 (d = 100); D_1 rebuilt to 1.2e-14 and 5.1e-14,
route.py's streamed D_1 to 1.5e-13 and 3.3e-13; the shifted plant
misses by 0.25.
- Z1 holds: u V_u is 0.625, 0.500314, 0.486121, 0.484145, 0.483959,
0.483947 at u = 1, 10, 100, 1000, 10^4, 4e4 (0.65625 at u = 2),
falling onto kappa = 0.483941 from above past u = 2; the largest on
4000..4e4 is 0.483996.
- Z2 holds, tighter than frozen: D_1 - 0.968 ln P lies in [0.7780,
0.7790] at every d read (0.7778 at d = 800), D_1 = 6.575283 at d =
400, 7.247336 at 800: D_1 ~ 0.968 ln P + 0.778 (2 kappa = 0.96788),
an observation.
- Z3 fails on the tight side: D_1^+ - D_1 is 0.1382 at d = 400, not
0.5 to 3.5; the least tau_c is 1/16 at every d, the root tail there
0.1698 (by hand 0.17), the far bumps 1.1e-4. The by-hand loss of 1
to 3 was the root tail at small tau_c, which the grid does not pick.
- Z4 holds: at d = 400, t0 = 14.711, X_R = 6.119e-3 (the planned 1e-3
bar missed, as predicted), X = 2.049e-4 (1.003e-4 with D_1 itself).
- Z5 fails on its premise: X steps down as t1 = floor(t0) steps and
rises slightly between steps (eps_t1 fixed, D_1^+ growing), 8.66e-2
at d = 120, 2.05e-4 at 400; it falls only in envelope.
- The KILL does not fire: D_1^+ >= D_1 at every d read; R_maj + F + X
= 0.063716 at d = 400 against 0.250620.
- Past 400, V_u past 4e4 as 0.483996/u: D_1^+ - 0.968 ln P = 0.9159
at d = 800 (U = 39900 inside the table, so the exact read), 1600,
3200 and 0.9158 at 10^4; X = 5.5e-8, 7.2e-15, 1.8e-27, 1.7e-76.
RUN RECORD: 10/12 checks, the misses Z3 and Z5, as recorded; 116.2 s
under the 512 MB guard, peak commit 232.4 MB. The first run stopped in
its last section: heat.roots' Newton step for beta overflowed b^d at d
= 1600; heat.roots now takes beta by its fixed point. Each bump
is read on a fixed window, its ends asserted below 1e-30 (the top only
where the window stops short of j = u), monotone beyond; the lower end
below j_lo is unpriced in far(), 1.2e-38 at u = 1 and 2.7e-83 at u =
1000 by a call of bump() after the run.
"""
import math
import os
import sys
import time
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np # noqa: E402
from heat import roots, coef # noqa: E402
from real_support import inverse_fft # noqa: E402
from route import stream # noqa: E402
CHECKS, NAMES = [], []
MISSED = ("Z3 D_1^+ - D_1 in [0.5, 3.5] at d = 400",
"Z5 X falls with d over the d = 0 mod 4 read")
KAPPA = 2 / math.sqrt(2 * math.pi * math.e)
TAUS = (1 / 400, 1 / 200, 1 / 100, 1 / 64, 1 / 48, 1 / 32, 1 / 24, 1 / 16,
1 / 12, 1 / 8)
UMAX = 40000
RMAJ_F = 0.063511 # uniform.py, P = 399
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
NAMES.append(name)
print(f" [{'PASS' if ok else 'FAIL'}] {name}" +
(f" -- {detail}" if detail else ""), flush=True)
def bump(u):
"""(j, T_u(j)) over j in [j_lo, j_hi], j_hi = min(u, 30 sqrt u + 100),
j_lo = -(30 sqrt u + 100): log C(n, u) by lgamma at the top, then
C(n + 1, u) = C(n, u)(n + 1)/(n + 1 - u) as j falls."""
w = int(math.ceil(30 * math.sqrt(u))) + 100
j_hi = min(u, w)
j = np.arange(j_hi, -w - 1, -1)
n = 2 * u - j
n0 = int(n[0])
lb0 = math.lgamma(n0 + 1) - math.lgamma(u + 1) - math.lgamma(n0 - u + 1) \
- n0 * math.log(2)
step = np.log((n[1:]) / (2.0 * (n[1:] - u)))
lb = lb0 + np.concatenate([[0.0], np.cumsum(step)])
return j, np.exp(lb) * j / (2.0 * n), j_hi
def V_table(umax):
"""V_u, u = 1..umax: the total variation of T_u over j, its end jumps
included (both below 1e-30 but at j = u, asserted)."""
V = np.zeros(umax + 1)
for u in range(1, umax + 1):
j, T, j_hi = bump(u)
end_lo = abs(T[-1])
assert end_lo < 1e-30, (u, end_lo)
if j_hi < u:
assert abs(T[0]) < 1e-30, (u, T[0])
V[u] = abs(T[0]) + np.abs(np.diff(T)).sum() + end_lo
return V
def walk_h(P, M, shift=0):
"""h_m on m = 0..M by the bumps (shift moves each bump's j: the plant)."""
h = np.zeros(M + 1)
h[0] += 0.5
u = 1
while True:
j, T, _ = bump(u)
m = u * P + j + shift
if m.min() > M:
break
keep = (m >= 0) & (m <= M)
np.add.at(h, m[keep], T[keep])
u += 1
return h
def far(P, U):
"""2 sum_(u > U) sum_(m <= M) |T_u(m)|, M = U P."""
M, s, u = U * P, 0.0, U + 1
while True:
j, T, _ = bump(u)
m = u * P + j
if m.min() > M:
break
s += np.abs(T[m <= M]).sum()
u += 1
return 2 * s
def root_tail(d, M, rt=None):
x, wt, b, _ = rt or roots(d)
c = coef(x, d)
ax = np.abs(x)
return float((wt * np.abs(c) * np.abs(1 - x) * np.exp(M * np.log(ax)) /
(1 - ax)).sum())
def edge(d, rt=None):
"""|h_-1| + |h_0 - h_-1| + 1/2, and beta, c_0."""
x, wt, b, _ = rt or roots(d)
c0 = float(coef(b, d))
hm1, h0 = -c0 / b, 1 - c0
return abs(hm1) + abs(h0 - hm1) + 0.5, b, c0
def D1_bound(d, Vc, kappa_x=None):
"""min over tau_c of the bound (Z3), with its pieces; V_u past the
table by kappa_x/u when kappa_x is given."""
P = d - 1
rt = roots(d)
A, b, c0 = edge(d, rt)
best = None
for tc in TAUS:
U = int(tc * P * P)
if U > UMAX:
if kappa_x is None:
continue
bumps = Vc[UMAX] + kappa_x * math.log(U / UMAX) # >= sum 1/u
else:
bumps = Vc[U]
fa = far(P, U)
rtl = root_tail(d, U * P, rt)
tot = A + bumps + fa + rtl
if best is None or tot < best[0]:
best = (tot, tc, A, bumps, fa, rtl)
return best
def excess(d, L, h0, D1p, b):
Ls = L - abs(h0) + abs(h0 - 1)
B = sum(b ** -s for s in range(d))
t0 = (Ls - 1) / D1p
t1 = math.floor(t0)
eps = b ** -t1 / B
XR = eps * (2 * L + 1)
X = (Ls - L) + eps * min(2 * L + 1 - Ls, (2 / b) * D1p)
return X, XR, t0, t1
def section_control():
print("Z0 CONTROL", flush=True)
for d in (40, 100):
P = d - 1
U = int(P * P / 16)
M = U * P
hf, hb, _ = inverse_fft(d)
hw = walk_h(P, M)
err = np.abs(hw - hf[:M + 1]).max()
check(f"Z0a the bumps' h equals the FFT's on [0, M] at d = {d}",
err < 1e-12, f"max {err:.1e}, M = {M}")
full = np.concatenate([hb[:0:-1], hf, [0.0]])
D1_fft = np.abs(np.diff(full)).sum()
A, b, c0 = edge(d)
rest = np.abs(np.diff(np.concatenate([hf[M:], [0.0]]))).sum()
D1_re = A - 0.5 + np.abs(np.diff(hw)).sum() + rest
L, E, h0s, D, tail, bs, integ = stream(d)
check(f"Z0b D_1 rebuilt and route.py's streamed D_1 equal the FFT's "
f"at d = {d}",
abs(D1_re - D1_fft) < 1e-10 and abs(D[1] - D1_fft) < 1e-9,
f"FFT {D1_fft:.12f}, rebuilt {D1_re - D1_fft:+.1e}, "
f"streamed {D[1] - D1_fft:+.1e}")
if d == 40:
hp = walk_h(P, M, shift=1)
errp = np.abs(hp - hf[:M + 1]).max()
check("Z0c PLANT: j shifted by one fails Z0a at d = 40",
errp >= 1e-12, f"max {errp:.1e}")
def section_V():
print("Z1 V_u, the bumps' total variation", flush=True)
t0 = time.time()
V = V_table(UMAX)
Vc = np.cumsum(V)
devs = []
for u in (1, 2, 5, 10, 100, 1000, 10000, 40000):
print(f" u = {u}: V_u {V[u]:.6e}, u V_u {u * V[u]:.6f}, "
f"sum to u {Vc[u]:.6f}", flush=True)
if u in (10, 100, 1000, 10000):
devs.append(abs(u * V[u] - KAPPA))
kx = float((np.arange(4000, UMAX + 1) * V[4000:]).max())
print(f" kappa {KAPPA:.6f}; largest u V_u on u = 4000..{UMAX}: {kx:.6f}"
f" ({time.time() - t0:.0f} s)", flush=True)
check("Z1 u V_u within 2% of kappa at u = 1000, the gap falling",
abs(1000 * V[1000] / KAPPA - 1) <= 0.02 and
all(b_ < a_ for a_, b_ in zip(devs, devs[1:])),
" ".join(f"{v:.2e}" for v in devs))
return Vc, kx
def section_read(Vc):
print("Z2-Z5 and KILL: D_1, D_1^+ and the excess at d <= 400 and 800",
flush=True)
ds = [101] + list(range(110, 401, 10)) + [800]
rows = {}
for d in ds:
t0 = time.time()
P = d - 1
L, E, h0, D, tail, b, integ = stream(d)
L = L + tail
D1 = D[1] + 2 * tail
tot, tc, A, bumps, fa, rtl = D1_bound(d, Vc)
X, XR, t0_, t1 = excess(d, L, h0, tot, b)
Xa, _, _, _ = excess(d, L, h0, D1, b)
rows[d] = (D1, tot, X, XR, Xa, L, b)
print(f" d = {d}: D_1 {D1:.6f} (- 0.968 ln P {D1 - 0.968 * math.log(P):.4f}),"
f" D_1^+ {tot:.6f} at tau_c 1/{1 / tc:.0f} (edge {A:.4f}, "
f"bumps {bumps:.4f}, far {fa:.1e}, roots {rtl:.4f}); t0 {t0_:.3f},"
f" X {X:.3e}, X_R {XR:.3e}, X by D_1 itself {Xa:.3e}; beta "
f"{b:.12f} ({time.time() - t0:.1f} s)", flush=True)
low = [d for d in ds if rows[d][1] < rows[d][0]]
check("KILL K1 not fired: D_1^+ >= D_1 at every d read", not low,
f"below at {low}")
sp = [rows[d][0] - 0.968 * math.log(d - 1) for d in ds if d <= 400]
check("Z2 D_1 - 0.968 ln P in a band of width 0.6 over d <= 400, D_1 in "
"[5.5, 7.5] at 400", max(sp) - min(sp) <= 0.6 and
5.5 <= rows[400][0] <= 7.5,
f"band [{min(sp):.4f}, {max(sp):.4f}], D_1(400) {rows[400][0]:.4f}")
g = rows[400][1] - rows[400][0]
check("Z3 D_1^+ - D_1 in [0.5, 3.5] at d = 400", 0.5 <= g <= 3.5,
f"{g:.4f}")
X, XR = rows[400][2], rows[400][3]
check("Z4 at d = 400 X_R in [3e-3, 0.1] and X in [1e-4, 5e-3]",
3e-3 <= XR <= 0.1 and 1e-4 <= X <= 5e-3,
f"X_R {XR:.3e} (1e-3 bar: {'met' if XR < 1e-3 else 'missed'}), "
f"X {X:.3e}")
q = [rows[d][2] for d in ds if d <= 400 and d % 4 == 0]
check("Z5 X falls with d over the d = 0 mod 4 read",
all(b_ < a_ for a_, b_ in zip(q, q[1:])),
" ".join(f"{v:.2e}" for v in q))
P = 399
budget = 1 / 4 + 1 / (4 * P) - 1 / P ** 2
check("KILL K2 not fired: R_maj + F + X below the budget at d = 400",
RMAJ_F + X < budget, f"{RMAJ_F + X:.6f} against {budget:.6f}")
return rows
def section_past(Vc, kx):
print("PAST 400 (V_u past the table as kappa'/u): an observation",
flush=True)
for d in (800, 1600, 3200, 10000):
P = d - 1
tot, tc, A, bumps, fa, rtl = D1_bound(d, Vc, kx)
L = (P * P - 1) / (4 * P) + 1 / P ** 2 # W at odd P; E > 0 left out
rt = roots(d)
h0 = 1 - float(coef(rt[2], d))
X, XR, t0_, t1 = excess(d, L, h0, tot, rt[2])
print(f" d = {d}: D_1^+ {tot:.4f} at tau_c 1/{1 / tc:.0f} (- 0.968 ln"
f" P {tot - 0.968 * math.log(P):.4f}); t0 {t0_:.2f}, X {X:.3e}, "
f"X_R {XR:.3e}", flush=True)
def main():
t0 = time.time()
section_control()
if not all(CHECKS):
print("control failed; stopping")
sys.exit(1)
Vc, kx = section_V()
section_read(Vc)
section_past(Vc, kx)
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed, {time.time() - t0:.0f} s")
sys.exit(as_recorded(MISSED))
def as_recorded(missed):
"""0 iff the checks that failed are exactly the record's misses, by name."""
got = {n for n, ok in zip(NAMES, CHECKS) if not ok}
print("misses %s; the record's %s: %s" % (
"; ".join(sorted(got)) or "none", "; ".join(sorted(missed)) or "none",
"as recorded" if got == set(missed) else "NOT AS RECORDED"))
return 0 if got == set(missed) else 1
if __name__ == "__main__":
main()