Example code

"""Find productive autocatalytic cores and independently verify their certificates."""
# EDIT THESE INPUTS. Rows are species; columns are reversible reaction identities.
SPECIES=('A','B','C')
REACTIONS=('A_to_2B','B_to_2C','C_to_2A')
LEFT=((1,0,0),(0,1,0),(0,0,1))
RIGHT=((0,0,2),(2,0,0),(0,2,0))
SEARCH_CANDIDATE_BUDGET=100000
CNF_CLAUSES=((1,),(-1,2))       # signed one-based literals; repeated occurrences retained
CNF_VARIABLES=2
ASSIGNMENT_BUDGET=65536

import argparse,hashlib,json,platform,csv
from pathlib import Path
from itertools import combinations
import numpy as np
import sympy as sp
from network import LiteralSource,SquareSearch,BinarySourceCodec
from linkage import CNF,DirectedLinkage,switch_check
from check_certificate import check
MANUSCRIPT_SHA256='a3ac0ec8661bb0553397f5bc881993d054eeb3b6917f2dcdec4569b0d87bf590'

def examples():
    return {
        'literal_self_amplification':LiteralSource(('A',),('A_to_2A',),((1,),),((2,),)),
        'reverse_flow_needed':LiteralSource(('A',),('2A_to_A',),((2,),),((1,),)),
        'balanced_conversion':LiteralSource(('A','B'),('A_to_B',),((1,),(0,)),((0,),(1,))),
        'external_resource':LiteralSource(('A','F'),('A_plus_F_to_2A',),((1,),(1,)),((2,),(0,))),
        'net_only_cancellation_error':LiteralSource(('A',),('cancelled',),((0,),),((1,),)),
        'closed_conservative_cycle':LiteralSource(('A','B','C'),('A_B','B_C','C_A'),((1,0,0),(0,1,0),(0,0,1)),((0,0,1),(1,0,0),(0,1,0)))
    }

def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));out=parser.parse_args().output;out.mkdir(parents=True,exist_ok=True)
    def dump(name,x):(out/name).write_text(json.dumps(x,indent=2,default=str)+'\n')
    source=LiteralSource(SPECIES,REACTIONS,LEFT,RIGHT);result=SquareSearch(SEARCH_CANDIDATE_BUDGET).search(source);dump('configured_search.json',result)
    if result['status']=='yes':dump('configured_certificate.json',result['witness']['certificate'])
    encoded=BinarySourceCodec.encode(source);(out/'configured_source.bits.txt').write_text(encoded+'\n');decoded=BinarySourceCodec.decode(encoded)
    if decoded.left!=source.left or decoded.right!=source.right:raise ArithmeticError('Literal binary round trip failed.')
    rows=[];details={}
    for name,s in examples().items():
        answer=SquareSearch(SEARCH_CANDIDATE_BUDGET).search(s);details[name]=answer;rows.append(dict(example=name,status=answer['status'],candidates=answer['attempted'],selected_entities=answer.get('witness',{}).get('certificate',{}).get('entities'),signed_flow=answer.get('witness',{}).get('certificate',{}).get('flow')))
    dump('literal_and_signed_flow_examples.json',details)
    dump('explicit_budget_exhaustion.json',SquareSearch(0).search(source))
    fresh_switch=switch_check();dump('complete_switch_check.json',fresh_switch)
    pairing={}
    for label,arcs in [('desired',(('a0','b0'),('a1','b1'))),('crossed',(('a0','b1'),('a1','b0')))]:
        g=DirectedLinkage(('a0','a1','b0','b1'),arcs,('a0','a1','b0','b1'));s=g.source();answer=SquareSearch(SEARCH_CANDIDATE_BUDGET).search(s);pairing[label]=dict(net=s.net,boundary_determinant=int(sp.Matrix(s.net).det()),search=answer)
    dump('desired_versus_crossed.json',pairing)
    census=[];arcs=(('a0','b0'),('a0','b1'),('a1','b0'),('a1','b1'))
    for mask in range(16):
        selected=tuple(e for i,e in enumerate(arcs) if mask&(1<<i));g=DirectedLinkage(('a0','a1','b0','b1'),selected,('a0','a1','b0','b1'));answer=SquareSearch(1000).search(g.source());expected=bool(mask&1 and mask&8)
        if (answer['status']=='yes')!=expected:raise ArithmeticError('Small complete linkage/PAC equivalence check failed.')
        census.append(dict(mask=mask,desired_paths=expected,pac_status=answer['status']))
    dump('all_direct_linkage_graphs.json',census)
    cnf=CNF(CNF_CLAUSES,CNF_VARIABLES);graph=cnf.compile();converted=graph.source();sat=cnf.assignments(ASSIGNMENT_BUDGET)
    formula=dict(clauses=cnf.clauses,variables=cnf.variable_count,occurrences=sum(map(len,cnf.clauses)),vertices=len(graph.vertices),arcs=len(graph.arcs),source_entities=len(converted.entities),source_reactions=len(converted.reactions),max_stoichiometry=max([0]+[v for matrix in (converted.left,converted.right) for row in matrix for v in row]),assignment_search=sat)
    dump('formula_source.json',dict(entities=converted.entities,reactions=converted.reactions,left=converted.left,right=converted.right))
    if sat['status']=='sat':
        paths=cnf.satisfying_paths(sat['assignment']);certificate=graph.path_certificate(paths);dump('formula_productivity_certificate.json',certificate);dump('formula_paths.json',dict(P=paths[0],Q=paths[1]));formula['certificate_verified']=certificate['verification']['accepted'];formula['square_size']=len(certificate['certificate']['entities']);formula['integer_flow_bits']=certificate['max_flow_bits']
    else:formula['scope']='No direct PAC search on the large compiled network. UNSAT implies no PAC by the paper reduction theorem; budget exhaustion implies neither.'
    dump('sat_to_network_summary.json',formula)
    boundary=[]
    for clauses,V in [((),0),(((),),0),(((1,),(-1,)),1),(((1,1),),1)]:
        f=CNF(clauses,V);g=f.compile();answer=f.assignments();record=dict(clauses=clauses,variables=V,status=answer['status'],terminals=g.terminals,distinct_terminals=len(set(g.terminals))==4)
        if answer['status']=='sat':record['certificate_verified']=g.path_certificate(f.satisfying_paths(answer['assignment']))['verification']['accepted']
        else:record['no_pac_scope']='Inference from the published reduction, not exhaustive support enumeration of this large instance.'
        boundary.append(record)
    dump('formula_boundary_cases.json',boundary)
    damaged=None
    if result['status']=='yes':
        c=json.loads(json.dumps(result['witness']['certificate']));c['flow']=[-v for v in c['flow']];damaged=check(c)
    dump('rejected_reversed_certificate.json',damaged)
    plot(out,source,result,pairing,rows)
    print('Configured network:',result['status'],'after',result['attempted'],'square candidates.',flush=True)
    print('Desired/crossed boundary determinants:',pairing['desired']['boundary_determinant'],pairing['crossed']['boundary_determinant'],flush=True)
    print('Complete 20-vertex switch replay:',fresh_switch,flush=True)
    print('Compiled formula:',formula,flush=True)
    print('Certificates establish literal structural productivity and PAC existence. They do not prescribe concentrations, kinetics or conserved elemental chemistry. Exhaustive NO and budget-limited UNKNOWN are distinct. Lean and machine-time proofs are not rerun.',flush=True)
    here=Path(__file__).parent;digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,source_sha256=digest(Path(__file__)),python=platform.python_version(),sympy=sp.__version__,numpy=np.__version__,module_sha256={p.name:digest(p) for p in sorted(here.glob('*.py')) if p.name not in ['example.py','test_example.py']},output_sha256={p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}))

def plot(out,source,result,pairing,rows):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,3,figsize=(11,4),layout='constrained')
    for ax,matrix,title in zip(axs,[source.left,source.right,source.net],['Literal reactants L','Literal products P','Net change N = P - L']):
        a=np.array(matrix);ax.imshow(a,cmap='coolwarm',vmin=-max(1,np.max(abs(a))),vmax=max(1,np.max(abs(a))))
        for (i,j),v in np.ndenumerate(a):ax.text(j,i,str(v),ha='center',va='center',color='black')
        ax.set(xticks=range(len(source.reactions)),xticklabels=[f'r{i}' for i in range(len(source.reactions))],yticks=range(len(source.entities)),yticklabels=source.entities,title=title,xlabel='Reversible reaction identity')
    fig.savefig(out/'literal_network.png',dpi=180);fig.savefig(out/'literal_network.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    for ax,(name,d) in zip(axs,pairing.items()):
        matrix=np.array(d['net']).T;ax.imshow(matrix,cmap='coolwarm',vmin=-4,vmax=4)
        for (i,j),v in np.ndenumerate(matrix):ax.text(j,i,str(v),ha='center',va='center',fontsize=20)
        ax.set(xticks=[0,1],xticklabels=['Path from a0','Path from a1'],yticks=[0,1],yticklabels=['Source boundary','Sink boundary'],title=f'{name.capitalize()} pairing: determinant {d["boundary_determinant"]}')
    fig.savefig(out/'pairing_obstruction.png',dpi=180);fig.savefig(out/'pairing_obstruction.svg');plt.close(fig)

if __name__=='__main__':main()