"""Exponential-space #Dominating Set: the analysis of van Rooij's thesis, Section 8.5, rerun with our pathwidth bounds (Section 13 of the paper). Source: J. M. M. van Rooij, Exact Exponential-Time Algorithms for Domination Problems in Graphs, PhD thesis, Utrecht University, 2011, Section 8.5: Algorithm 8.2 (p. 148), Overview 8.1 (p. 150), the weights and constraints (p. 151), Lemma 8.11 (p. 152), Lemma 8.12 (pp. 153-154), Theorem 8.13 (p. 154). The transcription follows the thesis. The journal version (Nederlof, van Rooij, van Dijk, Algorithmica 69 (2014), Section 7.1, pp. 724-731: Algorithm 4, Lemmas 29-30, Theorem 31) was compared later; it has the same algorithm, weights and analysis (results page, Section 2). Notation as in the thesis: v(i) is the weight of an element of frequency i, w(i) the weight of a set of size i; weights are constant from i = K = 7 on. A. The published analysis. The linear program of Lemma 8.12 is reproduced exactly (0.325782). The recurrences of Lemma 8.11 give alpha = 1.260893 for every case except branching on a set of size 4, where the published weights give 1.27583 (a set of size 4 whose four elements have frequency 4). Algorithm 8.2 branches on such sets (line 13). B. The model: Algorithm 8.2 with switching thresholds (dps, dpe) and a list of exceptional elements; the branching vectors of Lemma 8.11; the linear program of Lemma 8.12 with any set of per-degree pathwidth forms; and a weight optimizer (Kelley cutting planes for the branching constraints, LP duality for the DP condition, which is linear once alpha is fixed; the objective alpha^(v_max + w_max) does not depend on the scale, so alpha = 1.3 is fixed). C. Floating-point optima: the published algorithm repaired (all cases, FGSS coefficients), and variants with our pathwidth bounds (Theorems 2-4 and Lemma F1 of the results page). D. Certificates: rational weights; linear rows exact; every branching vector up to degree 9 by mpmath interval arithmetic (degree >= 10 reduces to degree 9, see the results page); the DP linear program by an exact rational dual certificate. Theorem numbers in the comments and labels follow the results page; in the manuscript: Theorem 4 (S) is thm:main, Theorem 3 (R), (P) are thm:deg5 with the triples R, P, Theorem 2(a), (b) are thm:deg4(a), (b), Theorem 1 is thm:cubic, and Lemma F1 is lem:n3; the third certificate of part D is the one of thm:ds. Run from the repository root: .venv/bin/python experiments/algorithmic-consequences/nrd_dominating_set.py """ import itertools import math from fractions import Fraction as Fr import numpy as np from mpmath import iv, mpf from scipy.optimize import linprog, brentq iv.prec = 120 fails = [] # ---------------------------------------------------------------- published data (thesis pp. 150-151) V_PUB = [0.000000, 0.275603, 0.551206, 0.658484, 0.719009, 0.745509, 0.749538] # v(1..7) W_PUB = [0.354816, 0.625964, 0.785112, 0.923325, 0.986711, 1.000000, 1.000000] # w(1..7) # Exceptional elements of frequency 5 (Overview 8.1): (r1,...,r5), r_i = number of neighbouring sets of size i. EXCEPTIONAL = [(1, 4, 0, 0, 0), (1, 3, 1, 0, 0), (1, 2, 2, 0, 0), (1, 1, 3, 0, 0), (1, 0, 4, 0, 0), (1, 3, 0, 1, 0), (1, 2, 1, 1, 0), (0, 5, 0, 0, 0), (0, 4, 1, 0, 0), (0, 3, 2, 0, 0), (0, 2, 3, 0, 0), (0, 1, 4, 0, 0), (0, 4, 0, 1, 0), (0, 3, 1, 1, 0)] # A larger list found by a greedy search for the design dps = dpe = 4 (part C); all tuples have r1 = 1 or are # published ones, so every exceptional set has size at most 5. TUNED = [(0, 1, 4, 0, 0), (0, 2, 3, 0, 0), (0, 3, 1, 1, 0), (0, 3, 2, 0, 0), (0, 4, 0, 1, 0), (0, 4, 1, 0, 0), (0, 5, 0, 0, 0), (1, 0, 0, 4, 0), (1, 0, 1, 3, 0), (1, 0, 2, 1, 1), (1, 0, 2, 2, 0), (1, 0, 3, 0, 1), (1, 0, 3, 1, 0), (1, 0, 4, 0, 0), (1, 1, 0, 2, 1), (1, 1, 0, 3, 0), (1, 1, 1, 1, 1), (1, 1, 1, 2, 0), (1, 1, 2, 0, 1), (1, 1, 2, 1, 0), (1, 1, 3, 0, 0), (1, 2, 0, 1, 1), (1, 2, 0, 2, 0), (1, 2, 1, 0, 1), (1, 2, 1, 1, 0), (1, 2, 2, 0, 0), (1, 3, 0, 0, 1), (1, 3, 0, 1, 0), (1, 3, 1, 0, 0), (1, 4, 0, 0, 0)] # Pathwidth forms {degree: coefficient}, each an upper bound (up to o(n)) for every graph of maximum degree <= 6 # after deleting the sets of size one. Key "n": coefficient on every remaining vertex. Coefficients are rationals # at least the true constants. F_UP = Fr(140159, 10 ** 6) # F ~ 0.1401583 (Theorem 2(b), triple P) F4_UP = Fr(24577197475, 10 ** 11) # F_4 <= 0.24577197475 (Theorem 2(a)) C3_UP = Fr(139823, 10 ** 6) # c_3 ~ 0.1398219 (Theorem 1 via Lemma F1) FGSS = [{3: Fr(1, 6), 4: Fr(1, 3), 5: Fr(13, 30), 6: Fr(23, 45)}] OURS = [{3: Fr(5, 32), 4: Fr(1, 4), 5: Fr(11, 32), 6: Fr(7, 16)}, # Theorem 4 (S) {3: Fr(1455, 10 ** 4), 4: Fr(2530, 10 ** 4), 5: Fr(3600, 10 ** 4), 6: Fr(1)}, # Theorem 3 (R); degree 6 deleted {3: F_UP, 4: Fr(2575, 10 ** 4), 5: Fr(3667, 10 ** 4), 6: Fr(1)}, # Theorem 3 (P) {"n": F4_UP, 3: Fr(0), 4: Fr(0), 5: 1 - F4_UP, 6: 1 - F4_UP}, # Theorem 2(a); degrees >= 5 deleted {3: F_UP, 4: Fr(2575, 10 ** 4), 5: Fr(1), 6: Fr(1)}, # Theorem 2(b) {3: C3_UP, 4: Fr(1), 5: Fr(1), 6: Fr(1)}] # Lemma F1; degrees >= 4 deleted def tau(a, b): """Root x > 1 of x^-a + x^-b = 1.""" if a <= 0 or b <= 0: return math.inf return brentq(lambda x: x ** -a + x ** -b - 1, 1 + 1e-12, 64) def multisets(n_items, types): """All vectors r indexed by types with sum r = n_items, as dicts.""" for combo in itertools.combinations_with_replacement(types, n_items): r = {} for t in combo: r[t] = r.get(t, 0) + 1 yield r # ---------------------------------------------------------------- B. the model and the weight optimizer def evalf(form, p): return sum(p[k] * x for k, x in form.items()) class Design: """Algorithm 8.2 with switching thresholds: dynamic programming when every non-exceptional set has size <= dps and every non-exceptional element has frequency <= dpe; exceptional elements have frequency dpe + 1 and a neighbourhood in `exc`; exceptional sets are sets of size > dps that contain an exceptional element. Branch on a set s if d(s) > d(e) or d(e) <= dpe, else on an element e (s, e of maximum degree, non-exceptional). forms: per-degree pathwidth forms {degree: coefficient}; `n` as a key means a coefficient on every vertex of degree >= 2 (after deleting sets of size one).""" def __init__(self, dps, dpe, exc, forms, K=7, dmax=9): self.dps, self.dpe, self.exc, self.forms, self.K, self.dmax = dps, dpe, list(exc), forms, K, dmax self.exc_set_sizes = sorted({i + 1 for c in self.exc for i, k in enumerate(c) if k > 0 and i + 1 > dps}) # weights are a vector p = (v(2..K), w(1..K)); helpers return linear forms as dicts {index: coeff} def vi(self, i): i = min(i, self.K) return {} if i <= 1 else {i - 2: 1.0} def wi(self, i): i = min(i, self.K) return {} if i <= 0 else {self.K - 1 + i - 1: 1.0} @property def nw(self): return 2 * self.K - 1 @staticmethod def add(*terms): out = {} for c, f in terms: for k, x in f.items(): out[k] = out.get(k, 0.0) + c * x return out def dv(self, i): return self.add((1, self.vi(i)), (-1, self.vi(i - 1))) def dw(self, i): return self.add((1, self.wi(i)), (-1, self.wi(i - 1))) def vectors(self): """Branching vectors as pairs of linear forms in p, with labels.""" out = [] dps, dpe = self.dps, self.dpe for d in range(dpe + 1, self.dmax + 1): # element branching, d(e) = d >= d(s) for r in multisets(d, range(1, d + 1)): if r.get(1, 0) > 1: continue if d == dpe + 1 and tuple(r.get(i, 0) for i in range(1, d + 1)) in self.exc: continue opt = self.add((1, self.vi(d)), *[(k, self.dw(i)) for i, k in r.items()]) inc = sum(k * (i - 1) for i, k in r.items()) forb = self.add((1, self.vi(d)), *[(k, self.wi(i)) for i, k in r.items()], (inc, self.dv(d))) out.append((opt, forb, ("element", d, dict(sorted(r.items()))))) for d in range(2, self.dmax + 1): # set branching, d(s) = d if d <= dps: if d - 1 <= dpe: continue emax = d - 1 else: emax = max(d - 1, dpe) D = max([d] + self.exc_set_sizes) for r in multisets(d, range(2, emax + 1)): disc = self.add((1, self.wi(d)), *[(k, self.dv(i)) for i, k in r.items()], (r.get(2, 0), self.wi(1))) inc = sum(k * (i - 1) for i, k in r.items()) take = self.add((1, self.wi(d)), *[(k, self.vi(i)) for i, k in r.items()], (inc, self.dw(D))) out.append((disc, take, ("set", d, dict(sorted(r.items()))))) return out def linear_constraints(self): """Rows g with g.p <= 0 (thesis p. 151, constraints 2-6): Delta v, Delta w >= 0; Delta v nonincreasing from i = 2 and Delta w from i = 1; 2 Delta v(dpe + 1) <= v(2), since elements are branched at frequency > dpe.""" rows = [] K = self.K for i in range(2, K + 1): rows.append(self.add((-1, self.dv(i)))) for i in range(1, K + 1): rows.append(self.add((-1, self.dw(i)))) for i in range(2, K): rows.append(self.add((1, self.dv(i + 1)), (-1, self.dv(i)))) for i in range(1, K): rows.append(self.add((1, self.dw(i + 1)), (-1, self.dw(i)))) rows.append(self.add((2, self.dv(self.dpe + 1)), (-1, self.vi(2)))) return rows # ---- the DP linear program (Lemma 8.12, generalized), in primal form for given weights def dp_columns(self): """Columns of the primal: sets x_i, elements y_i (non-exceptional, i <= dpe), exceptional e_c. Returns for each column: (measure form, size/frequency incidences (+ for sets, - for elements), degree after deleting size-one sets, kind).""" cols = [] smax = max([self.dps] + self.exc_set_sizes) for i in range(1, smax + 1): cols.append(("x", i, self.wi(i), i, 0 if i == 1 else i)) for i in range(2, self.dpe + 1): cols.append(("y", i, self.vi(i), -i, i)) for j, c in enumerate(self.exc): deg = self.dpe + 1 - (1 if c[0] > 0 else 0) cols.append(("e", j, self.vi(self.dpe + 1), -(self.dpe + 1), deg)) return cols def dp_rows(self): """Homogeneous rows: equalities (edge balance) and inequalities (<= 0) on the columns.""" cols = self.dp_columns() eq = [[c[3] for c in cols]] ub = [] smax = max([self.dps] + self.exc_set_sizes) for j in range(1, smax + 1): cnt = [0.0] * len(cols) for k, c in enumerate(cols): if c[0] == "e" and j <= len(self.exc[c[1]]): cnt[k] = self.exc[c[1]][j - 1] xj = next(k for k, c in enumerate(cols) if c[0] == "x" and c[1] == j) r = list(cnt); r[xj] -= j # sum_c r_cj e_c <= j x_j ub.append(r) if j > self.dps: # x_j <= sum_c r_cj e_c r = [-x for x in cnt]; r[xj] += 1 ub.append(r) return cols, eq, ub def form_value(self, form, col): """Coefficient of a form on one column (vertices of degree <= 2 cost nothing except under 'n').""" deg = col[4] if deg <= 1: return Fr(0) if deg == 2: return form.get("n", Fr(0)) val = form[deg] val += form.get("n", Fr(0)) return val def dp_value(self, p): cols, eq, ub = self.dp_rows() N = len(cols) meas = [sum(p[k] * x for k, x in c[2].items()) for c in cols] A_eq = [r + [0.0] for r in eq] + [meas + [0.0]] b_eq = [0.0] * len(eq) + [1.0] A_ub = [r + [0.0] for r in ub] for f in self.forms: A_ub.append([-float(self.form_value(f, c)) for c in cols] + [1.0]) b_ub = [0.0] * len(A_ub) res = linprog([0.0] * N + [-1.0], A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=[(0, None)] * N + [(None, None)], method="highs") sol = {f"{c[0]}{c[1]}": res.x[k] for k, c in enumerate(cols) if res.x[k] > 1e-9} return -res.fun, sol def optimize(design, alpha, iters=500, tol=1e-9, vecs=None, verbose=False): """Minimize v(K) + w(K) subject to all branching vectors at alpha, the steepness rows, and the DP condition (LP value <= log2 alpha), by Kelley's cutting planes on the convex constraints alpha^-a + alpha^-b <= 1. The DP condition enters through LP duality and is linear for fixed alpha.""" L = math.log2(alpha) la = math.log(alpha) if vecs is None: vecs = design.vectors() nw = design.nw VA = np.array([[a.get(k, 0.0) for k in range(nw)] for a, b, lab in vecs]) VB = np.array([[b.get(k, 0.0) for k in range(nw)] for a, b, lab in vecs]) cols, eq, ub = design.dp_rows() nf = len(design.forms) ne, nu_ = len(eq), len(ub) NV = nw + ne + nu_ + nf A_ub, b_ub, A_eq, b_eq = [], [], [], [] for g in design.linear_constraints(): r = [0.0] * NV for k, x in g.items(): r[k] += x A_ub.append(r); b_ub.append(0.0) # dual feasibility of the DP program: for each column k, # L*m_k(p) + sum_e mu_e eq_ek + sum_b nu_b ub_bk - sum_j lam_j F_jk >= 0, sum_j lam_j = 1 for k, c in enumerate(cols): r = [0.0] * NV for i, x in c[2].items(): r[i] -= L * x for e in range(ne): r[nw + e] -= eq[e][k] for b in range(nu_): r[nw + ne + b] -= ub[b][k] for j, f in enumerate(design.forms): r[nw + ne + nu_ + j] += float(design.form_value(f, c)) A_ub.append(r); b_ub.append(0.0) r = [0.0] * NV for j in range(nf): r[nw + ne + nu_ + j] = 1.0 A_eq.append(r); b_eq.append(1.0) bounds = [(0, 3)] * nw + [(None, None)] * ne + [(0, None)] * nu_ + [(0, None)] * nf obj = np.zeros(NV) obj[design.K - 2] = 1.0 # v(K) obj[design.K - 1 + design.K - 1] = 1.0 # w(K) A_ub = np.array(A_ub); b_ub = np.array(b_ub) cuts, rhss = [], [] p = np.ones(nw) obj_val = None def finish(p): """Scale p up until every branching vector holds; scaling keeps the linear rows and the DP dual valid.""" lo_s, hi_s = 1.0, 1.0 while (alpha ** -(VA @ (hi_s * p)) + alpha ** -(VB @ (hi_s * p))).max() > 1: hi_s *= 1 + 1e-6 return float(obj @ np.concatenate([hi_s * p, np.zeros(NV - nw)])), hi_s * p prev = None for it in range(iters): a, b = VA @ p, VB @ p ga, gb = alpha ** -a, alpha ** -b g = ga + gb bad = np.where(g > 1 + tol)[0] if it > 0 and (len(bad) == 0 or (prev is not None and abs(prev - obj_val) < 1e-12 and g.max() < 1 + 1e-6)): return finish(p) prev = obj_val bad = bad[np.argsort(-g[bad])][:80] for i in bad: r = np.zeros(NV) r[:nw] = -la * (ga[i] * VA[i] + gb[i] * VB[i]) cuts.append(r) rhss.append(1 - g[i] - la * (ga[i] * a[i] + gb[i] * b[i])) AU = np.vstack([A_ub, np.array(cuts)]) if cuts else A_ub res = linprog(obj, A_ub=AU, b_ub=np.concatenate([b_ub, rhss]), A_eq=A_eq, b_eq=b_eq, bounds=bounds, method="highs") if res.status != 0: return None, None p = res.x[:nw] obj_val = res.fun if verbose: print(it, obj_val, len(bad)) return None, None # ---------------------------------------------------------------- D. certificates def round_weights(design, p, grid=10 ** 7): """Rational weights from a float vector: round the increments, make them nonincreasing (Delta v from i = 2, Delta w from i = 1) and nonnegative, and sum them up.""" K = design.K v = [0.0, 0.0] + list(p[:K - 1]) # v(0..K) w = [0.0] + list(p[K - 1:]) # w(0..K) dv = [Fr(round((v[i] - v[i - 1]) * grid), grid) for i in range(2, K + 1)] dw = [Fr(round((w[i] - w[i - 1]) * grid), grid) for i in range(1, K + 1)] for seq in (dv, dw): for i in range(len(seq)): seq[i] = max(Fr(0), min(seq[i], seq[i - 1]) if i else seq[i]) vq = [sum(dv[:i], Fr(0)) for i in range(1, K)] # v(2..K) wq = [sum(dw[:i], Fr(0)) for i in range(1, K + 1)] return vq + wq def evalq(form, q): return sum((q[k] * Fr(x).limit_denominator(10 ** 6) for k, x in form.items()), Fr(0)) def certify(name, design, p, alpha=Fr(13, 10), dmax_check=9): """Certify alpha^(v(K) + w(K)) for rational weights derived from p. Returns the certified base (upper end).""" ok = True q = round_weights(design, p) dcheck = Design(design.dps, design.dpe, design.exc, design.forms, K=design.K, dmax=dmax_check) vecs = dcheck.vectors() # all coefficients of the branching forms are small integers; evaluate exactly ab = [(evalq(a, q), evalq(b, q)) for a, b, lab in vecs] # smallest scale s = 1 + j/10^6 for which every vector holds in floating point, then exactly by intervals fa = float(alpha) j = 0 while max(fa ** -(float(a) * (1 + j / 1e6)) + fa ** -(float(b) * (1 + j / 1e6)) for a, b in ab) > 1 - 1e-12: j += 1 s = Fr(10 ** 6 + j, 10 ** 6) q = [s * x for x in q] A = iv.mpf(alpha.numerator) / alpha.denominator worst_iv = 0 for (a, b) in ab: a, b = s * a, s * b g = A ** (-(iv.mpf(a.numerator) / a.denominator)) + A ** (-(iv.mpf(b.numerator) / b.denominator)) if not g.b <= 1: ok = False worst_iv = max(worst_iv, float(g.b)) # linear rows (homogeneous, so scaling keeps them) for g in design.linear_constraints(): if evalq(g, q) > 0: ok = False print(" linear row fails:", g) # DP: exact dual certificate. Solve the dual numerically, round, and compute the exact L it certifies. cols, eq, ub = design.dp_rows() eq = [[Fr(round(x)) for x in r] for r in eq] # all entries are integers ub = [[Fr(round(x)) for x in r] for r in ub] nf = len(design.forms) ne, nu_ = len(eq), len(ub) meas = [sum((q[k] * int(x) for k, x in c[2].items()), Fr(0)) for c in cols] Fm = [[design.form_value(f, c) for c in cols] for f in design.forms] # variables: L, mu (ne), nu (nu_), lam (nf); minimize L s.t. for each column k: # -L m_k - sum mu eq_k - sum nu ub_k + sum lam F_k <= 0; sum lam = 1 NV = 1 + ne + nu_ + nf A_ub, b_ub = [], [] for k in range(len(cols)): r = [-float(meas[k])] + [-float(eq[e][k]) for e in range(ne)] + [-float(ub[b][k]) for b in range(nu_)] + \ [float(Fm[j][k]) for j in range(nf)] A_ub.append(r); b_ub.append(0.0) A_eq = [[0.0] * (1 + ne + nu_) + [1.0] * nf] res = linprog([1.0] + [0.0] * (NV - 1), A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=[1.0], bounds=[(None, None)] * (1 + ne) + [(0, None)] * (nu_ + nf), method="highs") y = res.x mu = [Fr(y[1 + e]).limit_denominator(10 ** 9) for e in range(ne)] nu = [max(Fr(0), Fr(y[1 + ne + b]).limit_denominator(10 ** 9)) for b in range(nu_)] lam = [max(Fr(0), Fr(y[1 + ne + nu_ + j]).limit_denominator(10 ** 9)) for j in range(nf)] tot = sum(lam) lam = [x / tot for x in lam] Lc = max((sum(lam[j] * Fm[j][k] for j in range(nf)) - sum(mu[e] * eq[e][k] for e in range(ne)) - sum(nu[b] * ub[b][k] for b in range(nu_))) / meas[k] for k in range(len(cols))) dp_base = iv.mpf(2) ** (iv.mpf(Lc.numerator) / Lc.denominator) if not dp_base.b <= A.a: ok = False expo = q[design.K - 2] + q[2 * design.K - 2] # v(K) + w(K) base = A ** (iv.mpf(expo.numerator) / expo.denominator) vq = [0, 0] + q[:design.K - 1] wq = [0] + q[design.K - 1:] print(f" {name}: {'certified' if ok else 'FAILED'}: alpha = {alpha}, v_max + w_max = {float(expo):.7f}, " f"base < {float(base.b):.6f}") print(f" scale s = {s}; v(2..7) = {[str(x) for x in vq[2:]]}") print(f" w(1..7) = {[str(x) for x in wq[1:]]}") print(f" largest alpha^-a + alpha^-b (upper end) = {worst_iv:.9f} over {len(ab)} vectors up to degree " f"{dmax_check}; DP: 2^L < {float(dp_base.b):.9f} with L = {float(Lc):.9f} (exact dual certificate)") print(f" dual weights lambda of the forms (exact, in the order of the list of forms) = {[str(x) for x in lam]}") if not ok: fails.append(name) return float(base.b), q # ---------------------------------------------------------------- main if __name__ == "__main__": print("A. Published analysis (thesis Section 8.5)") pub = Design(3, 4, EXCEPTIONAL, FGSS, dmax=9) p_pub = np.array(V_PUB[1:] + W_PUB) scored = sorted(((tau(evalf(a, p_pub), evalf(b, p_pub)), lab) for a, b, lab in pub.vectors()), key=lambda t: -t[0]) rest = [t for t in scored if not (t[1][0] == "set" and t[1][1] == 4)] print(f" all cases except sets of size 4: alpha = {rest[0][0]:.6f} (published 1.26089); tight: " f"{[t[1] for t in rest[:3]]}") print(f" sets of size 4: alpha = {scored[0][0]:.6f} at {scored[0][1]}; next {scored[1][0]:.6f} at {scored[1][1]}") if not abs(rest[0][0] - 1.26089) < 5e-6: fails.append("published alpha (other cases)") if not scored[0][0] > 1.275: fails.append("size-4 set case") val, sol = pub.dp_value(list(p_pub)) print(f" dynamic programming (Lemma 8.12): {val:.6f}, 2^value = {2 ** val:.5f} (published 0.325782, 1.25334); " f"solution {({k: round(x, 6) for k, x in sol.items()})}") if not abs(val - 0.325782) < 2e-6: fails.append("published LP") print(f" Theorem 8.13: 1.26089^(v_max + w_max) = {1.26089 ** (V_PUB[-1] + W_PUB[-1]):.5f}") print("C. Floating-point optima (alpha fixed at 1.3; the objective does not depend on the scale)") AL = 1.3 designs = [("published algorithm, all cases, FGSS coefficients", Design(3, 4, EXCEPTIONAL, FGSS, dmax=8)), ("published algorithm, all cases, our bounds", Design(3, 4, EXCEPTIONAL, OURS, dmax=8)), ("published algorithm without size-4 set cases, FGSS (not a valid analysis)", None), ("published case range (sets of size d with elements of frequency < d), FGSS (not a valid analysis)", "published-range"), ("DP for sets <= 4, elements <= 4, published exceptions, FGSS", Design(4, 4, EXCEPTIONAL, FGSS, dmax=8)), ("DP for sets <= 4, elements <= 4, no exceptions, our bounds", Design(4, 4, [], OURS, dmax=8)), ("DP for sets <= 4, elements <= 4, published exceptions, Theorem 4 only", Design(4, 4, EXCEPTIONAL, OURS[:1], dmax=8)), ("DP for sets <= 4, elements <= 4, published exceptions, our bounds", Design(4, 4, EXCEPTIONAL, OURS, dmax=8)), ("DP for sets <= 4, elements <= 4, tuned exceptions, our bounds", Design(4, 4, TUNED, OURS, dmax=8)), ("DP for sets <= 3, elements <= 5, no exceptions, our bounds", Design(3, 5, [], OURS, dmax=8)), ("DP for sets <= 4, elements <= 5, no exceptions, our bounds", Design(4, 5, [], OURS, dmax=8)), ("DP for sets <= 5, elements <= 4, no exceptions, our bounds", Design(5, 4, [], OURS, dmax=8))] sols = {} for name, d in designs: if d == "published-range": # the case range of the published recurrences, which omits the set branchings where an element has # frequency at least the size of the set (see Appendix D.5 (`sec:nrd-gap`) of the paper: reproduces 1.50016) d0 = Design(3, 4, EXCEPTIONAL, FGSS, dmax=8) vecs = [x for x in d0.vectors() if not (x[2][0] == "set" and max(x[2][2]) >= x[2][1])] o, p = optimize(d0, AL, vecs=vecs) elif d is None: d0 = Design(3, 4, EXCEPTIONAL, FGSS, dmax=8) vecs = [x for x in d0.vectors() if not (x[2][0] == "set" and x[2][1] == 4)] o, p = optimize(d0, AL, vecs=vecs) else: o, p = optimize(d, AL) print(f" {AL ** o:.6f} {name}") sols[name] = (d, p) print("D. Certificates") for name in ["published algorithm, all cases, FGSS coefficients", "DP for sets <= 4, elements <= 4, published exceptions, our bounds", "DP for sets <= 4, elements <= 4, tuned exceptions, our bounds"]: d, p = sols[name] certify(name, d, p) print("E. Further values printed in Appendix D (floating point, for the certified weights of thm:ds as printed)") # The rational weights of the certificate in the committed output and in the paper: v(2..7), then w(1..7). # They are fixed here, since the optimizer of part C may return slightly different weights on other # versions of scipy. CERT_TUNED = [Fr(x) for x in ["2175171/10000000", "2175171/5000000", "2781121/5000000", "388797/625000", "6539177/10000000", "6655099/10000000", "1531123/5000000", "1067673/2000000", "1350949/2000000", "7616247/10000000", "8178813/10000000", "8324929/10000000", "8324929/10000000"]] d = Design(4, 4, TUNED, OURS, dmax=8) q = [float(x) for x in CERT_TUNED] val_s, _ = Design(d.dps, d.dpe, d.exc, OURS[:1], K=d.K, dmax=d.dmax).dp_value(q) print(f" DP value with the form of Theorem 4 (S) alone: {val_s:.7f} (log2 1.3 = {math.log2(1.3):.7f})") if not val_s > math.log2(1.3): fails.append("S alone above log2 alpha") d12 = Design(d.dps, d.dpe, d.exc, d.forms, K=d.K, dmax=12) qa = np.array(q) worst12, n12 = 0.0, 0 for a, b, lab in d12.vectors(): worst12 = max(worst12, 1.3 ** -evalf(a, qa) + 1.3 ** -evalf(b, qa)) n12 += 1 print(f" recheck of all branching vectors up to degree 12: {n12} vectors, largest alpha^-a + alpha^-b = " f"{worst12:.8f}") if not worst12 <= 1: fails.append("degree-12 recheck") print("FAILS:", fails if fails else "none")