"""Enzyme-conversion models, exact child reductions, core spectra and saturating kinetics."""
# EDITABLE STUDY INPUTS. Concentrations and time are nondimensional.
SITES = 3
CENSUS_SITES = (2, 3)
MAX_CENSUS_CHILDREN = 250000
BACKGROUND_REACTIVITY = '1/100'
OWNER_REACTIVITY = '1'
PERTURBATION = 1e-5
TIME_HORIZON = 40.
TIME_SAMPLES = 801
FLOWER_RATES = ('1', '2', '3', '1', '2', '1', '2')  # n=3: long cycle, then leaf
FLOWER_LOSSES = ('1/10',)*7
ASYMPTOTIC_SITES = (2,3,5,10,25,100,1000,10000)

from dataclasses import dataclass
from fractions import Fraction as Q
from itertools import combinations,product
from collections import Counter
import argparse
import csv
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 scipy.optimize import brentq
from scipy.special import lambertw,logsumexp

MANUSCRIPT_SHA256='370eea21e788907581cc3c677783a0bd1581afe53307babfdc2325e77120f2aa'


@dataclass(frozen=True)
class ConversionArm:
    substrate:str
    product:str
    enzyme:str
    intermediate:str
    name:str


class ConversionSystem:
    def __init__(self,substrates,enzymes,arms):
        self.substrates=tuple(substrates);self.enzymes=tuple(enzymes);self.arms=tuple(arms)
        arms=self.arms
        intermediates=tuple(a.intermediate for a in arms);self.names=self.substrates+self.enzymes+intermediates
        if len(set(self.names))!=len(self.names) or len({a.name for a in arms})!=len(arms):raise ValueError('Species roles must be disjoint; intermediates and arm names unique.')
        self.index={s:i for i,s in enumerate(self.names)};reactions=[];Y=[];P=[]
        for arm in arms:
            if arm.substrate not in self.substrates or arm.product not in self.substrates or arm.enzyme not in self.enzymes:raise ValueError('Arm role mismatch.')
            for suffix,inputs,outputs in [('b',(arm.substrate,arm.enzyme),(arm.intermediate,)),('u',(arm.intermediate,),(arm.substrate,arm.enzyme)),('v',(arm.intermediate,),(arm.product,arm.enzyme))]:
                reactions.append(arm.name+suffix);a=[0]*len(self.names);b=a.copy()
                for s in inputs:a[self.index[s]]+=1
                for s in outputs:b[self.index[s]]+=1
                Y.append(a);P.append(b)
        self.reactions=tuple(reactions);self.rindex={r:j for j,r in enumerate(reactions)};self.Y=np.array(Y,dtype=int).T;self.P=np.array(P,dtype=int).T;self.S=self.P-self.Y
        self.intermediate_rows=frozenset(self.index[s] for s in intermediates);self.enzyme_rows=frozenset(self.index[s] for s in self.enzymes)
    @classmethod
    def futile(cls,n):
        if type(n) is not int or not 1<=n<=30:raise ValueError('Explicit futile-cycle model requires 1 <= n <= 30.')
        arms=[]
        for i in range(n):
            arms += [ConversionArm(f'S{i}',f'S{i+1}','K',f'C{i}',f'K{i}'),ConversionArm(f'S{i+1}',f'S{i}','F',f'D{i+1}',f'F{i+1}')]
        return cls([f'S{i}' for i in range(n+1)],['K','F'],arms)
    def child(self,owners,reactions):
        if len(owners)!=len(reactions) or len(set(owners))!=len(owners) or len(set(reactions))!=len(reactions):raise ValueError('Distinct owners and injective reaction assignment required.')
        try:rows=tuple(self.index[s] for s in owners);cols=tuple(self.rindex[r] for r in reactions)
        except KeyError as e:raise ValueError('Unknown owner or reaction.') from e
        if not all(self.Y[i,j]>0 for i,j in zip(rows,cols)):raise ValueError('Every owner must be a reactant of its assigned reaction.')
        return ChildSelection(self,rows,cols)
    def selections(self,limit=MAX_CENSUS_CHILDREN):
        options=[(None,)+tuple(int(j) for j in np.flatnonzero(self.Y[i])) for i in range(len(self.names))];count=0
        def visit(i,rows,cols,used):
            nonlocal count
            if i==len(options):
                count+=1
                if count>limit:raise ValueError('Child census budget exceeded; partial counts are not a complete census.')
                yield ChildSelection(self,tuple(rows),tuple(cols));return
            yield from visit(i+1,rows,cols,used)
            for j in options[i][1:]:
                if j not in used:yield from visit(i+1,rows+[i],cols+[j],used|{j})
        yield from visit(0,[],[],set())
    def conservation(self):
        pools=[{e}|{a.intermediate for a in self.arms if a.enzyme==e} for e in self.enzymes]
        pools.append(set(self.substrates)|{a.intermediate for a in self.arms})
        return np.array([[int(s in pool) for s in self.names] for pool in pools],int)


def incidence_determinant(H):
    """Cofactor induction; each column has at most one entry of each sign."""
    H=np.asarray(H,int);n=len(H);sign=1;steps=[]
    if H.shape!=(n,n):raise ValueError('Square incidence matrix required.')
    if any(np.sum(H[:,j]==1)>1 or np.sum(H[:,j]==-1)>1 or np.any(abs(H[:,j])>1) for j in range(n)):raise ArithmeticError('Not a signed-incidence matrix.')
    while n:
        chosen=next((j for j in range(n) if np.count_nonzero(H[:,j])<=1),None)
        if chosen is None:
            if np.any(H.sum(axis=0)):raise ArithmeticError('Incidence zero-sum argument failed.')
            return 0,steps+[{'reason':'all columns sum to zero','dimension':n}]
        rows=np.flatnonzero(H[:,chosen])
        if len(rows)==0:return 0,steps+[{'reason':'zero column','column':chosen,'dimension':n}]
        i=int(rows[0]);value=int(H[i,chosen]);sign*=(-1)**(i+chosen)*value;steps.append({'row':i,'column':chosen,'entry':value,'dimension':n})
        H=np.delete(np.delete(H,i,axis=0),chosen,axis=1);n-=1
    return sign,steps


@dataclass(frozen=True)
class ChildSelection:
    system:ConversionSystem
    rows:tuple
    columns:tuple
    @property
    def matrix(self):return self.system.S[np.ix_(self.rows,self.columns)]
    def reduction(self,details=True):
        system=self.system;free=[i for i,r in enumerate(self.rows) if r not in system.intermediate_rows];bound=[i for i,r in enumerate(self.rows) if r in system.intermediate_rows];order=free+bound
        A=self.matrix[np.ix_(order,order)];p=len(free);q=len(bound)
        if not np.array_equal(A[p:,p:],-np.eye(q,dtype=int)):raise ArithmeticError('Selected intermediate block is not -I.')
        reduced=A[:p,:p]+A[:p,p:]@A[p:,:p];enzyme_count=0;H=reduced.copy()
        for i,original in enumerate(free):
            if self.rows[original] in system.enzyme_rows:H[i]*=-1;enzyme_count+=1
        detH,steps=incidence_determinant(H);value=(-1)**(q+enzyme_count)*detH
        if not details:return value
        return {'owners':[system.names[i] for i in self.rows],'reactions':[system.reactions[j] for j in self.columns],'determinant':value,
            'selected_intermediates':q,'enzyme_row_signs':enzyme_count,'reduced_matrix':reduced.tolist(),'signed_incidence_matrix':H.tolist(),'cofactor_certificate':steps}


def negative_child(system):return system.child(['S1','S3','K','F','C2','D1'],['K1b','F3b','K2b','F1b','K2v','F1v'])
def positive_child(system,n):return system.child(['F']+[s for i in range(1,n) for s in (f'S{i}',f'C{i}')]+[f'S{n}','D1'],['F1b']+[r for i in range(1,n) for r in (f'K{i}b',f'K{i}v')]+[f'F{n}b','F1v'])


def exact_routh(coefficients):
    coefficients=list(map(sp.Rational,coefficients));n=len(coefficients)-1;width=(n+2)//2;table=[[sp.Rational(0)]*width for _ in range(n+1)]
    for j,c in enumerate(coefficients):table[j%2][j//2]=c
    for i in range(2,n+1):
        if table[i-1][0]==0:raise ValueError('Degenerate Routh pivot requires a separate analysis; result unresolved.')
        for j in range(width-1):table[i][j]=(table[i-1][0]*table[i-2][j+1]-table[i-2][0]*table[i-1][j+1])/table[i-1][0]
    first=[r[0] for r in table]
    if any(v==0 for v in first):raise ValueError('Zero Routh first-column entry; result unresolved.')
    signs=[int(sp.sign(v)) for v in first]
    return {'first_column':list(map(str,first)),'signs':signs,'right_half_plane_roots':sum(a!=b for a,b in zip(signs,signs[1:]))}


def cycle_block_certificate(M):
    """SCCs of this restriction are isolated vertices or unit positive cycles.

    After arbitrary positive column scaling each cycle block has zero column
    sums and nonnegative off-diagonal entries, hence is nonunstable.
    """
    M=np.asarray(M,int);n=len(M);reach=[{i}|{j for j in range(n) if j!=i and M[j,i]} for i in range(n)]
    for k in range(n):
        for i in range(n):
            if k in reach[i]:reach[i]|=reach[k]
    remaining=set(range(n));components=[]
    while remaining:
        i=min(remaining);component=sorted(j for j in remaining if j in reach[i] and i in reach[j]);remaining-=set(component);B=M[np.ix_(component,component)]
        if np.any(np.diag(B)!=-1):raise ArithmeticError('Unexpected diagonal.')
        if len(component)>1:
            off=B+np.eye(len(B),dtype=int)
            if np.any(off<0) or not np.all(off.sum(axis=0)==1) or not np.all(off.sum(axis=1)==1):raise ArithmeticError('SCC is not a unit positive cycle.')
        components.append(component)
    return components


def negative_certificates(system):
    A=sp.Matrix(negative_child(system).matrix);B=A[:5,:5];D=sp.diag(1,1,2,1,2);z=sp.Symbol('z')
    energy=sp.Matrix([[1443,770,-1334,-1169,-1041],[770,1133,-371,-1024,153],[-1334,-371,1627,895,1518],[-1169,-1024,895,1278,456],[-1041,153,1518,456,1780]])
    minors=[energy[:k,:k].det() for k in range(1,6)]
    if B.T*energy+energy*B!=-218*sp.eye(5) or any(v<=0 for v in minors):raise ArithmeticError('Exceptional restriction Lyapunov certificate failed.')
    for mask in range(1,63):
        rows=[i for i in range(6) if mask>>i&1]
        if rows==list(range(5)):continue
        cycle_block_certificate(np.array(A.extract(rows,rows),int))
    for mask in range(1,31):cycle_block_certificate(np.array(B.extract([i for i in range(5) if mask>>i&1],[i for i in range(5) if mask>>i&1]),int))
    p=(B*D).charpoly(z).as_expr();value=sp.simplify(sp.re(sp.diff(p,z).subs(z,3*sp.I/5)/p.subs(z,3*sp.I/5)))
    if value!=-sp.Rational(24183425,1544506):raise ArithmeticError('Logarithmic derivative certificate failed.')
    return {'A_minus':{'matrix':A.tolist(),'determinant':str(A.det()),'routh':exact_routh(A.charpoly(z).all_coeffs()),'proper_restrictions_checked':62},
        'B':{'matrix':B.tolist(),'determinant':str(B.det()),'routh':exact_routh(B.charpoly(z).all_coeffs()),'lyapunov_matrix':energy.tolist(),'lyapunov_leading_minors':list(map(str,minors))},
        'B_scaled':{'scaling':[1,1,2,1,2],'routh':exact_routh((B*D).charpoly(z).all_coeffs()),'log_derivative_real_at_3i_over_5':str(value),'proper_D_nonunstable_restrictions_checked':30}}


class FlowerSpectrum:
    def __init__(self,n,rates=None,losses=None):
        if type(n) is not int or not 2<=n<=10000:raise ValueError('Flower size outside scalar-solver budget.')
        self.n=n;self.rates=tuple(Q(v) for v in (rates if rates is not None else [1]*(2*n+1)));self.losses=tuple(Q(v) for v in (losses if losses is not None else [0]*(2*n+1)))
        if len(self.rates)!=2*n+1 or len(self.losses)!=2*n+1 or min(self.rates)<=0 or min(self.losses)<0:raise ValueError('Positive column rates and nonnegative independent losses required.')
    def matrix(self):
        if self.n>30:raise ValueError('Use the scalar equation instead of a large dense matrix.')
        n=2*self.n;M=-np.eye(n+1)
        for i in range(n-1):M[i+1,i]=1
        M[0,n-1]=1;M[n,0]=1;M[0,n]=1
        return M*np.array(self.rates,float)-np.diag(np.array(self.losses,float))
    def gain_at_zero(self):
        f=[d/(d+eta) for d,eta in zip(self.rates,self.losses)];return math.prod(f[:-1])+f[0]*f[-1]
    def abscissa(self):
        d=np.array(self.rates,float);a=d+np.array(self.losses,float);logs=np.log(d)
        def loggain(beta):
            denominators=beta+a
            if np.any(denominators<=0):return math.inf
            f=logs-np.log(denominators);return float(logsumexp([np.sum(f[:-1]),f[0]+f[-1]]))
        lower=np.nextafter(-float(min(a)),math.inf);upper=max(1.,float(max(d)))
        while loggain(upper)>0:upper*=2
        return brentq(loggain,lower,upper,xtol=1e-14)
    def eigenvector(self,beta):
        d=np.array(self.rates,float);a=d+np.array(self.losses,float);u=np.ones(len(d))
        for i in range(1,len(d)-1):u[i]=d[i]*u[i-1]/(beta+a[i])
        u[-1]=d[-1]/(beta+a[-1]);v=u/d;return v/max(v)
    @staticmethod
    def unit_alpha(n):return brentq(lambda a:math.exp(-2*n*math.log1p(a))+math.exp(-2*math.log1p(a))-1,0.,1.,xtol=1e-15)
    @staticmethod
    def lambert_bracket(n):
        m=n-1;r=float(lambertw(m).real);low=r/m;high=r/(m-r)
        return low/(math.sqrt(1+low)+1),high/(math.sqrt(1+high)+1),r/(2*m)


class SaturatingKinetics:
    def __init__(self,system,reference,flux,reactivity):
        self.system=system;self.reference=sp.Matrix(reference);self.flux=sp.Matrix(flux);self.R=sp.Matrix(reactivity);S=sp.Matrix(system.S)
        if self.reference.shape!=(len(system.names),1) or self.flux.shape!=(len(system.reactions),1) or self.R.shape!=(len(system.reactions),len(system.names)):raise ValueError('Kinetic dimensions mismatch.')
        if any(v<=0 for v in list(self.reference)+list(self.flux)) or S*self.flux!=sp.zeros(len(system.names),1):raise ValueError('Positive reference and balanced positive flux required.')
        self.theta={}
        for j in range(len(system.reactions)):
            for i in range(len(system.names)):
                if system.Y[i,j]:
                    th=self.R[j,i]*self.reference[i]/self.flux[j]
                    if not 0<th<1:raise ValueError('Saturating factors require 0 < R[j,i]*xstar[i]/v[j] < 1.')
                    self.theta[j,i]=th
                elif self.R[j,i]!=0:raise ValueError('Reactivity outside reactant support.')
        self.S=np.array(S,float);self.xstar=np.array(self.reference,float).ravel();self.v=np.array(self.flux,float).ravel()
    @classmethod
    def paper(cls,system,background=BACKGROUND_REACTIVITY,owner=OWNER_REACTIVITY):
        R=sp.Matrix(system.Y.T)*sp.Rational(background);child=negative_child(system)
        for i,j in zip(child.rows,child.columns):R[j,i]=sp.Rational(owner)
        flux=[4 if r.endswith('b') else 2 for r in system.reactions]
        return cls(system,[1]*len(system.names),flux,R)
    @property
    def exact_jacobian(self):return sp.Matrix(self.system.S)*self.R
    def rates_and_jacobian(self,eta):
        eta=np.asarray(eta,float);relative=eta/self.xstar
        if np.any(relative<=-1):raise ArithmeticError('Nonpositive concentration; no clipping.')
        ratios=np.ones(len(self.v));delta=np.zeros(len(self.v));derivatives=np.zeros(self.R.shape)
        factors={}
        for (j,i),theta in self.theta.items():
            th=float(theta);den=1+(1-th)*relative[i];change=th*relative[i]/den;factors[j,i]=1+change;delta[j]=delta[j]*(1+change)+change;ratios[j]*=1+change
        for (j,i),theta in self.theta.items():
            th=float(theta);den=1+(1-th)*relative[i]
            derivatives[j,i]=self.v[j]*ratios[j]/factors[j,i]*th/(den*den*self.xstar[i])
        return self.v*ratios,self.S@(self.v*delta),self.S@derivatives
    def integrate(self,initial,times,method='Radau'):
        times=np.asarray(times,float)
        if times[0]!=0 or np.any(np.diff(times)<=0):raise ValueError('Increasing times starting at zero required.')
        solution=solve_ivp(lambda t,e:self.rates_and_jacobian(e)[1],(0,times[-1]),initial,method=method,
            jac=lambda t,e:self.rates_and_jacobian(e)[2],t_eval=times,rtol=2e-10,atol=2e-13)
        if not solution.success or np.any(solution.y.T+self.xstar<=0):raise ArithmeticError('Kinetic integration failed.')
        return solution.y.T+self.xstar


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
    if SITES<3:raise ValueError('The full negative-core kinetic illustration requires at least three sites.')
    system=ConversionSystem.futile(SITES);negative=negative_certificates(system);positive=positive_child(system,SITES);reduction=positive.reduction();size=2*SITES+1;T=np.diag([-1]+[1]*(size-2)+[-1])
    flower=FlowerSpectrum(SITES)
    if not np.array_equal(T@positive.matrix@T,flower.matrix()):raise ArithmeticError('Source-to-flower identity failed.')
    proper_count=0
    if size<=11:
        for mask in range(1,(1<<size)-1):
            rows=[i for i in range(size) if mask>>i&1];cycle_block_certificate((T@positive.matrix@T)[np.ix_(rows,rows)]);proper_count+=1
    census=[]
    for n in CENSUS_SITES:
        network=ConversionSystem.futile(n);counts=Counter(selection.reduction(False) for selection in network.selections())
        census.append({'n':n,'children_including_empty':sum(counts.values()),'determinants':{str(k):counts[k] for k in (-1,0,1)}})
    selected_flower=FlowerSpectrum(3,FLOWER_RATES,FLOWER_LOSSES);beta=selected_flower.abscissa();gain=selected_flower.gain_at_zero();v=selected_flower.eigenvector(beta)
    if np.max(abs(selected_flower.matrix()@v-beta*v))>1e-10:raise ArithmeticError('Return-equation eigenvector failed.')
    zero_loss=FlowerSpectrum(3,FLOWER_RATES);alpha=zero_loss.abscissa();loss_sweep=[(delta,alpha-delta) for delta in np.linspace(0,2*alpha,101)]
    asymptotic=[]
    for n in ASYMPTOTIC_SITES:
        a=FlowerSpectrum.unit_alpha(n);lower,upper,leading=FlowerSpectrum.lambert_bracket(n);asymptotic.append((n,2*n+1,lower,a,upper,leading,math.log(n)/(2*n)))
        if not lower<a<upper:raise ArithmeticError('Lambert bracket failed numerically.')
    model=SaturatingKinetics.paper(system);J=model.exact_jacobian;L=sp.Matrix(system.conservation());S=sp.Matrix(system.S)
    if L*S!=sp.zeros(3,len(system.reactions)) or S.rank()!=len(system.names)-3:raise ArithmeticError('Conservation rank failed.')
    coefficients=(100*J).charpoly().all_coeffs();zeros=0
    while coefficients[-1]==0:zeros+=1;coefficients.pop()
    full_routh=exact_routh(coefficients);p10=sp.Poly((sp.Symbol('z')+1)**21-(sp.Symbol('z')+1)**19-(sp.Symbol('z')+1),sp.Symbol('z'))
    eig,vectors=np.linalg.eig(np.array(J,float));index=max((i for i,z in enumerate(eig) if z.imag>1e-8),key=lambda i:eig[i].real);mode=vectors[:,index];mode/=np.max(abs(mode))
    initial=PERTURBATION*mode.real;Ln=np.array(L,float);initial-=Ln.T@np.linalg.solve(Ln@Ln.T,Ln@initial)
    times=np.linspace(0,TIME_HORIZON,TIME_SAMPLES);states=model.integrate(initial,times);independent=model.integrate(initial,times,'BDF');linear=PERTURBATION*np.real(np.exp(eig[index]*times[:,None])*mode)
    totals=states@Ln.T;conservation_error=float(np.max(abs(totals-totals[0])));linear_error=float(np.max(abs(states-1-linear)))
    spectra={'A_minus':np.linalg.eigvals(np.array(negative_child(system).matrix,float)),'B_scaled':np.linalg.eigvals(np.array(negative_child(system).matrix[:5,:5]@np.diag([1,1,2,1,2]),float)),'full_J':eig}
    result={'negative_certificates':negative,'positive_child':{'dimension':size,'reduction':reduction,'proper_restrictions_checked':proper_count,'scope':'All-scaling minimality follows from the flower cycle structure; finite restrictions replayed within the stated size guard.'},'census':census,
        'n10_routh':exact_routh(p10.all_coeffs()),'selected_flower':{'n':3,'rates':FLOWER_RATES,'independent_losses':FLOWER_LOSSES,'exact_gain_at_zero':str(gain),'spectral_abscissa':beta,'zero_loss_abscissa':alpha},
        'full_kinetics':{'species':system.names,'reactions':system.reactions,'stoichiometry':system.S.tolist(),'reference_flux':list(map(str,model.flux)),'reactivity':[[str(x) for x in row] for row in model.R.tolist()],
            'conservation_rows':L.tolist(),'reference_totals':list(map(str,L*model.reference)),'stoichiometric_rank':S.rank(),'exact_conservation_zeros':zeros,'nonzero_scaled_characteristic_coefficients':list(map(str,coefficients)),'routh':full_routh,
            'conservation_max_error':conservation_error,'independent_solver_difference':float(np.max(abs(states-independent))),'linearization_max_error':linear_error,'minimum_concentration':float(np.min(states))},
        'spectra':{name:[[float(z.real),float(z.imag)] for z in values] for name,values in spectra.items()},
        'scope':'Exact determinants and root counts plus numerical rate/trajectory illustrations. Full kinetics are saturating, not mass action; no Hopf bifurcation or periodic orbit is established. Lean not rerun.'}
    def write_json(name,data):(out/name).write_text(json.dumps(data,indent=2,default=lambda x:int(x) if isinstance(x,sp.Integer) else str(x))+'\n')
    write_json('results.json',result)
    def table(name,header,rows):
        with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    table('trajectories.csv',['time',*system.names,*[s+'_linear_departure' for s in system.names]],np.column_stack([times,states,linear]))
    table('unit_rate_scaling.csv',['sites','core_dimension','Lambert_lower','alpha','Lambert_upper','Lambert_leading','log_equivalent'],asymptotic)
    table('uniform_loss.csv',['independent_uniform_loss','spectral_abscissa'],loss_sweep)
    table('child_census.csv',['sites','children','det_minus_one','det_zero','det_plus_one'],[(r['n'],r['children_including_empty'],r['determinants']['-1'],r['determinants']['0'],r['determinants']['1']) for r in census])
    lines=[f'Determinant censuses: {census}.', 'Six-species negative core: two unstable roots; five-species restriction: stable at unit rates, two unstable roots after scaling (1,1,2,1,2).',
        f'Positive child dimension {size}; n=10 has three unstable roots, despite a unique positive real root.',
        f'Selected flower gain at zero={gain}; spectral abscissa={beta:.9g}; uniform-loss threshold without other losses={alpha:.9g}.',
        f'Full saturating model: {zeros} exact conservation zeros and {full_routh["right_half_plane_roots"]} unstable roots; conservation error={conservation_error:.3g}.',
        'Selected matrix rates/losses are not automatically realizable interventions on the full chemical model. No oscillatory orbit claimed.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,3,figsize=(11,3.8),layout='constrained')
    for ax,(name,values) in zip(axs,spectra.items()):
        values=np.array(values);ax.scatter(values.real,values.imag,c=['#bd5a24' if z.real>1e-8 else '#417d8c' for z in values]);ax.axvline(0,color='gray',lw=.7);ax.set(xlabel='Real part',ylabel='Imaginary part',title=name.replace('_',' '));ax.grid(alpha=.2)
    fig.savefig(out/'spectra.png',dpi=180);fig.savefig(out/'spectra.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained');data=np.array(asymptotic)
    axs[0].loglog(data[:,0],data[:,3],'o-',label='Exact scalar equation (numerical root)');axs[0].loglog(data[:,0],data[:,5],'--',label='Leading Lambert term');axs[0].loglog(data[:,0],data[:,6],':',label='log(n)/(2n)')
    axs[0].set(xlabel='Phosphorylation sites n',ylabel='Unit-rate spectral abscissa',title='Dominant growth rate versus\nphosphorylation-site count');axs[0].legend(fontsize=7)
    axs[1].plot(*np.array(loss_sweep).T);axs[1].axhline(0,color='gray',lw=.7);axs[1].axvline(alpha,color='#bd5a24',ls='--');axs[1].set(xlabel='Independent uniform loss',ylabel='Selected-matrix spectral abscissa',title='Selected-matrix stability versus uniform\nloss')
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'rates_and_losses.png',dpi=180);fig.savefig(out/'rates_and_losses.svg');plt.close(fig)
    fig,ax=plt.subplots(figsize=(8,4.2),layout='constrained')
    for name in ('S1','S3','K','F'):
        i=system.index[name];line,=ax.plot(times,(states[:,i]-1)/PERTURBATION,label=name);ax.plot(times,linear[:,i]/PERTURBATION,'--',color=line.get_color(),lw=.8)
    ax.set(xlabel='Time (nondimensional)',ylabel='Departure / initial perturbation scale',title='Full saturating network and its unstable linear mode');ax.legend();ax.grid(alpha=.2)
    fig.savefig(out/'kinetic_departure.png',dpi=180);fig.savefig(out/'kinetic_departure.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    write_json('run_metadata.json',{'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'}})


if __name__=='__main__':main()