Example code
In a small autocatalytic reactor, production must first wait for a catalyst molecule to form. A continuous concentration model alone does not give the probability of meeting that startup deadline. This example models the complete ordered binary-polymer catalogue and compares self-catalysis, food catalysis and catalytic deletion with the same basal chemistry and operating protocol.
The task requires catalytic stock at times 1, 100 and 199, positive export in two consecutive windows, material inventory within allowed limits throughout the run and a gross supply allowance. The stochastic simulator follows one uninterrupted trajectory, retains reverse-reaction losses and checks an exact signed synthesis ledger.
Reusable components separate the random catalytic source, marked environment, reactor, observation counters and reliability bounds. The conditional sampler retains the capped Zipf tail and the degree bias caused by requiring an incidence. It also preserves arbitrary nonfood background rather than planting a single edge into an otherwise empty chemistry.
The code reproduces the 224 catalyst–reaction assignments that support the specified production mechanism, the necessary 99% startup scale of 14,272,222, and a sufficient witness scale near 5 × 10⁴⁹. The enormous gap is traced to one absolute product-noise tolerance. At the small default catalogue, the finite source converse is uninformative; the example does not turn that bound into a success estimate.
Download the package for matched deterministic trajectories, fresh small-count startup records, optional full-mission simulation, source sampling, exact and directed scalar checks, and seven scientific test groups. The guide explains corrections to the paper's written compensator calculations. General probability theorems remain manuscript inputs, numerical illustrations supply no proof premise, and Lean is not rerun.
Python source
"""Two consecutive productive windows, with startup and material accounts.
Binary words and normalized rates are designed mathematical inputs. No chemical
species, clinical use, or physically achievable theorem-scale reactor is implied.
"""
from fractions import Fraction as F
from dataclasses import asdict
from pathlib import Path
import argparse
import csv
import hashlib
import json
import math
import platform
import numpy as np
import mpmath as mp
from scipy.stats import beta
from polymer import Catalogue, Environment, PolymerReactor, FOOD, SELECTED
from source import CappedZipf, source_operation_bounds, screening_trials
from certificates import TwoWindowCertificate, first_birth_ceiling, necessary_startup_scale, deterministic_noise_checks, C, D
# EDITABLE INPUTS ------------------------------------------------------------
MAX_WORD_LENGTH = 4
ZIPF_EXPONENT = '1.5' # >1; default critical exponent 2-2/n at n=4
RELIABILITY_COUNT_SCALE = 5001*10**46 # 5.001e49; theorem evaluation, never SSA
FAILURE_TARGET = '0.01'
TOLERANCE_FRACTION = '0.99' # fraction of each of the six allowed tolerances
STARTUP_REPLICATES = 16 # small diagnostic; does not measure rare success
STARTUP_SETTINGS = (('self',100),('self',200),('food',100),('deleted',100))
SEED = 59092026
EVENT_BUDGET = 100000
FULL_PATH_COUNT_SCALE = 20 # optional --full-path; one uninterrupted 199 run
CONDITIONAL_SOURCE_SEED = 59123
CONCENTRATION_MOLAR = '0.000001' # illustrative 1 micromolar dimensionalization
RESIDENCE_TIME_HOURS = '1'
MANUSCRIPT_SHA256 = '662a47f66e9145138f080525227f6d6da0e9af0e8a13ac22b0d27df83f296851'
# Basal epsilon=2e-9 and mark support {1,1.5,2} are fixed by the theorem model.
# Main ODEs use all-one marks; random-source sampling is a separate experiment.
# ---------------------------------------------------------------------------
def main():
parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'))
parser.add_argument('--full-path',action='store_true',help='Also simulate one small-count uninterrupted mission path.')
args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
catalogue=Catalogue(MAX_WORD_LENGTH);source=CappedZipf(ZIPF_EXPONENT);moments=source.moments(MAX_WORD_LENGTH)
reactors={kind:PolymerReactor(catalogue,Environment.matched(catalogue,kind)) for kind in ['self','food','deleted']}
times=np.unique(np.r_[np.linspace(0,199,700),1,100,199]);end=[int(np.flatnonzero(times==t)[0]) for t in [1,100,199]]
traces={};ode=[];X=len(catalogue.words);selected=catalogue.index['0011']
def table(name,header,rows):
with (out/name).open('w',newline='',encoding='utf-8') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
for kind,reactor in reactors.items():
y=reactor.deterministic(times);traces[kind]=y
mass=y[:,:X]@reactor.length;nonfood=y[:,:X]@reactor.nonfood;Q=y[:,X];J=y[:,X+3]
row=dict(environment=kind,selected_at_observations=y[end,selected].tolist(),window_exports=np.diff(Q[end]).tolist(),
gross_food=float(y[-1,X+1]),gross_monomer_feed=float(y[-1,X+2]),maximum_mass_error=float(max(abs(mass-10))),
maximum_nonfood_ledger_error=float(max(abs(nonfood+Q-J))),recovered_fraction=float(Q[-1]/(10+y[-1,X+2])),
evidence='Concentration ODE illustration; does not establish a finite-count success probability')
ode.append(row)
table(f'ode_{kind}.csv',['time','selected_product','total_mass','nonfood_mass','cumulative_export','food_arrivals','food_monomers','signed_synthesis'],
zip(times,y[:,selected],mass,nonfood,Q,y[:,X+1],y[:,X+2],J))
raw=[];startup=[]
if type(STARTUP_REPLICATES)!=int or STARTUP_REPLICATES<1:raise ValueError('Positive startup replicate count required.')
for kind,V in STARTUP_SETTINGS:
runs=[reactors[kind].simulate(V,1.,SEED+j,EVENT_BUDGET) for j in range(STARTUP_REPLICATES)]
raw.extend(dict(environment=kind,**r) for r in runs)
complete=all(r['completed'] for r in runs)
births=sum(r['first_nonfood'] is not None for r in runs);exits=sum(r['first_mass_exit'] is not None for r in runs)
before=sum(r['first_nonfood'] is not None and (r['first_mass_exit'] is None or r['first_nonfood']<r['first_mass_exit']) for r in runs)
k=sum(r['first_selected'] is not None for r in runs);n=len(runs)
ci=([0. if births==0 else float(beta.ppf(.025,births,n-births+1)),1. if births==n else float(beta.ppf(.975,births+1,n-births))] if complete else None)
startup.append(dict(environment=kind,V=V,runs=n,completed=complete,nonfood_births=births,selected_births=k,births_before_exit=before,corridor_exits=exits,
birth_frequency_interval95_numeric=ci,analytic_mission_ceiling=first_birth_ceiling(V) if kind!='food' else None))
certificate=TwoWindowCertificate(MAX_WORD_LENGTH,RELIABILITY_COUNT_SCALE,F(TOLERANCE_FRACTION));reliability=certificate.evaluate()
scales=certificate.sufficient_scale(F(FAILURE_TARGET));necessary=necessary_startup_scale(F(FAILURE_TARGET));scalar=deterministic_noise_checks()
residual=source.small_cap_residual(MAX_WORD_LENGTH)
source_bounds=None
# Both original finite converse scale conditions are checked independently.
converse_valid=RELIABILITY_COUNT_SCALE>=2*MAX_WORD_LENGTH/D and F(MAX_WORD_LENGTH,RELIABILITY_COUNT_SCALE)<=F(1,200000)
if reliability['failure_upper'] is not None and converse_valid:
from certificates import exp_negative_bounds
tail=8*exp_negative_bounds(2*C*RELIABILITY_COUNT_SCALE/MAX_WORD_LENGTH)[1]
with mp.workdps(60):
residual_upper=min(mp.mpf(1),residual['probability_at_least_two']+mp.mpf(tail.numerator)/tail.denominator)
error=mp.mpf(reliability['failure_upper'].numerator)/reliability['failure_upper'].denominator
source_bounds=source_operation_bounds(moments,error,residual_upper)
source_bounds['independent_fresh_experiments_sufficient_numeric']=screening_trials(source_bounds['success_lower'])
source_bounds['finite_converse_informative']=source_bounds['success_upper']<1
incidence=[]
for n in [4,8,16,32,64]:
for label,a in [('a=1.5','1.5'),('a=2','2'),('a=2.5','2.5'),('critical',str(2-2/n))]:
m=CappedZipf(a).moments(n);incidence.append(dict(label=label,n=n,incidence=float(m['incidence']),scaled_incidence=float(m['scaled_incidence']),same_row_pair_ratio=float(m['pair_ratio'])))
conditional=source.sample_environment(catalogue,CONDITIONAL_SOURCE_SEED,True)
selected_degree=sum(z=='0011' for z,r,H in conditional.incidences)
# Preserve the whole sampled environment for reuse, not just its planted edge.
conditional_record=dict(seed=CONDITIONAL_SOURCE_SEED,law='Six food rows empty and selected self-incidence required; correct size-biased row, other nonfood rows unrestricted',
witness_class=conditional.witness_class(catalogue),selected_row_degree=selected_degree,environment=asdict(conditional),
numerical_sampling='Analytic conditional degree law evaluated in floating point; uniform subsets sampled without replacement')
full=reactors['self'].simulate(FULL_PATH_COUNT_SCALE,199.,SEED,EVENT_BUDGET) if args.full_path else None
concentration=F(CONCENTRATION_MOLAR);residence=F(RESIDENCE_TIME_HOURS)
if concentration<=0 or residence<=0:raise ValueError('Positive dimensionalization scales required.')
avogadro=F(602214076)*10**15
dimensions=dict(concentration_molar=concentration,residence_hours=residence,mission_hours=199*residence,
necessary_startup_litres=F(necessary['necessary_integer_scale'])/(avogadro*concentration),
sufficient_witness_litres=F(scales['count_scale'])/(avogadro*concentration),
note='Illustrative units only; equal dimensionless forward/reverse coefficients have different physical molecularity units')
results=dict(inputs=dict(n=MAX_WORD_LENGTH,exponent=ZIPF_EXPONENT,theorem_V=RELIABILITY_COUNT_SCALE,startup_replicates=STARTUP_REPLICATES,seed=SEED),
catalogue=dict(molecules=X,ordered_splits=len(catalogue.splits),productive_labels=len(catalogue.productive_labels())),
deterministic=ode,startup=startup,noise_certificate=reliability,sufficient_scale=scales,uniform_certificate=certificate.uniform(),necessary_startup=necessary,
scalar_checks=scalar,source_moments=moments,small_cap_residual=residual,source_operation=source_bounds,
dimensionalization=dimensions,optional_full_path=full,
asymptotic_constants=dict(gamma=float((6/mp.pi**2)**6),posterior_liminf=float((6/mp.pi**2)**6/224),critical_incidence_limit=float(9/(2*mp.pi**2))),
evidence='Fresh exact/directed scalar checks and numeric source evaluations; imported trajectory and asymptotic theorems; actual ODE/SSA diagnostics; Lean not rerun')
def serializer(v):
if isinstance(v,F):return str(v)
if isinstance(v,mp.mpf):return mp.nstr(v,40)
return float(v)
def dump(name,obj):(out/name).write_text(json.dumps(obj,indent=2,default=serializer,allow_nan=False)+'\n',encoding='utf-8')
dump('results.json',results);dump('startup_paths.json',raw);dump('conditional_environment.json',conditional_record)
table('productive_census.csv',['catalyst','left','right','product'],[(z,*r) for z,r in catalogue.productive_labels()])
table('incidence.csv',list(incidence[0]),[list(r.values()) for r in incidence])
lines=[f'Catalogue: {X} molecules, {len(catalogue.splits)} ordered splits, {len(catalogue.productive_labels())} productive incidences.',
*[f'{r["environment"]} ODE exports: {r["window_exports"][0]:.8g}, {r["window_exports"][1]:.8g}; concentration illustration only.' for r in ode],
*[f'{r["environment"]}, V={r["V"]}: {r["nonfood_births"]}/{r["runs"]} startup births; {r["corridor_exits"]} corridor exits; complete={r["completed"]}.' for r in startup],
f'Necessary integer startup scale for the chosen reliability: {necessary["necessary_integer_scale"]:,}.',
f'Sufficient witness count scale: {scales["count_scale"]:.5e}; dominant term: {scales["dominant_term"]}.',
'The finite source converse is uninformative at n=4. No empirical ensemble success probability is estimated.',
'Scalar checks correct two compensator/forcing slips in the written derivation; theorem sources remain unchanged.']
(out/'console.txt').write_text('\n'.join(lines)+'\n',encoding='utf-8');print('\n'.join(lines))
plot(out,times,traces,selected,X,incidence,startup)
digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest();here=Path(__file__).parent
dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,python=platform.python_version(),source_sha256=digest(Path(__file__)),
module_sha256={n:digest(here/n) for n in ['polymer.py','source.py','certificates.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,times,traces,selected,X,incidence,startup):
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
for kind,y in traces.items():
axs[0].semilogy(times[1:],np.maximum(y[1:,selected],1e-30),label=kind);axs[1].plot(times,y[:,X],label=kind)
for ax in axs:
ax.axvline(100,color='gray',ls=':');ax.set_xlabel('Model time; one uninterrupted trajectory');ax.grid(alpha=.2);ax.legend(fontsize=8)
axs[0].set(title='Concentration trajectories: initiation differs',ylabel='Selected product concentration',ylim=(1e-10,1))
axs[1].set(title='Cumulative aggregate nonfood export',ylabel='Monomer equivalents per count scale')
fig.savefig(out/'trajectories.png',dpi=180);fig.savefig(out/'trajectories.svg');plt.close(fig)
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained');volumes=np.geomspace(10,1e8,200)
axs[0].loglog(volumes,-np.expm1(-(484/3)*2e-9*volumes),label='Food-silent mission ceiling')
axs[0].axhline(.99,color='gray',ls=':',label='99% reliability target');axs[0].set(xlabel='Count scale V',ylabel='Upper bound on success probability',title='Success ceiling from the first-catalyst\ndeadline');axs[0].legend(fontsize=8)
for label in ['a=1.5','a=2','a=2.5','critical']:
rows=[r for r in incidence if r['label']==label];axs[1].plot([r['n'] for r in rows],[r['same_row_pair_ratio'] for r in rows],'o-',label=label)
axs[1].axhline(2/3,color='gray',ls=':',label='a=1.5 limit');axs[1].set(xlabel='Maximum polymer length n',ylabel='P(second incidence | first, same row)',title='Conditional catalytic-assignment\nprobability');axs[1].legend(fontsize=8)
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'startup_source.png',dpi=180);fig.savefig(out/'startup_source.svg');plt.close(fig)
if __name__=='__main__':main()
Run output
Catalogue: 30 molecules, 68 ordered splits, 224 productive incidences. self ODE exports: 102.78205, 109.00104; concentration illustration only. food ODE exports: 136.27595, 136.27595; concentration illustration only. deleted ODE exports: 2.2093594e-05, 2.2175999e-05; concentration illustration only. self, V=100: 0/16 startup births; 0 corridor exits; complete=True. self, V=200: 0/16 startup births; 0 corridor exits; complete=True. food, V=100: 16/16 startup births; 1 corridor exits; complete=True. deleted, V=100: 0/16 startup births; 0 corridor exits; complete=True. Necessary integer startup scale for the chosen reliability: 14,272,222. Sufficient witness count scale: 5.00016e+49; dominant term: selected_product. The finite source converse is uninformative at n=4. No empirical ensemble success probability is estimated. Scalar checks correct two compensator/forcing slips in the written derivation; theorem sources remain unchanged.