Example Code
Code and saved outputs: RAF emergence across catalysis levels
This example studies when randomly assigned catalysts allow a polymer reaction network to sustain its own catalytic support, following the manuscript. The complete Python program appears below. The figures and numbers on this page were produced by running that program; readers can inspect them without installing anything.
The paper studies a precise random model: molecules are binary words, reactions join or split them, and the six words of length one or two are always supplied as food. Every molecule independently catalyzes every reaction with the same probability. Each reaction's two directions share its catalyst assignments, and different split positions remain different reaction identities.
A RAF is a nonempty reaction set whose reactants can be made from food and whose catalysts are available in that food-generated pool. Write n for maximum molecule length, f_n for the average number of reactions catalyzed by one molecule, and λ = f_n/n. The paper proves that when this ratio approaches any positive finite value, the probability of a RAF approaches a continuous value strictly between zero and one. The program demonstrates the model, checks finite examples and calculates finite probabilities. The general proof of convergence and continuity remains in the paper.
Read the saved figures
Left: finite systems containing a RAF, with maximum lengths 4, 6, 8 and 10. Right: independently marked reversible-reaction fields that reach beyond lengths 2, 3, 4 or 5. The right-hand probabilities are upper approximations to infinite survival; the plotted frequencies estimate those probabilities. The horizontal axis is enlarged to intensities from 0 to 0.35.
Each point uses 400 samples. Lines join sampled values, and the uncertainty marks are pointwise 95% Wilson intervals. They describe sampling error; they do not bound the difference between a finite probability and the paper's limiting probability. No curve here is labeled as the calculated infinite-size answer.
The finite simulation curve can become very steep. That observation cannot determine whether the limiting probability is exactly zero or one. The paper's theorem addresses that separate question.
This is a finite-size probability map. Columns are separately sampled intensity values, with their numerical values printed beneath them. Color and labels show the observed fraction of systems with a RAF. A zero or a one in a finite sample does not establish an exact probability of zero or one. The CSV tables provide unrounded counts and confidence intervals.
Selected numerical results
The following values are actual output from the saved run. At intensity 0.1 and maximum length 10, the average number of reactions catalyzed per molecule is one. The universe contains 2,046 molecule types and 16,388 reaction identities.
| Maximum length | Catalysis intensity | RAFs found / samples | Observed frequency | 95% interval | |---|---|---|---|---| | 4 | 0.1 | 317 / 400 | 79.25% | 75.01%–82.94% | | 6 | 0.1 | 103 / 400 | 25.75% | 21.71%–30.25% | | 8 | 0.1 | 24 / 400 | 6.00% | 4.06%–8.77% | | 10 | 0.1 | 5 / 400 | 1.25% | 0.54%–2.89% | | 10 | 0.5 | 400 / 400 | 100.00% | 99.05%–100.00% |
The paper provides no finite-size convergence rate. In particular, the first four rows do not establish that the limiting probability at intensity 0.1 is zero. The final row does not establish that the limiting probability at intensity 0.5 is one.
Why failure remains possible
A nonempty RAF must contain at least one reaction that can be used directly from food. The six-by-six table lists all 36 such identities. Four produce dimers already supplied as food; the remaining 32 can enlarge the food closure.
If all 36 reactions lack catalysts anywhere in the system, a RAF cannot exist. The probability of this particular failure event is exactly (1 − p)^(36 U_n), where U_n = 2^(n+1) − 2 counts molecule types, and n ≥ 4. With the paper's scaling this approaches exp(−36λ). Consequently the limiting RAF probability is at most 1 − exp(−36λ).
At intensity 0.03, the limiting probability of this particular failure event is approximately 33.96%, giving a ceiling of approximately 66.04% on the limiting RAF probability. Other ways to fail remain possible, so the ceiling is not the actual probability formula.
A RAF can use only food molecules. The script verifies the example 0 + 0 ↔ 00, catalyzed by 0. The correct necessary condition is a reaction usable from food, rather than production of a non-food molecule.
A catalyst that cannot be made can cause both reactions to disappear
In the full maximum-length-four model, give catalysts to only these two reactions:
| Reaction | Catalyst |
|---|---|
| 0 + 00 ↔ 000 | 1111 |
| 000 + 0 ↔ 0000 | 000 |
The pruning calculation first allows all 30 molecule types to provide catalysis, then restricts that pool to molecules actually obtainable from food using currently catalyzed reactions. It repeats until the pool stops changing.
| Round | Candidate molecule types | Reactions with a candidate catalyst | Molecule types obtainable from food | |---|---:|---:|---:| | 1 | 30 | 2 | 8 | | 2 | 8 | 1 | 6 | | 3 | 6 | 0 | 6 |
Neither reaction can produce 1111. Removing that unavailable catalyst removes support for the first reaction. Without the first reaction, the second cannot obtain 000 from food. No RAF remains.
Now replace the first reaction's catalyst by 0000. Both reactions remain: the reachable pool contains eight molecule types, including 000 and 0000, and each reaction has a catalyst in that pool. The script returns the two reaction identities and independently checks the RAF conditions on this returned set.
Closure is a construction of available molecule types. Its steps are not time steps of chemical kinetics, and RAF membership does not require each catalyst to appear before any use of its reaction.
What the infinite process says
In the separate infinite reaction model, every reaction is independently available with probability a = 1 − exp(−λ). The paper calls the probability of generating arbitrarily long words S(a), and proves that the limiting finite RAF probability is S(1 − exp(−λ)).
To inspect this using finite calculations, the script asks whether the closure reaches a word longer than a boundary B, allowing all intermediate lengths up to 2B. The factor of two is justified by Lemma 3.1: the first crossing is a join of two words each of length at most B. Reaching beyond B therefore has a finite witness within cap 2B.
At intensity 0.1 the saved escape counts are 374, 366, 353 and 350 out of 400 for boundaries 2, 3, 4 and 5. The events are nested for every sampled field. Their probabilities decrease to infinite survival as the boundary grows. The displayed finite frequencies do not certify infinite survival or specify how close the last boundary is to its limit.
Theorem 8.3 gives a stronger conclusion on survival: almost surely, a closure that reaches arbitrarily long words contains every finite binary word. One small piece of the proof is this construction from an already available 000 to target 101:
- Join
000and food word1to make0001. - Join
0001and food word0to make00010. - Join
00010and food word1to make000101. - Split
000101after its third bit to obtain000and101.
All four selected reaction identities have products longer than 000. When their marks have not yet been examined, their probability of all being open is a^4. The code checks the actual reaction identities and target generation. The paper uses carefully chosen previously reached words and a probability argument to obtain the almost-sure conclusion. The finite example alone does not prove it.
A strictly positive bound written without rounding
At a = 1/2, corresponding to intensity λ = ln 2, the program checks an integer inequality used by the paper's growth estimate. It finds k = 458 and L = 4580, satisfying
7 × (2 × 7^64 × 81) × 3^458 ≤ 4^458.
Applying Theorem 7.1 and Proposition 7.2 then gives the bounds
2^(−(4578 × 2^4581 + 5)) ≤ Θ(ln 2) ≤ 1 − 2^(−36).
The lower expression is strictly positive and the upper expression is strictly less than one. They are not estimates of the actual value. The lower bound is far too small for ordinary floating-point display, so the code keeps it symbolic. The integer calculation checks a sufficient condition in the paper's proof; it does not replace that proof.
How the program checks its simplified simulation
The main simulation uses the exact distributional replacement in Lemma 5.1. For each reaction it draws a rank representing the position of the first successful catalyst assignment in an ordered list of candidate molecules. A pool of size k activates that reaction when the rank is at most k. The same ranks are retained throughout pruning.
This saves storing a full molecule-by-reaction Boolean matrix. It preserves the distribution of entire reaction histories, as proved in the paper. It does not construct a literal catalytic witness for a particular matrix. That is why the program includes a separate routine using explicit catalytic assignments for the worked examples.
Two checks make the distinction visible:
- For a small three-species teaching system, all 64 possible catalysis matrices are enumerated. The literal and replacement processes have identical history probabilities, checked separately by the number of present coordinates, so the equality holds for every Bernoulli parameter in that example. Two individual matrices have different histories under the two processes: equality in distribution does not mean equality for each matrix.
- Independent samples from the full maximum-length-six model compare the literal and replacement implementations at intensity 0.2. The two observed frequencies naturally differ because they come from different random samples.
| Implementation | RAFs / samples | 95% interval | |---|---|---| | literal catalysis | 276 / 400 | 64.30%–73.33% | | cardinality replacement | 257 / 400 | 59.44%–68.79% |
The static closure-size condition in Lemma 6.1 is also checked at maximum length 6, intensity 0.2 and candidate-pool size 63. In 400 samples, 53 met that condition and 282 contained a RAF. Every sample meeting the condition contained a RAF; there were zero violations. The general implication is supplied by the lemma.
Diagnostic transcript
This transcript is copied directly from the program's saved output.
CRITICAL WINDOW: ACTUAL SAVED OUTPUT
Seed 20260908; 400 independent samples per displayed point.
PASS: U_5=62, R_5=196; 36 food gateways, of which 32 enlarge closure.
PASS: all 64 small-example matrices give equal history polynomials.
PASS: exact literal witnesses, failed-support example and food-only RAF.
PASS: append-then-split produces target 101 using four fresh reactions.
PASS: common-mark RAF monotonicity and nested static escape events.
PASS: static barrier: 53 barrier successes, 282 RAFs; zero implication violations in 400 samples.
PASS: exact positive-bound condition at a=1/2: k=458, L=4580.
Paper-backed lower bound: Theta(ln 2) >= 2^(-((4580-2)*2^(4580+1)+5)) > 0.
Paper-backed upper bound: Theta(ln 2) <= 1 - 2^(-36) < 1.
FINITE RAF FREQUENCIES (95% pointwise Wilson intervals)
n= 4, lambda=0.1: 317/400 = 0.7925 [0.7501, 0.8294]
n= 4, lambda=0.5: 400/400 = 1.0000 [0.9905, 1.0000]
n= 4, lambda=1: 400/400 = 1.0000 [0.9905, 1.0000]
n= 6, lambda=0.1: 103/400 = 0.2575 [0.2171, 0.3025]
n= 6, lambda=0.5: 400/400 = 1.0000 [0.9905, 1.0000]
n= 6, lambda=1: 400/400 = 1.0000 [0.9905, 1.0000]
n= 8, lambda=0.1: 24/400 = 0.0600 [0.0406, 0.0877]
n= 8, lambda=0.5: 400/400 = 1.0000 [0.9905, 1.0000]
n= 8, lambda=1: 400/400 = 1.0000 [0.9905, 1.0000]
n=10, lambda=0.1: 5/400 = 0.0125 [0.0054, 0.0289]
n=10, lambda=0.5: 400/400 = 1.0000 [0.9905, 1.0000]
n=10, lambda=1: 400/400 = 1.0000 [0.9905, 1.0000]
INDEPENDENT FINITE IMPLEMENTATION COMPARISON (n=6, lambda=0.2)
literal catalysis: 276/400, 95% interval [0.6430, 0.7333]
cardinality replacement: 257/400, 95% interval [0.5944, 0.6879]
Static escape frequencies estimate upper approximations to S(a), not S(a) itself.
The paper proves the infinite-size statements. This run did not compile Lean.
Run or download the program
Download critical_window.py. The complete program follows below, so it can also be read or copied directly from this page.
python -m pip install numpy matplotlib
python critical_window.py --output results --trials 400 --seed 20260908
This run used Python 3.11.0, NumPy 2.4.1 and Matplotlib 3.10.5. The source requires Python 3.10 or later with compatible versions of those packages. The tested environment was the version combination just listed. There is no input file, account, API key or external data requirement.
The seed and separate random streams are recorded in run_metadata.json. Exact random-number reproduction should use the recorded package versions. Increasing --trials reduces sampling uncertainty but does not make a finite cap an infinite system. Memory grows with the sample count and number of reactions.
Download the numerical data: finite_raf.csv, static_escape.csv, literal_comparison.csv, and worked_examples.json. PNG images are included above; SVG versions in the same results directory are available for website publication and resizing.
The source SHA-256 is recorded in the downloadable run metadata.
The displayed source below is the same file that produced this run. The paper reports a Lean verification in Appendix B; no Lean source was compiled in this companion task.
Website reproduction
The source and manuscript are unchanged copies of the supplied files. The figures, CSV tables, worked examples and diagnostic transcript were independently regenerated with seed 20260908 and 400 samples per point. Run metadata records package versions, source hash and sampling rules. The author and repository placeholders in the manuscript are not publication metadata; no author name or repository URL has been invented. The manuscript reports Lean verification; no Lean source or verification manifest was supplied or compiled here.
Python source
#!/usr/bin/env python3
"""Executable companion to A Nontrivial Critical Window for RAF Emergence.
Run: python critical_window.py --output results --trials 400 --seed 20260908
Requires Python 3.10+, NumPy and Matplotlib. No input files or network needed.
The finite RAF simulations use the exact distributional reduction in Lemma 5.1.
They do not approximate the infinite-size limit by a fitted curve. Independent
reaction marks are used separately for the infinite model's finite escape events.
The paper, not a finite simulation, proves positivity, continuity and convergence.
"""
from __future__ import annotations
import argparse
import csv
import hashlib
import json
import math
import platform
from collections import Counter
from dataclasses import dataclass
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
@dataclass
class Polymer:
n: int
def __post_init__(self):
if not 2 <= self.n <= 16:
raise ValueError("Use an integer cap between 2 and 16.")
self.words = [format(j, f"0{length}b")
for length in range(1, self.n + 1)
for j in range(2**length)]
self.index = {w: i for i, w in enumerate(self.words)}
self.reactions = [(w, cut) for w in self.words
for cut in range(1, len(w))]
self.rindex = {r: i for i, r in enumerate(self.reactions)}
self.U, self.R = len(self.words), len(self.reactions)
self.u = np.array([self.index[w[:i]] for w, i in self.reactions])
self.v = np.array([self.index[w[i:]] for w, i in self.reactions])
self.w = np.array([self.index[w] for w, i in self.reactions])
self.length = np.array([len(w) for w in self.words])
assert self.U == 2**(self.n + 1) - 2
assert self.R == (self.n - 2)*2**(self.n + 1) + 4
def closure(self, active):
"""Food closure: join when BOTH factors exist; split when product exists.
Catalysts are not prerequisites during this construction. They determine
the active reaction set before the closure is computed.
"""
u, v, w = self.u[active], self.v[active], self.w[active]
have = self.length <= 2
while True:
before = int(have.sum())
can_join = have[u] & have[v]
can_split = have[w]
have[w[can_join]] = True
have[u[can_split]] = True
have[v[can_split]] = True
if int(have.sum()) == before:
return have
def usable(self, have):
return have[self.w] | (have[self.u] & have[self.v])
def catalysis_probability(model, intensity):
return min(1.0, max(0.0, intensity*model.n/model.R))
def first_catalyst_ranks(uniforms, p, U):
"""One fixed random rank per reaction, reused throughout pruning.
P(rank <= k) = 1-(1-p)**k. Rank U+1 means no catalyst. These are
first successes in an ordered Bernoulli column, not actual molecule labels.
Reusing marks is essential; redrawing each pruning round is incorrect.
"""
if p <= 0:
return np.full(uniforms.shape, U + 1, dtype=np.int64)
if p >= 1:
return np.ones(uniforms.shape, dtype=np.int64)
value = np.floor(np.log1p(-uniforms)/math.log1p(-p)) + 1
return np.minimum(value, U + 1).astype(np.int64)
def canonical_raf(model, ranks):
"""Lemma 5.1: exact finite RAF-existence law using only pool sizes.
A sampled rank history is NOT a witness in a particular literal catalysis
matrix. The separate literal_raf routine returns such witnesses.
"""
k = model.U
while True:
active = ranks <= k
have = model.closure(active)
new_k = int(have.sum())
assert new_k <= k
if new_k == k:
# A food-only RAF counts: six retained molecules need not mean failure.
return bool(np.any(active & model.usable(have))), new_k
k = new_k
def literal_raf(model, catalyst_molecules, catalyst_reactions):
"""Prune an explicit list of (molecule, reaction) catalytic assignments."""
pool = np.ones(model.U, dtype=bool)
history = []
while True:
active = np.zeros(model.R, dtype=bool)
active[catalyst_reactions[pool[catalyst_molecules]]] = True
have = model.closure(active)
history.append({"candidate_molecules": int(pool.sum()),
"active_reactions": int(active.sum()),
"reachable_molecules": int(have.sum())})
assert np.all(~have | pool)
if np.array_equal(pool, have):
retained = np.flatnonzero(active & model.usable(have))
# Recompute from the returned set, then check Definition 2.1 literally.
keep = np.zeros(model.R, dtype=bool)
keep[retained] = True
verified = model.closure(keep)
assert np.array_equal(have, verified)
for r in retained:
assert model.usable(verified)[r]
assert np.any(verified[catalyst_molecules[catalyst_reactions == r]])
return retained, have, history
pool = have
def sample_literal(model, p, rng):
# Conditional on their count, iid Bernoulli successes form a uniform subset.
count = int(rng.binomial(model.U*model.R, p))
positions = rng.choice(model.U*model.R, count, replace=False)
return positions // model.R, positions % model.R
def wilson(successes, trials):
"""Pointwise approximate 95% binomial interval, including 0/N and N/N."""
z = 1.959963984540054
phat = successes/trials
denominator = 1 + z*z/trials
center = (phat + z*z/(2*trials))/denominator
radius = z*math.sqrt(phat*(1-phat)/trials + z*z/(4*trials**2))/denominator
return min(phat, max(0.0, center-radius)), max(phat, min(1.0, center+radius))
def probability_row(successes, trials, **fields):
lo, hi = wilson(successes, trials)
return dict(fields, successes=int(successes), trials=trials,
estimate=successes/trials, ci_low=lo, ci_high=hi)
def exact_history_check():
"""Exhaustively check Lemma 5.1 on a small monotone closure operator.
Teaching test only: three abstract species, food {0}, reactions 0->2, 2->1.
All 64 catalysis matrices are enumerated. Counts agree separately for each
number of present coordinates, proving equality for every Bernoulli p in
THIS small example. This is not an enumeration of the full polymer model.
"""
def phi(active):
result = {0}
if active & 1:
result.add(2)
if active & 2:
result.add(1)
return result
distributions = [Counter(), Counter()]
different_paths = 0
for bits in range(64):
histories = []
for canonical in (False, True):
pool = {0, 1, 2}
trace = []
for _ in range(4):
active = sum(1 << r for r in range(2)
if any(bits & (1 << (2*x+r)) for x in pool))
trace.append(active)
actual = phi(active)
pool = set(range(len(actual))) if canonical else actual
key = (tuple(trace), bits.bit_count())
distributions[int(canonical)][key] += 1
histories.append(tuple(trace))
different_paths += histories[0] != histories[1]
assert distributions[0] == distributions[1]
return {"matrices_checked": 64, "history_polynomials_equal": True,
"individual_matrix_histories_differ": different_paths}
def worked_examples():
model = Polymer(4)
def run(assignments):
xs = np.array([model.index[x] for x, r in assignments], dtype=int)
rs = np.array([model.rindex[r] for x, r in assignments], dtype=int)
retained, have, trace = literal_raf(model, xs, rs)
return {"assignments": [{"catalyst": x, "product": r[0], "cut": r[1]}
for x, r in assignments],
"RAF_exists": bool(len(retained)), "history": trace,
"retained_reactions": [list(model.reactions[r]) for r in retained],
"reachable_words": [model.words[i] for i in np.flatnonzero(have)]}
# All unlisted catalysis coordinates are absent; full n=4 reaction universe.
lost = run([("1111", ("000", 1)), ("000", ("0000", 3))])
mutual = run([("0000", ("000", 1)), ("000", ("0000", 3))])
food_only = run([("0", ("00", 1))])
assert not lost["RAF_exists"] and mutual["RAF_exists"] and food_only["RAF_exists"]
assert [r["reachable_molecules"] for r in lost["history"]] == [8, 6, 6]
assert [r["reachable_molecules"] for r in mutual["history"]] == [8, 8]
# Lemma 8.2: make target 101 from an already generated word 000.
repair_model = Polymer(6)
repairs = [("0001", 3), ("00010", 4), ("000101", 5), ("000101", 3)]
active = np.zeros(repair_model.R, dtype=bool)
for r in [("000", 1)] + repairs:
active[repair_model.rindex[r]] = True
assert all(len(w) > 3 for w, cut in repairs)
assert len(set(repairs)) == 4
assert repair_model.closure(active)[repair_model.index["101"]]
return {"lost_catalyst": lost, "mutual_support": mutual,
"food_only": food_only,
"append_then_split": {"starting_word": "000", "target": "101",
"reactions": repairs, "target_verified": True,
"all_open_probability": "a^4"}}
def exact_positive_bound():
"""Verify the numerical condition of Theorem 7.1 exactly at a=1/2, m=2.
theta <= 1/7 implies theta/(1-theta) <= 1/6. Proposition 7.2 then
gives S(1/2) >= 2^(-(R_L+1)). The contour theorem itself is in the paper.
Keep this tiny positive number symbolic to avoid floating-point underflow.
"""
A = 2*7**64*81
k = 1
while 7*A*3**k > 4**k:
k += 1
assert 7*A*3**k <= 4**k
assert k == 1 or 7*A*3**(k-1) > 4**(k-1)
L = 10*k
return {"a": "1/2", "lambda": "ln(2)", "m": 2, "k": k, "L": L,
"exact_condition_verified": "7 * (2*7^64*81) * 3^k <= 4^k",
"positive_lower_bound": f"2^(-(({L}-2)*2^({L}+1)+5))",
"upper_bound_from_36_gateways": "1 - 2^(-36)",
"meaning": "Bounds use Theorem 7.1 and Proposition 7.2; not a numerical estimate."}
def simulate_finite(intensities, caps, trials, seed):
rows = []
for n in caps:
model = Polymer(n)
rng = np.random.default_rng(np.random.SeedSequence(seed, spawn_key=(1, n)))
marks = rng.random((trials, model.R))
previous = np.zeros(trials, dtype=bool)
for intensity in intensities:
p = catalysis_probability(model, intensity)
ranks = first_catalyst_ranks(marks, p, model.U)
outcomes = np.array([canonical_raf(model, r)[0] for r in ranks])
assert np.all(~previous | outcomes) # common marks preserve inclusion
previous = outcomes
log_failure = -math.inf if p == 1 else 36*model.U*math.log1p(-p)
rows.append(probability_row(int(outcomes.sum()), trials, n=n,
intensity=float(intensity), p=p, molecule_count=model.U,
reaction_count=model.R, finite_gateway_ceiling=-math.expm1(log_failure)))
print(f"Finished finite RAF simulations: n={n}", flush=True)
return rows
def simulate_escape(intensities, boundaries, trials, seed):
"""H_B: produce a word longer than B, allowing intermediate length <=2B.
Lemma 3.1 says this exactly captures escape beyond B in the infinite field.
H_B decreases to infinite survival as B grows. Its probability is an UPPER
bound on S(a); our measured frequency estimates that bound, with sampling error.
"""
largest = Polymer(2*max(boundaries))
rng = np.random.default_rng(np.random.SeedSequence(seed, spawn_key=(2,)))
marks = rng.random((trials, largest.R))
rows = []
for intensity in intensities:
a = -math.expm1(-intensity)
previous = np.ones(trials, dtype=bool)
for boundary in boundaries:
model = Polymer(2*boundary)
outcomes = np.array([np.any(model.closure(u[:model.R] < a)
& (model.length > boundary)) for u in marks])
assert np.all(~outcomes | previous)
previous = outcomes
rows.append(probability_row(int(outcomes.sum()), trials,
boundary=boundary, cap=2*boundary,
intensity=float(intensity), openness=a))
print(f"Finished static escape simulations: lambda={intensity:g}", flush=True)
return rows
def reference_comparison(trials, seed):
model = Polymer(6)
rng = np.random.default_rng(np.random.SeedSequence(seed, spawn_key=(3,)))
successes = [0, 0]
p = catalysis_probability(model, 0.2)
for _ in range(trials):
xs, rs = sample_literal(model, p, rng)
successes[0] += bool(len(literal_raf(model, xs, rs)[0]))
ranks = first_catalyst_ranks(rng.random(model.R), p, model.U)
successes[1] += canonical_raf(model, ranks)[0]
return [probability_row(s, trials, method=method, n=6, intensity=0.2)
for method, s in zip(("literal catalysis", "cardinality replacement"), successes)]
def static_barrier_check(trials, seed):
"""Check Lemma 6.1 in its canonical coupling, using k=U/2 > six food words."""
model = Polymer(6)
k = model.U // 2
assert k > 6
p = catalysis_probability(model, .2)
rng = np.random.default_rng(np.random.SeedSequence(seed, spawn_key=(4,)))
barrier_count = raf_count = 0
for _ in range(trials):
ranks = first_catalyst_ranks(rng.random(model.R), p, model.U)
barrier = int(model.closure(ranks <= k).sum()) >= k
raf, _ = canonical_raf(model, ranks)
assert not barrier or raf
barrier_count += barrier
raf_count += raf
return {"n": 6, "intensity": .2, "retained_pool_size_k": k,
"static_openness": -math.expm1(k*math.log1p(-p)), "trials": trials,
"barrier_successes": barrier_count, "RAF_successes": raf_count,
"violations": 0,
"scope": "Samplewise implication in canonical coupling; law transfers by Lemma 5.1."}
def save_csv(path, rows):
with path.open("w", newline="", encoding="utf-8") as stream:
writer = csv.DictWriter(stream, fieldnames=list(rows[0]))
writer.writeheader()
writer.writerows(rows)
def make_figures(out, finite, escape, examples):
plt.rcParams.update({"font.family": "DejaVu Sans", "font.size": 11,
"axes.spines.top": False, "axes.spines.right": False,
"axes.titleweight": "bold", "axes.labelcolor": "#263444",
"text.color": "#263444", "axes.edgecolor": "#AAB4BF",
"figure.facecolor": "#FAFAF7", "axes.facecolor": "#FAFAF7",
"savefig.facecolor": "#FAFAF7"})
colors = ["#305879", "#147D83", "#B2742B", "#845577"]
def save(fig, name):
fig.savefig(out / f"{name}.png", dpi=180, bbox_inches="tight")
fig.savefig(out / f"{name}.svg", bbox_inches="tight")
plt.close(fig)
fig, axes = plt.subplots(1, 2, figsize=(13.8, 5.6))
for color, n in zip(colors, sorted({r["n"] for r in finite})):
data = [r for r in finite if r["n"] == n and r["intensity"] <= .35]
axes[0].errorbar([r["intensity"] for r in data], [r["estimate"] for r in data],
yerr=[[r["estimate"]-r["ci_low"] for r in data],
[r["ci_high"]-r["estimate"] for r in data]],
label=f"Maximum length {n}", color=color, marker="o", ms=3, lw=1.5, capsize=2)
for color, boundary in zip(colors, sorted({r["boundary"] for r in escape})):
data = [r for r in escape if r["boundary"] == boundary and r["intensity"] <= .35]
axes[1].plot([r["intensity"] for r in data], [r["estimate"] for r in data],
label=f"Reach length > {boundary}", color=color, marker="o", ms=3, lw=1.5)
axes[1].fill_between([r["intensity"] for r in data], [r["ci_low"] for r in data],
[r["ci_high"] for r in data], color=color, alpha=.08)
axes[0].set_title("Does a finite system contain a RAF?", loc="left", pad=14)
axes[1].set_title("Does the open-reaction closure escape?", loc="left", pad=14)
axes[0].set_ylabel("Fraction of samples containing a RAF")
axes[1].set_ylabel("Fraction of samples reaching beyond the boundary")
for ax in axes:
ax.set_xlabel(r"Catalysis intensity $\lambda$ (average catalysis / maximum length)")
ax.set_ylim(-.025, 1.025)
ax.set_xlim(-.005, .36)
ax.grid(axis="y", alpha=.2)
ax.legend(loc="lower right", fontsize=9, frameon=False)
fig.text(.07, .01, "Finite samples with pointwise 95% intervals; view enlarged to 0–0.35. Neither panel plots the unknown limiting curve.", fontsize=10)
fig.subplots_adjust(bottom=.19, wspace=.28)
save(fig, "probability_comparison")
caps = sorted({r["n"] for r in finite})
intensities = sorted({r["intensity"] for r in finite})
grid = np.array([[next(r["estimate"] for r in finite if r["n"] == n
and r["intensity"] == x) for x in intensities] for n in caps])
fig, ax = plt.subplots(figsize=(13, 4.7))
im = ax.imshow(grid, vmin=0, vmax=1, cmap="cividis", aspect="auto", origin="lower")
ax.set_xticks(range(len(intensities)), [f"{x:g}" for x in intensities])
ax.set_yticks(range(len(caps)), caps)
ax.set_xlabel(r"Catalysis intensity $\lambda$ (displayed columns are separate sampled values)")
ax.set_ylabel("Maximum molecule length")
ax.set_title("Finite RAF probability across size and catalysis", loc="left", pad=18)
for i in range(len(caps)):
for j in range(len(intensities)):
value = grid[i, j]
label = "<1%" if 0 < value < .01 else ">99%" if .99 < value < 1 else f"{100*value:.0f}%"
ax.text(j, i, label, ha="center", va="center",
color="white" if grid[i,j] < .48 else "#172530", fontsize=9)
fig.colorbar(im, ax=ax, label="Fraction of samples containing a RAF", pad=.025)
fig.text(.1, .01, "A finite-size probability map, not a proved phase boundary. A displayed 0% or 100% is a sample result.", fontsize=10)
fig.subplots_adjust(bottom=.2)
save(fig, "finite_size_map")
food = ["0", "1", "00", "01", "10", "11"]
fig, axes = plt.subplots(1, 2, figsize=(12.5, 5.6), gridspec_kw={"width_ratios": [1, 1.1]})
grid = np.ones((6, 6)); grid[:2, :2] = 0
axes[0].imshow(grid, cmap=matplotlib.colors.ListedColormap(["#DDD8CB", "#C9E6E1"]), vmin=0, vmax=1)
for i, u in enumerate(food):
for j, v in enumerate(food):
axes[0].text(j, i, f"{u}|{v}", ha="center", va="center", fontsize=10)
axes[0].set_xticks(range(6), food); axes[0].set_yticks(range(6), food)
axes[0].set_xlabel("Second food word"); axes[0].set_ylabel("First food word")
axes[0].set_title("36 reactions usable directly from food", loc="left", pad=15)
x = np.linspace(0, .2, 300)
axes[1].plot(x, -np.expm1(-36*x), color=colors[0], lw=2.5)
axes[1].fill_between(x, -np.expm1(-36*x), 1, color="#E9DACA", alpha=.7)
axes[1].text(.003, .98, "Excluded\nprobabilities", fontsize=10, va="top")
axes[1].text(.075, .32, "The actual limiting probability\ncan be below the ceiling.", fontsize=11)
axes[1].set_title("A proved ceiling on the limiting probability", loc="left", pad=15)
axes[1].set_xlabel(r"Catalysis intensity $\lambda$"); axes[1].set_ylabel("Probability")
axes[1].set_ylim(0, 1.02); axes[1].grid(axis="y", alpha=.2)
axes[1].text(.085, .75, r"$1-e^{-36\lambda}$", fontsize=16)
fig.text(.075, .01, "The bar marks the split. Four beige cells produce food already supplied; the other 32 can enlarge the closure.", fontsize=10)
fig.subplots_adjust(bottom=.17, wspace=.36)
save(fig, "food_boundary")
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--output", type=Path, default=Path("results"))
parser.add_argument("--trials", type=int, default=400)
parser.add_argument("--seed", type=int, default=20260908)
args = parser.parse_args()
if args.trials < 1 or args.seed < 0:
parser.error("trials must be positive and seed must be nonnegative")
args.output.mkdir(parents=True, exist_ok=True)
intensities = [0, .01, .03, .06, .1, .12, .14, .16, .18, .2, .25, .35, .5, 1, 2, 3]
model = Polymer(4)
gateway = (model.length[model.u] <= 2) & (model.length[model.v] <= 2)
assert int(gateway.sum()) == 36
assert int((gateway & (model.length[model.w] > 2)).sum()) == 32
assert (Polymer(5).U, Polymer(5).R) == (62, 196)
assert not canonical_raf(model, np.full(model.R, model.U+1))[0]
assert canonical_raf(model, np.ones(model.R))[0]
history = exact_history_check()
examples = worked_examples()
positive = exact_positive_bound()
finite = simulate_finite(intensities, [4, 6, 8, 10], args.trials, args.seed)
escape = simulate_escape(intensities, [2, 3, 4, 5], args.trials, args.seed)
comparison = reference_comparison(args.trials, args.seed)
barrier = static_barrier_check(args.trials, args.seed)
save_csv(args.output / "finite_raf.csv", finite)
save_csv(args.output / "static_escape.csv", escape)
save_csv(args.output / "literal_comparison.csv", comparison)
(args.output / "worked_examples.json").write_text(json.dumps(examples, indent=2)+"\n")
make_figures(args.output, finite, escape, examples)
metadata = {"paper": "Hordijk_Steel_Critical_Window_arxiv.pdf, September 2026",
"seed": args.seed, "trials_per_point": args.trials,
"python": platform.python_version(), "numpy": np.__version__,
"matplotlib": matplotlib.__version__,
"source_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
"exact_history_test": history, "positive_bound": positive,
"static_barrier_check": barrier,
"intervals": "Pointwise 95% Wilson; neither simultaneous nor finite-size error bounds.",
"randomness": "NumPy default_rng, SeedSequence with documented separate spawn keys.",
"coupling": "Common marks across intensities and static escape boundaries.",
"limits": "No infinite survival, continuity, convergence rate or Lean proof was numerically verified."}
(args.output / "run_metadata.json").write_text(json.dumps(metadata, indent=2)+"\n")
lines = ["CRITICAL WINDOW: ACTUAL SAVED OUTPUT", "",
f"Seed {args.seed}; {args.trials} independent samples per displayed point.",
"PASS: U_5=62, R_5=196; 36 food gateways, of which 32 enlarge closure.",
"PASS: all 64 small-example matrices give equal history polynomials.",
"PASS: exact literal witnesses, failed-support example and food-only RAF.",
"PASS: append-then-split produces target 101 using four fresh reactions.",
"PASS: common-mark RAF monotonicity and nested static escape events.",
f"PASS: static barrier: {barrier['barrier_successes']} barrier successes, "
f"{barrier['RAF_successes']} RAFs; zero implication violations in {args.trials} samples.",
f"PASS: exact positive-bound condition at a=1/2: k={positive['k']}, L={positive['L']}.",
f"Paper-backed lower bound: Theta(ln 2) >= {positive['positive_lower_bound']} > 0.",
"Paper-backed upper bound: Theta(ln 2) <= 1 - 2^(-36) < 1.", "",
"FINITE RAF FREQUENCIES (95% pointwise Wilson intervals)"]
for r in finite:
if r["intensity"] in (.1, .5, 1.0):
lines.append(f"n={r['n']:2d}, lambda={r['intensity']:g}: {r['successes']}/{args.trials} "
f"= {r['estimate']:.4f} [{r['ci_low']:.4f}, {r['ci_high']:.4f}]")
lines += ["", "INDEPENDENT FINITE IMPLEMENTATION COMPARISON (n=6, lambda=0.2)"]
for r in comparison:
lines.append(f"{r['method']}: {r['successes']}/{args.trials}, "
f"95% interval [{r['ci_low']:.4f}, {r['ci_high']:.4f}]")
lines += ["", "Static escape frequencies estimate upper approximations to S(a), not S(a) itself.",
"The paper proves the infinite-size statements. This run did not compile Lean."]
transcript = "\n".join(lines)+"\n"
(args.output / "diagnostics.txt").write_text(transcript)
print(transcript)
if __name__ == "__main__":
main()