integral.py
Checks: READING: The carry bound rule.
Cited on: READING.
To run it: save integral.py, with real_support.py beside it; needs mpmath, numpy, scipy.
"""
integral.py -- the carry bound at the d the real support and the pair
graph leave open (11, 12, 16, 20), by the carries' integrality.
QUESTION. Over every bounded relaxed path from a start (x integral,
delta_t = x_t - (x_(t+1) + ... + x_(t+d)) in {-1, 0, 1}, first window
sum 0), does a carry x_t or a window sum x_(t+1) + ... + x_(t+d) reach
C_d + 1, C_d = max(1, floor((d - 1)/4))? The real support bounds both
by L1 norms H_x, H_s, which pass C_d + 1 at d = 11, 12, 16, 20 (3.04,
3.24, 4.11, 5.03); there the integrality of x is put back near t.
BY HAND. x_t = g^(t) . delta (real_support.py), g^(t)_k = h_(k - t) +
beta^(-t) gam_k, gam = -sum_(s < d) h_(. - s)/(1 + ... + beta^(-(d-1))).
For each t a feasibility MILP asks for z_t >= C_d + 1 (or z_t - delta_t,
the window sum). Its columns: delta_m integral on [lo_m, F), real in
[-1, 1] on [F, F2), F = 10d, F2 past where h's forward tail falls below
1e-3; z_u integral in [-X, X], X = floor(H_x) = C_d + 1 at every open d.
Its rows: |row_u . delta - z_u| <= tau_u + 1e-8, tau_u the mass of row u
outside the columns (and its coefficients below 1e-9, moved there for
conditioning), and z_u - (z_(u+1) + ... + z_(u+d)) = delta_u wherever
the window lies inside. h is cut where its terms reach 1e-16 forward
and at 200d back, a lost mass below 1e-13, inside the 1e-8.
STATIONARY (every t >= B = 5 at once): coordinates relative to t, m in
[-B, F2), rows u in [-B + 1, F - d), row u being h_(m - u) plus
beta^(-u) w, w = beta^(-t) gam . delta the beta-mode as one real column,
cut loose from gam . delta in |w| <= beta^(-B) d L/(1 + ... +
beta^(-(d-1))), since ||gam||_1 <= d L/(1 + ... + beta^(-(d-1))) and
|delta| <= 1. At t >= B its rows sit at absolute times t + u >= 1 and
its deltas at t + m >= 0, so every path's window near t is a point of
it. B is a design variable: a B-point restricted to m >= -B', rows u >=
-B' + 1 is a B'-point for B' < B (the dropped columns' mass goes to the
larger tau, and |w| shrinks with B), so a smaller B is weaker; B = 5
keeps the verdicts of B = 10 and halves the absolute arms.
ABSOLUTE (one t < B at a time): rows g^(u), u in [0, F - d), k in [0,
F2), the first window sum 0. The mass of g^(u) past F2 is h's forward
tail T(F2 - u) plus beta^(-u) gam's, at most beta^(-u) d T(F2 - d + 1),
T(n) = sum_(m >= n) |h_m|; both are at most T(F2 - u - d), so tau_u =
T(F2 - u - d)(1 + d beta^(-u)) at u >= 0.
Only the positive side is asked: x -> -x flips delta, keeps the first
window sum 0, and every box is symmetric. Infeasible (status 2) at every
arm of d proves the bound at d.
PREDICTIONS, CONTROL AND KILL (fixed before the runs).
R1 at B = 10, d = 11, 12, the stationary carry and window sum
infeasible (TRANSPLANT from the stationary optimum without w at
B = 30, 2.0001); R2 every absolute t < 10 infeasible there; R3
likewise at d = 16, 20 (TRANSPLANT from 11, 12: the real
support's slack there 0.11 and 0.03).
KILL of the route at d: a feasible print at that d.
CONTROL, read first: at d = 13, target 3, the stationary carry and
window sum, the absolute t = 3 carry and t = 2 window sum all
feasible, the triangle wave being a point of each (carry 3 at t =
3 + 12j, window sum 3 at t = 2 + 12j).
Fixed after R1..R3's run, for the bare run:
T1 at B = 5, presolve off, the stationary arms infeasible at d = 11,
12, 16, 20 (TRANSPLANT from B = 10); T2 the absolute t = 0..4
infeasible there, one solve each; T3 the CONTROL at B = 5,
presolve off; T4 at d = 11 the stationary carry at target 3 as an
LP (integrality dropped) feasible (property, from B = 10's).
KILL of B = 5: a stationary arm feasible or unread at any of the four
d (then B = 10 kept); KILL of the fit: the bare run past 900 s.
FINDINGS (copied from the runs' prints).
- CONTROL holds: all four arms feasible at d = 13, presolve on at
B = 10 and presolve off at B = 5 (T3; the bare run's 26.3, 63.3,
30.1, 11.1 s).
- R1, R2, R3 hold at B = 10, every arm infeasible on two solves
(HiGHS's presolve on and off). T1, T2 hold at B = 5, one solve each.
The bound C_d holds at d = 11, 12, 16, 20 (C_d + 1 = 3, 3, 4, 5)
over every bounded relaxed path, legal or not, to HiGHS's floating
verdicts (below). The kills are silent; the bare run takes 621.4 s.
- Integrality refutes what the reals admit: at d = 11 the stationary
carry's LP at target 3 is feasible (T4, at B = 5 and B = 10), the
MILP is not; at t = 0 the LP alone is infeasible. At B = 10 the same
rows at target 2 are feasible (stationary carry and sum, t = 5 carry
and sum; at t = 0 not).
- A first form, one MILP for every t < 69 with a binary per (t, kind)
and their sum >= 1, found no point of d = 13 in 600 s (its big-M
relaxation is weak); one t at a time finds it. A search for
eventually periodic paths (one MILP per period) found the wave's C_d
at d = 9, 13, and found no lasso of 3 at d = 12 to n = 48, its 30 s
time-outs unread; the stationary optimum without w (valid at t >=
69) reads 2.000116 at d = 11. Neither is needed by the proof, and
both are cut.
- Float: at B = 10 HiGHS shifts bounds (values to 1.6e14) in 3 arms
at d = 12, each in the presolve-on solve (presolve off: no shift),
and in 3 arms at d = 16, the solve and size unread; at B = 5, in
the presolve-off solve of d = 16's t = 4 window sum (values to
-8.3e10; presolve on: no shift, infeasible), and none at 11, 12 or
20. Every arm solved twice agrees. The verdicts are HiGHS's
floating ones (its bound shifts at feasibility 1e-6 times 1 plus
a random fraction, as the shift prints read; the rows' tau +
1e-8), no exact certificate.
RUN RECORD: the bare run (B = 5, presolve off) 621.4 s, peak commit
349.7 MB, 6/6 checks, each arm within 0.3 s of the runs before it,
d = 16's shift as recorded. Before it, at B = 5, one process each
under a 512 MB ceiling: d = 13 130.4 s, 214.9 MB; d = 11, 12 246.4 s,
185.8 MB (d = 12's stationary arms 114.5 and 78.2 s); d = 16 110.5 s,
341.9 MB; d = 20 134.0 s, 297.6 MB. At B = 10, each arm solved twice
but d = 13's (presolve on), the four runs under a 512 MB ceiling: d =
13 130.1 s, 211.9 MB; TLIM=300 d = 11, 12 343.5 s, 179.0 MB; TLIM=300 d = 16 705.3 s, 300.4 MB;
TLIM=600 d = 20 520.4 s, 447.4 MB; the d = 11 LP and target-2 checks
684.3 s, 192.2 MB; the d = 12 presolve-off and B = 15, 20 checks
(infeasible, B = 20 without a shift) 272.1 s, 193.8 MB.
"""
import math
import os
import sys
import time
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np # noqa: E402
from scipy.optimize import Bounds, LinearConstraint, milp # noqa: E402
from scipy.sparse import coo_matrix # noqa: E402
from real_support import by_roots, norms # noqa: E402
TOL = 1e-8
B = 5
CHECKS = []
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
print(f" [{'PASS' if ok else 'FAIL'}] {name}" +
(f" -- {detail}" if detail else ""), flush=True)
def C_of(d):
return max(1, (d - 1) // 4)
def model(d, t, B=B, F=None, F2=None):
"""The rows at time t (None: stationary, every t >= B), sparse: a dict
with the coordinate list A of (row, col, value), lo, hi, the delta
count M (the first F - lo_m integral), the z columns zi, the box X,
the row u0 of z_t, and ext, the beta-mode w's column (stationary)."""
hf, hb = by_roots(d)
L, Ls, nx, ns = norms(d, hf, hb)
X = int(math.floor(max(nx.max(), L + 1e-12)))
F = F or 10 * d
if F2 is None: # forward tail past F2 - F below 1e-3
tail = np.cumsum(np.abs(hf)[::-1])[::-1]
F2 = min(len(hf) - 1, F + int(np.argmax(tail < 1e-3)))
beta = hb[1] / hb[2]
w = t is None
if w: # h_(m - u), m in [-B, F2)
lo_m = -B
def row(u):
m = np.arange(lo_m, F2)
j = m - u
v = np.where(j >= 0, hf[np.clip(j, 0, len(hf) - 1)],
hb[np.clip(-j, 0, len(hb) - 1)])
v[j >= len(hf)] = 0
v[-j >= len(hb)] = 0
full = np.abs(hf).sum() + np.abs(hb[1:]).sum()
return v, full - np.abs(v).sum()
us = range(-B + 1, F - d)
u0 = 0
else: # g^(u)_k, k in [0, F2)
lo_m = 0
K = np.arange(F2)
def shifted(s):
j = K - s
return np.where(j >= 0, hf[np.clip(j, 0, len(hf) - 1)],
hb[np.clip(-j, 0, len(hb) - 1)])
gam = -sum(shifted(s) for s in range(d)) / sum(beta ** -s
for s in range(d))
def row(u):
v = shifted(u) + beta ** -u * gam
# mass of g^(u) beyond F2: h's forward tail past F2 - u, plus
# gam's, bounded by d times h's tail past F2 - u - d
tail = np.abs(hf[max(0, F2 - u - d):]).sum() * (1 + d *
beta ** -u)
return v, tail
us = range(0, F - d)
u0 = t
M = F2 - lo_m # delta variables
us = list(us)
zi = {u: M + i for i, u in enumerate(us)}
wc = M + len(zi) if w else None
wb = beta ** -B * d * L / sum(beta ** -s for s in range(d))
A, lo, hi = [], [], []
r = 0
for u in us:
v, tau = row(u)
small = np.abs(v) < 1e-9 # into the slack, for conditioning
tau += np.abs(v[small]).sum()
A += [(r, int(c), float(v[c])) for c in np.flatnonzero(~small)]
A.append((r, zi[u], -1.0))
if w:
A.append((r, wc, beta ** -u))
lo.append(-tau - TOL)
hi.append(tau + TOL)
r += 1
for u in us: # z_u - sum z_(u+i) = delta_u
if u + d in zi and lo_m <= u < F:
A.append((r, zi[u], 1.0))
A += [(r, zi[u + i], -1.0) for i in range(1, d + 1)]
A.append((r, u - lo_m, -1.0))
lo.append(0)
hi.append(0)
r += 1
if not w: # the first window sum 0
A += [(r, zi[u], 1.0) for u in range(d)]
lo.append(0)
hi.append(0)
r += 1
return dict(A=A, lo=lo, hi=hi, M=M, F=F, lo_m=lo_m, zi=zi, X=X, u0=u0,
ext=[(-wb, wb, 0)] if w else [])
def solve(mo, c, tlim, presolve, lp):
"""HiGHS on the model's rows and its ext columns (lp: no integrality)."""
M, zi, X = mo["M"], mo["zi"], mo["X"]
nz = len(zi)
ext = mo["ext"]
nv = M + nz + len(ext)
rr, cc, vv = zip(*mo["A"])
A = coo_matrix((vv, (rr, cc)), shape=(len(mo["lo"]), nv)).tocsr()
integ = np.zeros(nv)
if not lp:
integ[:mo["F"] - mo["lo_m"]] = 1
integ[M:M + nz] = 1
lb = [-1] * M + [-X] * nz + [e[0] for e in ext]
ub = [1] * M + [X] * nz + [e[1] for e in ext]
return milp(c=c, integrality=integ, bounds=Bounds(lb, ub),
constraints=LinearConstraint(A, mo["lo"], mo["hi"]),
options={"time_limit": tlim, "presolve": presolve})
def feasible(d, t, kind, target=None, B=B, tlim=600, presolve=False,
lp=False):
"""A carry (or window sum) of target at t: t None reads every t >= B
at once, else the absolute problem at t. Zero objective; HiGHS's
status, 2 (infeasible) proving the bound there."""
target = target or C_of(d) + 1
mo = model(d, t, B=B)
r = len(mo["lo"])
mo["A"].append((r, mo["zi"][mo["u0"]], 1.0))
if kind == "sum":
mo["A"].append((r, mo["u0"] - mo["lo_m"], -1.0))
mo["lo"].append(target)
mo["hi"].append(np.inf)
nv = mo["M"] + len(mo["zi"]) + len(mo["ext"])
res = solve(mo, np.zeros(nv), tlim, presolve, lp)
return res.status
def arm(d, t, kind, t0, **kw):
t1 = time.time()
st = feasible(d, t, kind, **kw)
print(f" d = {d}, t {'>= %d' % B if t is None else '= %d' % t}, "
f"{kind} >= {kw.get('target') or C_of(d) + 1}"
f"{' (LP)' if kw.get('lp') else ''}: status {st} "
f"{time.time() - t1:.1f} s [{time.time() - t0:.1f} s]",
flush=True)
return st
def section_control(t0):
print("CONTROL d = 13, target 3: the wave a point of every arm",
flush=True)
sts = [arm(13, t, kind, t0, target=3)
for t, kind in ((None, "carry"), (None, "sum"), (3, "carry"),
(2, "sum"))]
check("every control arm feasible", all(s == 0 for s in sts),
f"statuses {sts}")
def section_f(t0):
print(f"F every arm at d = 11, 12, 16, 20: t >= {B}, then t < {B}",
flush=True)
for d in (11, 12, 16, 20):
sts = [arm(d, t, kind, t0) for t in [None] + list(range(B))
for kind in ("carry", "sum")]
check(f"d = {d}: no carry or window sum of {C_of(d) + 1}",
all(s == 2 for s in sts), f"statuses {sorted(set(sts))}")
if d == 11:
st = arm(11, None, "carry", t0, lp=True)
check("d = 11: the stationary carry's LP at 3 feasible "
"(integrality refutes it, not the box)", st == 0)
def main():
t0 = time.time()
section_control(t0)
section_f(t0)
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed "
f"[{time.time() - t0:.1f} s]")
sys.exit(0 if all(CHECKS) else 1)
if __name__ == "__main__":
main()