Example code

"""Shared NADPH: retained-state kinetics, exact thresholds, and a usable design model.

Run: python example.py --output outputs. Concentrations uM; time seconds.
The main inputs are immediately below. The exact paper audit stays fixed and is
reported separately from any edited preparation. See README for evidence scope.
"""
from __future__ import annotations
import argparse
import csv
from dataclasses import dataclass, replace, asdict
from fractions import Fraction as F
import hashlib
import json
import math
from pathlib import Path
import platform
from typing import Protocol
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import brentq

# EDITABLE PREPARATION: paper Table 1, source-derived effective kinetics.
# a_G, a_T, b_T already incorporate a maintained peroxide concentration 0.01 uM.
# Editing that challenge requires changing the peroxide-dependent rates too.
GPX = dict(G=371.56, z0=1.78, E=50., a=.21, b=.04, c=10., k=3.2)
TRX = dict(T=.505, z0=.075, E=19.096, a=.4, b=.00072, c=.003, d=15., e=2.1, k=20.)
SOURCE = dict(P=30.3, P0=.3, K=57., V=375.)
QUOTAS = (10., 4.)  # design requirements, uM/s; not health thresholds
REFERENCE_SCALE = .1
SLACK_SCALE = .12
SOURCE_TOLERANCE = .05  # family of constant multiplicative settings
HORIZON = 2000.
CARRIER_PERTURBATION = (1., .02)  # increases to z and zT, uM
ENZYME_PERTURBATION = (.2, .05)  # increases to e1 and h, uM
MANUSCRIPT_SHA256 = 'b1f97e45d58820881a8f2d779660a9f6e27270fadf9df8596477b2a7c787becd'


def positive(**values):
    if any(not math.isfinite(float(v)) or v <= 0 for v in values.values()):
        raise ValueError('Finite positive values required: '+', '.join(values))


@dataclass(frozen=True)
class QuadraticResponse:
    """Au²+(B+C/x)u-Gamma=0; j=pu/(vu+w). B may be negative."""
    A: float
    B: float
    C: float
    Gamma: float
    p: float
    v: float
    w: float

    def __post_init__(self):
        positive(A=self.A, C=self.C, Gamma=self.Gamma, p=self.p, w=self.w)
        if not math.isfinite(float(self.B)) or not math.isfinite(float(self.v)) or self.v < 0:
            raise ValueError('B finite; v finite and nonnegative')

    def polynomial(self, x, u):
        return self.A*u*u+(self.B+self.C/x)*u-self.Gamma

    def carrier(self, x):
        if x == 0: return 0.
        positive(x=x)
        beta=float(self.B+self.C/x)
        disc=math.hypot(beta, 2*math.sqrt(float(self.A*self.Gamma)))
        # Avoid cancellation and also avoid overflow in beta+disc.
        return float(self.Gamma)/(beta/2+disc/2) if beta >= 0 else (-beta/2+disc/2)/float(self.A)

    def current_from_carrier(self, u):
        return self.p*u/(self.v*u+self.w)

    def current(self, x):
        return self.current_from_carrier(self.carrier(x))

    def invert(self, quota):
        positive(quota=quota)
        den=self.p-quota*self.v
        if den <= 0: raise ValueError('Quota reaches the kinetic asymptote')
        Y=quota*self.w/den
        delta=self.Gamma-self.A*Y*Y-self.B*Y
        if delta <= 0: raise ValueError('No finite positive NADPH floor')
        return Y, self.C*Y/delta

    def derivative(self, x):
        u=self.carrier(x); beta=self.B+self.C/x
        du=self.C*u/(x*x*(2*self.A*u+beta))
        return self.p*self.w*du/(self.v*u+self.w)**2

    def enclose(self, x, steps=100):
        """Fresh rational bisection, with no floating root used as evidence."""
        if not isinstance(x, F) or x <= 0: raise ValueError('Positive Fraction x required')
        lo,hi=F(0),F(1)
        while self.polynomial(x,hi) < 0: hi *= 2
        for _ in range(steps):
            mid=(lo+hi)/2
            if self.polynomial(x,mid) <= 0: lo=mid
            else: hi=mid
        return lo,hi


@dataclass(frozen=True)
class GlutathioneBranch:
    G: float; z0: float; E: float; a: float; b: float; c: float; k: float

    def __post_init__(self):
        positive(G=self.G,E=self.E,a=self.a,b=self.b,c=self.c,k=self.k)
        if not math.isfinite(float(self.z0)) or self.z0 < 0: raise ValueError('Invalid GSSG baseline')
        self.response  # enforce positive-pool condition at construction

    @property
    def rho(self): return self.a*(1/self.b+1/self.c)
    @property
    def S(self): return self.G-2*self.z0
    @property
    def response(self):
        return QuadraticResponse(1,self.rho-self.S,2*self.E*self.a/self.k,
            self.S*self.rho-self.E*self.a/self.c,self.E*self.a,1,self.rho)

    def reconstruct(self,x):
        positive(x=x); g=self.response.carrier(x); j=self.response.current_from_carrier(g)
        return dict(g=g,z=self.z0+j/(self.k*x),e0=j/self.a,e1=j/(self.b*g),e2=j/(self.c*g),jG=j)

    def redesign(self,x,quota):
        """Both algebraic candidates, with admissibility checked independently."""
        positive(x=x,quota=quota)
        B=self.S-2*quota/(self.k*x); C=quota/self.c; disc=B*B-4*C
        if B <= 0 or disc < 0: return []
        large=(B+math.sqrt(disc))/2
        roots=[large] if disc == 0 else [large,C/large]
        result=[]
        for g in roots:
            E=quota*(g+self.rho)/(self.a*g)
            try:
                branch=replace(self,E=E); state=branch.reconstruct(x)
                admissible=min(state.values()) > 0 and abs(state['jG']-quota) < 1e-8*quota
            except ValueError: admissible=False
            result.append(dict(g=g,E=E,admissible=admissible))
        return result


@dataclass(frozen=True)
class ThioredoxinBranch:
    T: float; z0: float; E: float; a: float; b: float; c: float; d: float; e: float; k: float

    def __post_init__(self):
        positive(T=self.T,E=self.E,a=self.a,c=self.c,d=self.d,e=self.e,k=self.k)
        if any(not math.isfinite(float(v)) or v < 0 for v in [self.z0,self.b]): raise ValueError('Invalid baseline/hyperoxidation rate')
        self.response

    @property
    def rho(self): return 1/self.a+1/self.d+self.b/(self.c*self.d)
    @property
    def S(self): return self.T-self.z0
    @property
    def response(self):
        v=self.rho*self.e
        return QuadraticResponse(v,1-self.S*v,self.E*self.e/self.k,self.S,self.E*self.e,v,1)

    def reconstruct(self,x):
        positive(x=x); y=self.response.carrier(x); j=self.response.current_from_carrier(y)
        return dict(y=y,zT=self.z0+j/(self.k*x),r=j/self.a,h=j/self.d,
            w=self.b*j/(self.c*self.d),v=j/(self.e*y),jT=j)


class DecreasingSource(Protocol):
    """Contract: strictly decreasing, positive below N, zero at N."""
    N: float
    def unit_current(self,x): ...


@dataclass(frozen=True)
class Regeneration:
    P: float; P0: float; K: float; V: float
    def __post_init__(self):
        positive(P=self.P,K=self.K,V=self.V)
        if not math.isfinite(float(self.P0)) or not 0 <= self.P0 < self.P: raise ValueError('Require 0 <= P0 < P')
    @property
    def N(self): return self.P-self.P0
    @property
    def D(self): return self.K+self.P
    def unit_current(self,x): return self.V*(self.N-x)/(self.D-x)
    def derivative(self,x): return -self.V*(self.D-self.N)/(self.D-x)**2


@dataclass(frozen=True)
class InhibitedSource:
    """Illustrative product inhibition of source; not inhibition of a branch."""
    base: Regeneration
    KI: float
    def __post_init__(self): positive(KI=self.KI)
    @property
    def N(self): return self.base.N
    def unit_current(self,x): return self.base.unit_current(x)/(1+x/self.KI)


@dataclass(frozen=True)
class ServiceDesign:
    responses: tuple[QuadraticResponse,...]
    source: DecreasingSource

    def __post_init__(self):
        if not self.responses: raise ValueError('At least one branch required')

    def demand(self,x): return sum(r.current(x) for r in self.responses)

    def equilibrium(self,scale):
        positive(scale=scale)
        return brentq(lambda x:scale*self.source.unit_current(x)-self.demand(x),0,float(self.source.N),xtol=1e-14)

    def decision(self,quotas):
        if len(quotas) != len(self.responses): raise ValueError('One positive quota per branch')
        for q in quotas: positive(quota=q)
        ceilings=[r.current(self.source.N) for r in self.responses]
        if any(q >= c for q,c in zip(quotas,ceilings)):
            return dict(repairable=False,ceilings=ceilings,reason='A quota reaches a finite-pool ceiling')
        inv=[r.invert(q) for r,q in zip(self.responses,quotas)]; L=max(x for _,x in inv)
        currents=[r.current(L) for r in self.responses]; R=self.source.unit_current(L)
        return dict(repairable=True,ceilings=ceilings,carrier_floors=[y for y,_ in inv],
            nadph_floors=[x for _,x in inv],L=L,currents=currents,
            minimum=sum(currents)/R,quota_only=sum(quotas)/R,
            overdelivery=[j-q for j,q in zip(currents,quotas)],
            isolated=[q/self.source.unit_current(x) for q,(_,x) in zip(quotas,inv)])

    @staticmethod
    def robust_command(scale,tolerance):
        if not 0 <= tolerance < 1: raise ValueError('Tolerance in [0,1) required')
        return scale/(1-tolerance)


@dataclass(frozen=True)
class MaintainedKinetics:
    """Literal eight-state ODE. Equilibrium response curves are NOT a transient closure."""
    gpx: GlutathioneBranch
    trx: ThioredoxinBranch
    source: Regeneration
    coordinates=('x','z','e1','e2','zT','h','w','v')

    @property
    def design(self): return ServiceDesign((self.gpx.response,self.trx.response),self.source)

    def full(self,u):
        x,z,e1,e2,zT,h,w,v=u
        return dict(zip(self.coordinates,u),g=self.gpx.G-2*z-e2,e0=self.gpx.E-e1-e2,
            y=self.trx.T-zT,r=self.trx.E-h-w-v,nadp=self.source.P-x)

    def state(self,x):
        s=dict(x=x,**self.gpx.reconstruct(x),**self.trx.reconstruct(x))
        return np.array([s[k] for k in self.coordinates])

    def currents(self,u,scale):
        s=self.full(u); G,T=self.gpx,self.trx
        return dict(source=scale*self.source.unit_current(s['x']),phi0=G.a*s['e0'],
            phi1=G.b*s['g']*s['e1'],phi2=G.c*s['g']*s['e2'],psiG=G.k*s['x']*(s['z']-G.z0),
            thetaR=T.a*s['r'],thetaW=T.b*s['h'],thetaH=T.c*s['w'],thetaV=T.d*s['h'],
            thetaY=T.e*s['y']*s['v'],psiT=T.k*s['x']*(s['zT']-T.z0))

    def rhs(self,t,u,scale):
        c=self.currents(u,scale)
        return np.array([c['source']-c['psiG']-c['psiT'],c['phi2']-c['psiG'],
            c['phi0']-c['phi1'],c['phi1']-c['phi2'],c['thetaY']-c['psiT'],
            c['thetaR']+c['thetaH']-c['thetaW']-c['thetaV'],c['thetaW']-c['thetaH'],c['thetaV']-c['thetaY']])

    def margins(self,u):
        s=self.full(u)
        return np.array([s[k] for k in self.coordinates]+[s[k] for k in ['g','e0','y','r']]+
            [self.source.N-s['x'],s['z']-self.gpx.z0,s['zT']-self.trx.z0])

    def jacobian(self,u,scale):
        return np.column_stack([self.rhs(0,np.asarray(u,dtype=complex)+1e-20j*np.eye(8)[i],scale).imag/1e-20 for i in range(8)])

    def integrate(self,initial,scale,times):
        positive(scale=scale)
        if self.margins(initial).min() < -1e-12: raise ValueError('Initial state outside conserved physical domain')
        times=np.asarray(times)
        if len(times)<2 or times[0]!=0 or not np.all(np.diff(times)>0): raise ValueError('Increasing times starting at zero required')
        sol=solve_ivp(lambda t,u:self.rhs(t,u,scale),(0,times[-1]),initial,t_eval=times,
            method='Radau',rtol=2e-9,atol=2e-11,jac=lambda t,u:self.jacobian(u,scale))
        if not sol.success: raise RuntimeError(sol.message)
        if min(self.margins(u).min() for u in sol.y.T) < -1e-7: raise RuntimeError('Numerical physical-domain violation')
        return sol.y.T


def nominal(exact=False):
    conv=lambda d:{k:F(str(v)) for k,v in d.items()} if exact else d.copy()
    # Literal paper parameters: deliberately independent of the editable block.
    G=GlutathioneBranch(**conv(dict(G=371.56,z0=1.78,E=50,a=.21,b=.04,c=10,k=3.2)))
    T=ThioredoxinBranch(**conv(dict(T=.505,z0=.075,E=19.096,a=.4,b=.00072,c=.003,d=15,e=2.1,k=20)))
    R=Regeneration(**conv(dict(P=30.3,P0=.3,K=57,V=375)))
    return MaintainedKinetics(G,T,R)


def exact_paper_audit():
    M=nominal(True); G,T,R=M.gpx,M.trx,M.source
    yg,xg=G.response.invert(F(10)); yt,xt=T.response.invert(F(4)); L=max(xg,xt)
    assert xg==F(3294375,138400918) and xt==F(460180,489387) and yt==F(5000,23009)
    assert T.response.polynomial(L,yt)==0 and T.response.current_from_carrier(yt)==4
    glo,ghi=F('361.1185'),F('361.1186')
    assert G.response.polynomial(L,glo)<0<G.response.polynomial(L,ghi)
    slo=(G.response.current_from_carrier(glo)+4)/R.unit_current(L)
    shi=(G.response.current_from_carrier(ghi)+4)/R.unit_current(L)
    assert F('.11371266')<slo<shi<F('.11371268')
    tight=G.response.enclose(L)
    sq=F(14)/R.unit_current(L)
    # Isolation succeeds at s=.1; joint failure is already forced at the floor.
    assert F(10)/R.unit_current(xg)<F('.1') and F(4)/R.unit_current(xt)<F('.1')<slo
    # Ceiling tests use the inverted polynomial, without rounded ceiling currents.
    assert G.response.polynomial(R.N,yg)<0 and T.response.polynomial(R.N,yt)<0
    # Larger redesign root in (361,362) uses less enzyme; check its pool condition.
    B=G.S-20/(G.k*L); C=F(10)/G.c
    poly=lambda u:u*u-B*u+C
    assert poly(F(361))<0<poly(F(362))
    emax=10*(361+G.rho)/(G.a*361)
    assert 0<emax<G.E and G.S*G.rho>emax*G.a/G.c
    assert sq>F('.1')  # no admitted GPx-only redesign can avoid this lower bound
    # Exact enclosure gives e0<=50<g, hence opposite Jacobian off-diagonal signs.
    assert glo>G.E
    return dict(XG=xg,XT=xt,carrier_T=yt,carrier_G_bracket=[glo,ghi],scale_bracket=[slo,shi],
        tighter_scale_bracket=[(G.response.current_from_carrier(u)+4)/R.unit_current(L) for u in tight],
        quota_only=sq,redesign_GSH_bracket=[361,362],redesign_enzyme_upper=emax,
        separate_success_joint_failure=True,ceiling_tests=True,opposite_sign_pair=True,
        evidence='Fresh rational checks under the paper equations; Lean is not rerun')


def declared_parameter_box():
    """Appendix B: rational enclosure of a continuum, not a corner optimization."""
    M=nominal(True); G,T,R=M.gpx,M.trx,M.source
    emax=G.E*F(101,100); emin=G.E*F(99,100)
    tmax=T.E*F(101,100); tmin=T.E*F(99,100)
    yglo=10*G.rho/(emax*G.a-10); yghi=10*G.rho/(emin*G.a-10)
    ytlo=4/(tmax*T.e-4*T.rho*T.e); ythi=4/(tmin*T.e-4*T.rho*T.e)
    dg=G.S*G.rho-emax*G.a/G.c-yghi*yghi-(G.rho-G.S)*yglo
    dt=T.S-T.rho*T.e*ythi**2-(1-T.S*T.rho*T.e)*ytlo
    assert dg>0 and dt>0 and G.rho-G.S<0 and 1-T.S*T.rho*T.e<0
    xg=(2*emax*G.a/(G.k*F(19,20)))*yghi/dg
    xt=(tmax*T.e/(T.k*F(19,20)))*ythi/dt
    target=F(3,2)
    assert max(xg,xt)<target<R.N
    upper=emax*G.a*G.S/(G.S+G.rho)+tmax*T.e*T.S/(1+T.rho*T.e*T.S)
    scale=upper/R.unit_current(target)
    assert scale==F(11621793767394811,92312851976765625)
    return dict(threshold_upper=[xg,xt],target=target,demand_upper=upper,scale=scale,
        command_with_5pct_tolerance=scale/F('.95'),scope='Fixed nominal box: E +/-1%, reductase k +/-5%; sufficient, not optimal')


class NonmonotoneExample:
    """Dimensionless analytic counterexample, kept outside ServiceDesign's contract."""
    @staticmethod
    def current(x): return x*(1.2-x)
    @staticmethod
    def equilibrium(scale):
        positive(scale=scale)
        return 2*scale/(1.2+scale+math.sqrt((1.2+scale)**2-4*scale))
    @staticmethod
    def scale_interval(): return ((12-3*math.sqrt(6))/10,(12+3*math.sqrt(6))/10)


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


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',default='outputs');args=parser.parse_args()
    out=Path(args.output);out.mkdir(parents=True,exist_ok=True)
    model=MaintainedKinetics(GlutathioneBranch(**GPX),ThioredoxinBranch(**TRX),Regeneration(**SOURCE))
    design=model.design;decision=design.decision(QUOTAS)
    result=dict(inputs=dict(GPX=GPX,TRX=TRX,SOURCE=SOURCE,quotas=QUOTAS),configured=decision,
        fixed_paper_exact=exact_paper_audit(),fixed_parameter_box=declared_parameter_box())
    # An edited impossible quota is a result, not a failed root solve.
    dump=lambda name,obj:(out/name).write_text(json.dumps(obj,indent=2,default=lambda v:str(v) if isinstance(v,F) else float(v))+'\n',encoding='utf-8')
    if not decision['repairable']:
        dump('results.json',result);print('No finite regeneration repair: a quota reaches its branch ceiling.');return
    L=decision['L'];sm=decision['minimum'];eq=model.state(L)
    result['boundary_state']=model.full(eq)
    result['boundary_residual']=float(np.max(np.abs(model.rhs(0,eq,sm))))
    result['source_capacity']=sm*model.source.V
    result['boundary_accounting']=dict(NADPH=sum(decision['currents']),
        peroxide=decision['currents'][0]+decision['currents'][1]*(1+model.trx.b/model.trx.d),
        external_Prx_repair=model.trx.b/model.trx.d*decision['currents'][1])
    result['constant_tolerance_command']=design.robust_command(sm,SOURCE_TOLERANCE)
    result['inhibited_source_example']=ServiceDesign(design.responses,InhibitedSource(model.source,5.)).decision(QUOTAS)
    isolated=[]
    for response in design.responses:
        x=ServiceDesign((response,),model.source).equilibrium(REFERENCE_SCALE)
        isolated.append(dict(x=x,current=response.current(x)))
    xjoint=design.equilibrium(REFERENCE_SCALE)
    result['reference_comparison']=dict(isolated=isolated,joint_x=xjoint,
        joint_currents=[r.current(xjoint) for r in design.responses])
    redesign=[]
    if decision['nadph_floors'][0] <= decision['nadph_floors'][1]:
        for candidate in model.gpx.redesign(L,QUOTAS[0]):
            if candidate['admissible']:
                other=replace(model,gpx=replace(model.gpx,E=candidate['E']))
                candidate['decision']=other.design.decision(QUOTAS)
                candidate['full_residual']=float(np.max(np.abs(other.rhs(0,other.state(L),decision['quota_only']))))
            redesign.append(candidate)
    result['redesign_candidates']=redesign
    traces=[];dynamics=[]
    times=np.unique(np.r_[np.linspace(0,min(5.,HORIZON),251),np.linspace(min(5.,HORIZON),HORIZON,601)])
    for scale,label in [(REFERENCE_SCALE,'reference'),(sm,'boundary'),(SLACK_SCALE,'slack')]:
        equilibrium=model.state(design.equilibrium(scale));jac=model.jacobian(equilibrium,scale)
        eigen=np.linalg.eigvals(jac)
        for kind in ['carrier oxidation','enzyme allocation']:
            initial=equilibrium.copy()
            if kind=='carrier oxidation': initial[[1,4]]+=CARRIER_PERTURBATION
            else: initial[[2,5]]+=ENZYME_PERTURBATION
            states=model.integrate(initial,scale,times);services=[];dist=[]
            for t,u in zip(times,states):
                c=model.currents(u,scale);err=float(max(abs(u-equilibrium)));dist.append(err)
                services.append([c['phi0'],c['thetaY']])
                traces.append([label,kind,t,*u,c['phi0'],c['thetaY'],c['psiG']+c['psiT'],
                    c['phi0']+c['thetaR']+c['thetaW'],c['thetaH'],err])
            services=np.array(services)
            dynamics.append(dict(setting=label,scale=scale,perturbation=kind,
                largest_eigen_real=float(max(eigen.real)),sampled_service_min=services.min(axis=0).tolist(),
                sampled_domain_min=float(min(model.margins(u).min() for u in states)),
                final_distance=dist[-1],equilibrium_residual=float(max(abs(model.rhs(0,equilibrium,scale))))))
    result['numerical_dynamics']=dynamics
    J=model.jacobian(eq,sm);result['orthant_obstruction_pair']=[J[2,3],J[3,2]]
    # Sweep from a low quota to near its finite-source ceiling, resolving the narrow GPx binding interval.
    sweep=[]
    for branch in [0,1]:
        switch=decision['currents'][branch];ceil=decision['ceilings'][branch]
        grid=np.unique(np.r_[np.linspace(.8*QUOTAS[branch],min(switch,ceil*.999999),90),
            np.linspace(switch,ceil-(ceil-switch)*.01,80)])
        for q in grid:
            qs=list(QUOTAS);qs[branch]=q;d=design.decision(qs)
            if d['repairable']: sweep.append([branch,q,d['L'],d['minimum'],d['quota_only']])
    nonmono=[[s,NonmonotoneExample.equilibrium(s),NonmonotoneExample.current(NonmonotoneExample.equilibrium(s))] for s in np.linspace(.05,3,180)]
    result['nonmonotone_dimensionless']=dict(success_scale_interval=NonmonotoneExample.scale_interval(),unique_equilibrium=True,
        explanation='h derivative = ((1-x)^2+1/5)/(1-x)^2 >0; more source eventually breaks the quota')
    write_table(out/'transients.csv',['setting','perturbation','time_s',*model.coordinates,'GPx_peroxide','Trx_oxidation','NADPH_demand','total_peroxide','Prx_repair','state_distance'],traces)
    write_table(out/'quota_sweep.csv',['branch','quota_uM_per_s','NADPH_floor_uM','minimum_scale','quota_only_scale'],sweep)
    write_table(out/'nonmonotone.csv',['scale','x','current'],nonmono)
    write_table(out/'boundary_species.csv',['species','uM'],result['boundary_state'].items())
    dump('results.json',result)
    lines=[f'Both separate currents at s={REFERENCE_SCALE}: '+', '.join(format(v['current'],'.9f') for v in isolated),
        f'Joint currents: '+', '.join(f'{r.current(xjoint):.9f}' for r in design.responses),
        f'Minimum scale {sm:.12f}; quota-only {decision["quota_only"]:.12f}; capacity {sm*model.source.V:.9f} uM/s.',
        f'Boundary NADPH demand {sum(decision["currents"]):.9f} uM/s; full ODE residual {result["boundary_residual"]:.3g}.',
        'Exact paper audit and parameter-box bounds passed. Dynamic stability and transient service are numerical only.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
    plot(out,model,decision,traces,sweep,nonmono)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    dump('run_metadata.json',dict(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'}))


def plot(out,model,d,traces,sweep,nonmono):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    xs=np.linspace(0,min(model.source.N,3*d['L']),400)
    for i,label in enumerate(['GPx alone','Prx/Trx alone']):axs[0].plot(xs,[model.design.responses[i].current(x) for x in xs],label=label)
    axs[0].plot(xs,[model.design.demand(x) for x in xs],label='Joint demand')
    for s,label in [(REFERENCE_SCALE,'Reference source'),(d['minimum'],'Minimum source')]:axs[0].plot(xs,[s*model.source.unit_current(x) for x in xs],'--',label=label)
    axs[0].axvline(d['L'],color='k',ls=':');axs[0].set(xlabel='NADPH (μM)',ylabel='Current (μM/s)',title='NADPH regeneration and branch consumption');axs[0].legend(fontsize=7)
    data=np.array(sweep);a=data[data[:,0]==0]
    # Use NADPH floor to expose the tiny quota interval without collapsing the switch and ceiling.
    axs[1].plot(a[:,1],a[:,3],label='Fixed kinetics');axs[1].plot(a[:,1],a[:,4],'--',label='Quota-only');axs[1].axvline(d['currents'][0],ls=':',color='k',label='Binding switch')
    axs[1].set(xlabel='GPx quota (μM/s)',ylabel='Required source scale',ylim=(.75*d['minimum'],1.5*d['minimum']),title='Required regeneration scale versus GPx\ntarget');axs[1].legend(fontsize=7,loc='lower left')
    inset=axs[1].inset_axes([.12,.53,.42,.38]);tail=a[a[:,1]>=d['currents'][0]];inset.plot(tail[:,1],tail[:,3]);inset.axvline(d['ceilings'][0],ls=':',color='r');inset.tick_params(labelsize=6);inset.set_xticks([d['currents'][0],d['ceilings'][0]],[format(d['currents'][0],'.6f'),format(d['ceilings'][0],'.6f')]);inset.set_yscale('log');inset.set_title('Switch → ceiling',fontsize=7)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'shared_source.png',dpi=180);fig.savefig(out/'shared_source.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    for setting,style in [('boundary','--'),('slack','-')]:
        a=np.array([row[2:] for row in traces if row[0]==setting and row[1]=='carrier oxidation'],float)
        axs[0].plot(a[:,0],a[:,10],style,label=f'{setting}: Trx oxidation')
    axs[0].axhline(QUOTAS[1],color='k',ls=':',label='Trx quota');axs[0].set(xlim=(0,min(30,HORIZON)),xlabel='Time (s)',ylabel='Current (μM/s)',title='Thioredoxin turnover after carrier\noxidation');axs[0].legend(fontsize=7)
    a=np.array(nonmono);lo,hi=NonmonotoneExample.scale_interval();axs[1].plot(a[:,0],a[:,2]);axs[1].axhline(.3,color='k',ls=':');axs[1].axvspan(lo,hi,color='green',alpha=.12,label='Successful scales');axs[1].set(xlabel='Source scale (dimensionless)',ylabel='Service (dimensionless)',title='Service versus regeneration in the\nnon-monotone model');axs[1].legend(fontsize=7)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'operating_limits.png',dpi=180);fig.savefig(out/'operating_limits.svg');plt.close(fig)


if __name__=='__main__': main()