Example code
Compartments in two inherited chemical states share a finite pool of growth precursor. The high state contains about three times as much growth-linked species per unit size as the low state. A growth event consumes one and one precursor, adds one unit of size, and triggers division when the compartment reaches twice its newborn size. Resident feeds and removals remain maintained throughout the batch.
The quantity being selected is the number of high-state cells relative to low-state cells. Faster growth alone does not settle that comparison: a lineage can hold more material while waiting for its next division. The example implements the literal count reactions, shared resource, complementary daughters and chemical stopping checks, alongside the paper's explicit probability bounds.
The first figure evaluates the proved expressions. At the manuscript setting , two founders, growth coupling , logarithmic odds gain and deadline , the retained raw chemical sum gives a success lower bound of about 0.996889. The simpler error bound exceeds one there and gives no success guarantee. These enormous sufficient copy numbers reflect conservative certificates; they are not estimates of practical operating requirements or actual failure rates.
The second figure is a separate deterministic diagnostic that holds each state's concentration fixed. Size odds rise smoothly while cell-count odds change at divisions. It explains the factor-two allowance for unequal progress toward division when converting relative lineage sizes into relative cell counts. It does not simulate molecular recovery or establish chemical fidelity. For exactly two founders, the paper's successful event implies a high-state frequency of at least by integer cell-count constraints.
The downloadable model separates reaction channels, regions defined by the chemical energy function, complementary partition, population transitions, ancestral observables and error budgets. Editable batch inputs are at the top. Exact rational state updates preserve precursor-plus-size accounting; interval energy checks retain precision at large molecule counts. If the final precursor-consuming event divides a cell, both daughters appear in the endpoint record, and a simultaneous chemical flag takes priority over success.
The default run produces bound tables and figures. An optional direct count simulator reports unfinished paths when its event or numerical limit is reached; it cannot feasibly simulate a full batch at the certified scale. Seven scientific test groups check the count semantics, endpoint ordering, source data, analytic expressions and scope limits. The example does not rerun the Lean proof, and the single-batch guarantee does not automatically extend to repeated transfers.
Python source
"""Finite-batch enrichment under one shared precursor pool.
Default: evaluate the manuscript bounds and inspect a literal count-model kernel.
Optional --simulate runs a budget-limited direct jump simulation, not a proof.
"""
from fractions import Fraction as Q
# EDITABLE INPUTS: manuscript units, not laboratory-calibrated concentrations.
NEWBORN_SIZE = 111 * 10**20
HIGH_FOUNDERS = 1
LOW_FOUNDERS = 1
GROWTH_COUPLING = Q(1,10**11)
LOG_ODDS_GAIN = Q(1,10)
SCALED_DEADLINE = Q(8) # gamma*T; physical T is this divided by gamma
TARGET_ERROR = Q(1,100)
COPY_NUMBER_SWEEP = (14*10**19,111*10**20,308*10**20)
RANDOM_SEED = 2101
EVENT_BUDGET = 1000
MANUSCRIPT_SHA256 = 'b93345aac9b134b314f9e6115f46585fd3e5fc60c2ece42f5da6eacab32199ed'
import argparse
import csv
import hashlib
import json
import math
import platform
from dataclasses import dataclass, replace
from pathlib import Path
import mpmath as mp
import numpy as np
import sympy as sp
A = Q(1,512000000)
B = 16*A
ALPHA = Q(1,10**12)
LAMBDA = Q(1,10**15)
DELTA = Q(1,10**6)
C = Q(1,1024*10**18)
THETA = Q(1,1000)
N_MIN = 14*10**19
GAMMA_MAX = Q(1,10**11)
SPECIES = ('A','B','z','H')
MATRICES = {
'L': ((1115801,-613708,2346047,3444885),(-613708,1111499,-2339313,-3434384),
(2346047,-2339313,6767757,10178806),(3444885,-3434384,10178806,15517433)),
'H': ((846084,-1159490,2352853,3510058),(-1159490,4333988,-6781643,-10233028),
(2352853,-6781643,11398763,17222276),(3510058,-10233028,17222276,26082111)),
}
@dataclass(frozen=True)
class Interval:
lo: Q
hi: Q | None = None
def __post_init__(self):
object.__setattr__(self,'lo',Q(self.lo))
object.__setattr__(self,'hi',Q(self.lo if self.hi is None else self.hi))
if self.lo>self.hi:raise ValueError('Inverted interval')
@staticmethod
def cast(v):return v if isinstance(v,Interval) else Interval(v)
def __add__(self,v):v=self.cast(v);return Interval(self.lo+v.lo,self.hi+v.hi)
__radd__=__add__
def __neg__(self):return Interval(-self.hi,-self.lo)
def __sub__(self,v):return self+-self.cast(v)
def __rsub__(self,v):return self.cast(v)+-self
def __mul__(self,v):
v=self.cast(v);values=[x*y for x in (self.lo,self.hi) for y in (v.lo,v.hi)]
return Interval(min(values),max(values))
__rmul__=__mul__
def __truediv__(self,v):
v=self.cast(v)
if v.lo<=0<=v.hi:raise ZeroDivisionError('Interval denominator includes zero')
return self*Interval(1/v.hi,1/v.lo)
def __rtruediv__(self,v):return self.cast(v)/self
def json(self):return [str(self.lo),str(self.hi)]
def reconstruction(z):
b=60/(z+2);k=(20004*z*z-159984*z)/20001
return (z*b+k,b,z,(16*z+2*z*z)/Q(20001,10000))
def residual(z):
aa,bb,_,_=reconstruction(z)
return Q(1,100000)*aa*aa+aa+Q(99999,100000)*bb-33
class ChemicalRegions:
"""Rational energy enclosures around tightly isolated stationary roots."""
def __init__(self,refinements=150):
self.roots={};self.centers={}
for tag,left,right in [('L',Q('0.99579401232'),Q('0.99579401233')),
('H',Q('2.97636724376'),Q('2.97636724377'))]:
fl=residual(left);assert fl*residual(right)<0
for _ in range(refinements):
mid=(left+right)/2;fm=residual(mid)
if fm==0:left=right=mid;break
if fl*fm<0:right=mid
else:left=mid;fl=fm
self.roots[tag]=Interval(left,right)
self.centers[tag]=reconstruction(self.roots[tag])
def energy(self,cell):
delta=[Q(n,cell.size)-center for n,center in zip(cell.counts,self.centers[cell.tag])]
value=sum((Q(MATRICES[cell.tag][i][j],10**6)*delta[i]*delta[j]
for i in range(4) for j in range(4)),Interval(0))
return Interval(max(Q(0),value.lo),max(Q(0),value.hi))
def admits(self,cell,threshold,closed=False):
energy=self.energy(cell)
if (energy.hi<=threshold if closed else energy.hi<threshold):return True
if (energy.lo>threshold if closed else energy.lo>=threshold):return False
raise ArithmeticError('Energy comparison unresolved; refine the stationary-root brackets.')
def newborn(self,tag,N):
centers=self.centers[tag]
def nearest(v):return (2*v.numerator+v.denominator)//(2*v.denominator)
cell=Cell(tag,N,tuple(nearest(N*(v.lo+v.hi)/2) for v in centers))
if not self.admits(cell,4*A,closed=True):
raise ValueError('Rounded preparation is outside the closed newborn region at this N.')
return cell
@dataclass(frozen=True)
class Cell:
tag: str
size: int
counts: tuple[int,int,int,int]
def __post_init__(self):
if self.tag not in ('H','L') or type(self.size) is not int or self.size<1:
raise ValueError('A cell needs a valid tag and positive integer size.')
if len(self.counts)!=4 or any(type(n) is not int or n<0 for n in self.counts):
raise ValueError('Resident counts must be four nonnegative Python integers.')
@property
def readout(self):
value=10*self.counts[2]-self.counts[1]
return 'H' if value>0 else 'L' if value<0 else 'ambiguous'
@dataclass(frozen=True)
class Channel:
name: str
inputs: tuple[int,...]
outputs: tuple[int,...]
coefficient: Q
def rate(self,cell):
value=self.coefficient*Q(cell.size)**(1-sum(self.inputs))
for count,order in zip(cell.counts,self.inputs):
if count<order:return Q(0)
for offset in range(order):value*=count-offset
return value
def apply(self,cell):
if self.rate(cell)<=0:raise ValueError('A zero-propensity channel cannot fire.')
return replace(cell,counts=tuple(n-out+new for n,out,new in zip(cell.counts,self.inputs,self.outputs)))
class ResidentChemistry:
def __init__(self):
self.channels=(
Channel('A to B+z',(1,0,0,0),(0,1,1,0),Q(1)),
Channel('B+z to A',(0,1,1,0),(1,0,0,0),Q(1)),
Channel('z to H',(0,0,1,0),(0,0,0,1),Q(16)),
Channel('H to z',(0,0,0,1),(0,0,1,0),Q(1)),
Channel('H to 2z',(0,0,0,1),(0,0,2,0),Q(1)),
Channel('2z to H',(0,0,2,0),(0,0,0,1),Q(2)),
Channel('feed A',(0,0,0,0),(1,0,0,0),Q(6)),
Channel('remove A',(1,0,0,0),(0,0,0,0),Q(1)),
Channel('feed B',(0,0,0,0),(0,1,0,0),Q(27)),
Channel('remove B',(0,1,0,0),(0,0,0,0),Q(1)),
Channel('B to 2A',(0,1,0,0),(2,0,0,0),Q(1,100000)),
Channel('2A to B',(2,0,0,0),(0,1,0,0),Q(1,100000)),
Channel('remove H',(0,0,0,1),(0,0,0,0),Q(1,10000)),
)
def generator(self,cell,observable):
base=observable(cell)
return sum((ch.rate(cell)*(observable(ch.apply(cell))-base)
for ch in self.channels if ch.rate(cell)),Q(0))
@dataclass(frozen=True)
class BatchParameters:
N: int = NEWBORN_SIZE
h: int = HIGH_FOUNDERS
low: int = LOW_FOUNDERS
gamma: Q = GROWTH_COUPLING
gain: Q = LOG_ODDS_GAIN
scaled_deadline: Q = SCALED_DEADLINE
def __post_init__(self):
if any(type(v) is not int or v<1 for v in (self.N,self.h,self.low)):
raise ValueError('N and both founder numbers must be positive integers.')
for key in ('gamma','gain','scaled_deadline'):
object.__setattr__(self,key,Q(str(getattr(self,key))))
if self.gamma<=0 or self.scaled_deadline<0:raise ValueError('Invalid growth rate or deadline.')
@property
def M(self):return self.h+self.low
@property
def horizon(self):return self.scaled_deadline/self.gamma
@property
def theorem_scope(self):return self.N>=N_MIN and self.gamma<=GAMMA_MAX
@dataclass(frozen=True)
class Population:
precursor: int
cells: tuple[Cell,...]
divisions: int = 0
status: str = 'active'
@property
def total_size(self):return sum(c.size for c in self.cells)
def tag_size(self,tag):return sum(c.size for c in self.cells if c.tag==tag)
def tag_count(self,tag):return sum(c.tag==tag for c in self.cells)
class ComplementaryPartition:
def draw(self,counts,rng):
if max(counts)>2**53:
raise ArithmeticError('Binomial sampling exceeds the configured exact-integer input range.')
first=tuple(int(rng.binomial(n,.5)) for n in counts)
return first,tuple(n-d for n,d in zip(counts,first))
@staticmethod
def allocation_probability(counts,first):
if max(counts)>20000:raise ValueError('Exact allocation probability size guard.')
if any(d<0 or d>n for n,d in zip(counts,first)):return Q(0)
return Q(math.prod(math.comb(n,d) for n,d in zip(counts,first)),2**sum(counts))
class PopulationKernel:
"""Physical event outcomes retained even on a simultaneous terminal flag."""
def __init__(self,parameters=BatchParameters(),regions=None,partition=None):
self.p=parameters;self.regions=ChemicalRegions() if regions is None else regions
self.partition=ComplementaryPartition() if partition is None else partition
self.chemistry=ResidentChemistry()
def initial(self):
cells=tuple(self.regions.newborn(tag,self.p.N) for tag,count in [('H',self.p.h),('L',self.p.low)] for _ in range(count))
return Population(4*self.p.N*self.p.M,cells)
def growth_rate(self,pop,cell):
return self.p.gamma*Q(pop.precursor,4*self.p.N*self.p.M)*cell.counts[2] if pop.precursor else Q(0)
def events(self,pop):
if pop.status!='active':return []
events=[]
for i,cell in enumerate(pop.cells):
for k,ch in enumerate(self.chemistry.channels):
rate=ch.rate(cell)
if rate:events.append((i,k,rate))
rate=self.growth_rate(pop,cell)
if rate:events.append((i,13,rate))
return events
def step(self,pop,index,channel,rng,allocation=None):
if pop.status!='active':raise ValueError('Terminal records are absorbing.')
if not 0<=index<len(pop.cells) or not 0<=channel<=13:raise ValueError('Unknown event.')
cell=pop.cells[index];Qnext=pop.precursor;divisions=pop.divisions
if channel<13:
updated=self.chemistry.channels[channel].apply(cell);newcells=(updated,)
status='active' if self.regions.admits(updated,B) else 'outer_failure'
else:
if not self.growth_rate(pop,cell):raise ValueError('Disabled growth event.')
n=list(cell.counts);n[2]-=1
updated=Cell(cell.tag,cell.size+1,tuple(n));Qnext-=1
division=updated.size==2*self.p.N
if updated.size>2*self.p.N:raise ValueError('Cell exceeded the division threshold.')
if division:
first,second=self.partition.draw(updated.counts,rng) if allocation is None else (tuple(allocation),tuple(n-d for n,d in zip(updated.counts,allocation)))
newcells=(Cell(cell.tag,self.p.N,first),Cell(cell.tag,self.p.N,second));divisions+=1
else:newcells=(updated,)
# Outer flags, parent division energy, daughter partition, endpoint.
if not self.regions.admits(updated,B):status='outer_failure'
elif division and not self.regions.admits(updated,2*A,closed=True):status='division_energy_failure'
elif division and any(not self.regions.admits(c,4*A) for c in newcells):status='partition_failure'
elif Qnext==self.p.N*self.p.M:status='nutrient_endpoint'
else:status='active'
outcome=Population(Qnext,pop.cells[:index]+newcells+pop.cells[index+1:],divisions,status)
assert outcome.precursor+outcome.total_size==pop.precursor+pop.total_size
assert len(outcome.cells)-len(pop.cells)==divisions-pop.divisions
assert all(self.p.N<=c.size<2*self.p.N for c in outcome.cells)
return outcome
def simulate(self,budget=EVENT_BUDGET,seed=RANDOM_SEED):
if type(budget) is not int or budget<1:raise ValueError('Use a positive integer event budget.')
rng=np.random.default_rng(seed);pop=self.initial();time=0.;trace=[]
status='event_budget';steps=0
for steps in range(budget):
events=self.events(pop);rates=np.array([float(e[2]) for e in events]);total=rates.sum()
wait=float(rng.exponential(1/total))
if time+wait>float(self.p.horizon):time=float(self.p.horizon);status='deadline';break
event=events[int(rng.choice(len(events),p=rates/total))]
try:following=self.step(pop,event[0],event[1],rng)
except ArithmeticError as exc:status='computational_limit: '+str(exc);break
pop=following;time+=wait
if steps%100==0 or pop.status!='active':
trace.append([time,pop.precursor,pop.total_size,pop.tag_count('H'),pop.tag_count('L'),pop.status])
if pop.status!='active':status=pop.status;break
determined=status=='deadline' or pop.status!='active'
success=self.success(pop) if determined else None
final=[time,pop.precursor,pop.total_size,pop.tag_count('H'),pop.tag_count('L'),pop.status]
if not trace or trace[-1]!=final:trace.append(final)
return {'time':time,'steps_attempted':steps+1,'status':status,'success':success,
'theorem_parameter_scope':self.p.theorem_scope,'trace':trace}
def success(self,pop):
if pop.status!='nutrient_endpoint':return False
if any(c.readout!=c.tag for c in pop.cells):return False
H=pop.tag_count('H');L=pop.tag_count('L')
if not H or not L:return False
with mp.workdps(80):return bool(mp.log(mp.mpf(H)*self.p.low/(L*self.p.h))>=mpq(self.p.gain))
def mpq(v):
v=Q(v);return mp.mpf(v.numerator)/v.denominator
class ErrorBudget:
"""High-precision evaluation of proved expressions, not empirical failure rates."""
def __init__(self,parameters=BatchParameters()):self.p=parameters
def logs(self):
p=self.p;N=mp.mpf(p.N);M=mp.mpf(p.M);T=mpq(p.horizon);u=mpq(ALPHA*A)*N
logT=mp.log(T) if T else -mp.inf
branches={
'outer_initial':mp.log(7*M)-12*u,
'outer_drift':mp.log(M/84)+logT+mp.log(u)-mp.mpf('15.5')*u,
'division_initial':mp.log(M)+2*u-mpq(LAMBDA)*N,
'division_drift':mp.log(M/84)+logT+mp.log(u)-mp.mpf('1.5')*u,
'partition':mp.log(24*M)-N*mpq(DELTA**2/35),
}
raw=mp.log(sum(mp.exp(v) for v in branches.values()))
selection=-mpq(THETA)*N*(mp.log(2)/5-mpq(p.gain))
deadline=-mpq(THETA)*N*(mpq(Q(9,40)*p.scaled_deadline)-mp.log(4))
simple=mp.log(M*(32+T/42))-mpq(C)*N
sharp=mp.log(M*(32+T*u/42))-3*u/2
double=mp.log(M*(32+T/21))-u
return branches,{'raw':raw,'simplified':simple,'retained':sharp,'double_exponent':double},selection,deadline
def evaluate(self):
with mp.workdps(80):
branches,chem,sel,deadline=self.logs()
totals={name:mp.log(sum(mp.exp(v) for v in (value,sel,deadline))) for name,value in chem.items()}
format=lambda v:mp.nstr(v,35)
return {'theorem_parameter_scope':self.p.theorem_scope,
'branch_log10':{k:format(v/mp.log(10)) for k,v in branches.items()},
'chemical_error':{k:format(mp.exp(v)) for k,v in chem.items()},
'total_error_log10':{k:format(v/mp.log(10)) for k,v in totals.items()},
'success_lower':{k:format(max(mp.mpf(0),-mp.expm1(v))) if self.p.theorem_scope else None for k,v in totals.items()},
'selection_log10':format(sel/mp.log(10)),'deadline_log10':format(deadline/mp.log(10)),
'selection_margin':format(mp.log(2)/5-mpq(self.p.gain)),
'deadline_margin':format(mpq(Q(9,40)*self.p.scaled_deadline)-mp.log(4))}
def sufficient_size(self,error=TARGET_ERROR,double_exponent=False):
p=self.p
if not 0<Q(error)<1 or p.gamma>GAMMA_MAX:raise ValueError('Use error in (0,1) and certified gamma.')
if p.gain!=Q(1,10) or p.scaled_deadline!=8:
raise ValueError('This closed-form sufficient size is for s=0.1 and gamma*T=8.')
with mp.workdps(80):
coefficient=p.M*(32+(8 if double_exponent else 4)/(21*mpq(p.gamma)))
exponent=mpq(C)*(2 if double_exponent else 1)
candidates=[mp.mpf(N_MIN),mp.log(3*coefficient/mpq(error))/exponent,
mp.mpf(500000)/19*mp.log(3/mpq(error)),2500*mp.log(3/mpq(error))]
# One extra integer avoids understating a numerically evaluated ceiling.
return int(mp.ceil(max(candidates)))+1
class AncestralObservables:
def __init__(self,parameters):self.p=parameters
def log_phi(self,pop):
p=self.p
with mp.workdps(80):
H=mp.mpf(pop.tag_size('H'));L=mp.mpf(pop.tag_size('L'));W=H+L
return mpq(THETA)*p.N*(-mp.log(H/(p.N*p.h))+mp.log(L/(p.N*p.low))+mp.mpf(3)/5*mp.log(W/(p.N*p.M)))
def size_to_count_ratio(self,pop):
H=pop.tag_count('H');L=pop.tag_count('L')
if not H or not L:raise ValueError('Both tags must be present.')
return {'size_odds':Q(pop.tag_size('H'),pop.tag_size('L')),'count_odds':Q(H,L),
'phase_factor':Q(pop.tag_size('H')*L,pop.tag_size('L')*H)}
def relative_odds_drift(self,pop,high_z,low_z):
"""L Phi / (beta Phi), useful for checking the source-safe rate sandwich."""
with mp.workdps(80):
H=mp.mpf(pop.tag_size('H'));L=mp.mpf(pop.tag_size('L'));W=H+L
thetaN=mpq(THETA)*self.p.N
vH=thetaN*(-mp.log1p(1/H)+mp.mpf(3)/5*mp.log1p(1/W))
vL=thetaN*(mp.log1p(1/L)+mp.mpf(3)/5*mp.log1p(1/W))
return mpq(high_z)*H*mp.expm1(vH)+mpq(low_z)*L*mp.expm1(vL)
def log_position_exponential(self,cell,regions):
return LAMBDA*(cell.size-2*self.p.N)+ALPHA*self.p.N*regions.energy(cell)
def scalar_certificates():
assert ALPHA*A==2*C
assert LAMBDA>=Q(7,2)*ALPHA*A and DELTA**2/35>=Q(3,2)*ALPHA*A
assert 32*GAMMA_MAX*LAMBDA<=ALPHA*A/1344
assert (4*ALPHA*A-LAMBDA)*N_MIN<=Q(-4479965,32)<-1
assert 42*16*DELTA**2<A
assert 200000000+8000000000*N_MIN*GAMMA_MAX**2<=N_MIN*A/672
assert -Q(3,20000)+Q(6,390625)<0
# Exact positive-definite matrix checks replace numerical eigenvalue checks.
for tag,upper in [('L',24),('H',42)]:
P=sp.Matrix(MATRICES[tag])/10**6
for matrix in (P-sp.eye(4)/200,upper*sp.eye(4)-P):
assert all(matrix[:k,:k].det()>0 for k in range(1,5))
candidates=[Q(h,h+l) for h in range(1,9) for l in range(1,9) if h>l and h+l<=8]
assert min(candidates)==Q(4,7)
return {'exact_scalar_margins':True,'exact_matrix_sandwiches':True,'two_founder_frequency':'4/7',
'scope':'does not rerun full nonlinear finite-count generator or Lean certificates'}
def frozen_composition_diagnostic(p,regions):
"""Deterministic phase illustration, explicitly NOT the resident count model.
Chemistry is clamped at each stationary z. Exposure q satisfies dE/dt=gamma*Q/Omega.
Total ancestral size/N is h exp(zH E)+low exp(zL E). Division conserves it,
and equal daughters give count_i=founders_i*2**floor(z_i E/log(2)).
"""
zH=float((regions.roots['H'].lo+regions.roots['H'].hi)/2)
zL=float((regions.roots['L'].lo+regions.roots['L'].hi)/2)
left=0.;right=10.
for _ in range(70):
mid=(left+right)/2
if p.h*math.exp(zH*mid)+p.low*math.exp(zL*mid)>4*p.M:right=mid
else:left=mid
end=(left+right)/2;rows=[]
for exposure in np.linspace(0,end,201):
BH=p.h*math.exp(zH*exposure);BL=p.low*math.exp(zL*exposure)
CH=p.h*2**math.floor(zH*exposure/math.log(2));CL=p.low*2**math.floor(zL*exposure/math.log(2))
rows.append([exposure,BH,BL,CH,CL,BH/BL,CH/CL,(5*p.M-BH-BL)/(4*p.M)])
return rows
def main():
parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'))
parser.add_argument('--simulate',action='store_true',help='Attempt a budget-limited literal count path.')
args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
p=BatchParameters();regions=ChemicalRegions();kernel=PopulationKernel(p,regions)
initial=kernel.initial();assert initial.precursor+initial.total_size==5*p.N*p.M
budget=ErrorBudget(p);evaluation=budget.evaluate();checks=scalar_certificates()
table=[{'N':N,**ErrorBudget(replace(p,N=N)).evaluate()} for N in COPY_NUMBER_SWEEP]
phase=frozen_composition_diagnostic(p,regions)
result={'parameters':{'N':p.N,'h':p.h,'low':p.low,'gamma':str(p.gamma),'gain':str(p.gain),'horizon':str(p.horizon)},
'budget':evaluation,'copy_number_table':table,'scalar_checks':checks,
'founders':[{'tag':c.tag,'size':c.size,'counts':c.counts,'energy_interval':regions.energy(c).json(),
'readout':c.readout} for c in initial.cells],
'stationary_root_intervals':{tag:v.json() for tag,v in regions.roots.items()},
'initial_total_event_rate':str(sum(e[2] for e in kernel.events(initial))),
'sufficient_N_original':budget.sufficient_size() if p.gain==Q(1,10) and p.scaled_deadline==8 and p.gamma<=GAMMA_MAX else None,
'sufficient_N_double_exponent':budget.sufficient_size(double_exponent=True) if p.gain==Q(1,10) and p.scaled_deadline==8 and p.gamma<=GAMMA_MAX else None,
'simulation':kernel.simulate() if args.simulate else None,
'diagnostic_scope':'frozen-composition deterministic phase illustration, not a count-model simulation'}
(out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
def csv_file(name,headers,rows):
with (out/name).open('w',newline='') as stream:
w=csv.writer(stream);w.writerow(headers);w.writerows(rows)
csv_file('bound_table.csv',['N','raw_chemical_error','simplified_chemical_error','raw_success_lower','simplified_success_lower'],
[(r['N'],r['chemical_error']['raw'],r['chemical_error']['simplified'],r['success_lower']['raw'],r['success_lower']['simplified']) for r in table])
csv_file('phase_diagnostic.csv',['exposure','high_size_per_N','low_size_per_N','high_count','low_count','size_odds','count_odds','precursor_fraction'],phase)
sweep=[]
for N in np.geomspace(float(N_MIN),4e22,180):
d=ErrorBudget(replace(p,N=int(N))).evaluate()
sweep.append([int(N),d['total_error_log10']['raw'],d['total_error_log10']['simplified'],d['success_lower']['raw'],d['success_lower']['simplified']])
csv_file('copy_number_sweep.csv',['N','raw_total_log10','simplified_total_log10','raw_success_lower','simplified_success_lower'],sweep)
lines=[f'Parameter scope for the theorem: {p.theorem_scope}; N={p.N}; M={p.M}.',
f'Raw-formula success lower bound: {evaluation["success_lower"]["raw"]}.',
f'Simplified-formula success lower bound: {evaluation["success_lower"]["simplified"]}.',
f'Sufficient N for error {TARGET_ERROR} (original): {result["sufficient_N_original"]}.',
f'Sufficient N using double exponent: {result["sufficient_N_double_exponent"]}.',
'Two-founder success implies high-state frequency at least 4/7.',
'The phase diagnostic clamps chemistry; it does not estimate chemical fidelity.',
'The count kernel is available; default outputs evaluate bounds rather than simulating an astronomical batch.']
(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),layout='constrained');data=np.array(sweep,dtype=float)
for col,color,label in [(1,'#287a99','Raw chemical sum + population tails'),(2,'#777777','Simplified chemical envelope + tails')]:
axes[0].plot(data[:,0]/1e22,data[:,col],color=color,label=label)
axes[1].plot(data[:,0]/1e22,data[:,col+2],color=color,label=label)
axes[0].axhline(0,color='black',lw=.7);axes[0].set(ylabel='log10 of total error bound',ylim=(-35,12))
axes[1].set(ylabel='Success probability lower bound',ylim=(-.02,1.02));axes[1].legend(fontsize=7)
for ax in axes:ax.set_xlabel('Newborn size N / 10^22');ax.grid(alpha=.2)
fig.suptitle('Success bounds versus newborn molecule\ncount')
fig.savefig(out/'bounds.png',dpi=180);fig.savefig(out/'bounds.svg');plt.close(fig)
fig,axes=plt.subplots(1,2,figsize=(10,4),layout='constrained');data=np.array(phase)
axes[0].plot(data[:,0],data[:,5],label='Ancestral size odds');axes[0].step(data[:,0],data[:,6],where='post',label='Cell-count odds')
axes[0].set(ylabel='High / low',xlabel='Integrated precursor exposure');axes[0].legend(fontsize=8)
axes[1].plot(data[:,0],data[:,7]);axes[1].axhline(.25,color='black',ls='--',label='Batch endpoint')
axes[1].set(ylabel='Remaining precursor / initial precursor',xlabel='Integrated precursor exposure');axes[1].legend(fontsize=8)
for ax in axes:ax.grid(alpha=.2)
fig.suptitle('Lineage size and cell count with fixed\ncompositions')
fig.savefig(out/'division_phase.png',dpi=180);fig.savefig(out/'division_phase.svg');plt.close(fig)
digest=lambda p:hashlib.sha256(p.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':{p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}},indent=2)+'\n')
if __name__=='__main__':main()
Run output
Parameter scope for the theorem: True; N=11100000000000000000000; M=2. Raw-formula success lower bound: 0.99688942456136952951624618556887646. Simplified-formula success lower bound: 0.0. Sufficient N for error 1/100 (original): 30788748905704567310306. Sufficient N using double exponent: 15749265808868895654117. Two-founder success implies high-state frequency at least 4/7. The phase diagnostic clamps chemistry; it does not estimate chemical fidelity. The count kernel is available; default outputs evaluate bounds rather than simulating an astronomical batch.