Example code
The reactor makes a covalent template from foods , uses that template to promote further ligation, and exports product through washout. A reversible cleavage reaction, powered by a maintained external drive, competes with the cycle by returning free to its food components. The example implements the twenty labelled count channels, their deterministic density limit, and separate ledgers for covalent export, food input and gross driving service.
The operational claim is joint: establish catalysis by time 500, keep both resources within specified concentration ranges, retain catalytic material after entry, export enough in every completed unit window, and stay within supply budgets. Each window starts from the actual continuing chemical state; neither chemistry nor cumulative supplies are restarted.
At over 100 completed windows, interval arithmetic that rounds bounds outwards confirms a joint success probability of at least 0.9238. The ordinary evaluation is about 0.923849. Removing only the templated ligation pair makes the same output schedule have probability below , even when every resource exit is credited as a possible success. The comparison keeps basal reactions, binding, duplex release, food inputs, washout and both driving directions.
The guaranteed output is at least 20000 covalent mass units per window, corresponding to 5000 template equivalents. Through the full 600 time units, each food allowance is events and the gross driving allowance is events. Opposite driving events cannot cancel this cost. Resident inventory stays of order ; replenishment over a long run grows with elapsed time.
Editable inputs select the kinetic parameters, copy scale, duration and target confidence. A requested failure probability of 0.01 does not pass at the default copy scale. The full expression individually certifies as sufficient in the example's bounded search; the simpler rule gives . Neither value is claimed to be minimal. Startup and shared failures are charged once, while window errors are added without assuming independence.
The package separates reaction channels, event monitoring, probability evaluation and production accounting. Tests check deadline entry, right-endpoint window assignment, return exactly to the residence floor, persistent failure history, supply saturation and material balances. An optional count simulator explicitly reports unfinished paths when its computational budget runs out.
The ODE curves illustrate the mechanism and ledger, not finite-copy success frequencies. At a reference concentration of one micromolar and dilution of one per hour, the default scale translates to about 166 pL and a conditional equivalent-product rate of 50 pM/hour. Those are unit conversions, without experimental calibration or a claim of purified free product. Seven scientific test groups accompany the code; the upstream Lean proof is not rerun, and exponential finite operation does not imply perpetual survival.
Python source
"""Driven-cleavage reactor: literal counts, repeated-window ledger, and bounds."""
from __future__ import annotations
# EDITABLE INPUTS -------------------------------------------------------------
COPY_SCALE = 100_000_000
KINETIC_BIAS = '10'
DUPLEX_RELEASE = '20'
DRIVEN_CLEAVAGE = '1/20'
POST_STARTUP_DURATION = '100'
TARGET_FAILURE = '1/100'
BASAL_RATE = '1/500000000' # fixed by the theorem; changing it exits its scope
REFERENCE_MOLAR = '1/1000000'
DILUTION_PER_HOUR = '1'
TEMPERATURE_K = '298'
SSA_EVENT_BUDGET = 20_000
SSA_SEED = 25092026
MANUSCRIPT_SHA256 = 'b03c9193c19090688e71a36a46ad31bf33dc9bde1d910c5b0cf3615ae803b05e'
# ---------------------------------------------------------------------------
from dataclasses import dataclass
from fractions import Fraction as Q
import argparse
import csv
import hashlib
import json
import math
from pathlib import Path
import platform
import numpy as np
from scipy.integrate import solve_ivp
from mpmath import mp, iv
mp.dps=70;iv.dps=70
SPECIES=('U','W','X','C1','C2','Z')
A_WEIGHTS=(1,0,1,2,2,2);B_WEIGHTS=(0,1,1,1,2,2)
Y_WEIGHTS=(Q(0),Q(0),Q(1),Q(9,8),Q(7,5),Q(9,5))
M_WEIGHTS=(2,2,4,6,8,8);S_WEIGHTS=(0,0,4,4,4,8)
LEDGER=('export','food_U','food_W','drive_F','drive_P','wash_A','wash_B','wash_mass')
def dot(a,b):return sum(x*y for x,y in zip(a,b))
def number(ctx,x):
x=Q(str(x));return ctx.mpf(x.numerator)/x.denominator
def logsum(ctx,values):
values=list(values);shift=values[0]
return shift+ctx.log(sum(ctx.exp(v-shift) for v in values))
@dataclass(frozen=True)
class Parameters:
K: object=KINETIC_BIAS
release: object=DUPLEX_RELEASE
delta: object=DRIVEN_CLEAVAGE
epsilon: object=BASAL_RATE
def __post_init__(self):
for name in ('K','release','delta','epsilon'):
q=Q(str(getattr(self,name)))
if q<=0:raise ValueError('Positive reaction coefficients required.')
object.__setattr__(self,name,q)
@property
def eta(self):return self.epsilon/16
def theorem_scope(self):
return 8<=self.K<=12 and 18<=self.release<=22 and Q(1,100)<=self.delta<=Q(1,20) and self.epsilon==Q(1,500000000)
@dataclass(frozen=True)
class Channel:
label: int
name: str
inputs: tuple[int,...]
outputs: tuple[int,...]
coefficient: Q
marks: tuple[int,...]=(0,)*8
@property
def jump(self):return tuple(self.outputs.count(i)-self.inputs.count(i) for i in range(6))
def count_rate(self,counts,V):
rate=self.coefficient*Q(V)**(1-len(self.inputs));seen={}
for i in self.inputs:
rate*=counts[i]-seen.get(i,0);seen[i]=seen.get(i,0)+1
return rate
class DrivenReactor:
"""Preserves original labels after deleting only channels 6 and 7."""
def __init__(self,parameters=Parameters(),enabled=True):
self.p=parameters;self.enabled=enabled;p=parameters
pairs=(('basal',(0,1),(2,),p.epsilon,p.epsilon/p.K),
('bind_U',(2,0),(3,),Q(20),Q(20)),
('bind_W',(3,1),(4,),Q(20),Q(20)),
('ligation',(4,),(5,),Q(20),Q(20)/p.K),
('release',(5,),(2,2),p.release,p.release))
channels=[]
for j,(name,left,right,kp,km) in enumerate(pairs):
channels.extend([Channel(2*j,name+'_forward',left,right,kp),Channel(2*j+1,name+'_reverse',right,left,km)])
for i in (0,1):
marks=[0]*8;marks[1+i]=1
channels.append(Channel(10+i,'feed_'+SPECIES[i],(),(i,),Q(1),tuple(marks)))
for i in range(6):
marks=(S_WEIGHTS[i],0,0,0,0,A_WEIGHTS[i],B_WEIGHTS[i],M_WEIGHTS[i])
channels.append(Channel(12+i,'wash_'+SPECIES[i],(i,),(),Q(1),marks))
channels.extend([Channel(18,'driven_cleavage',(2,),(0,1),p.delta,(0,0,0,1,0,0,0,0)),
Channel(19,'driven_reverse',(0,1),(2,),p.delta*p.eta,(0,0,0,0,1,0,0,0))])
self.channels=tuple(c for c in channels if enabled or c.label not in (6,7))
self.by_label={c.label:c for c in self.channels}
self.jumps=np.array([c.jump for c in self.channels],int);self.marks=np.array([c.marks for c in self.channels],int)
self.coefficients=np.array([float(c.coefficient) for c in self.channels])
self.orders=np.array([len(c.inputs) for c in self.channels])
self.inputs=np.full((len(self.channels),2),6,int);self.offsets=np.zeros((len(self.channels),2),int)
for j,c in enumerate(self.channels):
seen={}
for k,i in enumerate(c.inputs):
self.inputs[j,k]=i;self.offsets[j,k]=seen.get(i,0);seen[i]=seen.get(i,0)+1
def rates(self,state,V=None):
factors=np.r_[state,1][self.inputs]
if V is not None:factors=np.maximum(0,factors-self.offsets)
rate=self.coefficients*np.prod(factors,axis=1)
return rate if V is None else rate*np.power(float(V),1-self.orders)
def generator(self,counts,V,weights):
rates=[c.count_rate(counts,V) for c in self.channels]
jumps=[dot(weights,c.jump) for c in self.channels]
return sum(a*d for a,d in zip(rates,jumps)),sum(a*d*d for a,d in zip(rates,jumps))
def deterministic(self,H,method='Radau'):
H=float(Q(str(H)))
if not 1<=H<=10000:raise ValueError('ODE illustration is limited to duration 1..10000; bounds have a separate interface.')
def rhs(t,state):
rates=self.rates(state[:6]);return np.r_[[email protected],[email protected]]
def entry(t,state):return np.array(Y_WEIGHTS,float)@state[:6]-1/2500
entry.direction=1
sol=solve_ivp(rhs,(0,500+H),[1,1,0,0,0,0]+[0]*8,method=method,dense_output=True,events=entry,
rtol=1e-10,atol=1e-16,max_step=1)
if not sol.success:raise ArithmeticError(sol.message)
return sol
def stochastic(self,policy,limit=SSA_EVENT_BUDGET,seed=SSA_SEED):
if not 1<=policy.V<10**12 or limit<1:raise ValueError('Direct SSA needs V<10^12 and a positive event budget.')
rng=np.random.default_rng(seed);monitor=OperatingMonitor(self,policy);events=0
while monitor.status=='active' and events<limit:
rates=self.rates(monitor.counts,policy.V);total=float(sum(rates))
event_time=monitor.time+rng.exponential(1/total)
boundary=float(monitor.next_boundary())
if event_time>boundary:
monitor.advance(boundary);continue
j=int(np.searchsorted(np.cumsum(rates),rng.random()*total,side='right'))
monitor.event(event_time,self.channels[j].label);events+=1
report=monitor.report();report.update(events=events)
if monitor.status=='active':report.update(status='unfinished',event_verdict=None)
return report
@dataclass(frozen=True)
class ObservationPolicy:
V: int=COPY_SCALE
H: object=POST_STARTUP_DURATION
def __post_init__(self):
object.__setattr__(self,'H',Q(str(self.H)))
if not isinstance(self.V,int) or self.V<1 or self.H<1:raise ValueError('Positive integer copy scale and H>=1 required.')
@property
def end(self):return 500+self.H
@property
def windows(self):return self.H.numerator//self.H.denominator
@property
def export_cap(self):return math.ceil(Q(self.V,5000))
@property
def supply_cap(self):return math.ceil(2*self.V*self.end)+1
class OperatingMonitor:
"""Tracks one continuing path. Boundaries reset only the window counter.
Events exactly at a window's right endpoint belong to that window. Entry
at time 500 is allowed. Return at Y<=V/5000 fails (including equality).
Natural ledgers provide a second audit of the capped measurement counters.
"""
def __init__(self,reactor,policy):
if policy.windows>100000:raise ValueError('Monitor window limit; use ProbabilityBudget for huge horizons.')
self.reactor=reactor;self.policy=policy;self.counts=[policy.V,policy.V,0,0,0,0]
self.ledger=[0]*8;self.capped_supplies=[0]*4;self.counter=0
self.time=0.;self.phase='pre_entry';self.entry_time=None;self.status='active'
self.started=False;self.completed=0;self.all_windows_pass=True;self.window_exports=[]
def next_boundary(self):
if not self.started:return Q(500)
return Q(501+self.completed) if self.completed<self.policy.windows else self.policy.end
def corridor(self):
V=self.policy.V
return all(9*V<=10*dot(w,self.counts)<=11*V for w in (A_WEIGHTS,B_WEIGHTS))
def _boundary(self):
if not self.started:
self.started=True
if self.reactor.enabled and self.phase=='pre_entry':self.status='missed_deadline'
elif self.completed<self.policy.windows:
self.window_exports.append(self.counter)
self.all_windows_pass &= self.counter>=self.policy.export_cap
self.completed+=1;self.counter=0
if self.status=='active' and self.time>=float(self.policy.end):self.status='complete'
def advance(self,t):
if t<self.time or t>float(self.policy.end):raise ValueError('Invalid monitor time.')
while self.status=='active' and float(self.next_boundary())<=t:
self.time=float(self.next_boundary());self._boundary()
if self.status=='active':self.time=t
def event(self,t,label):
if self.status!='active' or t<self.time or t>float(self.policy.end):raise ValueError('Event outside active observation interval.')
while self.status=='active' and float(self.next_boundary())<t:
self.time=float(self.next_boundary());self._boundary()
if self.status!='active':return
channel=self.reactor.by_label[label]
if channel.count_rate(self.counts,self.policy.V)<=0:raise ValueError('Reaction is not supported by molecule counts.')
self.time=t;self.counts=[n+d for n,d in zip(self.counts,channel.jump)]
if min(self.counts)<0:raise AssertionError('Count transition left the orthant.')
self.ledger=[n+d for n,d in zip(self.ledger,channel.marks)]
self.capped_supplies=[min(self.policy.supply_cap,n+d) for n,d in zip(self.capped_supplies,channel.marks[1:5])]
if t>500 and self.completed<self.policy.windows:
self.counter=min(self.policy.export_cap,self.counter+channel.marks[0])
V=self.policy.V;Y40=dot((0,0,40,45,56,72),self.counts)
if self.reactor.enabled:
if self.phase=='pre_entry' and 2500*Y40>=40*V:self.phase='entered';self.entry_time=t
elif self.phase=='entered' and 5000*Y40<=40*V:self.phase='returned';self.status='return_failure'
if not self.corridor():self.status='resource_exit'
self.audit_ledger()
# Apply deterministic gates after a reaction on their exact endpoint.
if self.status=='active' and float(self.next_boundary())==t:self._boundary()
def audit_ledger(self):
V=self.policy.V;E,IU,IW,QF,QP,OA,OB,MW=self.ledger
if dot(A_WEIGHTS,self.counts)+OA!=V+IU or dot(B_WEIGHTS,self.counts)+OB!=V+IW:
raise AssertionError('Moiety ledger failed.')
if dot(M_WEIGHTS,self.counts)+MW!=4*V+2*(IU+IW):raise AssertionError('Mass ledger failed.')
if OA*4<E or OB*4<E:raise AssertionError('Export exceeded washed moiety accounting.')
if self.capped_supplies!=[min(self.policy.supply_cap,n) for n in self.ledger[1:5]]:
raise AssertionError('Capped counters no longer match natural histories.')
def report(self):
V=self.policy.V;T=self.policy.end;IU,IW,QF,QP=self.capped_supplies
supplies=IU<=2*V*T and IW<=2*V*T and QF+QP<=V*T/8
verdict=None
if self.status!='active':
verdict=(self.status=='complete' and self.all_windows_pass and supplies) if self.reactor.enabled else (
self.status=='resource_exit' or (self.status=='complete' and self.all_windows_pass))
return {'status':self.status,'time':self.time,'phase':self.phase,'entry_time':self.entry_time,
'counts':self.counts,'natural_ledger':dict(zip(LEDGER,self.ledger)),
'capped_supply_counts':self.capped_supplies,'current_window_counter':self.counter,
'completed_windows':self.completed,'capped_window_exports':self.window_exports,
'all_windows_pass':bool(self.all_windows_pass),'supplies_pass':bool(supplies),'event_verdict':verdict,
'comparison_scope':'Enabled joint event' if self.reactor.enabled else 'Disabled output schedule OR any resource exit credited as success; no entry/residence/supply requirement.'}
class ProbabilityBudget:
def __init__(self,V=COPY_SCALE,H=POST_STARTUP_DURATION,parameters=Parameters()):
self.V=V;self.H=Q(str(H));self.p=parameters
if not isinstance(V,int) or V<100000000 or self.H<1 or not parameters.theorem_scope():
raise ValueError('Published probability scope: integer V>=1e8, H>=1, K in [8,12], r in [18,22], delta in [1/100,1/20], fixed epsilon.')
def logs(self,ctx=mp):
V=ctx.mpf(self.V);H=number(ctx,self.H);T=500+H;m=self.H.numerator//self.H.denominator
return {'entry':-ctx.mpf(8047)/312500000000*V,'startup_clock':-3*V/500,
'resource_initial':ctx.log(4)-V/1000,
'resource_leak':ctx.log(12000)+ctx.log(V)+ctx.log(T)+number(ctx,'1/50')-V/2000,
'return_initial':-V/62500,
'return_leak':ctx.log(3000)+ctx.log(V)+ctx.log(T)+number(ctx,'18/125')-V/15625,
'all_windows':ctx.log(m)+number(ctx,'1/8')-ctx.mpf(13)*V/1080000,
'two_foods':ctx.log(2)-(2*ctx.log(2)-1)*V*T,'gross_service':-V*T/200}
def interval_pass(self,rho):
rho=Q(str(rho))
if not 0<rho<1:raise ValueError('Failure target must lie between zero and one.')
return bool(logsum(iv,self.logs(iv).values()).b<=iv.log(number(iv,rho)).a)
def evaluate(self):
logs=self.logs();intervals=self.logs(iv);errorlog=logsum(mp,logs.values())
in_simple=bool(iv.log(number(iv,self.H)).b<=(iv.mpf(self.V)/1000000).a)
result={'V':self.V,'H':str(self.H),'complete_windows':self.H.numerator//self.H.denominator,
'component_log10':{k:mp.nstr(v/mp.log(10),25) for k,v in logs.items()},
'failure_formula':mp.nstr(mp.exp(errorlog),35),
'enabled_success_lower_ordinary':mp.nstr(max(mp.mpf(0),1-mp.exp(errorlog)),35),
'log_failure_outward_interval':str(logsum(iv,intervals.values())),
'simplified_horizon_verified':in_simple,
'simplified_success_lower':mp.nstr(1-2*mp.exp(-mp.mpf(self.V)/50000000),35) if in_simple else None,
'log10_exponential_horizon':mp.nstr(mp.mpf(self.V)/1000000/mp.log(10),25),
'disabled_log10_upper':None,'scope':'Outward mpmath interval evaluation is numerical certificate checking, not Lean evidence; no independent-window assumption.'}
if self.H.denominator==1:
off=[-mp.mpf(self.V)*self.H.numerator/25000,logs['resource_initial'],logs['resource_leak']]
offiv=[-iv.mpf(self.V)*self.H.numerator/25000,intervals['resource_initial'],intervals['resource_leak']]
result.update(disabled_log10_upper=mp.nstr(logsum(mp,off)/mp.log(10),35),
disabled_log_upper_interval=str(logsum(iv,offiv)),
disabled_below_10_power_minus_21699=bool(logsum(iv,offiv).b<(-21699*iv.log(10)).a))
return result
@classmethod
def sufficient_scale(cls,H,rho,parameters=Parameters()):
H=Q(str(H));rho=Q(str(rho))
if H<1 or not 0<rho<1:raise ValueError('H>=1 and failure target in (0,1) required.')
simple=int(mp.ceil(max(100000000,1000000*mp.log(number(mp,H)),50000000*mp.log(2/number(mp,rho)))))
for _ in range(8):
if bool(iv.log(number(iv,H)).b<=(iv.mpf(simple)/1000000).a) and bool(iv.log(2/number(iv,rho)).b<=(iv.mpf(simple)/50000000).a):break
simple+=1
else:raise ArithmeticError('Sizing rounding check unresolved.')
candidate=100000000
for _ in range(100):
if cls(candidate,H,parameters).interval_pass(rho):
return {'simple_sufficient_V':simple,'individually_interval_checked_V':candidate,'minimum_claimed':False}
candidate=(11*candidate+9)//10
raise RuntimeError('No checked geometric candidate in 100 steps.')
class ProductionLedger:
def __init__(self,policy,parameters=Parameters()):self.policy=policy;self.p=parameters
def guarantees(self):
V=self.policy.V;T=self.policy.end;m=self.policy.windows
return {'per_window_export_mass':str(Q(V,5000)),'total_completed_export_mass':str(Q(V*m,5000)),
'covalent_equivalents_per_window':str(Q(V,20000)),
'each_food_budget':str(2*V*T),'gross_drive_budget':str(V*T/8),
'resident_mass_interval':[str(Q(18,5)*V),str(Q(22,5)*V)],
'service_per_export_upper':str(625*T/m),'recovery_ratio_lower':str(Q(m,5000)/(4+8*T)),
'scope':'Deterministic consequences conditional on the joint event; not expected physical efficiency.'}
def aligned_interval(self,length):
length=Q(str(length))
if not 0<=length<=self.policy.windows:raise ValueError('Interval must fit within the completed-window span.')
return Q(self.policy.V,5000)*max(0,math.floor(length)-1)
def physical(self,Cref=REFERENCE_MOLAR,dilution=DILUTION_PER_HOUR,temperature=TEMPERATURE_K):
Cref=number(mp,Cref);dilution=number(mp,dilution);temperature=number(mp,temperature)
if min(Cref,dilution,temperature)<=0:raise ValueError('Positive dimensional scales required.')
NA=mp.mpf('6.02214076e23');volume=self.policy.V/(NA*Cref)
affinity=mp.log(number(mp,self.p.K/self.p.eta))
return {'volume_pL':mp.nstr(volume*10**12,25),'startup_hours':mp.nstr(500/dilution,25),
'post_startup_hours':mp.nstr(number(mp,self.policy.H)/dilution,25),
'covalent_equivalent_rate_pM_per_hour':mp.nstr(Cref*dilution/20000*10**12,25),
'reservoir_affinity_kBT':mp.nstr(affinity,25),
'reservoir_affinity_kJ_per_mol':mp.nstr(mp.mpf('8.314462618')*temperature*affinity/1000,25),
'scope':'Unit translation and maintained-reservoir requirement, not calibration; collected complexes are not purified free active product.'}
def algebraic_certificate():
"""Recompute finite source identities and exact worst-box coefficient slack."""
import sympy as sp
counts=sp.symbols('u w x c1 c2 z');V=100000000
reactor=DrivenReactor();off=DrivenReactor(enabled=False)
identities={}
for name,weights in (('A',A_WEIGHTS),('B',B_WEIGHTS)):
drift,noise=reactor.generator(counts,V,weights)
identities[name+'_restoring_drift']=sp.expand(drift-(V-dot(weights,counts)))==0
# Resource quadratic rate is bounded by V+2D because washout marks are 0,1,2.
identities[name+'_quadratic_bound']=all(Q(c.marks[5 if name=='A' else 6])**2<=2*c.marks[5 if name=='A' else 6]
for c in reactor.channels if 12<=c.label<=17)
p=reactor.p;u,w,x,c1,c2,z=counts
drift,noise=reactor.generator(counts,V,Y_WEIGHTS)
expansion=(p.epsilon+p.delta*p.eta)*u*w/V-(p.epsilon/p.K+p.delta)*x+(Q(5,2)*u/V-1)*x+(-Q(29,8)+Q(11,2)*w/V)*c1+Q(11,10)*c2+(p.release/5-Q(9,5)-8/p.K)*z-p.release/(5*V)*(x*x-x)
identities['Y_drift_expansion']=sp.expand(drift-expansion)==0
identities['mass_balance_each_label']=all(dot(M_WEIGHTS,c.jump)+c.marks[7]==2*(c.marks[1]+c.marks[2]) for c in reactor.channels)
identities['internal_moiety_conservation']=all(dot(A_WEIGHTS,c.jump)==dot(B_WEIGHTS,c.jump)==0 for c in reactor.channels if c.label<10 or c.label>=18)
identities['disabled_pathwise_covalent_creation_bound']=all(Q(dot(S_WEIGHTS,c.jump)+c.marks[0],4)<=int(c.label in (0,19)) for c in off.channels)
offdrift,_=off.generator(counts,V,S_WEIGHTS)
export=sum(c.count_rate(counts,V)*c.marks[0] for c in off.channels)
identities['disabled_generator_identity']=sp.expand(offdrift+export-4*((p.epsilon+p.delta*p.eta)*u*w/V-(p.epsilon/p.K+p.delta)*x))==0
drift_coefficients=(1-Q(1,20)-Q(22,5000)-Q(1,8000000),Q(31,40),Q(11,10),Q(4,5))
noise_coefficients=(1+Q(1,20)+Q(11,32)+Q(22,25000)+Q(1,8000000),Q(10374,3200),Q(2669,400),Q(113,25))
weights=Y_WEIGHTS[2:]
drift_slack=[a-Q(39,100)*w for a,w in zip(drift_coefficients,weights)]
noise_slack=[5*w-a for a,w in zip(noise_coefficients,weights)]
cE=399*Q(14,25)*Q(2,25)*Q(1,500000000)-Q(1,2500*40000)
cW=Q(1,27000)-Q(1,40000);lam=Q(121,100)*Q(321,320)*Q(1,500000000)
disabled_slack=Q(1,21000)-1004*lam-Q(1,25000)
if not all(identities.values()) or min(*drift_slack,*noise_slack,disabled_slack)<=0:raise AssertionError('Algebraic certificate failed.')
return {'source_identities':identities,'drift_coefficient_slack':list(map(str,drift_slack)),
'noise_coefficient_slack':list(map(str,noise_slack)),'entry_exponent':str(cE),'window_exponent':str(cW),
'disabled_creation_intensity_coefficient':str(lam),'disabled_chernoff_slack':str(disabled_slack),
'scope':'Source identities and rational coefficient checks, not a replay of the finite-kernel or Lean proof.'}
def main():
parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));parser.add_argument('--simulate',action='store_true')
args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
parameters=Parameters();policy=ObservationPolicy();reactor=DrivenReactor(parameters);off=DrivenReactor(parameters,False)
budget=ProbabilityBudget(policy.V,policy.H,parameters);bounds=budget.evaluate();ledger=ProductionLedger(policy,parameters)
result={'parameters':{n:str(getattr(parameters,n)) for n in parameters.__dataclass_fields__},
'probability':bounds,'requested_failure':TARGET_FAILURE,'requested_failure_interval_pass':budget.interval_pass(TARGET_FAILURE),
'sufficient_scales':budget.sufficient_scale(policy.H,TARGET_FAILURE,parameters),
'event_guarantees':ledger.guarantees(),'physical_reading':ledger.physical(),
'startup_no_seed_probability_lower_at_time_one':mp.nstr(mp.exp(-number(mp,parameters.epsilon+parameters.delta*parameters.eta)*policy.V),25),
'claimed_0_9238_success_interval_check':budget.interval_pass('0.0762'),
'algebraic_certificate':algebraic_certificate()}
result['expected_completed_export_lower_ordinary']=mp.nstr(number(mp,Q(policy.V*policy.windows,5000))*mp.mpf(bounds['enabled_success_lower_ordinary']),35)
if args.simulate:result['limited_count_path']=reactor.stochastic(policy)
if policy.H>10000:raise ValueError('Default illustration requires H<=10000. Use ProbabilityBudget directly for larger durations.')
enabled=reactor.deterministic(policy.H);disabled=off.deterministic(policy.H)
independent=reactor.deterministic(policy.H,method='BDF')
times=np.unique(np.r_[np.linspace(0,20,201),np.linspace(20,float(policy.end),581)])
values=enabled.sol(times).T;offvalues=disabled.sol(times).T
result['deterministic_illustration']={'entry_time':float(enabled.t_events[0][0]) if len(enabled.t_events[0]) else None,
'max_independent_solver_difference':float(np.max(np.abs(values-independent.sol(times).T))),
'scope':'Density-limit ODE with ledgers per V; does not estimate finite-copy startup or event probabilities.'}
windows=[]
for i in range(policy.windows):
end=501+i;start=end-1
windows.append([i,float(enabled.sol(end)[6]-enabled.sol(start)[6]),float(disabled.sol(end)[6]-disabled.sol(start)[6]),1/5000])
def write_csv(name,headers,rows):
with (out/name).open('w',newline='') as f:
writer=csv.writer(f);writer.writerow(headers);writer.writerows(rows)
write_csv('deterministic_trajectory.csv',['time',*SPECIES,*[n+'_per_V' for n in LEDGER]],[[t,*row] for t,row in zip(times,values)])
write_csv('window_exports.csv',['window','enabled_export_per_V','disabled_export_per_V','certified_threshold_per_V'],windows)
sweep=[]
for V in np.linspace(1e8,4e8,81).astype(np.int64):
b=ProbabilityBudget(int(V),policy.H,parameters).evaluate()
sweep.append([int(V),float(b['enabled_success_lower_ordinary']),float(b['simplified_success_lower']) if b['simplified_success_lower'] is not None else '',b['disabled_log10_upper']])
write_csv('copy_scale_bounds.csv',['V','full_success_lower','simplified_success_lower','disabled_log10_upper'],sweep)
write_csv('failure_components.csv',['failure_mode','log10_upper_component'],bounds['component_log10'].items())
(out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
lines=[f'Finite-duration joint success lower bound: {bounds["enabled_success_lower_ordinary"]}.',
f'Simplified bound: {bounds["simplified_success_lower"]}; disabled log10 upper: {bounds["disabled_log10_upper"]}.',
f'Sufficient scales for failure {TARGET_FAILURE}: {result["sufficient_scales"]}.',
f'Output and supply ledger: {result["event_guarantees"]}.',
'Entry and shared failures are paid once; window errors are added without independence.',
'ODE trajectories are illustrative; direct stochastic paths, if requested, may remain unfinished.']
(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.2),layout='constrained');rows=np.array([r[:3] for r in sweep],float)
axes[0].plot(rows[:,0]/1e8,rows[:,1],label='Full finite-duration bound')
axes[0].plot(rows[:,0]/1e8,rows[:,2],'--',label='Simplified bound')
axes[0].set(xlabel='Copy scale V / 1e8',ylabel='Joint success probability lower bound',title='Joint operating-success lower bound');axes[0].legend(fontsize=8)
if policy.H.denominator==1:
axes[1].plot(rows[:,0]/1e8,[float(r[3]) for r in sweep],color='#b95024')
axes[1].set(xlabel='Copy scale V / 1e8',ylabel='log10 disabled probability upper bound',title='Disabled-reactor output upper bound')
else:
axes[1].axis('off');axes[1].text(.1,.5,'Disabled finite-duration theorem shown\nonly for integer H.',transform=axes[1].transAxes)
for ax in axes:ax.grid(alpha=.2)
fig.savefig(out/'reliability.png',dpi=180);fig.savefig(out/'reliability.svg');plt.close(fig)
fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
axes[0].semilogy(times,np.maximum(values[:,:6]@np.array(Y_WEIGHTS,float),1e-14),label='Enabled Y/V')
axes[0].semilogy(times,np.maximum(offvalues[:,:6]@np.array(Y_WEIGHTS,float),1e-14),label='Disabled Y/V')
axes[0].axhline(1/2500,color='black',ls='--',lw=.8,label='Entry threshold')
axes[0].set(xlabel='Dimensionless time',ylabel='Weighted catalytic concentration',title='Catalyst buildup in the deterministic\nreactor');axes[0].legend(fontsize=8)
w=np.array(windows,float)
axes[1].semilogy(w[:,0]+1,w[:,1],label='Enabled')
axes[1].semilogy(w[:,0]+1,w[:,2],label='Disabled')
axes[1].axhline(1/5000,color='black',ls='--',label='Theorem threshold')
axes[1].set(xlabel='Completed unit window after time 500',ylabel='Deterministic covalent export / V',title='Export in consecutive operating windows');axes[1].legend(fontsize=8)
for ax in axes:ax.grid(alpha=.2)
fig.savefig(out/'reactor.png',dpi=180);fig.savefig(out/'reactor.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()
Run output
Finite-duration joint success lower bound: 0.92384922329337936362655799013895908.
Simplified bound: 0.72932943352677461621200101005503119; disabled log10 upper: -21699.858076776522049059661650761192.
Sufficient scales for failure 1/100: {'simple_sufficient_V': 264915869, 'individually_interval_checked_V': 194871710, 'minimum_claimed': False}.
Output and supply ledger: {'per_window_export_mass': '20000', 'total_completed_export_mass': '2000000', 'covalent_equivalents_per_window': '5000', 'each_food_budget': '120000000000', 'gross_drive_budget': '7500000000', 'resident_mass_interval': ['360000000', '440000000'], 'service_per_export_upper': '3750', 'recovery_ratio_lower': '1/240200', 'scope': 'Deterministic consequences conditional on the joint event; not expected physical efficiency.'}.
Entry and shared failures are paid once; window errors are added without independence.
ODE trajectories are illustrative; direct stochastic paths, if requested, may remain unfinished.