Example code

A reactor can have only one positive steady state, yet small concentration disturbances can grow instead of returning to it. This example builds the paper's full seven-species reversible reactor, including every reverse reaction and linear degradation channel, and verifies exact eigenvalues 1±8i1\pm8i: their positive real part describes growing disturbances, and their imaginary part describes oscillation.

The model separates reaction structure, stationary fluxes, rate reconstruction, dynamics and exact mode verification into reusable components. Editable study inputs appear at the top of the Python file. The adjacent JSON preserves the paper's integer witness; a positive stationary-coordinate parameterization lets readers construct and sweep new reactors without repeatedly solving for their equilibria.

Four species oscillate with growing departures from the stationary state. The nonlinear and linear trajectories nearly overlap; halving the perturbation reduces the linearization error by approximately four.
Solid curves show the full nonlinear reactor; dashed curves show the exact mode's linear prediction. The error plot adds 1e-16 for logarithmic display, not as an error certificate.
The full physical Jacobian has an unstable complex pair while its static contraction is stable. The internal response develops a substantial phase lag as frequency increases.
Zero-frequency contraction loses the internal species' relaxation. The phase plot evaluates s on the imaginary axis; the exact growing-mode identity is checked separately at s=1+8i.

The exact certificate checks two rational matrix identities, rather than trusting rounded eigenvalues. All 126 proper nonempty species restrictions also fail the source's autocatalysis condition, verifying its stated notion of minimality. The source family's uniqueness theorem comes from the paper; the code does not infer uniqueness from a numerical root search.

A perturbation of relative amplitude 10−510^{-5} follows the growing linear mode over the displayed interval. Its maximum linearization error is approximately 7.34×10−87.34\times10^{-8}; halving the amplitude reduces this to 1.84×10−81.84\times10^{-8}. Independent Radau and BDF integrations agree within approximately 1.22×10−111.22\times10^{-11}. These are numerical diagnostics, separate from the exact instability certificate.

The internal species CC has a finite relaxation time. Accounting for its response to a mode of growth rate and frequency ss introduces (s−L)−1(s-L)^{-1}; assuming instantaneous equilibration instead uses its value at s=0s=0 produces a stable contracted matrix and misses the growing mode. The package checks the full response at the actual eigenvalue s=1+8is=1+8i, and plots its phase along the imaginary axis to show what static contraction discards.

An important distinction is preserved: the matrix H=DF(x)diag⁡(x)H=DF(x)\operatorname{diag}(x), obtained by multiplying each Jacobian column by its stationary concentration, is stable in this example, but it is not the Jacobian of the log-coordinate dynamics. That Jacobian is similar to the physical Jacobian and retains the unstable spectrum. Changing coordinates cannot remove the instability.

An optional exact linear solver constructs new stationary rate assignments for a prescribed relative mode. Every successful candidate is checked independently with rational arithmetic; timeouts remain unresolved and infeasibility concerns only that prescribed mode. Seven scientific test groups and a reuse guide accompany the package. The witness is mathematical, not experimentally calibrated, and oscillatory instability alone does not certify a periodic orbit or Hopf bifurcation.

Python source

"""An exact growing mode in a literal seven-species reversible reactor.

Exact stationary coordinates, full nonlinear dynamics, and dynamic elimination.
"""
# EDITABLE STUDY INPUTS. The adjacent JSON holds the paper's exact integer table.
INITIAL_RELATIVE_AMPLITUDE = 1e-5
TIME_HORIZON = 3.
TIME_SAMPLES = 1201
DESIGN_TARGET_REAL_PART = '1'
DESIGN_TARGET_IMAGINARY_PART = '8'
MODAL_NORMALIZATION = 200
DESIGN_TIMEOUT_MS = 5000
FREQUENCIES = (0., 0.1, 0.3, 1., 3., 6., 8., 12., 30., 100.)

import argparse
from dataclasses import dataclass
from fractions import Fraction as F
import hashlib
import json
from pathlib import Path
import platform
import numpy as np
import sympy as sp
from scipy.integrate import solve_ivp
import z3

MANUSCRIPT_SHA256='7d9c5c80b49130dad4e00d307cde35fe4ac9fd33868c108b76a6d628ebdb8996'
SPECIES=('A0','B0','A1','B1','A2','B2','C')


class DilutedNetwork:
    """Reaction r consumes species r; columns of P give its products."""
    def __init__(self,names,products):
        self.names=tuple(names);self.P=sp.Matrix(products);n=len(names)
        if self.P.shape!=(n,n) or any(not v.is_Integer or v<0 for v in self.P):raise ValueError('Square nonnegative integer product matrix required.')
        if len(set(names))!=n:raise ValueError('Distinct species names required.')
        self.N=self.P-sp.eye(n)
        if self.N.det()==0:raise ValueError('Stationary-coordinate parametrization requires invertible N.')
    @classmethod
    def paper(cls):
        products=[(1,5),(6,),(3,1),(4,),(5,3),(0,),(2,)]
        P=sp.zeros(7)
        for r,targets in enumerate(products):
            for i in targets:P[i,r]+=1
        return cls(SPECIES,P)
    def monomials(self,x):return sp.Matrix([sp.prod(x[i]**self.P[i,r] for i in range(len(x))) for r in range(len(x))])
    def has_top(self,retained):
        """Exact induced split-graph SCC test for a productive fork."""
        S=set(retained);reach={s:{s}|{i for i in S if self.P[i,s]>0} for s in S}
        for k in S:
            for i in S:
                if k in reach[i]:reach[i]|=reach[k]
        for r in S:
            component={i for i in S if i in reach[r] and r in reach[i]}
            if sum(self.P[i,r] for i in component)>=2:return True
        return False
    def source_audit(self):
        n=len(self.names)
        if n>14:raise ValueError('Exhaustive source-minimality audit limited to fourteen species.')
        proper=[]
        for mask in range(1,2**n-1):
            S=[i for i in range(n) if mask>>i&1]
            if self.has_top(S):proper.append([self.names[i] for i in S])
        return {'full_has_top':self.has_top(range(n)),'proper_subsets_checked':2**n-2,'proper_top_witnesses':proper,
            'source_minimal_for_this_finite_network':self.has_top(range(n)) and not proper}


@dataclass(frozen=True)
class StationaryCoordinates:
    network:DilutedNetwork
    p:sp.Matrix
    q:sp.Matrix
    e:sp.Matrix
    x:sp.Matrix

    def __post_init__(self):
        n=len(self.network.names)
        for key in ('p','q','e','x'):
            value=sp.Matrix([sp.Rational(v) for v in getattr(self,key)]);object.__setattr__(self,key,value)
            if value.shape!=(n,1) or any(v<=0 for v in value):raise ValueError('All stationary coordinates must be positive rational vectors.')
        if self.network.N*(self.p-self.q)!=self.e:raise ValueError('Stationary flux balance N(p-q)=e failed.')
    @classmethod
    def from_positive_parameters(cls,network,x,e,excess_reverse):
        x,e,t=map(sp.Matrix,(x,e,excess_reverse));j=network.N.inv()*e
        if any(v<=0 for v in t):raise ValueError('Positive excess reverse flux required.')
        q=sp.Matrix([max(0,-v)+margin for v,margin in zip(j,t)])
        return cls(network,q+j,q,e,x)
    def rescale(self,c):
        c=sp.Rational(c)
        if c<=0:raise ValueError('Positive concentration scale required.')
        return StationaryCoordinates(self.network,self.p*c,self.q*c,self.e*c,self.x*c)
    @property
    def H(self):
        N=self.network.N;P=self.network.P
        return N*sp.diag(*self.p)-N*sp.diag(*self.q)*P.T-sp.diag(*self.e)
    @property
    def A(self):return self.H*sp.diag(*[1/v for v in self.x])
    @property
    def relative_jacobian(self):return sp.diag(*[1/v for v in self.x])*self.H
    def record(self):return {key:list(map(str,getattr(self,key))) for key in ('p','q','e','x')}


class ReversibleReactor:
    """Full polynomial reactor, with all seven species and every degradation."""
    def __init__(self,network,forward,reverse,degradation):
        self.network=network;self.forward=sp.Matrix(forward);self.reverse=sp.Matrix(reverse);self.degradation=sp.Matrix(degradation)
        n=len(network.names)
        if any(v.shape!=(n,1) for v in (self.forward,self.reverse,self.degradation)) or any(v<=0 for vec in (self.forward,self.reverse) for v in vec) or any(v<0 for v in self.degradation):raise ValueError('Positive reversible rates and nonnegative losses required.')
    @classmethod
    def at_stationary(cls,coordinates):
        c=coordinates;m=c.network.monomials(c.x)
        return cls(c.network,[p/x for p,x in zip(c.p,c.x)],[q/v for q,v in zip(c.q,m)],[e/x for e,x in zip(c.e,c.x)])
    def field(self,z):
        z=sp.Matrix(z);m=self.network.monomials(z)
        return self.network.N*sp.Matrix([k*x-b*y for k,x,b,y in zip(self.forward,z,self.reverse,m)])-sp.diag(*self.degradation)*z
    def jacobian(self,z):
        symbols=sp.Matrix(sp.symbols('x:'+str(len(self.network.names))))
        return self.field(symbols).jacobian(symbols).subs(dict(zip(symbols,z)))
    def simulate_relative(self,coordinates,initial_relative,times,method='Radau'):
        """Integrate y=z/x*, with the original rates held fixed; no elimination.

        Evaluating flux differences around y=1 reduces cancellation near x*.
        The resulting dynamical Jacobian is X^-1 H, similar to A, not H.
        """
        c=coordinates
        if self.field(c.x)!=sp.zeros(len(c.x),1):raise ValueError('Scaling state must be stationary for these rates.')
        reference=ReversibleReactor.at_stationary(c)
        if any(a!=b for v,w in ((self.forward,reference.forward),(self.reverse,reference.reverse),(self.degradation,reference.degradation)) for a,b in zip(v,w)):raise ValueError('Coordinates do not encode this rate assignment.')
        N=np.array(c.network.N,float);P=np.array(c.network.P,int);x=np.array(c.x,float).ravel();p=np.array(c.p,float).ravel();q=np.array(c.q,float).ravel();e=np.array(c.e,float).ravel()
        y0=np.asarray(initial_relative,float);times=np.asarray(times,float)
        if len(y0)!=len(x) or np.any(y0<=0) or times[0]!=0 or np.any(np.diff(times)<=0):raise ValueError('Positive state and increasing times beginning at zero required.')
        def rhs(t,eta):
            # N(p-q)-e=0 exactly, so subtract that constant before evaluation.
            # Accumulate product(1+eta)-1 without subtracting nearby unit values.
            difference=np.zeros(len(x))
            for i in range(len(x)):
                for power in range(int(max(P[i]))):
                    mask=P[i]>power;difference[mask]=difference[mask]*(1+eta[i])+eta[i]
            return (N@(p*eta-q*difference)-e*eta)/x
        sol=solve_ivp(rhs,(0,times[-1]),y0-1,t_eval=times,method=method,rtol=1e-10,atol=2e-15)
        if not sol.success:raise RuntimeError(sol.message)
        if np.any(sol.y<=-1) or not np.all(np.isfinite(sol.y)):raise RuntimeError('Nonpositive numerical trajectory; no clipping applied.')
        return 1+sol.y.T


@dataclass(frozen=True)
class Mode:
    real:sp.Matrix
    imag:sp.Matrix
    growth:sp.Rational=sp.Rational(1)
    frequency:sp.Rational=sp.Rational(8)

    def verify(self,c):
        u,v=sp.Matrix(self.real),sp.Matrix(self.imag);X=sp.diag(*c.x)
        first=c.H*u-X*(self.growth*u-self.frequency*v);second=c.H*v-X*(self.frequency*u+self.growth*v)
        nonzero=any(z!=0 for z in list(u)+list(v));valid=nonzero and first==sp.zeros(len(c.x),1) and second==sp.zeros(len(c.x),1)
        return {'valid':bool(valid),'growth':str(self.growth),'angular_frequency':str(self.frequency),'nonzero_mode':nonzero,
            'real_pencil_residual':list(map(str,first)),'imaginary_pencil_residual':list(map(str,second)),
            'scope':'Exact rational stationary pencil identities. Nonlinear instability follows from the smooth-system linearized-instability theorem; no periodic orbit or Hopf bifurcation is certified.'}
    def relative_solution(self,times,amplitude=1):
        u=np.array(self.real,float).ravel();v=np.array(self.imag,float).ravel();t=np.asarray(times,float)
        return amplitude*np.exp(float(self.growth)*t[:,None])*(np.cos(float(self.frequency)*t[:,None])*u-np.sin(float(self.frequency)*t[:,None])*v)


class InternalResponse:
    def __init__(self,A,internal=6):
        self.A=sp.Matrix(A);self.internal=internal;self.retained=[i for i in range(A.rows) if i!=internal]
        self.Jrr=A.extract(self.retained,self.retained);self.Jri=A.extract(self.retained,[internal]);self.Jir=A.extract([internal],self.retained);self.L=A[internal,internal]
    def effective(self,s):
        if s==self.L:raise ValueError('Internal resolvent has a pole.')
        return self.Jrr+self.Jri*self.Jir/(s-self.L)
    def verify_mode(self,c,mode):
        xi=sp.diag(*c.x)*(mode.real+sp.I*mode.imag);lam=mode.growth+sp.I*mode.frequency;ret=xi.extract(self.retained,[0])
        residual=self.effective(lam)*ret-lam*ret
        return all(sp.simplify(v)==0 for v in residual)
    def diagnostics(self):
        full=np.linalg.eigvals(np.array(self.A,float));static=np.linalg.eigvals(np.array(self.effective(0),float));corner=-float(self.L)
        return {'full_eigenvalues':[[float(v.real),float(v.imag)] for v in full],
            'static_schur_eigenvalues':[[float(v.real),float(v.imag)] for v in static],
            'internal_relaxation_rate':corner,'phase_at_8i_degrees':float(-np.degrees(np.arctan2(8,corner))),
            'resolvent_magnitude_at_8i':float(abs(1/(corner+8j))),
            'scope':'Floating-point spectra and imaginary-axis response diagnostics; the exact eigenpair uses s=1+8i, not s=8i.'}


class ModalDesigner:
    """Exact linear feasibility after fixing a RELATIVE mode and eigenvalue."""
    def __init__(self,timeout_ms=DESIGN_TIMEOUT_MS):self.timeout_ms=timeout_ms
    def solve(self,network,mode):
        n=len(network.names)
        if len(mode.real)!=n or len(mode.imag)!=n or not any(v!=0 for v in list(mode.real)+list(mode.imag)):raise ValueError('A nonzero matching relative mode is required.')
        e=sp.Matrix(sp.symbols('e:'+str(n)));q=sp.Matrix(sp.symbols('q:'+str(n)));x=sp.Matrix(sp.symbols('x:'+str(n)));j=network.N.inv()*e;p=q+j
        H=-network.N*sp.diag(*q)*network.N.T+network.N*sp.diag(*j)-sp.diag(*e);X=sp.diag(*x)
        equations=list(H*mode.real-X*(mode.growth*mode.real-mode.frequency*mode.imag))+list(H*mode.imag-X*(mode.frequency*mode.real+mode.growth*mode.imag))
        symbols=list(e)+list(q)+list(x);variables={s:z3.Real(str(s)) for s in symbols}
        def zq(value):
            value=sp.Rational(value);return z3.RealVal(f'{value.p}/{value.q}')
        def linear(expr):
            poly=sp.Poly(sp.expand(expr),*symbols)
            if poly.total_degree()>1:raise ValueError('Modal design must be linear.')
            return sum(zq(coefficient)*sp_product_z3(powers) for powers,coefficient in poly.terms())
        def sp_product_z3(powers):
            value=z3.RealVal(1)
            for symbol,power in zip(symbols,powers):
                if power:value*=variables[symbol]
            return value
        solver=z3.SolverFor('QF_LRA');solver.set(timeout=self.timeout_ms)
        solver.add(*[linear(v)==0 for v in equations],*[linear(v)>0 for v in list(e)+list(q)+list(x)+list(p)],sum(variables[v] for v in x)==1)
        answer=solver.check()
        if answer==z3.unknown:return {'status':'UNKNOWN','reason':solver.reason_unknown()}
        if answer==z3.unsat:return {'status':'UNSAT_FOR_PRESCRIBED_MODE','scope':'This does not prove stability or absence of a different unstable mode.'}
        model=solver.model()
        def extract(vector):
            values=[]
            for s in vector:
                v=model.eval(variables[s]);values.append(sp.Rational(v.numerator_as_long(),v.denominator_as_long()))
            return sp.Matrix(values)
        ee,qq,xx=map(extract,(e,q,x));coordinates=StationaryCoordinates(network,qq+network.N.inv()*ee,qq,ee,xx)
        if not mode.verify(coordinates)['valid']:raise ArithmeticError('Rational modal-design output failed exact replay.')
        return {'status':'SAT','coordinates':coordinates.record(),'mode_certificate':mode.verify(coordinates),'normalization':'sum(x)=1; rates differ from the published witness'}


def load_certificate():
    data=json.loads(Path(__file__).with_name('certificate_inputs.json').read_text());network=DilutedNetwork.paper()
    coordinates=StationaryCoordinates(network,*[sp.Matrix(data[key]) for key in ('p','q','e','x')])
    mode=Mode(sp.Matrix(data['u'])/MODAL_NORMALIZATION,sp.Matrix(data['v'])/MODAL_NORMALIZATION)
    return data,coordinates,mode


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)
    data,integer,mode=load_certificate();c=integer.rescale(sp.Rational(1,data['D']));reactor=ReversibleReactor.at_stationary(c)
    exact=mode.verify(integer)
    if not exact['valid'] or reactor.field(c.x)!=sp.zeros(7,1) or reactor.jacobian(c.x)!=c.A:raise ArithmeticError('Published stationary/eigenpair certificate failed.')
    response=InternalResponse(c.A);source=c.network.source_audit()
    target_mode=Mode(mode.real,mode.imag,sp.Rational(DESIGN_TARGET_REAL_PART),sp.Rational(DESIGN_TARGET_IMAGINARY_PART))
    design=ModalDesigner().solve(c.network,target_mode)
    times=np.linspace(0,TIME_HORIZON,TIME_SAMPLES);linear=mode.relative_solution(times,INITIAL_RELATIVE_AMPLITUDE);initial=1+linear[0]
    nonlinear=reactor.simulate_relative(c,initial,times);independent=reactor.simulate_relative(c,initial,times,method='BDF')
    half=reactor.simulate_relative(c,1+linear[0]/2,times);half_linear=linear/2
    np.savetxt(out/'trajectories.csv',np.column_stack([times,nonlinear,1+linear,half]),delimiter=',',header=','.join(['time',*[s+'_nonlinear_relative' for s in SPECIES],*[s+'_linear_relative' for s in SPECIES],*[s+'_half_amplitude_relative' for s in SPECIES]]),comments='')
    import csv
    with (out/'exact_rates.csv').open('w',newline='') as f:
        w=csv.writer(f);w.writerow(['reactant_species','x_normalized','forward_rate','reverse_rate','degradation_rate','forward_flux','reverse_flux','degradation_flux'])
        for i,s in enumerate(SPECIES):w.writerow([s,*map(str,(c.x[i],reactor.forward[i],reactor.reverse[i],reactor.degradation[i],c.p[i],c.q[i],c.e[i]))])
    diagnostics=response.diagnostics();H_spectrum=np.linalg.eigvals(np.array(c.H,float));relative_spectrum=np.linalg.eigvals(np.array(c.relative_jacobian,float))
    result={'exact_mode':exact,'source_structure':source,'frequency_response_exact_modal_identity':response.verify_mode(c,mode),
        'numerical_spectrum':diagnostics,'H_spectral_abscissa_numerical':float(max(H_spectrum.real)),
        'true_relative_ODE_jacobian_abscissa_numerical':float(max(relative_spectrum.real)),
        'nonlinear_illustration':{'amplitude':INITIAL_RELATIVE_AMPLITUDE,'horizon':TIME_HORIZON,
            'maximum_relative_departure':float(np.max(abs(nonlinear-1))),'linearization_max_error':float(np.max(abs(nonlinear-1-linear))),
            'half_amplitude_linearization_error':float(np.max(abs(half-1-half_linear))),
            'independent_solver_max_difference':float(np.max(abs(nonlinear-independent)))},
        'design':design,'scope':'Exact instability of the published positive stationary state. Uniqueness invokes the paper family theorem; it is not inferred from a root search. H is a column-scaled derivative, not the Jacobian of the log-coordinate ODE. No periodic orbit or Hopf point is certified; Lean not rerun.'}
    (out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
    with (out/'frequency_response.csv').open('w',newline='') as f:
        w=csv.writer(f);w.writerow(['angular_frequency','scalar_resolvent_real','scalar_resolvent_imag','magnitude','phase_degrees'])
        for omega in FREQUENCIES:
            value=1/(1j*omega-float(response.L));w.writerow([omega,value.real,value.imag,abs(value),np.degrees(np.angle(value))])
    lines=[f'Exact stationary balance and eigenpair 1 +/- 8i: {exact["valid"]}.',
        f'Source minimality: {source["proper_subsets_checked"]} proper species restrictions checked, {len(source["proper_top_witnesses"])} retain (Top).',
        f'Numerical H abscissa: {max(H_spectrum.real):.8g}; actual relative-ODE Jacobian abscissa: {max(relative_spectrum.real):.8g}.',
        f'Internal relaxation rate: {diagnostics["internal_relaxation_rate"]:.8g}; phase at 8i: {diagnostics["phase_at_8i_degrees"]:.6g} degrees.',
        f'Nonlinear/linear maximum difference: {result["nonlinear_illustration"]["linearization_max_error"]:.6g}; half-amplitude difference: {result["nonlinear_illustration"]["half_amplitude_linearization_error"]:.6g}.',
        f'Independent exact modal design: {design["status"]}. Oscillatory instability does not certify a periodic orbit.']
    (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')
    for i in (0,2,3,6):
        line,=axs[0].plot(times,(nonlinear[:,i]-1)/INITIAL_RELATIVE_AMPLITUDE,label=SPECIES[i]);axs[0].plot(times,linear[:,i]/INITIAL_RELATIVE_AMPLITUDE,'--',color=line.get_color(),lw=.8)
    axs[0].set(xlabel='Time (source units)',ylabel='Relative departure / initial amplitude',title='Full nonlinear reactor and linear prediction');axs[0].legend(fontsize=8)
    axs[1].semilogy(times,np.max(abs(nonlinear-1-linear),axis=1)+1e-16,label='Initial amplitude 1e-5');axs[1].semilogy(times,np.max(abs(half-1-half_linear),axis=1)+1e-16,label='Half initial amplitude')
    axs[1].set(xlabel='Time (source units)',ylabel='Maximum linearization error + 1e-16',title='Smaller perturbations follow the linear mode');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'growing_mode.png',dpi=180);fig.savefig(out/'growing_mode.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
    full=np.array(diagnostics['full_eigenvalues']);static=np.array(diagnostics['static_schur_eigenvalues'])
    axs[0].scatter(full[:,0],full[:,1],label='Full physical Jacobian');axs[0].scatter(static[:,0],static[:,1],marker='x',label='Zero-frequency contraction');axs[0].axvline(0,color='gray',lw=.7)
    axs[0].set(xscale='symlog',xlabel='Real part (symmetric logarithmic scale)',ylabel='Imaginary part',title='Eigenvalues of the full and statically\nreduced models');axs[0].legend(fontsize=8)
    omega=np.geomspace(.1,100,300);phase=-np.degrees(np.arctan2(omega,-float(response.L)));axs[1].semilogx(omega,phase);axs[1].axvline(8,color='#bd5a24',ls='--',label='Mode angular frequency 8')
    axs[1].set(xlabel='Angular frequency on the imaginary axis',ylabel='Phase of internal resolvent (degrees)',title='Internal-species response phase versus\nfrequency');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'dynamic_response.png',dpi=180);fig.savefig(out/'dynamic_response.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest();base=Path(__file__).resolve().parent
    (out/'run_metadata.json').write_text(json.dumps({'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':digest(Path(__file__)),'input_sha256':{'certificate_inputs.json':digest(base/'certificate_inputs.json')},'python':platform.python_version(),'z3':z3.get_version_string(),
        'output_sha256':{p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}},indent=2)+'\n')


if __name__=='__main__':main()
Run output
Exact stationary balance and eigenpair 1 +/- 8i: True.
Source minimality: 126 proper species restrictions checked, 0 retain (Top).
Numerical H abscissa: -0.0001179731; actual relative-ODE Jacobian abscissa: 1.
Internal relaxation rate: 6.1230329; phase at 8i: -52.5703 degrees.
Nonlinear/linear maximum difference: 7.34128e-08; half-amplitude difference: 1.83537e-08.
Independent exact modal design: SAT. Oscillatory instability does not certify a periodic orbit.