Example code
In this persister-cell model, treatment erases inherited molecular marks that affect survival. The example follows activating and repressive marks through chemical changes, death and division, separating treatment duration from its maximum erasure rate. Daughters share the mother's old marks through a correlated allocation law; treating them as independent draws changes extinction probabilities.
Three distinct checks bound what the intervention can achieve: a remaining-exposure barrier, an amplitude barrier valid against every admitted feedback policy, and a finite-time death-rate floor. At two sites, a freshly replayed Taylor enclosure bounds survival below 0.009819 for one fully activating founder, with exposure 11.6 and a deadline of 100 model-time units. At six sites, amplitude 0.29 instead leaves a certified survival floor of 0.45962721, regardless of budget or duration.
Reusable source, action, schedule and certificate classes expose editable molecular and demographic rates. They support arbitrary site count, predetermined schedules, an optional population-feedback simulator with a pathwise budget, separate withdrawal endpoints and an added-death actuator. Saved exact residuals explain why a certificate passes; changed sources are checked again rather than inheriting old guarantees.
The example reconstructs and checks the published witnesses for two through eight sites, sharpens endpoint weight vectors where their original rounding cannot establish a sign, and reproduces the preparation-dependent failure of erasure monotonicity. It distinguishes mean-growth brackets from a proven universal actuator threshold.
All rates are synthetic. Exposure means additional reaction rate integrated over time; it is neither drug concentration nor administered amount. Numerical schedule curves, rigorous finite computations and the conventional probability theorems remain labeled separately. Lean is not rerun.
Python source
"""Exposure, amplitude and duration are distinct limits in an inherited population."""
# ------------------------- EDITABLE INPUTS -------------------------
SITE_COUNT=2
BASELINE_ERASURE='1/100'
BASAL_WRITING='1/100'
RECRUITED_WRITING='1'
RECRUITED_ERASURE='1/2'
DIVISION_RATE='1/10'
PROTECTED_DIVISION_RATE='1/10'
PROTECTED_DEATH_RATE='1/100'
UNPROTECTED_DEATH_RATE='3/10'
FOUNDERS=1
INITIAL_MARK_COUNTS=(2,0)
TARGET_RISK=.01
ERASER_AMPLITUDE='29/100'
EXPOSURE_BUDGET='58/5'
HORIZON='100'
ADDED_DEATH_AMPLITUDE='29/100'
RUN_FEEDBACK_DEMONSTRATION=False
RANDOM_SEED=20260923
OUTPUT_DIRECTORY='outputs'
# These are synthetic reaction/demographic rates, not a treatment calibration.
# Fixed paper certificate replays below remain separate from edited scenarios.
# ------------------------------------------------------------------
import argparse,csv,hashlib,json,platform
from dataclasses import replace
from fractions import Fraction as F
from pathlib import Path
import numpy as np
from population import MolecularRates,InheritedPopulation,Action,Phase,population_survival,FeedbackSimulator
from certificates import ExposureCertificate,AmplitudeCertificate,WeightedDrift,collatz_bounds,m_matrix_test,mv,population_barrier_generator_identity,finite_time_floor,concentration_lower_bound
MANUSCRIPT_SHA256='1be8ca10a66dcac8fbc1175361298ce9ac2b729529221e7c6a94a8e28687b77c'
def configured_source():
return InheritedPopulation(MolecularRates(SITE_COUNT,F(BASAL_WRITING),F(RECRUITED_WRITING),F(RECRUITED_ERASURE),F(BASELINE_ERASURE),F(DIVISION_RATE),F(PROTECTED_DIVISION_RATE),F(PROTECTED_DEATH_RATE),F(UNPROTECTED_DEATH_RATE)))
def phases_for(exposure,amplitude,horizon,placement):
B,v,T=map(F,(exposure,amplitude,horizon))
if B<0 or v<=0 or T<=0 or B>v*T:raise ValueError('Schedule exceeds the amplitude/horizon exposure capacity.')
if B==0:return [Phase(T,Action())]
pulse=B/v
if placement=='spread':return [Phase(T,Action(B/T))]
on=Phase(pulse,Action(v))
if pulse==T:return [on]
off=Phase(T-pulse,Action())
if placement=='early':return [on,off]
if placement=='late':return [off,on]
raise ValueError('Unknown schedule placement.')
def replay_sites(data):
results=dict(mean_sign_endpoints={},exposure={},amplitude={},added_death={},spectral_at_point={})
for key,row in data['sites'].items():
N=int(key);source=InheritedPopulation(MolecularRates(sites=N));lo,hi=F(row['ec_lower']),F(row['ec_upper'])
al,ah=Action(lo-F(1,100)),Action(hi-F(1,100));wl=tuple(map(F,row['witness_lower']));wh=tuple(map(F,row['witness_upper']))
stored_cwl=collatz_bounds(source,al,wl);stored_cwh=collatz_bounds(source,ah,wh)
# Most released six-decimal Perron vectors straddle zero this close to
# the threshold. Propose finer vectors and accept ONLY exact signs.
def sharpen(action):
A=source.mean_matrix(action);ev,vec=np.linalg.eig(np.array(A,dtype=float))
v=np.real(vec[:,np.argmax(ev.real)]);v*=np.sign(v.sum());v/=v.min()
for digits in [10,12,14]:
w=tuple(F(str(round(float(x),digits))) for x in v)
if min(w)<=0:continue
bracket=collatz_bounds(source,action,w)
if bracket[0]>0 or bracket[1]<0:return w,bracket
raise RuntimeError('Cannot certify endpoint sign with the proposed weights.')
wl,cwl=sharpen(al);wh,cwh=sharpen(ah)
assert cwl[0]>0 and cwh[1]<0 and hi-lo==F(1,10**6)
lowtest=m_matrix_test(source,al);hitest=m_matrix_test(source,ah)
assert not lowtest['subcritical'] and hitest['subcritical']
w=tuple(map(F,row['w_half']));drift=WeightedDrift(w,F(row['gamma_half']),Action(F(49,100)));slack=drift.verify(source)
results['mean_sign_endpoints'][key]=dict(intrinsic_erasure_bracket=(lo,hi),amplitude_bracket=(lo-F(1,100),hi-F(1,100)),
lower_collatz=cwl,upper_collatz=cwh,lower_weights=wl,upper_weights=wh,
published_weight_collatz=(stored_cwl,stored_cwh),lower_M_matrix=lowtest,upper_M_matrix=hitest,
contraction_gamma=drift.gamma,contraction_slack=slack,founder_weight_ratio=w[-1]/min(w),
scope='Fresh endpoint sign certificates. A global first-crossing/uniqueness assertion is not inferred from bisection or endpoint checks.')
for key,row in data['exposure'].items():
source=InheritedPopulation(MolecularRates(sites=int(key)));cert=ExposureCertificate(tuple(map(F,row['q'])),F(row['c']))
results['exposure'][key]=dict(q=cert.q,c=cert.c,eta=1-cert.q[-1],residuals=cert.verify(source),necessary_one_founder_point01=cert.necessary_exposure(1,.01,source.size-1))
for key,row in data['amplitude'].items():
source=InheritedPopulation(MolecularRates(sites=row['N']));cert=AmplitudeCertificate(tuple(map(F,row['q'])),F(row['e_hi'])-F(1,100))
z=[0]*source.size;z[-1]=1
results['amplitude'][key]=dict(N=row['N'],maximum=cert.maximum,q=cert.q,residuals=cert.verify(source),fully_marked_floor=cert.survival_floor_exact(z,source),any_type_floor=1-max(cert.q))
for key,row in data['spectral_at_three_tenths'].items():
source=InheritedPopulation(MolecularRates(sites=int(key)));cw=collatz_bounds(source,Action(F(29,100)),tuple(map(F,row['witness'])))
assert cw==tuple(map(F,[row['lower'],row['upper']]))
results['spectral_at_point'][key]=dict(bounds=cw,scope='Mean growth at one constant action only.')
for key,pair in data['added_death'].items():
N=int(key.split('_')[0][1:]);source=InheritedPopulation(MolecularRates(sites=N));eraser=F(0) if key.endswith('no_eraser') else F(29,100)
lo,hi=map(F,pair);low=m_matrix_test(source,Action(eraser,lo));high=m_matrix_test(source,Action(eraser,hi))
assert not low['subcritical'] and high['subcritical'] and hi<=F(9,100)
results['added_death'][key]=dict(bracket=(lo,hi),lower=low,upper=high,scope='Added death makes the Metzler matrix componentwise smaller, so this endpoint bracket has a monotone threshold interpretation.')
return results
def encode(v):
if isinstance(v,F):return str(v)
if isinstance(v,np.ndarray):return v.tolist()
if isinstance(v,np.generic):return v.item()
raise TypeError(type(v).__name__)
def main():
parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',default=OUTPUT_DIRECTORY);args=parser.parse_args()
out=Path(args.output);out.mkdir(exist_ok=True);here=Path(__file__).parent
def dump(name,data):(out/name).write_text(json.dumps(data,indent=2,default=encode)+'\n',encoding='utf-8')
def table(name,header,rows):
with (out/name).open('w',newline='',encoding='utf-8') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
data=json.loads((here/'site_certificates.json').read_text());sites=replay_sites(data);dump('rechecked_site_certificates.json',sites)
source=InheritedPopulation();q=tuple(map(F,['.94','.982','.987','.27','.735','.235']));exposure=ExposureCertificate(q,F(67,73))
weight=WeightedDrift(tuple(map(F,[1419,1105,1000,8423,4305,10765])),F(21,200),Action(F(29,100)))
exposure_res=exposure.verify(source);weight_slack=weight.verify(source)
assert min(-v for v in exposure_res['baseline_residual'])==F(1,5000)
assert min(exposure_res['control_slack'])==0
rows=[]
for count in [1,10,100,1000,1000000]:
z=[0]*6;z[-1]=count
rows.append(dict(founders=count,necessary=exposure.necessary_exposure(count,TARGET_RISK,5),**weight.sufficient(source,z,TARGET_RISK)))
table('founder_exposure_bounds.csv',['founders','necessary_exposure','sufficient_exposure','sufficient_duration'],[(r['founders'],r['necessary'],r['eraser_exposure'],r['duration']) for r in rows])
# Directly check the exact population generator, not a product assumption.
identity=population_barrier_generator_identity(source,q,[2,0,1,3,0,4],Action(F(13,100)))
segment=[]
for a in [F(0),F(1,3),F(1)]:
qa=[1-a*(1-x) for x in q];h=[1-x for x in q];lhs=source.phi(qa);D=source.daughter(h,h)
rhs=[a*p-a*(1-a)*b*dd for p,b,dd in zip(source.phi(q),source.birth,D)];assert lhs==rhs
segment.append(dict(a=a,residual=lhs))
dump('two_site_barriers.json',dict(exposure_q=q,c=exposure.c,residuals=exposure_res,weight=weight.weights,weight_slack=weight_slack,founder_bounds=rows,
population_generator_identity=identity,segment_checks=segment,scope='Common predictable action with a pathwise rate-time budget. Exact source inequalities feed the conventional probability theorem; no independence is assumed under shared feedback.'))
print('Exact exposure, amplitude and mean-drift witnesses replayed against reconstructed sources N=2..8.',flush=True)
from validated_pgf import certify
policy=certify([(F(3,10),F(40)),(F(1,100),F(60))]);assert F(policy['survival_AA_upper'])<F(1,100)
policy['evidence']='Fresh outward-rounded dyadic Taylor enclosure of the fixed two-site source. Degree 11, step 1/20, invariant cube and one-sided Lipschitz bound; not a floating ODE tolerance.'
policy['exposure']=F(58,5);policy['founders']=1;policy['endpoint']='Any living descendant at time 100; eraser withdrawn at 40 while original killing continues.'
eps=F(policy['survival_AA_upper']);policy['multiple_founder_upper']={str(n):1-(1-eps)**n for n in [1,10,100]}
dump('validated_policy_replay.json',policy)
# All custom configurations use their own source, never the fixed witnesses.
configured=configured_source();idx=configured.index[INITIAL_MARK_COUNTS];z=[0]*configured.size;z[idx]=FOUNDERS
configured_records=[]
for placement in ['early','late','spread']:
phases=phases_for(EXPOSURE_BUDGET,ERASER_AMPLITUDE,HORIZON,placement);u=configured.schedule(phases)
configured_records.append(dict(placement=placement,extinction_vector=u,population_survival=population_survival(u,z),
phases=[dict(duration=p.duration,eraser=p.action.eraser) for p in phases]))
dump('configured_schedules.json',dict(states=configured.states,rates=vars(configured.rates),initial_configuration=z,schedules=configured_records,evidence='Numerical PGF for predetermined schedules; founder independence is valid only here.'))
sweep=[]
for B in np.linspace(0,20,21):
row=dict(exposure=float(B))
for placement in ['early','late','spread']:
u=source.schedule(phases_for(F(str(B)),F(29,100),100,placement));row[placement]=1-u[5]
sweep.append(row)
table('schedule_sweep.csv',['exposure','early','late','spread'],[[r[k] for k in ['exposure','early','late','spread']] for r in sweep])
# Endpoints distinguish live-at-100 from eventual survival after withdrawal.
# Long baseline integration from zero approximates the minimal fixed point;
# this is a numerical diagnostic, not a certified root enclosure.
off=source.schedule([Phase(2000,Action())]);terminal_full=np.full(6,.1)
early=phases_for(F(58,5),F(29,100),100,'early');on40=[Phase(40,Action(F(29,100)))]
endpoints={name:1-source.schedule(phases,terminal)[np.array([5,2])] for name,phases,terminal in [
('live_at_100',early,np.zeros(6)),('eventual_background_retained',early,off),('full_withdrawal_at_100',early,terminal_full),('full_withdrawal_at_40',on40,terminal_full)]}
lo=source.taylor_coefficients(Action());hi=source.taylor_coefficients(Action(F(29,100)))
difference=hi[4][2]-lo[4][2];assert difference==-F(49619,600000000) and all(hi[k][2]==lo[k][2] for k in range(4))
rr_low=source.schedule([Phase(1,Action())])[2];rr_high=source.schedule([Phase(1,Action(F(29,100)))])[2]
dump('preparation_and_withdrawal.json',dict(founder_order=['AA','RR'],endpoints=endpoints,baseline_terminal_extinction=off,
RR_fourth_order_extinction_difference=difference,RR_t1_numeric_extinction_difference=rr_high-rr_low,
scope='The negative exact fourth-order coefficient disproves universal erasure monotonicity near time zero. t=1 and withdrawal curves are numerical. Full withdrawal stipulates all-state death .01 and division .1, with extinction vector .1.'))
# Protected-division uncertainty: affine endpoint checks establish each interval.
robust=[]
for row in json.loads((here/'source_certificates.json').read_text()):
results=[]
for b in map(F,row['birth_interval']):
ss=InheritedPopulation(MolecularRates(sites=row['N'],protected_division=b))
ec=ExposureCertificate(tuple(map(F,row['q'])),F(row['c']));wc=WeightedDrift(tuple(map(F,row['w'])),F(row['gamma']),Action(F(29,100)))
results.append(dict(protected_division=b,exposure=ec.verify(ss),weighted_slack=wc.verify(ss)))
robust.append(dict(N=row['N'],interval=row['birth_interval'],checks=results,scope='Same witness at both endpoints; affine dependence extends inequalities through this fixed-parameter interval.'))
dump('demographic_sensitivity.json',robust)
added=[]
for N in range(2,9):
ss=InheritedPopulation(MolecularRates(sites=N));drift=WeightedDrift((F(1),)*ss.size,F(1,5),Action(F(0),F(29,100)));zz=[0]*ss.size;zz[-1]=FOUNDERS
added.append(dict(N=N,slack=drift.verify(ss),**drift.sufficient(ss,zz,TARGET_RISK)))
# .025 repairs N=6 at e=.3; exact positive weight proposed numerically.
ss=InheritedPopulation(MolecularRates(sites=6));act=Action(F(29,100),F(1,40));A=ss.mean_matrix(act)
eig,vec=np.linalg.eig(np.array(A,dtype=float));v=np.real(vec[:,np.argmax(eig.real)]);v*=np.sign(v.sum());v/=v.min()
w=tuple(F(str(round(float(x),8))) for x in v);bounds=collatz_bounds(ss,act,w);assert bounds[1]<0
configured_gamma=min(d+F(ADDED_DEATH_AMPLITUDE)*p-b for d,p,b in zip(configured.death,configured.protected,configured.birth))
configured_second=dict(amplitude=F(ADDED_DEATH_AMPLITUDE),gamma=configured_gamma,status='unit_weight_drift_not_negative')
if configured_gamma>0:
configured_weight=WeightedDrift((F(1),)*configured.size,configured_gamma,Action(F(0),F(ADDED_DEATH_AMPLITUDE)))
configured_second.update(status='certified',slack=configured_weight.verify(configured),bound=configured_weight.sufficient(configured,z,TARGET_RISK))
dump('second_actuator.json',dict(configured=configured_second,uniform_policy=added,six_site_repair=dict(action=vars(act),weights=w,collatz=bounds),
scope='Uniform sufficient amplitude is >.09, not an exact N-independent critical threshold. At .29 the sufficient exposure is 1.45*log(n/risk). Costs of distinct actuators must be declared separately.'))
Breq=exposure.necessary_exposure(FOUNDERS,TARGET_RISK,5)
dump('duration_and_units.json',dict(finite_time_floor_at_100=finite_time_floor(FOUNDERS,100,.3),
conditional_illustration=concentration_lower_bound(Breq,100,.29,1.,.1,1.),
scope='Illustrative uncalibrated saturating engagement v(C)=vmax*C/(K+C), and one-compartment clearance. AUC/amount outputs are conditional necessary bounds, not sufficient schedules; exposure is dimensionless additional reaction-rate*time. Time and exposure lower bounds are simultaneous, not multiplied.'))
if RUN_FEEDBACK_DEMONSTRATION:
sim=FeedbackSimulator(configured,F(ERASER_AMPLITUDE),F(EXPOSURE_BUDGET),RANDOM_SEED)
def feedback_policy(t,count,spent,history):return float(F(ERASER_AMPLITUDE)) if sum(count[i] for i,p in enumerate(configured.protected) if p)>0 else 0.
dump('feedback_path.json',sim.run(z,float(F(HORIZON)),feedback_policy))
plot(out,rows,sweep,sites,policy)
print(f'Fresh finite-time enclosure: one AA founder, exposure 11.6, survival upper {float(eps):.15f} at time 100.',flush=True)
print('Amplitude .29 leaves certified survival floors at N=5..8; a second actuator restores contraction. Exact and numerical endpoints remain distinct.',flush=True)
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(),
module_sha256={p.name:digest(p) for p in sorted(here.glob('*.py')) if p.name not in ['example.py','test_example.py']},
input_sha256={p.name:digest(p) for p in sorted(here.glob('*certificates.json'))},
output_sha256={p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}))
def plot(out,rows,sweep,sites,policy):
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
axs[0].semilogx([r['founders'] for r in rows],[r['necessary'] for r in rows],'o-',label='Necessary, any feedback')
axs[0].semilogx([r['founders'] for r in rows],[r['eraser_exposure'] for r in rows],'s-',label='Sufficient weighted policy')
axs[0].set(xlabel='AA founder count',ylabel='Reaction exposure',title='Exposure bounds versus founder count');axs[0].legend(fontsize=8)
for name in ['early','late','spread']:axs[1].semilogy([r['exposure'] for r in sweep],[r[name] for r in sweep],label=name.capitalize())
axs[1].scatter([11.6],[float(F(policy['survival_AA_upper']))],marker='*',s=100,color='black',label='Certified upper bound')
axs[1].axhline(.01,color='gray',ls=':');axs[1].set(xlabel='Reaction exposure',ylabel='One AA founder survival at time 100',title='Survival under equal-exposure treatment\nschedules');axs[1].legend(fontsize=8)
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'exposure_and_timing.png',dpi=180);fig.savefig(out/'exposure_and_timing.svg');plt.close(fig)
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
Ns=list(range(2,9));vals=[float(sites['mean_sign_endpoints'][str(n)]['intrinsic_erasure_bracket'][1]) for n in Ns]
axs[0].plot(Ns,vals,'o-',label='Certified mean-sign endpoint');axs[0].axhline(.3,color='black',ls='--',label='Baseline + .29 actuator')
axs[0].set(xlabel='Memory sites',ylabel='Intrinsic erasure rate',title='Mean-growth erasure thresholds versus\nmemory size');axs[0].legend(fontsize=8)
selected=[v for v in sites['amplitude'].values() if v['N']>=5 and v['maximum']==F(29,100)]
axs[1].bar([v['N'] for v in selected],[float(v['fully_marked_floor']) for v in selected])
axs[1].set(xlabel='Memory sites',ylabel='Certified survival floor, one founder',ylim=(0,1),title='All-policy survival lower bounds by memory\nsize')
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'amplitude_floor.png',dpi=180);fig.savefig(out/'amplitude_floor.svg');plt.close(fig)
if __name__=='__main__':main()
Run output
Exact exposure, amplitude and mean-drift witnesses replayed against reconstructed sources N=2..8. Exact determinant, root bracket, certificate slacks, and Taylor checks passed. Fresh finite-time enclosure: one AA founder, exposure 11.6, survival upper 0.009818264214474 at time 100. Amplitude .29 leaves certified survival floors at N=5..8; a second actuator restores contraction. Exact and numerical endpoints remain distinct.