Example code

"""Planar mass-action cusps: reusable kinetics, corrected jets and exact catalogue checks."""
# EDITABLE INPUTS. Concentrations and time are nondimensional.
EPSILON = '1/10'  # rational bistable family: mu=0, nu=epsilon^2, 0<epsilon<1
CUSTOM_RATES = None  # or five positive rational strings for template 02
INITIAL_STATES = ((.85,.7),(1.15,1.4))
TIME_HORIZON = 30000.
TIME_SAMPLES = 601
SWEEP_NU = '1/100'
SWEEP_MU_MAX = '1/1000'
SWEEP_SAMPLES = 201
QUERY_REACTIONS = (0,5,12,18,22)  # zero-based IDs; complete table exported
RUN_EXHAUSTIVE_CENSUS = True
VERIFY_ALL_POSITIVE_TEMPLATES = True
MAX_LITERAL_SETS = 142506

from collections import Counter
from dataclasses import dataclass
from fractions import Fraction as Q
from itertools import combinations
import argparse
import csv
import hashlib
import json
import math
from pathlib import Path
import platform
import numpy as np
import sympy as sp
from scipy.integrate import solve_ivp

MANUSCRIPT_SHA256='9c35a17e2f132b80f00d057d1ee52f1ba89391f78fd752be2d120f43b0ae086c'
COMPLEXES=((0,0),(1,0),(0,1),(2,0),(1,1),(0,2))
NAMES=('0','X','Y','2X','X+Y','2Y')
REACTIONS=tuple((a,b) for a in COMPLEXES for b in COMPLEXES if a!=b)
INDEX={r:i for i,r in enumerate(REACTIONS)}
t=sp.Symbol('t');x,y=sp.symbols('x y')


class PlanarNetwork:
    def __init__(self,indices):
        self.indices=tuple(indices)
        if len(self.indices)!=5 or len(set(self.indices))!=5 or any(type(i) is not int or not 0<=i<30 for i in self.indices):raise ValueError('Five distinct reaction IDs from 0 through 29 required.')
        self.reactions=tuple(REACTIONS[i] for i in self.indices)
        self.sources=tuple(a for a,b in self.reactions);self.vectors=tuple((b[0]-a[0],b[1]-a[1]) for a,b in self.reactions)
        self.N=sp.Matrix(self.vectors).T;self.A=sp.Matrix(self.sources)
    def key(self):
        rows=[]
        for a,v in zip(self.sources,self.vectors):
            g=math.gcd(*v);rows.append((*a,v[0]//g,v[1]//g))
        return min(tuple(sorted(rows)),tuple(sorted((r[1],r[0],r[3],r[2]) for r in rows)))
    def eligible(self):
        if len(set(self.sources))<4:return False
        for a in self.vectors:
            crosses=[a[0]*b[1]-a[1]*b[0] for b in self.vectors]
            if min(crosses)>=0 or max(crosses)<=0:return False
        return True  # surrounding condition already implies rank two
    def swapped(self):return tuple(sorted(INDEX[((a[1],a[0]),(b[1],b[0]))] for a,b in self.reactions))
    def symbolic_field(self,rates):return self.N*sp.Matrix([k*x**a[0]*y**a[1] for k,a in zip(rates,self.sources)])
    def hessian_at_one(self,rates,u,v):
        return self.N*sp.Matrix([k*(a*(a-1)*u[0]*v[0]+a*b*(u[0]*v[1]+u[1]*v[0])+b*(b-1)*u[1]*v[1]) for k,(a,b) in zip(rates,self.sources)])
    def transport_from(self,source):
        """Return source-index matching, species exchange and rate factors.

        target rate = source rate * factor; reactant labels remain fixed.
        """
        for swap in (False,True):
            unused=set(range(5));mapping=[];factors=[]
            for a,v in zip(self.sources,self.vectors):
                match=None
                for j in sorted(unused):
                    aa=source.sources[j][::-1] if swap else source.sources[j];vv=source.vectors[j][::-1] if swap else source.vectors[j]
                    if a!=aa or v[0]*vv[1]!=v[1]*vv[0] or sum(c*d for c,d in zip(v,vv))<=0:continue
                    coordinate=next(i for i in range(2) if v[i]);match=(j,sp.Rational(vv[coordinate],v[coordinate]));break
                if match is None:break
                j,factor=match;unused.remove(j);mapping.append(j);factors.append(factor)
            if not unused:return dict(source_indices=mapping,rate_factors=factors,species_exchange=swap)
        raise ValueError('Networks are not equivalent by labelled positive rays and optional species exchange.')


class AlgebraicTemplate:
    def __init__(self,row):
        self.row=row;self.network=PlanarNetwork(tuple(row['indices']));self.modulus=sp.Poly(row['minimal'],t,domain=sp.QQ)
        self.lo,self.hi=map(sp.Rational,row['interval'])
    def red(self,value):
        n,d=sp.cancel(value).as_numer_denom();m=self.modulus.as_expr()
        return sp.rem(sp.rem(n,m,t)*sp.invert(sp.rem(d,m,t),m,t),m,t)
    def interval(self,expression):
        lo=hi=sp.Rational(0)
        for c in sp.Poly(expression,t,domain=sp.QQ).all_coeffs():
            candidates=[lo*self.lo,lo*self.hi,hi*self.lo,hi*self.hi];lo=min(candidates)+c;hi=max(candidates)+c
        return lo,hi
    def verify(self):
        row=self.row;m=self.modulus.as_expr()
        if self.lo==self.hi:
            if self.modulus.degree()!=1 or self.modulus.eval(self.lo)!=0:raise ArithmeticError('Invalid rational root witness.')
        elif self.modulus.count_roots(self.lo,self.hi)!=1 or self.modulus.eval(self.lo)*self.modulus.eval(self.hi)>=0:raise ArithmeticError('Root is not uniquely isolated.')
        parse=lambda key:sp.Matrix([sp.sympify(v,locals={'t':t}) for v in row[key]])
        rates=parse('rates');q=parse('right_kernel');p=parse('left_kernel').T;h=parse('center');N=self.network.N;A=self.network.A;J=N*sp.diag(*rates)*A
        B=lambda u,v:self.network.hessian_at_one(rates,u,v).applyfunc(self.red)
        identities=list(N*rates)+list(J*q)+list(p*J)+[(p*q)[0]-1,(p*h)[0],(p*B(q,q))[0]]+list(J*h+B(q,q))
        if any(self.red(v)!=0 for v in identities):raise ArithmeticError('Equilibrium or cusp jet identity failed.')
        tau=self.red(sp.trace(J));c=self.red((p*B(q,h))[0]/2)
        first=[self.red((p*N[:,r])[0]) for r in range(5)]
        raw=[self.red(first[r]*(A[r,0]*q[0]+A[r,1]*q[1])) for r in range(5)]
        corrected_scaled=[self.red(tau*raw[r]-(p*B(q,N[:,r]))[0]) for r in range(5)]
        i,j=row['pair'];D=self.red(first[i]*corrected_scaled[j]-first[j]*corrected_scaled[i]);raw_minor=self.red(first[i]*raw[j]-first[j]*raw[i])
        U=sp.sympify(row['U'],locals={'t':t});V=sp.sympify(row['V'],locals={'t':t})
        if sp.expand(U*D+V*m-1)!=0:raise ArithmeticError('Recomputed corrected minor failed its Bezout identity.')
        for label,value in [('tau',tau),('c',c),('D',D)]:
            if self.red(value-sp.sympify(row[label],locals={'t':t}))!=0:raise ArithmeticError('Stored certificate disagrees with literal derivatives.')
        enclosures={name:self.interval(value) for name,value in [('tau',tau),('c',c),('D',D)]}
        if enclosures['tau'][1]>=0 or enclosures['c'][1]>=0 or enclosures['D'][0]<=0 or any(self.interval(k)[0]<=0 for k in rates):raise ArithmeticError('Strict sign enclosure failed.')
        self.rates=rates
        return dict(ordinal=row['ordinal'],degree=self.modulus.degree(),tau=str(tau),c=str(c),trace_scaled_corrected_minor=str(D),raw_minor=str(raw_minor),corrected_minor=str(self.red(D/tau)),sign_enclosures={k:list(map(str,v)) for k,v in enclosures.items()},root_interval=row['interval'],rates=list(map(str,rates)))
    def numerical_reactor(self,target=None,bisections=80):
        """Instantiate positive rational approximations for ODE exploration.

        Exact cusp certification applies at the algebraic root, not its midpoint.
        """
        if not hasattr(self,'rates'):self.verify()
        lo,hi=self.lo,self.hi
        for _ in range(bisections):
            if lo==hi:break
            mid=(lo+hi)/2;value=self.modulus.eval(mid)
            if value==0:lo=hi=mid;break
            if self.modulus.eval(lo)*value<0:hi=mid
            else:lo=mid
        rates=[k.subs(t,(lo+hi)/2) for k in self.rates]
        if target is not None:
            match=target.transport_from(self.network);rates=[rates[j]*factor for j,factor in zip(match['source_indices'],match['rate_factors'])]
        return MassActionReactor(target or self.network,rates)


class Catalogue:
    def __init__(self,path=Path(__file__).with_name('catalogue.json')):
        self.data=json.loads(Path(path).read_text());self.records={};self.templates={r['ordinal']:r for r in self.data['positive_templates']}
        for label,networks in self.data['partition'].items():
            for indices in networks:
                network=PlanarNetwork(tuple(indices));key=network.key()
                if key in self.records:raise ArithmeticError('Duplicate mechanism key.')
                self.records[key]=(label,tuple(indices))
    def classify(self,network):
        if not network.eligible():return dict(admits_cusp=False,evidence='Necessary rank/positive-flux/four-source filter fails; paper theorem.',layer='structural')
        if network.key() not in self.records:raise ArithmeticError('Eligible key missing; catalogue coverage unresolved.')
        label,indices=self.records[network.key()];mapping=network.transport_from(PlanarNetwork(indices))
        return dict(admits_cusp=label=='cusp',layer=label,representative=indices,transport=mapping,evidence='Published exact classification lookup; negative obstruction payloads and Lean proof are not replayed here.')
    def census(self,limit=MAX_LITERAL_SETS):
        keys=set();orbits=set();literal=eligible=0;sources=Counter()
        for indices in combinations(range(30),5):
            literal+=1
            if literal>limit:raise ValueError('Census budget exceeded; partial counts are not exhaustive.')
            # Cheap filter avoids constructing symbolic matrices for most candidates.
            pairs=[REACTIONS[i] for i in indices]
            if len({a for a,b in pairs})<4:continue
            vs=[(b[0]-a[0],b[1]-a[1]) for a,b in pairs]
            if any(not (min(a[0]*b[1]-a[1]*b[0] for b in vs)<0<max(a[0]*b[1]-a[1]*b[0] for b in vs)) for a in vs):continue
            network=PlanarNetwork(indices);eligible+=1;keys.add(network.key());orbits.add(min(indices,network.swapped()))
        if keys!=set(self.records):raise ArithmeticError('Imported partition does not cover exactly the regenerated eligible keys.')
        for key,(label,indices) in self.records.items():sources[label,len({REACTIONS[i][0] for i in indices})]+=1
        return dict(literal_sets=literal,eligible_sets=eligible,species_orbits=len(orbits),mechanism_classes=len(keys),layers=dict(Counter(label for label,ids in self.records.values())),source_counts={f'{label}:{n}':v for (label,n),v in sources.items()},scope='Enumeration and partition coverage replayed; imported obstruction tags remain dependent on the paper theorem.')


class MassActionReactor:
    def __init__(self,network,rates):
        self.network=network;self.rates=tuple(sp.Rational(k) for k in rates)
        if len(self.rates)!=5 or min(self.rates)<=0:raise ValueError('Five independent positive rates required.')
        self.N=np.array(network.N,float);self.A=np.array(network.A,int);self.k=np.array(self.rates,float)
    def field(self,state):return self.N@(self.k*np.prod(np.asarray(state,float)**self.A,axis=1))
    def jacobian(self,state):
        state=np.asarray(state,float);v=self.k*np.prod(state**self.A,axis=1)
        return self.N@(v[:,None]*self.A/state)
    def integrate(self,initial,times,method='Radau'):
        initial=np.asarray(initial,float);times=np.asarray(times,float)
        if initial.shape!=(2,) or min(initial)<=0 or times[0]!=0 or np.any(np.diff(times)<=0):raise ValueError('Positive initial concentrations and increasing times from zero required.')
        def fun(t,z):
            state=np.exp(z);return self.field(state)/state
        def jac(t,z):
            state=np.exp(z);return self.jacobian(state)*state[None,:]/state[:,None]-np.diag(self.field(state)/state)
        solution=solve_ivp(fun,(0,times[-1]),np.log(initial),t_eval=times,jac=jac,method=method,rtol=2e-10,atol=2e-12)
        states=np.exp(solution.y.T)
        if not solution.success or np.any(states<=0) or not np.all(np.isfinite(states)):raise ArithmeticError('Positive integration failed.')
        return states


class RationalCuspReactor(MassActionReactor):
    def __init__(self,rates):super().__init__(PlanarNetwork((0,5,12,18,22)),rates)
    @classmethod
    def unfolding(cls,mu,nu):
        mu,nu=map(sp.Rational,(mu,nu));return cls([(1+mu-nu)/11,(3-nu)/11,sp.Rational(3,11),sp.Rational(3,11),sp.Rational(1,11)])
    @classmethod
    def bistable(cls,epsilon):
        epsilon=sp.Rational(epsilon)
        if not 0<epsilon<1:raise ValueError('The explicit bistable family requires 0 < epsilon < 1.')
        return cls.unfolding(0,epsilon**2)
    def equilibrium_polynomial(self):
        k1,k2,k3,k4,k5=self.rates
        return sp.Poly(k1-k2*x+k4*x*x-k5*k4/k3*x**3,x,domain=sp.QQ)
    def equilibria(self):
        polynomial=self.equilibrium_polynomial();k1,k2,k3,k4,k5=self.rates;rows=[]
        for (lo,hi),multiplicity in polynomial.intervals(eps=sp.Rational(1,10**14)):
            if hi<=0:continue
            if lo<=0:raise ArithmeticError('Positive root isolation unresolved.')
            xx=(lo+hi)/2;yy=k4/k3*xx**2;trace=-k2-2*k4*xx-k5*yy-k3;det=-k3*polynomial.diff().eval(xx)
            stability='nonhyperbolic' if multiplicity>1 else ('stable' if det>0 else 'saddle')
            # Exact interval derivative sign must agree; midpoint alone is insufficient.
            if multiplicity==1:
                derivative=sp.Poly(polynomial.diff(),x);low=high=sp.Rational(0)
                for coefficient in derivative.all_coeffs():
                    candidates=[low*lo,low*hi,high*lo,high*hi];low=min(candidates)+coefficient;high=max(candidates)+coefficient
                if low<=0<=high:raise ArithmeticError('Derivative sign needs a tighter root enclosure.')
            rows.append(dict(x_interval=[str(lo),str(hi)],x=float(xx),y=float(yy),multiplicity=multiplicity,stability=stability,trace_midpoint=str(trace),determinant_midpoint=str(det)))
        return rows


def rational_identities():
    mu,nu,z,epsilon=sp.symbols('mu nu z epsilon');a,b=sp.symbols('a b')
    network=PlanarNetwork((0,5,12,18,22));F=network.symbolic_field([a,b,sp.Rational(3,11),sp.Rational(3,11),sp.Rational(1,11)])
    scalar=sp.expand(11*F[0].subs(y,x*x).subs({x:1+z,a:(1+mu-nu)/11,b:(3-nu)/11}))
    if sp.expand(scalar-(mu+nu*z-z**3))!=0:raise ArithmeticError('Exact equilibrium reduction failed.')
    field=F.subs({a:(1-epsilon**2)/11,b:(3-epsilon**2)/11});J=field.jacobian([x,y]);rows=[]
    for root in (1-epsilon,sp.Integer(1),1+epsilon):
        at={x:root,y:root**2}
        if any(sp.simplify(v) for v in field.subs(at)):raise ArithmeticError('Explicit equilibrium failed.')
        rows.append(dict(x=str(root),y=str(root**2),determinant=str(sp.factor(J.subs(at).det())),trace=str(sp.expand(sp.trace(J.subs(at))))))
    aa,bb,dd,ee,ff=sp.symbols('a b d e f');negative=PlanarNetwork((0,5,12,20,24));G=negative.symbolic_field([aa,bb,dd,ee,ff]);det=sp.factor(G.jacobian([x,y]).det().subs(x,dd/(ff-ee)))
    if sp.simplify(det+dd*(ff-3*ee)*y)!=0 or any(sp.simplify(v) for v in G.subs({ff:3*ee,x:dd/(2*ee),aa:bb*dd/(2*ee)})):raise ArithmeticError('Same-vector negative example failed.')
    return dict(scalar=str(scalar),equilibria=rows,negative_collision_determinant=str(det),scope='Scalar identity describes equilibria on the y-nullcline, not a one-dimensional time evolution or invariant manifold.')


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
    catalogue=Catalogue();census=catalogue.census() if RUN_EXHAUSTIVE_CENSUS else {'scope':'Census skipped by input.'}
    checked=[]
    for ordinal,row in catalogue.templates.items():
        if VERIFY_ALL_POSITIVE_TEMPLATES or ordinal==2:checked.append(AlgebraicTemplate(row).verify())
    model=RationalCuspReactor(CUSTOM_RATES) if CUSTOM_RATES is not None else RationalCuspReactor.bistable(EPSILON);equilibria=model.equilibria();identities=rational_identities();query=catalogue.classify(PlanarNetwork(tuple(QUERY_REACTIONS)))
    times=np.unique(np.r_[0,np.geomspace(1e-3,min(100.,TIME_HORIZON),120),np.linspace(0,TIME_HORIZON,TIME_SAMPLES)])
    trajectories=[model.integrate(initial,times) for initial in INITIAL_STATES]
    differences=[float(np.max(abs(a-model.integrate(initial,times,'BDF')))) for initial,a in zip(INITIAL_STATES,trajectories)]
    sweep=[];maximum=sp.Rational(SWEEP_MU_MAX)
    for i in range(SWEEP_SAMPLES):
        mu=-maximum+2*maximum*i/(SWEEP_SAMPLES-1);sample=RationalCuspReactor.unfolding(mu,SWEEP_NU)
        # Exact discriminant identifies root-count regime; plotted roots are numerical.
        disc=sp.discriminant(sample.equilibrium_polynomial().as_expr(),x)
        for root in np.roots([float(c) for c in sample.equilibrium_polynomial().all_coeffs()]):
            if abs(root.imag)<1e-9:
                xx=float(root.real);yy=xx*xx;det=np.linalg.det(sample.jacobian((xx,yy)))
                sweep.append((float(mu),float(sp.Rational(SWEEP_NU)),xx,yy,'stable' if det>0 else 'saddle',int(sp.sign(disc))))
    def table(name,header,rows):
        with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    def write_json(name,value):(out/name).write_text(json.dumps(value,indent=2,default=str)+'\n')
    table('reactions.csv',['id','reactant','product','dx','dy'],[(i,NAMES[COMPLEXES.index(a)],NAMES[COMPLEXES.index(b)],b[0]-a[0],b[1]-a[1]) for i,(a,b) in enumerate(REACTIONS)])
    table('trajectories.csv',['initial_index','time','X','Y'],[(i,tt,*state) for i,values in enumerate(trajectories) for tt,state in zip(times,values)])
    table('unfolding_sweep.csv',['mu','nu','X','Y','linear_stability','exact_discriminant_sign'],sweep)
    table('positive_templates.csv',['ordinal','reaction_ids','degree','tau','c','trace_scaled_corrected_minor','raw_minor','corrected_minor'],[(r['ordinal'],','.join(map(str,catalogue.templates[r['ordinal']]['indices'])),*[r[k] for k in ('degree','tau','c','trace_scaled_corrected_minor','raw_minor','corrected_minor')]) for r in checked])
    result=dict(census=census,positive_templates_checked=checked,query=query,rational_identities=identities,reactor_rates=list(map(str,model.rates)),equilibria=equilibria,numerical_solver_differences=differences,
        scope='Exact enumeration, all selected positive jet/sign/corrected-minor checks, and rational bistable model. Negative-layer classification uses the imported paper theorem; its 9,947 obstruction payloads and Lean/native evaluation are not rerun. No global basin or arbitrary-supernetwork claim.')
    write_json('results.json',result)
    lines=[f'Census: {census}.',f'Positive templates freshly verified: {len(checked)}; root degree counts: {dict(Counter(r["degree"] for r in checked))}.',
        f'Template 02: tau={next(r for r in checked if r["ordinal"]==2)["tau"]}, c={next(r for r in checked if r["ordinal"]==2)["c"]}, corrected minor={next(r for r in checked if r["ordinal"]==2)["corrected_minor"]}.',
        f'Configured reactor equilibria: {[(r["x"],r["y"],r["stability"]) for r in equilibria]}.',f'Independent solver differences: {differences}.',
        'The nullcline cubic is an equilibrium reduction, not an invariant one-dimensional dynamical model. Imported negative tags are not fresh obstruction proofs.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained');nu=np.linspace(0,.04,301);fold=2*(nu/3)**1.5
    axs[0].fill_between(nu,-fold,fold,color='#417d8c',alpha=.2);axs[0].plot(nu,fold,color='#417d8c');axs[0].plot(nu,-fold,color='#417d8c');axs[0].set(xlabel='nu',ylabel='mu',title='Three-equilibrium region between folds')
    for label,color in [('stable','#417d8c'),('saddle','#bd5a24')]:
        points=np.array([(r[0],r[2]) for r in sweep if r[4]==label]);axs[1].scatter(points[:,0],points[:,1],s=5,color=color,label=label) if len(points) else None
    axs[1].set(xlabel='mu',ylabel='Equilibrium X',title=f'Equilibrium branches at nu={SWEEP_NU}');axs[1].legend()
    axs[1].ticklabel_format(axis='x',style='sci',scilimits=(0,0))
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'unfolding.png',dpi=180);fig.savefig(out/'unfolding.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
    for i,values in enumerate(trajectories):
        axs[0].plot(values[:,0],values[:,1],label=f'Initial state {i+1}');axs[1].plot(times,values[:,0],label=f'X, initial {i+1}')
    for row in equilibria:axs[0].scatter(row['x'],row['y'],color='#bd5a24' if row['stability']=='saddle' else '#417d8c',marker='x' if row['stability']=='saddle' else 'o',zorder=5)
    axs[0].set(xlabel='X',ylabel='Y',title='Full two-dimensional mass-action trajectories');axs[1].set(xlabel='Time',ylabel='X concentration',title='Slow relaxation near the cusp')
    for ax in axs:ax.legend(fontsize=8);ax.grid(alpha=.2)
    fig.savefig(out/'bistability.png',dpi=180);fig.savefig(out/'bistability.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    write_json('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,source_sha256=digest(Path(__file__)),catalogue_sha256=digest(Path(__file__).with_name('catalogue.json')),python=platform.python_version(),output_sha256={p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}))


if __name__=='__main__':main()