cell_exact.py
Checks: SIZE: The cell tent bound theorem; property.
Cited on: SIZE.
To run it: save cell_exact.py, with erring_cells.py, erring_label.py, erring_lower.py, erring_widths.py, erring_zones.py beside it; needs numpy.
"""
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()