Example code

"""Construct globally optimal stationary mass-action production; keep scope explicit.

Run python example.py --output outputs. All concentrations and times are
dimensionless source-model quantities, not fitted laboratory parameters.
"""
from __future__ import annotations

# EDITABLE INPUTS -------------------------------------------------------------
RECYCLING_ORDER = 2
RESPONSE_RATIO = '3/2'          # 1 < r < 2; reaction response f=(1,m*r)
REQUIRED_CURRENT = '1'
REVERSE_BUDGETS = ('2', '6')    # reverse one-way fluxes at the design state
EQUILIBRIUM_CONSTANTS = ('9/4', '14/9')
INITIAL_INTERNAL_CONCENTRATIONS = (0.1, 0.5, 2.0, 4.0)
RELAXATION_END = 1.5
CURVE_POINTS = 241
MANUSCRIPT_SHA256 = 'bcdd5807347a31ab417f017f3dc921ae79bbf33d3cc9a184b3ed089b73e0720b'
# ---------------------------------------------------------------------------

import argparse
import csv
from dataclasses import dataclass
import hashlib
import json
import math
from pathlib import Path
import platform

import numpy as np
import sympy as sp
from scipy.integrate import solve_ivp
from sympy.solvers.simplex import InfeasibleLPError, lpmin


def rational(x):
    """Avoid the binary-float interpretation of user-entered decimals."""
    return sp.Rational(str(x))


def vector(xs):
    return sp.Matrix([rational(x) for x in xs])


def strings(xs):
    return [str(sp.simplify(x)) for x in xs]


@dataclass(frozen=True)
class ResponseProfile:
    f: sp.ImmutableMatrix
    q: sp.ImmutableMatrix
    u: sp.ImmutableMatrix
    phi: sp.Expr


class SquareSource:
    """Species rows, reaction columns. X is chemostatted; others are internal.

    This implements only the paper's square source class. It does not silently
    apply the result to rectangular networks or mixed-sign production modes.
    """
    def __init__(self, reactants, products, controlled=0, names=None):
        self.Sp = sp.ImmutableMatrix(reactants)
        self.Sm = sp.ImmutableMatrix(products)
        self.n = self.Sp.rows
        if self.Sp.shape != (self.n, self.n) or self.Sm.shape != self.Sp.shape:
            raise ValueError('Both stoichiometric arrays must be square and equal-sized.')
        if any(not x.is_Integer or x < 0 for x in (*self.Sp, *self.Sm)):
            raise ValueError('Mass-action stoichiometries must be nonnegative integers.')
        if not 0 <= controlled < self.n:
            raise ValueError('Invalid controlled species.')
        self.X = controlled
        self.names = tuple(names or [f'z{i}' for i in range(self.n)])
        self.S = self.Sm-self.Sp
        if self.Sp.det() == 0 or self.S.det() == 0:
            raise ValueError('Reactant matrix and net stoichiometry must be invertible.')
        self.T = self.Sp.inv()*self.Sm
        self.g = self.S.inv()*sp.eye(self.n)[:, self.X]
        if any(x < 0 for x in self.T) or any(x <= 0 for x in self.g):
            raise ValueError('Source theorem requires T>=0 and g>0.')
        self.N = sp.ilcm(*[int(x.q) for x in self.g]) if self.n > 1 else self.g[0].q
        self.weights = self.N*self.g

    def profile(self, raw_f):
        f = vector(raw_f)
        if len(f) != self.n or any(x <= 0 for x in f):
            raise ValueError('Every reaction response must be positive.')
        h = self.T.T*f
        if any(h[i] <= f[i] for i in range(self.n)):
            raise ValueError('Forward response requires T.T*f > f componentwise.')
        u = self.Sp.T.inv()*f
        if u[self.X] <= 0:
            raise AssertionError('Source-mode identity violated.')
        f = f/u[self.X]
        u = u/u[self.X]
        q = sp.Matrix([h[i]/raw for i, raw in enumerate(vector(raw_f))])
        phi = sp.prod(q[i]**self.weights[i] for i in range(self.n))
        return ResponseProfile(sp.ImmutableMatrix(f), sp.ImmutableMatrix(q),
                               sp.ImmutableMatrix(u), sp.factor(phi))

    def components(self):
        """SCCs of i -> j iff T[j,i]>0; exact transitive-closure implementation."""
        reach = [[i == j or self.T[j, i] > 0 for j in range(self.n)] for i in range(self.n)]
        for k in range(self.n):
            for i in range(self.n):
                for j in range(self.n):
                    reach[i][j] = reach[i][j] or (reach[i][k] and reach[k][j])
        pending = set(range(self.n)); components = []
        while pending:
            i = min(pending)
            c = tuple(j for j in sorted(pending) if reach[i][j] and reach[j][i])
            components.append(c); pending.difference_update(c)
        terminal = [c for c in components if not any(self.T[j, i] > 0
                    for i in c for j in range(self.n) if j not in c)]
        return components, terminal

    def uniqueness_certificate(self):
        components, terminal = self.components()
        if len(terminal) != 1:
            return {'certified': False, 'reason': 'Graph is not rooted; criterion not applicable.'}
        blocks = []
        for c in components:
            if c == terminal[0]:
                continue
            matrix = sp.eye(len(c))-self.T.extract(c, c)
            cert = HomogeneousCone(matrix).solve()
            blocks.append({'indices': list(c), 'subcritical_certificate': cert})
        ok = all(b['subcritical_certificate']['feasible'] for b in blocks)
        return {'certified': ok, 'terminal': list(terminal[0]), 'transient_blocks': blocks,
                'reason': 'Rooted with positive subcritical block vectors.' if ok else
                          'Sufficient uniqueness criterion fails; this alone does not prove nonuniqueness.'}


class HomogeneousCone:
    """Decide D*f>=0, f>=1 with a rational witness on either side.

    A solver verdict alone is never emitted as the mathematical certificate.
    The infeasible side gives y>=0, y.T*D<=0 with a strictly negative sum.
    """
    def __init__(self, matrix):
        self.D = sp.Matrix(matrix)

    def solve(self):
        x = sp.Matrix(sp.symbols(f'x0:{self.D.cols}'))
        try:
            _, solution = lpmin(0, [v >= 1 for v in x]+[v >= 0 for v in self.D*x])
            f = sp.Matrix([solution.get(v, 1) for v in x])
            residual = self.D*f
            if not (all(v >= 1 for v in f) and all(v >= 0 for v in residual)):
                raise ArithmeticError('Exact primal verification failed.')
            return {'feasible': True, 'f': strings(f), 'D_f': strings(residual)}
        except InfeasibleLPError:
            y = sp.Matrix(sp.symbols(f'y0:{self.D.rows}'))
            row = y.T*self.D
            _, solution = lpmin(0, [v >= 0 for v in y]+[v <= 0 for v in row]+[sum(row) <= -1])
            witness = sp.Matrix([solution.get(v, 0) for v in y])
            row = witness.T*self.D
            if not (all(v >= 0 for v in witness) and all(v <= 0 for v in row) and sum(row) < 0):
                raise ArithmeticError('Exact dual verification failed.')
            return {'feasible': False, 'dual_y': strings(witness), 'y_D': strings(row)}


class ReverseBudget:
    def __init__(self, source, budgets):
        self.source = source
        self.beta = vector(budgets)
        if len(self.beta) != source.n or any(x <= 0 for x in self.beta):
            raise ValueError('A positive reverse one-way budget is required for each reaction.')

    def certificate(self, required_current):
        J = rational(required_current)
        if J <= 0:
            raise ValueError('Required current must be positive.')
        D = self.source.T.T-sp.diag(*[1+J*self.source.g[i]/self.beta[i] for i in range(self.source.n)])
        result = HomogeneousCone(D).solve()
        result['D'] = [[str(x) for x in D.row(i)] for i in range(D.rows)]
        if result['feasible']:
            profile = self.source.profile(result['f'])
            reverse = [J*self.source.g[i]/(profile.q[i]-1) for i in range(self.source.n)]
            if any(reverse[i] > self.beta[i] for i in range(self.source.n)):
                raise AssertionError('Budget witness did not reconstruct admissible fluxes.')
            result.update(q=strings(profile.q), reverse_one_way=strings(reverse))
        return result


class MassActionRealization:
    """Composable rate model with exact expressions and numerical ODE access."""
    def __init__(self, source, profile, current=1, equilibrium_constants=None):
        self.source = source
        self.profile = profile
        self.J0 = rational(current)
        if self.J0 <= 0:
            raise ValueError('Design current must be positive.')
        self.reference = sp.ones(source.n, 1)
        if equilibrium_constants is not None:
            K = vector(equilibrium_constants)
            if len(K) != source.n or any(k <= 0 for k in K):
                raise ValueError('Equilibrium constants must be positive.')
            logz = source.S.T.inv()*sp.Matrix([sp.log(K[i]/profile.q[i]) for i in range(source.n)])
            self.reference = logz.applyfunc(lambda v: sp.simplify(sp.exp(sp.expand_log(v, force=True))))
        self.reverse = sp.Matrix([self.J0*source.g[i]/(profile.q[i]-1) for i in range(source.n)])
        self.forward = sp.matrix_multiply_elementwise(profile.q, self.reverse)
        self.kp = sp.Matrix([sp.simplify(self.forward[i]/self.monomial(self.reference, source.Sp[:, i])) for i in range(source.n)])
        self.km = sp.Matrix([sp.simplify(self.reverse[i]/self.monomial(self.reference, source.Sm[:, i])) for i in range(source.n)])
        self.z = sp.Matrix(sp.symbols(f'z0:{source.n}', positive=True))
        self.symbolic_currents = self.currents(self.z)
        self.symbolic_field = source.S*self.symbolic_currents
        self._numeric = sp.lambdify([tuple(self.z)], self.symbolic_field, 'numpy')

    @staticmethod
    def monomial(z, exponents):
        return sp.prod(z[j]**exponents[j] for j in range(len(z)))

    def currents(self, z):
        s = self.source
        return sp.Matrix([self.kp[i]*self.monomial(z, s.Sp[:, i])-
                          self.km[i]*self.monomial(z, s.Sm[:, i]) for i in range(s.n)])

    def field(self, z):
        a = np.asarray(z, dtype=float)
        if a.shape != (self.source.n,) or not np.all(np.isfinite(a)):
            raise ValueError('Supply one finite concentration per species.')
        return np.asarray(self._numeric(a), dtype=float).reshape(-1)

    def bottleneck(self, z):
        """Numerical audit of the all-state inequality; only stationarity equalizes C_i."""
        z = np.asarray(z, float)
        ref = np.array(self.reference, float).reshape(-1)
        if z.shape != ref.shape or np.any(z <= 0) or not np.all(np.isfinite(z)):
            raise ValueError('Bottleneck comparison requires a positive finite state.')
        xi = np.array(self.source.Sp.T, float)@np.log(z/ref)
        f = np.array(self.profile.f, float).reshape(-1)
        q = np.array(self.profile.q, float).reshape(-1)
        h = np.array(self.source.T.T, float)@xi
        C = (q*np.exp(xi)-np.exp(h))/(q-1)
        k = int(np.argmin(xi/f))
        delta = h[k]-q[k]*xi[k]
        convex = (math.expm1(q[k]*xi[k])-q[k]*math.expm1(xi[k]))/(q[k]-1)
        excess = math.exp(q[k]*xi[k])*math.expm1(delta)/(q[k]-1)
        return {'reaction': k, 'normalized_currents': C.tolist(), 'delta': float(delta),
                'gap': float(1-C[k]), 'convex_gap': convex, 'coupling_gap': excess}

    def local_curvature(self):
        """Exact current-branch curvature at a rooted regular maximum.

        The nullspace and bordered rank are explicitly checked. This is not a
        claim about the stability of the controlled dynamical system.
        """
        s = self.source; p = self.profile
        M = sp.diag(*self.reverse)*(sp.diag(*p.q)-s.T.T)*s.Sp.T
        left = M.T.nullspace()
        if len(left) != 1 or M.row_join(-s.g).rank() != s.n:
            raise ValueError('The one-dimensional regular-branch hypothesis failed.')
        lam = left[0]
        if (lam.T*s.g)[0] < 0:
            lam = -lam
        denom = (lam.T*s.g)[0]
        if denom <= 0 or any(x < 0 for x in lam):
            raise ValueError('No nonnegative regular left-null certificate.')
        value = -sum(lam[i]*self.reverse[i]*p.q[i]*(p.q[i]-1)*p.f[i]**2 for i in range(s.n))/denom
        return {'J_second_log_control': str(sp.factor(value)), 'left_null': strings(lam),
                'regular': True}

    def controlled_trajectory(self, internal_initial, end=RELAXATION_END, points=201):
        s = self.source; internal = [i for i in range(s.n) if i != s.X]
        initial = np.asarray(internal_initial, float)
        if initial.shape != (len(internal),) or np.any(initial <= 0) or end <= 0:
            raise ValueError('Positive initial internal concentrations and end time required.')
        fixed = np.array(self.reference, float).reshape(-1)
        def rhs(_, y):
            z = fixed.copy(); z[internal] = y
            return self.field(z)[internal]
        sol = solve_ivp(rhs, (0, end), initial, t_eval=np.linspace(0, end, points),
                        method='Radau', rtol=2e-10, atol=2e-12)
        if not sol.success or np.min(sol.y) <= 0:
            raise ArithmeticError('Numerical controlled trajectory failed.')
        return sol.t, sol.y.T


class RecyclingFamily:
    """X <-> Y and mY <-> 2X+(m-1)Y, with X controlled."""
    def __init__(self, m=2):
        if not isinstance(m, int) or m < 2:
            raise ValueError('Integer recycling order m>=2 required.')
        self.m = m
        self.source = SquareSource([[1, 0], [0, m]], [[0, 2], [1, m-1]], names=('X', 'Y'))
        self.capacity = 1+sp.Rational(1, m)

    def realize(self, r, current=1, equilibrium_constants=None):
        r = rational(r)
        if not 1 < r < 2:
            raise ValueError('The source cone is exactly 1<r<2.')
        return MassActionRealization(self.source, self.source.profile([1, self.m*r]), current, equilibrium_constants)

    def budget_optimum(self, current, budgets):
        J = rational(current); b1, b2 = vector(budgets)
        if min(J, b1, b2) <= 0:
            raise ValueError('Positive current and budgets required.')
        lo = 1+J/b1; hi = 2*b2/(b2+self.m*J)
        discriminant = (b2+self.m*b1)**2+4*self.m*b1*b2
        ceiling = 2*b1*b2/(sp.sqrt(discriminant)+b2+self.m*b1)
        return {'feasible': bool(lo <= hi), 'r_min': str(lo), 'r_max': str(hi),
                'minimum_phi': str(self.capacity+(1-sp.Rational(1, self.m))*J/b1) if lo <= hi else None,
                'maximum_current': str(sp.simplify(ceiling)), 'maximum_current_numeric': float(ceiling)}

    def stationary_curve(self, model, ys):
        """All positive stationary states parameterized by internal Y>0.

        The quadratic has exactly one positive root. No branch is chosen by
        continuing from the design point; some positive states have negative J.
        """
        kp = np.array(model.kp, float).reshape(-1); km = np.array(model.km, float).reshape(-1)
        ys = np.asarray(ys, float)
        if np.any(ys <= 0):
            raise ValueError('Positive Y required.')
        a = km[1]*ys**(self.m-1); b = kp[0]; c = km[0]*ys+kp[1]*ys**self.m
        xs = 2*c/(b+np.sqrt(b*b+4*a*c))
        js = kp[0]*xs-km[0]*ys
        return np.column_stack([xs, ys, js])


def rooted_example():
    """Exact counterexample: rootedness alone permits two global maxima."""
    source = SquareSource([[0, 1, 0], [2, 1, 0], [4, 0, 1]],
                          [[1, 0, 0], [1, 4, 0], [1, 8, 2]], names=('X', 'Y', 'Z'))
    model = MassActionRealization(source, source.profile([1, 1, 1]))
    second = sp.Matrix([5/sp.sqrt(3), sp.sqrt(3), 1])
    ratios = sp.Matrix([sp.simplify(model.kp[i]*model.monomial(second, source.Sp[:, i])/
                        (model.km[i]*model.monomial(second, source.Sm[:, i]))) for i in range(3)])
    second_response = source.Sp.T*sp.Matrix([1, 17, -6])
    internal_jac = model.symbolic_field.jacobian(model.z).extract([1, 2], [1, 2])
    unit_jac = internal_jac.subs(dict(zip(model.z, [1, 1, 1])))
    second_jac = internal_jac.subs(dict(zip(model.z, second)))
    p1, p2, p3, t = sp.symbols('p1 p2 p3 t')
    equations = [2*p1-p2*p3-t, 2*p2-p1**2-t, 2*p3-p3**2-t]
    factor_at_max = sp.factor(equations[1].subs({p2: 2*p1-1, t: 1}))
    return model, {
        'T': [[str(x) for x in source.T.row(i)] for i in range(3)],
        'g': strings(source.g), 'weights': strings(source.weights),
        'uniqueness_criterion': source.uniqueness_certificate(),
        'maximizers': [strings(sp.ones(3, 1)), strings(second)],
        'maximizer_currents': [strings(model.currents(sp.ones(3, 1))), strings(model.currents(second))],
        'unit_phi': str(model.profile.phi), 'unattained_forward_capacity': '4',
        'second_flux_ratios': strings(ratios),
        'second_phi': str(sp.prod(ratios[i]**source.weights[i] for i in range(3))),
        'second_response_u': ['1', '17', '-6'], 'second_response_f': strings(second_response),
        'capacity_scope': 'The second maximum has a negative reaction response and is outside the forward-response cone.',
        'stationary_monomial_equations': strings(equations),
        'factor_at_maximum': str(factor_at_max),
        'unit_internal_jacobian': [[str(x) for x in unit_jac.row(i)] for i in range(2)],
        'unit_internal_determinant': str(unit_jac.det()),
        'second_internal_determinant': str(sp.simplify(second_jac.det())),
        'local_curvature': model.local_curvature(),
        'dynamical_scope': 'The unit design is an internal saddle. Maxima occur at different X; this is not bistability at one control.'}


def write_csv(out, name, headers, rows):
    with (out/name).open('w', newline='') as stream:
        writer = csv.writer(stream); writer.writerow(headers); writer.writerows(rows)


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--output', type=Path, default=Path('outputs'))
    out = parser.parse_args().output; out.mkdir(parents=True, exist_ok=True)
    family = RecyclingFamily(RECYCLING_ORDER)
    model = family.realize(RESPONSE_RATIO, REQUIRED_CURRENT)
    fixed = family.realize(RESPONSE_RATIO, REQUIRED_CURRENT, EQUILIBRIUM_CONSTANTS)
    budget = ReverseBudget(family.source, REVERSE_BUDGETS)
    optimum = family.budget_optimum(REQUIRED_CURRENT, REVERSE_BUDGETS)
    _, rooted = rooted_example()
    result = {
        'inputs': {'m': RECYCLING_ORDER, 'r': RESPONSE_RATIO, 'J': REQUIRED_CURRENT,
                   'reverse_budgets': REVERSE_BUDGETS, 'equilibrium_constants': EQUILIBRIUM_CONSTANTS},
        'unit_design': {'kp': strings(model.kp), 'km': strings(model.km), 'g': strings(family.source.g),
                        'f': strings(model.profile.f), 'u': strings(model.profile.u), 'q': strings(model.profile.q),
                        'phi': str(model.profile.phi), 'unattained_capacity': str(family.capacity),
                        'uniqueness': family.source.uniqueness_certificate(), 'curvature': model.local_curvature()},
        'budget_certificate': budget.certificate(REQUIRED_CURRENT), 'budget_optimum': optimum,
        'higher_current_certificate': budget.certificate(rational(REQUIRED_CURRENT)*sp.Rational(6, 5)),
        'fixed_equilibrium_design': {'reference': strings(fixed.reference), 'kp': strings(fixed.kp), 'km': strings(fixed.km)},
        'rooted_counterexample': rooted,
        'scope': 'Exact finite algebra/LP witnesses plus numerical illustrations. Global theorem is explained, not proved by sampling. Lean is not rerun.'}
    (out/'results.json').write_text(json.dumps(result, indent=2)+'\n')
    curve = family.stationary_curve(model, np.linspace(.03, 3, CURVE_POINTS))
    write_csv(out, 'stationary_curve.csv', ['X', 'Y', 'stationary_current'], curve)
    turnover = []
    for gap in np.geomspace(1e-5, .95, 101):
        r = 1+gap
        q2 = 2/(family.m*r)+(family.m-1)/family.m
        phi = r*q2
        bounded_J = min(1, gap, q2-1)
        turnover.append([r, phi, phi-float(family.capacity), float(model.J0)/gap,
                         float(model.J0)/(q2-1), bounded_J, bounded_J/gap, bounded_J/(q2-1)])
    write_csv(out, 'capacity_turnover.csv', ['r', 'phi', 'capacity_gap', 'fixed_J_reverse1', 'fixed_J_reverse2',
                                           'bounded_rate_J', 'bounded_rate_reverse1', 'bounded_rate_reverse2'], turnover)
    sweep = []
    for J in np.linspace(.01, float(optimum['maximum_current_numeric'])*1.2, 121):
        b = family.budget_optimum(str(J), REVERSE_BUDGETS)
        sweep.append([J, float(rational(b['r_min'])), float(rational(b['r_max'])), b['feasible'],
                      float(rational(b['minimum_phi'])) if b['feasible'] else ''])
    write_csv(out, 'budget_sweep.csv', ['required_current', 'r_min', 'r_max', 'feasible', 'minimum_phi'], sweep)
    trajectories = []
    for initial in INITIAL_INTERNAL_CONCENTRATIONS:
        times, values = model.controlled_trajectory([initial])
        trajectories.extend([[initial, t, row[0]] for t, row in zip(times, values)])
    write_csv(out, 'controlled_trajectories.csv', ['initial_Y', 'time', 'Y'], trajectories)
    lines = [f'Recycling m={family.m}: capacity infimum {family.capacity}; design Phi={model.profile.phi}.',
             f'Unit design kp={strings(model.kp)}, km={strings(model.km)}, current={model.J0}.',
             f'Reverse-budget feasible: {result["budget_certificate"]["feasible"]}; optimum={optimum["minimum_phi"]}.',
             f'Fixed equilibrium constants give reference state {strings(fixed.reference)}.',
             f'Rooted counterexample: two global maxima; unit internal determinant={rooted["unit_internal_determinant"]}.',
             f'Second maximum Phi={rooted["second_phi"]}<4 but f={rooted["second_response_f"]} is not forward.',
             'Stationary production, uniqueness, and dynamical attraction are separate questions.',
             'Exact algebraic witnesses and numerical trajectories; no Lean proof rerun.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n'); print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig, axes = plt.subplots(1, 2, figsize=(10, 4), layout='constrained')
    axes[0].plot(curve[:, 0], curve[:, 2], color='#216b91')
    axes[0].scatter([1], [float(model.J0)], color='#b95024', zorder=3, label='Constructed optimum')
    axes[0].axhline(float(model.J0), color='black', ls='--', lw=1)
    axes[0].set(xlabel='Controlled concentration X', ylabel='Stationary production current J', title='Global stationary optimum')
    axes[0].legend(fontsize=8)
    for initial in INITIAL_INTERNAL_CONCENTRATIONS:
        rows = np.array([r for r in trajectories if r[0] == initial])
        axes[1].plot(rows[:, 1], rows[:, 2], label=f'Y(0)={initial:g}')
    axes[1].axhline(1, color='black', ls='--', lw=1)
    axes[1].set(xlabel='Time with X held at 1', ylabel='Internal concentration Y', title='Relaxation with X held at the production\noptimum')
    axes[1].legend(fontsize=8)
    for ax in axes: ax.grid(alpha=.2)
    fig.savefig(out/'production.png', dpi=180); fig.savefig(out/'production.svg'); plt.close(fig)
    fig, axes = plt.subplots(1, 2, figsize=(10, 4), layout='constrained')
    rows = np.array(turnover)
    axes[0].loglog(rows[:, 2], rows[:, 3], label='Reverse reaction 1')
    axes[0].loglog(rows[:, 2], rows[:, 4], label='Reverse reaction 2')
    axes[0].set(xlabel='Capacity gap Phi - L', ylabel='Reverse one-way flux at fixed J', title='Reverse flux required near the affinity\nlimit')
    axes[0].legend(fontsize=8)
    xx = [r[0] for r in sweep]
    axes[1].plot(xx, [r[1] for r in sweep], label='Lower response bound')
    axes[1].plot(xx, [r[2] for r in sweep], label='Upper response bound')
    axes[1].fill_between(xx, [r[1] for r in sweep], [r[2] for r in sweep],
                         where=[r[3] for r in sweep], alpha=.15, color='#216b91', label='Feasible responses')
    axes[1].axvline(optimum['maximum_current_numeric'], color='black', ls='--', lw=1)
    axes[1].set(xlabel='Required production current J', ylabel='Response ratio r', title='Finite reverse-flux budgets')
    axes[1].legend(fontsize=8)
    for ax in axes: ax.grid(alpha=.2)
    fig.savefig(out/'budgets.png', dpi=180); fig.savefig(out/'budgets.svg'); plt.close(fig)
    digest = lambda p: hashlib.sha256(p.read_bytes()).hexdigest()
    (out/'run_metadata.json').write_text(json.dumps({'manuscript_sha256': MANUSCRIPT_SHA256,
        'source_sha256': digest(Path(__file__)), 'python': platform.python_version(),
        'output_sha256': {p.name: digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name != 'run_metadata.json'}}, indent=2)+'\n')


if __name__ == '__main__':
    main()