bandheight.py
Checks: FLATTEN: Past the last failure the champion attains the least height rule; rule; property; pattern.
Cited on: FLATTEN.
To run it: save bandheight.py, with bandend.py, flatten.py beside it; needs numpy, scipy, sympy.
"""bandheight.py -- does the Gaussian limit of a class's least height
predict the last failing depth J_hi(r) of each rank's band, and its growth?
QUESTION. J_hi = 30, 33, 58, 60, 88, 106, 122, 152, 168, 195 at ranks
5..14 (gap.py, bandgrowth.py) is the last depth where some vector off the
champion's line beats the champion's height H_ch (J_ch = J_hi + 1 at every
rank read). bandscale.py's limit of the dual bound gives J_c; does the same
limit of the class's own least height give J_hi?
THE LIMIT (by hand, from bandscale.py's). Write the cofactor of class s as
(1 - x)^s c', c' = sum_k e_k (1 - x)^k, e_0 = c'(1) = 1. The coefficients of
(1 + x)^J (1 - x)^n tend to 2^J sigma^(-1-n) He_n phi / sqrt(2 pi), sigma =
sqrt(J)/2, phi = e^(-y^2/2), so b_k = e_k sigma^(-k) and the class's real
least height over H_ch tends to sigma^d mu_s / M_(r-1), d = r - 1 - s,
mu_s = min over real b of max |(He_s + sum_(m<=d) b_m He_(s+m)) phi|, M_n =
max |He_n phi|. Class s beats the champion below J_hi*_s = 4 (M_(r-1) /
mu_s)^(2/d); J_hi* is the largest over s <= r - 2.
Parity (proved): a term He_(s+m) with m odd has the parity opposite to
He_s, and max(|E + O|, |E - O|) = |E| + |O| >= |E| at y and -y, so mu_s
is a problem in the even m only; at d = 1, mu_s = M_s.
Holder: kappa*_s <= mu_s, so J_hi* <= J* is forced and reads nothing.
The lattice: the integer least height is at least the real one, so the
exact real relaxation's last crossing J_hi^R >= J_hi is forced.
THE CONSTANT (by hand). At d = 2, in the leading form of bandscale.py's
THE GROWTH, f ~ A_s delta y sin(w y - s pi/2) e^(-y^2/4), both minimisers
alike, the oscillation reaching its envelope: mu ~ A_s delta sqrt(2)
e^(-1/2), kappa ~ A_s delta pi sqrt(2 pi)/16, so mu/kappa -> 1.7428 and
J_hi* ~ 0.5738 J* ~ 4.663 r^(3/2). At d = 4 the envelope He_2(y) e^(-y^2/4)
gives a constant near 3.73, so d = 2 binds asymptotically. TRANSPLANT: the
oscillation reaching its envelope needs sqrt(s) large; at s ~ 5 it is not.
THE ARMS.
(G) mu_s by an exchange LP on a fine grid, the even m only, every
s <= r - 2, ranks 3..40: J_hi*_s, J_hi*, the binding order.
(R) the exact real relaxation at finite J: per class an LP over e_1..e_d
on the J + r coefficients of the columns (1 + x)^J (1 - x)^(s+k), exact
in integers then scaled; J_hi^R the first depth, scanning down from
J_c - 1, where some class beats H_ch, its beat rechecked exactly in
rationals at the LP's e. Ranks 5..14.
(A) mu_s at d = 1..4 to rank 801 from He_n/sqrt(n!), factorials in logs,
by golden-section search on the even parameters.
PREDICTIONS (frozen before the engine).
P1 J_hi*/J_hi in [0.9, 1.2] at every rank 5..14 (TRANSPLANT from J*/J_c =
1.019..1.059; the lattice pushes J_hi below the real crossing).
P2 J_hi^R/J_hi in [1.0, 1.15] and J_hi*/J_hi^R in [0.9, 1.1], ranks 5..14.
P3 the binding order of J_hi* is d = 2 (s = r - 3) at ranks 5..14.
P4 mu/kappa at d = 2 within 3% of 1.7428 at r = 800; J_hi*(800)/800^1.5
within 5% of 4.663.
KILLS. K1 J_hi*/J_hi outside [0.8, 1.25] at some rank 5..14. K2
J_hi*(800)/800^1.5 outside [4.2, 5.1], or log2(J_hi*(800)/J_hi*(400))
outside [1.45, 1.55]. P1..P4 read, not gated.
CONTROLS, read first. C1 the LP over every m (odd included) at d <= 4
equals the even-only LP to relative 1e-6, ranks 3..40 (at d = 1 that is
mu_s = M_s). C2 at rank 5, J = 30, 36, ..., 60, the (R) LP's least height at d = 1, 2
against a golden-section search (nested at d = 2: the parity is the
limit's, not exact at finite J, so both parameters are searched), relative
1e-6. C3 (R) over its limit at J = 16 r^2, s = r - 3, ranks 5..14, within
2% of 1. C4 (A) against (G) at r = 5..40, d = 1..4, relative 1e-4. C5
J_hi^R >= J_hi at every rank. C6 mu_s >= kappa*_s at every s read.
INSTRUMENT, before the counted run. A first run was stopped past twice its
estimate in (G), and the grid cut from 120001 to 30001 points (the cost is
the LP, not the grid). A second read C1 4.81e-06 and C6 False: kappa by
the normal equations loses the projection at large d, now by least
squares; HiGHS's absolute tolerances on rows of size ~1e-2, now each
column scaled to max 1 (HiGHS here returned no solution at tolerances
1e-9 and 1e-10). Its J_hi*, (R) and (A) rows equal the counted run's;
(G) took 227.7 s of the counted run.
FINDINGS (copied from the printed output).
(1) Controls: C1 worst 2.56e-07; C2 1.11e-15; C3 1.0009..1.0035 at
J = 16 r^2, ranks 5..14; C4 7.64e-07; C5 J_hi^R >= J_hi at every
rank, every beat rechecked exactly in rationals; C6 holds.
(2) K1 misses, P1 holds: J_hi* = 35.2, 36.8, 63.0, 69.5, 93.5, 112.4,
131.3, 161.6, 176.1, 201.2 at ranks 5..14 against J_hi = 30, 33, 58,
60, 88, 106, 122, 152, 168, 195: J_hi*/J_hi = 1.172, 1.115, 1.086,
1.158, 1.063, 1.060, 1.076, 1.063, 1.048, 1.032, falling with the rank
but not monotone. The binding order is d = 2 (s = r - 3) at every rank
4..40 (P3); at rank 3 it is d = 1, J_hi* = 10.9, and 14.3 at rank 4,
where no cell fails and the champion attains h from J = 7 and 13.
(3) P2 holds. The exact real relaxation crosses last at J_hi^R = 34, 35,
60, 68, 88, 109, 130, 156, 170, 199, always in class s = r - 3, no
near-tie: the lattice's share J_hi^R/J_hi = 1.133, 1.061, 1.034,
1.133, 1.000, 1.028, 1.066, 1.026, 1.012, 1.021, the limit's own
J_hi*/J_hi^R = 1.034, 1.051, 1.049, 1.022, 1.063, 1.031, 1.010, 1.036,
1.036, 1.011.
(4) K2 misses, P4 holds: J_hi*/r^1.5 = 4.3771, 4.5680, 4.6024, 4.6254,
4.6466 at r = 50, 100, 200, 400, 800 against the hand 4.663, the
local exponent 400..800 1.5066, d = 2 binding among d = 1..4 at every
r read; mu/kappa at d = 2 is 1.7415 at r = 800 and 1.7428 at 801
against the hand 1.7428, ragged at small r (1.40 at rank 5, 1.52 to
1.77 over ranks 6..14).
RUN RECORD. 248.5 s, peak commit 139.5 MB under a 512 MB guard. Its exit
code is red on C1..C6, on a beat not rechecked exactly, or on K1 or K2;
P1..P4 are printed history.
"""
import math
import os
import sys
import time
from fractions import Fraction
from math import comb
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np # noqa: E402
from scipy.optimize import linprog # noqa: E402
import bandend # noqa: E402
JHI = {5: 30, 6: 33, 7: 58, 8: 60, 9: 88, 10: 106,
11: 122, 12: 152, 13: 168, 14: 195} # gap.py, bandgrowth.py
JC = {5: 47, 6: 57, 7: 97, 8: 108, 9: 157, 10: 170,
11: 227, 12: 241, 13: 305, 14: 319} # bandend.py, bandgrowth.py
YG = np.linspace(-15.0, 15.0, 30001)
PHI = np.exp(-YG ** 2 / 2)
W2 = PHI * PHI
NMAX = 41
BIG = (50, 100, 200, 201, 400, 401, 800, 801)
def htable(needed):
"""{n: h_n phi} on YG for n in needed, h_n = He_n / sqrt(n!)."""
out, prev, cur = {}, None, np.ones_like(YG)
for n in range(max(needed) + 1):
if n in needed:
out[n] = cur * PHI
nxt = (YG * cur - (n ** 0.5 * prev if n else 0)) / (n + 1) ** 0.5
prev, cur = cur, nxt
return out
def cheb(F0, Fs):
"""min over x of max |F0 + sum x_i Fs[i]| on YG, by an exchange LP:
solve on an active set, add the grid's violating local maxima, repeat.
Returns (value on the full grid at x, x)."""
k = len(Fs)
if k == 0:
return float(np.max(np.abs(F0))), np.zeros(0)
act = np.unique(np.concatenate([np.arange(0, len(YG), 15),
[int(np.argmax(np.abs(F0)))]]))
# scaled so every column has max 1: HiGHS's tolerances are absolute
s0 = float(np.max(np.abs(F0)))
A = np.array(Fs).T
cs = np.max(np.abs(A), axis=0)
A, F0 = A / cs, F0 / s0
top, x = cheb_scaled(F0, A, act, k)
return top * s0, x * s0 / cs
def cheb_scaled(F0, A, act, k):
for _ in range(40):
Aa, f0 = A[act], F0[act]
ones = np.ones((len(act), 1))
A_ub = np.vstack([np.hstack([Aa, -ones]), np.hstack([-Aa, -ones])])
b_ub = np.concatenate([-f0, f0])
c = np.zeros(k + 1)
c[-1] = 1
res = linprog(c, A_ub=A_ub, b_ub=b_ub,
bounds=[(None, None)] * k + [(0, None)], method="highs")
if res.status != 0:
raise RuntimeError("exchange LP: " + res.message)
x, t = res.x[:k], res.x[-1]
g = np.abs(F0 + A @ x)
top = float(np.max(g))
if top <= t * (1 + 1e-10):
return top, x
loc = np.where((g[1:-1] >= g[:-2]) & (g[1:-1] >= g[2:])
& (g[1:-1] > t * (1 + 1e-10)))[0] + 1
act = np.unique(np.concatenate([act, loc, [int(np.argmax(g))]]))
return top, x
def kappa_h(H, s, d):
"""kappa*_s in the h basis: the L2(dy) projection, int f^2 / int |f|."""
b = np.linalg.lstsq(np.array([H[s + m] for m in range(1, d + 1)]).T,
-H[s], rcond=None)[0]
f = H[s] + sum(bi * H[s + i] for i, bi in enumerate(b, 1))
return float(np.sum(f * f) / np.sum(np.abs(f)))
def mu_h(H, s, d, odd=False):
ms = [m for m in range(1, d + 1) if odd or m % 2 == 0]
return cheb(H[s], [H[s + m] for m in ms])[0]
def jstar_from(r, s, Mh, val):
"""4 (M_(r-1)/val_s)^(2/d) for He, from the h-basis values."""
d = r - 1 - s
logfac = 0.5 * sum(math.log(j) for j in range(s + 1, r))
return 4 * math.exp((logfac + math.log(Mh / val)) * 2 / d)
# ---------------------------------------------------------------- (R)
def columns(J, r, s):
"""exact integer columns (1 + x)^J (1 - x)^(s+k), k = 0..d, width J + r."""
base = bandend.Q(J, s)
cols, cur = [], base
for k in range(r - s):
cols.append(cur + [0] * (J + r - len(cur)))
cur = bandend.pmul(cur, [1, -1])
return cols
def real_least(J, r, s, want_x=False):
"""least height over real c' with c'(1) = 1 in class s, over its scale
S_0 = max |column 0|; returns (t, S_0, cols, x scaled)."""
cols = columns(J, r, s)
S = [max(abs(v) for v in c) for c in cols]
F = [np.array([v / S[k] for v in c]) for k, c in enumerate(cols)]
k = len(F) - 1
if k == 0:
return 1.0, S[0], cols, np.zeros(0)
A = np.array(F[1:]).T
n = A.shape[0]
ones = np.ones((n, 1))
A_ub = np.vstack([np.hstack([A, -ones]), np.hstack([-A, -ones])])
b_ub = np.concatenate([-F[0], F[0]])
c = np.zeros(k + 1)
c[-1] = 1
res = linprog(c, A_ub=A_ub, b_ub=b_ub,
bounds=[(None, None)] * k + [(0, None)], method="highs")
return float(res.x[-1]), S[0], cols, res.x[:k]
def exact_height(cols, S0, x):
"""height of col_0 + sum (x_k S_0 / S_k) col_k, in rationals."""
S = [max(abs(v) for v in c) for c in cols]
coef = [Fraction(1)] + [Fraction(float(xk)) * S0 / S[k + 1]
for k, xk in enumerate(x)]
return max(abs(sum(cf * c[i] for cf, c in zip(coef, cols)))
for i in range(len(cols[0])))
def golden(f, lo, hi, it=120):
g = (5 ** 0.5 - 1) / 2
a, b = lo, hi
c, d = b - g * (b - a), a + g * (b - a)
fc, fd = f(c), f(d)
for _ in range(it):
if fc < fd:
b, d, fd = d, c, fc
c = b - g * (b - a)
fc = f(c)
else:
a, c, fc = c, d, fd
d = a + g * (b - a)
fd = f(d)
return min(fc, fd)
def main():
t0 = time.time()
H = htable(set(range(NMAX + 1)))
# C1 and (G)
print("(G) mu_s by exchange LP, ranks 3..40", flush=True)
c1w, c6 = 0.0, True
G = {}
for r in range(3, 41):
Mh = float(np.max(np.abs(H[r - 1])))
js, ratio2 = {}, None
for s in range(r - 1):
d = r - 1 - s
mu = mu_h(H, s, d)
if d <= 4:
full = mu_h(H, s, d, odd=True)
c1w = max(c1w, abs(full / mu - 1))
kap = kappa_h(H, s, d)
c6 &= mu >= kap * (1 - 1e-9)
if d == 2:
ratio2 = mu / kap
js[s] = jstar_from(r, s, Mh, mu)
top = max(js, key=js.get)
G[r] = (js[top], r - 1 - top, js, ratio2)
row = " r = %2d J_hi* = %8.1f at d = %d; mu/kappa (d = 2) %.4f;" \
" J_hi*/r^1.5 = %.3f [%.1f s]" % (
r, js[top], r - 1 - top, ratio2, js[top] / r ** 1.5,
time.time() - t0)
if r in JHI:
row += "; J_hi = %d, J_hi*/J_hi = %.3f" % (JHI[r], js[top] / JHI[r])
print(row, flush=True)
c1 = c1w < 1e-6
print("C1 every-m LP against the even-m LP, d <= 4: worst relative "
"%.2e (%s)" % (c1w, c1))
print("C6 mu_s >= kappa*_s at every s: %s [%.1f s]"
% (c6, time.time() - t0), flush=True)
# C2
worst = 0.0
r = 5
for J in range(30, 61, 6):
for s in (3, 2):
t, S0, cols, _ = real_least(J, r, s)
Fl = [np.array([v / max(abs(u) for u in c) for v in c])
for c in cols]
if s == 3:
g = golden(lambda a: float(np.max(np.abs(Fl[0] + a * Fl[1]))),
-50, 50)
else:
g = golden(lambda b: golden(
lambda a: float(np.max(np.abs(Fl[0] + a * Fl[1]
+ b * Fl[2]))), -50, 50,
80), -50, 50, 80)
worst = max(worst, abs(g / t - 1))
c2 = worst < 1e-6
print("C2 the (R) LP against golden search, rank 5, J = 30..60, d = 1, 2:"
" worst relative %.2e (%s) [%.1f s]"
% (worst, c2, time.time() - t0), flush=True)
# C3
c3 = True
row = "C3 (R) over its limit at J = 16 r^2, s = r - 3:"
for r in range(5, 15):
J, s = 16 * r * r, r - 3
t, S0, _, _ = real_least(J, r, s)
ex = t * (S0 / bandend.champ_height(r, J))
Mh = float(np.max(np.abs(H[r - 1])))
mu = mu_h(H, s, 2)
lim = (J ** 0.5 / 2) ** 2 * math.sqrt(math.factorial(s)
/ math.factorial(r - 1)) * mu / Mh
c3 &= abs(ex / lim - 1) < 0.02
row += " %d: %.4f" % (r, ex / lim)
print(row + " (%s) [%.1f s]" % (c3, time.time() - t0), flush=True)
# (R)
print("(R) the exact real relaxation, scanning down from J_c - 1")
R, c5, exact_ok = {}, True, True
for r in range(5, 15):
H_at = {}
found = None
for J in range(JC[r] - 1, 1, -1):
Hc = bandend.champ_height(r, J)
best = None
for s in range(r - 1):
t, S0, cols, x = real_least(J, r, s)
q = t * (S0 / Hc)
if best is None or q < best[0]:
best = (q, s, cols, S0, x)
H_at[J] = best[0]
if best[0] < 1 - 1e-6:
found = J
hx = exact_height(best[2], best[3], best[4])
exact_ok &= hx < Hc
break
near = [J for J, q in H_at.items() if abs(q - 1) < 1e-5]
R[r] = found
c5 &= found >= JHI[r]
print(" r = %2d J_hi^R = %d (class s = r - %d, least/H_ch %.5f),"
" J_hi = %d, J_hi^R/J_hi = %.3f, J_hi*/J_hi^R = %.3f;"
" near-ties %s [%.1f s]"
% (r, found, r - best[1], best[0], JHI[r], found / JHI[r],
G[r][0] / found, near, time.time() - t0), flush=True)
print("C5 J_hi^R >= J_hi at every rank: %s; every beat rechecked exactly:"
" %s" % (c5, exact_ok))
# (A)
need = set()
for r in list(range(5, 41)) + list(BIG):
need |= set(range(max(0, r - 5), r))
HA = htable(need)
def mus(r, Hx):
out = {}
for d in range(1, 5):
s = r - 1 - d
F0 = Hx[s]
if d == 1:
out[d] = float(np.max(np.abs(F0)))
elif d in (2, 3):
F2 = Hx[s + 2]
out[d] = golden(lambda b: float(np.max(np.abs(F0 + b * F2))),
-5, 5)
else:
F2, F4 = Hx[s + 2], Hx[s + 4]
out[d] = golden(lambda b4: golden(
lambda b2: float(np.max(np.abs(F0 + b2 * F2 + b4 * F4))),
-5, 5, 70), -5, 5, 70)
return out
worst = 0.0
for r in range(5, 41):
m = mus(r, HA)
Mh = float(np.max(np.abs(HA[r - 1])))
for d, v in m.items():
a = jstar_from(r, r - 1 - d, Mh, v)
worst = max(worst, abs(a / G[r][2][r - 1 - d] - 1))
c4 = worst < 1e-4
print("(A) C4 against (G) at r = 5..40, d = 1..4: worst relative %.2e"
" (%s) [%.1f s]" % (worst, c4, time.time() - t0), flush=True)
top = {}
for r in BIG:
m = mus(r, HA)
Mh = float(np.max(np.abs(HA[r - 1])))
js = {d: jstar_from(r, r - 1 - d, Mh, v) for d, v in m.items()}
s = r - 3
kap = kappa_h(HA, s, 2)
bind = max(js, key=js.get)
top[r] = js[bind]
print(" r = %3d J_hi*_d, d = 1..4: %s; binding d = %d;"
" J_hi*/r^1.5 = %.4f; mu/kappa (d = 2) = %.4f [%.1f s]"
% (r, ["%.1f" % js[d] for d in range(1, 5)], bind,
top[r] / r ** 1.5, m[2] / kap, time.time() - t0), flush=True)
# verdicts
p1 = [r for r in JHI if not 0.9 <= G[r][0] / JHI[r] <= 1.2]
k1 = [r for r in JHI if not 0.8 <= G[r][0] / JHI[r] <= 1.25]
p2 = [r for r in JHI if not (1.0 <= R[r] / JHI[r] <= 1.15
and 0.9 <= G[r][0] / R[r] <= 1.1)]
p3 = [r for r in JHI if G[r][1] != 2]
expo = math.log2(top[800] / top[400])
k2 = not (4.2 <= top[800] / 800 ** 1.5 <= 5.1 and 1.45 <= expo <= 1.55)
print("P1 ranks with J_hi*/J_hi outside [0.9, 1.2]: %s" % p1)
print("P2 ranks off P2's bands: %s" % p2)
print("P3 ranks whose binding order is not d = 2: %s" % p3)
print("K1 ranks with J_hi*/J_hi outside [0.8, 1.25]: %s (K1 %s)"
% (k1, "fires" if k1 else "misses"))
print("K2 J_hi*(800)/800^1.5 = %.4f, local exponent 400..800 = %.4f"
" (K2 %s)" % (top[800] / 800 ** 1.5, expo,
"fires" if k2 else "misses"))
print("%.1f s" % (time.time() - t0))
return (not (c1 and c2 and c3 and c4 and c5 and c6 and exact_ok)
or bool(k1) or k2)
if __name__ == "__main__":
sys.exit(1 if main() else 0)