erring_order.py

"""
erring_order.py -- does putting the widths out of order buy the erring
labelling a constant: the least continuum error over every partition
against the least over nondecreasing widths, at the same m.

QUESTION. P cuts [0, 1) into m intervals, E(P) its continuum error
(erring_cells.py's cell_errs, overflow [1, 2)). Over nondecreasing
widths E(P) >= (1 - o(1))/(2 log m) is a theorem; lim sup kappa_m log m
<= 1/2. At the same m, does the least E over all P fall below the least
over sorted P by a ratio that stays away from 1 as m grows? The same
question for LB2'(P) = sum_a T_aa + 2 sum_{a<b} T_ba, the tent bound
with each pair charged at its later interval, which never sorts:
E >= LB2 >= LB2' at every P.

DERIVATION (the route this search tests; proved through (3), open
after). P in position order, c_0 = 0 not a boundary, c_m = 1 one.
  (1) LB2 >= LB2' = sum_a T_aa + 2 sum_{a<b} T_ba: max >= either side,
      so each pair is charged at its later interval b; nothing sorted.
  (2) sum_{a<b} T_ba integrates h_b over t in [mid_b, mid_b + c_(b-1)),
      and on a gap I_j (j >= 2, both ends boundaries) h_b is the larger
      end tent, integral phi = w_b^2/4 if w_j >= w_b, else
      w_j (2 w_b - w_j)/4 (erring_cells.py's tents closed form); either
      way phi >= (w_b/4) min(w_b, w_j).
  (3) I_j lies in b's sweep when j > b and c_j <= 2 c_(b-1), so
      E >= sum_a T_aa + B/2,
      S0_j = sum over b < j with 2 c_(b-1) >= c_j of w_b min(w_b, w_j),
      B = sum_j (S0_j + o_j), o_j = w_j^2/2 when 2 c_(j-1) >= w_j, else
      0: then the left half of j's own tent at c_j lies on [mid_j, c_j),
      inside j's sweep and before every gap counted, integral w_j^2/8.
      Sorted, S0_j + o_j is at least the sorted proof's S'_j - w_j^2/2.
  (4) OPEN: B >= (1 - o(1))/log m once max w and the prefix [0, delta)
      are handled as in erring_cells.py (6). In the continuum (w(t) the
      width of t's interval) B = int S0(t)/w(t) dt, S0(t) = int over
      [t/2, t] of min(w(s), w(t)). A flat Cauchy-Schwarz is expected
      to fail (an unrun estimate: a sawtooth of steep ramps keeps B a
      constant while int w/S0 grows past log m); the threshold form
      was hoped not to: with theta = 2/log m, over [delta, 1], B >= theta
      (1 - delta - |H|) + |H|^2/Lambda_H >= (1 - delta)^2/L when
      Lambda_H = int over H of w/S0 is at most L and L >= (1 - delta)
      log m, H the steep landings S0(t) < theta w(t); the target needs
      L = (1 + o(1)) log m.
      Against the potential S(t) = int_0^t w, S0 = S(t) - S(t/2) - X(t),
      X(t) = int over [t/2, t] of (w(s) - w(t))_+: the S(t/2) part is
      the sorted dyadic lemma; X is a spike's excess (one interval of
      width W above the ramp around it adds W^2 to S but W w(t) to S0).
      Owed: the damage X adds to Lambda_H is o(log m) whenever
      B < 1/log m. Heuristically, over a ramp w ~ e^(Kt), a spike at s
      saves about W/K as a landing and costs about W s as a source, a
      loss once s > 1/K. (Superseded: b_slow.py's slow staircase, each
      member c' times its window's square mass and every member in H
      (X = 0 on it), raises Lambda_H to (2 - o(1)) log m at B < 1/log m wherever
      the exponential profile keeps B below 1/log m by order 1/log^2 m,
      so neither form reaches the constant and the X lemma is moot; E >= sum T_aa + B/2 stands.)
      The prefix may be delta = 1/sqrt(log m): (1 - delta)^2 -> 1, the
      window term 8 ceil(log2(1/delta)) is O(log(log m)), and with
      eta = 2/sqrt(log m) and max w < eta (else E >= S_m/8 spends it)
      every interval past delta has 2 c_(j-1) >= eta > w_j, so carries
      o_j and w_j^2/B_j <= 2. In u = log w,
      S0(t)/w(t) = int over [t/2, t] of exp(-(u(t) - u(s))_+) ds: scale
      free, and a point below an earlier peak in its window pays the
      peak's length.

DESIGN. A compass search per m on a width vector, each coordinate moved
by +-s times the mean, s halved when a sweep finds nothing, to s below
1e-7; then 12 hops, each a random multiplicative jolt of the best and a
descent, the best kept. SORT arm: the coordinates are the increments
d_1 = w_1 > 0, d_j = w_j - w_(j-1) >= 0, so ties are reachable; starts
geometric widths at ratios 1, 1.15, 1.4 and H_n's sorted widths. ALL
arm: the coordinates are the widths; starts the SORT optimum, three
shuffles of it and three random partitions (Dirichlet(1)). Then the
SORT arm restarts from the ALL optimum sorted, so a weak SORT search
cannot inflate the gap (4 hops there). Steps start at s = 0.25, at
0.05 in a hop; a jolt multiplies by exp(N(0, 0.3)); a move is kept when
it lowers f by more than 1e-15. The frozen design named a bracketed 1D
search on the cuts and the exponential cuts as SORT starts; the compass
on widths and the geometric starts replaced them before the first run.
Beside it, the integer lattice: every
composition of N = 14 into m parts, m = 2..7, least E over all of them
against least over the nondecreasing ones.

PREDICTIONS (written before the run). (P1) rho_m = E_all/E_sort < 1 at
small m: the overflow end favours a narrow last interval. (P2) The ALL
optimum's descents (w_(j+1) < w_j) sit within the last two positions,
m - j <= 2. (P3) rho_m rises with m. (P4) the LB2' arm's ratio >= 0.9.

KILL (of "order buys no constant"). rho_m <= 0.95 at each of the three
largest m of the E arm, rising across them by no more than 0.01, AND an
ALL optimum among them with a descent at m - j >= 3. A descent at the
top end alone kills nothing.

CONTROLS (read before any verdict). C1 cell_errs against E(H_n) =
(n^2 + 3n - 1)/(8n^2), n = 2..6. C2 the search recovers a planted
minimum: f = |w - w*|^2 at a random unsorted w*, m = 8, within 1e-6.
C3 at m = 2 both arms equal a 1D grid's least (20001 points, both
orders) to 1e-6. C4 LB2' <= E at every E-arm optimum (the LB2'
arm's optima are not checked). C5 E and LB2' at least sum T_aa + B/2 at
3000 random partitions, m = 2..29, Dirichlet with parameter 0.2, 1 or 5.

FINDINGS (copied from the printed run; toy scale, m <= 16). Controls
5/5: C1 max |diff| 6.66e-16, C2 2.63e-08, C3 0.21428571 on both arms and
the grid, C4 max LB2' - E -1.35e-03, C5 (DERIVATION (3), added after the
search ran, on its own generator) least margin 4.05e-15 over 3000
random partitions. KILL missed: rho at m = 8, 10, 12
is 1.0000, 0.9995, 0.9996. E arm, rho_m from m = 3 to 12: 1.0000,
0.9999, 0.9999, 0.9999, 1.0000, 0.9995, 0.9996; the ALL optimum's
descents are drops under 10% inside runs of similar widths (m = 10:
0.695 then 0.688, 0.756 then 0.703; m = 12: 0.552 then 0.499; widths x
m); the whole gap E_sort - E_all is at most 6.6e-5 (m = 10). P1 holds in
sign only (rho is 1 at m = 3), its mechanism wrong: the ALL optima at m
= 8, 10, 12 end in their widest interval. P2 is wrong: no top-end
effect. E_all * 2 log m climbs 0.3852 to 0.5893. The lattice ratio is 1
at m = 2..6; its 0.9672 at m = 7 is the grid's: the least E over every
composition of 14 into 7 parts equals m = 6's, 0.150510. LB2' arm: rho
0.9517, 0.9688, 0.9688, 0.9700, 0.9698, 0.9714 at m = 3, 4, 6, 8, 12,
16, flat from m = 4, its ALL optimum's one descent at the bottom end:
the first interval wider than the second (1.028/m to 2.494/m; at m = 16
the second, 2.969/m, wider than the third). Charging each pair at its
later interval leaves the first interval charged by its own cell alone;
whether the derivation's delta-prefix removes the whole gap is unrun. P4
holds. P3 is not readable at these m. Run: about 7 min wall-clock (the
output file's times; the per-m times printed sum to 413 s), about 30 MB
resident read mid-run by tasklist, neither printed by the script.
"""

import os
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")

import itertools  # noqa: E402
import math  # noqa: E402
import time  # noqa: E402

import numpy as np  # noqa: E402

from erring_cells import cell_errs, tents  # 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 cuts(w):
    w = np.asarray(w, dtype=float)
    c = np.concatenate([[0.0], np.cumsum(w / w.sum())])
    c[-1] = 1.0
    return c


def E(w):
    return float(cell_errs(cuts(w)).sum())


def lb2p(w):
    T = tents(cuts(w))
    return float(np.trace(T) + 2 * np.tril(T, -1).sum())


def bsum(w):
    c = cuts(w)
    w = np.diff(c)
    tot = 0.0
    for b in range(len(w)):
        for j in range(b + 1, len(w)):
            if c[j + 1] <= 2 * c[b]:
                tot += w[b] * min(w[b], w[j])
        if 2 * c[b] >= w[b]:
            tot += w[b] ** 2 / 2
    return tot


def widths_all(x):
    x = np.maximum(x, 1e-9)
    return x / x.sum()


def widths_sort(x):
    d = np.maximum(x, 0.0)
    d[0] = max(d[0], 1e-9)
    w = np.cumsum(d)
    return w / w.sum()


def incr(w):
    w = np.sort(np.asarray(w, dtype=float))
    return np.concatenate([[w[0]], np.diff(w)])


def compass(f, x0, to_w, s0=0.25, tol=1e-7):
    x = np.array(x0, dtype=float)
    fx = f(to_w(x))
    s = s0
    m = len(x)
    while s > tol:
        moved = False
        for i in range(m):
            for sg in (1.0, -1.0):
                y = x.copy()
                y[i] = y[i] + sg * s * x.mean()
                fy = f(to_w(y))
                if fy < fx - 1e-15:
                    x, fx, moved = y, fy, True
                    break
        if not moved:
            s /= 2
    return x, fx


def search(f, starts, to_w, rng, hops=12):
    best = None
    for x0 in starts:
        x, fx = compass(f, x0, to_w)
        if best is None or fx < best[1]:
            best = (x, fx)
    for _ in range(hops):
        jolt = best[0] * np.exp(rng.normal(0, 0.3, len(best[0])))
        x, fx = compass(f, jolt, to_w, s0=0.05)
        if fx < best[1]:
            best = (x, fx)
    return to_w(best[0]), best[1]


def descents(w):
    m = len(w)
    return [m - j for j in range(1, m) if w[j] < w[j - 1] - 1e-9]


def h_sorted(n):
    return np.sort(np.concatenate([[0.5], np.ones(n - 1), [0.5]]))


def arms(f, m, rng):
    geo = [np.array([q ** j for j in range(m)]) for q in (1.0, 1.15, 1.4)]
    sstarts = [incr(g / g.sum()) for g in geo]
    if m >= 3:
        sstarts.append(incr(h_sorted(m - 1)))
    ws, fs = search(f, sstarts, widths_sort, rng)
    astarts = [ws.copy()] + [rng.permutation(ws) for _ in range(3)] + \
        [rng.dirichlet(np.ones(m)) for _ in range(3)]
    wa, fa = search(f, astarts, widths_all, rng)
    ws2, fs2 = search(f, [incr(wa)], widths_sort, rng, hops=4)
    if fs2 < fs:
        ws, fs = ws2, fs2
    return ws, fs, wa, fa


def section_controls(rng):
    print("CONTROLS")
    worst = 0.0
    for n in range(2, 7):
        w = np.concatenate([[0.5], np.ones(n - 1), [0.5]])
        worst = max(worst, abs(E(w) - (n * n + 3 * n - 1) / (8 * n * n)))
    check("C1 cell_errs against E(H_n), n = 2..6", worst < 1e-12,
          f"max |diff| {worst:.2e}")
    tgt = rng.dirichlet(np.ones(8))
    w, fx = search(lambda v: float(((v - tgt) ** 2).sum()),
                   [np.ones(8) / 8], widths_all, rng, hops=0)
    check("C2 the search recovers a planted unsorted minimum, m = 8",
          np.abs(w - tgt).max() < 1e-6, f"max |diff| {np.abs(w - tgt).max():.2e}")
    s = np.linspace(1e-6, 1 - 1e-6, 20001)
    grid = min(E([t, 1 - t]) for t in s)
    gsort = min(E([t, 1 - t]) for t in s if t <= 0.5)
    ws, fs, wa, fa = arms(E, 2, rng)
    check("C3 m = 2: both arms against a 1D grid",
          abs(fa - grid) < 1e-6 and abs(fs - gsort) < 1e-6,
          f"ALL {fa:.8f} grid {grid:.8f}, SORT {fs:.8f} grid {gsort:.8f}")
    return ws, fs, wa, fa


def section_reduction(rng):
    print("REDUCTION: E >= LB2' >= sum T_aa + B/2 at random partitions")
    worst, n = math.inf, 0
    for _ in range(3000):
        m = int(rng.integers(2, 30))
        w = rng.dirichlet(np.ones(m) * rng.choice([0.2, 1.0, 5.0]))
        lo = float(np.trace(tents(cuts(w)))) + bsum(w) / 2
        worst = min(worst, lb2p(w) - lo, E(w) - lo)
        n += 1
    check("C5 LB2' and E at least sum T_aa + B/2, m = 2..29",
          worst >= -1e-12, f"{n} partitions, least margin {worst:.2e}")


def section_lattice():
    print("LATTICE: compositions of N = 14 into m parts")
    N = 14
    for m in range(2, 8):
        best_all, best_sort, arg = math.inf, math.inf, None
        for bars in itertools.combinations(range(1, N), m - 1):
            w = np.diff(np.array((0,) + bars + (N,))) / N
            e = E(w)
            if e < best_all:
                best_all, arg = e, w
            if np.all(np.diff(w) >= 0) and e < best_sort:
                best_sort = e
        print(f"  m {m}: E_all {best_all:.6f} at {np.round(arg * N).astype(int)}"
              f", E_sort {best_sort:.6f}, ratio {best_all / best_sort:.4f}")


def section_arm(name, f, ms, rng, guard=None):
    print(f"ARM {name}")
    rows = []
    for m in ms:
        t0 = time.time()
        ws, fs, wa, fa = arms(f, m, rng)
        L = 2 * math.log(m)
        rows.append((m, fs, fa, fa / fs, descents(wa)))
        print(f"  m {m:3d}: sort {fs:.6f} all {fa:.6f} rho {fa / fs:.4f}"
              f"  gap {fs - fa:.1e}"
              f"  all*2log m {fa * L:.4f}  descents(m-j) {descents(wa)}"
              f"  [{time.time() - t0:.1f} s]")
        print("         all widths x m: "
              + " ".join(f"{v * m:.3f}" for v in wa))
        if guard is not None:
            guard.append((wa, ws))
    return rows


def main():
    rng = np.random.default_rng(1693)
    section_controls(rng)
    section_reduction(np.random.default_rng(5))
    section_lattice()
    opt = []
    rows = section_arm("E", E, [3, 4, 5, 6, 8, 10, 12], rng, opt)
    worst = max(lb2p(w) - E(w) for pair in opt for w in pair)
    check("C4 LB2' <= E at every E-arm optimum", worst <= 1e-12,
          f"max LB2' - E {worst:.2e}")
    section_arm("LB2'", lb2p, [3, 4, 6, 8, 12, 16], rng)
    last = rows[-3:]
    kill = (all(r[3] <= 0.95 for r in last)
            and last[-1][3] - last[0][3] <= 0.01
            and any(any(d >= 3 for d in r[4]) for r in last))
    print(f"KILL of 'order buys no constant': {'FIRED' if kill else 'missed'}"
          f" (rho at m = {[r[0] for r in last]}: "
          f"{[round(r[3], 4) for r in last]})")
    print(f"checks: {sum(CHECKS)}/{len(CHECKS)}")
    raise SystemExit(0 if all(CHECKS) else 1)


if __name__ == "__main__":
    main()