Example code
Coupled autocatalytic reactors can keep producing through repeated harvesting even when the model includes an additional reaction intermediate. Its stored material must be tracked through recovery and collection. This example implements the paper's reversible seven-species exporter, common molecular exchange, withdrawal/loss/refill pulses, and repeated recovery from actual returned states.
The code distinguishes productive stock, collected inventory, free product, carried intermediate, food input, and gross reaction service. Product marks are existing washout events in the last unit of each four-unit cycle. Intermediate washout is not credited as product, and free-X collection is already included in inventory collection.
For the paper's connected two-reactor, 100-cycle instance, the finite-mission theorem gives a success probability above 99.5% at the stated count scale. Fresh rational arithmetic separately checks a sufficient exponential bound proving more than 99%. The example exposes the scale calculation and integer quotas instead of attempting to simulate that enormous molecular population.
The seven-species model carries the intermediate through every reaction, pulse, and cycle. A separate six-species donor solver demonstrates the exact reduction driven by the recorded material inventories: four productive phases plus the actual material histories recover the full trajectory. An explicit pair of ready states shows why simply deleting the refined intermediate cannot give the same closed dynamics or projected jump law.
Every successful mission of at least 55 cycles must create net new inventory beyond its initial supply. The code checks the inventory balance summed across all cycles, including intermediate storage, all washout, withdrawal, and pulse loss. For the worked mission the worst-ready lower bound is 3.652 × 10¹⁴ newly created inventory equivalents.
Editable inputs and composable classes support new symmetric graphs, node rates, admitted policies, deterministic missions, and optional small-count molecular simulations. Event-budget exhaustion is reported as unfinished, and failed cycles remain in the simulated law. The package includes source, scientific checks, full-state CSV outputs, theorem sizing, and reusable accounting components.
The probability theorem is imported from the manuscript; Lean is not rerun. Numerical trajectories illustrate the literal equations and do not estimate failure probabilities. The scale is sufficient rather than optimized or experimentally calibrated, and the result concerns finite missions beginning in readiness.
Python source
"""Coupled refined reactors: useful operation, exact ledgers, and honest certificates."""
# EDITABLE INPUTS: normalized paper example, not calibrated chemistry.
COUNT_SCALE=224_000_000_000_000
MISSION_CYCLES=100
FAILURE_ALLOWANCE='1/100'
EXCHANGE=((0.,1.),(1.,0.)) # same symmetric graph for every species
RELEASE=('20','20')
CLEAVAGE=('3/100','3/100')
THETA='1/100' # stochastic theorem fixes this value
INITIAL=(('19/20','19/20','1/20','0','0','0','0'),
('13/14','13/14','0','0','1/28','0','0'))
DENSITY_CYCLES=10
FEEDBACK_POLICY=False # True retains more after a small actual prior collection
SMALL_COUNT_SCALE=1400 # optional SSA, far below useful theorem scale
SMALL_CYCLES=2
EVENT_BUDGET=200000
RANDOM_SEED=51092026
CONCENTRATION_MOLAR='1/1000' # illustrative unit conversion
TIME_UNIT_SECONDS=60.
import argparse
from fractions import Fraction as F
import csv
import hashlib
import json
import math
from pathlib import Path
import platform
import numpy as np
from mpmath import iv
import sympy as sp
from reactor import (Chemistry,CommonExchange,Intervention,HistoryPolicy,ReactorNetwork,
DeterministicMission,MolecularMission,ready,A,B,I,J,Y,EPS,ETA,ACCOUNTS,SPECIES)
iv.dps=60
MANUSCRIPT_SHA256='6939f07afd2c656577771d525062a8375f3ec59a0a20e453f06c736e638d0820'
def interval(q):
q=F(q);return iv.mpf(q.numerator)/q.denominator
def lower(x):return math.nextafter(float(x.a),-math.inf)
def upper(x):return math.nextafter(float(x.b),math.inf)
class MissionCertificate:
"""Imported theorem interface; numerical tail evaluated with outward bounds."""
def __init__(self,network):self.network=network
def bound(self,V,cycles,initial):
if type(V) is not int or V<10000 or type(cycles) is not int or cycles<0:raise ValueError('Theorem requires integer V>=10000 and cycles>=0')
net=self.network
if not all(c.scope(True) for c in net.chemistries):raise ValueError('Rates or theta outside stochastic theorem')
if not ready(initial,net.chemistries[0].theta,V):raise ValueError('Every initial node must be ready')
refined=net.size==7;pref=29 if refined else 27;C=2*10**12 if refined else 2*10**14
degree=net.exchange.exact_degree;exponent=F(V)/(C*(1+degree)**3)
if cycles==0:prob=1.;tail=0.
else:
tail=upper(interval(pref*net.n*cycles)*iv.exp(-interval(exponent)))
prob=max(0.,math.nextafter(1-tail,-math.inf))
synthesis=net.n*cycles*(-(-V//56))+net.n*F(V,28)-sum(int(row@J[:net.size]) for row in initial)
worst=net.n*V*(F(cycles,56)+F(1,28)-F(161,160)) if cycles else None
return dict(lower_probability=prob,failure_union_upper=tail,exponent=exponent,
node_count=net.n,degree_bound=degree,cycles=cycles,V=V,refined=refined,
required_QI_per_node_cycle=-(-V//56),required_QX_per_node_cycle=-(-V//1080),
max_each_food=5*V,max_gross_service=V//5,
net_synthesis_lower_using_initial=synthesis if cycles else None,
net_synthesis_lower_worst_ready=worst,
scope='Theorem applied to ready initial counts and admitted literal history policy; not estimated from trajectories')
def sufficient_scale(self,cycles,delta):
if type(cycles) is not int or cycles<1 or not 0<F(delta)<1:raise ValueError('Positive mission and failure allowance in (0,1)')
net=self.network;pref=29 if net.size==7 else 27;C=2*10**12 if net.size==7 else 2*10**14
q=F(1)+net.exchange.exact_degree
val=interval(C*q**3)*iv.ln(interval(F(pref*net.n*cycles)/F(delta)))
return max(10000,math.ceil(upper(val)))
def exact_checks():
chemistry=Chemistry();s=sp.symbols('u w x c1 c2 z h');V=sp.Symbol('V',positive=True)
rates=[sp.Rational(c.k.numerator,c.k.denominator)*sp.prod(s[i] for i in c.inputs) for c in chemistry.channels]
f=[sp.expand(sum(rate*c.jump[i] for rate,c in zip(rates,chemistry.channels))) for i in range(7)]
assert sp.expand(sum(int(a)*v for a,v in zip(A,f))-(1-sum(int(a)*x for a,x in zip(A,s)) ))==0
assert sp.expand(sum(int(a)*v for a,v in zip(B,f))-(1-sum(int(a)*x for a,x in zip(B,s)) ))==0
count_rates=[]
for c in chemistry.channels:
count_rates.append(sp.Rational(c.k.numerator,c.k.denominator)*V**(1-len(c.inputs))*sp.prod(V*s[i]-c.inputs[:k].count(i) for k,i in enumerate(c.inputs)))
correction=[sp.simplify(sum(rate*c.jump[i]/V for rate,c in zip(count_rates,chemistry.channels))-f[i]) for i in range(7)]
assert correction==[0,0,40*s[2]/V,0,0,-20*s[2]/V,0]
yw=[F(0),F(0),F(1),F(9,8),F(7,5),F(9,5),F(0)]
assert sp.simplify(sum(w*c for w,c in zip(yw,correction))-4*s[2]/V)==0
increments=[sum(int(w)*v for w,v in zip(J,c.jump)) for c in chemistry.channels if c.direction==1]
assert increments==[1,0,0,1,0,0,-1]
# Finite adjoint: polynomial row, rational exponential upper bound, and noise margins.
H=sp.Matrix([[0,20,0,38],[0,5,20,0],[0,0,7,2],[0,0,20,24]])
row=(sum((H**j*sp.Rational(1,28)**j/math.factorial(j) for j in range(5)),sp.zeros(4)))[0,:]
expect=[F(1),F(961355,1229312),F(1847785,1843968),F(665611,307328)]
assert list(row)==[sp.Rational(q.numerator,q.denominator) for q in expect]
expupper=sum(F(6,7)**j/math.factorial(j) for j in range(8))+F(6,7)**8*F(9,8)/math.factorial(8)
assert expupper<=F(2357,1000) and F(2357,1000)**2<=F(961355,172872)
assert min(a/b for a,b in zip(expect,yw[2:6]))==F(961355,1382976)
freeX=F(5145,100000)/8-F(8,10000)
assert freeX==F(901,160000)
assert F(3575,100000)>F(1,56) and freeX-F(1,1000)>F(1,1080) and F(181,1000)<F(1,5)
# Multiplying by (1+Delta) makes the graph-noise inequality affine; coefficients suffice.
D=sp.Symbol('D',nonnegative=True)
eps=1/(10000*(1+D));eta=sp.Rational(11,2)*eps
surplus=sp.expand((sp.Rational(2,3)-sp.Rational(11,20))*sp.Rational(59,5000)*(1+D)-(2*D+sp.Rational(2,3))*sp.Rational(11,20000))
assert all(c>0 for c in sp.Poly(surplus,D).all_coeffs())
assert F(361,5120000*10**8)>=F(1,2*10**12)
taylor=sum(F(14)**j/math.factorial(j) for j in range(16))
assert taylor>600000 and 1-F(5800,600000)>F(99,100)
assert F(54,56)+F(1,28)-F(161,160)<0<F(55,56)+F(1,28)-F(161,160)
# Same six coordinates, distinct h: exact projected derivatives and projected event rates.
h=F(1,200);d=F(3,100);beta=theta=F(1,100)
projection_difference=[d*h/theta,d*h/theta,d*beta*h/theta,F(0),F(0),F(0)]
assert projection_difference[:3]==[F(3,200),F(3,200),F(3,20000)]
return dict(material_identities=True,falling_factorial_correction=[str(v) for v in correction],
inventory_forward_increments=increments,phase_polynomial_row=expect,exp_6_over_7_upper=expupper,
phase_free_X_floor=freeX,noise_margin_polynomial=str(surplus),witness_exp14_lower=taylor,
witness_probability_rational_lower=1-F(5800,600000),net_synthesis_positive_from_cycle=55,
nonclosure_projected_derivative_difference=projection_difference,
evidence='Fresh symbolic/rational identities; probability theorem imported, Lean not rerun')
def literal_initial(V,refined=True):
vals=[[F(x)*V for x in row] for row in INITIAL]
if any(v.denominator!=1 for row in vals for v in row):raise ValueError('Initial density is not integral at this V')
return np.array([[int(v) for v in row[:7 if refined else 6]] for row in vals],dtype=np.int64)
def main():
parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',default='outputs');parser.add_argument('--simulate',action='store_true');args=parser.parse_args()
out=Path(args.output);out.mkdir(parents=True,exist_ok=True)
graph=CommonExchange(EXCHANGE);net=ReactorNetwork([Chemistry(r,d,THETA) for r,d in zip(RELEASE,CLEAVAGE)],graph)
policy=HistoryPolicy(7,FEEDBACK_POLICY);initial=np.array([[float(F(v)) for v in row] for row in INITIAL])
certificate=MissionCertificate(net);bound=certificate.bound(COUNT_SCALE,MISSION_CYCLES,literal_initial(COUNT_SCALE))
run=DeterministicMission(net,policy).run(initial,DENSITY_CYCLES)
donor=ReactorNetwork([Chemistry(r,d,refined=False) for r,d in zip(RELEASE,CLEAVAGE)],graph)
post=np.array([Intervention(survival=('49/50',)*6).deterministic(row[:6])[0] for row in initial])
full=donor.evolve(post);reduced=donor.donor_reduced(post)
reduction_error=float(np.max(abs(full['end']-reduced)))
refined_first=run['traces'][0]
# Finite theta is retained: compare the illustrative endpoints, without convergence claims.
theta_rows=[]
for theta in ['1/100','1/1000','1/10000']:
other=ReactorNetwork([Chemistry(r,d,theta) for r,d in zip(RELEASE,CLEAVAGE)],graph)
after=np.array([Intervention().deterministic(row)[0] for row in initial]);flow=other.evolve(after)
theta_rows.append(dict(theta=theta,max_endpoint_difference_to_donor=float(np.max(abs(flow['end'][:,:6]-full['end']))),
carried_intermediate=flow['end'][:,6].tolist()))
result=dict(inputs=dict(V=COUNT_SCALE,mission_cycles=MISSION_CYCLES,theta=THETA,release=RELEASE,cleavage=CLEAVAGE,exchange=EXCHANGE,initial=INITIAL),
exact=exact_checks(),refined_certificate=bound,
donor_certificate=MissionCertificate(donor).bound(COUNT_SCALE,MISSION_CYCLES,literal_initial(COUNT_SCALE,False)),
sufficient_scale=certificate.sufficient_scale(MISSION_CYCLES,FAILURE_ALLOWANCE),
donor_reduction_endpoint_error=reduction_error,theta_illustrations=theta_rows,
numerical_mission={k:v for k,v in run.items() if k not in ('traces','final')},
units=dict(concentration_molar=CONCENTRATION_MOLAR,cycle_seconds=4*TIME_UNIT_SECONDS,
node_volume_liters=float(F(COUNT_SCALE)/(F('6.02214076e23')*F(CONCENTRATION_MOLAR)))))
if args.simulate:
result['small_count_path']=MolecularMission(net,policy,SMALL_COUNT_SCALE,RANDOM_SEED).run(literal_initial(SMALL_COUNT_SCALE),SMALL_CYCLES,EVENT_BUDGET)
dump=lambda name,obj:(out/name).write_text(json.dumps(obj,indent=2,default=lambda v:v.tolist() if isinstance(v,np.ndarray) else str(v) if isinstance(v,F) else float(v))+'\n',encoding='utf-8')
dump('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('cycles.csv',['cycle','node',*ACCOUNTS,'food_U','food_W',*SPECIES],
[(h['cycle'],i,*h['accounts'][i],*h['food'][i],*h['endpoint'][i]) for h in run['history'] for i in range(graph.n)])
table('first_cycle.csv',['t','node',*SPECIES,'material_A','material_B','stock_Y','inventory_I'],
[(t,i,*row,row@A,row@B,row@Y,row@I) for t,cs in zip(refined_first['times'],refined_first['states']) for i,row in enumerate(cs)])
sizes=[]
for degree in [0,1,4]:
other=ReactorNetwork(net.chemistries,CommonExchange([[0,degree],[degree,0]]));cert=MissionCertificate(other)
for m in [1,3,10,30,100,300,1000,3000,10000]:sizes.append([degree,m,cert.sufficient_scale(m,'1/100')])
table('sizing.csv',['maximum_weighted_degree','cycles','sufficient_V'],sizes)
lines=[f'Refined mission lower probability: {bound["lower_probability"]:.10f}; donor bound at the same scale: {result["donor_certificate"]["lower_probability"]:.10f}.',
f'Sufficient V for failure allowance {FAILURE_ALLOWANCE}: {result["sufficient_scale"]}.',
f'Worst-ready net synthesis lower bound: {bound["net_synthesis_lower_worst_ready"]} inventory equivalents.',
f'Donor full versus material-forced reduction endpoint gap: {reduction_error:.3g}.',
f'Deterministic ledger residuals: inventory {run["inventory_telescope_residual"]:.3g}, intermediate storage {run["storage_telescope_residual"]:.3g}.',
'Illustrations are deterministic; the probability theorem is imported. Lean is not rerun.']
(out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
plot(out,run,sizes)
digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,python=platform.python_version(),
source_sha256=digest(Path(__file__)),module_sha256={p.name:digest(p) for p in Path(__file__).parent.glob('*.py')},
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,run,sizes):
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
first=run['traces'][0];t=first['times'];c=first['states']
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
for i in range(c.shape[1]):
axs[0].plot(t,(c@Y)[:,i],label=f'Node {i+1}')
axs[1].plot(t,c[:,i,6],label=f'Node {i+1}')
axs[0].plot(t,np.minimum(.0118*np.exp(.55*t),.052)-.00055,'--',color='k',label='Good-event stock bound')
axs[1].plot(t,.009+.00201*np.exp(-3.02*t)+.0003,'--',color='k',label='Good-event intermediate bound')
for ax,title,yl in zip(axs,['Catalytic stock during recovery','Intermediate concentration during recovery'],['Stock Y','Intermediate D']):
ax.set(title=title,xlabel='Normalized time',ylabel=yl);ax.axvspan(3,4,alpha=.08);ax.grid(alpha=.2);ax.legend(fontsize=7)
fig.savefig(out/'refined_flow.png',dpi=180);fig.savefig(out/'refined_flow.svg');plt.close(fig)
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
h=run['history'];k=[v['cycle'] for v in h]
for i in range(c.shape[1]):
axs[0].plot(k,[v['accounts'][i][0] for v in h],'-o',label=f'Inventory, node {i+1}')
axs[0].plot(k,[v['accounts'][i][1] for v in h],'--',label=f'Free X, node {i+1}')
axs[0].set(title='Collected output over successive reactor\ncycles',xlabel='Cycle',ylabel='Collected normalized amount');axs[0].legend(fontsize=7)
for degree in [0,1,4]:
a=np.array([r for r in sizes if r[0]==degree]);axs[1].loglog(a[:,1],a[:,2],label=f'Degree ≤ {degree}')
axs[1].scatter([100],[224000000000000],color='k',s=20,label='Paper witness')
axs[1].set(title='A sufficient count scale',xlabel='Mission cycles',ylabel='Sufficient count scale V');axs[1].legend(fontsize=7)
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'mission.png',dpi=180);fig.savefig(out/'mission.svg');plt.close(fig)
if __name__=='__main__':main()
Run output
Refined mission lower probability: 0.9951771334; donor bound at the same scale: 0.0000000000. Sufficient V for failure allowance 1/100: 212332534120362. Worst-ready net synthesis lower bound: 365200000000000 inventory equivalents. Donor full versus material-forced reduction endpoint gap: 1.88e-12. Deterministic ledger residuals: inventory -1.67e-13, intermediate storage 0. Illustrations are deterministic; the probability theorem is imported. Lean is not rerun.