Example code

"""Static thermodynamic compatibility, exact triangle synthesis, and obstructions.

Linear programs propose witnesses/certificates; Fraction arithmetic checks them.
"""
from __future__ import annotations
import argparse
import csv
from dataclasses import dataclass
from fractions import Fraction as F
from itertools import combinations
import hashlib
import json
from math import log
from pathlib import Path
import platform
import time
import numpy as np
from scipy.optimize import linprog

# USER INPUTS ---------------------------------------------------------------
FACTORS = ('1','1','1','1','1') # shared, left1, left2, right1, right2
GAIN = 2                     # same product complex on both branches
A_ACTIVITY = '1/2'            # free construction choice, dimensionless, 0 < A < 1
K_SWEEP = tuple(F(i,100) for i in range(20,301,2))
TEMPERATURE_K = 298.0         # optional symmetric equal-prefactor energy reading
GAS_CONSTANT_J_MOL_K = 8.314462618
RATIONAL_DENOMINATOR_LIMIT = 1_000_000
MAX_PAC_ENUMERATION_BITS = 16 # species + reactions of a motif
PAPER_SHA256 = 'fcd130caf5eb188887d31e607de61c6ad6a65053acf1c6ed3a62de7b2e9006d0'
# Activities and factors use the paper's normalized current convention.
# These are static productive-current witnesses, not steady operating states.
# --------------------------------------------------------------------------


def matrix(rows): return np.array([[F(v) for v in row] for row in rows],dtype=object)
def subsets(values):
    values=tuple(values)
    for size in range(len(values)+1):
        for s in combinations(values,size): yield tuple(s)


@dataclass(frozen=True)
class LinearDecision:
    status: str
    vector: tuple[F,...]
    residuals: tuple[F,...]

    def record(self):
        return dict(status=self.status,vector=list(map(str,self.vector)),residuals=list(map(str,self.residuals)))


class StrictLinearFeasibility:
    """Solve A x > 0 over free x; return only exactly verified conclusions.

    Homogeneity allows A x >= 1. Infeasibility uses w>=0, sum(w)=1, A.T w=0.
    """
    @staticmethod
    def solve(rows):
        a=matrix(rows)
        if a.ndim!=2 or min(a.shape)<1: raise ValueError('A nonempty inequality matrix is required.')
        numeric=np.array(a,dtype=float);nr,nv=a.shape
        result=linprog(np.zeros(nv),A_ub=-numeric,b_ub=-np.ones(nr),bounds=[(None,None)]*nv,method='highs')
        if result.success:
            witness=tuple(F(float(v)).limit_denominator(RATIONAL_DENOMINATOR_LIMIT) for v in result.x)
            residuals=tuple(a@witness)
            if all(v>0 for v in residuals): return LinearDecision('feasible',witness,residuals)
        dual=linprog(np.zeros(nr),A_eq=np.vstack([numeric.T,np.ones(nr)]),b_eq=np.r_[np.zeros(nv),1],bounds=[(0,None)]*nr,method='highs')
        if dual.success:
            certificate=tuple(F(float(v)).limit_denominator(RATIONAL_DENOMINATOR_LIMIT) for v in dual.x)
            residuals=tuple(a.T@certificate)
            if all(w>=0 for w in certificate) and sum(certificate)>0 and all(v==0 for v in residuals):
                return LinearDecision('infeasible',certificate,residuals)
        return LinearDecision('unresolved',(),())


@dataclass(frozen=True)
class Reaction:
    reactants: tuple[int,...]
    products: tuple[int,...]
    factor: F

    def __post_init__(self):
        object.__setattr__(self,'reactants',tuple(self.reactants));object.__setattr__(self,'products',tuple(self.products))
        object.__setattr__(self,'factor',F(self.factor))
        if len(self.reactants)!=len(self.products) or not self.reactants or any(not isinstance(v,int) or v<0 for v in self.reactants+self.products):
            raise ValueError('Complexes require equally sized nonnegative integer vectors.')
        if self.factor<=0: raise ValueError('Factors must be strictly positive.')


@dataclass(frozen=True)
class Motif:
    species: tuple[int,...]
    reactions: tuple[int,...]


class CoreFamily:
    """Literal network and motif balances; PAC status can be audited separately."""
    def __init__(self,species,reactions,motifs,orientation=None):
        self.species=tuple(species);self.reactions=tuple(reactions);self.motifs=tuple(motifs)
        n,m=len(self.species),len(self.reactions)
        if n==0 or m==0 or len(set(self.species))!=n or not self.motifs: raise ValueError('Nonempty network and distinct species required.')
        if any(len(r.reactants)!=n for r in self.reactions): raise ValueError('Complex dimension mismatch.')
        self.orientation=tuple(orientation if orientation is not None else [1]*m)
        if len(self.orientation)!=m or any(e not in (-1,1) for e in self.orientation): raise ValueError('Orientation entries must be +1 or -1.')
        if any(not c.species or not c.reactions or len(set(c.species))!=len(c.species) or len(set(c.reactions))!=len(c.reactions)
               or not set(c.species)<=set(range(n)) or not set(c.reactions)<=set(range(m)) for c in self.motifs):
            raise ValueError('Invalid motif indices.')
        if set().union(*(set(c.reactions) for c in self.motifs))!=set(range(m)): raise ValueError('Every reaction must belong to a motif.')
        self.stoichiometry=matrix([[r.products[i]-r.reactants[i] for r in self.reactions] for i in range(n)])
        self.complexes=tuple(sorted({c for r in self.reactions for c in (r.reactants,r.products)}))
        self.incidence=matrix([[r.factor*(int(c==r.reactants)-int(c==r.products)) for c in self.complexes] for r in self.reactions])

    def balance_rows(self,motifs=None):
        return matrix([[self.stoichiometry[i,j] if j in c.reactions else 0 for j in range(len(self.reactions))]
                       for c in (self.motifs if motifs is None else motifs) for i in c.species])

    def flow_rows(self): return np.vstack([np.diag(self.orientation),self.balance_rows()])
    def direction_rows(self): return -np.diag(self.orientation)@self.stoichiometry.T
    def linear_complex_rows(self): return np.vstack([self.flow_rows()@self.incidence,np.eye(len(self.complexes),dtype=int)])
    def filters(self):
        return {name:StrictLinearFeasibility.solve(a) for name,a in [('productive_flow',self.flow_rows()),('direction',self.direction_rows()),('independent_complex',self.linear_complex_rows())]}

    def currents(self,activities):
        z=tuple(F(v) for v in activities)
        if len(z)!=len(self.species) or any(v<=0 for v in z): raise ValueError('A positive activity per species is required.')
        monomials=[]
        for c in self.complexes:
            value=F(1)
            for activity,power in zip(z,c): value*=activity**power
            monomials.append(value)
        return self.incidence@monomials

    def evaluate(self,activities):
        currents=self.currents(activities);residuals=self.balance_rows()@currents;directions=np.array(self.orientation)*currents
        return dict(activities=list(map(str,activities)),currents=list(map(str,currents)),
            motif_production=list(map(str,residuals)),whole_network_production=list(map(str,self.stoichiometry@currents)),
            strictly_productive=bool(all(v>0 for v in (*residuals,*directions))))

    def side_incident(self,motif):
        return bool(motif.species and motif.reactions) and all(any(self.reactions[j].reactants[i]>0 for i in motif.species)
            and any(self.reactions[j].products[i]>0 for i in motif.species) for j in motif.reactions)

    def audit_pac(self,motif):
        """Finite submotif check; currents can have either sign, as in the paper."""
        if len(motif.species)+len(motif.reactions)>MAX_PAC_ENUMERATION_BITS:
            raise ValueError('PAC submotif enumeration exceeds declared budget.')
        if not self.side_incident(motif): return dict(is_pac=False,reason='not side-incident')
        own=StrictLinearFeasibility.solve(self.balance_rows([motif])[:,motif.reactions])
        if own.status!='feasible': return dict(is_pac=False if own.status=='infeasible' else None,reason=own.status)
        checked=0
        for species in subsets(motif.species):
            for reactions in subsets(motif.reactions):
                sub=Motif(species,reactions)
                if sub==motif or not self.side_incident(sub): continue
                checked+=1;decision=StrictLinearFeasibility.solve(self.balance_rows([sub])[:,reactions])
                if decision.status!='infeasible': return dict(is_pac=False if decision.status=='feasible' else None,reason='proper submotif '+decision.status)
        return dict(is_pac=True,proper_side_incident_submotifs_checked=checked)

    def record(self):
        return dict(species=self.species,orientation=self.orientation,reactions=[dict(reactants=r.reactants,products=r.products,factor=str(r.factor)) for r in self.reactions],
            motifs=[dict(species=c.species,reactions=c.reactions) for c in self.motifs],complexes=self.complexes)


@dataclass(frozen=True)
class TriangleInterface:
    gain: F
    shared: F
    first: F
    second: F

    def __post_init__(self):
        for name in ('gain','shared','first','second'): object.__setattr__(self,name,F(getattr(self,name)))
        if self.gain<=1 or min(self.shared,self.first,self.second)<=0: raise ValueError('Gain > 1 and positive factors required.')

    @property
    def resistance(self): return self.shared*(1/self.first+1/self.second)
    @property
    def interval(self): return self.resistance/self.gain,self.resistance

    def currents(self,q):
        q=F(q);lower,upper=self.interval
        if not lower<q<upper: raise ValueError('Response is outside the strict open interval.')
        x,y=self.shared/self.first,self.shared/self.second;t=q/(x+y);ell=1/self.gain
        e=(t-ell)*(1-t)/(2*(x+y));u=t+e*y;v=t-e*x
        if not ell<v<u<1 or x*u+y*v!=q: raise ArithmeticError('Interface reconstruction failed.')
        return F(1),u,v

    def overlap(self,other):
        lower=max(self.interval[0],other.interval[0]);upper=min(self.interval[1],other.interval[1])
        return (lower,upper) if lower<upper else None


class TrianglePair:
    def __init__(self,factors,gain=2):
        factors=tuple(F(v) for v in factors)
        if len(factors)!=5 or not isinstance(gain,int) or gain<=1: raise ValueError('Five factors and one integer gain > 1 required.')
        self.gain=gain;self.factors=factors;b0,b1,b2,b3,b4=factors
        self.left=TriangleInterface(gain,b0,b1,b2);self.right=TriangleInterface(gain,b0,b3,b4)
        A=(1,0,0,0);B=(0,1,0,0);C=(0,0,1,0);D=(0,0,0,1);mA=(gain,0,0,0)
        self.family=CoreFamily(('A','B','C','D'),[Reaction(a,b,k) for (a,b),k in zip(((A,B),(B,C),(C,mA),(B,D),(D,mA)),factors)],
            [Motif((0,1,2),(0,1,2)),Motif((0,1,3),(0,3,4))])

    def witness(self,a=F(1,2)):
        a=F(a)
        if not 0<a<1: raise ValueError('A activity must lie strictly between zero and one.')
        overlap=self.left.overlap(self.right)
        if overlap is None: return None
        q=sum(overlap)/2;_,ul,vl=self.left.currents(q);_,ur,vr=self.right.currents(q)
        b=(a**self.gain+q*a)/(1+q);j=self.factors[0]*(a-b)
        c=b-ul*j/self.factors[1];d=b-ur*j/self.factors[3]
        result=self.family.evaluate((a,b,c,d))
        if not result['strictly_productive']: raise ArithmeticError('Species reconstruction failed literal network replay.')
        result.update(q=str(q),overlap=list(map(str,overlap)))
        return result


def toric_family():
    return CoreFamily(('A','B'),[Reaction((1,0),(0,1),1),Reaction((0,1),(2,0),1),Reaction((0,1),(3,0),2)],
        [Motif((0,1),(0,1)),Motif((0,1),(0,2))])


def direction_family():
    species=('e1','e2','e3','e4','e1p','e2p','eA','eB')
    def c(**kwargs): return tuple(kwargs.get(s,0) for s in species)
    pairs=[(c(e1=1),c(eB=1,e2=1)),(c(e2=1,e2p=1),c(e3=1)),(c(e3=1),c(e4=2)),
           (c(eA=1,e4=1),c(e1=1)),(c(e1p=1),c(eA=1,e2p=1)),(c(eB=1,e4=1),c(e1p=1))]
    return CoreFamily(species,[Reaction(a,b,1) for a,b in pairs],[Motif((0,1,2,3),(0,1,2,3)),Motif((4,5,2,3),(4,1,2,5))])


def dominance_family():
    return [dict(m=m,n=n,bm=bm,bn=bn) for m in range(2,9) for n in range(m+1,9) for bm in range(1,13) for bn in range(1,13) if F(1,m*bm)>=F(1,bn)]


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(sweep,toric,output):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    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');k=np.array([float(F(r['k'])) for r in sweep])
        ax.axhspan(1,2,color='#0072B2',alpha=.12,label='Left interval: 1 < q < 2')
        ax.plot(k,1/k,'--',color='#D55E00',label='Right boundaries: 1/k and 2/k');ax.plot(k,2/k,'--',color='#D55E00')
        lower=np.maximum(1,1/k);upper=np.minimum(2,2/k)
        ax.fill_between(k,lower,upper,where=lower<upper,color='#009E73',alpha=.4,label='Common open response interval')
        ax.axvline(.5,color='#666666',linewidth=.8);ax.axvline(2,color='#666666',linewidth=.8)
        ax.set(xlabel='Private-factor ratio k',ylabel='Response q = (B - A squared) / (A - B)',ylim=(0,3.3),xlim=(min(k),max(k)))
        ax.legend(fontsize=8,loc='upper right');ax.spines[['top','right']].set_visible(False)
        fig.savefig(output/'phase.png',dpi=220);fig.savefig(output/'phase.svg');plt.close(fig)
        fig,ax=plt.subplots(figsize=(7.2,4.5),layout='constrained')
        x=np.array([float(F(r['A'])) for r in toric]);lower=np.array([float(F(r['B_lower'])) for r in toric]);upper=np.array([float(F(r['B_upper'])) for r in toric])
        ax.plot(x,lower,color='#0072B2',label='B must be above this: core 1 produces A')
        ax.plot(x,upper,'--',color='#D55E00',label='B must be below this: core 2 produces B')
        ax.fill_between(x,upper,lower,color='#888888',alpha=.18,label='Incompatible gap for 0 < A < 1')
        ax.set(xlabel='Species activity A',ylabel='Required bounds on species activity B',xlim=(0,1),ylim=(0,1))
        ax.legend(fontsize=8);ax.spines[['top','right']].set_visible(False)
        fig.savefig(output/'monomial.png',dpi=220);fig.savefig(output/'monomial.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(factors=FACTORS,gain=GAIN,A_activity=A_ACTIVITY,k_sweep=list(map(str,K_SWEEP)),temperature_K=TEMPERATURE_K,
        gas_constant_J_mol_K=GAS_CONSTANT_J_MOL_K,rational_denominator_limit=RATIONAL_DENOMINATOR_LIMIT,max_pac_enumeration_bits=MAX_PAC_ENUMERATION_BITS)
    intro='Resolved inputs: '+json.dumps(inputs);print(intro,flush=True)
    chosen=TrianglePair(FACTORS,GAIN);witness=chosen.witness(F(A_ACTIVITY));sweep=[]
    for k in K_SWEEP:
        pair=TrianglePair((1,1,1,k,k));w=pair.witness(F(A_ACTIVITY));overlap=pair.left.overlap(pair.right)
        sweep.append(dict(k=str(k),compatible=w is not None,overlap_lower=str(overlap[0]) if overlap else '',overlap_upper=str(overlap[1]) if overlap else '',
            minimum_motif_production=str(min(map(F,w['motif_production']))) if w else '',activities=' '.join(w['activities']) if w else ''))
    write_csv(args.output/'phase_sweep.csv',sweep)
    families={'direction_obstruction':direction_family(),'magnitude_obstruction':TrianglePair((1,1,1,10,10)).family,'monomial_obstruction':toric_family(),'selected_triangles':chosen.family}
    records={};layer_rows=[]
    for name,family in families.items():
        filters=family.filters();audits=[family.audit_pac(c) for c in family.motifs]
        if any(d.status=='unresolved' for d in filters.values()) or any(a['is_pac'] is not True for a in audits):
            raise ArithmeticError('Numerical candidate could not be certified by rational arithmetic.')
        records[name]=dict(network=family.record(),filters={name:d.record() for name,d in filters.items()},pac_audits=audits)
        layer_rows.append(dict(example=name,**{key:d.status for key,d in filters.items()},species_activity='infeasible by manuscript analytic obstruction' if name!='selected_triangles' else 'feasible exact witness' if witness else 'infeasible exact interval separation'))
    write_csv(args.output/'layers.csv',layer_rows)
    direction=direction_family()
    if any(v!=0 for v in direction.direction_rows().sum(axis=0)): raise ArithmeticError('All-ones signed circuit failed.')
    toric=toric_family();complex_witness={(1,0):F(4),(0,1):F(3),(2,0):F(9,4),(3,0):F(11,4)}
    toric_linear=toric.incidence@[complex_witness[c] for c in toric.complexes]
    if tuple(toric_linear)!=(F(1),F(3,4),F(1,2)) or not all(v>0 for v in toric.flow_rows()@toric_linear):
        raise ArithmeticError('Paper linear-complex witness failed.')
    toric_rows=[]
    for i in range(1,100):
        a=F(i,100);lower=(a+2*a*a)/3;upper=(a+2*a**3)/3
        if lower-upper!=2*a*a*(1-a)/3 or not lower>upper: raise ArithmeticError('Monomial bound identity failed.')
        toric_rows.append(dict(A=str(a),B_lower=str(lower),B_upper=str(upper),gap=str(lower-upper)))
    write_csv(args.output/'monomial_bounds.csv',toric_rows);dominance=dominance_family();write_csv(args.output/'dominance_family.csv',dominance)
    (args.output/'models_and_certificates.json').write_text(json.dumps(records,indent=2)+'\n',encoding='utf-8')
    plot(sweep,toric_rows,args.output)
    if TEMPERATURE_K<=0: raise ValueError('Temperature must be positive.')
    summary=dict(inputs=inputs,selected_witness=witness,layers=layer_rows,dominance_count=len(dominance),
        barrier_window_J_per_mol=GAS_CONSTANT_J_MOL_K*TEMPERATURE_K*log(2),
        monomial_contradiction='Directions force 0<A<1; required B lower minus required B upper = 2 A^2 (1-A) / 3 > 0.',
        evidence='Exact rational witnesses and linear alternative certificates. Analytical species-level obstruction is from the manuscript. Static per-core productivity is not a steady state, stability, or whole-network-growth guarantee. 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 scipy,matplotlib
    metadata=dict(paper_sha256=PAPER_SHA256,source_sha256=hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),python=platform.python_version(),
        numpy=np.__version__,scipy=scipy.__version__,matplotlib=matplotlib.__version__,platform=platform.platform(),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()