erring_lower.py

"""
erring_lower.py -- the lower side of kappa_m: a width bound under the
continuum error, and the exponential widths that meet it up to a
constant.

QUESTION. Fill [0, c) with m intervals below one overflow label
[c, 2c - 1), e_m(c) the least error as in erring_label.py, kappa_m the
limsup over c of e_m(c)/c^2. Is kappa_m bounded below by a multiple of
1/log m, and above by one?

PROOF. P is 0 = c_0 < ... < c_m = 1, widths w_i, labels I_1..I_m and
overflow [1, 2), E(P) its continuum error (erring_zones.py's
cont_error: each cell's area less its largest label area).
  (1) kappa_m = inf_P E(P). Rounding cP to integers moves each count by
      O(c), so e_m(c) <= c^2 E(P) + O(m^2 c). Conversely a cut set at c
      scaled by 1/c is a P; in a cell, for each integer x the y-count of
      a label is an integer interval counted exactly, and its length is
      1-Lipschitz in x, so each label count is within the cell's x-width
      of c^2 times its area while the cell's count is exact: e >= c^2
      E(P) - O(m^2 c). So e_m(c)/c^2 converges to inf_P E(P).
  (2) E(P) >= LB(P) = (1/8) sum w_a^2 + (1/4) sum_{a<j} w_a min(w_a, w_j).
      A cell I_a x I_b errs at least the integral over y in I_b of
      e_a(y), I_a's measure off its largest label once shifted by y.
      A label boundary beta strictly inside I_a + y leaves every label
      class on one side, so e_a(y) >= w_a/2 - |beta - t|, t the shifted
      midpoint. Integrating t over a later label [c_{j-1}, c_j), both
      ends boundaries, gives at least w_a min(w_a, w_j)/4; the boundary
      1 gives w_a^2/8 on [1, 1 + w_a/2).
  (3) Each pair term w_a min(w_a, w_j) is at least min(w_a, w_j)^2, and
      the pairs' squared minima depend on the widths only as a
      multiset: with v sorted increasing and S_j = sum_{i<=j} v_i^2
      they sum to sum_i v_i^2 (m - i) = sum_{j<m} S_j. So LB >= G =
      S_m/8 + (1/4) sum_{j<m} S_j = R/4 - S_m/8 >= R/8, where
      R = S_1 + ... + S_m >= S_m.
  (4) R >= F(delta) = min(delta(2 - delta)/4,
      (1 - delta)^2/log(4m/delta^2)) at every mass delta in (0, 1):
      j0 the last index with v_1 + ... + v_j0 <= delta, A that sum
      (j0 < m as the widths sum to 1). A < delta/2 puts every later
      width above delta - A > delta/2, their sum 1 - A > 1 - delta/2, so
      R >= S_m > delta(2 - delta)/4. Else S_j0 >= A^2/j0 >=
      delta^2/(4m), and with D_j = log(S_j/S_{j-1}), v_j^2 <= S_j D_j
      (log(1 + x) >= x/(1 + x)); Cauchy-Schwarz:
      (1 - delta)^2 <= (sum_{j>j0} v_j)^2 <= R log(S_m/S_j0)
      <= R log(4m/delta^2), S_m <= 1. delta = 1/2 gives
      R >= 1/(4 log 16m), as 3/16 > 1/(4 log 16), so E(P) >= LB(P) >=
      R/8 >= 1/(32 log 16m) at every P, and kappa_m = inf_P E(P) >=
      1/(32 log 16m) at every m.
  (5) G >= F2(delta, eta) = min(delta(2 - delta)/32, eta^2/8,
      (1 - delta - eta)^2/(4 log(4m/delta^2))) for delta, eta > 0,
      delta + eta < 1, with j0 and A as in (4). A < delta/2: G >= S_m/8
      > delta(2 - delta)/32. Else S_j0 >= delta^2/(4m); if v_m >= eta,
      G >= S_m/8 >= eta^2/8; if not, the widths strictly between j0
      and m sum past 1 - delta - eta > 0, and Cauchy-Schwarz over them
      alone gives (1 - delta - eta)^2 <= log(S_{m-1}/S_j0)
      sum_{j<m} S_j. At delta = 8/log m, eta = (log m)^(-1/3) the three
      terms are (1 - o(1))/(2 log m), (log m)^(-2/3)/8 and
      (1 - o(1))/(4 log m), so kappa_m >= (1 - o(1))/(4 log m).
  (6) The upper side, by the zone construction ([0, t0) one interval,
      [t0, 1) in zones of width D, zone k a grid of width eps f(z_k) at
      an independent uniform phase; erring_zones.py's zones builds it
      with f(t) = t^p only), whose phases can be chosen to err at most
      the integral of (w(x) + w(y))/(4 w(x + y)) over x, y >= t0,
      x + y < 1, plus 2 t0 + 6(1/D + 2) eps f(1), at
      m* = 1 + sum_k (ceil(D/w_k) + 2) intervals; and kappa_m is
      nonincreasing in m: once c > m some interval [u, v) has
      v - u >= 2; split at v - 1, a split cell errs no more than the
      whole and only the sums equal to v - 1 change label, at most c
      pairs, so e_{m+1}(c) <= e_m(c) + c. f(t) = exp(K(t - 1)): the
      integrand is at most (e^{K(D - y)} + e^{K(D - x)})/4, integral at
      most 1/(2K) at t0 = D; t0 = D = K^-1.5, eps = K^-3 give
      (1/2 + O(K^-0.5))/K, the O about 8/sqrt(K), at
      m* = K^2 e^K (1 + o(1)) intervals. Given m, some K with
      m*(K) <= m has K = log m - O(log(log m)), and monotonicity carries
      the bound to m: kappa_m <= (1/2 + o(1))/log m.
  So kappa_m >= 1/(32 log 16m) at every m, and
  1/4 <= liminf kappa_m log m <= limsup kappa_m log m <= 1/2.

DESIGN. bounds computes LB, R/8 and the step-(4) floor at delta = 1/2,
floor_delta the floor F(delta), gain the bound G of step (3), f2 the
floor F2 of step (5). Partitions:
five width laws at m = 2..60 (Dirichlet; log-uniform over six decades;
many tiny widths then one wide; alternating wide and tiny; sorted
increasing), the growing-width profiles, rounding. The
discrete arm scales a P to c and reads erring_label's error against
c^2 E(P). The probe builds deterministic exponential cuts
t -> t + eps e^{K(t - 1)} from 0; the free-mass check reads F(delta)
against R, G against LB and F2 against G over the partitions and ten of
these cut sets (m up to 7673), F2 on the grid delta in {1/2, 1/4, 0.1,
8/log m when m > e^8}, eta in {1/4, 1/8, (log m)^(-1/3)}, delta + eta
< 1. delta = 8/log m enters no triple at the m run (delta + eta >= 1
up to m near 9,350 at eta = 1/8), so F2 log m at the limit's choice is
printed in closed form at m past every run, and the run prints its
count of triples at delta = 8/log m.

PREDICTIONS AND KILLS (written before the run).
  C1 min over a of E([0, a, 1]) is 3/14 at a = 2/7.
  C2 a planted break, LB with both constants set to 1, prints a ratio
     LB/E above 1 (the kill can fire).
  K1 max LB/E, max (R/8)/LB and max floor/(R/8), the floor
     1/(4 log 16m) over R, over every partition, all at most 1 + 1e-9;
     any above kills the step it checks. (Restated after the first run: the first
     run's code and this line read floor/R; the code was fixed and
     re-run.)
  K2 the discrete error at c = 200 and 400 at least c^2 E(P) - m^2 c.
  P3 E/LB bounded, 2 to 6, on the smooth profiles.
  P4 (a probe) the exponential cuts print E log m flat or falling in K.
  C3 a planted break, F(delta) with its 1/log factor dropped,
     min(delta(2 - delta)/4, (1 - delta)^2), prints a ratio F/R above 1
     on the same partitions (the kill can fire).
  K3 max F(delta)/R over every partition and delta in {1/2, 1/4, 0.1}
     at most 1 + 1e-9; above kills step (4).
  C4 a planted break, R/4 in place of G (gain with its S_m coefficient
     1/4), prints a ratio over LB above 1 (the kill can fire).
  C5 a planted break, G's S_m term dropped, the term that pays for
     the cases A < delta/2 and the wide last width (gain with S_m
     coefficient 0 against F2 without its eta term and with 1 - delta
     for 1 - delta - eta), prints a ratio above 1. The eta term dropped
     with the S_m term kept is no break this population catches.
  K4 max G/LB at most 1 + 1e-9; above kills step (3).
  K5 max F2/G at most 1 + 1e-9 over every (partition, delta, eta);
     above kills step (5).
  P5 G/(R/8) from about 1 (one wide width) to near 2 (many comparable
     widths), largest on the smooth profiles.
  P6 the best F2 on the grid times log m under 0.1 at every m run.

FINDINGS (from the printed run; C1 and C2 passed first: 0.2142857 at
a = 0.2857, and the break printed LB/E up to 3.333; C3 passed before K3,
max 1.502).
  K1 silent over 1186 partitions: the largest LB/E is 0.8325
     (rounding), every (R/8)/LB at most 1 (it nears 1 when one wide
     width carries S_m), floor/(R/8) at most 0.1618 (the p = 2
     profile), so step (4) is loose by about 6 against R/8 at these m.
     E/LB spans 1.20 to 4.00 over all of them, near 4 when one wide
     interval carries nearly all the length.
  K2 held: the least (e - c^2 E)/(m^2 c) over 24 runs is -0.0969.
  P3 failed as stated: E/LB 1.85 and 1.86 (p = 1) and 2.15
     (p = 2) on the profiles, p = 1 under the predicted 2.
  P4 failed as stated: E log m reads 0.373, 0.396, 0.394, 0.435, 0.419 at
     eps = 1/K, K = 3..7, and rises 0.401, 0.448, 0.464, 0.479 at
     eps = 1/K^2, K = 3..6, falling in neither arm, where the proof's o(1)
     is about 8/sqrt(K): no rate is read off it. The exponential cuts
     beat the growing-width profiles at nearby m: E 0.07247 at m = 406
     against p = 2's 0.0800 at m = 400, and 0.05978 at m = 1100 against
     0.0738 at m = 1000 (the profiles' limits, erring_zones.py's
     print), so kappa_1100 < 0.0598 (E bounds kappa_m above by step 1).
  K3 silent over 1196 partitions (m up to 7673), 3588 (partition,
     delta) pairs: max F/R 0.3806.
  C4 and C5 passed before K4 and K5: R/4 over LB up to 2.000, the
     dropped S_m term's F2/G up to 7.21e6.
  K4 silent over the 1196: max G/LB 1.0000.
  K5 silent over 9349 (partition, delta, eta) triples, 0 at
     delta = 8/log m: max F2/G 0.1924.
  P5 failed as stated: G/(R/8) runs 1.000 (tiny then wide, every one)
     to 1.978, its largest the exponential cuts', then alternating
     1.935 and Dirichlet 1.909; the profiles reach at most 1.788,
     rounding 1.879.
  P6 held: the best F2 log m is at most 0.0531 (exponential, m up to
     7673). At the limit's choice it reads 0.0075, 0.0519, 0.1364,
     0.1849 and 0.2111 at m = 10^8, 10^16, 10^64, 10^256 and 10^1024
     (the inline closed form equal to f2 within 2.2e-16 at m = 10^8 to
     10^256, where the log term binds), against its limit 1/4: the 1/4 is asymptotic, and no run reaches
     it.
"""

import math
import os

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np  # noqa: E402

from erring_label import error, rounding  # noqa: E402
from erring_widths import profile_at  # noqa: E402
from erring_zones import cont_error  # 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 bounds(cuts, k8=1 / 8, k4=1 / 4):
    """LB(P), R/8 and the floor 1/(32 log 16m) of a cut set."""
    w = np.diff(np.asarray(cuts, dtype=float))
    m = len(w)
    pair = sum(w[a] * np.minimum(w[a], w[a + 1:]).sum() for a in range(m))
    lb = k8 * (w * w).sum() + k4 * pair
    S = np.cumsum(np.sort(w) ** 2)
    return lb, S.sum() / 8, 1 / (32 * math.log(16 * m))


def cuts_of(w):
    w = np.asarray(w, dtype=float)
    return np.concatenate([[0.0], np.cumsum(w / w.sum())])


def laws(m, rng):
    yield "dirichlet", rng.dirichlet(np.ones(m))
    yield "log-uniform", 10.0 ** rng.uniform(-6, 0, m)
    yield "tiny then wide", np.append(np.full(m - 1, 1e-4), 1.0)
    yield "alternating", np.where(np.arange(m) % 2 == 0, 1.0, 1e-3)
    yield "sorted", np.sort(rng.exponential(1.0, m) ** 3)


def partitions(rng):
    for m in range(2, 61):
        for _ in range(4):
            for name, w in laws(m, rng):
                yield name, cuts_of(w)
    for c, p, m, u in [(4000, 1, 162, 0.8), (4000, 1, 400, 0.8),
                       (4000, 2, 400, 0.8)]:
        yield f"profile p = {p}, m = {m}", np.array(profile_at(c, p, m, u)) / c
    for m, u in [(18, 0.8), (66, 0.8), (162, 0.8)]:
        cu = rounding(4000, u * 4000 / (m - 2), u)
        yield f"rounding m = {len(cu) - 1}", np.array(cu) / 4000


def section_controls():
    print("controls")
    a = min(np.linspace(0.2, 0.36, 1601), key=lambda t: cont_error([0, t, 1]))
    e = cont_error([0, a, 1])
    check("C1: min E([0, a, 1]) is 3/14 near a = 2/7",
          abs(e - 3 / 14) < 1e-6 and abs(a - 2 / 7) < 1e-3,
          f"{e:.7f} at a = {a:.4f}, 3/14 = {3 / 14:.7f}")
    worst = max(bounds(cu, 1, 1)[0] / cont_error(cu)
                for cu in [[0, 1], [0, 0.5, 1], [0, 0.2, 0.6, 1]])
    check("C2: the planted break (constants 1, 1) prints LB/E above 1",
          worst > 1, f"max {worst:.3f}")


def section_bound():
    """K1 and P3 over every partition."""
    print("K1: the bound chain under E, over every partition")
    rng = np.random.default_rng(1567)
    rows = {}
    for name, cu in partitions(rng):
        E = cont_error(cu)
        lb, r8, floor = bounds(cu)
        key = name.split(",")[0] if name.startswith("profile") else name
        key = "rounding" if name.startswith("rounding") else key
        r = rows.setdefault(key, [0, 0.0, 0.0, 0.0, math.inf, 0.0])
        r[0] += 1
        r[1] = max(r[1], lb / E)
        r[2] = max(r[2], r8 / lb)
        r[3] = max(r[3], floor / r8)
        r[4] = min(r[4], E / lb)
        r[5] = max(r[5], E / lb)
    for key, (n, a, b, c, lo, hi) in rows.items():
        print(f"    {key:16s} n = {n:3d}: max LB/E {a:.4f}, max (R/8)/LB "
              f"{b:.4f}, max floor/(R/8) {c:.4f}; E/LB {lo:.2f}..{hi:.2f}")
    top = max(max(r[1], r[2], r[3]) for r in rows.values())
    check("K1: every ratio in the chain at most 1", top <= 1 + 1e-9,
          f"largest {top:.4f}")


def section_discrete():
    """K2: the discrete error at c against c^2 E(P) - m^2 c."""
    print("K2: step (1)'s lower count, discrete against continuum")
    rng = np.random.default_rng(7)
    worst = math.inf
    for m in (2, 3, 5, 8):
        for _ in range(3):
            P = np.sort(rng.choice(np.arange(1, 40), m - 1, replace=False))
            P = np.concatenate([[0], P, [40]]) / 40
            E = cont_error(P)
            for c in (200, 400):
                cu = tuple(int(round(x * c)) for x in P)
                d = error(cu + (2 * c - 1,), 2 * c - 1)
                slack = (d - c * c * E) / (m * m * c)
                worst = min(worst, slack)
    print(f"    least (e - c^2 E)/(m^2 c) over 24 runs: {worst:+.4f}")
    check("K2: e >= c^2 E(P) - m^2 c", worst >= -1)


def exp_cuts(K, eps):
    cuts, t = [0.0], 0.0
    while t < 1 - 1e-12:
        t = min(t + eps * math.exp(K * (t - 1)), 1.0)
        if 1 - t < 1e-9:
            t = 1.0
        cuts.append(t)
    return cuts


def section_exp():
    """P4: deterministic exponential widths, E log m against K."""
    print("P4: exponential cuts t -> t + eps e^(K (t - 1))")
    for K in (3, 4, 5, 6, 7):
        for eps in (1 / K, 1 / K ** 2):
            cu = exp_cuts(K, eps)
            m = len(cu) - 1
            if m > 3000:
                continue
            E = cont_error(cu)
            lb = bounds(cu)[0]
            print(f"    K = {K}  eps = 1/{round(1 / eps)}: m = {m:4d}, "
                  f"E {E:.5f}, E log m {E * math.log(m):.3f}, E/LB "
                  f"{E / lb:.2f}")


def floor_delta(m, delta):
    """F(delta): step (4) at a free mass delta, a lower bound on R."""
    return min(delta * (2 - delta) / 4,
               (1 - delta) ** 2 / math.log(4 * m / delta ** 2))


def sums(cu):
    """S_m, sum_{j<m} S_j and R of a cut set's sorted widths."""
    S = np.cumsum(np.sort(np.diff(np.asarray(cu, dtype=float))) ** 2)
    return S[-1], S[:-1].sum(), S.sum()


def gain(Sm, pair, k8=1 / 8):
    """G = S_m/8 + (1/4) sum_{j<m} S_j, step (3)'s bound under LB."""
    return k8 * Sm + pair / 4


def grid(m):
    L = math.log(m)
    ds = [0.5, 0.25, 0.1] + ([8 / L] if m > math.e ** 8 else [])
    es = [0.25, 0.125, L ** (-1 / 3)]
    return [(d, e) for d in ds for e in es if d + e < 1]


def f2(m, d, e, eta=True):
    """F2(delta, eta), step (5)'s floor under G (eta=False: the break)."""
    if not eta:
        return min(d * (2 - d) / 32,
                   (1 - d) ** 2 / (4 * math.log(4 * m / d ** 2)))
    return min(d * (2 - d) / 32, e * e / 8,
               (1 - d - e) ** 2 / (4 * math.log(4 * m / d ** 2)))


def population():
    rng = np.random.default_rng(1567)
    for name, cu in partitions(rng):
        key = name.split(",")[0] if name.startswith("profile") else name
        yield ("rounding" if name.startswith("rounding") else key), cu
    for K in (3, 4, 5, 6, 7):
        for eps in (1 / K, 1 / K ** 2):
            yield "exponential", exp_cuts(K, eps)


def pair_floors(cuts, k8=1 / 8, eta=True):
    """K4 and K5 over the cut sets: max G/LB, max F2/G and the rows."""
    k4, k5, n, at_limit, rows = 0.0, 0.0, 0, 0, {}
    for key, cu in cuts:
        m = len(cu) - 1
        Sm, pair, R = sums(cu)
        G = gain(Sm, pair, k8)
        k4 = max(k4, G / bounds(cu)[0])
        best = 0.0
        for d, e in grid(m):
            n += 1
            at_limit += d == 8 / math.log(m)
            v = f2(m, d, e, eta)
            best = max(best, v)
            k5 = max(k5, v / G if G > 0 else math.inf)
        r = rows.setdefault(key, [0, 0, math.inf, 0.0, 0.0])
        r[0] += 1
        r[1] = max(r[1], m)
        r[2] = min(r[2], G / (R / 8))
        r[3] = max(r[3], G / (R / 8))
        r[4] = max(r[4], best * math.log(m))
    return k4, k5, n, at_limit, rows


def section_constant():
    """C3 to C5, K3 to K5: F(delta) against R, G against LB, F2 against G."""
    print("K3: step (4) at a free mass delta, F(delta) against R")
    cuts = list(population())
    top, brk, n = 0.0, 0.0, 0
    for _, cu in cuts:
        m = len(cu) - 1
        R = sums(cu)[2]
        for d in (0.5, 0.25, 0.1):
            n += 1
            top = max(top, floor_delta(m, d) / R)
            brk = max(brk, min(d * (2 - d) / 4, (1 - d) ** 2) / R)
    check("C3: the planted break (no 1/log factor) prints F/R above 1",
          brk > 1, f"max {brk:.3f}")
    c4 = pair_floors(cuts, k8=1 / 4)[0]
    check("C4: R/4 in place of G prints a ratio over LB above 1", c4 > 1,
          f"max {c4:.3f}")
    c5 = pair_floors(cuts, k8=0, eta=False)[1]
    check("C5: G's S_m term dropped prints F2/G above 1", c5 > 1,
          f"max {c5:.3g}")
    if not all(CHECKS):
        print("    CONTROL FAILED: K3 to K5 are not read.")
        return
    print(f"    {len(cuts)} partitions (m up to "
          f"{max(len(cu) for _, cu in cuts) - 1}), {n} (partition, delta) "
          f"pairs")
    print("    the upper side's K B <= 1/2 + 8/sqrt(K) + 12/K^2 at K = 7: "
          f"{0.5 + 8 / math.sqrt(7) + 12 / 49:.2f}")
    check("K3: F(delta) <= R at every pair", top <= 1 + 1e-9,
          f"max F/R {top:.4f}")
    print("K4, K5: step (3)'s G against LB, step (5)'s F2 against G")
    k4, k5, n, at_limit, rows = pair_floors(cuts)
    print(f"    {n} (partition, delta, eta) triples, {at_limit} at "
          f"delta = 8/log m")
    for key, (k, mm, lo, hi, b) in rows.items():
        print(f"    {key:16s} n = {k:3d}, m <= {mm:4d}: G/(R/8) "
              f"{lo:.3f}..{hi:.3f}; best F2 log m {b:.4f}")
    check("K4: G <= LB at every partition", k4 <= 1 + 1e-9,
          f"max G/LB {k4:.4f}")
    check("K5: F2 <= G at every triple", k5 <= 1 + 1e-9,
          f"max F2/G {k5:.4f}")
    print("    the limit's choice delta = 8/log m, eta = (log m)^(-1/3), "
          "in closed form:")
    off = 0.0
    for ex in (8, 16, 64, 256, 1024):
        L = ex * math.log(10)
        d, e = 8 / L, L ** (-1 / 3)
        v = min(d * (2 - d) / 32, e * e / 8,
                (1 - d - e) ** 2 / (4 * (math.log(4) + L - 2 * math.log(d))))
        if ex <= 256:
            off = max(off, abs(v / f2(10.0 ** ex, d, e) - 1))
        print(f"    m = 10^{ex:<5d} F2 log m {v * L:.4f}")
    check("the closed form is f2 wherever m is a float (m <= 10^256)",
          off < 1e-12, f"max relative gap {off:.1e}")


def main():
    section_controls()
    if not all(CHECKS):
        print("\nCONTROL FAILED: nothing below is read.")
        raise SystemExit(1)
    section_bound()
    section_discrete()
    section_exp()
    section_constant()
    print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed")
    raise SystemExit(0 if all(CHECKS) else 1)


if __name__ == "__main__":
    main()