Example code
A RAF is a collection of reactions that can obtain their reactants from food and their catalysts from the resulting molecular closure. The example begins with two food-ready reactions, and , that catalyse each other. They form an elementary core. Two outside reactions, and , supply additional catalysts to that core.
The paper proves that a viable elementary core contains a reaction present in at least half the full RAF family, including the empty set. Here frequency means the fraction of allowed reaction sets containing a reaction, not how often it fires. The outside reactions may change which core reaction is most frequent. The code makes that distinction visible by fixing each exterior selection and counting its compatible core subsets.
There are 12 family members in the default system. Both and occur in 6, so each is globally abundant. With exterior selected, their conditional frequencies are and ; with selected, those values reverse. No single core reaction is abundant in both contexts. The theorem establishes this overall abundance by counting total core membership, without assuming that the conditional witness stays fixed.
The package reconstructs each context's support conditions and executes the proof's horizontal/vertical counting injection. It also accepts nonnegative exterior weights: the default reweighting gives core frequencies and , still guaranteeing some core witness. Every compatible row within a context must receive the same weight. These are structural counts, with no kinetic probability attached.
A second construction turns an explicit union-closed family into a producer–gate reaction system. Each coordinate gets two mutually catalysing reactions; marker reactants impose the family's membership conditions. The default benchmark has 8 reactions, 14 markers and 256 checked availability queries, including incomplete pairs. Its closure takes at most two growth rounds and it has no elementary RAF. This illustrates why sparse catalysis alone does not extend the abundance theorem to all systems: arbitrary union-closed families can be encoded in the reactant requirements.
Download the code for reusable reaction systems, closure and maxRAF queries, core and context analysis, projection-failure witnesses, shortest-cycle certificates, and both producer–gate and exact-ground marker constructions. Seven scientific test groups check the manuscript examples and exhaustive small cases. Polynomial query methods remain available independently of the exponential enumeration tools. The package does not rerun the Lean proofs or establish the unrestricted half-frequency conjecture.
Python source
"""Elementary cores force abundance, but arbitrary RAF abundance remains open.
Run python example.py --output outputs. Set-based RAF semantics have no kinetics.
"""
# EDITABLE INPUTS: the paper's changing-witness system, Table 1.
FOOD = ('f',)
REACTIONS = {
'a': {'reactants':('f',), 'products':('x',), 'catalysts':('y',)},
'b': {'reactants':('f',), 'products':('y',), 'catalysts':('x',)},
'g': {'reactants':('f',), 'products':('y',), 'catalysts':('f',)},
'h': {'reactants':('f',), 'products':('x',), 'catalysts':('f',)},
}
CORE = ('a','b')
MODULE = ('a','b')
# The same nonnegative weight applies to every compatible row in a context.
CONTEXT_WEIGHTS = {(): '1', ('g',): '3', ('h',): '1', ('g','h'): '1'}
UNION_GROUND = ('a','b','c','dead')
UNION_MEMBERS = ((), ('a','b'), ('a','c'), ('b','c'), ('a','b','c'))
ENUMERATION_LIMIT = 262144 # explicit rejection, never a silently partial count
MANUSCRIPT_SHA256 = '5d7518980a46f1080e39ca4631539a82ba2c8d91b4362c4b778535ef5fd1fd50'
import argparse
import csv
import hashlib
import json
import platform
from collections import deque
from dataclasses import dataclass
from fractions import Fraction as Q
from itertools import combinations
from pathlib import Path
def subsets(items, limit=ENUMERATION_LIMIT):
items=tuple(sorted(items))
if 2**len(items)>limit:
raise ValueError(f'Enumeration needs {2**len(items)} subsets; limit is {limit}.')
for size in range(len(items)+1):
for row in combinations(items,size):yield frozenset(row)
def rows_json(rows):return [sorted(s) for s in sorted(rows,key=lambda s:(len(s),tuple(sorted(s))))]
@dataclass(frozen=True)
class Reaction:
reactants: frozenset[str]
products: frozenset[str]
catalysts: frozenset[str]
def __post_init__(self):
for key in ('reactants','products','catalysts'):
object.__setattr__(self,key,frozenset(getattr(self,key)))
class CatalyticSystem:
"""Literal closure first; food generation and catalysis checked afterwards."""
def __init__(self,food,reactions):
self.food=frozenset(food)
self.reactions={k:v if isinstance(v,Reaction) else Reaction(**v) for k,v in reactions.items()}
self.names=frozenset(self.reactions)
def selected(self,selected):
selected=frozenset(selected)
if not selected<=self.names:raise ValueError('Unknown reaction selected.')
return selected
def closure_stages(self,selected):
selected=self.selected(selected);stages=[self.food]
while True:
current=stages[-1]
following=current.union(*(self.reactions[r].products for r in selected
if self.reactions[r].reactants<=current))
if following==current:return tuple(stages)
stages.append(following)
def closure(self,selected):return self.closure_stages(selected)[-1]
def food_generated(self,selected):
selected=self.selected(selected);available=self.closure(selected)
return all(self.reactions[r].reactants<=available for r in selected)
def fixed(self,selected):
selected=self.selected(selected);available=self.closure(selected)
return all(self.reactions[r].reactants<=available and bool(self.reactions[r].catalysts&available)
for r in selected)
def is_raf(self,selected):return bool(selected) and self.fixed(selected)
def max_raf(self,availability=None):
current=self.names if availability is None else self.selected(availability)
while True:
available=self.closure(current)
following=frozenset(r for r in current if self.reactions[r].reactants<=available
and self.reactions[r].catalysts&available)
if following==current:return current
current=following
def family(self):return tuple(s for s in subsets(self.names) if self.fixed(s))
@property
def food_ready(self):return frozenset(r for r in self.names if self.reactions[r].reactants<=self.food)
def product_union(self,selected):
return frozenset().union(*(self.reactions[r].products for r in selected))
def predecessors(self,head,module):
return frozenset(r for r in module if self.reactions[r].products&self.reactions[head].catalysts)
def as_json(self):
return {'food':sorted(self.food),'reactions':{r:{k:sorted(getattr(self.reactions[r],k))
for k in ('reactants','products','catalysts')} for r in sorted(self.names)}}
class SupportedFamily:
"""Executable horizontal/vertical cell injection from the counting proof."""
def __init__(self,ground,bodies,upward_predicate):
self.ground=frozenset(ground);self.bodies={r:frozenset(bodies[r]) for r in self.ground}
if any(not b or not b<=self.ground for b in self.bodies.values()):
raise ValueError('Every support body must be a nonempty subset of the ground.')
self.K=upward_predicate
def contains(self,S):
return self.K(S) and all(S&self.bodies[r] for r in S)
def audit_injection(self):
cube=tuple(subsets(self.ground));truth={s:bool(self.K(s)) for s in cube}
if any(truth[s] and not truth[s|{r}] for s in cube for r in self.ground-s):
raise ValueError('The additional constraint is not upward.')
bad={D for D in cube if not self.contains(self.ground-D)}
choices={r:min(self.bodies[r]) for r in self.ground}
records=[];targets=set()
for D in cube:
if D not in bad:continue
failed={r for r in self.ground-D if self.bodies[r]<=D}
for r in sorted(self.ground-D):
horizontal=not truth[self.ground-D] or bool(failed-{r})
target=(D|{r},r) if horizontal else (D,choices[r])
assert target[0] in bad and target[1] in target[0] and target not in targets
targets.add(target)
records.append({'source_deleted':sorted(D),'source_column':r,
'move':'horizontal' if horizontal else 'vertical',
'target_deleted':sorted(target[0]),'target_column':target[1]})
family=tuple(S for S in cube if self.contains(S))
zeros=sum(len(self.ground-D) for D in bad);ones=sum(len(D) for D in bad)
slack=sum(2*len(S)-len(self.ground) for S in family)
assert len(records)==zeros<=ones and slack==ones-zeros>=0
return {'family':rows_json(family),'slack':slack,'bad_zero_cells':zeros,
'bad_one_cells':ones,'injection':records}
class CoreAnalysis:
def __init__(self,system,core):
self.system=system;self.core=system.selected(core)
if not self.core<=system.food_ready or not system.is_raf(self.core):
raise ValueError('The designated core must itself be a nonempty food-ready RAF.')
def fibre(self,exterior):
q=self.system;U=self.core;T=q.selected(exterior)
if T&U:raise ValueError('Exterior context must be disjoint from the core.')
automatic=q.food|q.product_union(T)
bodies={r:frozenset({r}) if q.reactions[r].catalysts&automatic else q.predecessors(r,U) for r in U}
def constraint(S):
W=S|T
# Exterior products are NOT usable without this food-generation gate.
return q.food_generated(W) and all(q.reactions[t].catalysts&(q.food|q.product_union(W)) for t in T)
supported=SupportedFamily(U,bodies,constraint)
audit=supported.audit_injection()
direct=tuple(S for S in subsets(U) if q.fixed(S|T))
assert audit['family']==rows_json(direct)
return {'exterior':sorted(T),'bodies':{r:sorted(bodies[r]) for r in sorted(U)},
**audit,'counts':{r:sum(r in S for S in direct) for r in sorted(U)}}
def report(self,weights=None):
U=self.core;records=[self.fibre(T) for T in subsets(self.system.names-U)]
family=self.system.family();N=len(family)
counts={r:sum(r in W for W in family) for r in sorted(self.system.names)}
assert 2*sum(counts[r] for r in U)-len(U)*N==sum(r['slack'] for r in records)
supplied={} if weights is None else {frozenset(k):Q(v) for k,v in weights.items()}
if any(not k<=self.system.names-U or v<0 for k,v in supplied.items()):
raise ValueError('Context weights must be nonnegative and refer only to exterior reactions.')
Z=Q(0);weighted={r:Q(0) for r in sorted(U)}
for record in records:
h=Q(1) if weights is None else supplied.get(frozenset(record['exterior']),Q(0))
Z+=h*len(record['family'])
for r in U:weighted[r]+=h*record['counts'][r]
assert 2*sum(weighted.values())>=len(U)*Z
return {'family':rows_json(family),'N_including_empty':N,'counts':counts,
'abundant':[r for r in sorted(counts) if 2*counts[r]>=N], 'fibres':records,
'weighted_partition':str(Z),'weighted_core_counts':{r:str(v) for r,v in weighted.items()},
'weighted_core_frequencies':{r:str(v/Z) for r,v in weighted.items()} if Z else None,
'weighted_abundant':[r for r,v in weighted.items() if 2*v>=Z] if Z else []}
class ModuleDiagnostics:
def __init__(self,system,module):
self.system=system;self.module=system.selected(module)
if not self.module<=system.food_ready:raise ValueError('Module diagnostics require food-ready reactions.')
def projection_queries(self):
q=self.system;E=self.module;failures=[]
for r in sorted(E):
if q.reactions[r].catalysts&q.food:continue
W=q.max_raf(q.names-q.predecessors(r,E))
if r in W:
assert q.is_raf(W) and not q.fixed(W&E)
failures.append({'head':r,'witness':sorted(W),'projection':sorted(W&E)})
return {'preserves_projection':not failures,'failures':failures}
def defect(self):
q=self.system;E=self.module;U=q.max_raf(E)
if not U:return {'status':'no_viable_elementary_core'}
family=q.family();bad=[W for W in family if not q.fixed(W&E)]
lhs=2*sum(len(W&U) for W in family)
rhs=len(U)*(len(family)-len(bad))+2*sum(len(W&U) for W in bad)
assert lhs>=rhs
return {'core':sorted(U),'bad_rows':rows_json(bad),'lhs':lhs,'rhs':rhs}
class FoodReadyGraph:
def __init__(self,system):
self.system=system;self.vertices=system.food_ready
self.edges={p:frozenset(r for r in self.vertices if system.reactions[p].products&system.reactions[r].catalysts
or (p==r and system.food&system.reactions[r].catalysts)) for p in self.vertices}
def shortest_cycle(self,available=None):
available=self.vertices if available is None else frozenset(available)&self.vertices
candidates=[]
for start in sorted(available):
queue=deque([(start,(start,))]);seen={start}
while queue:
vertex,path=queue.popleft()
if start in self.edges[vertex]:
candidates.append(path);break
for following in sorted(self.edges[vertex]&available-seen):
seen.add(following);queue.append((following,path+(following,)))
return min(candidates,key=lambda c:(len(c),c)) if candidates else None
def greedy_packing(self):
remaining=self.vertices;cycles=[]
while cycle:=self.shortest_cycle(remaining):
assert self.system.is_raf(cycle)
cycles.append(cycle);remaining=remaining-set(cycle)
return cycles
def abundance_candidates(self):
"""Polynomial certificate locating a witness without counting the family."""
cycle=self.shortest_cycle()
if cycle:return {'reason':'viable_elementary_cycle','candidates':list(cycle)}
if not self.system.max_raf():return {'reason':'no_nonempty_raf','candidates':[]}
exceptional=self.system.names-self.vertices
if len(exceptional)<=2:
return {'reason':'at_most_two_non_food_ready_reactions','candidates':sorted(exceptional)}
return {'reason':'outside_proved_structural_cases','candidates':None}
class UnionClosedFamily:
def __init__(self,ground,members):
self.ground=frozenset(ground);self.members=frozenset(frozenset(s) for s in members)
if frozenset() not in self.members or any(not s<=self.ground for s in self.members):
raise ValueError('Members must lie in the ground set and include the empty set.')
if any(a|b not in self.members for a in self.members for b in self.members):
raise ValueError('The supplied family is not union-closed.')
def interior(self,A):
A=frozenset(A)
if not A<=self.ground:raise ValueError('Unknown coordinate.')
return frozenset().union(*(s for s in self.members if s<=A))
def blockers(self,v):
return tuple(B for B in subsets(self.ground-{v}) if v not in self.interior(self.ground-B))
@property
def accessible(self):
return all(not S or any(S-{v} in self.members for v in S) for S in self.members)
class FamilyRealization:
"""All-blocker construction; marker count may be exponential."""
def __init__(self,family):
self.family=family
self.blockers={v:family.blockers(v) for v in sorted(family.ground)}
self.markers={(v,B):f'marker:{v}:{i}' for v in sorted(family.ground)
for i,B in enumerate(self.blockers[v])}
def requirements(self,v):return frozenset(self.markers[v,B] for B in self.blockers[v])
def outputs(self,v):return frozenset(name for (u,B),name in self.markers.items() if v in B)
def paired(self):
reactions={}
for v in sorted(self.family.ground):
reactions['p:'+v]=Reaction({'food'},{'alpha:'+v}|self.outputs(v),{'beta:'+v})
reactions['g:'+v]=Reaction({'food'}|self.requirements(v),{'beta:'+v},{'alpha:'+v})
return CatalyticSystem({'food'},reactions)
@staticmethod
def encode(A):return frozenset(prefix+v for v in A for prefix in ('p:','g:'))
def complete_pairs(self,T):
return frozenset(v for v in self.family.ground if 'p:'+v in T and 'g:'+v in T)
def exact_ground(self,predecessors):
if not self.family.accessible:raise ValueError('Exact-ground food feasibility must be accessible.')
if set(predecessors)!=self.family.ground or any(not set(P)<=self.family.ground for P in predecessors.values()):
raise ValueError('Supply a predecessor subset for every coordinate.')
return CatalyticSystem({'food'},{v:Reaction({'food'}|self.requirements(v),
{'private:'+v}|self.outputs(v),{'private:'+p for p in predecessors[v]}) for v in sorted(self.family.ground)})
def audit_pairs(self):
q=self.paired();expected={self.encode(A) for A in self.family.members}
actual=set(q.family());assert actual==expected
queries=[];max_rounds=0
for T in subsets(q.names):
wanted=self.encode(self.family.interior(self.complete_pairs(T)))
found=q.max_raf(T);assert wanted==found
rounds=len(q.closure_stages(T))-1;assert rounds<=2;max_rounds=max(max_rounds,rounds)
queries.append({'available':sorted(T),'complete_pairs':sorted(self.complete_pairs(T)),
'max_raf':sorted(found),'closure_growth_rounds':rounds})
singletons=frozenset(v for v in self.family.ground if frozenset({v}) in self.family.members)
assert q.max_raf(q.food_ready)==self.encode(singletons)
for r in q.names:
assert len(q.reactions[r].catalysts)==1 and not q.reactions[r].products&q.reactions[r].catalysts
suppliers=q.predecessors(r,q.names)
assert suppliers==frozenset({'g:'+r[2:] if r.startswith('p:') else 'p:'+r[2:]})
return {'source_family':rows_json(self.family.members),'encoded_family':rows_json(actual),
'reaction_count':len(q.names),'marker_count':len(self.markers),
'max_closure_growth_rounds':max_rounds,'elementary_max_raf':sorted(q.max_raf(q.food_ready)),
'frequencies':{v:sum(v in A for A in self.family.members) for v in sorted(self.family.ground)},
'queries':queries,'system':q.as_json()}
def worked_system(name):
# Literal remaining manuscript examples; bracketed molecules are catalysts.
examples={
'feedback':{'a':(('f',),('ca',),('cb',)), 'b':(('f',),('cb','z'),('ca',)), 'g':(('z',),('cb','d'),('f',))},
'false_average':{'a':(('f',),('z',),('f',)), 'g':(('z',),('cg',),('cl',)), 'l':(('f',),('cl',),('cg',))},
'failed_projection':{'a':(('f',),('u',),('f',)), 'b':(('f',),('z',),('c',)), 'g':(('z',),('c',),('f',))},
}
return CatalyticSystem({'f'},{r:Reaction(*data) for r,data in examples[name].items()})
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)
system=CatalyticSystem(FOOD,REACTIONS)
report=CoreAnalysis(system,CORE).report(CONTEXT_WEIGHTS)
diagnostics=ModuleDiagnostics(system,MODULE)
graph=FoodReadyGraph(system);packing=graph.greedy_packing()
rare={r for r,count in report['counts'].items() if 2*count<report['N_including_empty']}
assert graph.shortest_cycle(rare) is None
source=FamilyRealization(UnionClosedFamily(UNION_GROUND,UNION_MEMBERS)).audit_pairs()
worked={}
for name in ('feedback','false_average','failed_projection'):
q=worked_system(name);U={'a','b'} if name=='feedback' else {'a'}
E={'a','b'} if name!='false_average' else {'a','l'}
worked[name]={'analysis':CoreAnalysis(q,U).report(),
'projection':ModuleDiagnostics(q,E).projection_queries(),
'defect':ModuleDiagnostics(q,E).defect()}
result={'input_system':system.as_json(),'core_report':report,
'module_projection':diagnostics.projection_queries(),'module_defect':diagnostics.defect(),
'greedy_disjoint_cycles':[list(c) for c in packing],
'abundance_candidate_certificate':graph.abundance_candidates(),
'abundant_reaction_lower_bound':len(packing),'rare_food_ready_subgraph_acyclic':True,
'paired_realization':source,'worked_examples':worked,
'scope':'finite exact audits; no kinetics, no Lean rerun, no unrestricted half-frequency proof'}
(out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
def csv_file(name,headers,rows):
with (out/name).open('w',newline='') as stream:
writer=csv.writer(stream);writer.writerow(headers);writer.writerows(rows)
csv_file('fibres.csv',['exterior','member_count',*sorted(CORE),'slack'],
[('|'.join(r['exterior']),len(r['family']),*(r['counts'][v] for v in sorted(CORE)),r['slack']) for r in report['fibres']])
csv_file('frequencies.csv',['reaction','count','N_including_empty','frequency'],
[(r,n,report['N_including_empty'],str(Q(n,report['N_including_empty']))) for r,n in report['counts'].items()])
csv_file('availability_queries.csv',['available','complete_pairs','max_raf','closure_growth_rounds'],
[('|'.join(r['available']),'|'.join(r['complete_pairs']),'|'.join(r['max_raf']),r['closure_growth_rounds']) for r in source['queries']])
csv_file('cell_injections.csv',['exterior','source_deleted','source_column','move','target_deleted','target_column'],
[('|'.join(f['exterior']),'|'.join(c['source_deleted']),c['source_column'],c['move'],'|'.join(c['target_deleted']),c['target_column'])
for f in report['fibres'] for c in f['injection']])
lines=[f'RAF family including empty: N={report["N_including_empty"]}; frequencies={report["counts"]}.',
f'Fibre slacks: {[f["slack"] for f in report["fibres"]]}.',
f'Weighted core frequencies: {report["weighted_core_frequencies"]}.',
f'Projection-preserving module: {result["module_projection"]["preserves_projection"]}.',
f'Greedy disjoint-cycle certificate: at least {len(packing)} distinct abundant reactions.',
f'Producer-gate benchmark: {source["reaction_count"]} reactions; {source["marker_count"]} markers; {len(source["queries"])} availability queries checked.',
f'Elementary maxRAF in benchmark: {source["elementary_max_raf"]}.',
'Abundance counts structural sets including the empty set; it is not a kinetic probability.']
(out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig,axes=plt.subplots(1,2,figsize=(10,4),layout='constrained')
contexts=report['fibres'];core=sorted(CORE);width=.8/max(1,len(core))
for i,r in enumerate(core):
axes[0].bar([j-.4+width*(i+.5) for j in range(len(contexts))],
[f['counts'][r]/len(f['family']) if f['family'] else float('nan') for f in contexts],width,label=r)
axes[0].set_xticks(range(len(contexts)),['empty' if not f['exterior'] else '+'.join(f['exterior']) for f in contexts])
axes[0].set(xlabel='Fixed exterior context',ylabel='Fraction of compatible core subsets',ylim=(0,1),title='Core-reaction frequencies by exterior\nselection')
axes[0].legend(title='Core reaction')
names=sorted(report['counts'])
axes[1].bar(names,[report['counts'][r]/report['N_including_empty'] for r in names],color=['#287a99' if r in CORE else '#999999' for r in names])
axes[1].set(xlabel='Reaction',ylabel='Fraction of all RAF-family members',ylim=(0,1),title='Reaction frequencies across the full RAF\nfamily')
for ax in axes:ax.axhline(.5,color='black',ls='--',lw=1);ax.grid(axis='y',alpha=.2)
fig.savefig(out/'abundance.png',dpi=180);fig.savefig(out/'abundance.svg');plt.close(fig)
digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
(out/'run_metadata.json').write_text(json.dumps({'manuscript_sha256':MANUSCRIPT_SHA256,
'source_sha256':digest(Path(__file__)),'python':platform.python_version(),
'output_sha256':{p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}},indent=2)+'\n')
if __name__=='__main__':main()
Run output
RAF family including empty: N=12; frequencies={'a': 6, 'b': 6, 'g': 7, 'h': 7}.
Fibre slacks: [0, 0, 0, 0].
Weighted core frequencies: {'a': '5/9', 'b': '4/9'}.
Projection-preserving module: False.
Greedy disjoint-cycle certificate: at least 3 distinct abundant reactions.
Producer-gate benchmark: 8 reactions; 14 markers; 256 availability queries checked.
Elementary maxRAF in benchmark: [].
Abundance counts structural sets including the empty set; it is not a kinetic probability.