Example code
A reaction graph is not uniquely determined by its concentration dynamics. This example builds a four-species reactor from thirteen literal mass-action reactions, then constructs a different graph with exactly the same polynomial vector field. At a suitable stationary state the alternative graph is complex balanced: incoming and outgoing flux agree at every reaction complex (the group of molecules on one side of a reaction).
The reactor couples cores on and through shared reactions and reservoir exchanges. Two stationary currents decide whether any equivalent complex-balanced graph exists: the fork current must be nonnegative, and the reverse-channel current must not exceed . These obstructions apply even when competing graphs use arbitrarily many auxiliary complexes of any molecularity.
The code evaluates the exact six-parameter criterion using rational arithmetic, with both equality boundaries included. It also constructs the paper's nineteen-entry flux table, removes zero-flux edges and reconstructs rate constants. Incoming/outgoing balance and every polynomial coefficient are checked exactly. The default example has 16 active realizing edges, , , and positive AB production margins of each.
Realizability, productivity and stability differ. On the paper's rational slice, the system is realizable exactly when and AB-productive when , yet its unique equilibrium stays locally asymptotically stable throughout . Crossing the realizability boundary does not cause a local bifurcation.
Where an equivalent complex-balanced graph exists, it supplies a Horn–Jackson entropy: a function of concentrations that decreases toward equilibrium. Together with the paper's permanence argument, this yields convergence from every positive initial state. The package computes conservative eventual concentration bounds and illustrates convergence from two initial states with independent integrators. The paper also establishes a nearby open region with global stability but no such complex-balanced realization, but gives no explicit radius; the example does not assign that conclusion to arbitrary off-locus parameters.
Reusable components include the literal reactor, exact criterion and root enclosure, supporting-plane certificates, realization builder and private catalytic attachments. An attachment preserves the base field and relaxes its private species toward a prescribed concentration, allowing composition without changing the realizability classification.
Editable inputs are nondimensional mathematical examples. The package includes exact rates, sweep and trajectory data, seven scientific test groups and a detailed reuse guide. It does not rerun Lean, and deterministic equivalence does not imply identical stochastic behavior.
Python source
"""A reusable mass-action assembly and exact complex-balanced realization certificates."""
# EDITABLE INPUTS: nondimensional mathematical examples, not fitted laboratory rates.
PARAMETERS = ('7/16','7/16','1','12/25','1/2','1') # a,b,u,v,e,d
EXACT_REFERENCE = ('1/2','1/2','5/8','13/48') # set None to solve numerically after changing parameters
INITIAL_STATES = ((.1,1.2,.2,.8),(1.4,.08,1.1,.04))
TIME_HORIZON = 35.
TIME_SAMPLES = 351
SLICE_ETAS = ('9/2','5','11/2','8')
ATTACHMENTS = (('W1',(1,0,1,0),'2','3','1/10'),('W2',(0,1,0,0),'1','2','2'))
GRID_E = ('1/20','4',81)
GRID_D = ('1/20','5/2',61)
ROOT_BISECTIONS = 100
from dataclasses import dataclass
from fractions import Fraction as Q
from pathlib import Path
import argparse
import csv
import hashlib
import json
import math
import platform
import numpy as np
import sympy as sp
from scipy.integrate import solve_ivp
MANUSCRIPT_SHA256='ef9c1c93dfe470c759488cd4ad2a5352e81cd75b04412a78501fd68b13036770'
SPECIES=('A','B','z','H')
COMPLEX_NAMES=('0','A','B+z','z','H','2z','B','2A')
COMPLEXES=((0,0,0,0),(1,0,0,0),(0,1,1,0),(0,0,1,0),(0,0,0,1),(0,0,2,0),(0,1,0,0),(2,0,0,0))
@dataclass(frozen=True)
class Reaction:
source:tuple
target:tuple
rate:Q
class MassActionNetwork:
"""Finite reactions, polynomial coefficients and positive-domain integration."""
def __init__(self,species,reactions):
self.species=tuple(species);self.reactions=tuple(reactions);n=len(species)
if len(set(species))!=n:raise ValueError('Distinct species names required.')
for r in self.reactions:
if len(r.source)!=n or len(r.target)!=n or any(type(v) is not int or v<0 for v in r.source+r.target) or r.rate<0:raise ValueError('Natural exponents and nonnegative rates required.')
self.Y=np.array([r.source for r in self.reactions],int);self.S=np.array([np.array(r.target)-r.source for r in self.reactions],float).T
self.k=np.array([float(r.rate) for r in self.reactions]);self.n=n
def coefficients(self):
rows={}
for r in self.reactions:
row=rows.setdefault(r.source,[Q(0)]*self.n)
for i in range(self.n):row[i]+=r.rate*(r.target[i]-r.source[i])
return {y:tuple(c) for y,c in rows.items() if any(c)}
def exact_field(self,state):
state=tuple(map(Q,state))
return tuple(sum(math.prod(x**y for x,y in zip(state,source))*c[i] for source,c in self.coefficients().items()) for i in range(self.n))
def field(self,state):
x=np.asarray(state,float)
return self.S@(self.k*np.prod(x**self.Y,axis=1))
def jacobian(self,state):
x=np.asarray(state,float)
if np.any(x<=0):raise ValueError('Positive state required.')
rates=self.k*np.prod(x**self.Y,axis=1)
return self.S@(rates[:,None]*self.Y/x)
def balance(self,state):
state=tuple(map(Q,state));imbalance={}
for r in self.reactions:
flux=r.rate*math.prod(x**y for x,y in zip(state,r.source))
imbalance[r.source]=imbalance.get(r.source,Q(0))+flux
imbalance[r.target]=imbalance.get(r.target,Q(0))-flux
return imbalance
def integrate(self,initial,times,method='Radau'):
initial=np.asarray(initial,float);times=np.asarray(times,float)
if initial.shape!=(self.n,) or np.any(initial<=0) or times[0]!=0 or np.any(np.diff(times)<=0):raise ValueError('Positive initial state and increasing times from zero required.')
def fun(t,y):
x=np.exp(y);return self.field(x)/x
def jac(t,y):
x=np.exp(y);return self.jacobian(x)*x[None,:]/x[:,None]-np.diag(self.field(x)/x)
sol=solve_ivp(fun,(0,times[-1]),np.log(initial),method=method,jac=jac,t_eval=times,rtol=2e-10,atol=2e-12)
values=np.exp(sol.y.T)
if not sol.success or not np.all(np.isfinite(values)) or np.any(values<=0):raise ArithmeticError('Positive-domain integration failed.')
return values
@dataclass(frozen=True)
class Parameters:
a:Q;b:Q;u:Q;v:Q;e:Q;d:Q
def __post_init__(self):
for name in ('a','b','u','v','e','d'):
object.__setattr__(self,name,Q(getattr(self,name)))
if getattr(self,name)<=0:raise ValueError('All six parameters must be positive.')
@classmethod
def slice(cls,eta):
eta=Q(eta)
return cls((19-2*eta)/25,(eta-4)/25,1,5,eta,3)
@classmethod
def from_stationary_design(cls,A,B,z,u,v,e):
"""Choose rational concentrations and three rates; solve the other balances.
Raises if this design needs a nonpositive feed or output rate.
"""
A,B,z,u,v,e=map(Q,(A,B,z,u,v,e))
if min(A,B,z,u,v,e)<=0:raise ValueError('Positive design inputs required.')
K=A-B*z;J=e*(B-A*A);H=(u*z+2*v*z*z-K)/3
if H<=0:raise ValueError('Design gives nonpositive H.')
p=cls(2*A-B*z-2*J,(1+z)*B-A+J,u,v,e,(u*z+v*z*z)/H-2)
return p,(A,B,z,H)
def literal(self):
edges=((1,2,1),(2,1,1),(3,4,self.u),(4,3,1),(5,4,self.v),(4,5,1),
(0,1,self.a),(1,0,1),(0,6,self.b),(6,0,1),(6,7,self.e),(7,6,self.e),(4,0,self.d))
return MassActionNetwork(SPECIES,[Reaction(COMPLEXES[i],COMPLEXES[j],Q(k)) for i,j,k in edges])
def reduced(self,z):
z=Q(z)
if z<=-2:raise ValueError('Reduced functions require z > -2.')
B=(self.a+2*self.b)/(z+2);K=(self.v*(1+2*self.d)*z*z-self.u*(1-self.d)*z)/(2+self.d)
A=z*B+K;H=(self.u*z+self.v*z*z)/(2+self.d)
E=self.b-(1+self.e)*B+K+self.e*A*A
return (A,B,z,H),K,E
def cut(self):return max(Q(0),self.u*(1-self.d)/(self.v*(1+2*self.d)))
def locus(self):
z0=self.cut();E0=self.reduced(z0)[2];s=self.a+self.b;t=(self.a+2*self.b)*(self.e-1)/(self.e*s*s)-2
if E0>0:status='no admissible nonnegative-fork-current equilibrium';inside=False;margin=None
elif self.e<=1 or t<=z0:status='disguised toric';inside=True;margin=None
else:
margin=s-self.reduced(t)[0][0];inside=margin>=0
status='disguised toric' if inside else 'reservoir current budget violated'
return dict(inside=inside,status=status,z0=z0,E_z0=E0,t=t,threshold_margin=margin)
def admissible_root(self,steps=ROOT_BISECTIONS):
"""Exact rational enclosure; no claim about roots below the admissible cut."""
lo=self.cut()
if self.reduced(lo)[2]>0:raise ValueError('No root at or above the nonnegative-current cut; other stationary states are not excluded.')
if self.reduced(lo)[2]==0:return lo,lo
hi=lo+1+(1+self.e)*(self.a+2*self.b)/self.b
if self.reduced(hi)[2]<=0:raise ArithmeticError('Analytic upper bracket failed.')
for _ in range(steps):
mid=(lo+hi)/2;residual=self.reduced(mid)[2]
if residual==0:return mid,mid
if residual<0:lo=mid
else:hi=mid
return lo,hi
def permanence_bounds(self):
C=(self.e+Q(3,2))*self.a+(2*self.e+1)*self.b+(self.u+2)**2/(8*self.v)
rho=min(1/(self.e+Q(3,2)),1/(2*self.e+1),Q(1),self.d);M=2*C/rho
ma=self.a/(2*(2+2*self.e*M));mb=self.b/(2*(1+self.e+M));mz=ma/(2*(M+self.u+2*self.v*M));mh=self.u*mz/(2*(2+self.d))
return dict(C=C,rho=rho,eventual_upper=M,eventual_lower=min(ma,mb,mz,mh),individual_lowers=(ma,mb,mz,mh))
def currents(p,x):
A,B,z,H=map(Q,x);K=A-B*z;J=p.e*(B-A*A)
return dict(K=K,J=J,B_minus_J=B-J,AB_margins=(-K+2*J,K-J),zH_margins=(-p.u*z+3*H-2*p.v*z*z,p.u*z-2*H+p.v*z*z))
class SupportingCertificate:
def __init__(self,kind):
if kind=='K':self.h=((0,0,1,0),)+((0,1,0,0),)*5+((0,0,1,0),)*2;self.pi=(0,0,-1,0,0,0,0,0)
elif kind=='J':self.h=((1,1,1,1),(1,1,0,1),(-1,-1,0,1),(1,1,0,1),(1,1,0,1),(1,1,-1,-1),(-1,-1,1,1),(-1,-1,1,1));self.pi=(0,-1,-1,0,-1,0,-1,0)
else:raise ValueError('Certificate kind must be K or J.')
self.kind=kind
def extension(self,extra=()):
vertices=list(COMPLEXES);h=list(self.h);pi=list(self.pi)
for point in extra:
if len(point)!=4:raise ValueError('Auxiliary exponent must have four coordinates.')
values=[p-sum(a*(b-c) for a,b,c in zip(slope,point,y)) for p,slope,y in zip(self.pi,self.h,COMPLEXES)]
active=max(range(8),key=lambda i:values[i]);vertices.append(tuple(point));h.append(self.h[active]);pi.append(values[active])
slacks=[[sum(h[i][k]*(y[k]-x[k]) for k in range(4))+pi[j]-pi[i] for j,y in enumerate(vertices)] for i,x in enumerate(vertices)]
if min(min(row) for row in slacks)<0:raise ArithmeticError('Supporting inequality failed.')
return slacks
def pairing(self,p,state):
x=tuple(map(Q,state));coeff=p.literal().coefficients()
return sum(math.prod(xi**yi for xi,yi in zip(x,y))*sum(hi*ci for hi,ci in zip(h,coeff[y])) for y,h in zip(COMPLEXES,self.h))
class BalancedRealization:
"""Exact witness at an exact rational stationary state; never round a failed budget."""
def __init__(self,p,state):
self.p=p;self.state=tuple(map(Q,state));A,B,z,H=self.state
if min(self.state)<=0 or any(p.literal().exact_field(self.state)):raise ValueError('Exact positive stationary state required; a rounded numerical root is not an exact witness.')
values=currents(p,self.state);K=values['K'];J=values['J']
if K<0 or J>B:raise ValueError('Supporting-plane budget excludes every finite complex-balanced realization.')
w=B*z;U=p.u*z;V=p.v*z*z;jp=max(J,Q(0));jm=max(-J,Q(0));tau=max(Q(0),V-Q(3,2)*H)
self.fluxes={(0,1):p.a,(0,6):p.b,(1,2):w,(2,1):w,(1,6):K,(1,3):K,(1,0):w+jm,(1,7):jm,
(6,7):p.e*B-jp,(6,1):2*jp,(6,0):B-jp,(7,6):p.e*A*A,(3,4):U,(3,5):tau,(3,0):tau,
(5,4):V,(4,3):3*H-2*V+2*tau,(4,5):V-tau,(4,0):K-tau}
if min(self.fluxes.values())<0:raise ArithmeticError('Negative realizing flux.')
self.network=MassActionNetwork(SPECIES,[Reaction(COMPLEXES[i],COMPLEXES[j],q/math.prod(x**y for x,y in zip(self.state,COMPLEXES[i]))) for (i,j),q in self.fluxes.items() if q>0])
if any(self.network.balance(self.state).values()) or self.network.coefficients()!=p.literal().coefficients():raise ArithmeticError('Complex balance or full polynomial equality failed.')
def entropy(self,x):
x=np.asarray(x,float);ref=np.array(self.state,float);return float(np.sum(x*np.log(x/ref)-x+ref))
def entropy_derivative(self,x):return float(np.log(np.asarray(x,float)/np.array(self.state,float))@self.network.field(x))
@dataclass(frozen=True)
class CatalyticAttachment:
name:str
catalyst:tuple
forward:Q
reverse:Q
def __post_init__(self):
object.__setattr__(self,'forward',Q(self.forward));object.__setattr__(self,'reverse',Q(self.reverse))
if len(self.catalyst)!=4 or any(type(i) is not int or i<0 for i in self.catalyst) or min(self.forward,self.reverse)<=0:raise ValueError('Base-only natural catalyst exponents and positive rates required.')
@property
def stationary(self):return self.forward/self.reverse
@staticmethod
def compose(base,attachments):
attachments=tuple(attachments);m=len(attachments)
if base.species!=SPECIES:raise ValueError('Attach to the four-species base model.')
reactions=[Reaction(r.source+(0,)*m,r.target+(0,)*m,r.rate) for r in base.reactions]
for j,a in enumerate(attachments):
source=a.catalyst+(0,)*m;target=a.catalyst+tuple(int(i==j) for i in range(m))
reactions.extend([Reaction(source,target,a.forward),Reaction(target,source,a.reverse)])
return MassActionNetwork(SPECIES+tuple(a.name for a in attachments),reactions)
def slice_certificate():
eta=sp.Symbol('eta',positive=True)
A,B,z,H=sp.symbols('A B z H');a=(19-2*eta)/25;b=(eta-4)/25
F=sp.Matrix([a-2*A+B*z+2*eta*(B-A*A),b+A-(1+z)*B-eta*(B-A*A),A-B*z-z-10*z*z+3*H,z+5*z*z-5*H])
state={A:sp.Rational(2,5),B:sp.Rational(1,5),z:sp.Rational(1,5),H:sp.Rational(2,25)}
if any(sp.expand(v) for v in F.subs(state)):raise ArithmeticError('Symbolic slice stationarity failed.')
J=F.jacobian([A,B,z,H]).subs(state)
coeff=J.charpoly().all_coeffs();_,c1,c2,c3,c4=coeff;quantities=coeff[1:]+[sp.expand(c1*c2-c3),sp.expand(c1*c2*c3-c3*c3-c1*c1*c4)]
if not all(v>0 for expression in quantities for v in sp.Poly(expression,eta).all_coeffs()):raise ArithmeticError('All-positive Hurwitz coefficients failed.')
D=sp.diag(sp.Rational(5,2),5,5,sp.Rational(25,2));J5=J.subs(eta,5);energy=-(D*J5+J5.T*D)/2;minors=[energy[:i,:i].det() for i in range(1,5)]
if any(v<=0 for v in minors):raise ArithmeticError('Boundary quadratic entropy certificate failed.')
return dict(characteristic=list(map(str,coeff)),hurwitz=list(map(str,quantities)),boundary_energy_minors=list(map(str,minors)),scope='Local stability for the entire positive slice; the paper proves a non-explicit open globally stable nontoric neighborhood, not every off-locus parameter.')
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)
p=Parameters(*PARAMETERS);network=p.literal();classification=p.locus();root=None;witness=None
if EXACT_REFERENCE is not None:
reference=tuple(map(Q,EXACT_REFERENCE))
if any(network.exact_field(reference)):raise ValueError('Update EXACT_REFERENCE or set it to None after changing PARAMETERS.')
if classification['inside']:witness=BalancedRealization(p,reference)
else:
root=p.admissible_root();reference=p.reduced(sum(root)/2)[0]
times=np.linspace(0,TIME_HORIZON,TIME_SAMPLES);trajectories=[];differences=[];realization_difference=None
for initial in INITIAL_STATES:
values=network.integrate(initial,times);independent=network.integrate(initial,times,'BDF');trajectories.append(values);differences.append(float(np.max(abs(values-independent))))
if witness:realization_difference=float(np.max(abs(witness.network.integrate(INITIAL_STATES[0],times)-trajectories[0])))
certificates={}
for kind in ('K','J'):
certificate=SupportingCertificate(kind);slacks=certificate.extension();expanded=certificate.extension(((9,4,7,11),(100,0,23,5),(0,0,0,12)))
certificates[kind]=dict(slacks=slacks,off_diagonal_zeros=sum(slacks[i][j]==0 for i in range(8) for j in range(8) if i!=j),pairing=certificate.pairing(p,reference),auxiliary_pairs_checked=len(expanded)**2)
slice_state=tuple(map(Q,('2/5','1/5','1/5','2/25')));slice_rows=[]
for eta in SLICE_ETAS:
q=Parameters.slice(eta);margins=currents(q,slice_state);slice_rows.append((eta,int(q.locus()['inside']),*map(str,margins['AB_margins']),str(margins['B_minus_J']),float(max(np.linalg.eigvals(q.literal().jacobian(np.array(slice_state,float))).real))))
grid=[];elo,ehi=map(Q,GRID_E[:2]);dlo,dhi=map(Q,GRID_D[:2])
for j in range(GRID_D[2]):
d=dlo+(dhi-dlo)*j/(GRID_D[2]-1)
for i in range(GRID_E[2]):
e=elo+(ehi-elo)*i/(GRID_E[2]-1);entry=Parameters(Q(1,10),Q(1,10),5,1,e,d).locus();category=0 if entry['inside'] else (1 if entry['E_z0']>0 else 2)
grid.append((float(e),float(d),category))
attachments=[CatalyticAttachment(name,tuple(catalyst),a,b) for name,catalyst,a,b,initial in ATTACHMENTS]
augmented=CatalyticAttachment.compose(network,attachments);attachment_initial=list(INITIAL_STATES[0])+[float(Q(row[4])) for row in ATTACHMENTS]
attachment_states=augmented.integrate(attachment_initial,times)
if witness:
extended=CatalyticAttachment.compose(witness.network,attachments);xstar=witness.state+tuple(a.stationary for a in attachments)
if extended.coefficients()!=augmented.coefficients() or any(extended.balance(xstar).values()):raise ArithmeticError('Attachment realization failed.')
def write_json(name,data):(out/name).write_text(json.dumps(data,indent=2,default=str)+'\n')
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',['initial_index','time',*SPECIES,'entropy','entropy_derivative'],[(i,t,*x,witness.entropy(x) if witness else '',witness.entropy_derivative(x) if witness else '') for i,values in enumerate(trajectories) for t,x in zip(times,values)])
table('attachments.csv',['time',*augmented.species],np.column_stack([times,attachment_states]))
table('stable_slice.csv',['eta','disguised_toric','AB_A_margin','AB_B_margin','B_minus_J','spectral_abscissa'],slice_rows)
table('parameter_locus.csv',['e','d','category_0_toric_1_fork_2_reservoir'],grid)
if witness:
table('realization.csv',['source','target','equilibrium_flux','rate','present'],[(COMPLEX_NAMES[i],COMPLEX_NAMES[j],str(q),str(q/math.prod(x**y for x,y in zip(witness.state,COMPLEXES[i]))),int(q>0)) for (i,j),q in witness.fluxes.items()])
elif (out/'realization.csv').exists():
(out/'realization.csv').unlink() # This run has no exact witness; remove only its stale optional export.
results=dict(parameters=PARAMETERS,classification=classification,reference=reference,reference_scope='exact rational stationary state' if EXACT_REFERENCE is not None else 'midpoint lift of exact root enclosure; not an exact witness',root_enclosure=root,
currents=currents(p,reference),supporting_certificates=certificates,permanence=p.permanence_bounds(),slice=slice_certificate(),
exact_realization_built=witness is not None,active_realization_edges=len(witness.network.reactions) if witness else None,
numerical_checks=dict(independent_solver_differences=differences,equivalent_field_solver_difference=realization_difference,attachment_base_difference=float(np.max(abs(attachment_states[:,:4]-trajectories[0]))),minimum_concentration=float(min(x.min() for x in trajectories))),
scope='Exact rational classification, supporting inequalities and witness reconstruction; numerical trajectories illustrate conventional permanence/global-stability theorems. Lean not rerun. Deterministic equivalence does not identify stochastic processes.')
write_json('results.json',results)
lines=[f'Parameter classification: {classification["status"]}.',f'Exact rational realization built: {witness is not None}; active edges: {results["active_realization_edges"]}.',
f'Current budgets: K={results["currents"]["K"]}, B-J={results["currents"]["B_minus_J"]}; AB production margins={results["currents"]["AB_margins"]}.',
'All 128 source-pair supporting inequalities verified; active-piece extension checked on three auxiliary complexes.',
'Slice: realizable iff eta <= 5, strictly AB-productive iff 9/2 < eta < 9; locally stable throughout 4 < eta < 19/2.',
f'Independent solver differences: {differences}.', 'No explicit globally stable nontoric neighborhood radius is supplied. No stochastic equivalence claimed.']
(out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
fig,axs=plt.subplots(1,2,figsize=(10,4.3),layout='constrained')
g=np.array(grid);axs[0].pcolormesh(g[:,0].reshape(GRID_D[2],GRID_E[2]),g[:,1].reshape(GRID_D[2],GRID_E[2]),g[:,2].reshape(GRID_D[2],GRID_E[2]),cmap=ListedColormap(['#3e8585','#cf8952','#866b9c']),shading='nearest',vmin=0,vmax=2)
axs[0].set(xlabel='Reverse-channel rate e',ylabel='Output rate d',title='Existence of an equivalent complex-balanced\nmodel')
from matplotlib.patches import Patch
axs[0].legend(handles=[Patch(color=c,label=s) for c,s in zip(['#3e8585','#cf8952','#866b9c'],['Realizable','Fork budget fails','Reservoir budget fails'])],fontsize=7)
etas=np.linspace(4.01,9.49,200);abscissae=[max(np.linalg.eigvals(Parameters.slice(str(e)).literal().jacobian(np.array(slice_state,float))).real) for e in etas]
axs[1].plot(etas,abscissae);axs[1].axvline(5,color='#866b9c',ls='--',label='Realizability boundary');axs[1].axhline(0,color='gray',lw=.7);axs[1].set(xlabel='Slice parameter eta',ylabel='Spectral abscissa',title='Equilibrium stability across the\nrealizability boundary');axs[1].legend(fontsize=7);axs[1].grid(alpha=.2)
fig.savefig(out/'classification.png',dpi=180);fig.savefig(out/'classification.svg');plt.close(fig)
fig,axs=plt.subplots(1,2,figsize=(10,4.3),layout='constrained')
for i,name in enumerate(SPECIES):
line,=axs[0].plot(times,trajectories[0][:,i],label=name)
for values in trajectories[1:]:axs[0].plot(times,values[:,i],ls='--',color=line.get_color(),lw=.8)
axs[0].set(xlabel='Time',ylabel='Concentration',title=f'{len(trajectories)} positive initial states');axs[0].legend()
if witness:
for i,values in enumerate(trajectories):axs[1].semilogy(times,np.maximum([witness.entropy(x) for x in values],1e-16),label=f'Initial state {i+1}')
axs[1].set(ylabel='Horn–Jackson entropy (display floor 1e-16)',title='Horn–Jackson entropy along reactor\ntrajectories')
else:
for i,a in enumerate(attachments):axs[1].plot(times,attachment_states[:,4+i],label=a.name)
axs[1].set(ylabel='Private-species concentration',title='Catalytic attachments')
axs[1].set_xlabel('Time');axs[1].legend()
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'trajectories.png',dpi=180);fig.savefig(out/'trajectories.svg');plt.close(fig)
digest=lambda path:hashlib.sha256(path.read_bytes()).hexdigest()
write_json('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'}))
if __name__=='__main__':main()
Run output
Parameter classification: disguised toric. Exact rational realization built: True; active edges: 16. Current budgets: K=3/16, B-J=3/8; AB production margins=(Fraction(1, 16), Fraction(1, 16)). All 128 source-pair supporting inequalities verified; active-piece extension checked on three auxiliary complexes. Slice: realizable iff eta <= 5, strictly AB-productive iff 9/2 < eta < 9; locally stable throughout 4 < eta < 19/2. Independent solver differences: [1.5610586157066564e-09, 1.4354611943545592e-09]. No explicit globally stable nontoric neighborhood radius is supplied. No stochastic equivalence claimed.