1"""Anytime primal-dual LP hierarchy MVP.
 2
 3Solves max c^T v subject to A v <= b and its nonnegative dual.
 4Nested primal bases S and nested dual columns T produce certified lower
 5and upper bounds whenever the corresponding reduced LPs are feasible.
 6"""
 7from dataclasses import dataclass
 8import numpy as np
 9from scipy.optimize import linprog
10
11
12@dataclass
13class Bound:
14    level: int
15    lower: float
16    upper: float
17    width: float
18    primal_residual: float
19    dual_eq_residual: float
20    dual_nonnegative: bool
21    covers_full: bool
22
23
24def solve_full(A, b, c):
25    A, b, c = map(lambda x: np.asarray(x, dtype=float), (A, b, c))
26    p = linprog(-c, A_ub=A, b_ub=b,
27                bounds=[(None, None)] * len(c), method="highs")
28    d = linprog(b, A_eq=A.T, b_eq=c,
29                bounds=[(0, None)] * len(b), method="highs")
30    if not p.success or not d.success:
31        raise RuntimeError(f"full LP failed: {p.message}; {d.message}")
32    return float(c @ p.x), float(b @ d.x), p.x, d.x
33
34
35def solve_reduced_primal(A, b, c, S):
36    """Return lower bound and full-space primal vector."""
37    S = np.asarray(S, dtype=float)
38    q = S.T @ c
39    sol = linprog(-q, A_ub=A @ S, b_ub=b,
40                  bounds=[(None, None)] * S.shape[1], method="highs")
41    if not sol.success:
42        raise RuntimeError(f"reduced primal failed: {sol.message}")
43    v = S @ sol.x
44    return float(c @ v), v, float(np.max(A @ v - b))
45
46
47def solve_reduced_dual(A, b, c, T):
48    """Return upper bound and full-space dual vector y=T alpha."""
49    T = np.asarray(T, dtype=float)
50    # alpha >= 0; equality is A.T T alpha = c.
51    sol = linprog(b @ T, A_eq=A.T @ T, b_eq=c,
52                  bounds=[(0, None)] * T.shape[1], method="highs")
53    if not sol.success:
54        raise RuntimeError(f"reduced dual failed: {sol.message}")
55    y = T @ sol.x
56    return float(b @ y), y, float(np.max(np.abs(A.T @ y - c)))
57
58
59def hierarchy(A, b, c, primal_bases, dual_bases, full_opt=None):
60    """Evaluate paired nested spaces and assert certificate properties."""
61    if full_opt is None:
62        full_opt = solve_full(A, b, c)[0]
63    out = []
64    for j, (S, T) in enumerate(zip(primal_bases, dual_bases)):
65        L, v, pres = solve_reduced_primal(A, b, c, S)
66        U, y, deres = solve_reduced_dual(A, b, c, T)
67        out.append(Bound(j, L, U, U - L, pres, deres,
68                         bool(np.min(y) >= -1e-8),
69                         bool(L - 1e-7 <= full_opt <= U + 1e-7)))
70    return out
71
72
73def stop_level(bounds, tau=0.02):
74    for x in bounds:
75        if x.width <= tau * max(1.0, abs(x.upper)):
76            return x.level
77    return None