Example code
This two-strain infection model tracks four compartments: susceptible individuals , individuals carrying only strain 1 or strain 2 (), and coinfected individuals . A strain disappears when both compartments carrying it are empty: or . These overlapping extinction sets are called siphons. The calculation counts their shared compartment once and tests a rare strain's growth in the actual resident population.
The example implements all fourteen mass-action reactions, editable rates, exact resident calculations and stiff numerical trajectories. It connects the structural overlap identity to the paper's permanence criterion: when both single-strain residents exist and each missing strain strictly invades the other resident, every positive solution eventually keeps every compartment above a parameter-dependent positive floor.
In the first sweep, the disease-free reproduction numbers stay fixed at , yet invasion of the strain-1 resident changes sign at . A separate exact example is stronger: changing the superinfection rate leaves the entire disease-free Jacobian unchanged while switching resident invasion from decay to growth. Disease-free data therefore cannot decide resident invasion.
The structural calculation verifies
Here each is the characteristic polynomial describing linear growth or decay for the indicated compartments. Every factor is evaluated at the same common boundary point. The code multiplies polynomials, so the check also works when the shared eigenvalue is zero. It retains additional factors if a supplied normal space includes missing compartments outside the siphons.
The model checks strict invasion with rational covectors and reproduces the paper's independent one-percent box for all fourteen rates. Exact relative growth margins are and . This establishes the theorem's hypotheses at every point in that box; it does not supply one numerical persistence floor shared across the box.
Reusable components also expose the normalized weighted masses used in the proof, the eight-dimensional normalized field, the corrected resident entropy and the conditional conversion from a justified product floor to species floors. The compactness-derived product floor remains unspecified. The plotted trajectory minima are finite observations, and neither convergence in these examples nor an observed minimum is promoted to a general theorem.
Seven scientific test groups check source balances, overlap, threshold cases, invasion certificates, interval bounds, growth identities and independent solver agreement. Inputs are dimensionless manuscript examples, without epidemiological calibration. The package invokes the paper's theorem only within its stated positive-rate and resident-existence hypotheses, and does not rerun Lean.
Python source
"""Coinfection source, overlap algebra, resident invasion and permanence scope."""
from __future__ import annotations
# EDITABLE INPUTS -------------------------------------------------------------
RATE_INPUTS = dict(recruitment='4', alpha1='2', alpha2='1', alpha3='1',
eta1='1/10', eta2='1/10', gamma1='1/10', gamma2='1/10', beta1='1/10', beta2='1',
mu0='1', mu1='1', mu2='1', mu3='1')
INITIAL_STATES = ((0.5, 3.5, 1e-5, 1e-5), (1, 1e-5, 3, 1e-5), (2, .4, .3, .1))
END_TIME = 180.0
SAMPLE_COUNT = 721
UNCERTAINTY_FRACTION = '1/100' # independent box around RATE_INPUTS
MANUSCRIPT_SHA256 = '3f0fc394b76fb14eae5f5e93292c0feea62e542904517384b4f69d86c36d3df3'
# Dimensionless reference rates, not epidemiological calibration or advice.
# No numerical permanence floor is assigned: the paper's floor is existential.
# ---------------------------------------------------------------------------
from dataclasses import dataclass
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
def Q(v): return sp.Rational(str(v))
def texts(v): return [str(sp.simplify(x)) for x in v]
@dataclass(frozen=True)
class Parameters:
recruitment: object = '4'
alpha1: object = '2'
alpha2: object = '1'
alpha3: object = '1'
eta1: object = '1/10'
eta2: object = '1/10'
gamma1: object = '1/10'
gamma2: object = '1/10'
beta1: object = '1/10'
beta2: object = '1'
mu0: object = '1'
mu1: object = '1'
mu2: object = '1'
mu3: object = '1'
def __post_init__(self):
for name in self.__dataclass_fields__:
value=Q(getattr(self,name))
if value<=0: raise ValueError('The literal fourteen-rate theorem requires strictly positive rates.')
object.__setattr__(self,name,value)
def as_json(self): return {name:str(getattr(self,name)) for name in self.__dataclass_fields__}
@dataclass(frozen=True)
class Reaction:
rate_name: str
reactant: tuple[int,...]
product: tuple[int,...]
class ReactionSource:
def __init__(self,species,reactions):
self.species=tuple(species);self.reactions=tuple(reactions);self.n=len(species)
for r in self.reactions:
if len(r.reactant)!=self.n or len(r.product)!=self.n or any(not isinstance(v,int) or v<0 for v in (*r.reactant,*r.product)):
raise ValueError('Complexes must be matching nonnegative integer vectors.')
self.Y=sp.Matrix.hstack(*[sp.Matrix(r.reactant) for r in self.reactions])
self.P=sp.Matrix.hstack(*[sp.Matrix(r.product) for r in self.reactions]);self.S=self.P-self.Y
def symbolic_field(self,rates,x):
flux=sp.Matrix([rates[r.rate_name]*sp.prod(x[i]**r.reactant[i] for i in range(self.n)) for r in self.reactions])
return self.S*flux
def is_siphon(self,indices):
S=set(indices)
if not S.issubset(range(self.n)): raise ValueError('Invalid species indices.')
return all(not any(r.product[i]>r.reactant[i] for i in S) or any(r.reactant[i]>0 for i in S) for r in self.reactions)
def siphons(self,budget=100000):
if 2**self.n>budget: raise RuntimeError('Siphon enumeration limit; no complete list produced.')
return [tuple(i for i in range(self.n) if mask&(1<<i)) for mask in range(2**self.n)
if self.is_siphon(i for i in range(self.n) if mask&(1<<i))]
def overlap_certificate(self,rates,point,S1,S2,normal_indices=None):
if any(rates[r.rate_name]<=0 for r in self.reactions): raise ValueError('Positive source rates required.')
S1,S2=set(S1),set(S2);union=S1|S2;intersection=S1&S2
if not self.is_siphon(S1) or not self.is_siphon(S2): raise ValueError('Both supplied sets must be siphons.')
point=sp.Matrix(point)
if len(point)!=self.n or any(v<0 for v in point) or any(point[i]!=0 for i in union):
raise ValueError('Use a common nonnegative boundary point where both siphons are absent.')
Z=sorted(union if normal_indices is None else set(normal_indices))
if not set(Z).issubset(range(self.n)) or not union.issubset(Z) or any(point[i]!=0 for i in Z): raise ValueError('Normal indices must cover the union and vanish.')
x=sp.Matrix(sp.symbols(f'x0:{self.n}'));z=sp.Symbol('z')
J=self.symbolic_field(rates,x).jacobian(x).subs(dict(zip(x,point)))
signature={i:frozenset(k for k,S in enumerate((S1,S2)) if i in S) for i in Z}
if any(J[i,j]!=0 and not signature[i].issubset(signature[j]) for i in Z for j in Z):
raise AssertionError('Membership containment failed.')
def char(indices):
I=sorted(indices)
return J.extract(I,I).charpoly(z).as_expr() if I else sp.Integer(1)
residual=sp.expand(char(union)*char(intersection)-char(S1)*char(S2))
groups={sig:[i for i in Z if signature[i]==sig] for sig in sorted(set(signature.values()),key=lambda s:(len(s),tuple(s)))}
factor=sp.prod(char(I) for I in groups.values())
if residual!=0 or sp.expand(char(Z)-factor)!=0: raise AssertionError('Characteristic factorization failed.')
return {'normal_indices':Z,'union':sorted(union),'intersection':sorted(intersection),
'normal_matrix':[[str(J[i,j]) for j in Z] for i in Z],
'chi_union':str(sp.factor(char(union))),'chi_intersection':str(sp.factor(char(intersection))),
'chi_S1':str(sp.factor(char(S1))),'chi_S2':str(sp.factor(char(S2))),
'cross_multiplied_residual':'0','chi_full_normal':str(sp.factor(char(Z))),
'signature_factors':[{'membership':sorted(sig),'species':I,'factor':str(sp.factor(char(I)))} for sig,I in groups.items()],
'scope':'All factors evaluated at the same boundary point; no determinant division.'}
class CoinfectionModel:
def __init__(self,parameters=Parameters()):
self.p=parameters
self.source=ReactionSource(('s','a','b','c'),(
Reaction('recruitment',(0,0,0,0),(1,0,0,0)),
Reaction('alpha1',(1,1,0,0),(0,2,0,0)),
Reaction('alpha2',(1,0,1,0),(0,0,2,0)),
Reaction('alpha3',(1,0,0,1),(0,0,0,2)),
Reaction('eta1',(0,1,0,1),(0,0,0,2)),
Reaction('eta2',(0,0,1,1),(0,0,0,2)),
Reaction('gamma1',(0,1,1,0),(0,0,1,1)),
Reaction('gamma2',(0,1,1,0),(0,1,0,1)),
Reaction('beta1',(1,0,0,1),(0,1,0,1)),
Reaction('beta2',(1,0,0,1),(0,0,1,1)),
Reaction('mu0',(1,0,0,0),(0,0,0,0)),
Reaction('mu1',(0,1,0,0),(0,0,0,0)),
Reaction('mu2',(0,0,1,0),(0,0,0,0)),
Reaction('mu3',(0,0,0,1),(0,0,0,0))))
self.x=sp.Matrix(sp.symbols('s a b c'))
self.rates={name:getattr(parameters,name) for name in parameters.__dataclass_fields__}
self.F=self.source.symbolic_field(self.rates,self.x)
self.J=self.F.jacobian(self.x)
self._F=sp.lambdify([tuple(self.x)],self.F,'numpy');self._J=sp.lambdify([tuple(self.x)],self.J,'numpy')
def residents(self):
p=self.p;s0=p.recruitment/p.mu0
rows=[sp.Matrix([s0,0,0,0])]
for i in (1,2):
mu=getattr(p,f'mu{i}');alpha=getattr(p,f'alpha{i}');s=mu/alpha;u=(p.recruitment-p.mu0*s)/mu
state=sp.Matrix([s,0,0,0]);state[i]=u
rows.append(state if u>0 else None)
return rows
def invasion_blocks(self):
_,E1,E2=self.residents()
if E1 is None or E2 is None: raise ValueError('Both positive single-strain residents must exist.')
return (self.J.subs(dict(zip(self.x,E1))).extract([2,3],[2,3]),
self.J.subs(dict(zip(self.x,E2))).extract([1,3],[1,3]))
def bounds(self):
p=self.p;R=p.recruitment/min(p.mu0,p.mu1,p.mu2,p.mu3)+1
rate=p.mu0+(p.alpha1+p.alpha2+p.alpha3+p.beta1+p.beta2)*R
return R,p.recruitment/(2*rate)
def field(self,x): return np.asarray(self._F(x),float).reshape(-1)
def simulate(self,initial,times,method='Radau'):
initial=np.asarray(initial,float);times=np.asarray(times,float)
if initial.shape!=(4,) or np.any(initial<=0) or not np.all(np.isfinite(initial)):
raise ValueError('Positive finite four-compartment initial state required.')
if len(times)<2 or times[0]!=0 or np.any(np.diff(times)<=0): raise ValueError('Times must increase from zero.')
# Log coordinates retain positivity without clipping compartments near extinction.
def rhs(t,z):
state=np.exp(z);return self.field(state)/state
def jac(t,z):
state=np.exp(z);J=np.asarray(self._J(state),float)
return J*state[None,:]/state[:,None]-np.diag(self.field(state)/state)
sol=solve_ivp(rhs,(0,times[-1]),np.log(initial),t_eval=times,jac=jac,method=method,rtol=2e-9,atol=2e-11)
states=np.exp(sol.y.T)
if not sol.success or not np.all(np.isfinite(states)) or np.any(states<=0): raise ArithmeticError('Log-state integration failed.')
return states
def permanence(self):
residents=self.residents()
if any(x is None for x in residents):
return {'status':'outside_resident_existence_hypotheses','numerical_floor':None}
invasion=[InvasionBlock(M).certificate() for M in self.invasion_blocks()]
statuses=[r['status'] for r in invasion]
status=('uniform_permanence_by_paper_theorem' if statuses==['strict_invasion']*2 else
'permanence_fails_resident_sink' if 'strict_decay' in statuses else 'threshold_not_decided')
return {'status':status,'residents':[texts(x) for x in residents],'invasion':invasion,
'numerical_floor':None,'scope':'Existential floor depends on fixed parameters, not initial state. No global interior convergence claim.'}
class InvasionBlock:
def __init__(self,M):
self.M=sp.Matrix(M)
if self.M.shape!=(2,2) or self.M[0,1]<=0 or self.M[1,0]<=0:
raise ValueError('Irreducible two-dimensional Metzler block required.')
def certificate(self):
A,B,C,D=list(self.M)
eigen=(A+D+sp.sqrt((A-D)**2+4*B*C))/2
if A>=0 or D>=0 or A*D<B*C:
status='strict_invasion'
r=1+max(0,-A/C) if D>=0 else (max(0,-A/C)+B/(-D))/2
elif A<0 and D<0 and B*C<A*D:
status='strict_decay';r=(B/(-D)+(-A)/C)/2
else:
return {'status':'threshold','matrix':[texts(self.M.row(i)) for i in range(2)],'lambda_plus':'0','numeric_lambda_plus':0.0}
col=sp.Matrix([[1,r]])*self.M
margin=min(col[0],col[1]/r) if status=='strict_invasion' else min(-col[0],-col[1]/r)
if margin<=0: raise AssertionError('Exact covector sign failed.')
return {'status':status,'matrix':[texts(self.M.row(i)) for i in range(2)],'left_weight':str(r),
'weighted_columns':texts(col),'relative_margin':str(margin),'lambda_plus':str(eigen),'numeric_lambda_plus':float(eigen)}
class NormalizedLift:
"""Polynomial extension of weighted relative growth to extinction faces.
Arbitrary directions above a boundary point need not lie in the physical
compact closure K. The certificate bounds larger direction simplices.
"""
def __init__(self,model):
self.model=model;cert=model.permanence()
if cert['status']!='uniform_permanence_by_paper_theorem': raise ValueError('Strict mutual invasion and resident existence required.')
self.ru=Q(cert['invasion'][1]['left_weight']);self.rv=Q(cert['invasion'][0]['left_weight'])
p=model.p;s0=p.recruitment/p.mu0
a0=p.alpha1*s0-p.mu1;b0=p.alpha2*s0-p.mu2;d0=p.alpha3*s0-p.mu3
self.rj=sp.Integer(1) if d0>=0 else min(1,(p.beta1+p.beta2)*s0/(-2*d0))
delta0=min(a0,b0,(p.beta1+p.beta2)*s0/self.rj+d0)
L0=min(a0,p.beta1*s0/self.ru+d0)+min(b0,p.beta2*s0/self.rv+d0)
self.k=max(1,int(sp.floor(-L0/delta0))+1)
self.fibre_margins=(L0+self.k*delta0,Q(cert['invasion'][0]['relative_margin']),Q(cert['invasion'][1]['relative_margin']))
if min(self.fibre_margins)<=0: raise AssertionError('Positive fibre margins failed.')
def lift(self,state):
s,a,b,c=state;U=a+self.ru*c;V=b+self.rv*c;J=a+b+self.rj*c
if min(U,V,J)<=0: raise ValueError('Undefined direction on extinction face; supply the extended coordinates explicitly.')
return list(state)+[a/U,b/V,a/J,b/J]
def growth(self,extended):
s,a,b,c,theta,phi,xi,zeta=extended;p=self.model.p;h=p.gamma1+p.gamma2
A=p.alpha1*s-p.gamma1*b-p.eta1*c-p.mu1
B=p.alpha2*s-p.gamma2*a-p.eta2*c-p.mu2
C=p.eta1*a+p.eta2*b+p.alpha3*s-p.mu3
HU=(A+self.ru*h*b)*theta+(p.beta1*s+self.ru*C)*(1-theta)/self.ru
HV=(B+self.rv*h*a)*phi+(p.beta2*s+self.rv*C)*(1-phi)/self.rv
HJ=A*xi+(B+self.rj*h*a)*zeta+((p.beta1+p.beta2)*s+self.rj*C)*(1-xi-zeta)/self.rj
return HU,HV,HJ,HU+HV+self.k*HJ
def field(self,extended):
s,a,b,c,theta,phi,xi,zeta=extended;p=self.model.p
HU,HV,HJ,_=self.growth(extended)
A=p.alpha1*s-p.gamma1*b-p.eta1*c-p.mu1;B=p.alpha2*s-p.gamma2*a-p.eta2*c-p.mu2
physical=self.model.F.subs(dict(zip(self.model.x,(s,a,b,c))))
return sp.Matrix(list(physical)+[A*theta+p.beta1*s/self.ru*(1-theta)-theta*HU,
B*phi+p.beta2*s/self.rv*(1-phi)-phi*HV,
A*xi+p.beta1*s/self.rj*(1-xi-zeta)-xi*HJ,
B*zeta+p.beta2*s/self.rj*(1-xi-zeta)-zeta*HJ])
def conditional_recovery(self,q):
"""Algebra conditional on an externally justified eventual product floor.
A supplied q is an assumption, not verified or estimated here.
"""
q=Q(q)
if q<=0: raise ValueError('Positive assumed product floor required.')
p=self.model.p;R,ell=self.model.bounds()
mu=q/((1+self.rv)*R*((2+self.rj)*R)**self.k)
mv=q/((1+self.ru)*R*((2+self.rj)*R)**self.k)
Ca=p.beta1*ell*mu/self.ru;Ka=p.mu1+(p.gamma1+p.eta1)*R+p.beta1*ell/self.ru
Cb=p.beta2*ell*mv/self.rv;Kb=p.mu2+(p.gamma2+p.eta2)*R+p.beta2*ell/self.rv
af=Ca/(2*Ka);bf=Cb/(2*Kb);cf=(p.gamma1+p.gamma2)*af*bf/(2*p.mu3)
return {'assumed_q':str(q),'conditional_only':True,'U_floor':str(mu),'V_floor':str(mv),
'species_floors':texts([ell,af,bf,cf]),'epsilon':str(min(ell,af,bf,cf))}
@dataclass(frozen=True)
class Interval:
lo: object
hi: object
def __post_init__(self):
object.__setattr__(self,'lo',Q(self.lo));object.__setattr__(self,'hi',Q(self.hi))
if self.lo>self.hi: raise ValueError('Reversed interval.')
@staticmethod
def coerce(x): return x if isinstance(x,Interval) else Interval(x,x)
def __add__(self,x):
x=self.coerce(x);return Interval(self.lo+x.lo,self.hi+x.hi)
__radd__=__add__
def __neg__(self):return Interval(-self.hi,-self.lo)
def __sub__(self,x):return self+-self.coerce(x)
def __mul__(self,x):
x=self.coerce(x);v=[a*b for a in (self.lo,self.hi) for b in (x.lo,x.hi)]
return Interval(min(v),max(v))
__rmul__=__mul__
def __truediv__(self,x):
x=self.coerce(x)
if x.lo<=0<=x.hi:raise ValueError('Division interval contains zero.')
return self*Interval(1/x.hi,1/x.lo)
def pair(self):return texts([self.lo,self.hi])
class ParameterBox:
def __init__(self,center,radius='1/100'):
radius=Q(radius)
if not 0<=radius<1:raise ValueError('Relative box radius must be in [0,1).')
self.box={n:Interval(getattr(center,n)*(1-radius),getattr(center,n)*(1+radius)) for n in center.__dataclass_fields__}
def certificate(self,weights=(2,1)):
p=self.box;R=p['recruitment'];mu0=p['mu0']
s1=p['mu1']/p['alpha1'];s2=p['mu2']/p['alpha2']
u1=(R-mu0*s1)/p['mu1'];u2=(R-mu0*s2)/p['mu2'];s0=R/mu0
rv,ru=map(Q,weights)
if min(rv,ru)<=0:raise ValueError('Positive covector weights required.')
# Expand the weighted columns BEFORE interval arithmetic to keep cancellation.
c11=p['alpha2']*s1-p['mu2']+(rv*p['gamma1']+(rv-1)*p['gamma2'])*u1
c12=p['beta2']*s1/rv+p['eta1']*u1+p['alpha3']*s1-p['mu3']
c21=p['alpha1']*s2-p['mu1']+((ru-1)*p['gamma1']+ru*p['gamma2'])*u2
c22=p['beta1']*s2/ru+p['eta2']*u2+p['alpha3']*s2-p['mu3']
margins=[min(c11.lo,c12.lo),min(c21.lo,c22.lo)]
existence=s1.hi<s0.lo and s2.hi<s0.lo and min(u1.lo,u2.lo)>0
return {'certified':bool(existence and min(margins)>0),'resident_s_intervals':[s1.pair(),s2.pair()],
'resident_u_intervals':[u1.pair(),u2.pair()],'weights':texts([rv,ru]),
'weighted_relative_column_intervals':[c11.pair(),c12.pair(),c21.pair(),c22.pair()],
'relative_margins':texts(margins),
'scope':'Independent continuous parameter box; permanence at each parameter, no common numerical floor over the box. Failure is inconclusive.'}
class ResidentEntropy:
def __init__(self,model,strain=1):
if strain not in (1,2):raise ValueError('Strain must be 1 or 2.')
resident=model.residents()[strain]
if resident is None:raise ValueError('Positive resident required.')
self.p=model.p;self.sbar=resident[0];self.ubar=resident[strain]
self.alpha=getattr(self.p,f'alpha{strain}');self.mu=getattr(self.p,f'mu{strain}')
self.R,self.ell=model.bounds();self.g=self.p.mu0+self.alpha*self.ubar
self.epsilon=min(self.g/(4*self.alpha*self.R**2),self.alpha*self.ell/(2*self.g*self.R))/2
def value_and_drift(self,s,u):
if min(s,u)<=0:raise ValueError('Positive face concentrations required.')
v=s-float(self.sbar);w=u-float(self.ubar);e=float(self.epsilon);g=float(self.g);alpha=float(self.alpha)
H=v-float(self.sbar)*math.log(s/float(self.sbar))+w-float(self.ubar)*math.log(u/float(self.ubar))
drift=-(g/s-e*alpha*u)*v*v-e*g*v*w-e*alpha*s*w*w
return H+e*v*w,drift
def write_csv(out,name,headers,rows):
with (out/name).open('w',newline='') as f:
writer=csv.writer(f);writer.writerow(headers);writer.writerows(rows)
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)
model=CoinfectionModel(Parameters(**RATE_INPUTS));p=model.p;residents=model.residents()
cert=model.permanence();lift=NormalizedLift(model) if cert['status']=='uniform_permanence_by_paper_theorem' else None
overlap=model.source.overlap_certificate(model.rates,residents[0],(1,3),(2,3))
threshold_overlap=model.source.overlap_certificate(model.rates,[p.mu3/p.alpha3,0,0,0],(1,3),(2,3))
# Fixed published obstruction variants, independent of user-edited primary inputs.
low=CoinfectionModel(Parameters(beta2='1/10'))
eta=CoinfectionModel(Parameters(beta2='1/10',eta1='1/5'))
dfe_equal=low.J.subs(dict(zip(low.x,low.residents()[0])))==eta.J.subs(dict(zip(eta.x,eta.residents()[0])))
result={'parameters':p.as_json(),'siphons':[list(s) for s in model.source.siphons()],
'permanence':cert,'overlap':overlap,'zero_shared_mode_overlap':threshold_overlap,
'independent_parameter_box':ParameterBox(p,UNCERTAINTY_FRACTION).certificate(),
'fixed_obstruction':{'complete_DFE_jacobians_equal_under_eta1_change':dfe_equal,
'low_beta2':low.permanence(),'changed_eta1':eta.permanence(),'beta2_threshold':'51/140'},
'normalized_weights':None if lift is None else {'ru':str(lift.ru),'rv':str(lift.rv),'rj':str(lift.rj),'k':lift.k,
'fibre_growth_lower_bounds':texts(lift.fibre_margins),'product_floor':None},
'absorbing_total_bound_and_eventual_s_floor':texts(model.bounds()),
'scope':'Fixed literal positive-rate model; theorem floor existential, sampled minima never promoted to a guarantee. No Lean rerun.'}
(out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
times=np.linspace(0,END_TIME,SAMPLE_COUNT);trajectories=[];curves=[]
for i,initial in enumerate(INITIAL_STATES):
states=model.simulate(initial,times);curves.append(states)
trajectories.extend([[i,t,*row] for t,row in zip(times,states)])
write_csv(out,'trajectories.csv',['initial_index','time','s','a','b','c'],trajectories)
comparison=low.simulate(INITIAL_STATES[0],times)
write_csv(out,'resident_sink_trajectory.csv',['time','s','a','b','c'],[[t,*row] for t,row in zip(times,comparison)])
sweep=[]
for beta in np.linspace(.02,1.2,181):
variant=CoinfectionModel(Parameters(beta2=str(beta)));blocks=variant.invasion_blocks()
c1,c2=[InvasionBlock(M).certificate() for M in blocks]
sweep.append([beta,c1['numeric_lambda_plus'],c2['numeric_lambda_plus'],c1['status']])
write_csv(out,'resident_invasion_sweep.csv',['beta2','lambda_at_E1','lambda_at_E2','E1_invasion_status'],sweep)
# Unequal mortality face example uses an independent published-model parameter choice.
face=CoinfectionModel(Parameters(mu0='6/5',mu1='13/10'))
entropy=ResidentEntropy(face);ft=np.linspace(0,30,301)
def face_rhs(t,z):
s,u=np.exp(z);ds=float(face.p.recruitment)-float(face.p.mu0)*s-float(face.p.alpha1)*s*u
du=float(face.p.alpha1)*s*u-float(face.p.mu1)*u
return [ds/s,du/u]
sol=solve_ivp(face_rhs,(0,30),np.log([2,.2]),t_eval=ft,method='Radau',rtol=1e-10,atol=1e-12)
if not sol.success:raise ArithmeticError('Resident face integration failed.')
face_states=np.exp(sol.y.T)
write_csv(out,'resident_entropy.csv',['time','s','a','corrected_entropy','exact_formula_drift'],
[[t,*row,*entropy.value_and_drift(*row)] for t,row in zip(ft,face_states)])
lines=[f'Permanence classification: {cert["status"]}.',
f'Siphons (species indices s,a,b,c): {result["siphons"]}.',
f'Overlap identity residual: {overlap["cross_multiplied_residual"]}; at shared zero mode: {threshold_overlap["cross_multiplied_residual"]}.',
f'Independent parameter-box margins: {result["independent_parameter_box"]["relative_margins"]}.',
f'Complete disease-free Jacobian unchanged under eta1 obstruction: {dfe_equal}.',
f'Normalized weights: {result["normalized_weights"]}.',
'No numerical universal species floor is computed or inferred from trajectories.']
(out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
sweep_array=np.array([r[:3] for r in sweep])
axes[0].plot(sweep_array[:,0],sweep_array[:,1],label='Missing strain 2 at resident E1')
axes[0].plot(sweep_array[:,0],sweep_array[:,2],label='Missing strain 1 at resident E2')
axes[0].axhline(0,color='black',lw=.8);axes[0].axvline(51/140,color='#777',ls='--',label='Exact threshold 51/140')
axes[0].set(xlabel='Single-strain transmission from coinfection beta2',ylabel='Dominant invasion eigenvalue',title='Resident invasion at fixed disease-free\nreproduction numbers')
axes[0].legend(fontsize=7)
for i,name in enumerate(model.source.species):
axes[1].semilogy(times,curves[0][:,i],label=name)
axes[1].set(xlabel='Time',ylabel='Compartment density',title='Recovery from near the strain-1 resident')
axes[1].legend(fontsize=8)
for ax in axes:ax.grid(alpha=.2)
fig.savefig(out/'invasion.png',dpi=180);fig.savefig(out/'invasion.svg');plt.close(fig)
fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
for i,states in enumerate(curves):axes[0].semilogy(times,np.min(states,axis=1),label=f'Initial state {i+1}')
axes[0].semilogy(times,np.min(comparison,axis=1),'--',color='black',label='beta2=0.1 near resident sink')
axes[0].set(xlabel='Time',ylabel='Smallest observed compartment',title='Minimum compartment density along sample\ntrajectories');axes[0].legend(fontsize=7)
values=np.array([entropy.value_and_drift(*row) for row in face_states])
axes[1].plot(ft,values[:,0],label='Corrected relative entropy')
axes[1].set(xlabel='Time on the single-strain face',ylabel='Corrected entropy',title='Resident convergence with unequal loss rates')
axes[1].legend(fontsize=8)
for ax in axes:ax.grid(alpha=.2)
fig.savefig(out/'persistence.png',dpi=180);fig.savefig(out/'persistence.svg');plt.close(fig)
digest=lambda path:hashlib.sha256(path.read_bytes()).hexdigest()
(out/'run_metadata.json').write_text(json.dumps({'manuscript_sha256':MANUSCRIPT_SHA256,
'source_sha256':digest(Path(__file__)),'python':platform.python_version(),
'output_sha256':{f.name:digest(f) for f in sorted(out.iterdir()) if f.is_file() and f.name!='run_metadata.json'}},indent=2)+'\n')
if __name__=='__main__':main()
Run output
Permanence classification: uniform_permanence_by_paper_theorem.
Siphons (species indices s,a,b,c): [[], [1, 3], [2, 3], [1, 2, 3]].
Overlap identity residual: 0; at shared zero mode: 0.
Independent parameter-box margins: ['5601/101000', '8701/25250'].
Complete disease-free Jacobian unchanged under eta1 obstruction: True.
Normalized weights: {'ru': '1', 'rv': '191/84', 'rj': '1', 'k': 1, 'fibre_growth_lower_bounds': ['47/5', '267/3820', '2/5'], 'product_floor': None}.
No numerical universal species floor is computed or inferred from trajectories.