Example code
When designing an autocatalytic reaction network, repeatedly deleting dispensable reactions can leave a much larger network than necessary. This example builds the paper's Set Cover construction, checks its exact correspondence, and makes that gap visible through two deletion orders.
A RAF is a nonempty reaction set whose reactants can be generated from food and whose reactions have catalysts in that food closure. An irreducible RAF contains no smaller RAF. A minimum RAF has the fewest reactions among all RAFs. The difference matters when using deletion as a design or search method.
The input is a Set Cover problem: choose from a numbered collection of sets so that every one of elements is covered. The construction makes required gate reactions and a block of reactions for each set. Each block needs the final gate product; its own final product catalyses its block and the gates corresponding to its elements. A RAF must contain all gates and whole blocks that collectively cover every element. Thus a cover corresponds to a RAF of size .
The default reproduces the paper's three-element example with three sets and two reactions per block. Testing all 512 reaction subsets finds exactly three RAFs, of sizes 7, 7, and 9, matching the three covers. The first figure separates the common gate cost from the selected blocks.
For the second figure, one set covers the whole universe and the remaining sets each cover one element. With block length equal to universe size , two deletion orders return irreducible RAFs of sizes and . At , that is 16 versus 72 reactions. Both are inclusion-minimal; only the smaller one is minimum.
The package also sweeps block length and records exact rational approximation ratios. When , a -approximate RAF yields a cover with ratio at most . Smaller blocks preserve the size identity but do not supply this guarantee: the fixed gate scaffold can hide a poor cover ratio.
Download the package, install requirements.txt, run python example.py --output outputs, and check it with python -m unittest -v. Edit the incidence sets and integer inputs at the top of the source. Import SetCover, Network, and CoverReduction to build other instances, inspect closure stages, try deletion orders, or connect a cover heuristic to the checked RAF decoder. Exact enumeration has explicit size caps.
The tests use literal reaction-subset checks, a family of 98 small coverable systems, and exact ratio arithmetic. They illustrate the construction and its consequences; the general complexity-hardness claim remains the manuscript theorem. This code models structural food closure and catalysis, not concentration dynamics, and does not compile or extract the Lean proof.
Python source
"""Executable Set Cover -> Min-RAF reduction from the paper's Construction 3.1.
Finite checks teach the reduction; they do not prove complexity hardness.
Run python example.py --output outputs.
"""
from __future__ import annotations
import argparse
import csv
from dataclasses import dataclass
from fractions import Fraction
import hashlib
from itertools import combinations
import json
from pathlib import Path
import platform
import time
# USER INPUTS ---------------------------------------------------------------
UNIVERSE_SIZE = 3 # m >= 1, elements numbered 0 through m-1
SETS = ((0,1), (1,2), (2,)) # paper's Figure 1; indexed sets may repeat
BLOCK_LENGTH = 2 # M >= 1; Figure 1 uses 2; AP bound needs M >= m
AMPLIFICATION_SWEEP = (1,2,3,6,12)
GAP_UNIVERSE_SIZES = (2,3,4,5,6,7,8) # whole-universe set versus singleton sets
EXHAUSTIVE_REACTION_CAP = 18 # subset oracle budget; not a scientific input
EXHAUSTIVE_SET_CAP = 18
PAPER_SHA256 = 'd0e814c5b965bffccfcdb5e6fe542c9d6d64dfe199e9d0c9814a0c5194bc6b79'
# Inputs are combinatorial counts and incidence sets, with no physical units.
# --------------------------------------------------------------------------
def subsets(values):
values=tuple(values)
for size in range(len(values)+1):
for selected in combinations(values,size): yield frozenset(selected)
@dataclass(frozen=True)
class SetCover:
universe_size: int
sets: tuple[frozenset[int],...]
def __post_init__(self):
if not isinstance(self.universe_size,int) or isinstance(self.universe_size,bool) or self.universe_size<1:
raise ValueError('Universe size must be a positive integer.')
sets=tuple(frozenset(s) for s in self.sets);object.__setattr__(self,'sets',sets)
universe=frozenset(range(self.universe_size))
if any(any(not isinstance(i,int) or isinstance(i,bool) for i in s) or not s<=universe for s in sets):
raise ValueError('Every set must contain only valid integer universe elements.')
if frozenset().union(*sets)!=universe:
raise ValueError('The union of the indexed sets must cover the universe.')
def covers(self,selected):
selected=frozenset(selected)
if not selected<=set(range(len(self.sets))): raise ValueError('Unknown set index.')
return frozenset().union(*(self.sets[j] for j in selected))==frozenset(range(self.universe_size))
def all_covers(self):
if len(self.sets)>EXHAUSTIVE_SET_CAP: raise ValueError('Exact cover search exceeds declared set budget.')
return tuple(d for d in subsets(range(len(self.sets))) if self.covers(d))
def minimal_covers(self):
return tuple(d for d in self.all_covers() if all(not self.covers(d-{j}) for j in d))
@dataclass(frozen=True)
class Reaction:
identifier: str
reactants: frozenset[str]
product: str
catalysts: frozenset[str]
class Network:
"""Literal food-closure semantics, independent of cover membership."""
def __init__(self,reactions,food=frozenset({'f'})):
reactions=tuple(reactions)
self.reactions={r.identifier:r for r in reactions};self.food=frozenset(food)
if len(self.reactions)!=len(reactions): raise ValueError('Reaction identifiers must be unique.')
self.identifiers=frozenset(self.reactions)
def selected(self,ids):
ids=frozenset(ids)
if not ids<=self.identifiers: raise ValueError('Unknown reaction identifier.')
return ids
def closure_stages(self,ids):
ids=self.selected(ids);stages=[self.food]
while True:
current=stages[-1]
following=current|frozenset(self.reactions[r].product for r in ids if self.reactions[r].reactants<=current)
if following==current: return stages
stages.append(following)
def is_raf(self,ids):
ids=self.selected(ids);closure=self.closure_stages(ids)[-1]
return bool(ids) and all(self.reactions[r].reactants<=closure and self.reactions[r].catalysts&closure for r in ids)
def maximal_raf(self,ids=None):
ids=self.identifiers if ids is None else self.selected(ids)
while ids:
closure=self.closure_stages(ids)[-1]
kept=frozenset(r for r in ids if self.reactions[r].reactants<=closure and self.reactions[r].catalysts&closure)
if kept==ids: return ids
ids=kept
return frozenset()
def shrink(self,order):
"""One deletion/pruning sweep returns an irreducible, not necessarily smallest, RAF."""
order=tuple(order)
if len(order)!=len(self.identifiers) or set(order)!=self.identifiers:
raise ValueError('Deletion order must list every reaction exactly once.')
current=self.maximal_raf()
for r in order:
if r not in current: continue
candidate=self.maximal_raf(current-{r})
if candidate: current=candidate
return current
def all_rafs(self):
if len(self.identifiers)>EXHAUSTIVE_REACTION_CAP:
raise ValueError('Literal exhaustive reaction search exceeds declared budget.')
return tuple(s for s in subsets(sorted(self.identifiers)) if self.is_raf(s))
class CoverReduction:
def __init__(self,instance:SetCover,block_length:int):
if not isinstance(block_length,int) or isinstance(block_length,bool) or block_length<1:
raise ValueError('Block length must be a positive integer.')
self.instance,self.block_length=instance,block_length
m=instance.universe_size;M=block_length
self.gates=frozenset(f'g{i}' for i in range(m))
self.blocks={j:frozenset(f'b{j}:{k}' for k in range(M)) for j in range(len(instance.sets))}
reactions=[]
for i in range(m):
reactants={'f'}|({f'y{i-1}'} if i else set())
catalysts={f'z{j}:{M-1}' for j,s in enumerate(instance.sets) if i in s}
reactions.append(Reaction(f'g{i}',frozenset(reactants),f'y{i}',frozenset(catalysts)))
for j in self.blocks:
for k in range(M):
previous=f'y{m-1}' if k==0 else f'z{j}:{k-1}'
reactions.append(Reaction(f'b{j}:{k}',frozenset({'f',previous}),f'z{j}:{k}',frozenset({f'z{j}:{M-1}'})))
self.network=Network(reactions)
def canonical(self,selected):
selected=frozenset(selected)
if not selected<=self.blocks.keys(): raise ValueError('Unknown set index.')
return self.gates.union(*(self.blocks[j] for j in selected))
def decode(self,ids):
ids=self.network.selected(ids)
if not self.network.is_raf(ids): raise ValueError('Decoder requires a feasible RAF.')
selected=frozenset(j for j in self.blocks if f'b{j}:{self.block_length-1}' in ids)
if not self.instance.covers(selected) or self.canonical(selected)!=ids:
raise ArithmeticError('Literal RAF violates the construction correspondence.')
return selected
def audit(self):
covers=self.instance.all_covers()
predicted={self.canonical(d) for d in covers}
actual=set(self.network.all_rafs())
if actual!=predicted: raise ArithmeticError('Exhaustive correspondence failed.')
tau=min(map(len,covers)); optimum=min(map(len,actual))
if optimum!=self.instance.universe_size+self.block_length*tau:
raise ArithmeticError('Optimum identity failed.')
return dict(reaction_subsets_checked=2**len(self.network.identifiers),raf_count=len(actual),
minimum_cover_size=tau,minimum_raf_size=optimum)
def approximation(self,selected):
if not self.instance.covers(selected): raise ValueError('Approximation candidate must be a cover.')
tau=min(map(len,self.instance.all_covers()));m=self.instance.universe_size;M=self.block_length
c=Fraction(len(self.canonical(selected)),m+M*tau)
cover_ratio=Fraction(len(selected),tau)
return dict(block_length=M,raf_ratio=str(c),cover_ratio=str(cover_ratio),
transferred_bound=str(2*c-1) if M>=m else None,
bound_applicable=M>=m,
bound_holds=bool(cover_ratio<=2*c-1) if M>=m else None)
def record(self):
return dict(universe_size=self.instance.universe_size,sets=[sorted(s) for s in self.instance.sets],
block_length=self.block_length,food=['f'],reactions=[dict(id=r.identifier,
reactants=sorted(r.reactants),product=r.product,catalysts=sorted(r.catalysts))
for r in self.network.reactions.values()])
def gap_instance(q):
if not isinstance(q,int) or q<1: raise ValueError('Gap universe size must be positive integer.')
return SetCover(q,(frozenset(range(q)),*(frozenset({i}) for i in range(q))))
def compare_deletion_orders(q):
reduction=CoverReduction(gap_instance(q),q);network=reduction.network
universal=reduction.blocks[0]
singleton=frozenset().union(*(v for j,v in reduction.blocks.items() if j!=0))
large=network.shrink(sorted(universal)+sorted(singleton)+sorted(reduction.gates))
small=network.shrink(sorted(singleton)+sorted(universal)+sorted(reduction.gates))
if reduction.decode(small)!={0} or reduction.decode(large)!=set(range(1,q+1)):
raise ArithmeticError('Deletion witness differs from paper family.')
return dict(universe_size=q,block_length=q,reactions=len(network.identifiers),
small_irraf=len(small),large_irraf=len(large),ratio=str(Fraction(len(large),len(small))))
def write_csv(path,rows):
with path.open('w',newline='',encoding='utf-8') as f:
writer=csv.DictWriter(f,fieldnames=list(rows[0]));writer.writeheader();writer.writerows(rows)
def plot(covers,gaps,amplification,output):
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
with plt.rc_context({'font.size':11,'axes.spines.top':False,'axes.spines.right':False,
'figure.facecolor':'white','savefig.facecolor':'white','svg.fonttype':'none'}):
fig,ax=plt.subplots(figsize=(7.2,4.5),layout='constrained')
labels=[r['set_indices'] for r in covers]
ax.bar(labels,[UNIVERSE_SIZE]*len(covers),label='Required gate scaffold',color='#0072B2')
ax.bar(labels,[r['raf_size']-UNIVERSE_SIZE for r in covers],bottom=UNIVERSE_SIZE,label='Selected set blocks',color='#D55E00')
ax.set(xlabel='Selected cover indices',ylabel='Reactions in the corresponding RAF')
ax.legend(fontsize=9);ax.grid(axis='y',color='#e3e6e8');ax.set_axisbelow(True)
fig.savefig(output/'correspondence.png',dpi=220);fig.savefig(output/'correspondence.svg');plt.close(fig)
fig,ax=plt.subplots(figsize=(7.2,4.5),layout='constrained')
ax.plot([r['universe_size'] for r in gaps],[r['small_irraf'] for r in gaps],'o-',color='#0072B2',label='Keep the whole-universe block')
ax.plot([r['universe_size'] for r in gaps],[r['large_irraf'] for r in gaps],'s--',color='#D55E00',label='Keep the singleton blocks')
ax.set(xlabel='Universe size q (block length also q)',ylabel='Reactions in returned irreducible RAF',ylim=(0,None))
ax.legend(fontsize=9);ax.grid(axis='y',color='#e3e6e8')
fig.savefig(output/'irreducible_gap.png',dpi=220);fig.savefig(output/'irreducible_gap.svg');plt.close(fig)
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);start=time.perf_counter()
instance=SetCover(UNIVERSE_SIZE,SETS);reduction=CoverReduction(instance,BLOCK_LENGTH)
inputs=dict(universe_size=UNIVERSE_SIZE,sets=SETS,block_length=BLOCK_LENGTH,
amplification_sweep=AMPLIFICATION_SWEEP,gap_sizes=GAP_UNIVERSE_SIZES,
exhaustive_reaction_cap=EXHAUSTIVE_REACTION_CAP,exhaustive_set_cap=EXHAUSTIVE_SET_CAP)
intro='Resolved inputs: '+json.dumps(inputs);print(intro,flush=True)
audit=reduction.audit();covers=[]
for d in instance.all_covers():
selected=reduction.canonical(d)
covers.append(dict(set_indices=' '.join(map(str,sorted(d))),cover_size=len(d),raf_size=len(selected),
inclusion_minimal=all(not instance.covers(d-{j}) for j in d),
reaction_ids=' '.join(sorted(selected)),closure_depth=len(reduction.network.closure_stages(selected))-1))
# Fixed gap-family cover makes amplification's effect visible and comparable.
amplified=gap_instance(4);all_singletons=frozenset(range(1,5))
amplification=[CoverReduction(amplified,M).approximation(all_singletons) for M in AMPLIFICATION_SWEEP]
gaps=[compare_deletion_orders(q) for q in GAP_UNIVERSE_SIZES]
write_csv(args.output/'covers.csv',covers);write_csv(args.output/'amplification.csv',amplification)
write_csv(args.output/'irreducible_gap.csv',gaps)
(args.output/'network.json').write_text(json.dumps(reduction.record(),indent=2)+'\n',encoding='utf-8')
stages=reduction.network.closure_stages(reduction.canonical(instance.minimal_covers()[0]))
write_csv(args.output/'closure_stages.csv',[dict(stage=i,molecules=' '.join(sorted(s))) for i,s in enumerate(stages)])
plot(covers,gaps,amplification,args.output)
summary=dict(inputs=inputs,exhaustive_audit=audit,covers=covers,amplification=amplification,irreducible_gap=gaps,
evidence='Exact finite incidence/closure/subset calculations and rational ratio arithmetic. Complexity hardness remains the manuscript theorem under P != NP; no Lean compilation.')
transcript=json.dumps(summary,indent=2)
(args.output/'summary.json').write_text(transcript+'\n',encoding='utf-8')
(args.output/'console.txt').write_text(intro+'\n'+transcript+'\n',encoding='utf-8')
import matplotlib
metadata=dict(paper_sha256=PAPER_SHA256,source_sha256=hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
python=platform.python_version(),matplotlib=matplotlib.__version__,platform=platform.platform(),processor=platform.processor(),
elapsed_seconds=time.perf_counter()-start,command='python example.py --output outputs',seed_policy='No random sampling.',
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',encoding='utf-8');print(transcript)
if __name__=='__main__': main()
Run output
Resolved inputs: {"universe_size": 3, "sets": [[0, 1], [1, 2], [2]], "block_length": 2, "amplification_sweep": [1, 2, 3, 6, 12], "gap_sizes": [2, 3, 4, 5, 6, 7, 8], "exhaustive_reaction_cap": 18, "exhaustive_set_cap": 18}
{
"inputs": {
"universe_size": 3,
"sets": [
[
0,
1
],
[
1,
2
],
[
2
]
],
"block_length": 2,
"amplification_sweep": [
1,
2,
3,
6,
12
],
"gap_sizes": [
2,
3,
4,
5,
6,
7,
8
],
"exhaustive_reaction_cap": 18,
"exhaustive_set_cap": 18
},
"exhaustive_audit": {
"reaction_subsets_checked": 512,
"raf_count": 3,
"minimum_cover_size": 2,
"minimum_raf_size": 7
},
"covers": [
{
"set_indices": "0 1",
"cover_size": 2,
"raf_size": 7,
"inclusion_minimal": true,
"reaction_ids": "b0:0 b0:1 b1:0 b1:1 g0 g1 g2",
"closure_depth": 5
},
{
"set_indices": "0 2",
"cover_size": 2,
"raf_size": 7,
"inclusion_minimal": true,
"reaction_ids": "b0:0 b0:1 b2:0 b2:1 g0 g1 g2",
"closure_depth": 5
},
{
"set_indices": "0 1 2",
"cover_size": 3,
"raf_size": 9,
"inclusion_minimal": false,
"reaction_ids": "b0:0 b0:1 b1:0 b1:1 b2:0 b2:1 g0 g1 g2",
"closure_depth": 5
}
],
"amplification": [
{
"block_length": 1,
"raf_ratio": "8/5",
"cover_ratio": "4",
"transferred_bound": null,
"bound_applicable": false,
"bound_holds": null
},
{
"block_length": 2,
"raf_ratio": "2",
"cover_ratio": "4",
"transferred_bound": null,
"bound_applicable": false,
"bound_holds": null
},
{
"block_length": 3,
"raf_ratio": "16/7",
"cover_ratio": "4",
"transferred_bound": null,
"bound_applicable": false,
"bound_holds": null
},
{
"block_length": 6,
"raf_ratio": "14/5",
"cover_ratio": "4",
"transferred_bound": "23/5",
"bound_applicable": true,
"bound_holds": true
},
{
"block_length": 12,
"raf_ratio": "13/4",
"cover_ratio": "4",
"transferred_bound": "11/2",
"bound_applicable": true,
"bound_holds": true
}
],
"irreducible_gap": [
{
"universe_size": 2,
"block_length": 2,
"reactions": 8,
"small_irraf": 4,
"large_irraf": 6,
"ratio": "3/2"
},
{
"universe_size": 3,
"block_length": 3,
"reactions": 15,
"small_irraf": 6,
"large_irraf": 12,
"ratio": "2"
},
{
"universe_size": 4,
"block_length": 4,
"reactions": 24,
"small_irraf": 8,
"large_irraf": 20,
"ratio": "5/2"
},
{
"universe_size": 5,
"block_length": 5,
"reactions": 35,
"small_irraf": 10,
"large_irraf": 30,
"ratio": "3"
},
{
"universe_size": 6,
"block_length": 6,
"reactions": 48,
"small_irraf": 12,
"large_irraf": 42,
"ratio": "7/2"
},
{
"universe_size": 7,
"block_length": 7,
"reactions": 63,
"small_irraf": 14,
"large_irraf": 56,
"ratio": "4"
},
{
"universe_size": 8,
"block_length": 8,
"reactions": 80,
"small_irraf": 16,
"large_irraf": 72,
"ratio": "9/2"
}
],
"evidence": "Exact finite incidence/closure/subset calculations and rational ratio arithmetic. Complexity hardness remains the manuscript theorem under P != NP; no Lean compilation."
}