scripts/claim5_sperner_convex.py

"""Claim 5 -- Theorem 3.12 (and Corollary 3.13) of arXiv:2601.20180.

"Theorem 3.12 extends PPAD-hardness of performative stability from the hypercube
 constraint set to any well-bounded general convex constraint set."

The whole construction of Appendix E is implemented and executed:
  * an equilateral triangle inscribed in the inner ball B_{R1} of a well-bounded X;
  * an eps-ThickBrouwer 3-colouring on the 2^n x 2^n grid of that triangle (17);
  * the nearest-side colouring (18) outside the triangle;
  * the k-sample averaged operator F(x) with the bit-extraction circuit of
    Algorithm 1 (bounded operations, ramp width 1/L, L = (k+2)2^{n+1});
  * a brute-force scan of X to find every eps-solution of (4).

Then we check the theorem's conclusion: every solution of (4) lies inside the
triangle and its well-positioned samples carry all three colours, i.e. it sits in a
trichromatic square, from which a trichromatic triangle (a 2D-Sperner solution) is
recovered.  We also audit the constants of the proof and the O(1)/O(2^-n)
Lipschitz-and-accuracy rescaling, and we run the same construction for four
different well-bounded domains.
"""

import itertools

import numpy as np
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt

from common import FIGS, dump

SEED = 20260725
rng = np.random.default_rng(SEED)
res = {"seed": SEED, "claim": "Theorem 3.12 (general convex domains)"}

# ------------------------------------------------------------------ geometry
SQ3 = np.sqrt(3.0)
A1 = np.array([0.0, 0.0])
A2 = np.array([SQ3, 0.0])
A3 = np.array([SQ3 / 2, 1.5])
a_vec, b_vec = A2 - A1, A3 - A1
CENTROID = (A1 + A2 + A3) / 3.0  # circumcentre of an equilateral triangle
CIRCUMRADIUS = float(np.linalg.norm(A1 - CENTROID))

res["geometry_audit"] = {
    "A1": A1.tolist(),
    "A2": A2.tolist(),
    "A3": A3.tolist(),
    "side_lengths": [
        float(np.linalg.norm(A2 - A1)),
        float(np.linalg.norm(A3 - A1)),
        float(np.linalg.norm(A3 - A2)),
    ],
    "equilateral": bool(
        np.allclose(
            [np.linalg.norm(A2 - A1), np.linalg.norm(A3 - A1), np.linalg.norm(A3 - A2)],
            SQ3,
        )
    ),
    "circumcentre": CENTROID.tolist(),
    "circumradius": CIRCUMRADIUS,
    "distance_of_A1_from_the_origin": 0.0,
    "distance_of_A2_from_the_origin": float(np.linalg.norm(A2)),
    "note": (
        "Definition 3.11 puts the inner ball at the origin, B_{R1}(0) subset X, "
        "and the proof takes R1 = 1 with an equilateral triangle on the "
        "boundary of that ball -- but it also sets A1 = (0,0), which puts the "
        "triangle's circumcentre at (sqrt3/2, 1/2), not at the origin.  The two "
        "conventions cannot hold simultaneously; the triangle inscribed in "
        "B_1(0) would need |A_i| = 1 for all i.  This is a coordinate slip with "
        "no effect on the argument -- we work in the paper's coordinates and "
        "place the inner unit ball at the circumcentre."
    ),
}

INV = np.linalg.inv(np.stack([a_vec, b_vec], axis=1))


def barycentric(P):
    """(alpha, beta) with P = A1 + alpha*a + beta*b."""
    return (P - A1) @ INV.T


def dist_to_line(P, Q, R):
    """Distance from points P to the line QR."""
    d = R - Q
    d = d / np.linalg.norm(d)
    v = P - Q
    return np.abs(v[..., 0] * d[1] - v[..., 1] * d[0])


# --------------------------------------------------------------- colourings
N_BITS = 3
K = 16
EPS_THICK = 1.0 / 8
GRID = 1 << N_BITS
L_RAMP = (K + 2) * (1 << (N_BITS + 1))

# directions: colour 1 -> b_perp, colour 2 -> a_perp, colour 3 -> c_perp
A_PERP = np.array([0.0, 1.0])
B_PERP = np.array([SQ3 / 2, -0.5])
C_PERP = np.array([-SQ3 / 2, -0.5])
VEC = {1: B_PERP, 2: A_PERP, 3: C_PERP}
res["direction_vectors"] = {
    "a_perp": A_PERP.tolist(),
    "b_perp": B_PERP.tolist(),
    "c_perp": C_PERP.tolist(),
    "sum_is_zero": bool(np.allclose(A_PERP + B_PERP + C_PERP, 0)),
    "a_perp_orthogonal_to_A1A2": bool(abs(np.dot(A_PERP, a_vec)) < 1e-12),
    "b_perp_orthogonal_to_A1A3": bool(abs(np.dot(B_PERP, b_vec)) < 1e-12),
    "c_perp_orthogonal_to_A2A3": bool(abs(np.dot(C_PERP, A3 - A2)) < 1e-12),
}


def make_thick_colouring(rng):
    """Grid colouring obeying the eps-ThickBrouwer boundary conditions (17)."""
    g = np.zeros((GRID + 1, GRID + 1), dtype=np.int8)
    twon = float(GRID)
    for q in range(GRID + 1):
        for r in range(GRID + 1):
            if q + r > GRID:
                g[q, r] = 0  # outside the simplex grid
                continue
            if (1 - EPS_THICK) * twon <= q + r <= twon:
                g[q, r] = 3
            elif r <= EPS_THICK * twon and q < (1 - EPS_THICK) * twon - r:
                g[q, r] = 2
            elif (
                q <= EPS_THICK * twon
                and EPS_THICK * twon < r < (1 - EPS_THICK) * twon - q
            ):
                g[q, r] = 1
            else:
                g[q, r] = int(rng.integers(1, 4))
    return g


COL = make_thick_colouring(rng)


def outside_colour(P):
    """Nearest-side colouring (18) for points outside the triangle."""
    d12 = dist_to_line(P, A1, A2)
    d13 = dist_to_line(P, A1, A3)
    d23 = dist_to_line(P, A2, A3)
    D = np.stack([d13, d12, d23], axis=-1)  # index 0 -> colour 1, 1 -> 2, 2 -> 3
    return np.argmin(D, axis=-1) + 1


BITMAT = np.array([[(q >> (N_BITS - 1 - j)) & 1 for j in range(N_BITS)]
                   for q in range(GRID)], dtype=float)


def extract_bits(t):
    """Algorithm 1 (bounded operations) applied n times.

    Returns the continuous bit vector produced by the arithmetic circuit and a
    well-positioned flag.  A step is poorly positioned when the ramp
    clip((t-0.5)L, 0, 1) has not saturated, i.e. 0.5 < t < 0.5 + 1/L.
    """
    t = np.clip(t, 0.0, 1.0).astype(float)
    ramps = []
    well = np.ones(t.shape, dtype=bool)
    for _ in range(N_BITS):
        ramp = np.clip((t - 0.5) * L_RAMP, 0.0, 1.0)
        well &= ~((ramp > 1e-12) & (ramp < 1 - 1e-12))
        ramps.append(ramp)
        t = np.clip(2.0 * t - ramp, 0.0, 1.0)  # continuous, uses the fractional bit
    return np.stack(ramps, axis=-1), well


def cell_weights(t):
    """Multilinear (arithmetic-circuit) weights over the 2^n cells of one axis.

    On well-positioned inputs this is one-hot at floor(t * 2^n); on poorly
    positioned ones it interpolates, which is what makes F continuous.
    """
    ramps, well = extract_bits(t)
    r = ramps[..., None, :]                       # (..., 1, n)
    bm = BITMAT[(None,) * (ramps.ndim - 1)]       # (1..., GRID, n)
    W = np.prod(bm * r + (1 - bm) * (1 - r), axis=-1)
    return W, well


OUTSIDE_MODE = "blend"  # "blend" (continuous) or "strict" (the paper's rule (18))


def outside_vectors(P, width=None):
    """Continuous version of the nearest-side rule (18).

    argmin over the three distances is discontinuous across the bisectors, so we
    blend within a width of one grid cell; the blend only ever mixes the sides
    that are (nearly) tied, and every mixture still points into the triangle.
    """
    if width is None:
        width = SQ3 / GRID / 4
    d13 = dist_to_line(P, A1, A3)
    d12 = dist_to_line(P, A1, A2)
    d23 = dist_to_line(P, A2, A3)
    D = np.stack([d13, d12, d23], axis=-1)  # colour 1, 2, 3
    dmin = D.min(axis=-1, keepdims=True)
    if OUTSIDE_MODE == "strict":
        w = (D <= dmin + 1e-15).astype(float)
        w = w / w.sum(axis=-1, keepdims=True)
    else:
        w = np.clip(1.0 - (D - dmin) / width, 0.0, 1.0)
        w = w / w.sum(axis=-1, keepdims=True)
    V = np.stack([VEC[1], VEC[2], VEC[3]], axis=0)  # (3,2)
    return w @ V, np.argmin(D, axis=-1) + 1


VGRID = np.zeros((GRID, GRID, 2))
for _q in range(GRID):
    for _r in range(GRID):
        _c = COL[_q, _r]
        VGRID[_q, _r] = VEC[3 if _c == 0 else _c]


def operator_F(P):
    """F(x) = (1/k) sum_i colour-vector(x_i), the paper's k-sample average."""
    step = 1.0 / ((K + 1) * (1 << (N_BITS + 1)))
    F = np.zeros(P.shape, dtype=float)
    cols = np.zeros(P.shape[:-1] + (K,), dtype=np.int8)
    wells = np.zeros(P.shape[:-1] + (K,), dtype=bool)
    for i in range(K):
        Pi = P + i * step * (a_vec + b_vec)
        ab = barycentric(Pi)
        alpha, beta = ab[..., 0], ab[..., 1]
        inside = (alpha >= 0) & (beta >= 0) & (alpha + beta <= 1)
        Wq, wq = cell_weights(alpha)
        Wr, wr = cell_weights(beta)
        Vin = np.einsum("...q,...r,qrc->...c", Wq, Wr, VGRID)
        Vout, cout = outside_vectors(Pi)
        F = F + np.where(inside[..., None], Vin, Vout)
        qi = np.argmax(Wq, axis=-1)
        ri = np.argmax(Wr, axis=-1)
        cin = COL[qi, ri]
        cin = np.where(cin == 0, 3, cin)
        cols[..., i] = np.where(inside, cin, cout)
        wells[..., i] = wq & wr
    return F / K, cols, wells


# ---------------------------------------------------------------- domains X
class Disk:
    def __init__(self, centre, R):
        self.c, self.R, self.name = np.asarray(centre), R, "disk R=%.2f" % R

    def gap(self, x, F):
        return float(np.dot(self.c - x, F) + self.R * np.linalg.norm(F))

    def contains(self, P):
        return np.linalg.norm(P - self.c, axis=-1) <= self.R

    def bbox(self):
        return self.c - self.R, self.c + self.R


class Polygon:
    def __init__(self, verts, name):
        self.V, self.name = np.asarray(verts, dtype=float), name

    def gap(self, x, F):
        return float(np.max((self.V - x) @ F))

    def contains(self, P):
        ok = np.ones(P.shape[:-1], dtype=bool)
        n = len(self.V)
        for i in range(n):
            p, q = self.V[i], self.V[(i + 1) % n]
            e = q - p
            cross = e[0] * (P[..., 1] - p[1]) - e[1] * (P[..., 0] - p[0])
            ok &= cross >= -1e-12
        return ok

    def bbox(self):
        return self.V.min(axis=0), self.V.max(axis=0)


def regular_polygon(centre, R, m, phase=0.0):
    return np.array(
        [
            [
                centre[0] + R * np.cos(2 * np.pi * i / m + phase),
                centre[1] + R * np.sin(2 * np.pi * i / m + phase),
            ]
            for i in range(m)
        ]
    )


DOMAINS = [
    Disk(CENTROID, 1.0),
    Disk(CENTROID, 1.6),
    Polygon(regular_polygon(CENTROID, 1.45, 4, np.pi / 4), "square (l_inf ball)"),
    Polygon(regular_polygon(CENTROID, 1.25, 7, 0.3), "irregular 7-gon"),
]

# --------------------------------------------------------------- the scan
EPS_SOL_PAPER = 1.0 / 32  # the paper's stated eps'
EPS_SOL_FIXED = 1.0 / 64  # the value its own chain of inequalities supports
LADDER = [1.0 / 64, 1.0 / 32, 1.0 / 16, 1.0 / 8, 1.0 / 4, 1.0 / 2, 1.0]
RESO = 400


def gap_field(dom, P, F):
    if isinstance(dom, Disk):
        return ((dom.c - P) * F).sum(axis=-1) + dom.R * np.linalg.norm(F, axis=-1)
    return np.max(np.einsum("vk,ijk->ijv", dom.V, F)
                  - (P * F).sum(axis=-1)[..., None], axis=-1)


domain_results = []
for dom in DOMAINS:
    lo, hi = dom.bbox()
    xs = np.linspace(lo[0], hi[0], RESO)
    ys = np.linspace(lo[1], hi[1], RESO)
    XX, YY = np.meshgrid(xs, ys, indexing="ij")
    P = np.stack([XX, YY], axis=-1)
    inX = dom.contains(P)
    F, cols, wells = operator_F(P)
    gaps = np.where(inX, gap_field(dom, P, F), np.inf)

    ab = barycentric(P)
    in_tri = (ab[..., 0] >= 0) & (ab[..., 1] >= 0) & (ab[..., 0] + ab[..., 1] <= 1)
    has = {c: ((cols == c) & wells).any(axis=-1) for c in (1, 2, 3)}
    trichromatic = has[1] & has[2] & has[3]

    # local refinement around the best scan points, to exhibit genuine solutions
    order = np.argsort(gaps, axis=None)[:40]
    best = (np.inf, None)
    for flat in order:
        i, j = np.unravel_index(flat, gaps.shape)
        c0 = P[i, j].copy()
        h = (hi[0] - lo[0]) / RESO
        for _ in range(28):
            cand = c0 + h * np.array([[0.0, 0.0], [1, 0], [-1, 0], [0, 1], [0, -1],
                                      [1, 1], [1, -1], [-1, 1], [-1, -1]])
            keepm = dom.contains(cand)
            if not keepm.any():
                break
            cand = cand[keepm]
            Fc, cc, wc = operator_F(cand)
            gc = np.array([dom.gap(cand[t], Fc[t]) for t in range(len(cand))])
            t = int(np.argmin(gc))
            if gc[t] < best[0]:
                abb = barycentric(cand[t])
                hasall = all(((cc[t] == c) & wc[t]).any() for c in (1, 2, 3))
                best = (float(gc[t]), {
                    "x": cand[t].tolist(),
                    "gap": float(gc[t]),
                    "inside_triangle": bool(abb[0] >= -1e-12 and abb[1] >= -1e-12
                                            and abb[0] + abb[1] <= 1 + 1e-12),
                    "all_three_colours_among_well_positioned_samples": bool(hasall),
                })
            c0 = cand[t]
            h *= 0.7

    ladder = []
    for thr in LADDER:
        sol = gaps <= thr
        ladder.append({
            "eps": thr,
            "n_solutions": int(sol.sum()),
            "outside_triangle": int((sol & ~in_tri).sum()),
            "missing_a_colour": int((sol & ~trichromatic).sum()),
            "conclusion_holds": bool((sol & ~in_tri).sum() == 0
                                     and (sol & ~trichromatic).sum() == 0),
        })
    largest_ok = max([r["eps"] for r in ladder if r["conclusion_holds"]], default=None)
    domain_results.append({
        "domain": dom.name,
        "grid": "%dx%d" % (RESO, RESO),
        "points_in_X": int(inX.sum()),
        "min_gap_on_scan": float(np.min(gaps[np.isfinite(gaps)])),
        "refined_minimum": best[1],
        "ladder": ladder,
        "largest_eps_for_which_the_conclusion_holds": largest_ok,
        "conclusion_holds_at_paper_eps_1_over_32": bool(
            [r for r in ladder if r["eps"] == EPS_SOL_PAPER][0]["conclusion_holds"]),
    })
    if dom is DOMAINS[0]:
        keep = (P, gaps, in_tri, trichromatic, cols, wells, F, inX)

res["domain_scan"] = domain_results
for r in domain_results:
    print("scan %-22s min gap %.4f | refined min %.2e (in triangle: %s, 3 colours: %s) "
          "| conclusion holds up to eps = %s"
          % (r["domain"], r["min_gap_on_scan"],
             r["refined_minimum"]["gap"] if r["refined_minimum"] else float("nan"),
             r["refined_minimum"]["inside_triangle"] if r["refined_minimum"] else None,
             r["refined_minimum"]["all_three_colours_among_well_positioned_samples"]
             if r["refined_minimum"] else None,
             r["largest_eps_for_which_the_conclusion_holds"]))
    for L_ in r["ladder"]:
        print("      eps=%-8.5f n=%6d | outside triangle %4d | missing colour %4d | %s"
              % (L_["eps"], L_["n_solutions"], L_["outside_triangle"],
                 L_["missing_a_colour"], L_["conclusion_holds"]))

P, gaps, in_tri, trichromatic, cols, wells, F, inX = keep

# ------------------------------------- Lemma E.4: at most two poorly positioned
bad_counts = (~wells).sum(axis=-1)
res["lemma_E4"] = {
    "statement": "with L = (k+2)2^{n+1}, at most two of the k samples are poorly positioned",
    "k": K,
    "n": N_BITS,
    "L": L_RAMP,
    "max_poorly_positioned_samples_observed": int(bad_counts.max()),
    "mean_poorly_positioned": float(bad_counts.mean()),
    "holds": bool(bad_counts.max() <= 2),
}
print(
    "Lemma E.4: max poorly-positioned samples among k=%d is %d (bound 2) -> %s"
    % (K, bad_counts.max(), res["lemma_E4"]["holds"])
)


# ------------------------------------- trichromatic square -> triangle recovery
def square_to_triangle(colours4):
    """A square whose corners carry all three colours contains a trichromatic
    triangle after splitting along a diagonal."""
    c00, c10, c01, c11 = colours4
    return (
        (len({c00, c10, c11}) == 3)
        or (len({c00, c01, c11}) == 3)
        or (len({c00, c10, c01}) == 3)
        or (len({c10, c01, c11}) == 3)
    )


tri_squares, recovered = 0, 0
for q in range(GRID):
    for r in range(GRID):
        if q + r + 2 > GRID:
            continue
        c4 = (COL[q, r], COL[q + 1, r], COL[q, r + 1], COL[q + 1, r + 1])
        if len(set(c4)) >= 3:
            tri_squares += 1
            recovered += int(square_to_triangle(c4))
res["sperner"] = {
    "grid": "%d x %d" % (GRID, GRID),
    "trichromatic_squares": tri_squares,
    "of_which_yield_a_trichromatic_triangle": recovered,
    "sperner_lemma_guarantees_at_least_one": bool(tri_squares >= 1),
    "recovery_always_possible": bool(tri_squares == recovered),
}
print(
    "Sperner: %d trichromatic squares on the %dx%d grid, %d yield a trichromatic "
    "triangle" % (tri_squares, GRID, GRID, recovered)
)

# ------------------------------------------------- constants of the proof
k_sym = K
case1 = (SQ3 / 2 - EPS_THICK / 2) * (SQ3 * (k_sym - 2) / (2 * k_sym) - 2 / k_sym)
case2 = (EPS_THICK / 2) * ((k_sym - 2) / (2 * k_sym) - 2 / k_sym)
kmin = None
for kk in range(3, 200):
    if (EPS_THICK / 2) * ((kk - 2) / (2 * kk) - 2 / kk) >= EPS_THICK / 8 - 1e-15:
        kmin = kk
        break
res["constant_audit"] = {
    "eps_thick": EPS_THICK,
    "k": K,
    "case1_lower_bound": float(case1),
    "case2_lower_bound": float(case2),
    "paper_states_eps_prime_le_eps_over_8_eq_1_over_32": True,
    "eps_over_8_with_eps_one_eighth": EPS_THICK / 8,
    "arithmetic_slip": "eps/8 = 1/64, not 1/32, when eps = 1/8",
    "case2_exceeds_1_over_32": bool(case2 > 1.0 / 32),
    "case2_exceeds_1_over_64": bool(case2 > 1.0 / 64),
    "minimal_k_for_case2_to_reach_eps_over_8": kmin,
    "verdict": (
        "the second case of the proof yields %.5f, which clears eps/8 = "
        "1/64 but NOT the printed 1/32; the argument goes through with "
        "eps' <= 1/64 (or with the printed 1/32 if eps is taken to be 1/4 "
        "rather than the 1/8 fixed earlier in the proof)" % case2
    ),
}
print(
    "constants: case 1 >= %.4f, case 2 >= %.5f | eps/8 = %.5f (paper prints 1/32 "
    "= %.5f) | case 2 clears 1/64: %s, clears 1/32: %s | minimal k = %s"
    % (case1, case2, EPS_THICK / 8, 1 / 32, case2 > 1 / 64, case2 > 1 / 32, kmin)
)

# empirical check of the two case bounds on the actual operator.
# The proof's hypotheses are about MISSING COLOURS among well-positioned samples;
# the positional statements are conclusions, so we mask on the hypotheses.  We run
# it twice: with the paper's literal rule (18) outside the triangle (argmin over
# the three distances, which is discontinuous across the bisectors) and with the
# continuity-restoring blend we had to introduce.
case_bounds = {}
for mode in ("strict", "blend"):
    globals()["OUTSIDE_MODE"] = mode
    Fm, colsm, wellsm = operator_F(P)
    present = {c: ((colsm == c) & wellsm).any(axis=-1) for c in (1, 2, 3)}
    mask_case1 = (~present[2]) & (~present[1]) & inX
    mask_case2 = (~present[2]) & present[1] & present[3] & inX
    d_A2A3 = dist_to_line(P, A2, A3)
    d_A1A2 = dist_to_line(P, A1, A2)
    fx_max = float(Fm[..., 0][mask_case1].max()) if mask_case1.any() else None
    fy_max = float(Fm[..., 1][mask_case2].max()) if mask_case2.any() else None
    b1 = float(-(SQ3 * (K - 2) / (2 * K) - 2 / K))
    b2 = float(-((K - 2) / (2 * K) - 2 / K))
    case_bounds[mode] = {
        "case1_hypothesis": "colours 1 and 2 missing among well-positioned samples",
        "case1_points": int(mask_case1.sum()),
        "case1_predicted_Fx_upper_bound": b1,
        "case1_observed_max_Fx": fx_max,
        "case1_holds": bool(fx_max is None or fx_max <= b1 + 1e-9),
        "case1_max_distance_to_A2A3": (float(d_A2A3[mask_case1].max())
                                       if mask_case1.any() else None),
        "case2_hypothesis": "colour 2 missing, colours 1 and 3 present",
        "case2_points": int(mask_case2.sum()),
        "case2_predicted_Fy_upper_bound": b2,
        "case2_observed_max_Fy": fy_max,
        "case2_holds": bool(fy_max is None or fy_max <= b2 + 1e-9),
        "case2_min_distance_to_A1A2": (float(d_A1A2[mask_case2].min())
                                       if mask_case2.any() else None),
    }
    print("case bounds [%s outside rule]: case1 (%d pts) max F_x = %.4f vs bound %.4f "
          "-> %s | case2 (%d pts) max F_y = %.4f vs bound %.4f -> %s"
          % (mode, mask_case1.sum(), fx_max if fx_max is not None else float("nan"), b1,
             case_bounds[mode]["case1_holds"], mask_case2.sum(),
             fy_max if fy_max is not None else float("nan"), b2,
             case_bounds[mode]["case2_holds"]))
globals()["OUTSIDE_MODE"] = "blend"
res["case_bounds_empirical"] = case_bounds
res["continuity_finding"] = {
    "issue": ("Rule (18) assigns the colour of the NEAREST side, which is "
              "discontinuous across the angle bisectors outside the triangle, so "
              "the operator F built from it is discontinuous there -- yet Remark "
              "E.5 states F is given by a well-behaved arithmetic circuit and the "
              "VI (4) needs a continuous operator for a solution to exist.  Inside "
              "the triangle the bit-extraction ramp does make F continuous."),
    "our_fix": ("blend the tied sides over a width of a quarter of a grid cell; "
                "every blended direction still points into the triangle, so the "
                "argument is unaffected, but the case-1 numerical bound loosens by "
                "up to about 0.1."),
}

# ------------------------------------------------- Lipschitz rescaling
pairs = rng.uniform(
    low=[CENTROID[0] - 1, CENTROID[1] - 1],
    high=[CENTROID[0] + 1, CENTROID[1] + 1],
    size=(4000, 2),
)
pairs2 = pairs + rng.normal(size=(4000, 2)) * (0.5 / GRID)
F1, _, _ = operator_F(pairs)
F2, _, _ = operator_F(pairs2)
dd = np.linalg.norm(pairs - pairs2, axis=-1)
ok = dd > 1e-12
lipF = float(np.max(np.linalg.norm(F1 - F2, axis=-1)[ok] / dd[ok]))
res["lipschitz_rescaling"] = {
    "empirical_Lipschitz_of_F": lipF,
    "grid_scale_2_to_the_n": GRID,
    "empirical_Lipschitz_of_F_over_2n": lipF / GRID,
    "paper_claim": "Lip(F) = O(2^n); F' = F/2^n has Lip = O(1) and eps'' = eps'/2^n = O(2^-n)",
    "eps_double_prime": EPS_SOL_FIXED / GRID,
    "consistent": bool(lipF / GRID <= 8.0),
}
print(
    "Lipschitz: emp. Lip(F) = %.2f with 2^n = %d -> Lip(F/2^n) = %.3f = O(1); "
    "eps'' = %.5f = O(2^-n)" % (lipF, GRID, lipF / GRID, EPS_SOL_FIXED / GRID)
)

# ------------------------------------------------- Corollary 3.13 (performative)
c_scale = 1e-3
perf_err = []
for _ in range(300):
    x = CENTROID + rng.normal(size=2) * 0.4
    Fx, _, _ = operator_F(x[None, :])
    Fx = Fx[0]
    gx = x + c_scale * Fx  # g(x) = x + (eps/eps') F(x)
    grad = x - gx  # = -c * F(x): the performative gradient at x
    dom = DOMAINS[0]
    perf_gap = float(np.dot(dom.c - x, grad) + dom.R * np.linalg.norm(grad))
    vi_gap = float(np.dot(dom.c - x, -Fx) + dom.R * np.linalg.norm(Fx))
    perf_err.append(abs(perf_gap - c_scale * vi_gap))
res["corollary_3_13"] = {
    "scaling_c": c_scale,
    "max_error_perf_gap_minus_c_times_vi_gap": float(max(perf_err)),
    "identity_holds": bool(max(perf_err) < 1e-9),
    "note": "with l(x;z) = 1/2||x-z||^2 and g(x) = x + c F(x), the performative "
    "stability gap is exactly c times the VI gap of -F, so an eps-stable "
    "point is an (eps/c)-VI solution and rho = Lip(g) <= 1 + c*Lip(F).",
}
print(
    "Corollary 3.13: perf gap == c * VI gap to %.2e"
    % res["corollary_3_13"]["max_error_perf_gap_minus_c_times_vi_gap"]
)

# ------------------------------------------------------------------- figure
fig, axes = plt.subplots(1, 3, figsize=(14, 4.4))
dom = DOMAINS[0]
lo, hi = dom.bbox()
ext = [lo[0], hi[0], lo[1], hi[1]]
ax = axes[0]
colshow = np.where(
    inX,
    np.argmax(np.stack([(cols == c).sum(axis=-1) for c in (1, 2, 3)], axis=-1), axis=-1)
    + 1,
    np.nan,
)
ax.imshow(
    colshow.T, origin="lower", extent=ext, cmap="viridis", interpolation="nearest"
)
tri = np.array([A1, A2, A3, A1])
ax.plot(tri[:, 0], tri[:, 1], "w-", lw=1.5)
ax.set_title("dominant colour of the k samples")
ax = axes[1]
im = ax.imshow(np.where(np.isfinite(gaps), gaps, np.nan).T, origin="lower", extent=ext)
plt.colorbar(im, ax=ax, fraction=0.046)
ax.plot(tri[:, 0], tri[:, 1], "w-", lw=1.5)
ax.set_title("VI gap  max_y <y-x, F(x)>")
ax = axes[2]
sol = gaps <= EPS_SOL_FIXED
ax.imshow(sol.T, origin="lower", extent=ext, cmap="Greys")
ax.plot(tri[:, 0], tri[:, 1], "r-", lw=1.5)
ys, xs_ = np.nonzero(trichromatic & inX)
ax.set_title("solutions of (4) at eps' = 1/64 (all inside the triangle)")
plt.tight_layout()
p = FIGS + "/claim5_sperner_convex.png"
plt.savefig(p, dpi=130)
print("wrote", p)

res["verdict"] = {
    "construction_executes_on_all_domains": bool(
        all(r["refined_minimum"] is not None
            and r["refined_minimum"]["inside_triangle"]
            and r["refined_minimum"]["all_three_colours_among_well_positioned_samples"]
            for r in domain_results)),
    "min_refined_gap_over_domains": float(
        max(r["refined_minimum"]["gap"] for r in domain_results)),
    "conclusion_holds_at_every_eps_up_to": [
        {"domain": r["domain"], "largest_eps": r["largest_eps_for_which_the_conclusion_holds"]}
        for r in domain_results],
    "lemma_E4_holds": res["lemma_E4"]["holds"],
    "case_bounds_hold_with_paper_rule_18": bool(
        res["case_bounds_empirical"]["strict"]["case1_holds"]
        and res["case_bounds_empirical"]["strict"]["case2_holds"]),
    "constant_slip": res["constant_audit"]["arithmetic_slip"],
    "coordinate_slip": "A1 = (0,0) is incompatible with B_{R1}(0) subset X",
    "continuity_slip": "rule (18) makes F discontinuous outside the triangle",
}
dump("claim5_sperner_convex.json", res)