Example code
Serial transfer tests whether a growth advantage persists when only some compartments seed the next batch. Here an inherited chemical state changes growth on a shared precursor, while every intact compartment has the same chance of transfer. The example connects the resident chemistry, complementary division, uniform transfer, precursor-free recovery and refill in reusable components.
The four internal species have two stationary compositions, with growth-linked concentration near 0.996 and 2.976. No separate fitness parameter is assigned to their ancestry labels. Growth consumes shared precursor and internal ; division partitions each molecule between complementary daughters. Transfer selects exactly cells without replacement, and recovery evolves the chemical states of those selected cells at their existing sizes.
The code freshly checks the material completion, compatible rate ratios, stationary reconstruction and exact bounds on the quadratic functions that test whether a chemical state is ready for the next cycle. It enumerates all 70 ways to retain four of eight cells. At the tight sampling tolerance, the exact failure fraction rises from for homogeneous sizes within each type to for mixed sizes. This shows why unbiased transfer still needs an explicit error allowance.
The certificate combines batch, transfer, recovery and service errors on the same two-cycle law. For the smaller paper witness, and , the original rational argument certifies 0.991998 success. Evaluating the full formulas with the distinct first- and second-cycle count floors gives the downward-rounded lower bound 0.993959999. These are existence-witness scales, not laboratory predictions.
On success, the high-state fraction increases from one half to more than 0.6930453. The key conversion is from total compartment size to compartment number: their log odds differ by the ratio of mean cell sizes. Intermediate division phases cancel across cycles, so newborn founders require a single correction for final cell sizes. Faster size growth alone would not establish selection in cell counts.
The package includes an exact molecule-count event model, shared-resource growth, complementary partition, intact-cell sampling, saturating service counters, a separate hard expenditure cutoff, probability-budget sweeps and seven scientific test groups. A bounded optional stochastic run keeps partial populations and reports unfinished stages. The displayed recovery curves are a separate density illustration from constructed admissible endpoint cells; they do not estimate finite-count probabilities.
Only two cycles are certified at the stated eligibility floors. Maintained reservoirs, effective growth/division and ready founders are assumptions; the formal probability and local dissipation results are imported without rerunning Lean. Fresh growth supply, discarded cells, removed medium and resident reservoir exchange are kept in separate accounts. The quotas do not represent autonomous reservoir depletion or total energy cost.
Python source
"""Resident chemistry, intact-cell transfer and two-cycle chemical-state selection."""
# EDITABLE INPUTS: theorem witness and separate, deliberately uncertified small diagnostics.
NEWBORN_SIZE = 65536 * 10**18
RETAINED_CELLS = 4 * 10**9
GROWTH_COUPLING = '1/100000000000'
SERVICE_FAILURE_ALLOWANCE = '1/1000'
DIAGNOSTIC_TRANSFER_TOLERANCE = '1/50'
DIAGNOSTIC_SIZE = 10
DIAGNOSTIC_CELLS = 4
DIAGNOSTIC_GROWTH = '2' # Accelerated exploration; outside theorem scope.
DIAGNOSTIC_RECOVERY = '.05' # The theorem instead requires 5376.
EVENT_BUDGET = 100000
RANDOM_SEED = 44092026
from dataclasses import dataclass,replace,asdict
from fractions import Fraction as Q
from itertools import combinations
from pathlib import Path
import argparse,csv,hashlib,json,math,platform
import numpy as np
import sympy as sp
from mpmath import mp,iv
from scipy.integrate import solve_ivp
from resident import Cell,ResidentChemistry,ChemicalRegions,Interval,MATRICES,A,reconstruction,residual
mp.dps=80;iv.dps=80
MANUSCRIPT_SHA256='d4d09105a0cfd39230861a7911081421a89f536afc2caf28d9297d00cab282ca'
ALPHA=Q(1,10**12);MIN_SIZE=14*10**19
SPECIES=('A','B','z','H');RESERVOIRS=('F','G','RA','RB','WH')
# Gross exchange for each directed resident label, including supplied and collected molecules.
EXCHANGE=(None,None,0,0,None,None,2,2,3,3,1,1,4)
def number(ctx,x):
x=Q(str(x));return ctx.mpf(x.numerator)/x.denominator
def dot(a,b):return sum(x*y for x,y in zip(a,b))
class ChemicalModel(ResidentChemistry):
def __init__(self):
super().__init__();self.jumps=np.array([np.subtract(c.outputs,c.inputs) for c in self.channels]);self.coefficients=np.array([float(c.coefficient) for c in self.channels])
def rates(self,cell):
a,b,z,h=cell.counts;m=cell.size
return np.array([a,b*z/m,16*z,h,h,2*z*(z-1)/m,6*m,a,27*m,b,b/100000,a*(a-1)/(100000*m),h/10000],float)
def density(self,x):
a,b,z,h=x
return np.array([6-2*a+b*z+2e-5*(b-a*a),27+a-(1+z)*b-1e-5*(b-a*a),a-b*z-16*z-4*z*z+3*h,16*z+2*z*z-2.0001*h])
def certificates(self):
# Full species A,B,z,H,F,G,RA,RB,WH, with maintained reservoirs projected out.
left=((0,),(2,4),(3,),(6,),(7,),(1,5));right=((1,2),(3,),(2,2),(0,),(1,),(0,0))
element=(2,1,1,2,1,3,2,1,2);boltzmann=(Q(1),Q(1),Q(1),Q(2),Q(1,8),Q(1),Q(1,6),Q(1,27),Q(1))
for j,(l,r) in enumerate(zip(left,right)):
if sum(element[i] for i in l)!=sum(element[i] for i in r):raise ArithmeticError('Full material balance failed.')
ratio=math.prod(boltzmann[i] for i in r)/math.prod(boltzmann[i] for i in l)
if ratio!=self.channels[2*j].coefficient/self.channels[2*j+1].coefficient:raise ArithmeticError('Common chemical potentials failed.')
for channel,reactants,products in [(self.channels[2*j],l,r),(self.channels[2*j+1],r,l)]:
if channel.inputs!=tuple(reactants.count(i) for i in range(4)) or channel.outputs!=tuple(products.count(i) for i in range(4)):raise ArithmeticError('Reservoir projection failed.')
if element[3]!=element[8]:raise ArithmeticError('Sink material imbalance.')
matrices={}
for tag,P in MATRICES.items():
P=sp.Matrix(P)/10**6
if any((P-sp.eye(4)/200)[:k,:k].det()<=0 or (42*sp.eye(4)-P)[:k,:k].det()<=0 for k in range(1,5)):raise ArithmeticError('Energy sandwich failed.')
matrices[tag]=P.tolist()
z=sp.Symbol('z');coords=reconstruction(z);a,b,zz,h=coords;eps=sp.Rational(1,100000)
f=(6-2*a+b*zz+2*eps*(b-a*a),27+a-(1+zz)*b-eps*(b-a*a),a-b*zz-16*zz-4*zz*zz+3*h,16*zz+2*zz*zz-sp.Rational(20001,10000)*h)
r=residual(z)
if any(sp.simplify(x)!=0 for x in (f[0]+2*r,f[1]-r,f[2],f[3])):raise ArithmeticError('Stationary reconstruction failed.')
return dict(element=element,boltzmann_weights=list(map(str,boltzmann)),energy_matrices=matrices,stationary_curve_identity=True,scope='Exact chemical, source projection and energy-matrix checks. The local stochastic exponential-generator inequality is imported from the paper.')
def recovery_density(self,cell,duration=5376,method='Radau'):
initial=np.array(cell.counts,float)/cell.size
sol=solve_ivp(lambda t,x:self.density(x),(0,duration),initial,method=method,rtol=2e-12,atol=1e-14,dense_output=True)
if not sol.success:raise ArithmeticError(sol.message)
return sol
@dataclass(frozen=True)
class Population:
cells:tuple
precursor:int
omega:int
endpoint:int
divisions:int=0
@property
def size(self):return sum(c.size for c in self.cells)
class ComplementaryDivision:
def split(self,counts,rng):
if max(counts)>=2**53:raise ValueError('Numerical division requires molecule counts below 2^53.')
first=tuple(int(rng.binomial(n,.5)) for n in counts)
return first,tuple(n-a for n,a in zip(counts,first))
class PopulationSource:
"""Physical source only: no type-dependent rates, reset, or analytical failure censoring."""
def __init__(self,N,gamma,chemistry=None,division=None):
self.N=N;self.gamma=Q(str(gamma));self.chemistry=chemistry or ChemicalModel();self.division=division or ComplementaryDivision()
if type(N) is not int or N<1 or self.gamma<=0:raise ValueError('Positive integer N and growth coefficient required.')
def refill(self,cells):
cells=tuple(cells)
if not cells or any(not self.N<=c.size<2*self.N for c in cells):raise ValueError('Cells must lie in the live-size interval.')
W=sum(c.size for c in cells);return Population(cells,4*W,4*W,W)
def growth_rate(self,pop,cell):return self.gamma*Q(pop.precursor,pop.omega)*cell.counts[2] if pop.precursor else Q(0)
def step(self,pop,index,label,rng,allocation=None):
if type(index) is not int or not 0<=index<len(pop.cells) or type(label) is not int or not 0<=label<=13:raise ValueError('Unknown cell or event label.')
cell=pop.cells[index];Qnext=pop.precursor;divisions=pop.divisions
if label<13:children=(self.chemistry.channels[label].apply(cell),)
elif label==13:
if self.growth_rate(pop,cell)<=0:raise ValueError('Disabled growth event.')
n=list(cell.counts);n[2]-=1;grown=Cell(cell.tag,cell.size+1,tuple(n));Qnext-=1
if grown.size==2*self.N:
if allocation is None:first,second=self.division.split(grown.counts,rng)
else:
first=tuple(allocation);second=tuple(n-a for n,a in zip(grown.counts,first))
if len(first)!=4 or min(first+second)<0:raise ValueError('Invalid complementary allocation.')
children=(Cell(cell.tag,self.N,first),Cell(cell.tag,self.N,second));divisions+=1
elif grown.size<2*self.N:children=(grown,)
else:raise ValueError('Division threshold exceeded.')
else:raise ValueError('Unknown event label.')
result=replace(pop,cells=pop.cells[:index]+children+pop.cells[index+1:],precursor=Qnext,divisions=divisions)
if result.size+result.precursor!=pop.size+pop.precursor:raise ArithmeticError('Precursor/size conservation failed.')
return result
def run(self,pop,duration,rng,event_budget,collect_endpoint=True,quota=None,hard_limit=None):
if duration<0 or type(event_budget) is not int or event_budget<0:raise ValueError('Nonnegative duration/event budget required.')
if quota is not None and (type(quota) is not int or quota<1):raise ValueError('Positive integer service quota required.')
if hard_limit is not None and (type(hard_limit) is not int or hard_limit<0):raise ValueError('Nonnegative integer hard budget required.')
time=0.;gross=[0]*5;steps=0;status='event_budget';counter=0
if duration==0:return pop,dict(status='duration',time=0.,events=0,gross_exchange=gross,service_counter=0,saturated=False,complete=True)
for _ in range(event_budget):
if collect_endpoint and pop.precursor==pop.endpoint:status='endpoint';break
arrays=[np.r_[self.chemistry.rates(c),float(self.growth_rate(pop,c))] for c in pop.cells];rates=np.concatenate(arrays);total=float(rates.sum())
proposed=time+float(rng.exponential(1/total))
if proposed>duration:time=float(duration);status='duration';break
if proposed<=time:raise ArithmeticError('Time resolution exhausted; use a smaller diagnostic.')
chosen=min(int(np.searchsorted(np.cumsum(rates),rng.random()*total,side='right')),len(rates)-1);index,label=divmod(chosen,14)
reservoir=EXCHANGE[label] if label<13 else None
if hard_limit is not None and reservoir is not None and gross[reservoir]>=hard_limit:status='hard_service_cutoff';break
pop=self.step(pop,index,label,rng);steps+=1;time=proposed
if reservoir is not None:gross[reservoir]+=1
counter=counter+1 if quota is None else min(quota,counter+1)
if collect_endpoint and pop.precursor==pop.endpoint:status='endpoint';break
# Saturating the counter never changes a physical transition or terminates the source.
return pop,dict(status=status,time=time,events=steps,gross_exchange=gross,service_counter=counter,saturated=quota is not None and counter==quota,complete=status in ('endpoint','duration'))
class UniformTransfer:
def sample(self,cells,M,rng):
cells=tuple(cells)
if type(M) is not int or not 1<=M<=len(cells):raise ValueError('Sample size must fit the endpoint population.')
chosen=set(map(int,rng.choice(len(cells),M,replace=False)))
return tuple(c for i,c in enumerate(cells) if i in chosen),tuple(c for i,c in enumerate(cells) if i not in chosen)
def enumerate(self,cells,M,tolerance):
cells=tuple(cells);eps=Q(tolerance);n=len(cells)
if not 0<eps<1 or not 1<=M<=n or math.comb(n,M)>100000:raise ValueError('Invalid transfer or exact enumeration exceeds 100,000 subsets.')
means={tag:Q(M,n)*sum(c.size for c in cells if c.tag==tag) for tag in ('H','L')};values=[];bad=0
for indices in combinations(range(n),M):
sizes={tag:sum(cells[i].size for i in indices if cells[i].tag==tag) for tag in ('H','L')}
failure=any(abs(sizes[tag]-means[tag])>eps*means[tag] for tag in means);bad+=failure;values.append(sizes)
moments={}
for tag in means:
mean=sum(Q(row[tag]) for row in values)/len(values);variance=sum((row[tag]-mean)**2 for row in values)/len(values)
moments[tag]=dict(mean=str(mean),variance=str(variance))
return dict(subsets=len(values),bad=bad,failure=str(Q(bad,len(values))),means={k:str(v) for k,v in means.items()},moments=moments)
@staticmethod
def weighted_moments(weights,M):
weights=tuple(map(Q,weights));n=len(weights)
if not 1<=M<=n or any(not 0<=w<=1 for w in weights):raise ValueError('Weights in [0,1] and valid sample size required.')
p=Q(M,n);joint=Q(M*(M-1),n*(n-1)) if n>1 else Q(0)
mean=p*sum(weights);variance=p*(1-p)*sum(w*w for w in weights)+2*(joint-p*p)*sum(weights[i]*weights[j] for i in range(n) for j in range(i+1,n))
if not 0<=variance<=mean:raise ArithmeticError('Without-replacement variance bound failed.')
return dict(mean=str(mean),variance=str(variance),pair_covariance=str(joint-p*p) if n>1 else None)
class SerialProtocol:
def __init__(self,source,M,recovery_duration=5376):
if type(M) is not int or M<1 or recovery_duration<0:raise ValueError('Positive retained count and nonnegative recovery duration required.')
self.source=source;self.M=M;self.recovery_duration=recovery_duration;self.transfer=UniformTransfer()
def run(self,cells,cycles,seed,event_budget):
if type(cycles) is not int or cycles<0 or type(event_budget) is not int or event_budget<0:raise ValueError('Nonnegative integer cycle count and event budget required.')
rng=np.random.default_rng(seed);pop=self.source.refill(cells);history=[];used=0
for cycle in range(cycles):
W0=pop.size;pop,batch=self.source.run(pop,float(8/self.source.gamma),rng,event_budget-used);used+=batch['events']
if batch['status']!='endpoint':return dict(complete=False,stage='batch',history=history,partial=batch,population=asdict(pop),events=used)
if len(pop.cells)<self.M:return dict(complete=False,stage='insufficient_cells',history=history,events=used)
selected,discarded=self.transfer.sample(pop.cells,self.M,rng);before=selected;residual=pop.precursor
recovering=replace(pop,cells=selected,precursor=0)
recovered,recovery=self.source.run(recovering,self.recovery_duration,rng,event_budget-used,collect_endpoint=False);used+=recovery['events']
record=dict(cycle=cycle,batch=batch,recovery=recovery,start_size=W0,consumed_precursor=3*W0,residual_discarded=residual,discarded_size=sum(c.size for c in discarded),selected_size=sum(c.size for c in selected),selected=[asdict(c) for c in before],recovered=[asdict(c) for c in recovered.cells])
history.append(record)
if recovery['status']!='duration':return dict(complete=False,stage='recovery',history=history,population=asdict(recovered),events=used)
pop=self.source.refill(recovered.cells)
return dict(complete=True,history=history,population=asdict(pop),events=used,scope='Physical diagnostic completion, not a theorem-success classification; analytical energy and odds marks are not applied.')
class TwoCycleCertificate:
def __init__(self,N=NEWBORN_SIZE,M=RETAINED_CELLS,gamma=GROWTH_COUPLING,service=SERVICE_FAILURE_ALLOWANCE):
self.N=N;self.M=M;self.gamma=Q(str(gamma));self.service=Q(str(service));self.u=Q(N)*ALPHA*A;self.T=8/self.gamma
if type(N) is not int or N<MIN_SIZE or type(M) is not int or M<2 or M%2 or not 0<self.gamma<=Q(1,10**11) or not 0<self.service<1:raise ValueError('Requires admitted N, even M>=2, 0<gamma<=1e-11 and service allowance in (0,1).')
self.qR=34000*N;self.qB=280000*N*M;self.JB=math.ceil(self.T*self.qB/self.service)+1;self.JR=math.ceil(M*5376*self.qR/self.service)+1
def logs(self,p,ctx=mp):
p=Q(str(p))
if not 0<p<=1:raise ValueError('Positive count floor <=1 required.')
N=ctx.mpf(self.N);M=ctx.mpf(self.M);u=number(ctx,self.u);T=number(ctx,self.T);ln=ctx.log
return dict(outer_initial=ln(M)-7*u,daughter_initial=ln(14*M)-4*u,outer_drift=ln(M*T*u/42)-15*u/2,division_initial=ln(M)-u,division_drift=ln(M*T*u/42)-3*u/2,partition=ln(56*M)-N/(35*10**12),deadline=-N/2500,size_odds=-19*N/500000,
batch_service=ln(T*self.qB/self.JB),transfer=ln(number(ctx,32/(Q(1,50)**2*p*self.M))),recovery_initial=ln(M)-u,recovery_return=ln(2*M)-u/2,recovery_exit=ln(M)-8*u,recovery_exit_drift=ln(16*M*u)-31*u/2,recovery_service=ln(M*5376*self.qR/self.JR))
def evaluate(self):
terms=[self.logs(p) for p in ('1/2','1/100')];intervals=[self.logs(p,iv) for p in ('1/2','1/100')]
error=sum(mp.exp(x) for d in terms for x in d.values());error_interval=sum(iv.exp(x) for d in intervals for x in d.values());lower=1-error_interval
digits=max(0,int(mp.floor(mp.mpf(lower.a)*10**9)))
while digits and not bool(number(iv,Q(digits,10**9)).b<=lower.a):digits-=1
return dict(N=self.N,M=self.M,u=str(self.u),success_lower_ordinary=mp.nstr(max(0,1-error),30),success_lower_certified=f'{digits/10**9:.9f}',error_interval=str(error_interval),cycle_component_log10=[{k:mp.nstr(v/mp.log(10),25) for k,v in d.items()} for d in terms],quotas=dict(qB=self.qB,qR=self.qR,JB=self.JB,JR=self.JR),fresh_precursor_bound=12*self.N*self.M,gross_each_reservoir_bound=2*(self.JB+self.M*self.JR),elapsed_model_time=str(2*(self.T+5376)))
def rational_witness(self):
if self.gamma!=Q(1,10**11) or self.u not in (128,256):raise ValueError('Paper rational envelopes implemented for u=128 or 256 at gamma=1e-11.')
env=Q(1,10**38);half=Q(1,2**64) if self.u==128 else env;N=self.N;M=self.M;u=self.u;T=self.T
if min(7*u,4*u,15*u/2,u,3*u/2,Q(N,35*10**12),Q(N,2500),Q(19*N,500000),8*u,31*u/2)<128:raise ArithmeticError('Envelope scope failed.')
chemical=(M+14*M+M*T*u/42+M+M*T*u/42+56*M+2)*env
recovery=M*((2+16*u)*env+2*half);transfer=Q(8000000,M)
bound=2*(chemical+recovery+transfer+2*self.service)
return dict(nonservice_uniform_bound=str(chemical+recovery+transfer),two_cycle_failure_bound=str(bound),success_lower=str(max(0,1-bound)),scope='Exact rational envelope uses e>=2, e^-128<=1e-38 and e^-64<=2^-64, with uniform p=1/100 in both cycles.')
class OddsAccounting:
@staticmethod
def coordinates(cells):
counts={t:sum(c.tag==t for c in cells) for t in ('H','L')};sizes={t:sum(c.size for c in cells if c.tag==t) for t in ('H','L')}
if not all(counts.values()):return dict(counts=counts,sizes=sizes,log_count=None,log_size=None,phase=None)
C=mp.log(mp.mpf(counts['H'])/counts['L']);S=mp.log(mp.mpf(sizes['H'])/sizes['L']);phase=mp.log(mp.mpf(sizes['H'])*counts['L']/(sizes['L']*counts['H']))
return dict(counts=counts,sizes=sizes,log_count=mp.nstr(C,30),log_size=mp.nstr(S,30),phase=mp.nstr(phase,30))
@staticmethod
def theorem_thresholds():
def values(ctx):
g=number(ctx,'3/5')*ctx.log(4)-number(ctx,'19/500')-ctx.log(number(ctx,'51/49'));one=g-ctx.log(2);two=2*g-ctx.log(2)
return g,one,two,1/(1+ctx.exp(-two))
g,one,two,f=values(mp);enclosures=values(iv)
if not bool(enclosures[3].a>number(iv,'0.6930453').b):raise ArithmeticError('Decimal fraction threshold unresolved.')
return dict(cycle_size_gain=mp.nstr(g,30),one_count_gain=mp.nstr(one,30),two_count_gain=mp.nstr(two,30),two_high_fraction=mp.nstr(f,30),two_fraction_certified_lower='0.6930453',two_count_interval=str(enclosures[2]),two_fraction_interval=str(enclosures[3]),eligibility_floors=['1/2','49/1600','2401/1280000'],exact_two_gain='(7/5)log(2)-19/250-2log(51/49)',scope='Two cycles from balanced ready newborns only. Intermediate division phases telescope; no separate log(2) penalty per cycle.')
def main():
parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));parser.add_argument('--simulate',action='store_true');args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
model=ChemicalModel();regions=ChemicalRegions();certificate=TwoCycleCertificate();founders={t:regions.newborn(t,NEWBORN_SIZE) for t in ('L','H')}
transfer=UniformTransfer();rows=[]
for name,sizes in [('equal',[10]*8),('high_larger',[19]*4+[10]*4),('low_larger',[10]*4+[19]*4),('mixed',[19,10,19,10,10,19,10,19])]:
cells=tuple(Cell('H' if i<4 else 'L',size,(20,20,30 if i<4 else 10,40)) for i,size in enumerate(sizes))
for eps in (DIAGNOSTIC_TRANSFER_TOLERANCE,'1/2'):rows.append(dict(endpoint=name,epsilon=eps,**transfer.enumerate(cells,4,eps)))
# Actual intact selected states, at large size, are followed by a density recovery illustration.
endpoint=tuple(replace(founders['H' if i<4 else 'L'],counts=tuple(n+(NEWBORN_SIZE//50000 if j==2 else 0) for j,n in enumerate(founders['H' if i<4 else 'L'].counts))) for i in range(8))
selected,discarded=transfer.sample(endpoint,4,np.random.default_rng(RANDOM_SEED));recovery=[];difference=0.
for index,cell in enumerate(selected):
sol=model.recovery_density(cell);other=model.recovery_density(cell,method='BDF');times=np.unique(np.r_[0,np.geomspace(.001,5376,180)])
values=sol.sol(times).T;difference=max(difference,float(np.max(abs(values-other.sol(times).T))))
recovery.extend((index,cell.tag,float(t),*x) for t,x in zip(times,values))
result=dict(chemistry=model.certificates(),founders={t:dict(cell=asdict(c),energy=regions.energy(c).json(),readout=c.readout,multiplicity=RETAINED_CELLS//2) for t,c in founders.items()},configured=certificate.evaluate(),paper_witnesses=[dict(evaluated=c.evaluate(),rational=c.rational_witness()) for c in (TwoCycleCertificate(65536*10**18,4*10**9,'1/100000000000','1/1000'),TwoCycleCertificate(131072*10**18,10**15,'1/100000000000','1/1000'))],thresholds=OddsAccounting.theorem_thresholds(),exact_transfer=rows,
selected_recovery=dict(selected=[asdict(c) for c in selected],discarded=[asdict(c) for c in discarded],selected_energy_enclosures=[regions.energy(c).json() for c in selected],solver_difference=difference,scope='Constructed admitted endpoint population, then intact sampling and full resident density ODE on the actual selected count/size inputs; fixed size and no precursor. Not a batch simulation or the finite-count recovery probability.'),scope='Bound evaluation imports the paper probability theorems; no Lean replay. Material, rate ratios, stationary identities, energy matrices and tiny transfer laws are freshly checked.')
if args.simulate:
cells=[]
for i in range(DIAGNOSTIC_CELLS):
tag='H' if i<DIAGNOSTIC_CELLS//2 else 'L';coords=regions.centers[tag];cells.append(Cell(tag,DIAGNOSTIC_SIZE,tuple(math.floor(DIAGNOSTIC_SIZE*(v.lo+v.hi)/2) for v in coords)))
result['uncertified_small_mission']=SerialProtocol(PopulationSource(DIAGNOSTIC_SIZE,DIAGNOSTIC_GROWTH),DIAGNOSTIC_CELLS,float(Q(DIAGNOSTIC_RECOVERY))).run(cells,2,RANDOM_SEED,EVENT_BUDGET)
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('recovery.csv',['selected_index','ancestry','time',*SPECIES],recovery)
table('transfer.csv',['endpoint','epsilon','subsets','bad','failure'],[(r['endpoint'],r['epsilon'],r['subsets'],r['bad'],r['failure']) for r in rows])
sweep=[]
for N in [int(x)*10**18 for x in np.linspace(140,70000,71)]:
c=TwoCycleCertificate(N,RETAINED_CELLS);d=c.logs('1/100');groups={k:mp.log(sum(mp.exp(v) for name,v in d.items() if name in names))/mp.log(10) for k,names in {'chemical':('outer_initial','daughter_initial','outer_drift','division_initial','division_drift','partition','deadline','size_odds'),'recovery':('recovery_initial','recovery_return','recovery_exit','recovery_exit_drift'),'transfer':('transfer',)}.items()};sweep.append((N,*map(float,groups.values())))
table('error_sweep.csv',['N','chemical_log10','recovery_log10','transfer_log10'],sweep)
lines=[f'Thirteen resident channels, balanced completion, stationary curves and exact energy sandwich checked.',f'Configured two-cycle success >= {result["configured"]["success_lower_certified"]}; measured high fraction strictly above {result["thresholds"]["two_fraction_certified_lower"]} on success.',f'Exact eight-cell transfer failures: {[r["failure"] for r in rows]}.',f'Actual selected states retained for precursor-free density recovery; independent solver difference {difference:.3g}.','The certificate covers two cycles from ready newborn founders, not indefinite selection. Small accelerated simulations are separately labeled and outside its hypotheses.']
(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.2),layout='constrained');data=np.array(sweep)
for j,label in enumerate(['Batch chemistry','Recovery','Transfer']):axs[0].plot(data[:,0]/1e22,data[:,j+1],label=label)
axs[0].set(xlabel='Newborn size N / 10^22',ylabel='log10 per-cycle error upper bound',ylim=(-50,20),title='Batch, transfer and recovery error bounds');axs[0].legend(fontsize=8)
x=np.arange(4);axs[1].bar(x-.18,[float(Q(rows[2*i]['failure'])) for i in range(4)],.36,label=f'Tolerance {DIAGNOSTIC_TRANSFER_TOLERANCE}');axs[1].bar(x+.18,[float(Q(rows[2*i+1]['failure'])) for i in range(4)],.36,label='Tolerance 1/2');axs[1].set_xticks(x,['Equal','High larger','Low larger','Mixed']);axs[1].set(ylabel='Exact transfer failure probability',title='Four of eight intact cells, all 70 subsets');axs[1].legend(fontsize=8)
for ax in axs:ax.grid(axis='y',alpha=.2)
fig.savefig(out/'selection.png',dpi=180);fig.savefig(out/'selection.svg');plt.close(fig)
fig,ax=plt.subplots(figsize=(8,4.2),layout='constrained')
for tag,color in [('L','#417d8c'),('H','#bd5a24')]:
for i,c in enumerate(selected):
if c.tag!=tag:continue
data=np.array([r[2:] for r in recovery if r[0]==i],float);center=float((regions.centers[tag][2].lo+regions.centers[tag][2].hi)/2)
ax.semilogx(np.maximum(data[:,0],1e-3),data[:,3]-center,color=color,label=f'Selected {tag}, index {i}')
ax.set(xlabel='Precursor-free recovery time',ylabel='z concentration minus stationary z',title='Chemical recovery after intact-cell\ntransfer');ax.legend(fontsize=8);ax.grid(alpha=.2);ax.ticklabel_format(axis='y',style='sci',scilimits=(0,0))
fig.savefig(out/'recovery.png',dpi=180);fig.savefig(out/'recovery.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__)),resident_sha256=digest(Path(__file__).with_name('resident.py')),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
Thirteen resident channels, balanced completion, stationary curves and exact energy sandwich checked. Configured two-cycle success >= 0.993959999; measured high fraction strictly above 0.6930453 on success. Exact eight-cell transfer failures: ['17/35', '1/35', '17/35', '1/35', '17/35', '1/35', '27/35', '13/35']. Actual selected states retained for precursor-free density recovery; independent solver difference 1.29e-10. The certificate covers two cycles from ready newborn founders, not indefinite selection. Small accelerated simulations are separately labeled and outside its hypotheses.