"""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()