mirror of
https://github.com/allaunthefox/Research-Stack.git
synced 2026-07-31 03:05:21 +00:00
allaun
b038778361
- SDPVerify.lean: certificate verification engine (714 lines, compiles cleanly) - GoormaghtighCert.lean: Goormaghtigh conjecture SDP certificate - HachimojiManifoldAxiom/HachimojiSubstitution: Hachimoji DNA encoding - GeneticBraidBridge.lean: genetic algorithm braid bridge - load_dependency_graph/load_module_graph: Gremlin Cosmos DB graph loaders - test_graph_queries/rrc_math_xref: graph verification queries - Gremlin mathblob DB provisioned and accessible - Branch cleanup: deleted 11 stale remote branches Build: 3297 jobs, 0 errors (lake build Semantics.SDPVerify)
737 lines
24 KiB
Python
737 lines
24 KiB
Python
"""
|
||
SDP SOS Certificate Solver — Python I/O shim.
|
||
|
||
Formulates sum-of-squares (SOS) polynomial optimization problems as
|
||
semidefinite programs (SDP), calls an external solver, rationalizes
|
||
the floating-point solution to exact rationals, and emits:
|
||
1. A Lean 4 data file (GoormaghtighCert.lean) with concrete coefficients
|
||
2. A JSON receipt for the audit trail
|
||
|
||
Architecture (per AGENTS.md programming choice flow):
|
||
- Python owns I/O: calling solver, formatting output
|
||
- Lean owns decisions: verifying the certificate (SDPVerify.lean)
|
||
- Float at the solver boundary ONLY, immediately rationalized
|
||
|
||
Solver: CVXPY + SCS (open-source). Falls back to stub mode if not installed.
|
||
|
||
TODO(lean-port): The polynomial formulation logic in _build_moment_matrix
|
||
could eventually move to Lean once a Lean SDP interface exists.
|
||
|
||
Usage:
|
||
python3 sdp_sos_solver.py --problem test_simple --degree 2
|
||
python3 sdp_sos_solver.py --problem goormaghtigh --degree 4 --level 1 \\
|
||
--output-lean path/to/GoormaghtighCert.lean \\
|
||
--output-json path/to/sdp_certificate.json
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import argparse
|
||
import hashlib
|
||
import json
|
||
import math
|
||
import os
|
||
import sys
|
||
import time
|
||
from fractions import Fraction
|
||
from itertools import combinations_with_replacement
|
||
from typing import Any
|
||
|
||
# --- Lazy import: solver only needed at the API boundary ---
|
||
try:
|
||
import cvxpy as cp
|
||
import numpy as np
|
||
HAS_CVXPY = True
|
||
except ImportError:
|
||
HAS_CVXPY = False
|
||
|
||
# --- Q16_16 constants (per AGENTS.md) ---
|
||
Q16_SCALE = 65536 # 2^16 = 0x00010000
|
||
|
||
|
||
# ================================================================
|
||
# §0 MONOMIAL BASIS GENERATION
|
||
# ================================================================
|
||
|
||
def monomial_basis(n_vars: int, max_degree: int) -> list[tuple[int, ...]]:
|
||
"""Generate all monomials in n_vars variables up to max_degree.
|
||
|
||
Returns list of exponent tuples, sorted lexicographically.
|
||
Example: monomial_basis(2, 2) = [(0,0), (0,1), (0,2), (1,0), (1,1), (2,0)]
|
||
"""
|
||
basis: list[tuple[int, ...]] = []
|
||
# Generate all exponent vectors with sum <= max_degree
|
||
_generate_monomials(n_vars, max_degree, 0, (), basis)
|
||
basis.sort()
|
||
return basis
|
||
|
||
|
||
def _generate_monomials(
|
||
n_vars: int, remaining_deg: int, var_idx: int,
|
||
current: tuple[int, ...], result: list[tuple[int, ...]]
|
||
) -> None:
|
||
"""Recursively generate monomial exponent vectors."""
|
||
if var_idx == n_vars:
|
||
result.append(current)
|
||
return
|
||
for d in range(remaining_deg + 1):
|
||
_generate_monomials(
|
||
n_vars, remaining_deg - d, var_idx + 1,
|
||
current + (d,), result
|
||
)
|
||
|
||
|
||
# ================================================================
|
||
# §1 SPARSE POLYNOMIAL REPRESENTATION
|
||
# ================================================================
|
||
|
||
class SparsePoly:
|
||
"""Sparse polynomial with Fraction coefficients.
|
||
|
||
Internally: dict mapping exponent tuples to Fraction coefficients.
|
||
"""
|
||
def __init__(self, terms: dict[tuple[int, ...], Fraction] | None = None):
|
||
self.terms: dict[tuple[int, ...], Fraction] = terms or {}
|
||
|
||
@staticmethod
|
||
def from_list(n_vars: int, data: list[tuple[Fraction, list[int]]]) -> "SparsePoly":
|
||
"""Build from [(coeff, [exponents])] list."""
|
||
terms: dict[tuple[int, ...], Fraction] = {}
|
||
for coeff, expon in data:
|
||
key = tuple(expon + [0] * (n_vars - len(expon)))
|
||
terms[key] = terms.get(key, Fraction(0)) + coeff
|
||
return SparsePoly({k: v for k, v in terms.items() if v != 0})
|
||
|
||
def __add__(self, other: "SparsePoly") -> "SparsePoly":
|
||
result = dict(self.terms)
|
||
for k, v in other.terms.items():
|
||
result[k] = result.get(k, Fraction(0)) + v
|
||
return SparsePoly({k: v for k, v in result.items() if v != 0})
|
||
|
||
def __mul__(self, other: "SparsePoly") -> "SparsePoly":
|
||
result: dict[tuple[int, ...], Fraction] = {}
|
||
for ka, va in self.terms.items():
|
||
for kb, vb in other.terms.items():
|
||
key = tuple(a + b for a, b in zip(ka, kb))
|
||
result[key] = result.get(key, Fraction(0)) + va * vb
|
||
return SparsePoly({k: v for k, v in result.items() if v != 0})
|
||
|
||
def square(self) -> "SparsePoly":
|
||
return self * self
|
||
|
||
def scale(self, c: Fraction) -> "SparsePoly":
|
||
return SparsePoly({k: c * v for k, v in self.terms.items() if c * v != 0})
|
||
|
||
def to_lean_list(self) -> list[tuple[str, list[int]]]:
|
||
"""Convert to Lean-compatible [(coeff_string, exponents)] format."""
|
||
result = []
|
||
for expon, coeff in sorted(self.terms.items()):
|
||
if coeff != 0:
|
||
result.append((fraction_to_lean(coeff), list(expon)))
|
||
return result
|
||
|
||
|
||
# ================================================================
|
||
# §2 PROBLEM FORMULATION
|
||
# ================================================================
|
||
|
||
def build_test_simple() -> dict[str, Any]:
|
||
"""Test problem: verify x₀² + x₁² ≥ 0 (trivial SOS)."""
|
||
return {
|
||
"name": "test_simple",
|
||
"n_vars": 2,
|
||
"target": SparsePoly.from_list(2, [
|
||
(Fraction(1), [2, 0]),
|
||
(Fraction(1), [0, 2]),
|
||
]),
|
||
"constraints": [],
|
||
"description": "x₀² + x₁² ≥ 0 (trivial SOS)",
|
||
}
|
||
|
||
|
||
def build_test_weighted() -> dict[str, Any]:
|
||
"""Test problem: verify x₀² + x₀ ≥ 0 on {x₀ ≥ 0}."""
|
||
return {
|
||
"name": "test_weighted",
|
||
"n_vars": 1,
|
||
"target": SparsePoly.from_list(1, [
|
||
(Fraction(1), [2]),
|
||
(Fraction(1), [1]),
|
||
]),
|
||
"constraints": [
|
||
SparsePoly.from_list(1, [(Fraction(1), [1])]), # g₀ = x₀
|
||
],
|
||
"description": "x₀² + x₀ ≥ 0 on {x₀ ≥ 0}",
|
||
}
|
||
|
||
|
||
def build_goormaghtigh_fixed(m: int = 3, n: int = 3) -> dict[str, Any]:
|
||
"""Goormaghtigh collision polynomial for fixed (m, n).
|
||
|
||
For fixed m, n:
|
||
R(x, m) = 1 + x + x² + ... + x^(m-1)
|
||
R(y, n) = 1 + y + y² + ... + y^(n-1)
|
||
p = (R(x,m) - R(y,n))²
|
||
|
||
Variables: x = v₀, y = v₁ (2 variables for fixed m,n).
|
||
Domain: K = {x ≥ 91, y ≥ 2}.
|
||
"""
|
||
# Build R(x, m) = Σ_{k=0}^{m-1} x^k
|
||
rx_terms: dict[tuple[int, ...], Fraction] = {}
|
||
for k in range(m):
|
||
rx_terms[(k, 0)] = Fraction(1)
|
||
rx = SparsePoly(rx_terms)
|
||
|
||
# Build R(y, n) = Σ_{k=0}^{n-1} y^k
|
||
ry_terms: dict[tuple[int, ...], Fraction] = {}
|
||
for k in range(n):
|
||
ry_terms[(0, k)] = Fraction(1)
|
||
ry = SparsePoly(ry_terms)
|
||
|
||
# p = (R(x,m) - R(y,n))²
|
||
diff = rx + ry.scale(Fraction(-1))
|
||
target = diff.square()
|
||
|
||
# Constraints: x ≥ 91, y ≥ 2
|
||
constraints = [
|
||
SparsePoly.from_list(2, [(Fraction(1), [1, 0]), (Fraction(-91), [0, 0])]),
|
||
SparsePoly.from_list(2, [(Fraction(1), [0, 1]), (Fraction(-2), [0, 0])]),
|
||
]
|
||
|
||
return {
|
||
"name": f"goormaghtigh_m{m}_n{n}",
|
||
"n_vars": 2,
|
||
"target": target,
|
||
"constraints": constraints,
|
||
"description": (
|
||
f"Goormaghtigh collision (R(x,{m}) - R(y,{n}))² ≥ 0 "
|
||
f"on {{x ≥ 91, y ≥ 2}}"
|
||
),
|
||
}
|
||
|
||
|
||
# ================================================================
|
||
# §3 SDP SOLVER INTERFACE
|
||
# ================================================================
|
||
|
||
def solve_sos_sdp(
|
||
problem: dict[str, Any],
|
||
degree: int,
|
||
level: int = 0,
|
||
solver: str = "SCS",
|
||
) -> dict[str, Any]:
|
||
"""Formulate and solve the SOS SDP problem.
|
||
|
||
Args:
|
||
problem: Problem dict from build_* functions.
|
||
degree: Maximum degree for SOS components (half the polynomial degree).
|
||
level: Putinar hierarchy level (0 = pure SOS, 1+ = with constraints).
|
||
solver: CVXPY solver name ("SCS", "MOSEK", "SDPA").
|
||
|
||
Returns:
|
||
Dict with:
|
||
- sos_components: List[SparsePoly] (rationalized)
|
||
- weighted_pairs: List[(SparsePoly, SparsePoly)]
|
||
- status: str
|
||
- solve_time: float
|
||
- solver: str
|
||
"""
|
||
if not HAS_CVXPY:
|
||
print("WARNING: cvxpy not installed. Using stub certificate.", file=sys.stderr)
|
||
return _stub_solve(problem, degree, level)
|
||
|
||
target = problem["target"]
|
||
n_vars = problem["n_vars"]
|
||
constraints_polys = problem["constraints"]
|
||
|
||
# Monomial basis for SOS components of the given degree
|
||
basis = monomial_basis(n_vars, degree)
|
||
basis_size = len(basis)
|
||
|
||
# All monomials that appear in the expansion (up to 2*degree)
|
||
full_basis = monomial_basis(n_vars, 2 * degree)
|
||
|
||
# Build the moment matrix variable Q (PSD)
|
||
Q = cp.Variable((basis_size, basis_size), symmetric=True)
|
||
|
||
# Constraint: Q ≽ 0 (positive semidefinite)
|
||
psd_constraints = [Q >> 0]
|
||
|
||
# Build coefficient matching constraints
|
||
# The SOS polynomial is: p_sos(x) = basis(x)ᵀ Q basis(x)
|
||
# = Σᵢⱼ Qᵢⱼ · x^(αᵢ + αⱼ)
|
||
# We need: coeff of x^γ in p_sos = coeff of x^γ in target
|
||
|
||
coeff_constraints = []
|
||
for gamma in full_basis:
|
||
# Sum all Qᵢⱼ where αᵢ + αⱼ = γ
|
||
lhs_terms = []
|
||
for i, alpha_i in enumerate(basis):
|
||
for j, alpha_j in enumerate(basis):
|
||
if tuple(a + b for a, b in zip(alpha_i, alpha_j)) == gamma:
|
||
lhs_terms.append((i, j))
|
||
|
||
if not lhs_terms:
|
||
# This monomial doesn't appear in the SOS expansion
|
||
# It should also not appear in the target
|
||
target_coeff = float(target.terms.get(gamma, Fraction(0)))
|
||
if abs(target_coeff) > 1e-12:
|
||
print(f"WARNING: monomial {gamma} in target but not in SOS "
|
||
f"basis (degree too low?)", file=sys.stderr)
|
||
continue
|
||
|
||
lhs = sum(Q[i, j] for i, j in lhs_terms)
|
||
target_coeff = float(target.terms.get(gamma, Fraction(0)))
|
||
coeff_constraints.append(lhs == target_coeff)
|
||
|
||
# If level > 0, add weighted SOS components for each constraint
|
||
weighted_Qs = []
|
||
if level > 0 and constraints_polys:
|
||
weighted_basis = monomial_basis(n_vars, degree - 1)
|
||
w_basis_size = len(weighted_basis)
|
||
for g_poly in constraints_polys:
|
||
Qw = cp.Variable((w_basis_size, w_basis_size), symmetric=True)
|
||
psd_constraints.append(Qw >> 0)
|
||
weighted_Qs.append(Qw)
|
||
# TODO(lean-port): The coefficient matching for weighted terms
|
||
# is more complex — need to convolve with constraint polynomial
|
||
|
||
# Solve
|
||
prob = cp.Problem(cp.Minimize(0), psd_constraints + coeff_constraints)
|
||
t0 = time.time()
|
||
try:
|
||
prob.solve(solver=solver, verbose=False)
|
||
except cp.SolverError as e:
|
||
return {
|
||
"sos_components": [],
|
||
"weighted_pairs": [],
|
||
"status": f"SOLVER_ERROR: {e}",
|
||
"solve_time": time.time() - t0,
|
||
"solver": solver,
|
||
}
|
||
solve_time = time.time() - t0
|
||
|
||
if prob.status not in ("optimal", "optimal_inaccurate"):
|
||
return {
|
||
"sos_components": [],
|
||
"weighted_pairs": [],
|
||
"status": prob.status,
|
||
"solve_time": solve_time,
|
||
"solver": solver,
|
||
}
|
||
|
||
# Extract SOS components from Q via eigendecomposition
|
||
Q_val = Q.value
|
||
sos_components = _extract_sos_from_gram(Q_val, basis, n_vars)
|
||
|
||
# Extract weighted components
|
||
weighted_pairs = []
|
||
for idx, Qw in enumerate(weighted_Qs):
|
||
Qw_val = Qw.value
|
||
weighted_basis = monomial_basis(n_vars, degree - 1)
|
||
w_components = _extract_sos_from_gram(Qw_val, weighted_basis, n_vars)
|
||
# For each weighted SOS component sⱼ, pair with constraint gⱼ
|
||
for s_comp in w_components:
|
||
weighted_pairs.append((s_comp, constraints_polys[idx]))
|
||
|
||
return {
|
||
"sos_components": sos_components,
|
||
"weighted_pairs": weighted_pairs,
|
||
"status": prob.status,
|
||
"solve_time": solve_time,
|
||
"solver": solver,
|
||
}
|
||
|
||
|
||
def _extract_sos_from_gram(
|
||
Q: "np.ndarray",
|
||
basis: list[tuple[int, ...]],
|
||
n_vars: int,
|
||
tol: float = 1e-8,
|
||
) -> list[SparsePoly]:
|
||
"""Extract SOS components from a Gram matrix via eigendecomposition.
|
||
|
||
Q = Σ λᵢ vᵢvᵢᵀ where λᵢ ≥ 0 (PSD).
|
||
Each eigenvector vᵢ with λᵢ > tol gives an SOS component:
|
||
qᵢ(x) = √λᵢ · Σⱼ vᵢⱼ · x^αⱼ
|
||
|
||
Returns list of SparsePoly (rationalized coefficients).
|
||
"""
|
||
eigenvalues, eigenvectors = np.linalg.eigh(Q)
|
||
components = []
|
||
|
||
for k in range(len(eigenvalues)):
|
||
lam = eigenvalues[k]
|
||
if lam < tol:
|
||
continue
|
||
sqrt_lam = math.sqrt(float(lam))
|
||
v = eigenvectors[:, k]
|
||
|
||
# Build the polynomial: qₖ(x) = √λₖ · Σⱼ vⱼ · x^αⱼ
|
||
terms: dict[tuple[int, ...], Fraction] = {}
|
||
for j, alpha in enumerate(basis):
|
||
coeff_float = sqrt_lam * float(v[j])
|
||
if abs(coeff_float) < tol:
|
||
continue
|
||
# Float → rational at the boundary (immediately rationalized)
|
||
coeff_rat = Fraction(coeff_float).limit_denominator(10**8)
|
||
key = tuple(alpha)
|
||
terms[key] = terms.get(key, Fraction(0)) + coeff_rat
|
||
|
||
if terms:
|
||
components.append(SparsePoly(
|
||
{k: v for k, v in terms.items() if v != 0}
|
||
))
|
||
|
||
return components
|
||
|
||
|
||
def _stub_solve(
|
||
problem: dict[str, Any], degree: int, level: int
|
||
) -> dict[str, Any]:
|
||
"""Stub solver for when CVXPY is not installed.
|
||
|
||
Returns a hand-constructed certificate for known test problems.
|
||
"""
|
||
name = problem["name"]
|
||
n_vars = problem["n_vars"]
|
||
|
||
if name == "test_simple":
|
||
# x₀² + x₁² = x₀² + x₁² (trivial: q₀ = x₀, q₁ = x₁)
|
||
return {
|
||
"sos_components": [
|
||
SparsePoly.from_list(2, [(Fraction(1), [1, 0])]),
|
||
SparsePoly.from_list(2, [(Fraction(1), [0, 1])]),
|
||
],
|
||
"weighted_pairs": [],
|
||
"status": "stub_optimal",
|
||
"solve_time": 0.0,
|
||
"solver": "stub",
|
||
}
|
||
elif name == "test_weighted":
|
||
# x₀² + x₀ on {x₀ ≥ 0}: sos=[x₀], weighted=[(1, x₀)]
|
||
return {
|
||
"sos_components": [
|
||
SparsePoly.from_list(1, [(Fraction(1), [1])]),
|
||
],
|
||
"weighted_pairs": [
|
||
(SparsePoly.from_list(1, [(Fraction(1), [0])]), # s₀ = 1
|
||
SparsePoly.from_list(1, [(Fraction(1), [1])])), # g₀ = x₀
|
||
],
|
||
"status": "stub_optimal",
|
||
"solve_time": 0.0,
|
||
"solver": "stub",
|
||
}
|
||
else:
|
||
return {
|
||
"sos_components": [],
|
||
"weighted_pairs": [],
|
||
"status": f"stub_unsupported:{name}",
|
||
"solve_time": 0.0,
|
||
"solver": "stub",
|
||
}
|
||
|
||
|
||
# ================================================================
|
||
# §4 RATIONALIZATION
|
||
# ================================================================
|
||
|
||
def fraction_to_lean(f: Fraction) -> str:
|
||
"""Convert a Fraction to a Lean ℚ literal string.
|
||
|
||
Examples:
|
||
Fraction(3, 1) → "3"
|
||
Fraction(1, 2) → "(1/2)"
|
||
Fraction(-5, 3) → "(-5/3)"
|
||
"""
|
||
if f.denominator == 1:
|
||
return str(f.numerator)
|
||
else:
|
||
return f"({f.numerator}/{f.denominator})"
|
||
|
||
|
||
def rationalize_float(x: float, max_denom: int = 10**8) -> Fraction:
|
||
"""Convert a float to the closest rational with bounded denominator.
|
||
|
||
This is the BOUNDARY CONVERSION: float (solver output) → ℚ (exact).
|
||
Per AGENTS.md: Float at external boundary only, immediately bracketed.
|
||
"""
|
||
return Fraction(x).limit_denominator(max_denom)
|
||
|
||
|
||
# ================================================================
|
||
# §5 LEAN CERTIFICATE EMITTER
|
||
# ================================================================
|
||
|
||
def emit_lean_certificate(
|
||
result: dict[str, Any],
|
||
problem: dict[str, Any],
|
||
degree: int,
|
||
level: int,
|
||
output_path: str,
|
||
) -> str:
|
||
"""Generate a Lean 4 data file with the SOS certificate.
|
||
|
||
Returns the SHA256 of the generated file.
|
||
"""
|
||
n_vars = problem["n_vars"]
|
||
target = problem["target"]
|
||
sos_components = result["sos_components"]
|
||
weighted_pairs = result["weighted_pairs"]
|
||
|
||
lines = [
|
||
"/-",
|
||
" Machine-generated SOS certificate.",
|
||
f" Problem: {problem['name']}",
|
||
f" Description: {problem['description']}",
|
||
f" Solver: {result['solver']}",
|
||
f" Status: {result['status']}",
|
||
f" Solve time: {result['solve_time']:.4f}s",
|
||
f" Generated: {time.strftime('%Y-%m-%dT%H:%M:%SZ', time.gmtime())}",
|
||
"",
|
||
" DO NOT EDIT — regenerate with:",
|
||
f" python3 sdp_sos_solver.py --problem {problem['name']} "
|
||
f"--degree {degree} --level {level}",
|
||
"-/",
|
||
"",
|
||
"import Semantics.SDPVerify",
|
||
"",
|
||
"open Semantics.SDPVerify",
|
||
"",
|
||
f"namespace Semantics.SDPCert.{_lean_name(problem['name'])}",
|
||
"",
|
||
]
|
||
|
||
# Emit SOS components
|
||
lines.append(f"/-- SOS certificate for: {problem['description']} -/")
|
||
lines.append(f"def certificate : SDPCertificate {n_vars} :=")
|
||
lines.append(f" mkCertificate {n_vars}")
|
||
|
||
# SOS components
|
||
lines.append(" -- SOS components (qᵢ polynomials)")
|
||
lines.append(" [" + ",\n ".join(
|
||
_poly_to_lean(q, n_vars) for q in sos_components
|
||
) + "]")
|
||
|
||
# Weighted pairs
|
||
lines.append(" -- Weighted pairs (sⱼ, gⱼ)")
|
||
if weighted_pairs:
|
||
lines.append(" [" + ",\n ".join(
|
||
f"({_poly_to_lean(s, n_vars)}, {_poly_to_lean(g, n_vars)})"
|
||
for s, g in weighted_pairs
|
||
) + "]")
|
||
else:
|
||
lines.append(" []")
|
||
|
||
# Target polynomial
|
||
lines.append(" -- Target polynomial p")
|
||
lines.append(" " + _poly_to_lean_flat(target, n_vars))
|
||
|
||
# Degree and level
|
||
lines.append(f" -- degree = {degree}, level = {level}")
|
||
lines.append(f" {degree} {level}")
|
||
lines.append("")
|
||
|
||
# Verification witness
|
||
lines.append("/-- #eval witness: certificate is valid. -/")
|
||
lines.append("#eval verifyCertificate certificate -- Expected: true")
|
||
lines.append("")
|
||
lines.append(f"end Semantics.SDPCert.{_lean_name(problem['name'])}")
|
||
lines.append("")
|
||
|
||
content = "\n".join(lines)
|
||
|
||
os.makedirs(os.path.dirname(output_path) or ".", exist_ok=True)
|
||
with open(output_path, "w") as f:
|
||
f.write(content)
|
||
|
||
sha = hashlib.sha256(content.encode()).hexdigest()
|
||
print(f"Wrote Lean certificate to {output_path} (SHA256: {sha})",
|
||
file=sys.stderr)
|
||
return sha
|
||
|
||
|
||
def _lean_name(name: str) -> str:
|
||
"""Convert a problem name to a Lean-valid namespace component."""
|
||
return "".join(c if c.isalnum() else "_" for c in name).capitalize()
|
||
|
||
|
||
def _poly_to_lean(p: SparsePoly, n_vars: int) -> str:
|
||
"""Convert a SparsePoly to a Lean list literal."""
|
||
terms = p.to_lean_list()
|
||
if not terms:
|
||
return "[]"
|
||
parts = []
|
||
for coeff_str, expon in terms:
|
||
parts.append(f"({coeff_str}, {list(expon)})")
|
||
return "[" + ", ".join(parts) + "]"
|
||
|
||
|
||
def _poly_to_lean_flat(p: SparsePoly, n_vars: int) -> str:
|
||
"""Convert a SparsePoly to a Lean list literal on a single line."""
|
||
return _poly_to_lean(p, n_vars)
|
||
|
||
|
||
# ================================================================
|
||
# §6 JSON RECEIPT EMITTER
|
||
# ================================================================
|
||
|
||
def emit_json_receipt(
|
||
result: dict[str, Any],
|
||
problem: dict[str, Any],
|
||
degree: int,
|
||
level: int,
|
||
lean_sha: str,
|
||
output_path: str,
|
||
) -> None:
|
||
"""Write a machine-readable JSON receipt for the SDP computation."""
|
||
receipt = {
|
||
"schema": "sdp_sos_certificate_v1",
|
||
"problem": {
|
||
"name": problem["name"],
|
||
"description": problem["description"],
|
||
"n_vars": problem["n_vars"],
|
||
"n_target_terms": len(problem["target"].terms),
|
||
"n_constraints": len(problem.get("constraints", [])),
|
||
},
|
||
"solver": {
|
||
"name": result["solver"],
|
||
"status": result["status"],
|
||
"solve_time_seconds": result["solve_time"],
|
||
},
|
||
"certificate": {
|
||
"n_sos_components": len(result["sos_components"]),
|
||
"n_weighted_pairs": len(result["weighted_pairs"]),
|
||
"degree": degree,
|
||
"putinar_level": level,
|
||
"lean_file_sha256": lean_sha,
|
||
},
|
||
"claim_boundary": (
|
||
"sdp-computed;python-rationalized;lean-verified"
|
||
),
|
||
"generated_utc": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
|
||
}
|
||
|
||
os.makedirs(os.path.dirname(output_path) or ".", exist_ok=True)
|
||
with open(output_path, "w") as f:
|
||
json.dump(receipt, f, indent=2)
|
||
print(f"Wrote JSON receipt to {output_path}", file=sys.stderr)
|
||
|
||
|
||
# ================================================================
|
||
# §7 CLI ENTRY POINT
|
||
# ================================================================
|
||
|
||
PROBLEMS = {
|
||
"test_simple": build_test_simple,
|
||
"test_weighted": build_test_weighted,
|
||
}
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser(
|
||
description="SDP SOS Certificate Solver",
|
||
formatter_class=argparse.RawDescriptionHelpFormatter,
|
||
epilog=__doc__,
|
||
)
|
||
parser.add_argument(
|
||
"--problem", required=True,
|
||
help="Problem name (test_simple, test_weighted, goormaghtigh_mM_nN)",
|
||
)
|
||
parser.add_argument(
|
||
"--degree", type=int, default=2,
|
||
help="Maximum degree for SOS components (default: 2)",
|
||
)
|
||
parser.add_argument(
|
||
"--level", type=int, default=0,
|
||
help="Putinar hierarchy level (default: 0 = pure SOS)",
|
||
)
|
||
parser.add_argument(
|
||
"--solver", default="SCS",
|
||
help="CVXPY solver name (default: SCS)",
|
||
)
|
||
parser.add_argument(
|
||
"--output-lean",
|
||
help="Path to write the Lean certificate file",
|
||
)
|
||
parser.add_argument(
|
||
"--output-json",
|
||
help="Path to write the JSON receipt file",
|
||
)
|
||
parser.add_argument(
|
||
"--m", type=int, default=3,
|
||
help="Goormaghtigh parameter m (default: 3)",
|
||
)
|
||
parser.add_argument(
|
||
"--n", type=int, default=3,
|
||
help="Goormaghtigh parameter n (default: 3)",
|
||
)
|
||
|
||
args = parser.parse_args()
|
||
|
||
# Build the problem
|
||
if args.problem in PROBLEMS:
|
||
problem = PROBLEMS[args.problem]()
|
||
elif args.problem.startswith("goormaghtigh"):
|
||
problem = build_goormaghtigh_fixed(args.m, args.n)
|
||
else:
|
||
print(f"ERROR: Unknown problem '{args.problem}'", file=sys.stderr)
|
||
print(f"Available: {list(PROBLEMS.keys())} + goormaghtigh",
|
||
file=sys.stderr)
|
||
sys.exit(1)
|
||
|
||
print(f"Problem: {problem['name']}", file=sys.stderr)
|
||
print(f" {problem['description']}", file=sys.stderr)
|
||
print(f" {len(problem['target'].terms)} target terms", file=sys.stderr)
|
||
print(f" {len(problem.get('constraints', []))} constraints",
|
||
file=sys.stderr)
|
||
print(f" Degree: {args.degree}, Level: {args.level}", file=sys.stderr)
|
||
|
||
# Solve
|
||
result = solve_sos_sdp(
|
||
problem, args.degree, args.level, args.solver
|
||
)
|
||
|
||
print(f"Status: {result['status']}", file=sys.stderr)
|
||
print(f"Solve time: {result['solve_time']:.4f}s", file=sys.stderr)
|
||
print(f"SOS components: {len(result['sos_components'])}", file=sys.stderr)
|
||
print(f"Weighted pairs: {len(result['weighted_pairs'])}", file=sys.stderr)
|
||
|
||
# Verify locally (Python-side sanity check)
|
||
if result["sos_components"]:
|
||
recon = SparsePoly()
|
||
for q in result["sos_components"]:
|
||
recon = recon + q.square()
|
||
for s, g in result["weighted_pairs"]:
|
||
recon = recon + (s * g)
|
||
|
||
# Compare with target
|
||
diff = recon + problem["target"].scale(Fraction(-1))
|
||
max_residual = max(
|
||
(abs(v) for v in diff.terms.values()), default=Fraction(0)
|
||
)
|
||
print(f"Max residual: {float(max_residual):.2e}", file=sys.stderr)
|
||
|
||
# Emit outputs
|
||
lean_sha = ""
|
||
if args.output_lean:
|
||
lean_sha = emit_lean_certificate(
|
||
result, problem, args.degree, args.level, args.output_lean
|
||
)
|
||
if args.output_json:
|
||
emit_json_receipt(
|
||
result, problem, args.degree, args.level,
|
||
lean_sha, args.output_json
|
||
)
|
||
|
||
# Summary
|
||
if result["status"] in ("optimal", "optimal_inaccurate", "stub_optimal"):
|
||
print("SUCCESS: SOS certificate computed.", file=sys.stderr)
|
||
sys.exit(0)
|
||
else:
|
||
print(f"FAILED: solver status = {result['status']}", file=sys.stderr)
|
||
sys.exit(1)
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|