position_charge.py

"""
position_charge.py -- does a one-way charge in POSITION order carry
LB_low(P) >= G2(sort P)?

QUESTION. cell_exact.py proves LB_low >= G2(sort P) at m <= 3 by an
exact case split, and prints that no one-way charge in RANK order pays
the sorted shares: interval a's row D_a = U_aa + 2 sum_(b narrower)
U_ab against G_(rank a). A charge in both directions along a path is
vacuous (any surplus vector with sum >= 0 is carried by the flow whose
edge (r, r + 1) holds the sum past r), so the one natural rule left
is the line itself: sigma_a = D_a - G_(rank a) with a in position
order; "earlier pays later" holds iff every prefix sum is >= 0, "later
pays earlier" iff every suffix sum is.

BY HAND (frozen before the code).
  S1 later pays earlier fails at m = 2 with the narrow interval last
     (widths (1 - s, s), s <= 1/2): I_2 + y starts at 2(1 - s) >= 1, so
     U_22 = 0, while the sorted share is s^2/4: the last suffix is
     -s^2/4 on the whole open range, least -1/16 at s = 1/2.
  S2 earlier pays later fails at (1, 3, 2, 4)/10, the sorted tie C_2 =
     C_3/2: cell_exact prints the gaps by rank (0, 3/2, -3/4, 0) x
     10^-2, positions holding ranks (1, 3, 2, 4), so the prefix of two
     is -3/4. Off the tie, at (1, 3, 2.001, 4), predicted negative (the
     tail from rank 3 is -0.7498 there and the first row pays its own
     share exactly); m = 3 unpredicted.
KILL, printed: the family dies iff each direction goes negative at an
  open face (no sorted tie inside) or at an off-tie point.
CONTROL, read first: the same sums in RANK order reproduce cell_exact's
  prints: the m = 3 tails from ranks 2 and 3 least 0 and 0, the tail
  from rank 3 at (1, 3, 2.001, 4) -0.7498, the rank 1 row at (2, 1)
  -2.7778 (x100).

DESIGN. cell_exact's engine: each face's D_a and shares are exact
quadratics in the free cuts, every order of tied widths read. The
least over an open top face (segments at m = 2, cells at m = 3) is the
least of its quadratic on the closure, and a negative there is a
negative at open points, which hold no sorted tie (each tie is one of
the compared forms' lines, so faces are split on it). Beside it,
cell_exact's least_over reads every open cell, segment and interior
vertex, every order; a width of 0 is no face, met only as a limit.

FINDINGS (copied from the printed run; controls C1..C3 matched first).
  - The kill fired, both directions. Later pays earlier: the last
    suffix least -1/16 over the open faces at m = 2 (widths (1/2,
    1/2), the open segment reaching it as -s^2/4) and at m = 3
    (suffixes of 1 and 2 both -1/16, at (0, 1/2, 1/2) and (1/2, 0,
    1/2), closure limits on the domain's boundary), and -1/12 at m = 3
    over every face (the tie (1/3, 1/3, 1/3)). Earlier pays later:
    every prefix least 0 over every cell, segment and vertex, every
    order, at m = 2 and 3, the least a limit at boundary corners
    ((1, 0), (0, 0, 1)), so >= 0, never shown attained; at the widths
    (1, 3, x, 4), normalized (x100), the prefix of two is +3.7508 at x
    = 1.999, -0.7500 at the tie x = 2 and -0.7498 at x = 2.001, where
    the sorted share of the width-3 interval counts the boundary C_3
    just inside 2 C_2; the prefix of one and the suffix of one print
    +0.0000 there (the rank-4 gap is 0 at 4 decimals).
  - So no one-way charge in position order carries LB_low >= G2(sort
    P), as none in rank order does (cell_exact.py's D2); earlier pays
    later holds at every P with m <= 3 and fails at m = 4.
  - Run: 5/5 checks, 2.6 s.
"""

import sys
import time
from fractions import Fraction as Fr

from cell_exact import (Quad, Run, at_point, fmt, least_over,
                        poly_min, seg_min, widths_of)

CHECKS = []


def check(name, ok, detail=""):
    CHECKS.append(bool(ok))
    print(f"  {'PASS' if ok else 'FAIL'}  {name}  {detail}")


def sigma(r, o):
    return [o["D"][a] - r["share"][o["rank"][a]] for a in range(len(o["D"]))]


def sums(r, o):
    """{('pre', k): sigma_1 + .. + sigma_k, ('suf', k): sigma_k + ..,
    ('tail', k): the rows of rank >= k}, k = 1 .. m - 1 (0-based)."""
    s = sigma(r, o)
    m = len(s)
    out = {}
    for k in range(1, m):
        out["pre", k] = sum(s[:k], Quad())
        out["suf", k] = sum(s[k:], Quad())
        out["tail", k] = sum((s[a] for a in range(m) if o["rank"][a] >= k),
                             Quad())
    return out


def open_least(run, top):
    """{label: (least, point, ranks)} over the open top faces only."""
    best = {}
    faces = run.cells if top == 2 else run.segs
    for F, r in faces:
        for o in r["orders"]:
            for lab, q in sums(r, o).items():
                c = poly_min(q, F) if top == 2 else seg_min(q, *F)
                for v, X, _ in c:
                    if lab not in best or v < best[lab][0]:
                        best[lab] = (v, X, o["rank"])
    return best


def at(w):
    r = at_point(w)
    return {lab: min(sums(r, o)[lab].at(()) for o in r["orders"])
            for lab in sums(r, r["orders"][0])}


def section_control(r3):
    print("C  the rank-order sums against cell_exact's prints")
    b = least_over(r3, sums)
    check("C1 m = 3 tails from ranks 2, 3 least 0, 0",
          b["tail", 1][0] == 0 and b["tail", 2][0] == 0,
          f"({b['tail', 1][0]}, {b['tail', 2][0]})")
    t = at((1, 3, Fr(2001, 1000), 4))["tail", 2] * 100
    check("C2 (1, 3, 2.001, 4) tail from rank 3 -0.7498 (x100)",
          f"{float(t):+.4f}" == "-0.7498", f"{float(t):+.4f}")
    r = at_point((2, 1))
    g = min(sigma(r, o)[o["rank"].index(0)].at(())
            for o in r["orders"]) * 100
    check("C3 (2, 1) rank 1 row -2.7778 (x100)",
          f"{float(g):+.4f}" == "-2.7778", f"{float(g):+.4f}")


def report(name, best, m):
    for kind in ("pre", "suf"):
        for k in range(1, m):
            v, X, rank = best[kind, k]
            n = k if kind == "pre" else m - k
            print(f"    {name}: least {kind}fix of {n} {float(v):+.6f}"
                  f" ({v}) at widths {fmt(widths_of(X))},"
                  f" ranks {tuple(x + 1 for x in rank)}")


def section_position(r2, r3):
    print("S  position-order sums, least over open top faces")
    b2, b3 = open_least(r2, 1), open_least(r3, 2)
    report("m = 2", b2, 2)
    report("m = 3", b3, 3)
    print("   over every cell, segment and vertex, every order:")
    for m, run in ((2, r2), (3, r3)):
        b = least_over(run, sums)
        print(f"    m = {m}: least prefixes "
              + ", ".join(f"{b[lab][0]} at {fmt(widths_of(b[lab][1]))}"
                          for lab in (("pre", k) for k in range(1, m)))
              + "; least suffixes "
              + ", ".join(f"{b[lab][0]} at {fmt(widths_of(b[lab][1]))}"
                          for lab in (("suf", k) for k in range(1, m))))
    print("   at the widths (1, 3, x, 4), normalized, x100:")
    for x in (Fr(1999, 1000), Fr(2), Fr(2001, 1000)):
        s = at((1, 3, x, 4))
        print(f"    x = {x}: prefixes "
              + ", ".join(f"{float(s['pre', k] * 100):+.4f}"
                          for k in (1, 2, 3))
              + "; suffixes "
              + ", ".join(f"{float(s['suf', k] * 100):+.4f}"
                          for k in (3, 2, 1)))
    s = at((1, 3, Fr(2001, 1000), 4))
    pre_neg = min(b2["pre", 1][0], *(b3["pre", k][0] for k in (1, 2)),
                  *(s["pre", k] for k in (1, 2, 3)))
    suf_neg = min(b2["suf", 1][0], *(b3["suf", k][0] for k in (1, 2)),
                  *(s["suf", k] for k in (1, 2, 3)))
    print("KILL")
    print(f"    least off-tie prefix {float(pre_neg):+.6f}, "
          f"least off-tie suffix {float(suf_neg):+.6f}")
    check("S1 later pays earlier goes negative off a tie", suf_neg < 0)
    check("S2 earlier pays later goes negative off a tie", pre_neg < 0)


def main():
    t0 = time.time()
    r2, r3 = Run(2), Run(3)
    r2.faces()
    r3.faces()
    section_control(r3)
    if not all(CHECKS):
        print("control red: nothing below is read")
        sys.exit(1)
    section_position(r2, r3)
    print(f"{sum(CHECKS)}/{len(CHECKS)} checks, {time.time() - t0:.1f} s")
    sys.exit(0 if all(CHECKS) else 1)


if __name__ == "__main__":
    main()