Example code

"""Exact Hopf/deletion/child certificates and a reusable high-order mass-action model."""
# EDITABLE STUDY INPUTS. This is an abstract high-order mathematical witness.
ORBIT_AMPLITUDE = 0.02
MIDPOINT_REFINEMENTS = (128, 256, 512)
REFERENCE_PARAMETER = '1'
SPECTRAL_RANGE = (0.5, 1.0)
SPECTRAL_SAMPLES = 201
CONTINUOUS_SAMPLES = 1201
SHOOTING_TOLERANCE = 2e-9

from dataclasses import dataclass
from fractions import Fraction as Q
from itertools import combinations, permutations, product
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
from scipy.optimize import root
from interval_certificate import certify

MANUSCRIPT_SHA256='9bdeae52d9666dd64bb1f20cacb59053d81494d3a2102b340139b8f6d9a3bb95'


class MassActionSource:
    def __init__(self,reactants,products):
        self.Y=sp.Matrix(reactants);self.P=sp.Matrix(products);self.S=self.P-self.Y
        if self.Y.shape!=self.P.shape or any(not x.is_Integer or x<0 for M in (self.Y,self.P) for x in M):raise ValueError('Matching nonnegative integer complexes required.')
        self.flux=sp.Matrix([2,3,2,2,2])
        if self.S.shape!=(4,5) or self.S*self.flux!=sp.zeros(4,1):raise ValueError('This example uses a four-species, five-reaction positive circuit.')
    @classmethod
    def paper(cls):return cls([[1,0,0,0,0],[0,4,0,0,0],[0,2,375,0,0],[400,0,0,2,0]],[[0,0,0,1,0],[0,0,0,6,0],[4,0,371,1,2],[398,0,2,0,2]])
    def skeleton(self):
        Y=self.S.applyfunc(lambda x:max(0,-x));return MassActionSource(Y,self.S+Y)
    def jet(self,order,inverse_state):
        tensor={}
        for indices in product(range(4),repeat=order):
            coefficients=[]
            for r in range(5):
                factor=1
                for i in range(4):
                    for k in range(indices.count(i)):factor*=int(self.Y[i,r])-k
                coefficients.append(factor)
            tensor[indices]=[inverse_state[i]*sum(int(self.S[i,r])*int(self.flux[r])*coefficients[r] for r in range(5)) for i in range(4)]
        return tensor
    def structural_certificate(self):
        minors=[self.S[:,[r for r in range(5) if r!=j]].det() for j in range(5)]
        if not all(v!=0 for v in minors) or self.S.rank()!=4:raise ArithmeticError('Positive circuit rank failed.')
        covectors=[]
        for j in range(5):
            retained=[r for r in range(5) if r!=j];w=self.S[:,retained].T.inv()*sp.ones(4,1)
            if w.T*self.S[:,retained]!=sp.ones(1,4):raise ArithmeticError('Deletion covector failed.')
            covectors.append(list(map(str,w)))
        return {'rank':4,'positive_kernel_vector':list(map(str,self.flux)),'four_column_minors':list(map(str,minors)),'deletion_covectors':covectors,
            'scope':'With at least one retained positive rate, each deleted system has a strictly increasing covector observable on the positive orthant. All-zero rates instead give constant solutions.'}


class ChildCertificates:
    def __init__(self,source):self.source=source
    def evaluate(self):
        S,Y=self.source.S,self.source.Y;children=[]
        for n in range(1,5):
            for rows in combinations(range(4),n):
                for cols in permutations(range(5),n):
                    if all(Y[i,j]>0 for i,j in zip(rows,cols)):children.append(frozenset(zip(rows,cols)))
        cover=[((2,3),(1,0),(1,2)),((0,2,3),(0,1,3),(10,2,7)),((0,1,2,3),(0,1,2,3),(100,35,40,140)),((1,3),(1,0),(1,2)),((2,3),(2,0),(1,2)),((1,2,3),(1,2,0),None)]
        sets=[frozenset(zip(rows,cols)) for rows,cols,_ in cover]
        if len(children)!=24 or not all(any(ch<=cs for cs in sets) for ch in children):raise ArithmeticError('Child cover failed.')
        records=[]
        for rows,cols,weights in cover:
            M=S.extract(rows,cols)
            if weights:
                Pi=sp.diag(*weights);H=-(Pi*M+M.T*Pi)
                minors=[H.extract(I,I).det() for size in range(1,len(rows)+1) for I in combinations(range(len(rows)),size)]
                if any(v<0 for v in minors):raise ArithmeticError('Child energy matrix is not positive semidefinite.')
                evidence={'weights':weights,'energy_matrix':[[str(v) for v in row] for row in H.tolist()],'all_principal_minors':list(map(str,minors))}
            else:
                lam=sp.Symbol('lambda');d=sp.symbols('d1:4');polynomial=(lam*sp.eye(3)-M*sp.diag(*d)).det()
                if sp.expand(polynomial-lam*(lam+4*d[0])*(lam+4*d[1]+2*d[2]))!=0:raise ArithmeticError('Singular child factorization failed.')
                evidence={'scaled_characteristic_factorization':'lambda*(lambda+4*d1)*(lambda+4*d2+2*d3)'}
            records.append({'species_zero_based':rows,'reactions_zero_based':cols,**evidence})
        return {'child_count':24,'covers':records,'conclusion':'Every positive diagonal scaling of every child has nonpositive spectral real parts; singular children are not claimed strictly stable.'}


class HopfCrossing:
    @staticmethod
    def h(t):return 2825760+t*(242612196+t*(3570069381-4001047375*t))
    @staticmethod
    def coefficients(t):return Q(386,125)+201*t,(562+5297*t)/250,(40+479*t)/250,6*t/25
    def isolate(self,steps=110):
        lo,hi=Q(1,2),Q(1)
        for _ in range(steps):
            mid=(lo+hi)/2
            if self.h(mid)>0:lo=mid
            else:hi=mid
        if not self.h(lo)>0>self.h(hi):raise ArithmeticError('Crossing root bracket failed.')
        return lo,hi
    def exact_certificate(self,source):
        t=sp.Symbol('t',positive=True);lam=sp.Symbol('lambda');inverse=[sp.Rational(1,2),sp.Rational(1,600),sp.Rational(1,1500),t/8]
        A=sp.diag(*inverse)*source.S*sp.diag(*source.flux)*source.Y.T
        coefficients=[sp.sympify(x) for x in self.coefficients(t)];a1,a2,a3,a4=coefficients
        polynomial=(lam*sp.eye(4)-A).det();expected=lam**4+a1*lam**3+a2*lam**2+a3*lam+a4
        delta=a1*a2*a3-a3*a3-a1*a1*a4;h=self.h(t)
        identities=[sp.factor(polynomial-expected),sp.factor(delta-sp.Rational(3,7812500)*h),
            sp.factor(expected-(lam*lam+a3/a1)*(lam*lam+a1*lam+a1*a4/a3)-lam*lam*delta/(a1*a3)),
            sp.expand(t*sp.diff(h,t)-3*h+8477280+485224392*t+3570069381*t*t)]
        alpha,beta,b,c=sp.symbols('alpha beta b c',real=True)
        generic=sp.Poly((lam*lam-2*alpha*lam+alpha*alpha+beta*beta)*(lam*lam+b*lam+c),lam).all_coeffs()
        _,g1,g2,g3,g4=generic;generic_delta=g1*g2*g3-g3*g3-g1*g1*g4
        identities += [sp.expand(generic_delta.subs(alpha,0)),sp.expand(sp.diff(generic_delta,alpha).subs(alpha,0)+2*b*((c-beta*beta)**2+b*b*beta*beta))]
        if any(x!=0 for x in identities):raise ArithmeticError('Literal vector field does not match the certified Hopf family.')
        if source.Y[:,:4].det()!=3000:raise ArithmeticError('Equilibrium determinant changed.')
        lo,hi=self.isolate();return {'identities_zero':True,'h_half':str(self.h(Q(1,2))),'h_one':str(self.h(Q(1))),
            'positive_root_bracket':[str(lo),str(hi)],'descartes_signs_ascending':[1,1,1,-1],
            'transversality':'t*h_prime=3*h-(8477280+485224392*t+3570069381*t^2), so h_prime<0 at the positive root; alpha_prime>0 by the exact Hurwitz derivative identity.',
            'equilibrium_det_numerator':'216000*k5^4','scope':'Exactly one positive crossing; stable before, two unstable roots after. Determinant of a positive equilibrium never vanishes.'}
    def numeric(self):
        lo,hi=self.isolate();t=float((lo+hi)/2);a1,a2,a3,a4=map(float,self.coefficients(t));omega=math.sqrt(a3/a1)
        hp=242612196+2*3570069381*t-3*4001047375*t*t;c=a1*a4/a3
        alpha_prime=-(3/7812500)*hp/(2*a1*((c-omega*omega)**2+a1*a1*omega*omega))
        return t,omega,alpha_prime


@dataclass(frozen=True)
class LogRates:
    values:tuple
    def __post_init__(self):
        if len(self.values)!=5 or not all(np.isfinite(self.values)):raise ValueError('Five finite log-rate constants required.')
    def equilibrium_log(self):
        k1,k2,k3,k4,k5=self.values
        x4=(k5-k4)/2;x3=(k5-k3)/375;x2=(math.log(1.5)+k5-k2-2*x3)/4;x1=k5-k1-400*x4
        return np.array([x1,x2,x3,x4])


class RelativeReactor:
    def __init__(self,source,state,flux):
        self.source=source;self.state=np.asarray(state,float);self.flux=np.asarray(flux,float)
        if self.state.shape!=(4,) or self.flux.shape!=(5,) or np.any(self.state<=0) or np.any(self.flux<=0) or not np.all(np.isfinite(self.state)) or not np.all(np.isfinite(self.flux)):raise ValueError('Positive finite stationary state and flux required.')
        self.S=np.array(source.S,float);self.Y=np.array(source.Y,float)
        if np.max(abs([email protected]))>1e-12*max(self.flux):raise ValueError('Stationary flux balance failed.')
        self.W=self.S*self.flux/self.state[:,None]
        self.log_rates=LogRates(tuple(np.log(self.flux)[email protected](self.state)))
    @classmethod
    def family(cls,source,t):
        if not np.isfinite(t) or t<=0:raise ValueError('Positive family parameter required.')
        return cls(source,[2,600,1500,8/t],[2,3,2,2,2])
    @classmethod
    def from_log_rates(cls,rates):
        state=np.exp(rates.equilibrium_log());k5=math.exp(rates.values[4]);return cls(MassActionSource.paper(),state,k5*np.array([1,1.5,1,1,1]))
    @property
    def jacobian(self):return [email protected]
    def field_jac(self,u,amplitude):
        if amplitude==0:return [email protected](u),self.jacobian
        eta=amplitude*np.asarray(u);z=1+eta
        if np.any(z<=0):raise ArithmeticError('Positive-domain exit; no clipping.')
        [email protected](eta)
        if np.max(logs)>700:raise ArithmeticError('Monomial overflow; reduce excursion or change numerical representation.')
        # Subtract the exactly balanced stationary flux before floating evaluation.
        return [email protected](logs)/amplitude,(self.W*np.exp(logs))@(self.Y.T/z)
    def integrate(self,initial,period,amplitude,method='Radau',samples=CONTINUOUS_SAMPLES):
        if amplitude<=0 or period<=0:raise ValueError('Positive amplitude and period required.')
        sol=solve_ivp(lambda s,u:self.field_jac(u,amplitude)[0],(0,period),initial,method=method,
            jac=lambda s,u:self.field_jac(u,amplitude)[1],rtol=2e-11,atol=2e-13,t_eval=np.linspace(0,period,samples))
        if not sol.success or np.any(1+amplitude*sol.y<=0):raise ArithmeticError('Continuous integration failed or exited the positive domain.')
        return sol.t,sol.y.T


class PeriodicShooter:
    def __init__(self,source,amplitude):
        if not 0<amplitude<=.05:raise ValueError('Local shooting amplitude must be in (0,0.05].')
        self.source=source;self.amplitude=amplitude;self.tH,self.omega,_=HopfCrossing().numeric()
        A=RelativeReactor.family(source,self.tH).jacobian;values,vectors=np.linalg.eig(A);j=next(i for i,z in enumerate(values) if z.imag>0)
        stable=sorted((i for i,z in enumerate(values) if abs(z.imag)<1e-9),key=lambda i:values[i].real,reverse=True)
        self.basis=np.column_stack([vectors[:,j].real,vectors[:,j].imag,vectors[:,stable[0]].real,vectors[:,stable[1]].real]);self.inverse=np.linalg.inv(self.basis)
    def initial(self,variables):return [email protected]([1,0,variables[0],variables[1]])
    def midpoint(self,variables,steps):
        reactor=RelativeReactor.family(self.source,variables[2]);dt=variables[3]/steps;u=self.initial(variables);trace=[u.copy()];midpoints=[]
        if dt<=0:raise ValueError('Positive trial period required.')
        for _ in range(steps):
            mid=u.copy()
            for _ in range(15):
                field,J=reactor.field_jac(mid,self.amplitude);correction=np.linalg.solve(np.eye(4)-dt*J/2,mid-u-dt*field/2);mid-=correction
                if np.linalg.norm(correction,np.inf)<2e-13:break
            else:raise ArithmeticError('Midpoint Newton solve failed.')
            u=2*mid-u;trace.append(u.copy());midpoints.append(mid.copy())
        return self.inverse@(u-trace[0]),np.array(trace),np.array(midpoints)
    def solve_midpoint(self,steps):
        if type(steps) is not int or not 32<=steps<=4096:raise ValueError('Midpoint steps outside budget.')
        variables=np.array([0.,0.,self.tH,2*steps*math.tan(math.pi/steps)/self.omega]);eps=[1e-5,1e-5,1e-5,1e-4]
        for _ in range(10):
            residual,trace,mids=self.midpoint(variables,steps)
            if np.linalg.norm(residual,np.inf)<SHOOTING_TOLERANCE:return variables,trace,mids,residual
            J=np.column_stack([(self.midpoint(variables+np.eye(4)[i]*e,steps)[0]-self.midpoint(variables-np.eye(4)[i]*e,steps)[0])/(2*e) for i,e in enumerate(eps)])
            variables-=np.linalg.solve(J,residual)
        raise ArithmeticError('Discrete shooting did not converge.')
    def continuous(self,guess):
        def defect(variables):
            reactor=RelativeReactor.family(self.source,variables[2]);_,trace=reactor.integrate(self.initial(variables),variables[3],self.amplitude,samples=2)
            return self.inverse@(trace[-1]-trace[0])
        answer=root(defect,guess,tol=1e-8)
        residual=defect(answer.x)
        if np.max(abs(residual))>2e-8:raise ArithmeticError('Continuous shooting residual too large: '+str(residual))
        reactor=RelativeReactor.family(self.source,answer.x[2]);times,trace=reactor.integrate(self.initial(answer.x),answer.x[3],self.amplitude)
        _,independent=reactor.integrate(self.initial(answer.x),answer.x[3],self.amplitude,method='BDF')
        return answer.x,times,trace,residual,float(np.max(abs(self.amplitude*(trace-independent))))


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)
    source=MassActionSource.paper();crossing=HopfCrossing();exact=crossing.exact_certificate(source);structure=source.structural_certificate();children=ChildCertificates(source).evaluate();interval=certify(source,crossing)
    shooter=PeriodicShooter(source,ORBIT_AMPLITUDE);refinements=[]
    for steps in MIDPOINT_REFINEMENTS:
        variables,trace,mids,residual=shooter.solve_midpoint(steps);reactor=RelativeReactor.family(source,variables[2]);ratios=np.exp(np.log1p(ORBIT_AMPLITUDE*mids)@reactor.Y)
        refinements.append({'steps':steps,'parameter':variables[2],'period':variables[3],'return_residual':float(np.max(abs(residual))),'mean_flux_ratios_midpoint':np.mean(ratios,axis=0).tolist(),'mean_z4_midpoint':float(np.mean(1+ORBIT_AMPLITUDE*mids[:,3]))})
    variables,times,trace,residual,solver_difference=shooter.continuous(variables);reactor=RelativeReactor.family(source,variables[2]);eta=ORBIT_AMPLITUDE*trace;z=1+eta;ratios=np.exp(np.log1p(eta)@reactor.Y)
    means=np.trapezoid(ratios,times,axis=0)/variables[3];mean_eta4=np.trapezoid(eta[:,3],times)/variables[3];variance=np.trapezoid((eta[:,3]-mean_eta4)**2,times)/variables[3]
    tH,omega,alpha_prime=crossing.numeric();spectral=[]
    for t in np.linspace(*SPECTRAL_RANGE,SPECTRAL_SAMPLES):
        ev=np.linalg.eigvals(RelativeReactor.family(source,t).jacobian);skeleton=np.linalg.eigvals(RelativeReactor.family(source.skeleton(),t).jacobian)
        spectral.append((t,crossing.h(t),max(ev.real),max(skeleton.real)))
    result={'exact_crossing':exact,'structure':structure,'children':children,'crossing_orientation':{'tH':tH,'omegaH':omega,'periodH':2*math.pi/omega,'alpha_prime':alpha_prime},'midpoint_refinements':refinements,
        'continuous_shooting':{'amplitude':ORBIT_AMPLITUDE,'parameter':variables[2],'period':variables[3],'stable_coordinates':variables[:2].tolist(),'return_residual_in_basis':float(np.max(abs(residual))),
            'independent_BDF_relative_difference':solver_difference,'minimum_relative_concentration':float(np.min(z)),'mean_flux_ratios_trapezoidal':means.tolist(),
            'mean_z4_minus_one':float(mean_eta4),'variance_z4':float(variance),'mean_square_identity_defect':float(2*mean_eta4+mean_eta4**2+variance),
            'scope':'Numerical continuous shooting and quadrature, not a validated orbit at this decimal parameter.'},
        'scope':'Exact child/deletion/crossing identities and a rational interval Lyapunov sign; local Hopf attraction uses the conventional theorem. No global attraction, explicit certified parameter neighborhood, or laboratory realization. Lean not rerun.'}
    def write_json(name,data):(out/name).write_text(json.dumps(data,indent=2)+'\n')
    write_json('results.json',result);write_json('lyapunov_certificate.json',interval)
    def table(name,header,rows):
        with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    table('orbit.csv',['time',*[f'z{i+1}' for i in range(4)],*[f'flux_ratio_{i+1}' for i in range(5)]],np.column_stack([times,z,ratios]))
    table('spectral_crossing.csv',['parameter','h','full_spectral_abscissa','skeleton_spectral_abscissa'],spectral)
    parameter=sp.Rational(REFERENCE_PARAMETER)
    if parameter<=0:raise ValueError('Positive reference parameter required.')
    x=sp.Matrix([2,600,1500,8/parameter]);rates=[source.flux[j]/sp.prod(x[i]**source.Y[i,j] for i in range(4)) for j in range(5)]
    table('reference_rates.csv',['reaction','reference_parameter','exact_rate','log_rate_at_numerical_orbit'],[(j+1,str(parameter),str(rates[j]),reactor.log_rates.values[j]) for j in range(5)])
    table('deletion_covectors.csv',['deleted_reaction','w1','w2','w3','w4'],[(j+1,*w) for j,w in enumerate(structure['deletion_covectors'])])
    lines=[f'24 children certified D-nonunstable; all five deletion covectors checked exactly.',f'Unique crossing tH={tH:.15f}; omega={omega:.15f}; alpha_prime={alpha_prime:.9g}.',
        'Outward-rounded rational certificate: -23/1000 < unit-norm l1 < -22/1000.',f'Numerical continuous shooting: t={variables[2]:.15f}, period={variables[3]:.10f}, residual={np.max(abs(residual)):.3g}.',
        f'Mean z4 minus one={mean_eta4:.8g}; variance={variance:.8g}; independent solver difference={solver_difference:.3g}.',
        'Maximum reactant molecularity 401. Numerical orbit is not an interval-validated decimal witness; all-zero deleted rates give constant solutions.']
    (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');spectral=np.array(spectral)
    axs[0].plot(spectral[:,0],spectral[:,1]/1e9);axs[0].set(xlabel='Family parameter t',ylabel='h(t) / 1e9',title='Exact cubic determines the crossing')
    axs[1].plot(spectral[:,0],1e3*spectral[:,2],label='Padded source');axs[1].plot(spectral[:,0],1e3*spectral[:,3],label='Unpadded skeleton')
    axs[1].set(xlabel='Family parameter t',ylabel='1000 x spectral abscissa',title='Equilibrium stability with and without\ncatalytic padding');axs[1].legend(fontsize=8)
    for ax in axs:ax.axvline(tH,color='#bd5a24',ls='--');ax.axhline(0,color='gray',lw=.8);ax.grid(alpha=.2)
    fig.savefig(out/'crossing.png',dpi=180);fig.savefig(out/'crossing.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained');axs[0].plot(eta[:,0],1e6*eta[:,3]);axs[0].set(xlabel='Relative departure z1 - 1',ylabel='1e6 x (z4 - 1)',title='Numerically returned positive orbit')
    axs[0].set_xticks(np.linspace(-ORBIT_AMPLITUDE,ORBIT_AMPLITUDE,5))
    axs[1].plot(times/variables[3],eta[:,0],label='z1 - 1');axs[1].plot(times/variables[3],100*eta[:,2],label='100 x (z3 - 1)');axs[1].set(xlabel='Time / returned period',ylabel='Scaled relative departure',title='Concentration variations over one return');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'orbit.png',dpi=180);fig.savefig(out/'orbit.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest();base=Path(__file__).resolve().parent
    write_json('run_metadata.json',{'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':digest(Path(__file__)),'support_source_sha256':{'interval_certificate.py':digest(base/'interval_certificate.py')},'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()