Example code
In a growing cell population, inherited sensitive or resistant states affect the distribution of future cell counts, even when their average is predicted correctly. This example compares a two-type branching process, where descendants inherit the parent's current state, with a single-population model built from the exact expected type proportions. The scalar model matches both the mean count and the expected birth and death fluxes, but loses dependence between family size and inherited state.
Each founder receives one sensitive or resistant type draw. Cells divide, die or switch independently at fixed per-cell rates; division preserves the parent's type. The paper's synthetic witness makes switching rare enough that one founder's state strongly influences its family's size.
For one founder, the scalar closure assigns about 97.68% probability to counts 0–2. The inherited-state model gives only 94.76%. Fresh rational and polynomial checks recover the paper's analytic result: endpoint 2 fails 95% coverage, while endpoint 3 is the smallest successful one-sided upper endpoint for this preparation. That correction is instance-specific, not a general add-one rule.
The runnable count model uses equations for the probabilities of low terminal counts, obtained from generating functions. A degree cutoff does not stop the population: histories that grow beyond the reported count and later shrink remain represented. An independent joint forward calculation retains overflow and checks the result. Separate components handle rates, founder preparation, exact moments, scalar laws, count probabilities, prediction assessment and optional event simulation.
With many independent founders, relative fluctuations shrink but the variance ratio remains wrong. At 100 founders the scalar interval is and actual coverage is about 85.11%; the numerical limiting coverage is about 84.88%. Finite-founder coverage moves in jumps and can exceed 95%, so the asymptotic statement is not a bound on every finite case.
The example also integrates the covariance identity explaining the variance gap, sweeps switching rates, includes positive birth and death rates in both states, and retains a source-linked fast-switching preparation where the closure covers correctly. It therefore demonstrates a failure mechanism without claiming that inheritance always inflates variance or that scalar closures always fail.
Editable inputs support new rates, horizons, founder numbers and binomial plating. New thresholds computed by ordinary ODE solvers are numerical recommendations; only the fixed paper witness receives the stated analytic certificate. These are future-count prediction regions under known model parameters, not fitted experimental results or confidence intervals, and Lean is not rerun.
Python source
"""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@[email protected]+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()
Run output
Exact witness arithmetic: endpoint 2 fails 95%; endpoint 3 is the smallest successful one-sided endpoint. Configured coverage [0,2]: actual 0.9475977178, scalar 0.9767606151; actual [0,3] 0.9738271218. Configured mean 0.5186169458, actual/scalar variances 1.0875939888/0.5832556124. Variance-defect identity residual 1.24e-14; forward/backward gap 3.89e-15, retained overflow 7.38e-16. Independent-founder limit (numerical) 0.848799756; source-linked control coverage 0.953216498. Numerical repairs for new preparations are not certified thresholds. The exact witness is synthetic, not a fitted assay.