TJC-1814: Python 3.14 default; numba-jitted vertex enumeration (4.2x certify-sweep)

.github/workflows/ci.yml

@@ -39,7 +39,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python: ["3.11", "3.13"]
+        python: ["3.11", "3.14"]
     steps:
       - uses: actions/checkout@11d5960a326750d5838078e36cf38b85af677262 # v4
       - uses: astral-sh/setup-uv@c771a70e6277c0a99b617c7a806ffedaca235ff9
@@ -58,7 +58,7 @@ jobs:
       - uses: astral-sh/setup-uv@c771a70e6277c0a99b617c7a806ffedaca235ff9
         with:
           version: "0.11.6"
-          python-version: "3.13"
+          python-version: "3.14"
           enable-cache: true
       - run: uv run --with cffconvert==2.0.0 cffconvert --validate
 

.github/workflows/release.yml

@@ -54,7 +54,7 @@ jobs:
       - uses: astral-sh/setup-uv@c771a70e6277c0a99b617c7a806ffedaca235ff9
         with:
           version: "0.11.6"
-          python-version: "3.13"
+          python-version: "3.14"
           enable-cache: true
       - run: uv sync --frozen
       - run: uv run --with cffconvert==2.0.0 cffconvert --validate

.github/workflows/scientific.yml

@@ -23,7 +23,7 @@ jobs:
       - uses: astral-sh/setup-uv@c771a70e6277c0a99b617c7a806ffedaca235ff9
         with:
           version: "0.11.6"
-          python-version: "3.13"
+          python-version: "3.14"
           enable-cache: true
       - run: uv sync --frozen
       - run: uv run --frozen pytest -m slow

.python-version

@@ -1 +1 @@
-3.13
+3.14

README.md

@@ -37,7 +37,7 @@ exploratory computation.
 | `verify-goemans` | Exhaustive integer enumeration | Exact 58-versus-60 counterexample and exact violation-budget thresholds |
 | `amplify-family` | `fractions.Fraction` | Exact checks for the displayed integral families |
 | `revision-checks` | Exhaustive rational enumeration | Exact checks for the general-family examples, gate collapse, and dummy-terminal embedding |
-| `certify-sweep` | Seeded NumPy floating-point vertex enumeration and local search | Numerical stress test on $H^{\ast}$ plus a sampled four-terminal side search; boundary points and cases with more than nine admissible routing vectors are skipped, so this is not a proof certificate |
+| `certify-sweep` | Seeded floating-point vertex enumeration (numba-jitted) and local search | Numerical stress test on $H^{\ast}$ plus a sampled four-terminal side search; boundary points and cases with more than nine admissible routing vectors are skipped, so this is not a proof certificate |
 | `interval-search` | Exhaustive finite structure enumeration plus truncated assignment enumeration and floating-point grid/refinement | Exploratory records in the stated chain subclass; at most the first 16 minimal-transit assignments per structure are examined |
 | `laned-search` | Seeded floating-point heuristic over every fifth canonical structure | Exploratory records in the stated laned subclass; a deterministic one-fifth structure subsample |
 
@@ -107,7 +107,7 @@ breaks the build.
 ### Python
 
 [uv](https://docs.astral.sh/uv/) 0.11.6 or later is recommended; Python 3.11
-through 3.13 is supported.
+through 3.14 is supported, and 3.14 is the default.
 
 ```bash
 uv sync --frozen

paper/goemans_gap_paper.typ

@@ -316,7 +316,7 @@ _Proof._ The splits are integral and feasible: $(a - b) M - 1 >= 0$ since $a > b
 
 #text(weight: "bold")[Concurrent public family.] In the public discussion following Rybin's announcement, H. Juang posted a three-parameter integral counterexample family on the same digraph — Rybin's instance is its member $(b, m, g) = (10, 5, 1)$ up to scaling — and observed that parameters near $m\/b = sqrt(2) - 1$ push the ratio toward $(4 sqrt(2) + 2)\/7$; the reparametrization $a = b + m$ identifies that direction with $delta -> 1\/sqrt(2)$ here, and Juang's example $(70, 29, 1)$ has the demands and unsplittable optimum $2 a b$ of Proposition 4.1's Pell member $(a, b) = (99, 70)$ with a slightly different split. Beyond the family itself, the contributions here are the closed-form ratio, the exact limit, and — with Theorem A's upper bound — that the common value $rho^*$ is the exact, non-attained supremum.
 
-#text(weight: "bold")[Machine verification (verification, not proof).] The formulas of Proposition 4.1 were checked by exhaustive enumeration in exact rational arithmetic (#raw("fractions.Fraction")) for $(a, b)$ in ${(3, 2), (7, 5), (10, 7), (17, 12), (41, 29), (99, 70), (577, 408)}$ and $M$ in ${1, 1000}$. The complete path set of $H^*$ is kernel-checked in Lean (#raw("walk_s_terminal_iff"), Appendix A) and independently re-derived by depth-first search in the Rybin-instance verifier. Separately, a seeded floating-point sweep over $(d, f)$, using NumPy dual-vertex enumeration for each sampled cost problem (4000 random samples plus local search), never exceeded $rho^*$ and approached the exhibited family's boundary configuration $d = (1, 0.70711, 1)$, $f = (0.29289, 0.41421, 0.29289)$. Candidate points with any fractional arc load at or below $10^(-12)$ are excluded because the numerical dual formulation divides by those loads. This sweep is a numerical stress test, not an exact LP certificate. An exact integer violation-budget sweep shows the counterexamples analyzed here are knife-edge at budget $d_max$: at budget $d_max + 1$ (in Rybin's integers) the minimum cost collapses below the fractional cost. They therefore say nothing against the still-open $2 d_max$ or $O(d_max)$ cost conjectures, nor against the Morell-Skutella two-sided conjecture.
+#text(weight: "bold")[Machine verification (verification, not proof).] The formulas of Proposition 4.1 were checked by exhaustive enumeration in exact rational arithmetic (#raw("fractions.Fraction")) for $(a, b)$ in ${(3, 2), (7, 5), (10, 7), (17, 12), (41, 29), (99, 70), (577, 408)}$ and $M$ in ${1, 1000}$. The complete path set of $H^*$ is kernel-checked in Lean (#raw("walk_s_terminal_iff"), Appendix A) and independently re-derived by depth-first search in the Rybin-instance verifier. Separately, a seeded floating-point sweep over $(d, f)$, using jit-compiled dual-vertex enumeration (numba) for each sampled cost problem (4000 random samples plus local search), never exceeded $rho^*$ and approached the exhibited family's boundary configuration $d = (1, 0.70711, 1)$, $f = (0.29289, 0.41421, 0.29289)$. Candidate points with any fractional arc load at or below $10^(-12)$ are excluded because the numerical dual formulation divides by those loads. This sweep is a numerical stress test, not an exact LP certificate. An exact integer violation-budget sweep shows the counterexamples analyzed here are knife-edge at budget $d_max$: at budget $d_max + 1$ (in Rybin's integers) the minimum cost collapses below the fractional cost. They therefore say nothing against the still-open $2 d_max$ or $O(d_max)$ cost conjectures, nor against the Morell-Skutella two-sided conjecture.
 
 // ============================================================
 // Part II
@@ -509,7 +509,7 @@ where #raw("hStarRatios") is the parametrized ratio set and #raw("graphRatios")
 
 #raw("verify_goemans.py") exhaustively verifies Rybin's instance with exact integer arithmetic, derives all paths by depth-first search, and computes the exact violation-budget sweep. #raw("amplify_family.py") and #raw("revision_checks.py") use exact rational arithmetic for the displayed family members, the gate-gadget collapse, and the dummy-terminal embedding.
 
-The seeded floating-point scripts are evidence, not proof certificates, and their exact budgets and filters are part of the reported scope. #raw("certify_sweep.py") excludes points having any fractional arc load at or below $10^(-12)$. On $H^*$ it evaluates one boundary-near point with $epsilon = 10^(-6)$, 4000 random points, and two 4000-proposal hill climbs. Its four-terminal side search samples 120 of 1248 eligible designs and 120 parameter points per retained design, skips parameter points with more than nine admissible routing load vectors, and subjects only the best retained design to 6000 additional refinement proposals. #raw("interval_search.py") visits all 32, 112, 257, and 973 canonical structures in the four displayed classes, respectively, but retains at most the first 16 Cartesian-product assignments of inclusion-minimal transit sets per structure. It uses the four-value demand grid stated in Section 6 and refines only the best candidate in each class, with 4000 continuous and 300 exit-weight proposals. #raw("laned_search.py") enumerates 927 canonical structures, evaluates the 185 whose deterministic indices are divisible by five, and uses six starts totaling 3000 proposals for each retained structure. Complete outputs are committed with the repository, and the exact scripts are regression-tested byte for byte. The slow numerical tests compare normalized scientific records with those outputs and assert stable guardrails while allowing platform-dependent tie breaking and least-significant digits.
+The seeded floating-point scripts are evidence, not proof certificates, and their exact budgets and filters are part of the reported scope. #raw("certify_sweep.py") jit-compiles its dual-vertex enumeration with numba (partial-pivot LU with the same nonsingularity and feasibility thresholds as the batched formulation it replaces; the archived transcript is reproduced byte for byte) and excludes points having any fractional arc load at or below $10^(-12)$. On $H^*$ it evaluates one boundary-near point with $epsilon = 10^(-6)$, 4000 random points, and two 4000-proposal hill climbs. Its four-terminal side search samples 120 of 1248 eligible designs and 120 parameter points per retained design, skips parameter points with more than nine admissible routing load vectors, and subjects only the best retained design to 6000 additional refinement proposals. #raw("interval_search.py") visits all 32, 112, 257, and 973 canonical structures in the four displayed classes, respectively, but retains at most the first 16 Cartesian-product assignments of inclusion-minimal transit sets per structure. It uses the four-value demand grid stated in Section 6 and refines only the best candidate in each class, with 4000 continuous and 300 exit-weight proposals. #raw("laned_search.py") enumerates 927 canonical structures, evaluates the 185 whose deterministic indices are divisible by five, and uses six starts totaling 3000 proposals for each retained structure. Complete outputs are committed with the repository, and the exact scripts are regression-tested byte for byte. The slow numerical tests compare normalized scientific records with those outputs and assert stable guardrails while allowing platform-dependent tie breaking and least-significant digits.
 
 #sechead[References]
 

pyproject.toml

@@ -5,7 +5,7 @@ version = "0.1.0"
 description = "Reproducibility suite for the Rybin-digraph extremal theorem for Goemans's cost conjecture"
 readme = "README.md"
 authors = [{ name = "A. Richard Caputo III", email = "rcaputo3@tjclp.com" }]
-requires-python = ">=3.11,<3.14"
+requires-python = ">=3.11,<3.15"
 license = "MIT"
 license-files = ["LICENSE"]
 keywords = ["combinatorial-optimization", "formal-verification", "lean4", "unsplittable-flow"]
@@ -16,7 +16,11 @@ classifiers = [
     "Topic :: Scientific/Engineering :: Mathematics",
 ]
 dependencies = [
-    "numpy>=2.0",
+    # numba 0.66 requires numpy < 2.5; the jit-compiled vertex enumeration in
+    # certify_sweep reproduces the archived transcript byte-for-byte and runs
+    # the sweep about 4x faster than the batched-LAPACK formulation.
+    "numba>=0.66",
+    "numpy>=2.0,<2.5",
 ]
 
 [project.urls]

src/goemans_gap/certify_sweep.py

@@ -24,6 +24,7 @@
 import math
 import random
 
+import numba
 import numpy as np
 
 RHO = 2 / (2 * math.sqrt(2) - 1)
@@ -34,8 +35,103 @@
 MAX_TEMPLATE_ROUTING_VECTORS = 9
 
 
+@numba.njit(cache=True, fastmath=False)
+def _lp_value_kernel(matrix, rhs, routing_count):
+    """Fused dual-LP vertex enumeration.
+
+    Enumerates basis-row combinations in lexicographic order (matching the
+    previous ``itertools.combinations`` formulation), factors each
+    ``k x k`` basis system by partial-pivot LU with the same ``|det| > 1e-10``
+    nonsingularity screen, solves by back substitution, and applies the same
+    ``1e-9`` feasibility slack against every constraint row. Returns ``nan``
+    when no basis is nonsingular and ``inf`` when none is feasible.
+    """
+    rows, k = matrix.shape
+    m = k - 1
+    idx = np.empty(m, dtype=np.int64)
+    for i in range(m):
+        idx[i] = i
+    a = np.empty((k, k))
+    b = np.empty(k)
+    v = np.empty(k)
+    best = np.inf
+    found_nonsingular = False
+    while True:
+        for i in range(m):
+            row = idx[i]
+            for j in range(k):
+                a[i, j] = matrix[row, j]
+            b[i] = rhs[row]
+        for j in range(k):
+            a[m, j] = 1.0 if j < routing_count else 0.0
+        b[m] = 1.0
+
+        det = 1.0
+        singular = False
+        for col in range(k):
+            pivot_row = col
+            pivot_mag = abs(a[col, col])
+            for r in range(col + 1, k):
+                mag = abs(a[r, col])
+                if mag > pivot_mag:
+                    pivot_mag = mag
+                    pivot_row = r
+            if pivot_row != col:
+                for j in range(k):
+                    swap = a[col, j]
+                    a[col, j] = a[pivot_row, j]
+                    a[pivot_row, j] = swap
+                swap = b[col]
+                b[col] = b[pivot_row]
+                b[pivot_row] = swap
+                det = -det
+            pivot = a[col, col]
+            det *= pivot
+            if pivot == 0.0:
+                singular = True
+                break
+            inv = 1.0 / pivot
+            for r in range(col + 1, k):
+                factor = a[r, col] * inv
+                if factor != 0.0:
+                    for j in range(col, k):
+                        a[r, j] -= factor * a[col, j]
+                    b[r] -= factor * b[col]
+
+        if not singular and abs(det) > 1e-10:
+            found_nonsingular = True
+            for r in range(k - 1, -1, -1):
+                s = b[r]
+                for j in range(r + 1, k):
+                    s -= a[r, j] * v[j]
+                v[r] = s / a[r, r]
+            feasible = True
+            for r in range(rows):
+                s = 0.0
+                for j in range(k):
+                    s += matrix[r, j] * v[j]
+                if s > rhs[r] + 1e-9:
+                    feasible = False
+                    break
+            if feasible and v[routing_count] < best:
+                best = v[routing_count]
+
+        i = m - 1
+        while i >= 0 and idx[i] == rows - m + i:
+            i -= 1
+        if i < 0:
+            break
+        idx[i] += 1
+        for j in range(i + 1, m):
+            idx[j] = idx[j - 1] + 1
+
+    if not found_nonsingular:
+        return np.nan
+    return best
+
+
 def lp_value(x, routings):
-    """Approximate max-cost value by NumPy dual-LP vertex enumeration."""
+    """Approximate max-cost value by jit-compiled dual-LP vertex enumeration."""
     routing_count, arc_count = routings.shape
     # variables v = (lambda_1..K, t); constraints:
     #   G v <= h : for each arc a: sum_k Y[k,a] lam_k - t x_a <= 0
@@ -49,36 +145,10 @@ def lp_value(x, routings):
     matrix[arc_count : arc_count + routing_count, :routing_count] = -np.eye(
         routing_count
     )
-    equality = np.zeros(variable_count)
-    equality[:routing_count] = 1.0
-    combinations = np.array(
-        list(
-            itertools.combinations(
-                range(arc_count + routing_count), variable_count - 1
-            )
-        )
-    )
-    systems = np.concatenate(
-        [
-            matrix[combinations],
-            np.broadcast_to(equality, (len(combinations), 1, variable_count)),
-        ],
-        axis=1,
-    )
-    system_rhs = np.concatenate(
-        [rhs[combinations], np.ones((len(combinations), 1))], axis=1
-    )
-    determinants = np.abs(np.linalg.det(systems))
-    nonsingular = determinants > 1e-10
-    if not nonsingular.any():
-        return None
-    vertices = np.linalg.solve(
-        systems[nonsingular], system_rhs[nonsingular][..., None]
-    )[..., 0]
-    feasible = np.all(vertices @ matrix.T <= rhs + 1e-9, axis=1)
-    if not feasible.any():
+    result = _lp_value_kernel(matrix, rhs, routing_count)
+    if math.isnan(result) or math.isinf(result):
         return None
-    return float(vertices[feasible][:, routing_count].min())
+    return float(result)
 
 
 # ---------- Rybin digraph ----------