"""Complete structural irrRAF catalogues via deterministic supplier resolutions. Catalysts do NOT gate food closure. Catalysis is disjunctive. This is a finite structural calculation, not a kinetic or thermodynamic reactor simulation. """ from __future__ import annotations import argparse import csv from dataclasses import dataclass, asdict from fractions import Fraction import hashlib from itertools import combinations, product import json import math from pathlib import Path import platform import time # USER INPUTS --------------------------------------------------------------- FOOD = ('f',) # species freely available to every reaction # id, reactants, products, alternative catalysts; paper's k=2 gate example. REACTIONS = ( (0, ('f',), ('x1',), ('z2',)), (1, ('f',), ('x1',), ('z2',)), (2, ('f',), ('x2',), ('z2',)), (3, ('f',), ('x2',), ('z2',)), (4, ('x1',), ('z1',), ('z2',)), (5, ('z1', 'x2'), ('z2',), ('z2',)), ) DELETION_COSTS = {0: 2, 1: 3, 2: 4, 3: 1, 4: 7, 5: 6} # positive abstract costs AVAILABLE_PROBABILITY = '1/2' # independent static reaction availability REQUIRED_REACTIONS = (0,) FORBIDDEN_REACTIONS = (2,) DELETED_REACTIONS = (0,) SCALING_SIZES = (1, 2, 3, 4, 5, 6) # paper families; controlled size sweep LARGE_INDEPENDENT_MODULES = 64 DETERMINISTIC_CHAIN_LENGTH = 64 RESOLUTION_BUDGET = 4096 # fail explicitly before exceeding this budget EXHAUSTIVE_REACTION_CAP = 18 # small-instance verification/queries only PAPER_SHA256 = 'c92ad401d52719b68ab0118296ae5dc686c994cbaf9bea19a5bb8c4c4820523a' # All inputs are finite identifiers/sets, counts, costs or probabilities; # there are no concentration units or rate constants in the RAF predicate. # -------------------------------------------------------------------------- @dataclass(frozen=True) class Reaction: identifier: int reactants: frozenset[str] products: frozenset[str] catalysts: frozenset[str] def __post_init__(self): if not isinstance(self.identifier,int) or isinstance(self.identifier,bool): raise ValueError('Reaction identifiers must be integers.') for field in ('reactants','products','catalysts'): values=frozenset(getattr(self,field)) if any(not isinstance(x,str) or not x for x in values): raise ValueError('Species names must be nonempty strings.') object.__setattr__(self,field,values) class ReactionSystem: def __init__(self,food,reactions): self.food=frozenset(food) if any(not isinstance(x,str) or not x for x in self.food): raise ValueError('Food species names must be nonempty strings.') reactions=tuple(reactions) self.reactions={r.identifier:r for r in reactions} if len(self.reactions)!=len(reactions): raise ValueError('Reaction identifiers must be unique.') self.identifiers=frozenset(self.reactions) def selected(self,ids=None): chosen=self.identifiers if ids is None else frozenset(ids) if not chosen<=self.identifiers: raise ValueError('Unknown reaction identifier.') return chosen def closure_stages(self,ids=None): selected=self.selected(ids); stages=[self.food] while True: current=stages[-1] new=current.union(*(self.reactions[r].products for r in selected if self.reactions[r].reactants<=current)) if new==current: return stages stages.append(new) def closure(self,ids=None): return self.closure_stages(ids)[-1] def is_raf(self,ids): ids=self.selected(ids);closure=self.closure(ids) 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): current=self.selected(ids) while current: closure=self.closure(current) remaining=frozenset(r for r in current if self.reactions[r].reactants<=closure and self.reactions[r].catalysts&closure) if remaining==current: return current current=remaining return frozenset() def restrict(self,ids): return ReactionSystem(self.food,[self.reactions[r] for r in sorted(self.selected(ids))]) def record(self): return dict(food=sorted(self.food),reactions=[dict(id=r.identifier,reactants=sorted(r.reactants), products=sorted(r.products),catalysts=sorted(r.catalysts)) for r in self.reactions.values()]) @classmethod def from_record(cls,data): return cls(data['food'],[Reaction(r['id'],r['reactants'],r['products'],r['catalysts']) for r in data['reactions']]) @dataclass class SupplierOptions: maximal: frozenset[int] producers: dict[str,tuple[int,...]] catalysts: dict[int,tuple[str,...]] @classmethod def from_system(cls,system): maximal=system.maximal_raf() available=system.food.union(*(system.reactions[r].products for r in maximal)) producers={x:tuple(sorted(r for r in maximal if x in system.reactions[r].products)) for x in sorted(available-system.food)} catalysts={} for r in sorted(maximal): food=system.reactions[r].catalysts&system.food catalysts[r]=(min(food),) if food else tuple(sorted(system.reactions[r].catalysts&available)) return cls(maximal,producers,catalysts) @property def beta(self): return sum(len(v)-1 for v in (*self.producers.values(),*self.catalysts.values())) @property def count(self): return math.prod(len(v) for v in (*self.producers.values(),*self.catalysts.values())) def resolutions(self): keys=list(self.producers)+list(self.catalysts) options=list(self.producers.values())+list(self.catalysts.values()) boundary=len(self.producers) for values in product(*options): yield dict(zip(keys[:boundary],values[:boundary])),dict(zip(keys[boundary:],values[boundary:])) def resolve(system,options,producer,catalyst): return ReactionSystem(system.food,[Reaction(r,system.reactions[r].reactants, frozenset(x for x in system.reactions[r].products if x in system.food or producer[x]==r), frozenset((catalyst[r],))) for r in sorted(options.maximal)]) def sink_components(graph): """Iterative Kosaraju SCCs, avoiding recursion-depth limits on long chains.""" visited=set();order=[] for root in sorted(graph): if root in visited: continue visited.add(root); stack=[(root,iter(sorted(graph[root])))] while stack: node,children=stack[-1] child=next(children,None) if child is None: order.append(node);stack.pop() elif child not in visited: visited.add(child);stack.append((child,iter(sorted(graph[child])))) reverse={v:set() for v in graph} for v,edges in graph.items(): for w in edges: reverse[w].add(v) visited=set();components=[] for root in reversed(order): if root in visited: continue component=set();stack=[root];visited.add(root) while stack: node=stack.pop();component.add(node) for child in reverse[node]-visited: visited.add(child);stack.append(child) components.append(frozenset(component)) return [c for c in components if all(graph[v]<=c for v in c)] @dataclass class Catalogue: members: tuple[frozenset[int],...] beta: int global_resolution_count: int resolutions_inspected: int candidates: int rejected: int validation_calls: int component_count: int largest_local_beta: int def filter(self,required=(),forbidden=()): required,forbidden=frozenset(required),frozenset(forbidden) if required&forbidden: raise ValueError('Required and forbidden identifiers overlap.') return tuple(c for c in self.members if required<=c and not forbidden&c) def frequencies(self): return {r:sum(r in c for c in self.members) for r in sorted(set().union(*self.members))} class SupplierEnumerator: def __init__(self,resolution_budget=RESOLUTION_BUDGET): if not isinstance(resolution_budget,int) or resolution_budget<1: raise ValueError('Resolution budget must be a positive integer.') self.budget=resolution_budget def global_catalogue(self,system): options=SupplierOptions.from_system(system) if options.count>self.budget: raise ValueError(f'{options.count} resolutions exceed budget {self.budget}; use components or raise budget explicitly.') members=set();candidates=rejected=calls=inspected=0 if options.maximal: for producer,catalyst in options.resolutions(): inspected+=1 resolved=resolve(system,options,producer,catalyst) kept=resolved.maximal_raf() graph={r:{producer[x] for x in (resolved.reactions[r].reactants|resolved.reactions[r].catalysts)-system.food} for r in kept} if any(not edges<=kept for edges in graph.values()): raise ArithmeticError('Dependency escaped pruned system.') resolution_calls=0 for component in sink_components(graph): candidates+=1;valid=True for r in sorted(component): calls+=1;resolution_calls+=1 if system.maximal_raf(component-{r}): valid=False;break if valid: members.add(component) else: rejected+=1 if resolution_calls>len(options.maximal): raise ArithmeticError('Validation budget invariant failed.') ordered=tuple(sorted(members,key=lambda c:tuple(sorted(c)))) result=Catalogue(ordered,options.beta,options.count,inspected,candidates,rejected,calls, 1 if options.maximal else 0,options.beta) if sum(map(len,ordered))>len(options.maximal)*options.count: raise ArithmeticError('Output membership bound failed.') return result def components(self,system): options=SupplierOptions.from_system(system) # Every nonfood reactant, product and effective catalyst incidence counts. touching={} for r in options.maximal: reaction=system.reactions[r] for x in (reaction.reactants|reaction.products|frozenset(options.catalysts[r]))-system.food: touching.setdefault(x,set()).add(r) adjacency={r:set() for r in options.maximal} for group in touching.values(): first=min(group) for r in group: adjacency[first].add(r);adjacency[r].add(first) blocks=[];seen=set() for r in sorted(adjacency): if r in seen: continue block=set();stack=[r];seen.add(r) while stack: v=stack.pop();block.add(v) for w in adjacency[v]-seen: seen.add(w);stack.append(w) blocks.append(frozenset(block)) return blocks def catalogue(self,system): options=SupplierOptions.from_system(system);blocks=self.components(system) local=[system.restrict(b) for b in blocks] required=sum(SupplierOptions.from_system(q).count for q in local) if required>self.budget: raise ValueError(f'Total local resolutions {required} exceed budget {self.budget}.') results=[self.global_catalogue(q) for q in local] return Catalogue(tuple(sorted((c for result in results for c in result.members),key=lambda c:tuple(sorted(c)))), options.beta,options.count,sum(r.resolutions_inspected for r in results), sum(r.candidates for r in results),sum(r.rejected for r in results), sum(r.validation_calls for r in results),len(blocks),max((r.beta for r in results),default=0)) def coupled_gates(k): if not isinstance(k,int) or k<1: raise ValueError('Gate count must be a positive integer.') reactions=[] for i in range(k): for alternative in range(2): reactions.append(Reaction(2*i+alternative,{'f'},{f'x{i+1}'},{f'z{k}'})) for i in range(k): reactants={f'x{i+1}'}|({f'z{i}'} if i else set()) reactions.append(Reaction(2*k+i,reactants,{f'z{i+1}'},{f'z{k}'})) return ReactionSystem({'f'},reactions) def independent_modules(count,choices=True): if not isinstance(count,int) or count<1: raise ValueError('Module count must be positive integer.') reactions=[];per=3 if choices else 2 for i in range(count): for a in range(per-1): reactions.append(Reaction(per*i+a,{'f'},{f'x{i}'},{f'z{i}'})) reactions.append(Reaction(per*i+per-1,{f'x{i}'},{f'z{i}'},{f'z{i}'})) return ReactionSystem({'f'},reactions) def deterministic_chain(length): if not isinstance(length,int) or length<1: raise ValueError('Chain length must be positive integer.') return ReactionSystem({'f'},[Reaction(i,{'f' if i==0 else f'x{i}'},{f'x{i+1}'},{f'x{length}'}) for i in range(length)]) def minimal_cuts(system,catalogue): if len(system.identifiers)>EXHAUSTIVE_REACTION_CAP: raise ValueError('Exact cut search exceeds subset budget.') found=[] for size in range(len(system.identifiers)+1): for values in combinations(sorted(system.identifiers),size): cut=frozenset(values) if not any(c<=cut for c in found) and all(cut&member for member in catalogue.members): found.append(cut) return tuple(found) def greedy_cut(catalogue,costs): costs={r:Fraction(str(c)) for r,c in costs.items()} if any(c<0 for c in costs.values()): raise ValueError('Deletion costs must be nonnegative.') remaining=set(catalogue.members);chosen=[] while remaining: options=[(cost/sum(r in c for c in remaining),r) for r,cost in costs.items() if any(r in c for c in remaining)] if not options: raise ValueError('Eligible reactions cannot hit all irrRAFs.') _,r=min(options);chosen.append(r);remaining={c for c in remaining if r not in c} return chosen,sum((costs[r] for r in chosen),Fraction(0)) def exact_availability(system,catalogue,p=Fraction(AVAILABLE_PROBABILITY)): p=Fraction(p) if not 0<=p<=1: raise ValueError('Availability must lie in [0,1].') n=len(system.identifiers) if n>EXHAUSTIVE_REACTION_CAP: raise ValueError('Exact availability search exceeds subset budget.') total=Fraction(0) for size in range(n+1): for ids in combinations(sorted(system.identifiers),size): if any(c<=frozenset(ids) for c in catalogue.members): total+=p**size*(1-p)**(n-size) return total def zero_excess_reliability(catalogue,p=Fraction(AVAILABLE_PROBABILITY)): p=Fraction(p) if catalogue.beta!=0: raise ValueError('Disjoint zero-excess formula does not apply to this catalogue.') if not 0<=p<=1: raise ValueError('Availability must lie in [0,1].') return 1-math.prod(1-p**len(c) for c in catalogue.members) def family_record(result): return dict(members=[sorted(c) for c in result.members],beta=result.beta, global_resolution_count=result.global_resolution_count,resolutions_inspected=result.resolutions_inspected, candidates=result.candidates,rejected=result.rejected,validation_calls=result.validation_calls, component_count=result.component_count,largest_local_beta=result.largest_local_beta, frequencies=result.frequencies(),total_memberships=sum(map(len,result.members))) def figures(system,catalogue,scaling,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'}): if catalogue.members and len(catalogue.members)<=32 and len(system.identifiers)<=30: fig,ax=plt.subplots(figsize=(7.2,4.5),layout='constrained') ids=sorted(system.identifiers) ax.imshow([[int(r in c) for r in ids] for c in catalogue.members],cmap='Blues',vmin=0,vmax=1,aspect='auto',interpolation='nearest') ax.set(xticks=range(len(ids)),xticklabels=ids,yticks=range(len(catalogue.members)), yticklabels=[f'I{i+1}' for i in range(len(catalogue.members))], xlabel='Original reaction identifier',ylabel='Irreducible RAF') for i,c in enumerate(catalogue.members): for j,r in enumerate(ids): ax.text(j,i,'1' if r in c else '0',ha='center',va='center',color='white' if r in c else '#202124') fig.savefig(output/'catalogue.png',dpi=220);fig.savefig(output/'catalogue.svg');plt.close(fig) fig,ax=plt.subplots(figsize=(7.2,4.5),layout='constrained') for family,mode,color,style in [('independent','global','#D55E00','--'),('independent','components','#0072B2','-'),('coupled','components','#7B3294',':')]: rows=[r for r in scaling if r['family']==family and r['mode']==mode] ax.plot([r['size'] for r in rows],[r['resolutions_inspected'] for r in rows],marker='o',linestyle=style,color=color,label=f'{family.capitalize()}, {mode}') ax.set(xlabel='Number of modules or gates',ylabel='Supplier resolutions inspected',yscale='log',xticks=list(SCALING_SIZES)) ax.legend(fontsize=9);ax.grid(axis='y',color='#e3e6e8') fig.savefig(output/'scaling.png',dpi=220);fig.savefig(output/'scaling.svg');plt.close(fig) def main(): parser=argparse.ArgumentParser(description=__doc__) parser.add_argument('--output',type=Path,default=Path('outputs')) parser.add_argument('--input',type=Path,help='Optional JSON food/reactions model; see model.json.') args=parser.parse_args();args.output.mkdir(parents=True,exist_ok=True);start=time.perf_counter() system=ReactionSystem.from_record(json.loads(args.input.read_text())) if args.input else ReactionSystem(FOOD,[Reaction(*r) for r in REACTIONS]) inputs=dict(model=system.record(),resolution_budget=RESOLUTION_BUDGET,availability=AVAILABLE_PROBABILITY, deletion_costs=DELETION_COSTS,required=REQUIRED_REACTIONS,forbidden=FORBIDDEN_REACTIONS,deleted=DELETED_REACTIONS, scaling_sizes=SCALING_SIZES,large_modules=LARGE_INDEPENDENT_MODULES,chain_length=DETERMINISTIC_CHAIN_LENGTH) intro='Resolved inputs: '+json.dumps(inputs);print(intro,flush=True) engine=SupplierEnumerator();catalogue=engine.catalogue(system) summary=dict(inputs=inputs,catalogue=family_record(catalogue)) if not args.input: cuts=minimal_cuts(system,catalogue);chosen,cost=greedy_cut(catalogue,DELETION_COSTS) summary['worked_queries']=dict(required=list(REQUIRED_REACTIONS),forbidden=list(FORBIDDEN_REACTIONS), matching=[sorted(c) for c in catalogue.filter(REQUIRED_REACTIONS,FORBIDDEN_REACTIONS)], after_deletion=[sorted(c) for c in catalogue.filter(forbidden=DELETED_REACTIONS)], minimal_cuts=[sorted(c) for c in cuts],minimum_cut_cost=min(sum(DELETION_COSTS[r] for r in c) for c in cuts), greedy_cut=chosen,greedy_cost=str(cost),availability_probability=str(exact_availability(system,catalogue))) # Paper's catalyst-alternative example: resolved minimum is not original minimum. counterexample=ReactionSystem({'f'},[Reaction(0,{'f'},{'x'},{'x','y'}),Reaction(1,{'x'},{'y'},{'x'})]) summary['necessary_validation']=family_record(engine.global_catalogue(counterexample)) scaling=[] for size in SCALING_SIZES: for family,builder in [('independent',independent_modules),('coupled',coupled_gates)]: model=builder(size) for mode,enumerate_ in [('global',engine.global_catalogue),('components',engine.catalogue)]: result=enumerate_(model) scaling.append(dict(family=family,mode=mode,size=size,reactions=len(model.identifiers), beta=result.beta,global_resolution_count=result.global_resolution_count, resolutions_inspected=result.resolutions_inspected,outputs=len(result.members), validation_calls=result.validation_calls,total_memberships=sum(map(len,result.members)))) large=engine.catalogue(independent_modules(LARGE_INDEPENDENT_MODULES)) summary['large_independent']=family_record(large) chain=deterministic_chain(DETERMINISTIC_CHAIN_LENGTH) summary['deterministic_chain']=dict(**family_record(engine.global_catalogue(chain)),closure_depth=len(chain.closure_stages())-1) zero=engine.catalogue(independent_modules(4,choices=False)) summary['zero_excess_availability']=str(zero_excess_reliability(zero)) summary['evidence']='Exact finite set computations and rational availability probabilities. Bounded independent tests, not Lean extraction or compilation. Structural RAF membership does not assert kinetic function.' with (args.output/'scaling.csv').open('w',newline='',encoding='utf-8') as f: writer=csv.DictWriter(f,fieldnames=list(scaling[0]));writer.writeheader();writer.writerows(scaling) with (args.output/'catalogue.csv').open('w',newline='',encoding='utf-8') as f: writer=csv.writer(f);writer.writerow(['irrRAF','reaction_ids']);writer.writerows((i+1,' '.join(map(str,sorted(c)))) for i,c in enumerate(catalogue.members)) (args.output/'model.json').write_text(json.dumps(system.record(),indent=2)+'\n',encoding='utf-8') figures(system,catalogue,scaling,args.output) 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(), matplotlib=matplotlib.__version__, python=platform.python_version(),platform=platform.platform(),processor=platform.processor(), elapsed_seconds=time.perf_counter()-start,command='python example.py --output outputs', seed_policy='Deterministic enumeration; no stochastic 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()