Example code
In this sparse polymer model, a few active molecules can catalyse many reactions. The paper proves that a reaction budget proportional to maximum polymer length eventually becomes insufficient for a RAF, even when cleavage is allowed. This example supplies the literal random chemistry and explains the two counting mechanisms behind the result without treating small simulations as an asymptotic proof.
Molecules are binary words of length at most . Each reversible channel is a product word and a split position: . Food contains all words up to a fixed length . PolymerModel builds that catalogue and its reversible food closure; CatalysedNetwork finds maximum RAFs and exactly tests whether any single reachable molecule catalyses a nonempty RAF.
SparseCatalysisLaw activates each molecule once, with probability , then independently gives each active molecule catalytic activity for each channel with probability . Both reaction orientations share that assignment. The distinction matters: one molecule catalyses two channels with probability , not . The example preserves this dependence and reports any clipping of .
The default run uses fixed , and 64 trials at each length from 4 through 8. It finds single-catalyst RAFs in 62, 64, 64, 63 and 63 trials. These finite systems frequently exhibit an event whose probability the theorem says eventually vanishes. The figure reports pointwise Wilson intervals; it does not estimate minimum RAF size or assert an asymptotic trend.
For few catalysts, reaching a catalyst that cannot be made in a few steps requires a longer sequence of molecule-producing reactions. At a fixed number of catalysts, choosing depth leaves an eventual factor. TheoremBounds evaluates the literal constants in logarithmic form. Raw bounds above one are uninformative, and both raw and capped values are saved.
For many catalysts, the key saving is a factorial in the count of productive supports. DependencyCodec labels channels arbitrarily and encodes each ligation or cleavage using food references or earlier-produced endpoints. Acyclic decoding uniquely determines the labelled channels. The manuscript's three-channel example has only one productive firing order, but six distinct labelled codes. This is why the factorial counts labellings, not firing orders.
saturated_program extracts productive channels while retaining the full closure and all reachable catalysts. A food-only RAF can have an empty extraction; the code and probability bound handle that case separately. Exact small enumeration finds 4, 22 and 86 productive supports of sizes 1, 2 and 3 at , , and verifies their labelled-code bound.
The full linear-budget evaluator checks the large-rank premise before using it. For the default parameters and , its conservative catalyst cutoff is 31. At , the assembled raw bound is about : mathematically valid but practically only the trivial bound 1. The example makes that limitation explicit rather than suggesting the paper supplies a sharp finite-size threshold.
Download the package, install requirements.txt, run python example.py --output outputs, and check it with python -m unittest -v. Editable inputs appear first. Import the polymer model, catalysis law, RAF algorithms, dependency codec or bound evaluator independently. Simulation enumerates an exponentially growing catalogue; scalar bound calculations can explore much larger without constructing that catalogue.
Seven test groups compare RAF results with an independent exhaustive oracle, validate reversible closure and empty extraction, audit code injectivity, test shared catalytic activity, and compare the bounds with exact rational formulas. A guarded exact_minimum method is available for small sources. Maximum RAF size is never presented as minimum size. The example does not determine the open scaling law between the linear and quadratic budgets or reproduce the Lean build.
Python source
"""Sparse binary polymer RAFs: literal sampler, construction codes, theorem bounds.
Run python example.py --output outputs. Settings below are dimensionless.
"""
# EDITABLE INPUTS -------------------------------------------------------------
FOOD_HORIZON = 2
INTENSITY = 1.0 # fixed lambda; expected channels/molecule=lambda*n unless clipped
SIMULATION_LENGTHS = (4, 5, 6, 7, 8)
TRIALS_PER_LENGTH = 64
RANDOM_SEED = 20260922
CATALYST_RANKS = (1, 2, 3) # each fixed before n grows
BOUND_LENGTHS = tuple(10 ** k for k in range(1, 19)) # scalar formulas; no molecule enumeration
LINEAR_BUDGET_CONSTANT = 1 # integer C, held fixed in the theorem
# Small exact construction count; independent of simulation's food horizon.
COUNT_MAX_LENGTH = 3
COUNT_FOOD_HORIZON = 1
COUNT_DEPTH = 3
# ---------------------------------------------------------------------------
import argparse
from dataclasses import dataclass
from fractions import Fraction
import hashlib
from itertools import combinations, permutations, product
import json
import math
from pathlib import Path
import platform
import csv
import numpy as np
def molecule_count(n):
return 2 ** (n + 1) - 2 if n >= 0 else 0
def channel_count(n):
return (n - 2) * 2 ** (n + 1) + 4 if n >= 2 else 0
@dataclass(frozen=True, order=True)
class Channel:
"""One reversible channel: product word AND split position identify it."""
word: str
split: int
def __post_init__(self):
if not self.word or set(self.word) - {"0", "1"} or not 1 <= self.split < len(self.word):
raise ValueError("A binary product and a strictly internal split are required")
@property
def endpoints(self):
return self.word[:self.split], self.word[self.split:], self.word
def enabled(self, available):
u, v, w = self.endpoints
return w in available or (u in available and v in available)
def productive(self, available):
return self.enabled(available) and not set(self.endpoints) <= available
@dataclass(frozen=True)
class PolymerModel:
n: int
food_horizon: int
def __post_init__(self):
if not isinstance(self.n, int) or not isinstance(self.food_horizon, int) or self.n < max(2, self.food_horizon) or self.food_horizon < 0:
raise ValueError("Require integer n >= max(2,t) and t >= 0")
def molecules(self):
return tuple("".join(bits) for length in range(1, self.n+1) for bits in product("01", repeat=length))
def food(self):
return frozenset("".join(bits) for length in range(1, self.food_horizon+1) for bits in product("01", repeat=length))
def channels(self):
return tuple(Channel(w, i) for w in self.molecules() for i in range(1, len(w)))
def closure(self, channels):
channels = tuple(channels)
available = set(self.food())
changed = True
while changed:
changed = False
for channel in channels:
if channel.productive(available):
available.update(channel.endpoints)
changed = True
return frozenset(available)
def saturated_program(self, channels):
"""Keep productive firings only, preserving the entire original closure."""
channels = tuple(sorted(channels))
available, sequence = set(self.food()), []
while True:
channel = next((c for c in channels if c.productive(available)), None)
if channel is None:
return ProductiveProgram(self, tuple(sequence))
sequence.append(channel)
available.update(channel.endpoints)
def programs(self, depth):
"""Exhaustive small-depth construction enumeration; exponential."""
catalogue = self.channels()
def visit(sequence, available):
if len(sequence) == depth:
yield ProductiveProgram(self, sequence)
return
for channel in catalogue:
if channel.productive(available):
yield from visit(sequence+(channel,), available | set(channel.endpoints))
yield from visit((), set(self.food()))
@dataclass(frozen=True)
class ProductiveProgram:
model: PolymerModel
sequence: tuple[Channel, ...]
def states(self):
states = [self.model.food()]
for c in self.sequence:
if len(c.word) > self.model.n or not c.productive(states[-1]):
raise ValueError("Each firing must be enabled and add a new molecule")
states.append(states[-1] | frozenset(c.endpoints))
return tuple(states)
@dataclass(frozen=True)
class Reference:
food: str | None = None
label: int | None = None
endpoint: int | None = None
def __post_init__(self):
if self.food is not None:
if self.label is not None or self.endpoint is not None:
raise ValueError("Food and endpoint references are distinct")
elif self.label is None or self.label < 0 or self.endpoint not in (0, 1, 2):
raise ValueError("Endpoint reference needs a label and position")
def resolve(self, assignment, food):
if self.food is not None:
if self.food not in food:
raise ValueError("Reference is not food")
return self.food
return assignment[self.label].endpoints[self.endpoint]
@dataclass(frozen=True)
class Instruction:
kind: str
references: tuple[Reference, ...]
split: int | None = None
class DependencyCodec:
"""Acyclic references encode arbitrary labels, not just firing orders."""
def __init__(self, model):
self.model = model
def encode(self, program, labelled_channels):
if program.model != self.model:
raise ValueError("Program and codec must use the same source model")
program.states() # validate productivity
labels = {c: i for i, c in enumerate(labelled_channels)}
if len(labels) != len(program.sequence) or set(labels) != set(program.sequence):
raise ValueError("Labels must bijectively label the program's support")
refs = {w: Reference(food=w) for w in self.model.food()}
code = [None] * len(labels)
for c in program.sequence:
u, v, w = c.endpoints
label = labels[c]
if u in refs and v in refs:
code[label] = Instruction("ligation", (refs[u], refs[v]))
else:
code[label] = Instruction("cleavage", (refs[w],), c.split)
for endpoint, word in enumerate(c.endpoints):
refs.setdefault(word, Reference(label=label, endpoint=endpoint))
return tuple(code)
def decode(self, code):
assignment = {}
while len(assignment) < len(code):
ready = [i for i, ins in enumerate(code) if i not in assignment and
all(ref.label is None or ref.label in assignment for ref in ins.references)]
if not ready:
raise ValueError("Cyclic or unresolved dependency code")
for i in ready:
ins = code[i]
values = tuple(ref.resolve(assignment, self.model.food()) for ref in ins.references)
if ins.kind == "ligation" and len(values) == 2:
c = Channel(values[0]+values[1], len(values[0]))
elif ins.kind == "cleavage" and len(values) == 1:
c = Channel(values[0], ins.split)
else:
raise ValueError("Malformed instruction")
if len(c.word) > self.model.n:
raise ValueError("Decoded channel exceeds maximum length")
assignment[i] = c
result = tuple(assignment[i] for i in range(len(code)))
if len(set(result)) != len(result):
raise ValueError("Labelled assignments must be injective")
return result
@dataclass(frozen=True)
class CatalysedNetwork:
model: PolymerModel
channels: tuple[Channel, ...]
catalysts: tuple[frozenset[str], ...]
active: frozenset[str]
def __post_init__(self):
if len(self.channels) != len(self.catalysts) or len(set(self.channels)) != len(self.channels):
raise ValueError("Unique channels with one catalyst set each required")
if any(len(c.word) > self.model.n for c in self.channels):
raise ValueError("Channel outside model")
if any(not cs <= self.active for cs in self.catalysts):
raise ValueError("Every assigned catalyst must be active")
if any(not x or len(x) > self.model.n or set(x)-{"0", "1"} for x in self.active):
raise ValueError("Active molecule outside model")
def max_raf(self, allowed=None, catalyst_filter=None):
remaining = set(range(len(self.channels)) if allowed is None else allowed)
if any(i < 0 or i >= len(self.channels) for i in remaining):
raise ValueError("Channel index outside catalogue")
while remaining:
closure = self.model.closure(self.channels[i] for i in remaining)
permitted = closure if catalyst_filter is None else closure & frozenset(catalyst_filter)
keep = {i for i in remaining if self.channels[i].enabled(closure) and self.catalysts[i] & permitted}
if keep == remaining:
break
remaining = keep
return frozenset(remaining)
def single_catalyst_witnesses(self):
"""Exact existence for each active catalyst, including food-only RAFs."""
witnesses = {}
for x in sorted(self.active):
allowed = [i for i, cs in enumerate(self.catalysts) if x in cs]
result = self.max_raf(allowed, (x,))
if result:
witnesses[x] = result
return witnesses
def exact_minimum(self, max_channels=18):
"""Small-model oracle. Returns None iff no RAF; refuses oversized searches.
A maximum RAF or one irreducible RAF is never reported as a minimum.
"""
maximum = tuple(sorted(self.max_raf()))
if len(maximum) > max_channels:
raise ValueError("Exact subset search limit exceeded")
for size in range(1, len(maximum)+1):
for subset in combinations(maximum, size):
if self.max_raf(subset) == frozenset(subset):
return frozenset(subset)
return None
@dataclass(frozen=True)
class SparseCatalysisLaw:
intensity: float
def __post_init__(self):
if not math.isfinite(self.intensity) or self.intensity < 0:
raise ValueError("Finite nonnegative intensity required")
def parameters(self, n):
raw = self.intensity * n*n / channel_count(n)
p, q = min(1.0, raw), 1.0/n
return {"activity": p, "conditional": q, "clipped": raw > 1,
"expected_active": molecule_count(n)*p,
"expected_channels_per_molecule": p*q*channel_count(n)}
def sample(self, model, rng):
molecules, channels = model.molecules(), model.channels()
params = self.parameters(model.n)
# Binomial count + uniform subset is exactly iid Bernoulli membership.
# Activity is sampled ONCE per molecule, shared by both orientations and
# all its channel bits. Never replace it with independent p*q marginals.
count = rng.binomial(len(molecules), params["activity"])
active = frozenset(molecules[int(i)] for i in rng.choice(len(molecules), count, replace=False))
catalysts = [set() for _ in channels]
for word in sorted(active):
count = rng.binomial(len(channels), params["conditional"])
for i in rng.choice(len(channels), count, replace=False):
catalysts[int(i)].add(word)
return CatalysedNetwork(model, channels, tuple(map(frozenset, catalysts)), active)
def log_add(values):
values = tuple(values)
largest = max(values, default=-math.inf)
return largest if largest == -math.inf else largest+math.log(sum(math.exp(x-largest) for x in values))
@dataclass(frozen=True)
class TheoremBounds:
food_horizon: int
intensity: float
def __post_init__(self):
if self.food_horizon < 0 or not isinstance(self.food_horizon, int):
raise ValueError("Nonnegative integer food horizon required")
SparseCatalysisLaw(self.intensity)
@property
def food_size(self):
return molecule_count(self.food_horizon)
def prefix_constants(self, depth):
p, h = 1, 0
f = self.food_size
for j in range(depth):
size = f+3*j
h += p*size
p *= size*size + size*self.food_horizon*2**j
return p, h
def log_parameters(self, n):
"""Stable logs even for huge n; avoid subtracting two ~n*log(2) values."""
if n < max(2, self.food_horizon):
raise ValueError("n must contain food and at least one channel")
if self.intensity == 0:
return -math.inf, -math.inf
small = math.exp(-(n+1)*math.log(2))
divisor = n-2+4*small
log_raw_p = math.log(self.intensity)+2*math.log(n)-(n+1)*math.log(2)-math.log(divisor)
log_p = min(0.0, log_raw_p)
log_np = ((n+1)*math.log(2)+math.log1p(-2*small) if log_raw_p > 0 else
math.log(self.intensity)+2*math.log(n)+math.log1p(-2*small)-math.log(divisor))
return log_p, log_np
def log_rank_bound(self, n, k):
if k < 1:
raise ValueError("Positive catalyst rank required")
if self.food_size == 0 or self.intensity == 0:
return -math.inf
p, h = self.prefix_constants(k+1)
log_p, log_np = self.log_parameters(n)
return log_add((math.log(h)+log_p,
math.log(p)+(k+1)*math.log(k)+k*log_np-(k+1)*math.log(n)))
def log_single_bound(self, n):
if self.food_size == 0 or self.intensity == 0:
return -math.inf
lp, lnp = self.log_parameters(n)
t = self.food_horizon
return log_add((math.log(molecule_count(2*t))+lp,
math.log(channel_count(2*t))+math.log(channel_count(4*t))+lnp-2*math.log(n)))
def cutoff(self, budget):
"""Certified conservative cutoff: replace e by 3 and use integer powers.
For k>=4*C, C*log(D*k)-k*log(4/3) decreases, since log(4/3)>1/4.
Checking its first nonpositive value therefore covers all later ranks.
"""
if not isinstance(budget, int) or budget < 0:
raise ValueError("Use a nonnegative integer linear budget constant")
d = 3*(self.food_size+3)*((self.food_size+3)*budget+1)
k = max(1, 4*budget)
while (d*k)**budget * 3**k > 4**k:
k += 1
return k-1, d
def linear_bound(self, n, budget):
"""Assembled bound ONLY after checking the large-rank scalar premise.
Logs are numerical evaluations of analytic inequalities, not interval
certificates. Raw bounds often greatly exceed 1 at practical sizes.
"""
cutoff, d = self.cutoff(budget)
if budget == 0 or self.food_size == 0 or self.intensity == 0:
return {"applicable": True, "cutoff": cutoff, "log10_raw_bound": None, "capped_bound": 0.0}
lp, _ = self.log_parameters(n)
premise = math.log(self.food_size+3*budget*n)+lp <= n*math.log(9/16)
if not premise:
return {"applicable": False, "cutoff": cutoff, "reason": "Large-rank activity premise not yet satisfied"}
terms = [math.log(self.food_size)+lp]
terms += [self.log_rank_bound(n, k) for k in range(1, cutoff+1)]
terms += [math.log(budget*n)+math.log(self.food_size+3*budget*n)+n*math.log(3/4)]
bound = log_add(terms)
return {"applicable": True, "cutoff": cutoff, "conservative_D": d,
"log10_raw_bound": bound/math.log(10), "capped_bound": math.exp(min(0, bound))}
def log_support_bound(self, n, size):
if size < 1:
return 0.0
refs = self.food_size+3*size
return size*math.log(refs*(refs+n))-math.lgamma(size+1)
def wilson(successes, trials):
z = 1.959963984540054
p = successes/trials
denominator = 1+z*z/trials
center = (p+z*z/(2*trials))/denominator
radius = z*math.sqrt(p*(1-p)/trials+z*z/(4*trials*trials))/denominator
return max(0, center-radius), min(1, center+radius)
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--output", type=Path, default=Path("outputs"))
args = parser.parse_args()
args.output.mkdir(parents=True, exist_ok=True)
if TRIALS_PER_LENGTH < 1:
raise ValueError("At least one trial is required")
law, rng = SparseCatalysisLaw(INTENSITY), np.random.default_rng(RANDOM_SEED)
summaries, trial_rows = [], []
for n in SIMULATION_LENGTHS:
model = PolymerModel(n, FOOD_HORIZON)
raf_count, single_count, active_counts = 0, 0, []
for trial in range(TRIALS_PER_LENGTH):
network = law.sample(model, rng)
maximum = network.max_raf()
single = network.single_catalyst_witnesses()
raf_count += bool(maximum)
single_count += bool(single)
active_counts.append(len(network.active))
trial_rows.append((n, trial, len(network.active), len(maximum), len(single)))
summaries.append({"n": n, **law.parameters(n), "trials": TRIALS_PER_LENGTH,
"raf_count": raf_count, "single_catalyst_count": single_count,
"mean_active": sum(active_counts)/TRIALS_PER_LENGTH,
"raf_wilson95": wilson(raf_count, TRIALS_PER_LENGTH),
"single_wilson95": wilson(single_count, TRIALS_PER_LENGTH)})
with (args.output/"trials.csv").open("w", newline="") as stream:
writer = csv.writer(stream)
writer.writerow(("n", "trial", "active_molecules", "maximum_raf_channels_NOT_minimum", "single_catalyst_witness_count"))
writer.writerows(trial_rows)
bounds = TheoremBounds(FOOD_HORIZON, INTENSITY)
bound_rows = []
for n in BOUND_LENGTHS:
for k in CATALYST_RANKS:
log_bound = bounds.log_rank_bound(n, k)
bound_rows.append((n, k, log_bound/math.log(10), math.exp(min(0, log_bound))))
with (args.output/"rank_bounds.csv").open("w", newline="") as stream:
writer = csv.writer(stream)
writer.writerow(("n", "rank", "log10_raw_upper_bound", "capped_probability_upper_bound"))
writer.writerows(bound_rows)
# Paper's example: two ligations, then a productive cleavage. Only one order
# works, while every one of the 3! arbitrary labelings has a distinct code.
model = PolymerModel(5, 2)
sequence = (Channel("010", 2), Channel("01011", 3), Channel("01011", 4))
program, codec = ProductiveProgram(model, sequence), DependencyCodec(model)
codes = []
for labelled in permutations(sequence):
code = codec.encode(program, labelled)
assert codec.decode(code) == labelled
codes.append({"labels": [vars(c) for c in labelled], "instructions": [
{"kind": ins.kind, "split": ins.split, "references": [vars(r) for r in ins.references]} for ins in code]})
counter = PolymerModel(COUNT_MAX_LENGTH, COUNT_FOOD_HORIZON)
counts = []
for depth in range(1, COUNT_DEPTH+1):
programs = tuple(counter.programs(depth))
supports = {frozenset(p.sequence) for p in programs}
f = len(counter.food())
instruction_count = (f+3*depth)**2+counter.n*(f+3*depth)
assert len(supports)*math.factorial(depth) <= instruction_count**depth
counts.append({"depth": depth, "productive_orders": len(programs), "supports": len(supports),
"labelled_supports": len(supports)*math.factorial(depth),
"instruction_code_bound": instruction_count**depth})
single_log = bounds.log_single_bound(max(BOUND_LENGTHS))/math.log(10)
results = {"simulation": summaries, "exact_small_counts": counts, "paper_example_codes": codes,
"paper_example_states": [sorted(s) for s in program.states()],
"linear_bound_at_n1000": bounds.linear_bound(1000, LINEAR_BUDGET_CONSTANT),
"single_catalyst_log10_bound_at_largest_n": single_log if math.isfinite(single_log) else None,
"zero_bound_log_convention": "null in JSON and -inf in CSV mean an exactly zero bound"}
(args.output/"results.json").write_text(json.dumps(results, indent=2, allow_nan=False)+"\n")
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
plt.rcParams.update({"font.size": 11, "axes.spines.top": False, "axes.spines.right": False,
"svg.hashsalt": "sparse-polymer-v1"})
fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), layout="constrained")
for field, interval, label, color in (("raf_count", "raf_wilson95", "Any RAF", "#187f91"),
("single_catalyst_count", "single_wilson95", "Single-catalyst RAF", "#ae6539")):
means = [s[field]/s["trials"] for s in summaries]
axes[0].errorbar([s["n"] for s in summaries], means,
yerr=[[p-s[interval][0] for s,p in zip(summaries,means)],
[s[interval][1]-p for s,p in zip(summaries,means)]],
fmt="o-", capsize=3, label=label, color=color)
axes[0].set(xlabel="Maximum polymer length n", ylabel="Sample fraction (95% Wilson intervals)",
ylim=(-.04, 1.04), title="Finite samples: existence is not a size theorem")
axes[0].legend(frameon=False)
axes[1].plot([s["n"] for s in summaries], [s["mean_active"] for s in summaries], "o-", label="Sample mean", color="#187f91")
axes[1].plot([s["n"] for s in summaries], [s["expected_active"] for s in summaries], "--", label="Exact law expectation", color="#ae6539")
axes[1].set(xlabel="Maximum polymer length n", ylabel="Number of active molecules", title="Active catalyst counts across polymer\nlengths")
axes[1].legend(frameon=False)
for ext in ("png", "svg"):
fig.savefig(args.output/f"sampling.{ext}", dpi=165)
plt.close(fig)
fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), layout="constrained")
for k in CATALYST_RANKS:
rows = [r for r in bound_rows if r[1] == k]
axes[0].plot([math.log10(r[0]) for r in rows], [r[2] for r in rows], label=f"Rank {k}")
axes[0].axhline(0, color="#555555", lw=1)
axes[0].set(xlabel="log10(n); analytic evaluation, no simulation", ylabel="log10(raw probability upper bound)",
title="Fixed catalyst ranks: eventually a 1/n bound")
axes[0].legend(frameon=False)
sizes = range(1, 41)
axes[1].plot(list(sizes), [math.lgamma(m+1)/math.log(10) for m in sizes], color="#187f91")
axes[1].set(xlabel="Productive support size m", ylabel="log10(m!) removed from the code count",
title="The factorial counts labels, not valid firing orders")
for ext in ("png", "svg"):
fig.savefig(args.output/f"proof_mechanisms.{ext}", dpi=165)
plt.close(fig)
console = "\n".join(["Sparse binary polymer model; reversible split-position channels",
f"Seed {RANDOM_SEED}; {TRIALS_PER_LENGTH} trials per n; fixed lambda={INTENSITY}, t={FOOD_HORIZON}",
"n / RAF samples / single-catalyst samples: " + str([(s["n"], s["raf_count"], s["single_catalyst_count"]) for s in summaries]),
"Paper example: one productive order, six distinct labelled dependency codes.",
"Small exact construction counts: " + str(counts),
"Full linear-budget bound at n=1000: " + str(results["linear_bound_at_n1000"]),
"The simulation does not estimate minimum RAF size or prove asymptotic absence."])
print(console)
(args.output/"console.txt").write_text(console+"\n")
metadata = {"paper_sha256": "3e04e1ca8aaafd81a2695e0b21b7a6c7b726e3016d80c00a355ab9e9c57ccc3a",
"source_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
"python": platform.python_version(), "numpy": np.__version__, "matplotlib": matplotlib.__version__,
"seed": RANDOM_SEED, "inputs": {"food_horizon": FOOD_HORIZON, "intensity": INTENSITY,
"simulation_lengths": SIMULATION_LENGTHS, "trials": TRIALS_PER_LENGTH,
"bound_lengths": BOUND_LENGTHS, "ranks": CATALYST_RANKS, "linear_budget": LINEAR_BUDGET_CONSTANT,
"count_n": COUNT_MAX_LENGTH, "count_food": COUNT_FOOD_HORIZON, "count_depth": COUNT_DEPTH},
"output_sha256": {p.name: hashlib.sha256(p.read_bytes()).hexdigest() for p in sorted(args.output.iterdir())
if p.is_file() and p.name != "run_metadata.json"}}
(args.output/"run_metadata.json").write_text(json.dumps(metadata, indent=2)+"\n")
if __name__ == "__main__":
main()
Run output
Sparse binary polymer model; reversible split-position channels
Seed 20260922; 64 trials per n; fixed lambda=1.0, t=2
n / RAF samples / single-catalyst samples: [(4, 63, 62), (5, 64, 64), (6, 64, 64), (7, 64, 63), (8, 64, 63)]
Paper example: one productive order, six distinct labelled dependency codes.
Small exact construction counts: [{'depth': 1, 'productive_orders': 4, 'supports': 4, 'labelled_supports': 4, 'instruction_code_bound': 40}, {'depth': 2, 'productive_orders': 28, 'supports': 22, 'labelled_supports': 44, 'instruction_code_bound': 7744}, {'depth': 3, 'productive_orders': 224, 'supports': 86, 'labelled_supports': 516, 'instruction_code_bound': 3652264}]
Full linear-budget bound at n=1000: {'applicable': True, 'cutoff': 31, 'conservative_D': 270, 'log10_raw_bound': 258.1902461138213, 'capped_bound': 1.0}
The simulation does not estimate minimum RAF size or prove asymptotic absence.