b_chain.py
Checks: SIZE: The cell tent bound theorem; property.
Cited on: SIZE.
To run it: save b_chain.py, with erring_cells.py, erring_label.py, erring_lower.py, erring_order.py, erring_widths.py, erring_zones.py, order_b.py beside it; needs numpy.
"""
b_chain.py -- the gap sum B telescoped: over nondecreasing widths B
carries the constant alone, and at every P a sliding-window chain
reduces the lemma on H to the window's exits and the level's net drops
(the lemma itself fails: b_slow.py).
QUESTION. E(P) >= sum_a T_aa + B/2 at every P, B = sum_j (S0_j + o_j)
(erring_order.py's DERIVATION (3)); lim kappa_m log m = 1/2 needs
B >= (1 - o(1))/log m once max w and the prefix [0, delta) are spent.
Does the sorted proof's telescoping (erring_cells.py's PROOF (5)) carry
B over nondecreasing widths, where B lacks half of each own tent, and
what does it leave owed at every P?
DERIVATION (proved; P in position order, 1-based here).
(L0) For 0 <= x < 1, x/(1 - x/2) <= -log(1 - x): the series are
sum x^n / 2^(n-1) and sum x^n / n, and 1/n >= 1/2^(n-1). So
w^2/(Q - w^2/2) <= log(Q/(Q - w^2)) whenever w^2 < Q: half of
the own tent costs the telescoping nothing.
(A) Nondecreasing widths, J the j with w_j <= c_(j-1), j0 the last j
with c_j <= delta, W = ceil(log2(1/delta)). On J, B_j = S'_j -
w_j^2/2 exactly (S0_j = S'_j - w_j^2 and o_j is there). With
S = S_j, X = S - S', a = S - w^2/2 >= S/2 and b = S' - w^2/2 >=
S'/2 >= S/8 (S' >= w^2 on J, S' >= S/4 by erring_cells.py's (4)),
1/b = 1/a + X/(ab) <= 1/a + 16 X/S^2, so
(A1) w_j^2/B_j <= log(S_j/S_(j-1)) + 16 w_j^2 X_j/S_j^2 by (L0);
(A2) over a dyadic window c_j in (U/2, U], sum w^2 X/S^2 <= 2
(erring_cells.py's (5));
(A3) Lambda^B = sum over J, j0 < j < m, of w_j^2/B_j
<= log(S_(m-1)/S_j0) + 32 W;
(A4) with erring_cells.py's (6) (A = c_j0 >= delta/2, w_m < eta,
the J widths strictly between j0 and m summing past
1 - delta - 3 eta, S_j0 >= delta^2/(4m)) and Cauchy-Schwarz,
B >= (1 - delta - 3 eta)^2 / (log(4m/delta^2) + 32 W).
At delta = 8/log m, eta = 2/sqrt(log m) that is (1 - o(1))/log m:
over nondecreasing widths B alone gives E >= (1 - o(1))/(2 log m)
unless S_m/8 does.
(C) Every P. set_j = {b < j : 2c_(b-1) >= c_j}, F(l; S) = sum over
b in S of w_b min(w_b, l), P_j = F(w_j; set_j) = S0_j, Q_j = P_j
+ w_j^2. Where o_j is there, B_j = Q_j - w_j^2/2, so by (L0)
w_j^2/B_j <= log(Q_j/P_j) at every j with S0_j > 0, with no
window error. max w < eta does not give S0_j > 0 past delta (an
empty set_j: c_(j-2) = 0.1 delta, w_(j-1) = 0.95 delta, w_j =
0.5 delta), so the chain below assumes S0_j > 0 at every j past
delta; a j with S0_j = 0 is a term w_j^2/B_j = 2 it must step
over. set_(j+1) lies in set_j + {j}, so the chain from
Q_j to P_(j+1) is an exit (the window forgets) and a level move
(w_j to w_(j+1) at fixed set), and exactly
sum_{j >= j1} log(Q_j/P_j) = log(Q_m/P*) + X + D - R,
j1 the first interval starting past delta, P* = F(max w over
set_j1; set_j1), X the exits log(F(w_j; set_j + {j}) /
F(w_j; set_(j+1))) >= 0, and the level moves log(F(w_j;
set_(j+1)) / F(w_(j+1); set_(j+1))), D their positive parts (the
drops, with the first move log(P*/P_j1)) and R their negative
parts (the rises). The lemma on H becomes a bound on the chain
total, range plus X + D - R, by (1 + o(1))log m whenever
B < 1/log m: B must pay for the window's exits and the level's net
drops. Neither part is bounded here: once w_j1 > c_(j1-1) the
count gives P* no lower bound, so the range log(Q_m/P*) no upper
one, and at two ramps (L1 = delta) it reads 0.08 to 0.10 log m
(m = 250 to 2000; the chain total less the last column), the first
drop carrying the second ramp's rise. Each j with S0_j = 0 costs 2
outside the chain, so the count of them in H must be o(log m).
DESIGN. Sorted: 2000 nondecreasing partitions, m = 3..300 (sorted
Dirichlet with parameter 0.2, 1 or 5, geometric widths at ratios 1.05
to 1.5, the sorted half-shifted grid), (delta, eta) = (1/sqrt(log m),
2/sqrt(log m)), (0.05, 0.1) and (0.1, 0.2); and the exponential profile
at order_b.py's
balanced K, m = 250, 500, 1000, 2000. Unsorted at those m, delta =
1/sqrt(log m), theta = 2/log m, H the j past delta with B_j < theta
w_j: (U1) the sorted profile; (U2) front(m, n, 0.02 m, 0.05), n = 8,
16, 32, 64 at m = 250 to 2000 (order_b.py's optimum shape); (U3)
sawtooths N = 2, 4, 8; (U4) two ramps, the profile's shape on [0, L1]
and again on [L1, 1], L1 in {delta/4, delta/2, delta}; (U5) 5 seeds,
10 random intervals of the widest tenth each moved to a random earlier
position past delta.
PREDICTIONS (written before the run). (Q1) (C)'s identity residual
below 1e-9, and (C)'s and (A1)'s per-term inequalities hold at every
term read (both are proofs). (Q2) the sorted profile: D = R = 0
exactly, X/log m below 0.3. (Q3) sawtooths N >= 2: B log m > 1. (Q4)
two ramps at L1 = delta/2: B log m < 1, Lambda_H/log m >= 0.6 and
(X + D - R)/log m < 0.3. (Q5) (A3)'s slack at the profile is at least
half its 32 W.
(Q6, at the re-freeze) two ramps at L1 = delta/2 and delta: B_past
below the sorted profile's at every m, and the chain total within
0.1 log m of range0.
KILL (of "exits and net drops are what B pays"; re-frozen after the
rehearsal). At m = 1000 and at m = 2000, some unsorted row with B_past
(B over the j past delta, the part Cauchy-Schwarz uses) below the
sorted profile's and net0/log m > 0.5, net0 = X + D - R with the first
move left out, the chain anchored at range0 = log(Q_m/P_j1). The
frozen kill (B log m < 1 and (X + D - R)/log m > 0.5 at m = 2000)
was retired: the rehearsal (m = 250, 500) read B log m 0.51 to 0.84 at
every row, sawtooths included, since a constant B (the two-ramp
shadow costs about 1/16) sits below 1/log m at these m; and P*
booked a second ramp's whole rise as its first drop. The rehearsal met
the frozen kill's condition (sawtooths N = 2 and 8, two ramps at
L1 = delta, seed 1, at m = 500), and the full run meets it at m = 2000
(two ramps at L1 = delta, 0.8997; sawtooth N = 2, 0.5335; the last
column), where "B log m < 1" holds at every row and so tells none
apart.
CONTROLS (read first). C1 the chain on the sorted profile: D = R = 0
to 1e-12 (the level moves at a tie are rounding: 3.6e-17 and 7.2e-17 read) and the
identity's residual. C2 plants, one per step, each must flag: (A1)
with 16 replaced by 0; (A2) with X_j replaced by S_j (the frozen 2 ->
0.5 never flagged, the largest window sum then read 0.41); (A3) with
32 W replaced by 0; (C) with o_j = w_j^2/4; the identity with its level
moves frozen at 0, at a sawtooth. C3 the B_j total against erring_order.bsum at 50 random
partitions. (A2)'s plant became, after the first full run, the error's
S_j^2 replaced by S_j^3: X_j -> S_j never flagged either, the window
sums sitting 1.382 below 2 at every read, so (A2)'s check catches only
gross faults. After the first full run, too: the delta = 8/log m arm
never ran (it needs m > e^8) and at delta = 1/sqrt(log m) (A4)'s base
1 - delta - 3 eta is not positive while log m <= 49 (the bound vacuous,
clipped to 0 in the code), so (A4) was read at no positive bound; the
fixed (delta, eta) arms replaced it, and B sits at least 59.5 times its
bound there, so (A4)'s check, like (A2)'s, catches only gross faults
and carries no plant. The identity's plant became the level moves
frozen at 0 (dropping X put it off by X, by construction).
FINDINGS (copied from the printed run; m = 250 to 2000). Controls:
C3 1.94e-16; C1 D 3.6e-17, R 7.2e-17, residual 0; the o_j = w^2/4
plant breaks (C) (least margin -0.350); the level moves frozen put the
identity off by 2.598 at a sawtooth; the sorted plants flag at 5903,
5865 and 5886 of 5903 reads. (A): 5903 (partition, delta) reads, 4958
meeting (A4)'s hypotheses, 2996 of those with its bound positive;
B_j = S'_j - w_j^2/2 on J to 2.1e-17; (A1) least margin 2.35e-02, (A2)
largest window sum 0.618, (A3) least slack 30.5, (A4) least B - bound
0.0866, least B/bound 59.5. At the profile Lambda^B 2.854,
3.362, 3.911, 4.438 at m = 250 to 2000, its slack 63.63 to 63.68 of
32 W = 64: the window term spends at most 0.37 (Q5 holds). (C): the
identity's residual at most 2.2e-14, the per-term margin at least
2.9e-13, no S0_j = 0 and no missing o_j past delta at any of the 52
rows (Q1 holds). KILL missed at m = 1000 and 2000. The rows with B_past
below the sorted profile's at m = 1000 and 2000 are front-loading and
the two ramps (seed 4 ties sorted to four places at m = 1000, its net0
0.037 log m; at m = 250 the sawtooth N = 2 is below too, 0.4561 against
0.4574): B_past log m at m = 2000 is 0.5098 (front), 0.5089, 0.4651,
0.4972 (L1 = delta/4, delta/2, delta) against sorted's 0.5497, and
their net0/log m is at most 0.031 (m = 1000, 2000; -0.029 and -0.036
at L1 = delta, where rises outpay the exits). (Q2) holds: D and R zero to
7e-17 at m = 500 (C1), 0.0000 printed at the rest, X/log m 0.031 to
0.058. (Q3) is wrong: B log m is 0.58 to 0.88
at the sawtooths, below 1 like every row, a constant still below
1/log m at these m. (Q4) holds at L1 = delta/2: B log m 0.51 to 0.59,
Lambda_H/log m 0.86 to 0.90, X + D - R with the first move 0.020 to
0.037 of log m (the last column). (Q6)
holds: the two ramps' B_past is below sorted's at every m (L1 = delta/2:
0.387, 0.417, 0.442, 0.465 against 0.457, 0.493, 0.524, 0.550), the
chain total within 0.037 log m of range0 (L1 = delta/2) and 0.036
log m (L1 = delta). Of the rows with B_past below sorted's, the two
ramps at L1 = delta read the largest chain total and Lambda_H: the
chain total 1.0072, 1.0010, 1.0010, 1.0001 of log m,
Lambda_H/log m 0.950, 0.946, 0.947, 0.947 (sorted's 0.72 to 0.76): by
hand, the second ramp spends the count's whole range past delta, where
sorted spends 1 - delta of it. The largest chain total read is 1.370
log m, a sawtooth N = 8 at m = 250 (Lambda_H/log m 1.333, B_past log m
0.491 above sorted's 0.457); the largest net0, 0.660 log m, a sawtooth
N = 8 at m = 500 (B_past 0.719 against 0.493). At the sawtooths N = 8,
D (without the first move) and R run 2.1 to 2.8 log m each and their
difference reads 0.05, 0.20, -0.26, -0.49 log m at m = 250 to 2000. Run: 10.3 s,
peak commit 88.8 MB.
"""
import os
import sys
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import math # noqa: E402
import time # noqa: E402
import numpy as np # noqa: E402
from erring_order import bsum, cuts, h_sorted # noqa: E402
from order_b import balanced, front, profile # noqa: E402
CHECKS = []
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
tail = f" ({detail})" if detail else ""
print(f" {'PASS' if ok else 'FAIL'} {name}{tail}")
def lagmat(c):
"""lag[b, j]: b < j and 2 c_(b-1) >= c_j (0-based: c[j+1] <= 2 c[b])."""
m = len(c) - 1
idx = np.arange(m)
return (idx[:, None] < idx[None, :]) & (c[1:][None, :] <= 2 * c[:-1][:, None])
def terms(c, own=0.5):
"""S0_j and o_j per interval, o_j = own * w_j^2 when 2 c[j] >= w_j."""
w = np.diff(c)
lag = lagmat(c)
s0 = (lag * (w[:, None] * np.minimum(w[:, None], w[None, :]))).sum(0)
o = np.where(2 * c[:-1] >= w, own * w * w, 0.0)
return w, s0, o
# ---------------------------------------------------------------- sorted
def sorted_chain(w, delta, eta, k16=16.0, k32=32.0, power=2):
"""(A1)-(A4) at one nondecreasing partition; returns a dict."""
c = cuts(w)
w = np.diff(c)
m = len(w)
w_, s0, o = terms(c)
B = s0 + o
S = np.cumsum(w * w)
lag = lagmat(c)
# S'_j: squares of b <= j with c[b] >= c[j+1]/2 (b = j included).
inwin = lag | np.eye(m, dtype=bool)
inwin &= (c[:-1][:, None] >= c[1:][None, :] / 2)
Sp = (inwin * (w * w)[:, None]).sum(0)
X = S - Sp
i0s = np.nonzero(c[1:] <= delta)[0]
if len(i0s) == 0:
return None
i0 = int(i0s[-1])
W = math.ceil(math.log2(1 / delta))
out = {"eq": 0.0, "a1": math.inf, "a2": -math.inf, "a3": math.inf,
"a4": None, "nJ": 0, "b4": 0.0}
js = [j for j in range(i0 + 1, m - 1) if w[j] <= c[j]]
lam = 0.0
win = {}
for j in js:
out["eq"] = max(out["eq"], abs(B[j] - (Sp[j] - w[j] ** 2 / 2)))
lhs = w[j] ** 2 / B[j]
err = w[j] ** 2 * X[j] / S[j] ** power
rhs = math.log(S[j] / S[j - 1]) + k16 * err
out["a1"] = min(out["a1"], rhs - lhs)
lam += lhs
k = math.floor(-math.log2(c[j + 1]))
win[k] = win.get(k, 0.0) + err
out["nJ"] = len(js)
if win:
out["a2"] = max(v - 2.0 for v in win.values())
if js:
out["a3"] = math.log(S[m - 2] / S[i0]) + k32 * W - lam
out["lam"] = lam
out["slack"] = out["a3"]
out["W"] = W
A = c[i0 + 1]
if A >= delta / 2 and w[-1] < eta:
bound = max(0.0, 1 - delta - 3 * eta) ** 2 / (
math.log(4 * m / delta ** 2) + 32 * W)
out["a4"] = float(B.sum()) - bound
out["b4"] = bound
return out
def sorted_population(rng):
pop = []
for _ in range(2000):
m = int(rng.integers(3, 301))
kind = int(rng.integers(0, 5))
if kind < 3:
w = np.sort(rng.dirichlet(np.full(m, (0.2, 1.0, 5.0)[kind])))
elif kind == 3:
w = float(rng.uniform(1.05, 1.5)) ** np.arange(m)
else:
w = np.sort(h_sorted(m - 1))
pop.append(np.maximum(w, 1e-120))
return pop
def section_sorted(rng):
print("\n(A) SORTED: the chain for B over nondecreasing widths")
pop = sorted_population(rng)
worst = {"eq": 0.0, "a1": math.inf, "a2": -math.inf, "a3": math.inf,
"a4": math.inf}
n4 = n4pos = nread = 0
plant = {"a1": 0, "a2": 0, "a3": 0}
ratio = math.inf
for w in pop:
m = len(w)
L = math.log(m)
for delta, eta in ((1 / math.sqrt(L), 2 / math.sqrt(L)),
(0.05, 0.1), (0.1, 0.2)):
r = sorted_chain(w, delta, eta)
if r is None or r["nJ"] == 0:
continue
nread += 1
for k in ("eq",):
worst[k] = max(worst[k], r[k])
worst["a1"] = min(worst["a1"], r["a1"])
worst["a2"] = max(worst["a2"], r["a2"])
worst["a3"] = min(worst["a3"], r["a3"])
if r["a4"] is not None:
n4 += 1
n4pos += r["b4"] > 0
worst["a4"] = min(worst["a4"], r["a4"])
if r["b4"] > 0:
ratio = min(ratio, (r["a4"] + r["b4"]) / r["b4"])
p1 = sorted_chain(w, delta, eta, k16=0.0)
p2 = sorted_chain(w, delta, eta, power=3)
p3 = sorted_chain(w, delta, eta, k32=0.0)
plant["a1"] += p1["a1"] < -1e-12
plant["a2"] += p2["a2"] > 1e-12
plant["a3"] += p3["a3"] < -1e-12
print(f" {nread} (partition, delta) reads; {n4} meet (A4)'s "
f"hypotheses, {n4pos} of them with its bound positive")
print(f" B_j - (S'_j - w_j^2/2) on J: max |diff| {worst['eq']:.2e}")
print(f" (A1) least margin {worst['a1']:.3e}; (A2) largest window sum "
f"minus 2 {worst['a2']:.3e}; (A3) least slack {worst['a3']:.3e}; "
f"(A4) least B - bound {worst['a4']:.3e}")
print(f" C2 plants flagged at: (A1, 16 -> 0) {plant['a1']}, "
f"(A2, S^2 -> S^3) {plant['a2']}, (A3, 32 W -> 0) {plant['a3']} "
f"of {nread} reads; (A4) least B / bound {ratio:.2f}")
check("B_j = S'_j - w_j^2/2 on J", worst["eq"] < 1e-12)
check("(A1) at every term", worst["a1"] > -1e-12)
check("(A2) at every window", worst["a2"] < 1e-12)
check("(A3) at every read", worst["a3"] > -1e-12)
check("(A4) at every read meeting its hypotheses",
n4pos > 0 and worst["a4"] > -1e-12)
check("C2 sorted plants each flag", min(plant.values()) > 0)
# ---------------------------------------------------------------- chain
def chain(w, delta, theta, own=0.5, still=False):
"""(C) at one partition; returns a dict of the chain's parts."""
c = cuts(w)
w = np.diff(c)
m = len(w)
_, s0, o = terms(c, own)
B = s0 + o
L = math.log(m)
lag = lagmat(c)
past = c[:-1] >= delta
js = np.nonzero(past)[0]
j1 = int(js[0])
H = past & (B < theta * w)
lamH = float((w[H] ** 2 / B[H]).sum())
noo = int((past & (o == 0)).sum())
zero = int((past & (s0 <= 0)).sum())
P = s0
Q = s0 + w * w
ok = past & (s0 > 0) & (o > 0)
margin = float(np.min(np.log(Q[ok] / P[ok]) - w[ok] ** 2 / B[ok]))
def F(l, setmask):
return float((w[setmask] * np.minimum(w[setmask], l)).sum())
s1 = lag[:, j1]
Pstar = F(w[s1].max(), s1) if s1.any() else 0.0
X = D = R = 0.0
if Pstar > 0 and P[j1] > 0:
D += math.log(Pstar / P[j1])
for j in range(j1, m - 1):
sn = lag[:, j + 1]
A = F(w[j], sn)
X += math.log(Q[j] / A)
lv = 0.0 if still else math.log(A / P[j + 1])
if lv > 0:
D += lv
else:
R -= lv
lhs = float(np.log(Q[j1:] / P[j1:]).sum())
rng_ = math.log(Q[m - 1] / Pstar)
resid = lhs - (rng_ + X + D - R)
d0 = math.log(Pstar / P[j1]) if Pstar > 0 and P[j1] > 0 else 0.0
range0 = math.log(Q[m - 1] / P[j1])
net0 = X + D - d0 - R
bpast = float(B[past].sum()) * L
return {"BL": float(B.sum()) * L, "H": int(H.sum()), "lamH": lamH / L,
"range": rng_ / L, "X": X / L, "D": D / L, "R": R / L,
"net": (X + D - R) / L, "resid": resid, "margin": margin,
"zero": zero, "noo": noo, "Xraw": X, "lhs": lhs,
"range0": range0 / L, "net0": net0 / L, "bpast": bpast,
"total": lhs / L, "D0": (D - d0) / L}
def two_ramps(m, K, L1):
m1 = max(2, int(round(m * L1)))
a = profile(m1, K) * L1
b = profile(m - m1, K) * (1 - L1)
return np.concatenate([a, b])
def section_chain(ms):
print("\n(C) EVERY P: the chain, delta = 1/sqrt(log m), theta = 2/log m")
print(" columns: B log m, B_past log m, |H|, then /log m: Lambda_H, "
"the chain total, range0, X, D (no first move), R, net0; "
"residual, least margin, S0 = 0, no o_j, X + D - R /log m with "
"the first move")
rows = []
rng = np.random.default_rng(20261020)
for m in ms:
t0 = time.time()
L = math.log(m)
delta, theta = 1 / math.sqrt(L), 2 / L
_, K, ws = balanced(m)
fam = [("U1 sorted", ws)]
a = sorted_chain(ws, delta, 2 / math.sqrt(L))
print(f" (A3) at the profile: Lambda^B {a['lam']:.4f}, slack "
f"{a['slack']:.4f} of 32 W = {32 * a['W']}")
n = {250: 8, 500: 16, 1000: 32, 2000: 64}[m]
fam.append((f"U2 front n={n}",
ws[front(m, n, int(round(0.02 * m)), 0.05)]))
for N in (2, 4, 8):
fam.append((f"U3 saw N={N}",
ws[[k for r in range(N) for k in range(r, m, N)]]))
for f in (0.25, 0.5, 1.0):
fam.append((f"U4 ramps L1={f:g}d", two_ramps(m, K, f * delta)))
cw = cuts(ws)
first = int(np.nonzero(cw[:-1] >= delta)[0][0])
for s in range(5):
o = list(range(m))
for _ in range(10):
k = int(rng.integers(int(0.9 * m), m))
x = o.index(k)
if x <= first + 1:
continue
o.pop(x)
o.insert(int(rng.integers(first, x)), k)
fam.append((f"U5 seed {s}", ws[o]))
print(f"m {m}: K {K:.3f}, log m {L:.3f}, delta {delta:.4f}")
for name, w in fam:
r = chain(w, delta, theta)
rows.append((m, name, r))
print(f" {name:17s} {r['BL']:6.4f} {r['bpast']:6.4f}"
f" {r['H']:5d} {r['lamH']:6.4f} {r['total']:6.4f}"
f" {r['range0']:6.4f} {r['X']:6.4f} {r['D0']:6.4f}"
f" {r['R']:6.4f} {r['net0']:7.4f} {r['resid']:8.1e}"
f" {r['margin']:8.1e} {r['zero']:2d} {r['noo']:2d}"
f" {r['net']:7.4f}")
print(f" [{time.time() - t0:.1f} s]")
return rows
def section_controls(rng):
print("\nCONTROLS")
worst = 0.0
for _ in range(50):
m = int(rng.integers(2, 40))
w = rng.dirichlet(np.ones(m))
c = cuts(w)
_, s0, o = terms(c)
worst = max(worst, abs(float((s0 + o).sum()) - bsum(w)))
check("C3 B_j total against bsum at 50 random partitions",
worst < 1e-12, f"max |diff| {worst:.2e}")
m = 500
L = math.log(m)
_, K, ws = balanced(m)
r = chain(ws, 1 / math.sqrt(L), 2 / L)
check("C1 sorted profile: D = R = 0 to 1e-12, identity exact",
r["D"] < 1e-12 and r["R"] < 1e-12 and abs(r["resid"]) < 1e-9,
f"D {r['D']:.2e}, R {r['R']:.2e}, residual {r['resid']:.1e}")
p = chain(ws, 1 / math.sqrt(L), 2 / L, own=0.25)
check("C2 plant o_j = w^2/4 breaks (C)'s per-term inequality",
p["margin"] < -1e-12, f"least margin {p['margin']:.2e}")
saw = ws[[k for q in range(4) for k in range(q, m, 4)]]
s = chain(saw, 1 / math.sqrt(L), 2 / L, still=True)
check("C2 plant level moves frozen breaks the identity at a sawtooth",
abs(s["resid"]) > 1e-6, f"residual {s['resid']:.3e}")
def main():
t0 = time.time()
rng = np.random.default_rng(1695)
section_controls(rng)
section_sorted(rng)
ms = [250, 500] if "--quick" in sys.argv else [250, 500, 1000, 2000]
rows = section_chain(ms)
check("(Q1) identity residual below 1e-9 at every row",
max(abs(r["resid"]) for _, _, r in rows) < 1e-9)
check("(Q1) (C)'s per-term inequality at every row",
min(r["margin"] for _, _, r in rows) > -1e-12)
srt = {mm: r["bpast"] for mm, n, r in rows if n == "U1 sorted"}
for mm in (1000, 2000):
if mm not in srt:
continue
hit = [(n, round(r["bpast"], 4), round(r["net0"], 4))
for m2, n, r in rows if m2 == mm and n != "U1 sorted"
and r["bpast"] < srt[mm] and r["net0"] > 0.5]
print(f"\nKILL at m = {mm}: rows with B_past below sorted's "
f"({srt[mm]:.4f}) and net0/log m > 0.5: "
f"{hit if hit else 'none'}")
print(f"\n{sum(CHECKS)}/{len(CHECKS)} checks pass "
f"[{time.time() - t0:.1f} s]")
return 0 if all(CHECKS) else 1
if __name__ == "__main__":
sys.exit(main())