Example code
A chemical mixture can store a label in a stable concentration state. Copying that label through division requires a parent within the specified operating region to divide and both daughters to recover into a region where copying can repeat. This example makes that complete cycle explicit and keeps the paper's three different chemical models separate.
Division is complementary. If one daughter receives of a species, the other receives from the parent's molecules. Their return probabilities cannot be multiplied as though the siblings were independent. The code sums the actual joint allocation law exactly and propagates state-dependent probabilities of daughter outcomes. Uniform return gives a bound along one lineage, but for every division in a binary family of depth .
For the constructed two-module chemistry, the package supplies every reversible resident reaction, shared consuming growth, cross-catalysis, domain monitor and event quota. Each of the four saved stochastic paths divides and returns both daughters to its specified ranges of molecule counts. These paths illustrate the protocol; they do not estimate its uniform 0.9901 guarantee. The advertised 874 molecules counts newborn residents per module, with controlled pools and other resources accounted for separately.
The four-species model addresses a different question: how redundancy grows with module count and copying horizon. Its reusable reaction and exchange generator implements the specified rate scaling. The analytic bound needs logarithmically increasing copy numbers for a lineage, versus linear growth with family depth. The sufficient constants are enormous: about copies even for the displayed one-division case.
The nominal eight-species model permits repair after split and refill. Both deterministic recovery curves return to the specified ellipsoidal recovery regions after 500 seconds. The thiol threshold decodes the label, but repeatable copying requires the stronger recovery-region test. The package reproduces the rational scalar budget underlying the 0.992 joint-return bound and shows how small rate changes can leave a nominal region while preserving the decoder's answer.
Download the complete package for editable inputs, exact allocation tools, reusable count and concentration models, region data, sensitivity outputs, seven scientific test groups, and worked import examples. The scalar budget checks use the paper's certified generator and polynomial inequalities as premises: this package does not recheck the large certificate tables. Its stochastic paths and deterministic integrations are clearly distinguished from those uniform probability guarantees.
Python source
"""Chemical copying as a complete cycle: return regions, allocation, recovery.
Three manuscript models remain distinct. Numerical trajectories do not replace
the paper's large generator certificates; scalar consequences are checked here.
"""
from __future__ import annotations
from dataclasses import dataclass
from fractions import Fraction as Q
from pathlib import Path
from functools import lru_cache
import argparse
import csv
import hashlib
import itertools
import json
import math
import platform
import mpmath as mp
import numpy as np
from scipy.integrate import solve_ivp
# EDITABLE INPUTS: constructed source (counts, seconds, prescribed volume).
BIRTH_MARKER = 53
CROSS_CATALYSIS = 0.05
REVERSE_GROWTH = 1e-7
DIVISION_DEADLINE = 20.0
PHYSICAL_EVENT_QUOTA = 100_000_000
COMPUTE_EVENT_LIMIT = 500_000 # unfinished is not a scientific failure
COPY_SEED = 17092026
PREPARATIONS = ((53,0),(530,53))
# Redundancy theorem: different four-species chemistry, different copy number.
MODULE_COUNT = 1
GROWTH_RATE = Q(1,10**11)
TARGET_ERROR = Q(1,100)
HORIZONS = tuple(range(1,41))
# Nominal eight-species chemistry: molar concentrations, seconds.
SPECIES = ('S','C','P','E','D','U','V','I')
FEED = (0.05,0,0,0,0.1,0,0,0.00231)
DILUTION = 0.002
EXCHANGE = 0.65
LIGATION = 0.41
HYDROLYSIS = 9.26e-6
INHIBITION = 150.0
RECOVERY_TIME = 500.0
OMEGA = 24*10**18 # molecules per molar; NOT directly simulated
PERTURBATION = 0.001 # nominal-region diagnostic, not uncertainty certification
MANUSCRIPT_SHA256 = '50450e40c8f43be86391b39a88460038b1d3d931ce7d2bba0be788916d1af07e'
REGION_SHA256 = 'e214ee720185c3e078d4880be7bfa247a473c4084165182477be5f7fae379e5a'
@dataclass(frozen=True)
class CountBox:
lower: tuple[int, ...]
upper: tuple[int, ...]
def __post_init__(self):
if len(self.lower) != len(self.upper) or any(not 0 <= lo <= hi for lo,hi in zip(self.lower,self.upper)):
raise ValueError('Nonempty nonnegative coordinate intervals required.')
def contains(self, counts):
return len(counts) == len(self.lower) and all(lo <= n <= hi for n,lo,hi in zip(counts,self.lower,self.upper))
class ComplementaryAllocation:
@staticmethod
@lru_cache(maxsize=4096)
def coordinate_return(total, lower, upper):
"""Exact P(lower <= D,total-D <= upper), D~Bin(total,1/2)."""
if not isinstance(total,int) or not 0 <= total <= 20000 or lower > upper:
raise ValueError('Exact allocation guard: integer total in [0,20000].')
lo,hi = max(0,lower,total-upper),min(total,upper,total-lower)
return Q(sum(math.comb(total,i) for i in range(lo,hi+1)),2**total)
@classmethod
def joint_return(cls, parent, box):
if len(parent) != len(box.lower):
raise ValueError('Allocation dimension differs from admission box.')
return math.prod(cls.coordinate_return(int(n),lo,hi) for n,lo,hi in zip(parent,box.lower,box.upper))
@staticmethod
def draw(parent,rng):
first = rng.binomial(np.asarray(parent,dtype=np.int64),.5)
return first,np.asarray(parent,dtype=np.int64)-first
class ConstructedSource:
boxes = (CountBox((6,0),(106,8)),CountBox((371,6),(742,132)))
def __init__(self,epsilon=CROSS_CATALYSIS,rho=REVERSE_GROWTH,birth=BIRTH_MARKER):
if epsilon < 0 or rho < 0 or birth < 1:
raise ValueError('Nonnegative interactions and positive marker required.')
self.epsilon,self.rho,self.birth = epsilon,rho,birth
jumps = []
kinds = []
resident = ((1,0),(-1,0),(2,-1),(-2,1),(-1,1),(1,-1),(-1,0),(1,0),(-1,0),(1,0))
for i in range(2):
for dx,dy in resident:
jump = [0]*5; jump[2*i:2*i+2] = [dx,dy]
jumps.append(jump);kinds.append('resident')
for dx,dm,kind in ((-1,1,'growth'),(1,-1,'reverse_growth')):
jump = [0]*5;jump[2*i]=dx;jump[4]=dm
jumps.append(jump);kinds.append(kind)
for dx,dy in ((-1,1),(1,-1)):
jump = [0]*5;jump[2*i:2*i+2]=[dx,dy]
jumps.append(jump);kinds.append('cross_catalysis')
jumps.append([0]*5);kinds.append('self_event')
self.jumps,self.kinds = np.array(jumps,dtype=np.int64),tuple(kinds)
def rates(self,state,fluid=False):
m = state[4]
if m <= 0:
raise ValueError('Positive marker required.')
rates = []
for i in range(2):
x,y = map(float,state[2*i:2*i+2]);other_y=float(state[2*(1-i)+1])
xx = x*x if fluid else x*max(0,x-1)
rates.extend((193*m/5,193*x/750000,1555*y,311*xx/(30000*m),
15*xx/m,15*x*y/m,100*x*y/m,y/1000,105*x/2,21*m/40000,
x/10,self.rho*m,self.epsilon*x*other_y/m,self.epsilon*y*other_y/m))
return np.array([*rates,1.0])
def admitted(self,state,word):
return state[4] == self.birth and all(self.boxes[label].contains(state[2*i:2*i+2]) for i,label in enumerate(word))
def run(self,word,seed=COPY_SEED,initial=None,limit=COMPUTE_EVENT_LIMIT,
deadline=DIVISION_DEADLINE,quota=PHYSICAL_EVENT_QUOTA):
if len(word) != 2 or any(label not in (0,1) for label in word) or limit < 1 or quota < 1 or deadline <= 0:
raise ValueError('Two binary labels and positive limits required.')
state = np.array((*PREPARATIONS[word[0]],*PREPARATIONS[word[1]],self.birth) if initial is None else initial,dtype=np.int64)
if not self.admitted(state,word):
raise ValueError('Initial parent must be admitted.')
rng,time = np.random.default_rng(seed),0.0
trace = [(time,*state.tolist())]
status = 'unfinished'
for step in range(limit):
rates = self.rates(state)
total = rates.sum()
time_next = time+rng.exponential(1/total)
if time_next > deadline:
time=deadline;status='deadline';break
time=time_next
event=int(np.searchsorted(np.cumsum(rates),rng.random()*total,side='right'))
state += self.jumps[event]
assert np.min(state) >= 0
if step+1 >= quota:
status='event_quota';break
if self.kinds[event] == 'reverse_growth':
status='reverse_growth_monitor';break
if any(state[2*i] > 16*state[4] or state[2*i+1] > 4*state[4] for i in range(2)):
status='domain_exit';break
if state[4] == 2*self.birth:
status='division';break
if step % 2000 == 0:
trace.append((time,*state.tolist()))
trace.append((time,*state.tolist()))
result={'word':list(word),'status':status,'time':time,'events':step+1,
'parent_at_stop':state.tolist(),'joint_partition_return':None,'both_return':None,'trace':trace}
if status == 'division':
payoff = math.prod(ComplementaryAllocation.joint_return(state[2*i:2*i+2],self.boxes[label]) for i,label in enumerate(word))
A,B = ComplementaryAllocation.draw(state[:4],rng)
result.update(joint_partition_return=str(payoff),
both_return=all(self.boxes[label].contains(d[2*i:2*i+2]) for i,label in enumerate(word) for d in (A,B)),
daughters=[[*A.tolist(),self.birth],[*B.tolist(),self.birth]])
elif status != 'unfinished':
result['both_return']=False
return result
@staticmethod
def paper_budget():
"""Scalar assembly conditional on the paper's checked generator rows.
This function does NOT regenerate/check the 22-million-row certificate.
"""
premonitor=Q(9911,10000)-20*Q(2,10**6)-Q(1,50000)/Q(5,8)
quota=Q(11,4)**80/Q(27,10)**100
assert Q(2718,1000)**25 > 250000**2
assert Q(101,100)**100 >= Q(27,10) and quota < Q(2,10**7)
final=Q(991,1000)-80*Q(2,10**7)*53-Q(2,10**7)
return {'premonitor':str(premonitor),'quota_upper_bound':str(quota),
'assembled_lower_bound':str(final),'reported_lower_bound':'9901/10000',
'generator_table_rechecked':False,'capacity_per_pool_molecules':4000*53+1001*10**8,
'maximum_newborn_residents_per_module':874}
class OffspringKernel:
"""Finite joint daughter kernel; missing row mass is failed copying."""
def __init__(self,states,rows):
self.states,self.rows=tuple(states),rows
if set(rows) != set(states):
raise ValueError('One row per admitted restart state required.')
for row in rows.values():
if any(a not in states or b not in states or p < 0 for (a,b),p in row.items()) or sum(row.values()) > 1:
raise ValueError('Invalid substochastic joint row.')
def success(self,generations,family=False):
if generations < 0:
raise ValueError('Nonnegative horizon required.')
values={s:Q(1) for s in self.states}
for _ in range(generations):
values={s:sum((p*values[a]*(values[b] if family else 1)
for (a,b),p in self.rows[s].items()),Q(0)) for s in self.states}
return values
@property
def uniform_one_cycle(self):
return min(sum(row.values()) for row in self.rows.values())
class RedundancyBounds:
def __init__(self,k=MODULE_COUNT,gamma=GROWTH_RATE):
if not isinstance(k,int) or k < 1 or not 0 < gamma <= Q(1,10**11):
raise ValueError('The four-species bound requires k >= 1 and 0 < gamma <= 1e-11.')
self.k,self.gamma=k,gamma
def required(self,depth,error=TARGET_ERROR,family=False):
if not isinstance(depth,int) or depth < 1 or not 0 < error < 1:
raise ValueError('Positive horizon and error in (0,1) required.')
with mp.workdps(80):
divisions=(1 << depth)-1 if family else depth
eta=mp.mpf(error.numerator)/error.denominator
C=2703+mp.mpf(11)*self.gamma.denominator/(5*self.gamma.numerator)
sufficient=max(14*10**19,int(mp.ceil(1024*10**18*mp.log(self.k*divisions*C/eta))))
necessary=mp.log(self.k*divisions/(-mp.log1p(-eta)))/(280*mp.log(2))
return {'depth':depth,'divisions':divisions,'sufficient_N':sufficient,
'necessary_real_lower_bound':mp.nstr(max(0,necessary),30),
'criterion':'immediate_return','family':family}
class FourSpeciesModule:
"""Reference chemistry for the redundancy theorem: A, B, z, H."""
def __init__(self,e=1e-5,d0=1e-4):
self.e,self.d0=e,d0
# Consumption, production, reference coefficient; source/sink channels
# represent the theorem's open chemistry, not the controlled pools above.
pairs=(((0,),(1,2),1.,1.),((2,),(3,),16.,1.),
((2,2),(3,),2.,1.),((),(0,),6.,1.),((),(1,),27.,1.),
((1,),(0,0),e,e))
self.reactions=tuple(r for a,b,k,l in pairs for r in ((a,b,k),(b,a,l)))+(((3,),(),d0),)
def fluid(self,u):
A,B,z,H=u;e,d0=self.e,self.d0
return np.array((6-2*A+z*B+2*e*(B-A*A),27+A-(1+z)*B-e*(B-A*A),
A-B*z-16*z-4*z*z+3*H,16*z+2*z*z-(2+d0)*H))
@staticmethod
def stationary_reduction(z):
e,d0=Q(1,100000),Q(1,10000)
B=60/(z+2);H=(16*z+2*z*z)/(2+d0)
K=(2*(1+2*d0)*z*z-16*(1-d0)*z)/(2+d0)
A=z*B+K
return (A,B,z,H),27-(1+e)*B+K+e*A*A
class CoupledModules:
"""Literal finite-count generator with shared consuming growth and exchange.
Physical coefficients scale with k exactly as in this theorem. Building a
generator does not assert that arbitrary parameters satisfy local quality.
"""
def __init__(self,k,weights,gamma=float(GROWTH_RATE),module=None):
weights=np.asarray(weights,dtype=float)
if (k < 1 or weights.shape != (k,k) or np.any(weights < 0)
or not np.array_equal(weights,weights.T) or np.any(np.diag(weights) != 0) or gamma <= 0):
raise ValueError('Symmetric nonnegative exchange without self-edges required.')
self.k,self.weights,self.gamma,self.module=k,weights,gamma,module or FourSpeciesModule()
def rates_and_jumps(self,counts,marker):
counts=np.asarray(counts,dtype=np.int64)
if counts.shape != (self.k,4) or np.min(counts) < 0 or marker <= 0:
raise ValueError('Nonnegative module counts and positive shared marker required.')
rates,jumps=[],[]
for i in range(self.k):
for inputs,outputs,coefficient in self.module.reactions:
order=len(inputs)
rate=coefficient*self.k**(order-2)*float(marker)**(1-order)
used={};jump=np.zeros(4*self.k+1,dtype=np.int64)
for x in inputs:
rate*=max(0,int(counts[i,x])-used.get(x,0));used[x]=used.get(x,0)+1;jump[4*i+x]-=1
for x in outputs:jump[4*i+x]+=1
rates.append(rate);jumps.append(jump)
jump=np.zeros(4*self.k+1,dtype=np.int64);jump[4*i+2]=-1;jump[-1]=1
rates.append(self.gamma*counts[i,2]/self.k);jumps.append(jump)
for j in range(self.k):
if self.weights[i,j]:
jump=np.zeros(4*self.k+1,dtype=np.int64);jump[4*i+2]=-1;jump[4*j+2]=1
rates.append(self.weights[i,j]*counts[i,2]/self.k);jumps.append(jump)
return np.array(rates),np.array(jumps)
def interaction_within_paper_allowance(self):
return bool(np.max(self.weights.sum(axis=1)) <= 1e-11 and self.gamma <= 1e-11)
@dataclass(frozen=True)
class ReturnEllipsoid:
label: str
center: tuple[Q, ...]
metric: tuple[tuple[Q, ...], ...]
eta: Q
omega: int
@classmethod
def load(cls,path=Path(__file__).with_name('terminal_regions.json')):
data=json.loads(path.read_text(encoding='utf-8'))
return tuple(cls(r['label'],tuple(map(Q,r['center_M'])),
tuple(tuple(map(Q,row)) for row in r['terminal_metric']),Q(r['eta']),r['volume_scale'])
for r in data['regions'])
def exact_count_energy(self,counts):
if len(counts) != len(self.center):
raise ValueError('Wrong region dimension.')
y=tuple(Q(int(n),self.omega)-z for n,z in zip(counts,self.center))
if len(y) != len(self.center):
raise ValueError('Wrong region dimension.')
return sum(y[i]*self.metric[i][j]*y[j] for i in range(8) for j in range(8))/self.eta
def energy(self,concentration):
y=np.asarray(concentration)-np.array(self.center,dtype=float)
return float([email protected](self.metric,dtype=float)@y/float(self.eta))
def cube_tolerance(self):
return math.sqrt(float(self.eta/sum(abs(v) for row in self.metric for v in row)))
@dataclass(frozen=True)
class SemenovModel:
exchange: float=EXCHANGE
ligation: float=LIGATION
hydrolysis: float=HYDROLYSIS
inhibition: float=INHIBITION
dilution: float=DILUTION
feed: tuple[float,...]=FEED
def __post_init__(self):
if (min(self.exchange,self.ligation,self.hydrolysis,self.inhibition,self.dilution) <= 0
or len(self.feed) != 8 or min(self.feed) < 0):
raise ValueError('Positive kinetic coefficients and eight nonnegative feeds required.')
def chemistry(self):
# Input/output species indices; inactive products remain outside tracking.
pairs=(((4,3),(1,5)),((4,2),(1,6)),((6,3),(2,5)))
reactions=[]
for left,right in pairs:
reactions.extend(((left,right,self.exchange),(right,left,self.exchange)))
reactions += [((0,1),(2,3),self.ligation),((0,),(3,),self.hydrolysis)]
reactions += [((x,7),(),self.inhibition) for x in (1,2,3)]
return tuple(reactions)
def rhs(self,time,u):
derivative=self.dilution*(np.array(self.feed)-u)
for inputs,outputs,k in self.chemistry():
rate=k*math.prod(u[i] for i in inputs)
for i in inputs: derivative[i]-=rate
for i in outputs: derivative[i]+=rate
return derivative
def recovery(self,region,rtol=2e-10,atol=2e-14):
initial=(np.array(region.center,dtype=float)+np.array(self.feed))/2
times=np.linspace(0,RECOVERY_TIME,251)
answer=solve_ivp(self.rhs,(0,RECOVERY_TIME),initial,t_eval=times,method='BDF',rtol=rtol,atol=atol)
if not answer.success: raise RuntimeError(answer.message)
return times,answer.y.T
def jump_rates(self,counts,omega):
"""Literal 11 chemical + 8 feed + 8 outflow rates and jumps.
This reusable kernel is not a practicable SSA at the certified volume.
"""
rates,jumps=[],[]
for inputs,outputs,k in self.chemistry():
used={};rate=k*float(omega)**(1-len(inputs));jump=[0]*8
for i in inputs:
rate*=max(0,int(counts[i])-used.get(i,0));used[i]=used.get(i,0)+1;jump[i]-=1
for i in outputs:jump[i]+=1
rates.append(rate);jumps.append(jump)
for i in range(8):
jump=[0]*8;jump[i]=1;rates.append(float(omega)*self.dilution*self.feed[i]);jumps.append(jump)
for i in range(8):
jump=[0]*8;jump[i]=-1;rates.append(self.dilution*int(counts[i]));jumps.append(jump)
return np.array(rates),np.array(jumps)
def semenov_scalar_budget(label):
"""Exact endpoint budget, conditional on the supplied paper's tube bounds."""
eta,L0,vs,r=((Q(6,10**12),1601,Q(115,1000),Q(9,100)),
(Q(9,10**14),3818,Q(114,1000),Q(26,100)))[label]
omega,K=24*10**18,12*10**18
m0=Q(omega)*Q(15231,100000)/2
lam=Q(omega)*Q(15231,100000)*Q(2,1000)
m1=L0*vs/(omega*eta)
m2=Q(L0**2)/eta**2*(24*vs**2/omega**2+8*vs/omega**3)
energy=r**6
terms={'inherited_parent':2*energy-energy**2,'second_moment':(4+48*r)*m1,
'fourth_moment':12*m2,'refill_inventory':64*m0/K**2,
'recovery_and_inventory_drift':1000*(Q(1,10**9)+4*lam/K**2)}
total=sum(terms.values())
assert total <= Q(1,250)
return {'terms':{k:str(v) for k,v in terms.items()},'joint_failure_upper_bound':str(2*total),
'polynomial_tube_certificates_rechecked':False}
def write_csv(path,headers,rows):
with path.open('w',newline='') as stream:
writer=csv.writer(stream);writer.writerow(headers);writer.writerows(rows)
def main():
parser=argparse.ArgumentParser(description=__doc__)
parser.add_argument('--output',type=Path,default=Path('outputs'))
args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
source=ConstructedSource()
copies=[]
for index,word in enumerate(itertools.product((0,1),repeat=2)):
result=source.run(word,COPY_SEED+index)
write_csv(out/f'copy_{word[0]}{word[1]}.csv',['time','X0','Y0','X1','Y1','marker'],result.pop('trace'))
copies.append(result)
regions=ReturnEllipsoid.load()
if any(region.omega != OMEGA for region in regions):
raise ValueError('Changing OMEGA requires compatible return-region data; nominal certification does not transfer automatically.')
model=SemenovModel()
from dataclasses import replace
recoveries=[]
trajectories=[]
sensitivity=[]
for label,region in enumerate(regions):
times,trajectory=model.recovery(region)
_,fine=model.recovery(region,rtol=2e-12,atol=2e-16)
preparation=[int(region.omega*z) for z in region.center]
recoveries.append({'label':region.label,'endpoint_energy':region.energy(fine[-1]),
'thiol_endpoint_molar':float(sum(fine[-1,1:4])),
'tolerance_endpoint_difference':float(np.max(np.abs(fine[-1]-trajectory[-1]))),
'cube_tolerance_molar':region.cube_tolerance(),'integer_preparation':preparation,
'exact_preparation_energy':str(region.exact_count_energy(preparation)),
'scalar_budget':semenov_scalar_budget(label)})
trajectories.append(fine)
write_csv(out/f'recovery_{region.label}.csv',['time',*SPECIES],zip(times,*fine.T))
for group in ('exchange','ligation','hydrolysis','inhibition','dilution'):
for sign in (-1,1):
changed=replace(model,**{group:getattr(model,group)*(1+sign*PERTURBATION)})
_,run=changed.recovery(region)
sensitivity.append({'label':region.label,'group':group,'relative_change':sign*PERTURBATION,
'terminal_energy_in_nominal_region':region.energy(run[-1]),
'decoded_high':bool(sum(run[-1,1:4]) > .005)})
redundancy=[RedundancyBounds().required(G,family=family) for family in (False,True) for G in HORIZONS]
kernel=OffspringKernel(('interior','edge'),{
'interior':{('interior','interior'):Q(9,10),('edge','edge'):Q(9,100)},
'edge':{('interior','edge'):Q(1,2),('edge','interior'):Q(12,25)}})
kernel_result={'role':'illustrative state-dependent kernel, not fitted chemical data',
'uniform_one_cycle':str(kernel.uniform_one_cycle),
'lineage_5':{s:str(p) for s,p in kernel.success(5).items()},
'family_depth_5':{s:str(p) for s,p in kernel.success(5,True).items()}}
root_brackets=[]
for left,right in (('0.99579401232','0.99579401233'),('2.97636724376','2.97636724377')):
_,fl=FourSpeciesModule.stationary_reduction(Q(left))
_,fr=FourSpeciesModule.stationary_reduction(Q(right))
assert fl*fr < 0
root_brackets.append({'z_interval':[left,right],'residual_signs':[1 if fl > 0 else -1,1 if fr > 0 else -1]})
result={'constructed_runs':copies,'constructed_scalar_budget':source.paper_budget(),
'four_species_stationary_brackets':root_brackets,
'nominal_recovery':recoveries,'sensitivity':sensitivity,'redundancy':redundancy,'joint_kernel_example':kernel_result}
(out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
lines=[f'Constructed word {r["word"]}: {r["status"]} at {r["time"]:.6g} s; both admitted={r["both_return"]}.' for r in copies]
lines += [f'Nominal {r["label"]} recovery: terminal energy {r["endpoint_energy"]:.8g}; cube tolerance {r["cube_tolerance_molar"]*1e9:.6g} nM.' for r in recoveries]
lines += ['Scalar probability budgets reproduced; the large generator/polynomial certificate tables are not rechecked.',
'Four SSA paths and deterministic recovery curves do not estimate the certified uniform copying probability.']
(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=(11,4.4),layout='constrained')
for region,trajectory in zip(regions,trajectories):
axes[0].plot(times,trajectory[:,1:4].sum(axis=1)*1000,label=region.label)
axes[0].axhline(5,color='black',ls='--',label='Thiol decoder threshold')
axes[0].set(xlabel='Recovery time (s)',ylabel='Free thiol (mM)',title='Nominal split, refill and recovery')
for family,label in ((False,'Designated lineage'),(True,'Complete binary family')):
rows=[r for r in redundancy if r['family']==family]
axes[1].plot([r['depth'] for r in rows],[r['sufficient_N']/1e22 for r in rows],label=label)
axes[1].set(xlabel='Lineage divisions / family depth',ylabel='Sufficient N / 10^22',title='Different model: four-species redundancy bound')
for ax in axes:ax.legend(fontsize=8);ax.grid(alpha=.2)
fig.savefig(out/'copying.png',dpi=180);fig.savefig(out/'copying.svg');plt.close(fig)
fig,axes=plt.subplots(1,2,figsize=(11,4.4),layout='constrained')
totals=range(241)
joint=[float(ComplementaryAllocation.coordinate_return(n,6,106)) for n in totals]
marginal=[sum(math.comb(n,d) for d in range(6,min(106,n)+1))/2**n for n in totals]
axes[0].plot(totals,joint,label='Both complementary daughters admitted')
axes[0].plot(totals,np.square(marginal),'--',label='Product of marginals (incorrect for siblings)')
axes[0].set(xlabel='Parent X count at division',ylabel='Probability for the X coordinate',
title='One coordinate: both daughters need 6 to 106')
horizons=np.arange(1,11)
axes[1].plot(horizons,.9901**horizons,label='Designated lineage: p^G')
axes[1].plot(horizons,.9901**(2**horizons-1),label='Complete family: p^(2^G - 1)')
axes[1].set(xlabel='Lineage divisions / family depth',ylabel='Guaranteed lower bound',
yscale='log',title='Uniform one-cycle bound p = 0.9901')
for ax in axes:ax.legend(fontsize=8);ax.grid(alpha=.2)
fig.savefig(out/'allocation.png',dpi=180);fig.savefig(out/'allocation.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__)),'region_source_sha256':digest(Path(__file__).with_name('terminal_regions.json')),
'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
Constructed word [0, 0]: division at 3.46622 s; both admitted=True. Constructed word [0, 1]: division at 0.574302 s; both admitted=True. Constructed word [1, 0]: division at 0.572689 s; both admitted=True. Constructed word [1, 1]: division at 0.44585 s; both admitted=True. Nominal low recovery: terminal energy 4.109884e-13; cube tolerance 16.3665 nM. Nominal high recovery: terminal energy 1.4076979e-09; cube tolerance 1.78818 nM. Scalar probability budgets reproduced; the large generator/polynomial certificate tables are not rechecked. Four SSA paths and deterministic recovery curves do not estimate the certified uniform copying probability.