b_slow.py

"""
b_slow.py -- a slow staircase against the lemma on H: a block that
is never steep, climbing Lambda by its whole climb at almost no cost in
the gap sum B.

QUESTION. The lemma on H, the Cauchy-Schwarz route to B's constant, is
Lambda_H <= (1 + o(1)) log m whenever B < 1/log m, H the j past delta
with B_j < theta w_j, theta = 2/log m; a restatement takes out Y_c =
{j past delta : w_j^2 > c S0_j}, or Y_c's closure (re-read against the window less Y
until stable). Does a block that is never steep raise Lambda off the
closure by its climb while B stays below 1/log m?

DERIVATION (by hand). (D1) Base: the exponential profile (order_b.py),
K = log m + log x, top width about x/K; continuum, S0_j about the
integral of w over (t/2, t]: B about 1/K - 2/K^2 + x/(2K^2), so B log m
about 1 - (log x + 2 - x/2)/log m (TRANSPLANT at the top, where widths
jump by e^x: printed). (D2) With u = 1 - delta - eta and B < theta/2,
Cauchy-Schwarz on H forces Lambda_H > 4(u - 1/2)/theta over the base's
intervals alone, the inserted ones counted nowhere. (D3) Block at
delta: w_1^2 = c' M, M the base's square mass in the first member's
window, w_k^2 = c'(M + w_1^2 + ... + w_(k-1)^2); every earlier member is
in each window (block width W << delta), so w_k^2 = c' S0_k up to the
window's edge, no member is in Y_c for c' < c, and the closure never
starts on the block. Each term w_k^2/B_k = c'/(1 + c'/2); n members
climb log(w_n^2/(c' M))/log(1 + c'). (D4) In H while w_top < theta
c'/(1 + c'/2). (D5) Cost: about delta W from the base j whose windows
hold the block, the block's own B below theta W, about 2W/K from the
shifted windows: B(P') <= B(P)(1 + O(W log m)). (D6) Asymptotically, at
w_top = 1/log^4 m and c' = c_m/2, 1/log m <= c_m -> 0 (so W <= 4/log^3
m, and the base's Y_c, its steep top, has width o(1), which (D2) needs:
u falls by |Y|), Lambda off the closure reaches (2 - o(1)) log m with
B(P') < 1/log m ((D2)'s 1 - 2 delta - 2 eta plus the block's 1 -
delta), if (D1)'s gap holds.

DESIGN. delta = 0.2, eta = 0.1 (fixed: the sorted proof's 8/log m
and 2/sqrt(log m) exceed 1/2 at every reachable m, so the run is toy scale;
eta = 0.1 keeps x/K below eta at m >= 1e5),
theta = 2/log m. Base rows: x in {0.5, 1, 2}, m = 1e3 to 1e6: B log m,
(1 - B log m) log m, max w. Block rows: x = 1, m = 1e5 and 1e6 (the
block's n taken from the base, so m counts both), c' = 0.05, c = 0.1,
w_top in {0.001, 0.002, 0.005, 0.01, 0.02}. Per row: n, W, B(P') log m,
|Y_c|, |closure| and its width, Lambda off the closure for P and P'
over log m, the block's own terms against n c'/(1 + c'/2), the base's
terms' move over log m, and the block's w^2/S0 range. S0 is computed
fast: every run is increasing, so only the block's members wider than
w_j need min(w_b, w_j).

PREDICTIONS (written before the run). (P1) base B log m below 1 at
x = 1 and 2, (1 - B log m) log m near 1 to 2, not falling. (P2) the block's
w^2/S0 within 1e-3 of c', none in Y_c or its closure. (P3) Lambda off
the closure rises by n c'/(1 + c'/2) to 1%, the base's terms moving
under 1% of log m. (P4) B(P') log m below 1 where W delta log m is
below (D1)'s gap.

KILL (of "a slow staircase raises Lambda off the closure by its climb at
B below 1/log m"). At m = 1e5 or 1e6, x = 1, the largest w_top on the
grid with B(P') log m < 1: the rise in Lambda off the closure below 0.5
n c'/(1 + c'/2), or no w_top keeps B(P') log m < 1.

CONTROLS (read first). C1 the fast S0 against b_chain.terms on P' at m =
1500, max |diff| < 1e-12. C2 plant: the block's min read as w_b^2 must
break C1. C3 (positive) a block at c' = 0.08 with w_1^2 = M, w_top =
0.005, m = 1e4: Y_c takes that member alone and the closure the whole block
(each later member is c' times a mass at least 2M, so steep against M
alone). C3 as frozen (c' = 0.05, w_1^2 = 10 M) built an empty block at
m = 1500, its first member past w_top: changed after the first run.

CHANGED AFTER THE FIRST RUN. The block was inserted at the first base
interval starting past delta, so after rescaling by 1 + W it started
below delta and its lower members were not past delta (46 of 76 in H at
m = 1e5, w_top 0.001); that run's kill FIRED (m = 1e5, w_top 0.005,
then the largest with B log m < 1: rise -0.0119 against bar 0.3008).
The predictions froze eta = 0.05; the run takes 0.1, the base's top x/K being
above 0.05 at every block row (eta enters only (D2), which no row prints).
Now the insertion point moves right until the rescaled
start is past delta. P2 is checked where W <= delta/2: at W >= 0.28
the window's edge passes the base mass and the ratio reaches 0.0550.

FINDINGS (copied from the printed run; toy scale, delta = 0.2). Controls:
C1 max |diff| 1.08e-17 (m = 1576, n = 76); C2 9.52e-04; C3 n = 45, Y_c
takes 1, the closure all 45. (P1) holds but for its band: B log m 0.770
to 0.894 (x =
0.5), 0.725 to 0.860 (x = 1), 0.694 to 0.834 (x = 2) at m = 1e3 to 1e6;
(1 - B log m) log m 1.588 to 1.466 (x = 0.5, falling), 1.897 to 1.938
(x = 1), 2.111 to 2.298 (x = 2), above the continuum's 1.06, 1.5 and
1.69: below 1 everywhere, near 1 to 2 at x = 1, not at x = 2, not falling at x =
1 and 2. (P2) holds at every row with W <= delta/2 and fails as frozen
at w_top 0.01 and 0.02 (up to 0.0550). (P3) half: the block's
own terms equal n c'/(1 + c'/2) at every row with W <= delta/2 (0.3135
at m = 1e5, w_top 0.001; 0.5014 at m = 1e6, w_top 0.002), every member
in H at those rows (146 of 155 and 174 of 191 at w_top 0.02), none in
Y_c or its closure at any row; but the base's terms fall 0.091 to 0.164
log m at W 0.033 to 0.073, not under 1%: the base j in (delta, 2 delta]
gain w_j W in B_j against about w_j/K, a factor 1 + WK, WK 0.38 to 1.01
here, and some leave H once B_j passes theta w_j; no row prints the
band's share of the falls (the asymptotic W = 4/log^3 m makes WK -> 0).
(P4) holds at w_top <= 0.002: B(P') log m 0.854, 0.885 (m = 1e5) and
0.899, 0.948 (m = 1e6), against the base's 0.832 and 0.860; at w_top
0.005 it passes 1 (1.015, 1.128). Lambda off the closure: the base
0.6559 and 0.6678 log m; P' 0.9210 log m at m = 1e5 and 1.0047 log m at
m = 1e6 (w_top 0.002). KILL MISSES: m = 1e5 rise 0.2651 against bar
0.2119 (0.5 of 0.4237); m = 1e6 rise 0.3369 against 0.2507. Run: 1.8
s, peak commit 204.6 MB.

VERDICT. The slow staircase does what (D3) says, term for term: it is
never steep, so neither Y_c nor its closure takes it; every member
sits in H; its terms sum to n c'/(1 + c'/2), its climb, while B stays
below 1/log m. So the restated lemma, on H less Y_c or less its closure,
fails asymptotically wherever (D1)'s gap holds (printed to m = 1e6, not
proved): Lambda off the closure reaches (2 - o(1)) log m at B <
1/log m. The toy rise is smaller than the climb by the base's lost
terms, in part a factor 1 + WK that the asymptotic choice of W
removes. Beyond the lemma: P' has B within 1 + o(1) of P and,
asymptotically, about twice its Lambda (1.40 and 1.50 times at the toy
rows at w_top 0.002), so
no bound of B through (sum of w_j over H)^2/Lambda_H reaches the
constant: a climb in Lambda can be all but free in B.
"""

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 b_chain import terms  # noqa: E402
from order_b import profile  # noqa: E402

CHECKS = []
DELTA = 0.2


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 block_at(base, delta, cp, wtop, first=None):
    """Widths of the base (summing to 1) with the block inserted at the
    first base interval starting at or past delta (1 + W), so that it
    starts past delta once the widths are rescaled; returns (w, s, n)."""
    at = delta
    for _ in range(40):
        w, i0, n = _block(base, at, cp, wtop, first)
        tot = float(w.sum())
        if float(base[:i0].sum()) >= delta * tot:
            return w, i0, n
        at = max(at * 1.001, delta * tot)
    raise RuntimeError("block start not past delta")


def _block(base, at, cp, wtop, first):
    c = np.concatenate([[0.0], np.cumsum(base)])
    i0 = int(np.searchsorted(c[:-1], at))
    st = c[:-1]
    sq = np.concatenate([[0.0], np.cumsum(base * base)])
    blk, mass = [], 0.0
    end = c[i0]
    while True:
        # the window of a member ending near end + w: base b < i0 with
        # start >= (end + w)/2, every block member before it included
        w = 0.0
        for _ in range(3):
            lo = int(np.searchsorted(st, (end + w) / 2))
            M = sq[i0] - sq[min(lo, i0)]
            w = math.sqrt(cp * (M + mass))
        if not blk and first is not None:
            w = math.sqrt(first * M)
        if w > wtop:
            break
        blk.append(w)
        mass += w * w
        end += w
    blk = np.array(blk)
    return np.concatenate([base[:i0], blk, base[i0:]]), i0, len(blk)


def s0_fast(w, s, n, plant=False):
    """S0_j and o_j; every run increasing, the block [s, s+n) the only
    run whose members may exceed w_j inside j's window."""
    c = np.concatenate([[0.0], np.cumsum(w)])
    m = len(w)
    st, en = c[:-1], c[1:]
    lo = np.searchsorted(st, en / 2, side="left")
    lo = np.minimum(lo, np.arange(m))
    p2 = np.concatenate([[0.0], np.cumsum(w * w)])
    s0 = p2[np.arange(m)] - p2[lo]
    if n:
        bw = w[s:s + n]
        b1 = np.concatenate([[0.0], np.cumsum(bw)])
        b2 = np.concatenate([[0.0], np.cumsum(bw * bw)])
        a = np.clip(np.maximum(lo, s) - s, 0, n)
        e = np.clip(np.minimum(np.arange(m), s + n) - s, 0, n)
        f = np.searchsorted(bw, w, side="right")
        a2 = np.maximum(a, f)
        ok = e > a2
        if not plant:
            corr = np.where(ok, (b2[e] - b2[np.minimum(a2, e)])
                            - w * (b1[e] - b1[np.minimum(a2, e)]), 0.0)
            s0 = s0 - corr
    o = np.where(2 * st >= w, 0.5 * w * w, 0.0)
    return s0, o, lo


def closure(w, s0, lo, past, c):
    """Y_c and its closure: the least Y holding Y_c with every j past
    delta and w_j^2 > c (S0_j less Y's members in j's window) in Y."""
    y0 = past & (w * w > c * s0)
    Y = y0.copy()
    s0p = s0.copy()
    done = np.zeros(len(w), dtype=bool)
    hi = np.searchsorted(lo, np.arange(len(w)), side="right") - 1
    while True:
        new = np.flatnonzero(Y & ~done)
        if not len(new):
            break
        for y in new:
            j1 = hi[y]
            if j1 > y:
                sl = slice(y + 1, j1 + 1)
                s0p[sl] -= w[y] * np.minimum(w[y], w[sl])
            done[y] = True
        Y |= past & (w * w > c * s0p * (1 + 1e-12))
    return y0, Y


def read(w, s, n, c, theta):
    w = w / w.sum()
    s0, o, lo = s0_fast(w, s, n)
    B = s0 + o
    st = np.concatenate([[0.0], np.cumsum(w)])[:-1]
    past = st >= DELTA
    H = past & (B < theta * w)
    y0, Y = closure(w, s0, lo, past, c)
    term = np.divide(w * w, B, out=np.zeros_like(w), where=H & ~Y)
    return dict(w=w, B=float(B.sum()), lam=float(term.sum()), term=term,
                y0=int(y0.sum()), Y=int(Y.sum()), Yw=float(w[Y].sum()),
                Yblk=int(Y[s:s + n].sum()) if n else 0, H=H,
                ratio=(w[s:s + n] ** 2 / s0[s:s + n]) if n else None)


def section_controls():
    print("CONTROLS")
    m = 1500
    K = math.log(m)
    base = profile(m, K)
    wb, s, n = block_at(base, DELTA, 0.05, 0.01)
    w = wb / wb.sum()
    _, s0b, ob = terms(np.concatenate([[0.0], np.cumsum(w)]))
    s0, o, _ = s0_fast(w, s, n)
    d = float(np.abs(s0 - s0b).max())
    check("C1 fast S0 matches b_chain.terms on P'", d < 1e-12,
          f"m = {len(w)}, n = {n}, max |diff| {d:.2e}")
    s0p, _, _ = s0_fast(w, s, n, plant=True)
    dp = float(np.abs(s0p - s0b).max())
    check("C2 plant: the block read without min breaks C1", dp > 1e-9,
          f"max |diff| {dp:.2e}")
    base = profile(10 ** 4, math.log(10 ** 4))
    wb, s, n = block_at(base, DELTA, 0.08, 0.005, first=1.0)
    r = read(wb, s, n, 0.1, 2 / math.log(len(wb)))
    y0b = int((r["w"][s:s + n] ** 2 > 0.1 * s0_fast(
        r["w"], s, n)[0][s:s + n]).sum())
    check("C3 Y_c takes the wide first member, the closure the block",
          n >= 2 and y0b == 1 and r["Yblk"] == n,
          f"n = {n}, Y_c on the block {y0b}, closure on it {r['Yblk']}")


def main():
    t0 = time.time()
    section_controls()
    print("\nBASE  exponential profile, K = log m + log x")
    gaps = {}
    for x in (0.5, 1.0, 2.0):
        for m in (10 ** 3, 10 ** 4, 10 ** 5, 10 ** 6):
            base = profile(m, math.log(m) + math.log(x))
            s0, o, _ = s0_fast(base / base.sum(), 0, 0)
            bl = float((s0 + o).sum()) * math.log(m)
            gaps[(x, m)] = bl
            print(f"  x {x:3.1f} m {m:8d}  B log m {bl:.5f}  "
                  f"(1 - B log m) log m {(1 - bl) * math.log(m):.4f}  "
                  f"max w {base.max():.4f}")
    check("P1 base B log m below 1 at every x and m",
          all(v < 1 for v in gaps.values()))
    print("\nBLOCK  x = 1, c' = 0.05, c = 0.1; Lambda off the closure "
          "(/log m) for P and P'")
    cp, c = 0.05, 0.1
    per = cp / (1 + cp / 2)
    kill = {}
    ratios_ok = True
    for m in (10 ** 5, 10 ** 6):
        lm = math.log(m)
        theta = 2 / lm
        base0 = profile(m, lm)
        rp = read(base0, 0, 0, c, theta)
        for wtop in (0.001, 0.002, 0.005, 0.01, 0.02):
            _, _, n0 = block_at(base0, DELTA, cp, wtop)
            base = profile(m - n0, lm)
            wb, s, n = block_at(base, DELTA, cp, wtop)
            r = read(wb, s, n, c, theta)
            W = float(r["w"][s:s + n].sum())
            own = float(r["term"][s:s + n].sum())
            rest = r["lam"] - own
            rb = read(base, 0, 0, c, theta)
            rat = r["ratio"]
            if W <= DELTA / 2:
                ratios_ok &= bool(np.abs(rat - cp).max() < 1e-3)
            inH = int(r["H"][s:s + n].sum())
            st = np.concatenate([[0.0], np.cumsum(r["w"])])[s]
            assert st >= DELTA, st
            print(f"  m {len(wb):8d} w_top {wtop:.3f} n {n:4d} (in H "
                  f"{inH}) W {W:.4f}  B log m {r['B'] * lm:.5f}  "
                  f"Y_c {r['y0']} closure {r['Y']} (block {r['Yblk']}, "
                  f"width {r['Yw']:.4f})")
            print(f"      Lambda P {rp['lam'] / lm:.4f} P' "
                  f"{r['lam'] / lm:.4f} (x{r['lam'] / rp['lam']:.2f}) "
                  f"rise {(r['lam'] - rp['lam']) / lm:+.4f}"
                  f"  block {own / lm:.4f} against n c'/(1 + c'/2) "
                  f"{n * per / lm:.4f}  base terms move "
                  f"{(rest - rb['lam']) / lm:+.4f}  w^2/S0 "
                  f"{rat.min():.5f}..{rat.max():.5f}")
            if r["B"] * lm < 1:
                kill[m] = (wtop, (r["lam"] - rp["lam"]) / lm, n * per / lm)
    check("P2 block w^2/S0 within 1e-3 of c' where W <= delta/2",
          ratios_ok)
    print("\nKILL: at each m the largest w_top with B(P') log m < 1 "
          "[w_top, rise /log m, n c'/(1 + c'/2) /log m]:")
    fired = False
    for m in (10 ** 5, 10 ** 6):
        if m not in kill:
            print(f"  m {m}: no w_top keeps B log m < 1")
            fired = True
            continue
        wtop, rise, pred = kill[m]
        print(f"  m {m}: {wtop:.3f}, {rise:.4f}, {pred:.4f}")
        fired |= rise < 0.5 * pred
    print(f"  {'FIRES' if fired else 'MISSES'}")
    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())