"""Inherited-state branching: exact witness arithmetic and reusable count models.

Constant rates, faithful same-type division, independent founder families.
Prediction coverage concerns future counts, not uncertainty in fitted parameters.
"""
from __future__ import annotations
import argparse
import csv
from dataclasses import dataclass,asdict,replace
from fractions import Fraction as Q
import hashlib
import json
import math
from pathlib import Path
import platform
import numpy as np
from scipy.integrate import solve_ivp
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import expm_multiply
from scipy.stats import binom,norm
import sympy as sp

# EDITABLE INPUTS -----------------------------------------------------------
BIRTH_RATES = (0., .1)                  # sensitive, resistant; per day
DEATH_RATES = (.3, 0.)
SWITCH_RATES = (.001, .00001)            # S->R, R->S; per day
RESISTANT_FOUNDER_PROBABILITY = 5/24      # drawn once per founder, not at each division
HORIZON_DAYS = math.log(2)/.10001
TARGET_COVERAGE = .95
FOUNDER_COUNTS = (1,2,3,5,10,20,50,100)
TERMINAL_COEFFICIENT_CAP = 192           # low PGF degrees; NOT a population cap
FORWARD_POPULATION_CAP = 48              # independent forward check retains overflow
ALPHA_GRID = (.0001,.001,.01,.1,1.)
BETA_GRID = (.00001,.0001,.001,.01,.1,1.)
RTOL,ATOL = 2e-12,2e-14
RANDOM_SEED = 480
EVENT_BUDGET = 100000
MANUSCRIPT_SHA256 = 'a30348c7b611dc86e0435ca9aa69eeaafaeb26bd6992488937114d9f71f78f1d'
# --------------------------------------------------------------------------
SWITCH_COV=np.array([[1.,-1.],[-1.,1.]])


@dataclass(frozen=True)
class Rates:
    birth: tuple = BIRTH_RATES
    death: tuple = DEATH_RATES
    switching: tuple = SWITCH_RATES
    def __post_init__(self):
        if any(len(v)!=2 for v in (self.birth,self.death,self.switching)) or not np.all(np.isfinite([self.birth,self.death,self.switching])) or min(*self.birth,*self.death,*self.switching)<0:
            raise ValueError('Two finite nonnegative constant rates for each event class required.')
    @property
    def drift(self):
        a,b=self.switching;g=np.array(self.birth)-self.death
        return np.array([[g[0]-a,b],[a,g[1]-b]])


@dataclass(frozen=True)
class FounderPreparation:
    potential: int = 1
    plating: float = 1.
    resistant: float = RESISTANT_FOUNDER_PROBABILITY
    def __post_init__(self):
        if not isinstance(self.potential,int) or isinstance(self.potential,bool) or self.potential<0 or not np.isfinite(self.plating) or not np.isfinite(self.resistant) or not 0<=self.plating<=1 or not 0<=self.resistant<=1:
            raise ValueError('Nonnegative integer potential founders and probabilities in [0,1] required.')
    def moments(self,single_mean,single_variance):
        return self.potential*self.plating*single_mean,self.potential*(self.plating*single_variance+self.plating*(1-self.plating)*single_mean**2)
    def compose(self,single,cap):
        """Truncated polynomial powers preserve low coefficients without renormalizing.

        Each potential founder is absent with probability 1-f; otherwise its type
        was drawn independently. This is not fixed counts of each initial type.
        """
        base=self.plating*np.asarray(single,float).copy();base[0]+=1-self.plating
        result=np.array([1.]);power=self.potential
        while power:
            if power&1:result=np.convolve(result,base)[:cap+1]
            power//=2
            if power:base=np.convolve(base,base)[:cap+1]
        return np.pad(result,(0,max(0,cap+1-len(result))))[:cap+1]


class MomentClosure:
    def __init__(self,rates,resistant):
        if not np.isfinite(resistant) or not 0<=resistant<=1:raise ValueError('Founder type probability in [0,1] required.')
        self.rates=rates;self.p=resistant
    def solve(self,duration):
        if not np.isfinite(duration) or duration<=0:raise ValueError('Positive finite horizon required.')
        rates=self.rates;M=rates.drift;birth=np.array(rates.birth);death=np.array(rates.death);growth=birth-death;switch=np.array(rates.switching)
        def rhs(t,y):
            mu=y[:2];cov=y[2:6].reshape(2,2);mean=mu.sum()
            if mean<=0:raise ArithmeticError('Mean underflow; rescale this numerical problem.')
            gbar=growth@mu/mean;noise=np.diag((birth+death)*mu)+(switch@mu)*SWITCH_COV
            covdot=M@cov+cov@M.T+noise
            excess=float(np.ones(2)@cov@(growth-gbar*np.ones(2)))
            # int(lambda/m), scalar variance, integrated covariance / m^2,
            # and an independently integrated scalar mean.
            return np.r_[M@mu,covdot.ravel(),(birth@mu)/mean**2,2*gbar*y[7]+(birth+death)@mu,
                         excess/mean**2,gbar*y[9]]
        y0=np.r_[[1-self.p,self.p],(self.p*(1-self.p)*SWITCH_COV).ravel(),0.,0.,0.,1.]
        t=np.linspace(0,duration,301);sol=solve_ivp(rhs,(0,duration),y0,t_eval=t,rtol=RTOL,atol=ATOL,method='DOP853')
        if not sol.success:raise RuntimeError(sol.message)
        final=sol.y[:,-1];mean=final[:2].sum();A=mean*final[6];vn=final[2:6].sum();vz=final[7]
        assert abs(vz-(mean*(1+2*A)-mean**2))<1e-8*max(1.,abs(vz))
        defect=2*mean**2*final[8];difference=abs(defect-(vn-vz))
        traces=[]
        for ti,y in zip(t,sol.y.T):
            mu=y[:2];cov=y[2:6].reshape(2,2);m=mu.sum();gb=growth@mu/m;excess=float(np.ones(2)@cov@(growth-gb*np.ones(2)))
            traces.append((ti,*mu,m,y[9],cov.sum(),y[7],excess,2*m*m*y[8]))
        return dict(mean=float(mean),A=float(A),variance_actual=float(vn),variance_scalar=float(vz),
            scalar_over_actual=float(vz/vn) if vn>0 else None,founder_limit=float(2*norm.cdf(norm.ppf((1+TARGET_COVERAGE)/2)*math.sqrt(vz/vn))-1) if vn>0 else None,
            defect_identity_error=float(difference),mean_identity_error=float(max(abs(sol.y[:2].sum(axis=0)-sol.y[9]))),trace=traces)


class ScalarTerminalLaw:
    """Kendall law for one founder of the mean-composition closure."""
    def __init__(self,mean,birth_parameter):
        self.m=float(mean);self.A=float(birth_parameter)
        if not np.isfinite(self.m+self.A) or self.m<0 or self.A<0 or self.m>1+self.A+1e-10:raise ValueError('Invalid mean / accumulated-birth parameter.')
    def coefficients(self,cap):
        p=np.empty(cap+1);p[0]=1-self.m/(1+self.A)
        p[1:]=self.m/(1+self.A)**2*(self.A/(1+self.A))**np.arange(cap)
        return p
    def tail_above(self,k):return self.m/(1+self.A)*(self.A/(1+self.A))**k


class BackwardCountLaw:
    """Closed PGF coefficients; grow-then-shrink histories remain represented."""
    def __init__(self,rates):self.rates=rates
    def solve(self,duration,cap):
        if not isinstance(cap,int) or cap<1 or duration<=0 or not np.isfinite(duration):raise ValueError('Degree cap >=1 and positive finite duration required.')
        def rhs(t,v):
            q=v.reshape(2,cap+1);out=np.empty_like(q)
            for i in range(2):
                out[i]=self.rates.birth[i]*(np.convolve(q[i],q[i])[:cap+1]-q[i])-self.rates.death[i]*q[i]+self.rates.switching[i]*(q[1-i]-q[i])
                out[i,0]+=self.rates.death[i]
            return out.ravel()
        initial=np.zeros((2,cap+1));initial[:,1]=1
        sol=solve_ivp(rhs,(0,duration),initial.ravel(),rtol=RTOL,atol=ATOL,method='DOP853')
        if not sol.success:raise RuntimeError(sol.message)
        p=sol.y[:,-1].reshape(2,cap+1)
        if p.min()<-1e-10:raise ArithmeticError('Negative numerical coefficient; no clipping applied.')
        return p


class ForwardOverflowLaw:
    """Independent joint forward solver, absorbing all exits and initial excess."""
    def __init__(self,rates,cap):self.rates=rates;self.cap=cap
    def solve(self,duration,preparation):
        if not isinstance(self.cap,int) or self.cap<1 or not np.isfinite(duration) or duration<0:raise ValueError('Positive integer population cap and nonnegative finite duration required.')
        states=[(s,n-s) for n in range(self.cap+1) for s in range(n+1)];index={v:i for i,v in enumerate(states)};overflow=len(states)
        rr=[];cc=[];vv=[];rates=self.rates
        for j,(s,r) in enumerate(states):
            events=[((s+1,r),rates.birth[0]*s),((s,r+1),rates.birth[1]*r),((s-1,r),rates.death[0]*s),((s,r-1),rates.death[1]*r),
                ((s-1,r+1),rates.switching[0]*s),((s+1,r-1),rates.switching[1]*r)]
            for dest,rate in events:
                if rate:rr.extend([index.get(dest,overflow),j]);cc.extend([j,j]);vv.extend([rate,-rate])
        G=coo_matrix((vv,(rr,cc)),shape=(overflow+1,overflow+1)).tocsc();p=preparation
        initial=np.array([binom.pmf(s+r,p.potential,p.plating)*binom.pmf(r,s+r,p.resistant) for s,r in states]+[binom.sf(self.cap,p.potential,p.plating)])
        final=expm_multiply(duration*G,initial)
        law=np.bincount([s+r for s,r in states],weights=final[:-1],minlength=self.cap+1)
        return dict(pmf=law,overflow=float(final[-1]),mass_residual=float(abs(final.sum()-1)),
            scope='Absorbing overflow supplies an exact-arithmetic enclosure mechanism; floating exponential and residuals are not validated roundoff bounds.')


class PredictionAssessment:
    def __init__(self,target=TARGET_COVERAGE):
        if not 0<target<1:raise ValueError('Coverage target in (0,1) required.')
        self.target=target
    @staticmethod
    def quantile(pmf,probability):
        k=int(np.searchsorted(np.cumsum(pmf),probability))
        if k>=len(pmf):raise ArithmeticError('Coefficient cap does not resolve this quantile; increase it.')
        return k
    def assess(self,actual,scalar):
        lo=self.quantile(scalar,(1-self.target)/2);hi=self.quantile(scalar,(1+self.target)/2)
        repair=self.quantile(actual,self.target)
        return dict(scalar_interval=[lo,hi],scalar_own_coverage=float(sum(scalar[lo:hi+1])),actual_coverage=float(sum(actual[lo:hi+1])),
            actual_one_sided_endpoint=repair,actual_one_sided_coverage=float(sum(actual[:repair+1])),
            actual_previous_coverage=float(sum(actual[:repair])),actual_mass_deficit=float(1-sum(actual)),scalar_mass_deficit=float(1-sum(scalar)),
            scope='Numerical quantiles and coverage, not a validated probability enclosure. One-sided repair is a distinct interval convention.')


class EventSimulator:
    def __init__(self,rates):self.rates=rates
    def run(self,preparation,duration,seed=RANDOM_SEED,event_budget=EVENT_BUDGET):
        if not np.isfinite(duration) or duration<0 or not isinstance(event_budget,int) or event_budget<0:raise ValueError('Nonnegative finite duration and integer event budget required.')
        rng=np.random.default_rng(seed);n=rng.binomial(preparation.potential,preparation.plating);r=rng.binomial(n,preparation.resistant);s=n-r;t=0.;events=0
        while t<duration:
            prop=np.array([self.rates.birth[0]*s,self.rates.birth[1]*r,self.rates.death[0]*s,self.rates.death[1]*r,self.rates.switching[0]*s,self.rates.switching[1]*r])
            rate=prop.sum()
            if not rate:return dict(finished=True,time=duration,S=int(s),R=int(r),events=events)
            wait=rng.exponential(1/rate)
            if t+wait>duration:return dict(finished=True,time=duration,S=int(s),R=int(r),events=events)
            if events>=event_budget:return dict(finished=False,time=t,S=int(s),R=int(r),events=events,reason='event budget; no terminal-count verdict')
            t+=wait;j=int(rng.choice(6,p=prop/rate));ds,dr=[(1,0),(0,1),(-1,0),(0,-1),(-1,1),(1,-1)][j];s+=ds;r+=dr;events+=1
        return dict(finished=True,time=t,S=int(s),R=int(r),events=events)


def exact_witness():
    c=Q(10000,10001);pr=Q(5,24);ps=1-pr
    L=pr/2*sum((c/2)**j for j in range(2,7));H=ps*Q(300,301)+pr/2*sum((c/2)**j for j in range(3))
    geo=pr/2*(c/2)**2/(1-c/2);variance=pr/2*sum(n*n*(c/2)**(n-1) for n in range(1,8))-Q(11,20)**2
    assert L==Q(151415045117539070312500,3001800450060004500180003)>Q(1,20)
    assert H==Q(467880212635,481696324816)>Q(9713,10000) and variance>Q(4,5)
    x=sp.Symbol('x');D=5*(1+x)**4+19;U=120*(1+x)**6+1368*(1+x)**2
    P=sp.Rational(25897,10000)+sp.Rational(3457,2000)*x-sp.Rational(32099,5000)*x*x+sp.Rational(179,50)*x**3
    coefficients=[Q(2292,625),Q(41292,625),Q(10987,625),Q(39713,625),Q(12938273,2500),Q(28405139,1250),Q(54190307,1250),Q(26139962,625),Q(174528933,10000),Q(131199,1250),Q(58443,2000),Q(836124,625)]
    assert all(v>0 for v in coefficients)
    assert sp.expand(P*D*D-U-sum(sp.Rational(v)*x**j*(1-x)**(11-j) for j,v in enumerate(coefficients)))==0
    assert sp.integrate(P,(x,0,1))==sp.Rational(132541,60000)<sp.Rational(9,4)
    meanbox=Q(19,24)*Q(1001,8000)+Q(1,20000)+Q(26,25)*Q(5,24)*2
    Abox=Q(11,20)*Q(26,25)/Q(993,1000)**2*Q(173,264)
    assert meanbox<Q(11,20) and Abox==Q(1124500,2958147)<Q(2,5)
    assert Q(11,20)*Q(2,5)**2/Q(7,5)**3==Q(11,343)
    assert Q(9,5)*Q(11,20)-Q(11,20)**2==Q(11,16) and Q(9,5)-2*Q(11,20)>0
    t=sp.Symbol('t',nonnegative=True);birth=sp.Rational(1,10);switch=sp.Rational(1,100000);cs=sp.Rational(c)
    qn=lambda n:sp.exp(-(birth+switch)*t)*(cs*(1-sp.exp(-(birth+switch)*t)))**(n-1)
    for n in range(1,5):
        rhs=-(birth+switch)*n*qn(n)+(birth*(n-1)*qn(n-1) if n>1 else 0)
        assert sp.simplify(sp.diff(qn(n),t)-rhs)==0
        assert sp.simplify(qn(n).subs(t,sp.log(2)/(birth+switch))-sp.Rational(1,2)*(cs/2)**(n-1))==0
    S,R,bs,br,ds,dr,alpha,beta=sp.symbols('S R bs br ds dr alpha beta');N=S+R
    def generator(f):
        return sum(rate*(f.subs({S:S+dS,R:R+dR},simultaneous=True)-f) for dS,dR,rate in [(1,0,bs*S),(0,1,br*R),(-1,0,ds*S),(0,-1,dr*R),(-1,1,alpha*S),(1,-1,beta*R)])
    growth=(bs-ds)*S+(br-dr)*R;turnover=(bs+ds)*S+(br+dr)*R
    assert sp.expand(generator(N)-growth)==sp.expand(generator(N*N)-2*N*growth-turnover)==0
    return dict(actual_tail_at_least_three=L,actual_coverage_at_most_two_upper=1-L,geometric_upper_at_two=1-geo,
        actual_coverage_at_most_three_lower=H,scalar_coverage_at_most_two_lower=Q(332,343),
        actual_variance_lower=variance,scalar_variance_upper=Q(11,16),variance_ratio_upper=Q(55,64),
        mean_box=Q(11,20),birth_parameter_box=Q(2,5),polynomial=str(P),bernstein_coefficients=coefficients,
        markov_only_endpoint=10,minimal_one_sided_95_endpoint=3,
        scope='Fixed published positive-switching witness only. Exact rational/polynomial replay plus imported analytic process arguments; no Lean replay.')


def scenario(rates,p,duration,founders=(1,),cap=TERMINAL_COEFFICIENT_CAP,plating=1.):
    moment=MomentClosure(rates,p).solve(duration);types=BackwardCountLaw(rates).solve(duration,cap)
    one=(1-p)*types[0]+p*types[1];scalar=ScalarTerminalLaw(moment['mean'],moment['A']).coefficients(cap);rows=[]
    for count in founders:
        prep=FounderPreparation(int(count),plating,p);actual=prep.compose(one,cap);approx=prep.compose(scalar,cap);row=PredictionAssessment().assess(actual,approx)
        mean,vn=prep.moments(moment['mean'],moment['variance_actual']);_,vz=prep.moments(moment['mean'],moment['variance_scalar'])
        row.update(founders=int(count),plating=plating,mean=mean,variance_actual=vn,variance_scalar=vz,scalar_over_actual=vz/vn if vn>0 else None)
        rows.append(row)
    return moment,one,scalar,rows


def main():
    ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('--output',type=Path,default=Path('outputs'));ap.add_argument('--simulate',action='store_true');args=ap.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
    rates=Rates();moment,actual,scalar,founders=scenario(rates,RESISTANT_FOUNDER_PROBABILITY,HORIZON_DAYS,FOUNDER_COUNTS)
    exact=exact_witness();trace=moment.pop('trace')
    forward=ForwardOverflowLaw(rates,FORWARD_POPULATION_CAP).solve(HORIZON_DAYS,FounderPreparation(resistant=RESISTANT_FOUNDER_PROBABILITY))
    comparison=float(abs(forward['pmf'][:3].sum()-actual[:3].sum()));forward.pop('pmf')
    assert comparison<forward['overflow']+1e-9 and moment['defect_identity_error']<1e-8 and moment['mean_identity_error']<1e-8
    grid=[]
    for a in ALPHA_GRID:
        for b in BETA_GRID:
            m,_,_,f=scenario(replace(rates,switching=(a,b)),RESISTANT_FOUNDER_PROBABILITY,HORIZON_DAYS,cap=48)
            grid.append(dict(alpha=a,beta=b,interval=f[0]['scalar_interval'],coverage=f[0]['actual_coverage'],founder_limit=m['founder_limit']))
    extensions=[]
    for a,b,e in [(.001,.00001,.001),(.001,.00001,.01),(.1,.1,.01)]:
        r=Rates((e,.1),(.3,e),(a,b));m,_,_,f=scenario(r,5/24,HORIZON_DAYS,(1,100));m.pop('trace');extensions.append(dict(rates=asdict(r),moments=m,founders=f))
    # Source-linked control as specified in the paper, with different plating/type preparation.
    p_ref=(.57-math.sqrt(.57**2-4*.05*.02))/(2*.05);plating=817*614/(math.pi*4500**2)
    ref_rates=Rates((0,.1),(.3,0),(1.,.01));refm,_,_,reff=scenario(ref_rates,p_ref,7.,(1000,),plating=plating);refm.pop('trace')
    refforward=ForwardOverflowLaw(ref_rates,60).solve(7.,FounderPreparation(1000,plating,p_ref));interval=reff[0]['scalar_interval'];forward_cov=float(refforward['pmf'][interval[0]:interval[1]+1].sum());refforward.pop('pmf')
    control=dict(source_commit='401c054dcf588262edfa1263767ffa7ce366533b',scope='Pinned Figure 8 two-type preparation only, not the continuous-phenotype inference.',
        resistant=p_ref,plating=plating,result=reff[0],forward_coverage=forward_cov,forward=refforward,
        actual_over_scalar_variance=reff[0]['variance_actual']/reff[0]['variance_scalar'])
    result=dict(manuscript_sha256=MANUSCRIPT_SHA256,configured_rates=asdict(rates),configured_moments=moment,exact_paper_witness=exact,
        configured_one_founder=dict(coverage_0_2=float(sum(actual[:3])),coverage_0_3=float(sum(actual[:4])),scalar_coverage_0_2=float(sum(scalar[:3]))),
        founders=founders,forward=forward,forward_backward_low_count_difference=comparison,grid=grid,positive_rate_extensions=extensions,source_linked_control=control,
        scope='Exact published-witness inequalities are separate from numerical quantiles for edited parameters. Overflow and solver residuals are diagnostics, not validated roundoff. Independent founders and constant rates only.')
    if args.simulate:result['optional_single_realization']=EventSimulator(rates).run(FounderPreparation(10,1,RESISTANT_FOUNDER_PROBABILITY),HORIZON_DAYS)
    def dump(name,data):(out/name).write_text(json.dumps(data,indent=2,default=str)+'\n')
    def table(name,header,rows):
        with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    dump('results.json',result);table('count_laws.csv',['count','actual','scalar'],[(i,a,b) for i,(a,b) in enumerate(zip(actual,scalar))])
    table('moments.csv',['time','mean_S','mean_R','total_mean','scalar_mean','actual_variance','scalar_variance','excess_growth_covariance','integrated_variance_defect'],trace)
    table('founders.csv',['F','scalar_lower','scalar_upper','actual_coverage','scalar_coverage','actual_one_sided_endpoint'],[(r['founders'],*r['scalar_interval'],r['actual_coverage'],r['scalar_own_coverage'],r['actual_one_sided_endpoint']) for r in founders])
    table('switching_grid.csv',['alpha','beta','lower','upper','coverage','founder_limit'],[(r['alpha'],r['beta'],*r['interval'],r['coverage'],r['founder_limit']) for r in grid])
    lines=['Exact witness arithmetic: endpoint 2 fails 95%; endpoint 3 is the smallest successful one-sided endpoint.',
        f'Configured coverage [0,2]: actual {sum(actual[:3]):.10f}, scalar {sum(scalar[:3]):.10f}; actual [0,3] {sum(actual[:4]):.10f}.',
        f'Configured mean {moment["mean"]:.10f}, actual/scalar variances {moment["variance_actual"]:.10f}/{moment["variance_scalar"]:.10f}.',
        f'Variance-defect identity residual {moment["defect_identity_error"]:.3g}; forward/backward gap {comparison:.3g}, retained overflow {forward["overflow"]:.3g}.',
        f'Independent-founder limit (numerical) {moment["founder_limit"]:.9f}; source-linked control coverage {reff[0]["actual_coverage"]:.9f}.',
        'Numerical repairs for new preparations are not certified thresholds. The exact witness is synthetic, not a fitted assay.']
    (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),layout='constrained');data=np.array(trace)
    for col,label in [(1,'Sensitive'),(2,'Resistant'),(3,'Total = scalar mean')]:axs[0].plot(data[:,0],data[:,col],label=label)
    axs[0].set(xlabel='Time (days)',ylabel='Expected count',title='Mean population count under both models');axs[0].legend(fontsize=8)
    n=np.arange(8);axs[1].bar(n-.18,actual[:8],.36,label='Inherited-state');axs[1].bar(n+.18,scalar[:8],.36,label='Scalar closure');axs[1].axvline(2.5,color='k',ls=':');axs[1].set(xlabel='Terminal count',ylabel='Probability',title='Terminal count probabilities with and\nwithout inherited state');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(axis='y',alpha=.2)
    fig.savefig(out/'counts.png',dpi=180);fig.savefig(out/'counts.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained');axs[0].plot([r['founders'] for r in founders],[r['actual_coverage'] for r in founders],'o-',label='Actual coverage of scalar interval');axs[0].axhline(TARGET_COVERAGE,color='k',ls=':',label='Target');axs[0].axhline(moment['founder_limit'],color='#bd5a24',ls='--',label='Numerical founder limit');axs[0].set(xscale='log',xlabel='Independent founders',ylabel='Actual coverage',title='Prediction-interval coverage versus founder\ncount');axs[0].legend(fontsize=7)
    axs[1].plot(data[:,0],data[:,5]-data[:,6],label='Actual variance − scalar variance');axs[1].plot(data[:,0],data[:,8],'--',label='Integrated covariance identity');axs[1].set(xlabel='Time (days)',ylabel='Variance defect',title='Variance difference and inherited-state\ncovariance');axs[1].legend(fontsize=7)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'coverage.png',dpi=180);fig.savefig(out/'coverage.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,source_sha256=digest(Path(__file__)),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()