Example code
In a random catalytic network, merging duplicate reaction descriptions requires updating their catalysis probabilities. A+AA AAA and AA+A AAA merge into one channel. If their catalytic marks were independent with probability , the merged mark has probability . Assigning a fresh probability- mark after merging changes the model.
The paper proves a critical-window law for this quotient catalogue: when the mean channels catalysed per molecule grows like , RAF probability tends to . Here is the probability that reversible food closure grows without bound in an infinite field of independently open channels. The example separates that theorem from exact finite calculations and simulations.
In the four-molecule fixture, at , the exact split RAF probability is . Marking a merged channel whenever either original is marked preserves that value, while assigning probability p to every merged channel give . The code checks all 4,096 split configurations, preserving RAF existence and lifting every successful quotient witness using an actually catalysed representative. Distinct products such as AB and BA remain distinct throughout.
Reusable components construct binary catalogues, merge channels, sample independent catalysis, compute reversible closure and peel unsupported reactions to a terminal RAF witness. One mark serves both reaction directions. Food-only RAFs count, and a catalysed channel whose endpoints cannot be generated is rejected.
The changing catalyst pool is adaptive. The code therefore checks the paper's probability law of the entire pruning history, rather than assuming independence after conditioning on a random pool. All six histories of the quotient fixture have matching exact count polynomials when pools are replaced by fixed prefixes of equal size.
Exact channel enumeration through word length 12 agrees with the paper's count formula. At and , full-pool channel openness is approximately 0.699150, compared with its limit 0.632121. That openness supplies the argument of the survival function; it is not itself a RAF probability.
Finite catalytic simulations retain the 34-channel gateway bound, including food-only RAFs. Separate static-field simulations estimate finite escape probabilities using cap for escape beyond length . The 30 productive gateways bound infinite survival, a different event. Pointwise sampling intervals describe Monte Carlo uncertainty and are not certified survival brackets.
Editable inputs, exact probability tables, saved histories, reusable classes and seven scientific test groups accompany the package. The paper supplies neither a practical numerical value of the infinite survival profile nor a finite-size convergence rate, and these examples make neither claim. They demonstrate structural autocatalysis, without adding a kinetic or thermodynamic interpretation.
Python source
"""Reaction identity, exact OR laws, reversible RAF peeling, and finite-size diagnostics."""
# EDITABLE STUDY INPUTS. Intensities are dimensionless, not kinetic rates.
CATALYTIC_CAP = 6
INTENSITIES = ('0.05', '0.1', '0.25', '0.5', '1')
MONTE_CARLO_TRIALS = 256
RANDOM_SEED = 36092026
CENSUS_ENUMERATION_MAX = 12
CENSUS_FORMULA_MAX = 40
TOY_PROBABILITIES = ('1/3', '1/2')
ESCAPE_CUTOFFS = (2, 3, 4)
ESCAPE_OPENNESS = '1/10'
MAX_EXACT_CATALYTIC_BITS = 16
MAX_EXPLICIT_CAP = 12
from dataclasses import dataclass
from fractions import Fraction as Q
from itertools import product
from collections import Counter
import argparse
import csv
import hashlib
import json
import math
from pathlib import Path
import platform
import numpy as np
MANUSCRIPT_SHA256='99c0aaf9ca053f4e56a901b84a96ce0c67bba99106b4616f3b5bc17474fd6b1b'
def probability(p):
if isinstance(p,float):raise ValueError('Use an integer or rational string for a model probability.')
p=Q(p)
if not 0<=p<=1:raise ValueError('Probability outside [0,1]; no clamp is applied.')
return p
def words(n):return tuple(''.join(w) for k in range(1,n+1) for w in product('AB',repeat=k))
@dataclass(frozen=True,order=True)
class Description:
u:str
v:str
def __post_init__(self):
if not self.u or not self.v:raise ValueError('Nonempty factors required.')
@property
def w(self):return self.u+self.v
def canonical(self):
return Description(self.v,self.u) if self.u+self.v==self.v+self.u and len(self.v)<len(self.u) else self
def label(self):return self.u+' + '+self.v+' <-> '+self.w
class ReversibleNetwork:
def __init__(self,molecules,food,channels):
self.molecules=tuple(molecules);self.channels=tuple(channels);self.index={m:i for i,m in enumerate(self.molecules)}
if len(self.index)!=len(self.molecules) or not set(food)<=set(self.molecules):raise ValueError('Unique molecules and a valid food set required.')
self.food=sum(1<<self.index[x] for x in set(food));self.all_molecules=(1<<len(self.molecules))-1
self.factors=[];self.products=[];self.endpoints=[]
for r in self.channels:
if any(x not in self.index for x in (r.u,r.v,r.w)):raise ValueError('Channel endpoint missing from molecule catalogue.')
factors=(1<<self.index[r.u])|(1<<self.index[r.v]);prod=1<<self.index[r.w]
self.factors.append(factors);self.products.append(prod);self.endpoints.append(factors|prod)
def closure(self,active,seed=None):
H=self.food if seed is None else seed
if H<0 or H&~self.all_molecules:raise ValueError('Invalid seed mask.')
active=tuple(active)
if any(type(j) is not int or not 0<=j<len(self.channels) for j in active):raise ValueError('Unknown channel.')
while True:
new=H
for j in active:
if H&self.products[j] or H&self.factors[j]==self.factors[j]:new|=self.endpoints[j]
if new==H:return H
H=new
def is_raf(self,active,marks):
active=tuple(active);H=self.closure(active)
return bool(active) and all(self.endpoints[j]&H==self.endpoints[j] and marks[j]&H for j in active)
def gateways(self,productive=False):
return tuple(j for j in range(len(self.channels)) if (self.products[j]&self.food or self.factors[j]&self.food==self.factors[j]) and (not productive or not self.products[j]&self.food))
def generated_words(self,mask):return [w for i,w in enumerate(self.molecules) if mask>>i&1]
class ChannelQuotient:
def __init__(self,split):
self.split=split;canonical=sorted(set(r.canonical() for r in split.channels));lookup={r:j for j,r in enumerate(canonical)}
self.mapping=tuple(lookup[r.canonical()] for r in split.channels)
fibres=[[] for _ in canonical]
for i,j in enumerate(self.mapping):fibres[j].append(i)
self.fibres=tuple(tuple(f) for f in fibres)
self.quotient=ReversibleNetwork(split.molecules,split.generated_words(split.food),canonical)
if any(len(f)>2 for f in self.fibres):raise ValueError('Input repeats identical descriptions; this binary-word quotient assumes each ordered pair appears once.')
if len(set(split.channels))!=len(split.channels):raise ValueError('Duplicate records require their own dependence model.')
@classmethod
def binary(cls,n):
if type(n) is not int or not 2<=n<=MAX_EXPLICIT_CAP:raise ValueError('Explicit catalogue cap outside implementation budget.')
molecules=words(n);return cls(ReversibleNetwork(molecules,words(2),[Description(w[:i],w[i:]) for w in molecules for i in range(1,len(w))]))
def project_marks(self,marks):
if len(marks)!=len(self.split.channels):raise ValueError('Wrong split matrix size.')
result=[]
for fibre in self.fibres:
value=0
for i in fibre:value|=marks[i]
result.append(value)
return tuple(result)
def or_parameters(self,p):
p=probability(p);return tuple(1-(1-p)**len(f) for f in self.fibres)
def lift_witness(self,selected,split_marks):
"""Choose representatives using catalytic witnesses, not arbitrary labels."""
projected=self.project_marks(split_marks);selected=tuple(selected)
if not self.quotient.is_raf(selected,projected):raise ValueError('Input is not a quotient RAF for the projected marks.')
H=self.quotient.closure(selected)
lifted=tuple(next(i for i in self.fibres[j] if split_marks[i]&H) for j in selected)
if not self.split.is_raf(lifted,split_marks) or self.split.closure(lifted)!=H:raise ArithmeticError('Witness lift failed.')
return lifted
class CatalyticPeeler:
def __init__(self,network):self.network=network
def run(self,marks,prefix=False):
q=self.network
if len(marks)!=len(q.channels) or any(type(m) is not int or m<0 or m&~q.all_molecules for m in marks):raise ValueError('One valid catalyst mask per reversible channel required.')
active=frozenset(j for j,m in enumerate(marks) if m);history=[tuple(sorted(active))]
for _ in range(len(q.channels)+1):
H=q.closure(active);pool=(1<<H.bit_count())-1 if prefix else H
next_active=frozenset(j for j,m in enumerate(marks) if m&pool)
if not next_active<=active:raise ArithmeticError('Peeling history is not decreasing.')
if next_active==active:
usable=tuple(j for j in sorted(active) if q.endpoints[j]&H==q.endpoints[j])
if not prefix and usable and not q.is_raf(usable,marks):raise ArithmeticError('Terminal RAF failed literal check.')
return {'history':tuple(history),'closure':H,'usable':usable,'has_raf':bool(usable)}
active=next_active;history.append(tuple(sorted(active)))
raise ArithmeticError('Peeling did not terminate within its finite bound.')
class IndependentCatalysis:
def __init__(self,rng):self.rng=rng
def sample(self,network,p):
"""Binomial count plus uniform subset = independent Bernoulli coordinates."""
p=float(probability(p));M=len(network.molecules);counts=self.rng.binomial(M,p,len(network.channels));marks=[]
for count in counts:
mask=0
for i in self.rng.choice(M,size=int(count),replace=False):mask|=1<<int(i)
marks.append(mask)
return tuple(marks)
class ExactCatalyticLaw:
def __init__(self,network):
self.network=network;self.M=len(network.molecules);self.bits=self.M*len(network.channels)
if self.bits>MAX_EXACT_CATALYTIC_BITS:raise ValueError('Exact catalytic enumeration budget exceeded.')
def configurations(self):
for mask in range(1<<self.bits):yield tuple((mask>>(j*self.M))&((1<<self.M)-1) for j in range(len(self.network.channels)))
def probability(self,parameters,event=None):
parameters=tuple(probability(p) for p in parameters)
if len(parameters)!=len(self.network.channels):raise ValueError('One parameter per channel required.')
peeler=CatalyticPeeler(self.network);event=event or (lambda marks:peeler.run(marks)['has_raf']);total=Q(0)
for marks in self.configurations():
if not event(marks):continue
weight=Q(1)
for p,m in zip(parameters,marks):k=m.bit_count();weight*=p**k*(1-p)**(self.M-k)
total+=weight
return total
def history_polynomials(self,prefix=False):
"""Counts by history and number of marked coordinates, valid at every p."""
histories={};peeler=CatalyticPeeler(self.network)
for marks in self.configurations():
history=peeler.run(marks,prefix)['history'];histories.setdefault(history,Counter())[sum(m.bit_count() for m in marks)]+=1
return histories
def history_probability(self,history,p):
"""Exact column-deletion law, including the final stability step."""
p=probability(p);history=tuple(frozenset(A) for A in history)
if not history or any(not A<=set(range(len(self.network.channels))) for A in history):raise ValueError('Invalid history.')
if any(not B<A for A,B in zip(history,history[1:])):return Q(0)
extended=history+(history[-1],)
sizes=[self.M]+[self.network.closure(A).bit_count() for A in history]
s=1-p;weight=Q(1)
for j in range(len(self.network.channels)):
drop=next((t for t,A in enumerate(extended) if j not in A),None)
weight*=1-s**sizes[-1] if drop is None else s**sizes[0] if drop==0 else s**sizes[drop]-s**sizes[drop-1]
return weight
class ChannelCensus:
@staticmethod
def formula(n):
if type(n) is not int or not 2<=n<=10000:raise ValueError('Count cap must be an integer from 2 to 10000.')
primitive=[0]*(n//3+1)
for d in range(1,len(primitive)):primitive[d]=2**d-sum(primitive[e] for e in range(1,d) if d%e==0)
deleted=sum(primitive[d]*((n//d-1)**2//4) for d in range(1,len(primitive)))
M=2**(n+1)-2;R=(n-2)*2**(n+1)+4
return {'n':n,'M':M,'split':R,'quotient':R-deleted,'deleted':deleted,'elementary_loss_bound':n*(2**((n-1)//2+1)-2)}
@staticmethod
def openness(p,k):
p=probability(p)
if type(k) is not int or k<0:raise ValueError('Nonnegative pool size required.')
if not k or not p:return 0.
if p==1:return 1.
return -math.expm1(k*math.log1p(-float(p)))
def seed_cutoff(a0,m):
"""Exact integer verification of the sufficient contour cutoff; no enumeration."""
a0=probability(a0)
if not a0 or type(m) is not int or m<2:raise ValueError('Positive openness and integer m >= 2 required.')
b=a0/2;r=1-b*b;C=2*7**64*81;D=m*(m+1)+1
decay=-math.log1p(-float(b*b))
if decay<=0 or math.log(C*D)/decay>1000000:raise ValueError('Exact seed-cutoff verification exceeds arithmetic budget.')
k=max(1,math.ceil(math.log(C*D)/decay))
holds=lambda j:C*D*pow(r.numerator,j)<=pow(r.denominator,j)
while not holds(k):k+=1
while k>1 and holds(k-1):k-=1
return {'a0':str(a0),'m':m,'k':k,'L':10*k,'exact_integer_inequality':True,'scope':'Sufficient supplied-seed cutoff; not a convergence rate or practical survival estimate.'}
def wilson(successes,trials):
if trials==0:return None
z=1.959963984540054;p=successes/trials;den=1+z*z/trials;center=(p+z*z/(2*trials))/den
radius=z*math.sqrt(p*(1-p)/trials+z*z/(4*trials*trials))/den
return [0. if successes==0 else max(0.,center-radius),1. if successes==trials else min(1.,center+radius)]
def toy():
return ChannelQuotient(ReversibleNetwork(('A','AA','AAA','AAAA'),('A','AA'),(Description('A','AA'),Description('AA','A'),Description('A','AAA'))))
def main():
parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
q=toy();split_law=ExactCatalyticLaw(q.split);quotient_law=ExactCatalyticLaw(q.quotient);toyrows=[]
for pstr in TOY_PROBABILITIES:
p=probability(pstr);split=split_law.probability([p]*3);projected=quotient_law.probability(q.or_parameters(p));homogeneous=quotient_law.probability([p]*2)
if split!=projected or not split>homogeneous:raise ArithmeticError('Toy OR-law comparison failed (use interior probabilities for this comparison).')
toyrows.append({'p':str(p),'split_raf':str(split),'projected_or_raf':str(projected),'homogeneous_quotient_raf':str(homogeneous)})
original=quotient_law.history_polynomials();prefix=quotient_law.history_polynomials(True)
if original!=prefix:raise ArithmeticError('Complete-history law failed.')
for history,counts in original.items():
for pstr in TOY_PROBABILITIES:
p=probability(pstr)
enumerated=sum(count*p**k*(1-p)**(quotient_law.bits-k) for k,count in counts.items())
if enumerated!=quotient_law.history_probability(history,p):raise ArithmeticError('Column-deletion history weights disagree with enumeration.')
lifts=0
for marks in split_law.configurations():
a=CatalyticPeeler(q.split).run(marks);b=CatalyticPeeler(q.quotient).run(q.project_marks(marks))
if a['has_raf']!=b['has_raf']:raise ArithmeticError('Configuration-wise RAF projection failed.')
if b['has_raf']:q.lift_witness(b['usable'],marks);lifts+=1
census=[]
for n in range(2,CENSUS_FORMULA_MAX+1):
row=ChannelCensus.formula(n)
if n<=CENSUS_ENUMERATION_MAX:
catalog=ChannelQuotient.binary(n);net=catalog.quotient
if len(catalog.split.channels)!=row['split'] or len(net.channels)!=row['quotient']:raise ArithmeticError('Count formula disagrees with enumeration.')
row.update(gateways=len(net.gateways()),productive_gateways=len(net.gateways(True)))
if n>=4 and (row['gateways'],row['productive_gateways'])!=(34,30):raise ArithmeticError('Gateway census mismatch.')
census.append(row)
normalization=[]
for row in census:
for lam in INTENSITIES:
p=Q(lam)*row['n']/row['quotient']
if not 0<=p<=1:continue
normalization.append({'n':row['n'],'lambda':lam,'p':str(p),'normalization_ratio':str(Q(row['n']*row['M'],row['quotient'])),
'full_pool_openness':ChannelCensus.openness(p,row['M']),'limiting_openness':-math.expm1(-float(Q(lam)))})
rng=np.random.default_rng(RANDOM_SEED);sampler=IndependentCatalysis(rng);net=ChannelQuotient.binary(CATALYTIC_CAP).quotient;finite=[]
if type(MONTE_CARLO_TRIALS) is not int or not 1<=MONTE_CARLO_TRIALS<=100000:raise ValueError('Trial count outside budget.')
for lam in INTENSITIES:
p=probability(Q(lam)*CATALYTIC_CAP/len(net.channels));hits=gateways=food_only=0
for _ in range(MONTE_CARLO_TRIALS):
marks=sampler.sample(net,p);result=CatalyticPeeler(net).run(marks);gateway=any(marks[j] for j in net.gateways())
if result['has_raf'] and not gateway:raise ArithmeticError('RAF without catalytic gateway.')
hits+=result['has_raf'];gateways+=gateway;food_only+=result['has_raf'] and result['closure']==net.food
finite.append({'n':CATALYTIC_CAP,'lambda':lam,'p':str(p),'trials':MONTE_CARLO_TRIALS,'rafs':hits,'gateway_samples':gateways,'food_only_rafs':food_only,
'raf_estimate':hits/MONTE_CARLO_TRIALS,'raf_wilson95':wilson(hits,MONTE_CARLO_TRIALS),
'gateway_probability_exact_formula_float':ChannelCensus.openness(p,len(net.gateways())*len(net.molecules)),
'conditional_estimate':hits/gateways if gateways else None,'conditional_wilson95':wilson(hits,gateways)})
openness=probability(ESCAPE_OPENNESS);escape=[]
for K in ESCAPE_CUTOFFS:
if type(K) is not int or K<2:raise ValueError('Escape cutoff must be an integer at least two.')
network=ChannelQuotient.binary(2*K).quotient;hits=0
for _ in range(MONTE_CARLO_TRIALS):
active=tuple(int(j) for j in np.flatnonzero(rng.random(len(network.channels))<float(openness)));H=network.closure(active)
hits+=any(len(w)>K for w in network.generated_words(H))
escape.append({'K':K,'cap':2*K,'openness':str(openness),'hits':hits,'trials':MONTE_CARLO_TRIALS,'estimate':hits/MONTE_CARLO_TRIALS,'wilson95':wilson(hits,MONTE_CARLO_TRIALS),
'exact_E2':str(1-(1-openness)**30) if K==2 else None,
'scope':'Estimate of finite escape probability e_K, an upper bound on infinite survival; interval is sampling uncertainty, not a certified survival bracket.'})
result={'toy':toyrows,'projection':{'configurations_checked':1<<split_law.bits,'successful_witness_lifts':lifts},
'history_identity':{'configurations_checked':1<<quotient_law.bits,'distinct_histories':len(original),'exact_count_polynomials_equal':True},
'seed_cutoffs':[seed_cutoff('1/10',10),seed_cutoff('1/2',10)],'finite_catalytic':finite,'finite_escape':escape,
'scope':'The theorem identifies the limit as S_qt(1-exp(-lambda)); these computations do not evaluate that unknown infinite profile or prove a finite-size convergence rate. Lean not rerun.'}
def write_json(name,data):(out/name).write_text(json.dumps(data,indent=2)+'\n')
write_json('results.json',result);write_json('census.json',census);write_json('history_polynomials.json',[{'history':[list(A) for A in history],'counts_by_mark_number':dict(sorted(counts.items()))} for history,counts in sorted(original.items())])
def table(name,header,rows):
with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
table('normalization.csv',list(normalization[0]),[list(r.values()) for r in normalization])
table('finite_catalytic.csv',['n','lambda','trials','rafs','food_only_rafs','estimate','wilson_low','wilson_high','conditional_estimate'],[(r['n'],r['lambda'],r['trials'],r['rafs'],r['food_only_rafs'],r['raf_estimate'],*r['raf_wilson95'],r['conditional_estimate']) for r in finite])
table('toy_probabilities.csv',list(toyrows[0]),[list(r.values()) for r in toyrows])
lines=[f'Toy p={r["p"]}: split={r["split_raf"]}, OR projection={r["projected_or_raf"]}, homogeneous quotient={r["homogeneous_quotient_raf"]}.' for r in toyrows]
lines += [f'All {1<<split_law.bits} split configurations preserve RAF existence; {lifts} witnesses lifted.',f'Complete peeling histories have identical count polynomials across {1<<quotient_law.bits} configurations.',
f'Finite catalytic estimates use n={CATALYTIC_CAP}, {MONTE_CARLO_TRIALS} trials per intensity; they are not the limiting survival curve.',
'34 gateways bound finite RAF existence; 30 productive gateways bound infinite survival. Food-only RAFs are retained.']
(out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained');positions=np.arange(len(toyrows));width=.25
for i,(key,label) in enumerate((('split_raf','Split independent marks'),('projected_or_raf','Quotient with OR law'),('homogeneous_quotient_raf','Homogeneous quotient'))):axs[0].bar(positions+(i-1)*width,[float(Q(r[key])) for r in toyrows],width,label=label)
axs[0].set(xticks=positions,xticklabels=[r['p'] for r in toyrows],xlabel='Catalytic probability p',ylabel='Exact toy RAF probability',ylim=(0,1.08),title='RAF probability under three\nchannel-assignment rules');axs[0].legend(fontsize=8,loc='lower right')
for lam in INTENSITIES:
points=[r for r in normalization if r['lambda']==lam];line,=axs[1].plot([r['n'] for r in points],[r['full_pool_openness'] for r in points],label='lambda='+lam)
axs[1].axhline(-math.expm1(-float(Q(lam))),color=line.get_color(),ls=':',lw=.8)
axs[1].set(xlabel='Maximum word length n',ylabel='Full-pool channel openness (not RAF probability)',title='Channel openness versus maximum polymer\nlength');axs[1].legend(fontsize=8)
for ax in axs:ax.grid(axis='y',alpha=.2)
fig.savefig(out/'identity_and_normalization.png',dpi=180);fig.savefig(out/'identity_and_normalization.svg');plt.close(fig)
fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained');x=[float(Q(r['lambda'])) for r in finite];y=[r['raf_estimate'] for r in finite]
axs[0].errorbar(x,y,yerr=[[v-r['raf_wilson95'][0] for v,r in zip(y,finite)],[r['raf_wilson95'][1]-v for v,r in zip(y,finite)]],fmt='o-',label='Finite RAF estimate, pointwise Wilson 95%')
axs[0].plot(x,[r['gateway_probability_exact_formula_float'] for r in finite],'--',label='Finite catalytic gateway upper bound')
axs[0].set(xlabel='Dimensionless intensity lambda',ylabel='Probability',ylim=(0,1.05),title=f'Actual quotient catalysis at n={CATALYTIC_CAP}');axs[0].legend(fontsize=7)
y=[r['estimate'] for r in escape]
axs[1].errorbar([r['K'] for r in escape],y,yerr=[[v-r['wilson95'][0] for v,r in zip(y,escape)],[r['wilson95'][1]-v for v,r in zip(y,escape)]],fmt='o',label='Finite escape estimate, pointwise Wilson 95%')
axs[1].axhline(float(1-(1-openness)**30),ls='--',color='gray',label='Exact productive-gateway escape E2')
axs[1].set(xlabel='Escape cutoff K (computed with cap 2K)',ylabel='Finite escape probability',ylim=(0,1.05),title=f'Static open channels at a={ESCAPE_OPENNESS}');axs[1].legend(fontsize=7)
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'finite_diagnostics.png',dpi=180);fig.savefig(out/'finite_diagnostics.svg');plt.close(fig)
digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
write_json('run_metadata.json',{'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':digest(Path(__file__)),'python':platform.python_version(),'numpy':np.__version__,'random_seed':RANDOM_SEED,'output_sha256':{p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}})
if __name__=='__main__':main()
Run output
Toy p=1/3: split=505585/531441, OR projection=505585/531441, homogeneous quotient=5137/6561. Toy p=1/2: split=4077/4096, OR projection=4077/4096, homogeneous quotient=239/256. All 4096 split configurations preserve RAF existence; 4077 witnesses lifted. Complete peeling histories have identical count polynomials across 256 configurations. Finite catalytic estimates use n=6, 256 trials per intensity; they are not the limiting survival curve. 34 gateways bound finite RAF existence; 30 productive gateways bound infinite survival. Food-only RAFs are retained.