edge.py
Checks: READING: The walk theorem.
Cited on: READING.
To run it: save edge.py, with heat.py, real_support.py beside it; needs mpmath, numpy, scipy.
"""
edge.py -- E at small tau by the walk's binomial series, and the
remainder the leading heat term leaves past it.
QUESTION Y. Is E's share at small tau, m < m1 = P^3/c, bounded
rigorously and exponentially small; and does the leading heat term phi =
-P^-2 K_x(r/P, m/P^3), as an oracle (h known), leave a remainder below
1/4 on m >= m1 at d = 400?
BY HAND.
(Y1) G_m = sum_n Pr(S_n = m), S the walk with steps -1, +d (each 1/2).
u up-steps in n steps reach u(d + 1) - n, so n = u(P + 2) - m,
down-steps u(P + 1) - m >= 0, and G_m = sum_u b(u, n_u), b(u, n) =
C(n, u) 2^-n. As C(n - 1, u) = C(n, u)(n - u)/n, b(u, n) - b(u, n -
1) = b(u, n)(2u - n)/n, and 2u - n = m - uP, so with h_m = (G_m -
G_(m+1))/2,
h_m = [m = 0]/2 + sum_(u >= 0, n_u >= max(u, 1)) T_u(m),
T_u(m) = b(u, n_u) (m - uP)/(2 n_u).
Read: h_0 = 1/2 - (P/(2(P + 2))) sum_(u>=1) b(u, u(P + 2)), so
1 - 2 h_0 ~ P 2^-(P+2), the u = 1 term; m < 0 has the u = 0 term
-2^(m-1), against -c_0 beta^m.
(Y2) T_u(m) has the sign of j = m - uP, j = m (mod P). The sawtooth's
sign s is + on r < (P - 1)/2, - on r > (P - 1)/2 and on r = (P -
1)/2 at odd P, so T_u is right-signed at j in [-(P + 1)/2, (P -
1)/2), and wrong only where |j| >= (P - 1)/2. Per m, |h| - s h <=
sum_u (|T_u| - s T_u), so E over any m-range is at most
Ebar = sum over that range and the u with T_u wrong of b |j|/n.
Hoeffding: b(u, n) <= exp(-j^2/(2n)) (one-sided tail at |u - n/2|
= |j|/2), so Ebar <= Ebar_H = the same with b replaced.
At m < m1 = P^3/c, u <= P^2/c, n ~ 2u: the exponent is at least
((P - 1)/2)^2/(4P^2/c) ~ c/16, over ~P^2/c bumps; c = 400 puts
Ebar near 400 e^-25 ~ 1e-8 at P = 399, uniformly c/16 in P.
(Y3) Oracle remainder on m >= m1: |h| - s h <= (|phi| - s phi) +
2|h - phi| at every m, zero where s h > 0, so
E_(>=m1) <= R = sum over Bad of (|phi| - s phi) + 2|h - phi|,
Bad = {m >= m1: s phi <= |h - phi|}. TRANSPLANT (heat.py Z4):
h - phi ~ P^-3 K_xx z(tau) near each zero x0, phi ~ -P^-3 K_xx
rho, rho classes from x0, so Bad is |rho| <= |z| and R ~ (2/P)
sum_(x0) int |K_xx| |z| (2|z| + 1) d tau: P R of order 10 to 50,
flat in P; z -> 0 at small tau, so tau0 = 1/c adds little.
PREDICTIONS, CONTROL AND KILL (fixed before the code).
Y0 CONTROL, read first. (a) The binomial series equals the root sum's
h to 1e-12 on m in [-80, 40P] at d = 13, 40, 100. (b) At c = 4
(tau0 = 1/4) E_(<m1)/E at d = 400 is 1 - 0.0045 within 0.001
(heat.py's share of [1/4, inf) at d = 400). (c) Ebar >= E_(<m1) and
R >= E_(>=m1) at every d read.
Y1 THE EDGE: at c = 400, Ebar <= 1e-6 and Ebar_H <= 1e-6 at d = 100,
200, 400.
Y2 THE REMAINDER: P R between 10 and 50 at d = 100, 200, 400.
KILL of the leading phi: R >= 1/4 at d = 400 (P = 399, d = 0 mod 4,
margin 1/4 + 1/(4P) - 1/P^2 less E and the excess); phi must then
carry the first order.
FINDINGS (copied from the run's prints).
- Y0 holds. (a) The binomial series equals the root sum on m in [-80,
40P] to 3.3e-16, 4.4e-16, 9.2e-16 at d = 13, 40, 100; 1 - 2 h_0 =
7.336e-4 against P 2^-(P+2) = 7.324e-4 at d = 13, 1.774e-11 against
1.774e-11 at d = 40. (b) At c = 4, E's share below m1 at d = 400 is
0.99545. (c) Ebar >= E below m1 and R >= E past m1 at every d read;
only Ebar's half at c = 400 can fail, R's holding term by term and
Ebar at c = 4 being 125 against E's 0.0118. The window W drops terms
with j^2/(2n) above about 64 (by hand), so truncation cannot move
Ebar at the printed digits.
- Ebar is tight: E below m1 against Ebar at c = 400 is 7.070e-7 /
7.070e-7 (d = 100), 1.503e-8 / 1.504e-8 (200), 9.835e-10 / 9.931e-10
(400): the edge's wrong-signed mass is single bumps' tails.
- Y1 fails as frozen at d = 100 only: Ebar_H is 1.68e-5 there (Ebar
7.07e-7 holds); 4.61e-7 and 4.48e-8 at d = 200, 400. Hoeffding drops
the local factor of b, about (pi n/2)^(-1/2), and the bump count
P^2/c multiplies e^(-c/16); the bar is met from d = 200.
- Y2 holds, and the KILL does not fire: P R = 27.29, 24.16, 19.87 at d =
100, 200, 400, falling, not flat; R = 0.0498 at d = 400, a fifth of
1/4. R's shares over [tau0, 1/16), [1/16, 1/4), [1/4, 1): 0.746 0.252
0.002 (d = 100), 0.714 0.283 0.003 (200), 0.655 0.342 0.004 (400).
RUN RECORD: 14/15 checks, the miss Y1 at d = 100, which the exit answers
for (0 iff the misses are the record's); 32.3 s, peak 218.3 MB.
"""
import math
import os
import sys
import time
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np # noqa: E402
from scipy.special import gammaln # noqa: E402
from heat import roots, coef # noqa: E402
CHECKS, NAMES = [], []
MISSED = ("Y1",) # the record's one miss
LN2 = math.log(2)
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 signs(P):
r = np.arange(P)
return np.sign(((P - 1) / 2 - r) / P - 1 / P ** 2)
def logb(u, n):
"""log of C(n, u) 2^-n."""
return gammaln(n + 1) - gammaln(u + 1) - gammaln(n - u + 1) - n * LN2
def h_binom(d, ms):
"""h_m = [m = 0]/2 + sum_u b(u, n)(m - uP)/(2n), n = u(P + 2) - m
(Y1), at each m of ms."""
P = d - 1
out = []
for m in ms:
u0 = max(0, -(-m // (P + 1)))
u = np.arange(u0, max(u0, m // P) + 12, dtype=float)
n = u * (P + 2) - m
ok = (n >= np.maximum(u, 1))
u, n = u[ok], n[ok]
t = np.exp(logb(u, n)) * (m - u * P) / (2 * n)
out.append(t.sum() + (0.5 if m == 0 else 0.0))
return np.array(out)
def ebar(d, m1, mlo):
"""Ebar and Ebar_H (Y2) over m in [mlo, m1): sum over the wrong-signed
T_u of b |j|/n, and the same with b <= exp(-j^2/(2n))."""
P = d - 1
s = signs(P)
tot = totH = 0.0
u = 0
while True:
W = int(130 + 16 * math.sqrt(u + 1))
jlo = max(mlo - u * P, -W)
jhi = min(u, m1 - 1 - u * P, W)
if u * P - W >= m1:
break
if jlo <= jhi:
j = np.arange(jlo, jhi + 1)
j = j[j != 0]
n = (2 * u - j).astype(float)
j, n = j[n >= max(u, 1)], n[n >= max(u, 1)]
wrong = s[j % P] != np.sign(j)
j, n = j[wrong], n[wrong]
if len(j):
aj = np.abs(j)
tot += (np.exp(logb(float(u), n)) * aj / n).sum()
totH += (np.exp(-j.astype(float) ** 2 / (2 * n)) * aj / n).sum()
u += 1
return tot, totH
def kx(x, tau):
"""K_x of the periodic heat kernel (K_tau = K_xx, period 1), by images
below tau = 0.05 and by Fourier above."""
out = np.zeros_like(tau)
lo = tau < 0.05
if lo.any():
t, xx = tau[lo], x[lo]
acc = np.zeros_like(t)
for k in range(-3, 4):
y = xx - k
acc += (4 * np.pi * t) ** -0.5 * np.exp(-y * y / (4 * t)) \
* (-y / (2 * t))
out[lo] = acc
hi = ~lo
if hi.any():
t, xx = tau[hi], x[hi]
q = np.exp(-4 * np.pi ** 2 * t)
qk, q2k = q.copy(), q ** 3 # q^(k^2), q^(2k + 1)
acc = np.zeros_like(t)
for k in range(1, 11):
acc += -4 * np.pi * k * np.sin(2 * np.pi * k * xx) * qk
qk = qk * q2k
q2k = q2k * q * q
out[hi] = acc
return out
def stream(d, m1):
"""E below m1 and at or past it, and the oracle remainder R past m1
(Y3) with its tau bands, h by the root sum over [-80 d, P^3)."""
P = d - 1
x, wt, b, resid = roots(d)
c = coef(x, d) * wt
lx = np.log(x)
absx = np.abs(x)
s = signs(P)
out = dict(Elt=0.0, Ege=0.0, R=0.0, Rb=np.zeros(3))
edges = np.array([1 / 16, 1 / 4])
def take(m, h):
e = np.abs(h) - s[m % P] * h
lt = m < m1
out["Elt"] += e[lt].sum()
out["Ege"] += e[~lt].sum()
if (~lt).any():
mm, hh = m[~lt], h[~lt]
ph = -kx((mm % P) / P, mm / P ** 3) / P ** 2
sm = s[mm % P]
bad = sm * ph <= np.abs(hh - ph)
rr = (np.abs(ph) - sm * ph + 2 * np.abs(hh - ph))[bad]
out["R"] += rr.sum()
band = np.searchsorted(edges, mm[bad] / P ** 3, side="right")
out["Rb"] += np.bincount(band, weights=rr, minlength=3)
mneg = np.arange(-80 * d, 0)
take(mneg, -coef(b, d) * b ** mneg.astype(float))
M = P ** 3
order = np.argsort(-absx)
m0 = 0
while m0 < M:
act = order[absx[order] ** m0 > 1e-20]
R = max(len(act), 1)
n = int(min(max(m0 // 8, 4096), 2 ** 20, 2 ** 21 // min(R, 64),
M - m0))
j = np.arange(n)
base = c[act] * np.exp(m0 * lx[act])
h = np.zeros(n)
for a in range(0, len(act), 64):
sl = slice(a, a + 64)
h += (np.exp(np.outer(j, lx[act][sl])) @ base[sl]).real
take(m0 + j, h)
m0 += n
ac = np.abs(c)
tail = (ac * absx ** M / (1 - absx)).sum() + (ac * 1e-20 /
(1 - absx)).sum()
out["tail"] = tail
return out
def section_control():
print("Y0 CONTROL", flush=True)
for d in (13, 40, 100):
P = d - 1
ms = np.arange(-80, 40 * P + 1)
hb = h_binom(d, ms)
x, wt, b, _ = roots(d)
c = coef(x, d) * wt
pos = ms[ms >= 0]
hr = np.concatenate([-coef(b, d) * b ** ms[ms < 0].astype(float),
(np.exp(np.outer(pos, np.log(x))) @ c).real])
err = np.abs(hb - hr).max()
check(f"Y0a binomial series = root sum on m in [-80, 40P] at d = {d}",
err < 1e-12, f"max difference {err:.1e}; 1 - 2 h_0 = "
f"{1 - 2 * hb[80]:.3e} against P 2^-(P+2) = "
f"{P * 2.0 ** -(P + 2):.3e}")
def section_edge():
print("Y0b, Y0c, Y1, Y2 at d = 100, 200, 400", flush=True)
res = {}
for d in (100, 200, 400):
t0 = time.time()
P = d - 1
for cc in ((4, 400) if d == 400 else (400,)):
m1 = P ** 3 // cc
o = stream(d, m1)
eb, ebH = ebar(d, m1, -80 * d)
E = o["Elt"] + o["Ege"]
R = o["R"] + 2 * o["tail"]
res[(d, cc)] = (o, eb, ebH, R)
print(f" d = {d}, c = {cc}: P E {P * E:.5f}, E below m1 "
f"{o['Elt']:.3e} (share {o['Elt'] / E:.5f}), Ebar {eb:.3e},"
f" Ebar_H {ebH:.3e}; E past m1 {o['Ege']:.4e}, R {R:.4e},"
f" P R {P * R:.3f}, R's bands [tau0, 1/16) [1/16, 1/4) "
f"[1/4, 1) {' '.join(f'{v / o['R']:.4f}' for v in o['Rb'])}"
f"; tail {o['tail']:.1e} ({time.time() - t0:.1f} s)",
flush=True)
check(f"Y0c Ebar >= E below m1 and R >= E past m1, d = {d}, "
f"c = {cc}", eb >= o["Elt"] - 1e-15 and R >= o["Ege"])
o, eb, ebH, R = res[(400, 4)]
sh = o["Elt"] / (o["Elt"] + o["Ege"])
check("Y0b at c = 4, E's share below m1 at d = 400 is 1 - 0.0045 "
"within 0.001", abs(sh - (1 - 0.0045)) <= 0.001, f"{sh:.5f}")
for d in (100, 200, 400):
o, eb, ebH, R = res[(d, 400)]
check(f"Y1 Ebar, Ebar_H <= 1e-6 at c = 400, d = {d}",
eb <= 1e-6 and ebH <= 1e-6, f"{eb:.2e}, {ebH:.2e}")
for d in (100, 200, 400):
o, eb, ebH, R = res[(d, 400)]
check(f"Y2 P R in [10, 50] at d = {d}", 10 <= (d - 1) * R <= 50,
f"P R {(d - 1) * R:.3f}")
o, eb, ebH, R = res[(400, 400)]
check("KILL not fired: R < 1/4 at d = 400", R < 0.25, f"R {R:.4f}")
def main():
section_control()
if not all(CHECKS):
print("control failed; stopping")
sys.exit(1)
section_edge()
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed")
sys.exit(as_recorded(MISSED))
def as_recorded(missed):
"""0 iff the checks that failed are exactly the record's misses."""
got = {n.split()[0] 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()