""" cell_exact.py -- LB_s >= G2(sort P) at m = 2 and m = 3, at every real partition, by an exact case split. QUESTION. P cuts [0, 1) into I_1, ..., I_m at 0 = c_0 < ... < c_m = 1, overflow [1, 2) one label. e_a(y) is w_a less the longest piece of I_a + y cut at the boundaries c_1, ..., c_m, U_ab its integral over y in I_b, and LB_s = sum_a U_aa + 2 sum over pairs of U at (the wider row, the narrower column), read at its least over the orders of equal widths (LB_low). G2(sort P) is erring_cells.py's g2 on the widths sorted increasing, ties c_(b-1) = c_j/2 in. Is LB_low >= G2(sort P) at every P with m = 2 or 3? With LB_low <= LB3 <= E, a yes gives E >= G2(sort P) there. BY HAND (fixed before the engine was written). At m = 2, with s the narrower width, L = 1 - s and M = 1 - 2s: narrow first, U_aa = T^Q_aa and LB_s - G2 = 2U_21 = s^2 (s <= 1/3), 3s - 7s^2/2 - 1/2 (1/3 <= s < 1/2); wide first, U_22 = 0 and LB_s - G2 = s/2 + s^2/2 (s <= 1/4), 5s/2 - 7s^2/2 - 1/4 (1/4 <= s < 1/2); equal widths, 3/16 = 3/16. DESIGN. Every quantity in U, T^Q and G2 is built from the cuts by affine maps, comparisons of affine forms and products of two affine forms. The engine carries each quantity as an affine form in the free cuts (Fractions) and takes each comparison's branch at an interior point of a face, recording the compared form. When no recorded form changes strict sign over the face's vertices, every point of the open face takes the same branches (an affine form >= 0 at the vertices and > 0 at one is > 0 on the open face), so each LB_pi - G2 is one quadratic there; else the face is split on that form's line. Faces: the open polygons of 0 < c_1 < c_2 < 1 (m = 3), then every polygon edge as an open segment, split alike, then every vertex, so G2's inclusive jumps and LB_s's order flips are read where they sit. At a face where widths tie, every order of the tied widths is read (LB_low >= G2 iff each is). A quadratic's least on a closed polygon is exact: the vertices, each edge's critical point, and the interior critical point when the Hessian is positive definite. U by segments of y in I_b, cut at c_(b-1), c_b, beta - c_(a-1), beta - c_a; on each the longest piece is max(p - y, q + y, r), cut at (p - q)/2, p - r, r - q. T^Q_jj by h(t) = max(0, r - dist(t, B)), cut at beta, beta +- r and the midpoints of consecutive boundaries. Both integrals by the trapezoid rule, exact on linear pieces. PREDICTIONS AND KILLS (fixed before the engine was written). X1 at m = 2 the engine's pieces equal the hand forms exactly. X2 the engine against erring_cells' shifts and sorted_g2 at a grid of widths, m = 2, 3, 4: |diff| < 1e-12; at (2, 4, 2, 5)/13, 169 LB_s = 169 G2 = 27 exactly. X3 break: the reverse charge (each pair to its narrower row) yields a certified negative face at m = 3 (a float probe read the reverse charge at 0.636 of G2 at widths (1, 1, 2)). X4 (added with the engine, not frozen) each final cell's and segment's G2 and LB, every order, equal a point evaluation at two interior points off its centroid. P LB_low >= G2 on every face at m = 2 and 3; zeros only at ties, at m = 3 (1, 1, 1) and (1, 1, 2). K a face whose exact least of LB_pi - G2 is below 0 kills it. D per sorted rank j, the least over the simplex of the row share D_a = U_aa + 2 sum_(b narrower) U_ab less the sorted share G_j = T^Q_jj + (v_j^2/2) #{i >= j: C_i <= 2C_(j-1)} (G2 is built as their sum, so they sum to it by construction); predicted: rank 1 fails, as at m = 2. FINDINGS (from the printed run; X2 max |dU| 1.1e-16 and |dG2| 8.3e-17 over 146 widths with 169 LB_s = 169 G2 = 27 at the tie, X1 0 mismatches, X4 0 off over 590, X3 least -0.062500 at (1/4, 1/4, 1/2) passed first). P held: LB_low - G2(sort P) is never below 0 at m = 2 or 3, so E >= LB3 >= LB_low >= G2(sort P) at every P with m <= 3. At m = 2, 5 segments cut at c_1 = 1/3, 1/2, 2/3, 3/4, the hand forms' seams, and one zero in the open range, c_1 = 1/2. At m = 3, 116 cells, 179 segments and 64 open vertices (474 evaluations). The zeros are not the two predicted points but the whole segment between them, widths (t, t, 1 - 2t) for 1/4 <= t <= 1/3 (w_1 = w_2, the sorted tie C_1 = C_2/2), in either order there; no other zero in the open simplex, and no cell's closure meets 0 off the domain's boundary, every order read. D (over cells, segments and vertices): the rows' shares fail at ranks 1 and 2, not 3: rank 1 least -0.027778 at (1/3, 1/3, 1/3) with ranks (2, 3, 1); rank 2 -0.083333 there with ranks (3, 1, 2); rank 3 least 0, the widest row always paying its own share. The locations are the first least met; the values are the reading. D2: at m = 3 every tail of the widest rows pays its shares (least 0 from rank 2 and from rank 3): surplus flows from wider rows to narrower. At m = 4 it flows the other way: at (1, 3, 2, 4)/10 the rank 3 row is short 3/4 (in 1/100) and the rank 2 row holds 3/2 over its share, at the tie C_2 = C_3/2 of the sorted shares, and on its open side the tail from rank 3 stays short (-0.7498 at (1, 3, 2.001, 4); +3.7508 at 1.999). At m = 2, wide first, the rank 1 row is short 2.7778 (x100) at widths (2, 1). So a charge moving surplus only from wider rows to narrower fails at m = 4, one moving it only from narrower to wider fails at m = 2 and 3: no one-way charge carries both. Run: 3.7 s, peak commit 21.9 MB, 6/6. """ import functools import itertools import os import time from fractions import Fraction as Fr os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") HALF = Fr(1, 2) CHECKS = [] def check(name, ok, detail=""): CHECKS.append((name, ok)) print(f" [{'PASS' if ok else 'FAIL'}] {name} {detail}") class Aff: """k[0] + sum k[i + 1] x_i, exact.""" __slots__ = ("k",) def __init__(self, k): self.k = tuple(k) @staticmethod def const(c, d): return Aff((Fr(c),) + (Fr(0),) * d) @staticmethod def var(i, d): k = [Fr(0)] * (d + 1) k[i + 1] = Fr(1) return Aff(k) def __add__(self, o): return Aff(a + b for a, b in zip(self.k, o.k)) def __sub__(self, o): return Aff(a - b for a, b in zip(self.k, o.k)) def __neg__(self): return Aff(-a for a in self.k) def sc(self, c): return Aff(a * c for a in self.k) def at(self, p): return self.k[0] + sum(a * x for a, x in zip(self.k[1:], p)) def mul(self, o): q = {} for i, a in enumerate(self.k): if a: for j, b in enumerate(o.k): if b: key = (min(i, j), max(i, j)) q[key] = q.get(key, 0) + a * b return Quad(q) class Quad: """sum q[(i, j)] x_i x_j over x_0 = 1, x_1, ..., x_d.""" __slots__ = ("q",) def __init__(self, q=None): self.q = {k: v for k, v in (q or {}).items() if v} def __add__(self, o): q = dict(self.q) for k, v in o.q.items(): q[k] = q.get(k, 0) + v return Quad(q) def __sub__(self, o): return self + o.sc(-1) def sc(self, c): return Quad({k: v * c for k, v in self.q.items()}) def at(self, p): x = (Fr(1),) + tuple(p) return sum(v * x[i] * x[j] for (i, j), v in self.q.items()) def along(self, P0, D): """(al, be, ga): q(P0 + t D) = al t^2 + be t + ga.""" x0 = (Fr(1),) + tuple(P0) dx = (Fr(0),) + tuple(D) al = be = ga = Fr(0) for (i, j), v in self.q.items(): al += v * dx[i] * dx[j] be += v * (x0[i] * dx[j] + dx[i] * x0[j]) ga += v * x0[i] * x0[j] return al, be, ga def __eq__(self, o): return not (self - o).q def qsum(items, zero): tot = zero for x in items: tot = tot + x return tot class Ctx: """Takes each comparison's branch at p, recording the compared form.""" def __init__(self, p): self.p = tuple(p) self.forms = {} def sign(self, f): v = f.at(self.p) if any(f.k[1:]): lead = next(a for a in f.k[1:] if a) key = tuple(a / abs(lead) for a in f.k) self.forms[key] = f return (v > 0) - (v < 0) def mx(self, f, g): return f if self.sign(f - g) >= 0 else g def sort(self, xs): xs = sorted(xs, key=functools.cmp_to_key(lambda f, g: self.sign(f - g))) out = [xs[0]] for x in xs[1:]: if self.sign(x - out[-1]) != 0: out.append(x) return out def inside(self, t, lo, hi): return self.sign(t - lo) > 0 and self.sign(hi - t) > 0 def trapezoid(pts, f, zero): return qsum(((f(u0) + f(u1)).sc(HALF).mul(u1 - u0) for u0, u1 in zip(pts, pts[1:])), zero) def U(ctx, C, a, b, zq): """U_ab: w_a less the longest piece of I_a + y, integrated over I_b.""" B = C[1:] L0, R0 = C[a], C[a + 1] w = R0 - L0 y0b, y1b = C[b], C[b + 1] cand = [t for beta in B for t in (beta - L0, beta - R0) if ctx.inside(t, y0b, y1b)] ys = ctx.sort([y0b, y1b] + cand) tot = zq for y0, y1 in zip(ys, ys[1:]): ym = (y0 + y1).sc(HALF) ins = [k for k, beta in enumerate(B) if ctx.inside(beta, L0 + ym, R0 + ym)] if not ins: continue lo, hi = ins[0], ins[-1] p, q = B[lo] - L0, R0 - B[hi] r = None for k in range(lo, hi): g = B[k + 1] - B[k] r = g if r is None else ctx.mx(r, g) brk = [(p - q).sc(HALF)] + ([p - r, r - q] if r is not None else []) pts = ctx.sort([y0, y1] + [t for t in brk if ctx.inside(t, y0, y1)]) def f(y): v = ctx.mx(p - y, q + y) return v if r is None else ctx.mx(v, r) tot = tot + w.mul(y1 - y0) - trapezoid(pts, f, zq) return tot def T(ctx, Cq, j, zq): """T^Q_jj: I_j's tent integrated over its own shifts.""" B = Cq[1:] d = len(Cq[0].k) - 1 zero = Aff.const(0, d) r = (Cq[j + 1] - Cq[j]).sc(HALF) t0, t1 = Cq[j] + r + Cq[j], Cq[j] + r + Cq[j + 1] cand = [x for beta in B for x in (beta - r, beta, beta + r)] cand += [(u + v).sc(HALF) for u, v in zip(B, B[1:])] pts = ctx.sort([t0, t1] + [x for x in cand if ctx.inside(x, t0, t1)]) def h(t): best = zero for beta in B: best = ctx.mx(best, r - ctx.mx(t - beta, beta - t)) return best return trapezoid(pts, h, zq) def evaluate(ctx, C): """U, G2 and the sorted shares; LB, the reverse charge and the row shares D for every order of the widths consistent with them.""" m = len(C) - 1 d = len(C[0].k) - 1 zq = Quad() Um = [[U(ctx, C, a, b, zq) for b in range(m)] for a in range(m)] w = [C[a + 1] - C[a] for a in range(m)] v = sorted(w, key=functools.cmp_to_key(lambda f, g: ctx.sign(f - g))) Cq = [Aff.const(0, d)] for x in v: Cq.append(Cq[-1] + x) Tq = [T(ctx, Cq, j, zq) for j in range(m)] share = [] for j in range(m): n = sum(1 for i in range(j + 1, m + 1) if ctx.sign(Cq[j].sc(2) - Cq[i]) >= 0) share.append(Tq[j] + v[j].mul(v[j]).sc(Fr(n, 2))) g2 = qsum(share, zq) less = {(a, b): ctx.sign(w[a] - w[b]) < 0 for a in range(m) for b in range(m)} diag = qsum((Um[a][a] for a in range(m)), zq) out = [] for rank in itertools.permutations(range(m)): if any(less[a, b] and rank[a] > rank[b] for a in range(m) for b in range(m)): continue D = [Um[a][a] + qsum((Um[a][b] for b in range(m) if rank[b] < rank[a]), zq).sc(2) for a in range(m)] rev = diag + qsum((Um[a][b] for a in range(m) for b in range(m) if rank[a] < rank[b]), zq).sc(2) out.append({"rank": rank, "lb": qsum(D, zq), "rev": rev, "D": D}) return {"U": Um, "g2": g2, "share": share, "orders": out} def cuts(p, m): d = m - 1 return ([Aff.const(0, d)] + [Aff.var(i, d) for i in range(d)] + [Aff.const(1, d)]) def point_cuts(widths): tot = sum(Fr(x) for x in widths) C, c = [Aff.const(0, 0)], Fr(0) for x in widths: c += Fr(x) / tot C.append(Aff.const(c, 0)) return C def at_point(widths): return evaluate(Ctx(()), point_cuts(widths)) # --- faces --------------------------------------------------------------- def centroid(P): n = len(P) return tuple(sum(x[i] for x in P) / n for i in range(len(P[0]))) def mixed(f, P): s = [f.at(x) for x in P] return any(v > 0 for v in s) and any(v < 0 for v in s) def split_poly(P, f): out = ([], []) n = len(P) for i in range(n): A, B = P[i], P[(i + 1) % n] fa, fb = f.at(A), f.at(B) if fa >= 0: out[0].append(A) if fa <= 0: out[1].append(A) if (fa > 0 and fb < 0) or (fa < 0 and fb > 0): t = fa / (fa - fb) X = tuple(a + t * (b - a) for a, b in zip(A, B)) out[0].append(X) out[1].append(X) return out def seg_min(q, A, B): """[(value, point, interior)]: the candidates for the least of q on [A, B], interior when off both ends.""" D = tuple(b - a for a, b in zip(A, B)) al, be, ga = q.along(A, D) c = [(ga, A, False), (al + be + ga, B, False)] if al > 0: t = -be / (2 * al) if 0 < t < 1: c.append((al * t * t + be * t + ga, tuple(a + t * x for a, x in zip(A, D)), True)) return c def poly_min(q, P): c = [] for i in range(len(P)): c += [(v, X, False) for v, X, _ in seg_min(q, P[i], P[(i + 1) % len(P)])] a, b, cc = q.q.get((1, 1), 0), q.q.get((1, 2), 0), q.q.get((2, 2), 0) det = 4 * a * cc - b * b if a > 0 and det > 0: dx, dy = q.q.get((0, 1), 0), q.q.get((0, 2), 0) X = (-2 * cc * dx + b * dy) / det Y = (b * dx - 2 * a * dy) / det if strictly_inside((X, Y), P): c.append((q.at((X, Y)), (X, Y), True)) return c def strictly_inside(X, P): n = len(P) s = [] for i in range(n): A, B = P[i], P[(i + 1) % n] s.append((B[0] - A[0]) * (X[1] - A[1]) - (B[1] - A[1]) * (X[0] - A[0])) return all(v > 0 for v in s) or all(v < 0 for v in s) def on_edge_of_domain(X, m): if m == 2: return X[0] in (0, 1) return X[0] == 0 or X[0] == X[1] or X[1] == 1 class Run: def __init__(self, m): self.m = m self.cells, self.segs, self.verts = [], [], {} self.evals = 0 def eval_at(self, p): self.evals += 1 ctx = Ctx(p) return ctx, evaluate(ctx, cuts(p, self.m)) def cell(self, P0): stack = [P0] while stack: P = stack.pop() ctx, r = self.eval_at(centroid(P)) f = next((f for f in ctx.forms.values() if mixed(f, P)), None) if f is None: self.cells.append((P, r)) else: stack.extend(Q for Q in split_poly(P, f) if len(Q) >= 3) def seg(self, A0, B0): stack = [(A0, B0)] while stack: A, B = stack.pop() ctx, r = self.eval_at(centroid([A, B])) f = next((f for f in ctx.forms.values() if mixed(f, [A, B])), None) if f is None: self.segs.append(((A, B), r)) for X in (A, B): self.verts.setdefault(X, None) else: fa, fb = f.at(A), f.at(B) t = fa / (fa - fb) X = tuple(a + t * (b - a) for a, b in zip(A, B)) stack += [(A, X), (X, B)] def faces(self): m = self.m if m == 2: self.seg((Fr(0),), (Fr(1),)) else: self.cell([(Fr(0), Fr(0)), (Fr(0), Fr(1)), (Fr(1), Fr(1))]) edges = set() for P, _ in self.cells: for i in range(len(P)): A, B = P[i], P[(i + 1) % len(P)] if (on_edge_of_domain(A, m) and on_edge_of_domain(B, m) and on_edge_of_domain(centroid([A, B]), m)): continue edges.add((min(A, B), max(A, B))) self.verts.setdefault(A, None) self.verts.setdefault(B, None) for A, B in sorted(edges): self.seg(A, B) for X in list(self.verts): if not on_edge_of_domain(X, m): self.verts[X] = self.eval_at(X)[1] def least_over(run, terms): """{label: (least, point, ranks)} of each quadratic terms(r, o) names, over every cell, segment and vertex and every order read there.""" best = {} def take(lab, c, rank): for v, X, _ in c: if lab not in best or v < best[lab][0]: best[lab] = (v, X, rank) for P, r in run.cells: for o in r["orders"]: for lab, q in terms(r, o).items(): take(lab, poly_min(q, P), o["rank"]) for (A, B), r in run.segs: for o in r["orders"]: for lab, q in terms(r, o).items(): take(lab, seg_min(q, A, B), o["rank"]) for X, r in run.verts.items(): if r is not None: for o in r["orders"]: for lab, q in terms(r, o).items(): take(lab, [(q.at(X), X, True)], o["rank"]) return best def on_seg(X, A, B): u = [x - a for x, a in zip(X, A)] v = [b - a for a, b in zip(A, B)] cross = u[0] * v[1] - u[1] * v[0] if len(X) == 2 else 0 dot = sum(x * y for x, y in zip(u, v)) return cross == 0 and 0 <= dot <= sum(y * y for y in v) def faces_report(run, key): """The least of key(order) - G2 over every face; the open segments on which it vanishes identically, joined; the other open-simplex points where it is 0; and the cells whose closure meets 0 off the domain's boundary.""" m = run.m least, at, _ = least_over(run, lambda r, o: {0: key(o) - r["g2"]})[0] zsegs, zeros, every = set(), [], True for (A, B), r in run.segs: D = tuple(b - a for a, b in zip(A, B)) flat = [(key(o) - r["g2"]).along(A, D) == (0, 0, 0) for o in r["orders"]] if any(flat): zsegs.add((A, B)) every &= all(flat) for o, f in zip(r["orders"], flat): if not f: q = key(o) - r["g2"] zeros += [X for v, X, i in seg_min(q, A, B) if v == 0 and i] for X, r in run.verts.items(): if r is not None and any((key(o) - r["g2"]).at(X) == 0 for o in r["orders"]): zeros.append(X) touch = 0 for P, r in run.cells: for o in r["orders"]: c = poly_min(key(o) - r["g2"], P) zeros += [X for v, X, i in c if v == 0 and i] touch += any(v == 0 and not on_edge_of_domain(X, m) for v, X, _ in c) joined = [] for A, B in sorted(zsegs): if (joined and joined[-1][1] == A and on_seg(A, joined[-1][0], B)): joined[-1] = (joined[-1][0], B) else: joined.append((A, B)) zeros = sorted(set(X for X in zeros if not on_edge_of_domain(X, m) and not any(on_seg(X, A, B) for A, B in joined))) return least, at, zeros, joined, touch, every def widths_of(X): c = (Fr(0),) + tuple(X) + (Fr(1),) return tuple(c[i + 1] - c[i] for i in range(len(c) - 1)) def fmt(X): return "(" + ", ".join(str(x) for x in X) + ")" # --- sections ------------------------------------------------------------ def hand_m2(x): """The hand forms of BY HAND: (LB_s, G2) at widths (x, 1 - x).""" s = min(x, 1 - x) L, M = 1 - s, 1 - 2 * s if s == HALF: return Fr(3, 16), Fr(3, 16) g2 = L * L / 4 - s * s / 4 if s <= Fr(1, 3) else s * s / 4 + M * M / 2 if x < HALF: gap = s * s if s <= Fr(1, 3) else 3 * s - Fr(7, 2) * s * s - HALF else: gap = (s / 2 + s * s / 2 if s <= Fr(1, 4) else Fr(5, 2) * s - Fr(7, 2) * s * s - Fr(1, 4)) return g2 + gap, g2 def section_controls(): print("CONTROLS") import numpy as np from erring_cells import shifts, sorted_g2 from erring_lower import cuts_of pts = ([(k, 23 - k) for k in range(1, 23)] + [(a, b, 17 - a - b) for a in range(1, 16) for b in range(1, 17 - a)] + [(3, 1, 4, 1), (5, 9, 2, 6), (2, 7, 1, 8), (1, 1, 2, 3)]) du = dg = 0.0 for wt in pts: r = at_point(wt) Uf = np.array([[float(x.at(())) for x in row] for row in r["U"]]) du = max(du, float(np.abs(Uf - shifts(cuts_of(wt))).max())) dg = max(dg, abs(float(r["g2"].at(())) - sorted_g2(wt))) check("X2 U and G2 against erring_cells", du < 1e-12 and dg < 1e-12, f"over {len(pts)} widths: max |dU| {du:.1e}," f" |dG2| {dg:.1e}") r = at_point((2, 4, 2, 5)) st = [o for o in r["orders"] if o["rank"] == (0, 2, 1, 3)][0] lb, g = st["lb"].at(()) * 169, r["g2"].at(()) * 169 check("X2 the tie (2, 4, 2, 5)/13", lb == 27 and g == 27, f"169 LB_s {lb}, 169 G2 {g}") def section_m2(): print("m = 2") t0 = time.time() run = Run(2) run.faces() bad = 0 for (A, B), r in run.segs: o = r["orders"] if len(o) != 1: bad += 1 continue for k in (1, 2, 3): x = A[0] + (B[0] - A[0]) * k / 4 lb, g2 = hand_m2(x) bad += (o[0]["lb"].at((x,)) != lb) + (r["g2"].at((x,)) != g2) for X, r in run.verts.items(): if r is not None: lb, g2 = hand_m2(X[0]) bad += (min(o["lb"].at(X) for o in r["orders"]) != lb or r["g2"].at(X) != g2) print(f" {len(run.segs)} segments, cut at " f"{', '.join(str(X[0]) for X in sorted(run.verts))}; " f"{run.evals} evaluations, {time.time() - t0:.1f} s") check("X1 the engine equals the hand forms", bad == 0, f"({bad} mismatches at 3 points a segment and every cut)") least, at, zeros, zsegs = faces_report(run, lambda o: o["lb"])[:4] print(f" least LB_low - G2 {least} at c_1 = {fmt(at)}; zeros in the " f"open range at {[fmt(X) for X in zeros]}, zero segments " f"{len(zsegs)}") return run, least def section_m3(): print("m = 3") t0 = time.time() run = Run(3) run.faces() print(f" {len(run.cells)} cells, {len(run.segs)} segments, " f"{sum(r is not None for r in run.verts.values())} open " f"vertices; {run.evals} evaluations, {time.time() - t0:.1f} s") off = n = 0 faces = [(P, r) for P, r in run.cells] + [((A, B), r) for (A, B), r in run.segs] for P, r in faces: Cn = centroid(P) for X in (P[0], P[1]): Y = tuple((3 * c + x) / 4 for c, x in zip(Cn, X)) q = at_point(widths_of(Y)) n += 1 off += (q["g2"].at(()) != r["g2"].at(Y) or sorted(o["lb"].at(()) for o in q["orders"]) != sorted(o["lb"].at(Y) for o in r["orders"])) check("X4 each face's quadratics against two point evaluations", off == 0, f"({off} off over {n}, cells and segments, every " "order)") least, at = faces_report(run, lambda o: o["rev"])[:2] check("X3 break: the reverse charge goes negative", least < 0, f"least {float(least):.6f} at widths {fmt(widths_of(at))}") least, at, zeros, zsegs, touch, every = faces_report( run, lambda o: o["lb"]) print(f" least LB_low - G2 {least} at widths {fmt(widths_of(at))}") print(f" zero segments (widths; every order zero there: {every}): " + "; ".join( f"{fmt(widths_of(A))} to {fmt(widths_of(B))}" for A, B in zsegs)) print(f" other zeros in the open simplex: " f"{[fmt(widths_of(X)) for X in zeros]}; cells whose closure " f"meets 0 off the domain's boundary: {touch}") return run, least def section_shares(run): print("D the row shares against the sorted shares, by rank") worst = least_over(run, lambda r, o: { o["rank"][a]: o["D"][a] - r["share"][o["rank"][a]] for a in range(run.m)}) for j in sorted(worst): v, X, rank = worst[j] print(f" rank {j + 1}: least D - G_j {float(v):+.6f} at widths " f"{fmt(widths_of(X))}, ranks {tuple(x + 1 for x in rank)}") def section_suffix(run): print("D2 the tails: the widest k rows against their shares") worst = least_over(run, lambda r, o: { k: sum((o["D"][a] - r["share"][o["rank"][a]] for a in range(run.m) if o["rank"][a] >= k), Quad()) for k in range(1, run.m)}) print(" m = 3, least over the simplex of the tail from rank " + ", ".join(f"{k + 1}: {worst[k][0]}" for k in sorted(worst))) w = (1, 3, 2, 4) r = at_point(w) o = r["orders"][0] gap = {o["rank"][a]: (o["D"][a] - r["share"][o["rank"][a]]).at(()) * 100 for a in range(4)} print(f" at {w}/10, D - G_j by rank (x100): " + ", ".join(str(gap[j]) for j in range(4)) + f"; the tail from rank 3: {gap[2] + gap[3]}") for x in (Fr(1999, 1000), Fr(2001, 1000)): r = at_point((1, 3, x, 4)) o = r["orders"][0] g = sum((o["D"][a] - r["share"][o["rank"][a]]).at(()) for a in range(4) if o["rank"][a] >= 2) print(f" at (1, 3, {x}, 4), the tail from rank 3 (x100) " f"{float(g * 100):+.4f}") r = at_point((2, 1)) for o in r["orders"]: a = o["rank"].index(0) print(f" m = 2 at widths (2, 1), the rank 1 row's D - G_j (x100) " f"{float((o['D'][a] - r['share'][0]).at(()) * 100):+.4f}") def main(): t0 = time.time() section_controls() if not all(ok for _, ok in CHECKS): print("CONTROL FAILED: nothing below is read.") raise SystemExit(1) _, l2 = section_m2() run, l3 = section_m3() section_shares(run) section_suffix(run) print("KILL") check("K least LB_low - G2 is not below 0", min(l2, l3) >= 0, f"(m = 2 {l2}, m = 3 {l3})") print(f"{sum(ok for _, ok in CHECKS)}/{len(CHECKS)} checks passed, " f"{time.time() - t0:.1f} s") raise SystemExit(0 if all(ok for _, ok in CHECKS) else 1) if __name__ == "__main__": main()