#!/usr/bin/env python3
"""Exact examples for Type II core unistationarity (paper dated 7 September 2026).

Run: python type_ii_demo.py
Optional plotting dependencies: numpy, matplotlib. For arithmetic only, use
    python type_ii_demo.py --no-figures

All certificates use fractions.Fraction, without a numerical tolerance.
The general uniqueness theorem is supplied by the paper, not proved by sampling.
The phase map is a derived consequence for the symmetric family defined here.
"""

import argparse
import csv
from dataclasses import dataclass
from fractions import Fraction as Q
import hashlib
import itertools
import json
import math
from pathlib import Path
import platform


def require(condition, message):
    """Keep checks active even if Python is run with -O."""
    if not condition:
        raise AssertionError(message)


def matvec(a, x):
    return [sum((v * w for v, w in zip(row, x)), Q(0)) for row in a]


def matmul(a, b):
    return [[sum((a[i][k] * b[k][j] for k in range(len(b))), Q(0))
             for j in range(len(b[0]))] for i in range(len(a))]


def solve(a, b):
    """Gaussian elimination with exact fractions; singular input raises an error."""
    n = len(a)
    m = [[Q(v) for v in row] + [Q(b[i])] for i, row in enumerate(a)]
    for c in range(n):
        pivot = next((r for r in range(c, n) if m[r][c]), None)
        require(pivot is not None, "Singular exact linear system")
        m[c], m[pivot] = m[pivot], m[c]
        p = m[c][c]
        m[c] = [v / p for v in m[c]]
        for r in range(n):
            if r != c:
                p = m[r][c]
                m[r] = [v - p * w for v, w in zip(m[r], m[c])]
    return [row[-1] for row in m]


def determinant(a):
    m = [[Q(v) for v in row] for row in a]
    result = Q(1)
    for c in range(len(m)):
        pivot = next((r for r in range(c, len(m)) if m[r][c]), None)
        if pivot is None:
            return Q(0)
        if pivot != c:
            m[c], m[pivot] = m[pivot], m[c]
            result = -result
        p = m[c][c]
        result *= p
        for r in range(c + 1, len(m)):
            t = m[r][c] / p
            for j in range(c + 1, len(m)):
                m[r][j] -= t * m[c][j]
    return result


def inverse(a):
    columns = [solve(a, [int(i == j) for i in range(len(a))])
               for j in range(len(a))]
    return [list(row) for row in zip(*columns)]


@dataclass
class Network:
    name: str
    weights: tuple
    backs: dict
    forward: list
    reverse: list
    degradation: list

    def __post_init__(self):
        self.n = len(self.weights)
        require(self.n >= 3, "At least three species required")
        require(all(isinstance(w, int) and w >= 1 for w in self.weights),
                "Product coefficients must be positive integers")
        require(all(0 <= r < self.n and 0 <= b < self.n
                    and b not in (r, (r + 1) % self.n)
                    for r, b in self.backs.items()), "Invalid back target")
        for key in ("forward", "reverse", "degradation"):
            values = [Q(v) for v in getattr(self, key)]
            require(len(values) == self.n and all(v > 0 for v in values),
                    f"{key} must contain one positive value per species/reaction")
            setattr(self, key, values)
        # Column r records the full product side of reaction r.
        self.P = [[0] * self.n for _ in range(self.n)]
        for r, w in enumerate(self.weights):
            self.P[(r + 1) % self.n][r] += w
            if r in self.backs:
                self.P[self.backs[r]][r] += 1
        self.N = [[self.P[i][r] - int(i == r) for r in range(self.n)]
                  for i in range(self.n)]

    def products(self, z):
        return [math.prod(Q(z[i]) ** self.P[i][r] for i in range(self.n))
                for r in range(self.n)]

    def balance(self, z):
        require(len(z) == self.n and all(v > 0 for v in z),
                "This check expects a strictly positive concentration vector")
        products = self.products(z)
        current = [self.forward[r] * z[r] - self.reverse[r] * products[r]
                   for r in range(self.n)]
        production = matvec(self.N, current)
        loss = [d * v for d, v in zip(self.degradation, z)]
        return current, production, loss, [p - d for p, d in zip(production, loss)]


def six_species_certificate():
    """Table 1 and equations (8.2)-(8.4), transcribed as integers/fractions."""
    net = Network("Six-species example from Theorem 8.1", (2, 2, 1, 3, 1, 3),
                  {0: 5, 2: 1, 4: 3},
                  [704180941994460, 619799798320080, 1064148026967810,
                   1379889114510795, 5001118197194658, 1928841702241720],
                  [18981846318750, 1952234297034000, 3747998870714400,
                   18981846318750, 937378238937273, 98424792637500],
                  [4806051633136950, 18981846318750, 18981846318750,
                   18981846318750, 18981846318750, 2918522144328875])
    x = [Q(31, 50), Q(57, 50), Q(22, 25), Q(3, 5), Q(14, 25), Q(3, 5)]
    require(net.products(x) == [Q(9747, 12500), Q(484, 625), Q(171, 250),
                                Q(2744, 15625), Q(9, 25), Q(29791, 125000)],
            "Table 1 product monomials disagree")
    return net, x


def seven_species_certificate():
    """Table 2: note that the final two concentrations are 4/3 and 12."""
    D = 3546790543
    net = Network("Seven-species example from Theorem 8.3", (1,) * 7,
                  {0: 6, 2: 1, 3: 2, 4: 3, 5: 4, 6: 5},
                  [Q(79873083, D), Q(4870030864, 11 * D), Q(6049869075, 11 * D),
                   Q(562292665, 2 * D), Q(560635312, D), Q(76005216, D),
                   Q(80869 * 43, D)],
                  [Q(150591284, D), Q(1809648, D), Q(401938000, 11 * D),
                   Q(1809648, D), Q(486297815, D), Q(1809648, D), Q(1809648, D)],
                  [Q(72385920, D), Q(1809648, D), Q(413617671, 2 * D),
                   Q(6785880025, 22 * D), Q(558389511, 2 * D),
                   Q(1809648, D), Q(1809648, D)])
    require(sum(net.forward + net.reverse + net.degradation) == 1,
            "Table 2 normalization disagrees")
    require(min(net.forward + net.reverse + net.degradation) == Q(1809648, D),
            "Table 2 smallest constant disagrees")
    return net, [Q(1), Q(1, 16), Q(1, 9), Q(1, 5), Q(1, 3), Q(4, 3), Q(12)]


def symmetric_core(a=Q(3), b=Q(4), d=Q(1)):
    """A teaching example: three identical forks, each separated by one species."""
    return Network("Symmetric minimal six-species core", (1,) * 6,
                   {0: 5, 2: 1, 4: 3}, [a, b] * 3, [1] * 6, [d] * 6)


def symmetric_state(a, b, d):
    """Derived in the plan. Uniqueness + cyclic symmetry makes this exhaustive."""
    a, b, d = Q(a), Q(b), Q(d)
    require(min(a, b, d) > 0, "Positive kinetic and degradation constants required")
    h = a * b - d * (a + b + 1) - d * d
    if h <= 0:
        return None
    tail = h / (1 + 2 * d)
    fork = (b - d) * tail / (1 + 2 * d)
    require(min(fork, tail) > 0, "Positive branch formula failed")
    return [fork, tail] * 3


def strongly_connected(net, species, reactions):
    """Reachability in the forward species graph, using only retained products."""
    adjacency = {i: set() for i in species}
    for r in reactions:
        if r in species:
            adjacency[r].update(i for i in species if net.P[i][r])
    for start in species:
        seen, todo = {start}, [start]
        while todo:
            for j in adjacency[todo.pop()] - seen:
                seen.add(j)
                todo.append(j)
        if seen != set(species):
            return False
    return True


def check_small_core_minimality(net):
    """Exhaust all retained species/reaction sets for this small unit-weight example.

    A productive strongly connected component in any restriction is itself one
    of these retained sets. Thus checking connected sets with a two-product fork
    exhausts (Top) witnesses. Reverse currents belong to the kinetic model;
    the classified diluted source graph uses the displayed forward reactions.
    This is not a general minimal-core recognizer for other source conventions.
    """
    require(all(w == 1 for w in net.weights), "This finite check uses unit weights")
    full = set(range(net.n))
    require(strongly_connected(net, full, full), "Full source is not connected")
    require(any(sum(net.P[i][r] for i in full) >= 2 for r in full), "No fork")
    checked = 0
    for k in range(1, net.n + 1):
        for species in itertools.combinations(range(net.n), k):
            # Reactions with omitted sources are external feeds, not autonomous
            # reactions of the retained diluted source, so do not retain them.
            for mask in range(1, 1 << k):
                reactions = [species[j] for j in range(k) if mask >> j & 1]
                if set(species) == full and set(reactions) == full:
                    continue
                checked += 1
                if not any(sum(net.P[i][r] for i in species) >= 2 for r in reactions):
                    continue
                require(not strongly_connected(net, species, reactions),
                        f"Proper (Top) restriction found: {species}, {reactions}")
    return checked


def restriction_witness(net, species, reactions, flux):
    s = [[net.N[i][r] for r in reactions] for i in species]
    gain = matvec(s, flux)
    require(all(v > 0 for v in gain), "Restriction is not strictly productive")
    return {"species": species, "reactions": reactions, "matrix": s,
            "flux": flux, "production": gain}


def current_difference_check(net, rho):
    """Equations (5.1)-(5.7), on the teaching core at reference state y = 1.

    rho is a candidate second concentration vector. This checks the reduction
    at that candidate; it does not turn finite tests into a universal proof.
    """
    n = net.n
    p, q, e = net.forward, net.reverse, net.degradation
    require(not any(net.balance([Q(1)] * n)[3]), "Reference state not stationary")
    C = [[Q(0)] * n for _ in range(n)]
    for r, w in enumerate(net.weights):
        successor = (r + 1) % n
        secant = sum((rho[successor] ** j for j in range(w)), Q(0))
        if r in net.backs:
            back = net.backs[r]
            C[r][back] = 1
            C[r][successor] = rho[back] * secant
        else:
            C[r][successor] = secant
    offset = [v - 1 for v in rho]
    require(matvec(C, offset) == [v - 1 for v in net.products(rho)],
            "Exact product difference identity failed")
    A = [[p[r] * int(r == i) - q[r] * C[r][i] for i in range(n)] for r in range(n)]
    scaled_N = [[Q(net.N[i][r]) / e[i] for r in range(n)] for i in range(n)]
    AN = matmul(A, scaled_N)
    K = [[int(r == s) - AN[r][s] for s in range(n)] for r in range(n)]
    base_current = [u - v for u, v in zip(p, q)]
    defect = matvec(K, base_current)
    require(defect == [q[r] * (sum(C[r]) - 1) for r in range(n)],
            "Stationary-current comparison identity failed")
    forks, gaps = sorted(net.backs), [i for i in range(n) if i not in net.backs]
    block = lambda rows, cols: [[K[i][j] for j in cols] for i in rows]
    H = block(gaps, gaps)
    H_inverse = inverse(H)
    require(all(v >= 0 for row in H_inverse for v in row), "Gap inverse has wrong sign")
    correction = matmul(matmul(block(forks, gaps), H_inverse), block(gaps, forks))
    R = [[K[i][j] - correction[r][s] for s, j in enumerate(forks)]
         for r, i in enumerate(forks)]
    m = len(forks)
    margins = [R[j][j] - R[j][(j + 1) % m] for j in range(m)]
    require(all(R[j][(j - 1) % m] <= 0 and R[j][(j + 1) % m] > 0
                and margins[j] >= 1 for j in range(m)), "Fork inequalities failed")
    require(all(v > 0 for v in matvec(R, [base_current[i] for i in forks])),
            "Positive fork comparison failed")
    det_K, det_R = determinant(K), determinant(R)
    require(det_K == determinant(H) * det_R and det_K > 0, "Determinant check failed")
    delta = matvec(A, offset)
    residual = net.balance(rho)[3]
    # If rho were stationary, residual would vanish and K delta would be zero.
    require(matvec(K, delta) == [-v for v in matvec(A, [r / d for r, d in zip(residual, e)])],
            "Current-difference equation failed")
    require(any(residual), "Intended nonstationary candidate unexpectedly stationary")
    return {"candidate": rho, "C": C, "K": K, "gap_block": H, "fork_matrix": R,
            "diagonal_margins": margins, "determinant_K": det_K,
            "determinant_fork_matrix": det_R, "comparison": defect,
            "candidate_balance": residual}


def passive_chain_check():
    """Lemma 4.1 in a two-internal-species example; coefficients are exact."""
    H = [[5, -1], [-3, 7]]
    left = solve(H, [2, 0])
    right = solve(H, [0, 2])
    c, beta = 4 * left[1], right[0]
    leak_left, leak_right = 2 - left[0] - c, 2 - 4 * right[1] - beta
    require(c > 0 and beta > 0 and min(leak_left, leak_right) >= 0, "Tail is not passive")
    X, Y = Q(2), Q(1)
    internal = solve(H, [2 * X, 2 * Y])
    j_left, j_right = 2 * X - internal[0], 4 * internal[1] - 2 * Y
    require(j_left == c * X - beta * Y + leak_left * X, "Left endpoint formula failed")
    require(j_right == c * X - beta * Y - leak_right * Y, "Right endpoint formula failed")
    require(j_left - j_right == internal[0] + 2 * internal[1], "Internal loss mismatch")
    return {"c": c, "beta": beta, "left_loss_coefficient": leak_left,
            "right_loss_coefficient": leak_right, "internal_state": internal,
            "entering_current": j_left, "leaving_current": j_right}


def save_csv(path, header, rows):
    with path.open("w", newline="", encoding="utf-8") as f:
        writer = csv.writer(f)
        writer.writerow(header)
        writer.writerows(rows)


def vector_text(values):
    return "(" + ", ".join(str(v) for v in values) + ")"


def make_figures(out, certificates, boundary):
    import numpy as np
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    plt.rcParams.update({"font.family": "DejaVu Sans", "font.size": 11,
                         "axes.spines.top": False, "axes.spines.right": False,
                         "axes.edgecolor": "#a7b2b5", "axes.labelcolor": "#23363c",
                         "text.color": "#23363c", "xtick.color": "#44565b",
                         "ytick.color": "#44565b", "savefig.facecolor": "white"})
    teal, orange = "#087e8b", "#cb6632"
    versions = {"numpy": np.__version__, "matplotlib": matplotlib.__version__}

    fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), layout="constrained")
    for ax, (net, x) in zip(axes, certificates):
        idx = np.arange(net.n)
        ax.bar(idx - .19, np.ones(net.n), .34, color=teal, label="State y")
        ax.bar(idx + .19, [float(v) for v in x], .34, color=orange, label="State x")
        if net.n == 7:
            ax.set_yscale("log", base=2)
            ax.set_ylim(1 / 32, 32)
            ax.set_yticks([1 / 16, 1 / 4, 1, 4, 16], ["1/16", "1/4", "1", "4", "16"])
        else:
            ax.set_ylim(0, 1.42)
        for i, value in enumerate(x):
            ax.annotate(str(value), (i + .19, float(value)), xytext=(0, 5),
                        textcoords="offset points", ha="center", fontsize=9)
        ax.set_xticks(idx, [f"$X_{i}$" for i in idx])
        ax.set_ylabel("Concentration / concentration in state y")
        if net.n == 7:
            ax.set_ylabel("Concentration / concentration in state y\n(logarithmic scale)")
        ax.set_title(f"{net.n} species · two exactly verified steady states", loc="left", fontsize=12, pad=14)
        ax.grid(axis="y", alpha=.14)
        ax.set_axisbelow(True)
        ax.legend(frameon=False, loc="upper left" if net.n == 7 else "upper right")
    fig.suptitle("The same rate constants can support two different concentration states", fontsize=16)
    fig.savefig(out / "two_stationary_states.png", dpi=190)
    fig.savefig(out / "two_stationary_states.svg")
    plt.close(fig)

    fig, ax = plt.subplots(figsize=(9, 6.1), layout="constrained")
    a = np.array([row[0] for row in boundary])
    critical = np.array([row[1] for row in boundary])
    ax.fill_between(a, .05, critical, color="#d6eeee")
    ax.fill_between(a, critical, 2.5, color="#f1ece3")
    ax.plot(a, critical, color=teal, lw=2.3, label="Exact boundary formula")
    ax.scatter([3], [1], c=teal, edgecolors="white", s=100, zorder=4)
    ax.annotate("Default: every concentration = 1", (3, 1), xytext=(3.45, .75),
                fontsize=10, arrowprops={"arrowstyle": "-", "color": teal})
    ax.scatter([3], [1.5], marker="x", c=orange, s=75, zorder=4)
    ax.annotate("Comparison: no positive steady state", (3, 1.5), xytext=(.5, 1.8),
                fontsize=10, arrowprops={"arrowstyle": "-", "color": orange})
    ax.text(4.65, .32, "ONE positive\nsteady state", color=teal, fontsize=14, ha="center", weight="bold")
    ax.text(2.6, 2.15, "NO positive steady state", color="#806249", fontsize=14, ha="center", weight="bold")
    ax.set(xlim=(.25, 6), ylim=(.05, 2.5), xlabel="Forward rate constant at each fork, a",
           ylabel="Degradation constant for every species, d")
    ax.set_title("A minimal core: where does a positive steady state exist?\n"
                 "Six species · b = 4 · all reverse rate constants = 1", loc="left", pad=18)
    ax.text(.02, .98, "Boundary: 4a − d(a + 5) − d² = 0", transform=ax.transAxes,
            va="top", fontsize=11)
    fig.savefig(out / "positive_state_phase_map.png", dpi=190)
    fig.savefig(out / "positive_state_phase_map.svg")
    plt.close(fig)

    # The segment connects two exact states; it is not a time trajectory.
    fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), layout="constrained")
    line_rows = []
    for ax, (net, x), color in zip(axes, certificates, [teal, orange]):
        scale = max(net.forward + net.reverse + net.degradation)
        normalized = Network(net.name, net.weights, net.backs,
                             [p / scale for p in net.forward],
                             [q / scale for q in net.reverse],
                             [d / scale for d in net.degradation])
        errors = []
        for j in range(101):
            t = Q(j, 100)
            z = [1 + t * (v - 1) for v in x]
            residual = normalized.balance(z)[3]
            value = max(abs(v) for v in residual)
            errors.append(float(value))
            line_rows.append([net.n, str(t), str(value), float(value)])
        ax.plot(np.linspace(0, 1, 101), errors, color=color, lw=2.4,
                label=f"{net.n}-species example")
        ax.scatter([0, 1], [0, 0], color=color, s=45, clip_on=False)
        ax.set(xlabel="Position t on the segment from y to x",
               ylabel="Largest absolute concentration derivative\n(after a common rescaling of all rates)")
        ax.set_title(f"{net.n}-species example", loc="left", fontsize=12, pad=12)
        ax.grid(alpha=.13)
    fig.suptitle("Both endpoints balance exactly; intermediate concentrations need not balance", fontsize=15)
    fig.savefig(out / "between_stationary_states.png", dpi=190)
    fig.savefig(out / "between_stationary_states.svg")
    plt.close(fig)
    save_csv(out / "between_states.csv", ["species_count", "t_exact", "max_derivative_exact", "max_derivative_decimal"], line_rows)
    return versions


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--output-dir", type=Path, default=Path(__file__).resolve().parent / "outputs")
    parser.add_argument("--no-figures", action="store_true", help="Run exact checks without NumPy or Matplotlib")
    args = parser.parse_args()
    out = args.output_dir
    out.mkdir(parents=True, exist_ok=True)
    report, data, balances = [], {}, []
    emit = report.append
    emit("TYPE II CORES: EXACT CHECKS AND WORKED EXAMPLES")
    emit("Paper: type_II_l_unistationarity_arxiv.pdf, 7 September 2026")
    emit("All reported balance zeros use exact rational arithmetic.")
    certificates = [six_species_certificate(), seven_species_certificate()]
    for net, x in certificates:
        states = {"y": [Q(1)] * net.n, "x": x}
        require(x != states["y"] and all(v > 0 for v in x), "States must be positive and distinct")
        emit(f"\n{net.name}")
        emit("  x = " + vector_text(x))
        emit("  All forward, reverse and degradation constants: strictly positive")
        state_data = {}
        for name, z in states.items():
            current, production, loss, residual = net.balance(z)
            require(residual == [0] * net.n, f"{net.name}, state {name}: nonzero balance")
            emit(f"  State {name}: exact balance residual = {vector_text(residual)}")
            state_data[name] = {"concentrations": z, "current": current,
                                "net_reaction_production": production, "degradation_loss": loss,
                                "balance_residual": residual}
            for i in range(net.n):
                balances.append([net.n, name, i, z[i], production[i], loss[i], residual[i]])
        data[f"certificate_{net.n}"] = {"name": net.name, "P": net.P, "N": net.N,
                                        "forward": net.forward, "reverse": net.reverse,
                                        "degradation": net.degradation, "states": state_data}
    six, seven = certificates[0][0], certificates[1][0]
    witnesses = [restriction_witness(six, [1, 2], [1, 2], [2, 3]),
                 restriction_witness(seven, [2, 3, 4], [2, 3, 4], [2, 3, 2])]
    data["proper_autocatalytic_restrictions"] = witnesses
    emit("\nWhy these examples are not minimal cores:")
    for witness in witnesses:
        emit("  Retain species " + ", ".join(f"X{i}" for i in witness['species'])
             + ": flux " + vector_text(witness['flux'])
             + " gives production " + vector_text(witness['production']))
    # Lemma 3.2's formula at M=2 is (3,5), a different valid witness.
    require(matvec(witnesses[0]["matrix"], [3, 5]) == [2, 1], "Lemma 3.2 witness failed")
    data["lemma_3_2_length_two_witness"] = {"flux": [3, 5], "production": [2, 1]}

    core = symmetric_core()
    count = check_small_core_minimality(core)
    require(core.balance([Q(1)] * 6)[3] == [0] * 6, "Teaching core default balance failed")
    core_gain = matvec(core.N, [2, 3] * 3)
    require(core_gain == [1] * 6, "Full core production witness failed")
    data["teaching_core"] = {"P": core.P, "N": core.N, "forward": core.forward,
                             "reverse": core.reverse, "degradation": core.degradation,
                             "proper_restrictions_checked": count, "productive_flux": [2, 3] * 3,
                             "production": core_gain}
    emit(f"\nTeaching core: {count} proper species/reaction restrictions checked; none retain (Top).")
    emit("  At a=3, b=4, d=1, the state (1,1,1,1,1,1) balances exactly.")
    emit("  Its uniqueness follows from Theorem 1.1 after checking the core assumptions.")
    presets = [(3, 4, 1), (5, 4, 1), (3, 4, Q(3, 2)), (2, 4, 1)]
    data["phase_presets"] = []
    for a, b, d in presets:
        z = symmetric_state(a, b, d)
        if z is not None:
            require(symmetric_core(a, b, d).balance(z)[3] == [0] * 6, "Symmetric formula failed")
        data["phase_presets"].append({"a": a, "b": b, "d": d, "positive_state": z})
        emit(f"  a={a}, b={b}, d={d}: " + ("positive state " + vector_text(z) if z else "no positive steady state"))

    kernel = current_difference_check(core, certificates[0][1])
    data["current_difference_example"] = kernel
    emit("\nExact current-difference calculation on the teaching core:")
    emit(f"  det(K) = {kernel['determinant_K']} > 0")
    emit(f"  det(reduced fork matrix) = {kernel['determinant_fork_matrix']} > 0")
    emit("  This rejects the selected candidate; the paper proves nonsingularity generally.")
    tail = passive_chain_check()
    data["passive_chain_example"] = tail
    emit(f"\nUnit-chain elimination: c={tail['c']}, beta={tail['beta']}, "
         f"left loss={tail['left_loss_coefficient']}, right loss={tail['right_loss_coefficient']}.")

    # Classify a rational grid exactly. Only drawing the boundary uses floats.
    grid = []
    for ia in range(25, 601, 5):
        for id_ in range(5, 251, 5):
            a, d = Q(ia, 100), Q(id_, 100)
            z = symmetric_state(a, 4, d)
            if z is not None:
                require(symmetric_core(a, 4, d).balance(z)[3] == [0] * 6, "Grid balance failed")
            grid.append([a, d, int(z is not None), z[0] if z else "", z[1] if z else ""])
    from math import sqrt
    boundary = []
    for j in range(576):
        a = (25 + j) / 100
        s = a + 5
        boundary.append([a, 8 * a / (sqrt(s * s + 16 * a) + s)])
    save_csv(out / "exact_balances.csv", ["species_count", "state", "species_index", "concentration", "net_reaction_production", "degradation_loss", "balance_residual"], balances)
    save_csv(out / "phase_grid.csv", ["a_exact", "d_exact", "positive_state_count", "fork_concentration_exact", "intermediate_concentration_exact"], grid)
    save_csv(out / "phase_boundary.csv", ["a_decimal", "critical_d_decimal"], boundary)
    versions = {} if args.no_figures else make_figures(out, certificates, boundary)
    emit(f"\nPhase map: {len(grid)} parameter pairs classified exactly; every displayed positive state checked.")
    emit("  Entire regions follow from the derived formula + the paper's uniqueness theorem.")
    emit("\nScope: two counterexamples verified, not a count of all their steady states.")
    emit("No stability, convergence, or thermodynamic-consistency conclusion is asserted.")
    emit("The Lean sources were not supplied or compiled by this program.")
    emit("ALL CHECKS PASSED")
    text_report = "\n".join(report) + "\n"
    (out / "run_output.txt").write_text(text_report, encoding="utf-8")
    data["provenance"] = {"paper": "type_II_l_unistationarity_arxiv.pdf", "paper_date": "2026-09-07",
                          "python": platform.python_version(), "plotting_versions": versions,
                          "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
                          "arithmetic": "exact rational for all checks; floating point for drawing only",
                          "phase_map_scope": "symmetric unit-weight six-species core; theorem + derived formula"}
    (out / "results.json").write_text(json.dumps(data, indent=2, default=str) + "\n", encoding="utf-8")
    print(text_report, end="")


if __name__ == "__main__":
    main()