Example code

"""Build literal RAF realizations from feasible families and catalyst graphs.

The paper's marker construction and paired construction; finite exact checks.
Run python example.py --output outputs. Importing this module performs no work.
"""
from __future__ import annotations
import argparse
import csv
from dataclasses import dataclass
import hashlib
from itertools import combinations
import json
from pathlib import Path
import platform
import time

# USER INPUTS ---------------------------------------------------------------
ELEMENTS = ('a', 'b', 'c')
FEASIBLE_SETS = ((), ('a',), ('a','b'), ('a','b','c'))
CATALYST_ARCS = (('a','a'), ('a','c'), ('c','b'))  # producer -> catalysed reaction
USE_MINIMAL_BLOCKERS = False  # True uses the paper's unformalized reduction remark
CHAIN_SIZES = (1,2,3,4,5,6,7,8)
PAIRED_THRESHOLD_SIZE = 5
THRESHOLD = 3
EXHAUSTIVE_ANTIMATROID_SIZE = 4  # all labelled antimatroids through this size
EXHAUSTIVE_GRAPH_SIZE = 3       # all graphs for the antimatroids through this size
MAX_GROUND_SIZE = 8            # explicit construction/enumeration budgets
MAX_REACTION_ENUMERATION = 16
PAPER_SHA256 = '0b534fbc24da9d3482699977e9f464b3fc7a5e5b48e1ae46a720b9ff4f9bcd44'
# Inputs are finite sets and directed incidences, without kinetic/physical units.
# --------------------------------------------------------------------------


def subsets(elements):
    elements=tuple(elements)
    for size in range(len(elements)+1):
        for selected in combinations(elements,size): yield frozenset(selected)


def ordered(family):
    return sorted(family,key=lambda s:(len(s),sorted(s)))


def label(s):
    return '{'+','.join(sorted(s))+'}' if s else 'empty'


class FixedFamily:
    """Union-closed fixed sets, including empty; owns the interior operator."""
    def __init__(self,elements,members):
        elements=tuple(elements)
        if len(set(elements))!=len(elements) or any(not isinstance(e,str) or not e for e in elements):
            raise ValueError('Elements must be distinct nonempty strings.')
        if len(elements)>MAX_GROUND_SIZE: raise ValueError('Ground-set budget exceeded.')
        self.elements=frozenset(elements);self.members=frozenset(frozenset(s) for s in members)
        if frozenset() not in self.members or any(not s<=self.elements for s in self.members):
            raise ValueError('Members must be subsets of the ground set and include empty.')
        if any(a|b not in self.members for a in self.members for b in self.members):
            raise ValueError('Fixed family must be union-closed.')

    def interior(self,s):
        s=frozenset(s)
        if not s<=self.elements: raise ValueError('Unknown element.')
        return frozenset().union(*(a for a in self.members if a<=s))

    def blockers(self,e,minimal=False):
        if e not in self.elements: raise ValueError('Unknown blocker target.')
        result=tuple(b for b in subsets(sorted(self.elements-{e})) if e not in self.interior(self.elements-b))
        return tuple(b for b in result if not any(a<b for a in result)) if minimal else result

    def accessible(self):
        return all(not s or any(s-{e} in self.members for e in s) for s in self.members)


class Antimatroid(FixedFamily):
    def __init__(self,elements,members):
        super().__init__(elements,members)
        if not self.accessible(): raise ValueError('An antimatroid must be accessible.')


class Digraph:
    """Loops permitted; arc direction is catalyst producer -> consumer."""
    def __init__(self,elements,arcs):
        self.elements=frozenset(elements);self.arcs=frozenset(tuple(a) for a in arcs)
        if any(len(a)!=2 or not set(a)<=self.elements for a in self.arcs):
            raise ValueError('Graph arcs must have two known endpoints.')
        self.predecessors={e:frozenset(u for u,v in self.arcs if v==e) for e in self.elements}

    def supported(self,s):
        s=frozenset(s)
        if not s<=self.elements: raise ValueError('Unknown graph vertex.')
        return all(self.predecessors[e]&s for e in s)


@dataclass(frozen=True)
class Reaction:
    identifier: str
    reactants: frozenset[str]
    products: frozenset[str]
    catalysts: frozenset[str]

    def __post_init__(self):
        for field in ('reactants','products','catalysts'):
            object.__setattr__(self,field,frozenset(getattr(self,field)))


class ReactionSystem:
    """Closure uses reactants alone. Catalysis is checked after closure."""
    def __init__(self,reactions,food=('f',)):
        reactions=tuple(reactions);self.reactions={r.identifier:r for r in reactions}
        if len(reactions)!=len(self.reactions): raise ValueError('Duplicate reaction identifiers.')
        self.food=frozenset(food);self.identifiers=frozenset(self.reactions)

    def selected(self,s):
        s=frozenset(s)
        if not s<=self.identifiers: raise ValueError('Unknown reaction identifier.')
        return s

    def closure_stages(self,s):
        s=self.selected(s);stages=[self.food]
        while True:
            current=stages[-1]
            following=current.union(*(self.reactions[e].products for e in s if self.reactions[e].reactants<=current))
            if following==current: return stages
            stages.append(following)

    def food_generated(self,s):
        s=self.selected(s);closure=self.closure_stages(s)[-1]
        return all(self.reactions[e].reactants<=closure for e in s)

    def is_raf(self,s):
        s=self.selected(s);closure=self.closure_stages(s)[-1]
        return bool(s) and all(self.reactions[e].reactants<=closure and self.reactions[e].catalysts&closure for e in s)

    def maximal_raf(self,s):
        s=self.selected(s)
        while s:
            closure=self.closure_stages(s)[-1]
            kept=frozenset(e for e in s if self.reactions[e].reactants<=closure and self.reactions[e].catalysts&closure)
            if kept==s: return s
            s=kept
        return s

    def fixed_family(self):
        if len(self.identifiers)>MAX_REACTION_ENUMERATION: raise ValueError('Reaction-subset enumeration budget exceeded.')
        return frozenset(s for s in subsets(sorted(self.identifiers)) if not s or self.is_raf(s))

    def product_graph(self):
        arcs={(u,e) for u,r in self.reactions.items() for e,t in self.reactions.items() if r.products&t.catalysts}
        arcs|={(e,e) for e,r in self.reactions.items() if self.food&r.catalysts}
        return Digraph(self.identifiers,arcs)

    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()])


class MarkerInventory:
    """Shared mechanical encoding; names use indices to avoid label collisions."""
    def __init__(self,family,minimal=False):
        self.elements=sorted(family.elements);self.index={e:i for i,e in enumerate(self.elements)}
        self.blockers={e:family.blockers(e,minimal) for e in self.elements}
        self.records=[(e,b,f'm:{self.index[e]}:{k}') for e in self.elements for k,b in enumerate(self.blockers[e])]

    def required(self,e): return frozenset(m for t,b,m in self.records if t==e)
    def produced(self,e): return frozenset(m for t,b,m in self.records if e in b)


class MarkerRealization:
    def __init__(self,feasible:Antimatroid,graph:Digraph,minimal=False):
        if not isinstance(feasible,Antimatroid): raise ValueError('Same-ground construction requires an Antimatroid.')
        if feasible.elements!=graph.elements: raise ValueError('Family and graph must use the same ground set.')
        self.feasible,self.graph=feasible,graph;self.markers=MarkerInventory(feasible,minimal)
        private={e:f'c:{i}' for e,i in self.markers.index.items()}
        self.system=ReactionSystem(Reaction(e,{'f'}|self.markers.required(e),{private[e]}|self.markers.produced(e),
                    {private[u] for u in graph.predecessors[e]}) for e in self.markers.elements)

    def audit(self):
        expected=FixedFamily(self.feasible.elements,(s for s in self.feasible.members if self.graph.supported(s)))
        rows=[]
        for s in subsets(sorted(self.feasible.elements)):
            fg=self.system.food_generated(s);support=self.graph.supported(s);raf=self.system.is_raf(s)
            maximum=self.system.maximal_raf(s)
            if fg!=(s in self.feasible.members) or (not s or raf)!=(s in expected.members) or maximum!=expected.interior(s):
                raise ArithmeticError('Literal realization differs from the two-layer specification.')
            rows.append(dict(selected=label(s),food_generated=fg,catalyst_supported=support,is_raf=raf,
                fixed_set=(not s or raf),maximal_raf=label(maximum),closure_stages=len(self.system.closure_stages(s))-1))
        if self.system.product_graph().arcs!=self.graph.arcs: raise ArithmeticError('Product catalyst graph mismatch.')
        return rows


class PairedRealization:
    def __init__(self,family:FixedFamily,minimal=False):
        self.family=family;self.markers=MarkerInventory(family,minimal)
        self.pairs={e:(f'p:{i}',f'g:{i}') for e,i in self.markers.index.items()};reactions=[]
        for e,i in self.markers.index.items():
            p,g=self.pairs[e];alpha,beta=f'alpha:{i}',f'beta:{i}'
            reactions.extend([Reaction(p,{'f'},{alpha}|self.markers.produced(e),{beta}),
                Reaction(g,{'f'}|self.markers.required(e),{beta},{alpha})])
        self.system=ReactionSystem(reactions)

    def encode(self,s):
        s=frozenset(s)
        if not s<=self.family.elements: raise ValueError('Unknown element.')
        return frozenset(r for e in s for r in self.pairs[e])

    def audit(self):
        actual=self.system.fixed_family();expected=frozenset(self.encode(s) for s in self.family.members)
        if actual!=expected: raise ArithmeticError('Paired construction has missing or spurious fixed sets.')
        # Includes incomplete pairs: projection means both partner reactions available.
        for t in subsets(sorted(self.system.identifiers)):
            available={e for e,pair in self.pairs.items() if set(pair)<=t}
            if self.system.maximal_raf(t)!=self.encode(self.family.interior(available)):
                raise ArithmeticError('Paired maxRAF differs on a restricted reaction set.')
        return dict(elements=len(self.family.elements),reactions=len(self.system.identifiers),
            reaction_subsets_checked=2**len(self.system.identifiers),nonempty_rafs=len(actual)-1,
            marker_molecules=len(self.markers.records))


def threshold_family(n,k):
    if not isinstance(n,int) or not 0<=n<=MAX_GROUND_SIZE or not isinstance(k,int) or k<1:
        raise ValueError('Threshold family requires a bounded nonnegative size and positive threshold.')
    elements=tuple(str(i) for i in range(n))
    return FixedFamily(elements,(s for s in subsets(elements) if not s or len(s)>=k))


def four_element_certificate():
    elements=('a','b','c','d')
    small=((),('a',),('b',),('a','b'),('a','c'),('b','c'))
    feasible=Antimatroid(elements,(*small,*(s for s in subsets(elements) if len(s)>=3)))
    graph=Digraph(elements,(('b','a'),('d','a'),('c','b'),('d','b'),('a','c'),('d','c'),('a','d'),('b','d')))
    return MarkerRealization(feasible,graph)


def all_antimatroids(n):
    if not isinstance(n,int) or not 0<=n<=4: raise ValueError('All-family enumeration is capped at four elements.')
    elements=tuple(str(i) for i in range(n));universe=tuple(subsets(elements))
    for mask in range(1<<(len(universe)-1)):
        family=frozenset([frozenset(),*(s for i,s in enumerate(universe[1:]) if mask>>i&1)])
        if any(s and not any(s-{e} in family for e in s) for s in family): continue
        if any(a|b not in family for a in family for b in family): continue
        yield Antimatroid(elements,family)


def exhaustive_checks(max_antimatroid=4,max_graph=3):
    if not 0<=max_graph<=min(3,max_antimatroid) or not 0<=max_antimatroid<=4:
        raise ValueError('Exhaustive checks permit antimatroids through 4 and graphs through 3.')
    rows=[]
    for n in range(max_antimatroid+1):
        count=fg_checks=raf_checks=0
        for feasible in all_antimatroids(n):
            count+=1;elements=sorted(feasible.elements)
            base=MarkerRealization(feasible,Digraph(elements,()))
            selections=tuple(subsets(elements))
            for s in selections:
                if base.system.food_generated(s)!=(s in feasible.members): raise ArithmeticError('Food family mismatch.')
                fg_checks+=1
            if n<=max_graph:
                arcs=tuple((u,v) for u in elements for v in elements)
                for mask in range(1<<len(arcs)):
                    graph=Digraph(elements,(a for i,a in enumerate(arcs) if mask>>i&1))
                    model=MarkerRealization(feasible,graph)
                    for s in selections:
                        if (not s or model.system.is_raf(s))!=((s in feasible.members) and graph.supported(s)):
                            raise ArithmeticError('Literal RAF factorization mismatch.')
                        raf_checks+=1
        rows.append(dict(elements=n,antimatroids=count,food_comparisons=fg_checks,raf_comparisons=raf_checks))
    return rows


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(truth,scaling,output):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    from matplotlib.colors import ListedColormap
    with plt.rc_context({'font.size':11,'figure.facecolor':'white','savefig.facecolor':'white','svg.fonttype':'none'}):
        fig,ax=plt.subplots(figsize=(7.2,4.5),layout='constrained')
        data=[[int(r[c]) for c in ('food_generated','catalyst_supported','fixed_set')] for r in truth]
        ax.imshow(data,cmap=ListedColormap(['#e8edf0','#0072B2']),vmin=0,vmax=1,aspect='auto')
        ax.set_xticks(range(3),['Food generated','Catalyst supported','Fixed set (intersection)'])
        ax.set_yticks(range(len(truth)),[r['selected'] for r in truth])
        for i,row in enumerate(data):
            for j,v in enumerate(row): ax.text(j,i,'yes' if v else 'no',ha='center',va='center',color='white' if v else '#30363a')
        ax.tick_params(length=0);fig.savefig(output/'layers.png',dpi=220);fig.savefig(output/'layers.svg');plt.close(fig)
        fig,ax=plt.subplots(figsize=(7.2,4.5),layout='constrained')
        ax.plot([r['elements'] for r in scaling],[r['all_blockers'] for r in scaling],'o-',label='All blockers (main construction)',color='#0072B2')
        ax.plot([r['elements'] for r in scaling],[r['minimal_blockers'] for r in scaling],'s--',label='Minimal blockers (paper remark)',color='#D55E00')
        ax.set(xlabel='Length of the chain antimatroid',ylabel='Auxiliary marker molecules')
        ax.spines[['top','right']].set_visible(False);ax.grid(axis='y',color='#e3e6e8');ax.legend(fontsize=9)
        fig.savefig(output/'markers.png',dpi=220);fig.savefig(output/'markers.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()
    inputs=dict(elements=ELEMENTS,feasible_sets=FEASIBLE_SETS,catalyst_arcs=CATALYST_ARCS,minimal_blockers=USE_MINIMAL_BLOCKERS,
        chain_sizes=CHAIN_SIZES,paired_threshold_size=PAIRED_THRESHOLD_SIZE,threshold=THRESHOLD,
        exhaustive_antimatroid_size=EXHAUSTIVE_ANTIMATROID_SIZE,exhaustive_graph_size=EXHAUSTIVE_GRAPH_SIZE,
        max_ground_size=MAX_GROUND_SIZE,max_reaction_enumeration=MAX_REACTION_ENUMERATION)
    intro='Resolved inputs: '+json.dumps(inputs);print(intro,flush=True)
    model=MarkerRealization(Antimatroid(ELEMENTS,FEASIBLE_SETS),Digraph(ELEMENTS,CATALYST_ARCS),USE_MINIMAL_BLOCKERS)
    truth=model.audit();write_csv(args.output/'layers.csv',truth)
    blockers=[dict(target=e,blocker=label(b),marker=m) for e,b,m in model.markers.records]
    (args.output/'blockers.json').write_text(json.dumps(blockers,indent=2)+'\n',encoding='utf-8')
    (args.output/'network.json').write_text(json.dumps(model.system.record(),indent=2)+'\n',encoding='utf-8')
    write_csv(args.output/'closure.csv',[dict(stage=i,molecules=' '.join(sorted(s))) for i,s in enumerate(model.system.closure_stages(ELEMENTS))])
    certificate=four_element_certificate();certificate.audit()
    if any(bool(s) and (certificate.system.is_raf(s)!=(len(s)>=3)) for s in subsets(sorted(certificate.feasible.elements))):
        raise ArithmeticError('Four-element threshold certificate failed.')
    paired=PairedRealization(threshold_family(PAIRED_THRESHOLD_SIZE,THRESHOLD),USE_MINIMAL_BLOCKERS);paired_audit=paired.audit()
    (args.output/'paired_network.json').write_text(json.dumps(paired.system.record(),indent=2)+'\n',encoding='utf-8')
    scaling=[]
    for n in CHAIN_SIZES:
        elements=tuple(str(i) for i in range(n));family=Antimatroid(elements,(elements[:i] for i in range(n+1)))
        all_model=MarkerRealization(family,Digraph(elements,((e,e) for e in elements)))
        minimal_model=MarkerRealization(family,all_model.graph,True)
        all_model.audit();minimal_model.audit()
        scaling.append(dict(elements=n,reactions=n,all_blockers=len(all_model.markers.records),minimal_blockers=len(minimal_model.markers.records)))
    write_csv(args.output/'marker_scaling.csv',scaling)
    checks=exhaustive_checks(EXHAUSTIVE_ANTIMATROID_SIZE,EXHAUSTIVE_GRAPH_SIZE);write_csv(args.output/'exhaustive_checks.csv',checks)
    plot(truth,scaling,args.output)
    summary=dict(inputs=inputs,worked_fixed_sets=[label(s) for s in ordered(model.system.fixed_family())],
        worked_markers=len(model.markers.records),four_element_threshold_rafs=len(certificate.system.fixed_family())-1,
        paired_threshold=paired_audit,marker_scaling=scaling,exhaustive_checks=checks,
        total_food_comparisons=sum(r['food_comparisons'] for r in checks),total_raf_comparisons=sum(r['raf_comparisons'] for r in checks),
        evidence='Exact finite constructions and exhaustive subset checks. Five-element same-ground impossibility is the manuscript theorem, not searched here. 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()