Example code
Treatment order can affect extinction differently in two cell-population models even when all their expected counts agree. Two seven-type branching models share the same molecular reactions, division, death, acquisition of resistance and one-daughter marginals. They differ only in how sister states are paired: complementary allocation conserves the mother's old marks; independent sampling does not conserve them across the pair.
For the declared eventual-extinction objective, fresh interval calculations prove that complementary sisters have greater extinction probability under B then A, while independent sisters favour A then B, across a small certified killing-rate interval. Following the independent approximation when complementary inheritance holds loses about 1.70134 × 10⁻⁶ in one-founder extinction probability. The effect is intentionally small and near a tie.
Reusable sister-law, acquisition, resistant-lineage, action and phase components support new comparisons. The example exposes editable inputs, full mean matrices, distinct survival endpoints, backward calculation of extinction through successive treatment phases and both leading-order and finite-amplitude dependence responses. Numerical exploration stays separate from the fixed reference certificates.
Eight freshly enclosed flows establish the central comparison and its continuum extension. Exact arithmetic also checks the shared means, clearing continuation, covariance-cone conditions, founder maximum at 32 and scalar feedback floor. Separate continuation calculations account for residual acquisition and the substantial extra exposure of finite clearing.
The reference continuation is part of the result: it clears molecular cells without creating new resistant cells. It is not drug-free. Rates and acquisition probabilities are synthetic; the example demonstrates a model-reduction failure, not a treatment recommendation or an optimal schedule. Lean is not rerun.
Python source
"""Same means, different sisters, opposite extinction decisions."""
# ------------------------- EDITABLE INPUTS -------------------------
ACQUISITION_SCALE='1/10'
ACQUISITION_PROFILE=('1/10','1/10','1/10','2/5','1/10','2/5')
ACTION_A_ERASURE='3/10'
ACTION_A_EXTRA_DEATH='0'
ACTION_B_ERASURE='1/100'
ACTION_B_EXTRA_DEATH='120553/1000000'
PULSE_DURATION='10'
RESISTANT_DIVISION='1/10'
RESISTANT_DEATH='1/50'
INITIAL_COUNTS=(0,0,0,0,0,1,0) # UU,UR,RR,AU,AR,AA,M
CORRECTION_ERROR_ALLOWANCE='4/1000000'
OUTPUT_DIRECTORY='outputs'
# Synthetic model rates and per-division probabilities, not measured dosing.
# The fixed reference certificate is kept separate from custom scenarios.
# ------------------------------------------------------------------
import argparse,csv,hashlib,json,platform
from fractions import Fraction as F
from pathlib import Path
import numpy as np
import sympy as sp
from branching import *
from validated import certify
from consequence_kernel import scalar_floor,mean_matrix,linear_flow
MANUSCRIPT_SHA256='96c8eb87625ee95044617d81f9a11ed08d44cf309c63d8036f84ea0775fbf3ff'
def subtract(a,b):return (a[0]-b[1],a[1]-b[0])
def bounds(run):return tuple(map(F,run['bounds'][5]))
def winner(gap):return 'AB' if gap[0]>0 else 'BA' if gap[1]<0 else 'unresolved'
def reference_certificates():
A=(F(3,10),F(0),F(10));B=(F(1,100),F(120553,10**6),F(10))
runs={law+order:certify([A,B] if order=='AB' else [B,A],law) for law in ['J','I'] for order in ['AB','BA']}
q={k:bounds(v) for k,v in runs.items()};gaps={law:subtract(q[law+'AB'],q[law+'BA']) for law in ['J','I']}
assert gaps['J'][1]<-F(16,10**7) and gaps['I'][0]>F(16,10**7)
left,center,right=map(F,['.120545','.120553','.120561']);K=F(6912);boundary={}
for law,end in [('J',left),('I',right)]:
BB=(F(1,100),end,F(10));rr={order:certify([A,BB] if order=='AB' else [BB,A],law) for order in ['AB','BA']}
gap=subtract(bounds(rr['AB']),bounds(rr['BA']))
num=subtract(gaps[law],gap) if law=='J' else subtract(gap,gaps[law]);distance=abs(center-end)
secant=tuple(v/distance for v in num);derivative=(secant[0]-K*(right-left),secant[1]+K*(right-left))
assert derivative[1]<0 and (gap[1]<0 if law=='J' else gap[0]>0)
boundary[law]=dict(endpoint=end,gap=gap,secant=secant,derivative=derivative,runs=rr)
return dict(runs=runs,AA_extinction=q,order_gaps=gaps,preferred={law:winner(gap) for law,gap in gaps.items()},
continuum=dict(interval=(left,right),curvature_bound=K,endpoints=boundary,scope='Four fresh endpoint flows plus the analytic curvature and mean-value argument prove the whole interval. Not a maximal reversal set.'),
evidence='Fresh exact outward dyadic Taylor runs; no saved certificate values loaded. Conventional soundness/probability proofs, not Lean.'),q,gaps
def decision_consequences(q,gaps):
mJ=(-gaps['J'][1],-gaps['J'][0]);mI=gaps['I']
errors={s:subtract(q['I'+s],q['J'+s]) for s in ['AB','BA']};paired=subtract(errors['AB'],errors['BA']);corrected=subtract(gaps['I'],paired)
regret=tuple(mJ[i]*mI[i]/(mJ[i]+mI[i]) for i in range(2))
founders={}
for law,h,l in [('J',q['JBA'],q['JAB']),('I',q['IAB'],q['IBA'])]:
increasing=l[0]**31*(1-l[1])-h[1]**31*(1-h[0]);decreasing=l[1]**32*(1-l[0])-h[0]**32*(1-h[1])
assert increasing>0 and decreasing<0
founders[law]=dict(unique_maximum=32,increment31_lower=increasing,increment32_upper=decreasing,advantage32=(h[0]**32-l[1]**32,h[1]**32-l[0]**32))
symmetric=subtract(gaps['I'],(-F(CORRECTION_ERROR_ALLOWANCE),F(CORRECTION_ERROR_ALLOWANCE)))
return dict(direct_misspecification_cost=mJ,independent_margin=mI,dependence_errors=errors,paired_correction=paired,
corrected_gap=corrected,corrected_decision=winner(corrected),symmetric_corrected_gap=symmetric,symmetric_decision=winner(symmetric),
minimax_regret=regret,founders=founders,
scope='AA founders, deterministic schedules and this continuation. No arbitrary-mixture ranking or global optimization. Minimax randomizes a complete course under a two-model information restriction; it is not the direct mistake cost.')
def cone_check():
states=source.STATES;pre=lambda i,j:states[i][0]<=states[j][0] and states[i][1]>=states[j][1]
covers=[(i,j) for i in range(6) for j in range(6) if i!=j and pre(i,j) and not any(k not in [i,j] and pre(i,k) and pre(k,j) for k in range(6))]
upsets=[[F((mask>>i)&1) for i in range(6)] for mask in range(1,64) if all(not ((mask>>i)&1) or not pre(i,j) or (mask>>j)&1 for i in range(6) for j in range(6))]
assert len(covers)==6 and len(upsets)==7 and all(source.KAPPA[i]<=source.KAPPA[j] for i,j in covers)
values={}
for label,(e,c) in [('A',(F(3,10),F(0))),('B',(F(1,100),F(120553,10**6)))]:
A=mean_matrix(e,c,F(0));vv=[sum((A[j][k]-A[i][k])*u[k] for k in range(6)) for i,j in covers for u in upsets if u[i]==u[j]]
assert len(vv)==32 and min(vv)>=0;values[label]=vv
return dict(covers=covers,upsets=upsets,checks=values,scope='Exact finite cone-invariance checks. The Harris covariance and positive transport conclusions use the conventional argument; an asymptotic sign alone does not prove the finite reversal.')
def continuation_bounds(decisions):
A=sp.Matrix(mean_matrix(F(3,10),F(0),F(0)));f=sp.Matrix([F(1,10)*k for k in source.KAPPA]);v=[F(x) for x in -A.inv()*f]
assert min(v)>0 and all(sum(F(A[i,j])*v[j] for j in range(6))==-F(f[i]) for i in range(6))
a=(F(3,10),F(0),F(10));b=(F(1,100),F(120553,10**6),F(10));w=list(map(F,[14,11,10,84,43,107]));U={};C={}
for label,phases in [('AB',[a,b]),('BA',[b,a])]:
U[label]=tuple(F(4,5)*x for x in linear_flow(phases,v));C[label]=tuple(x/10 for x in linear_flow(phases,w))
mJ=decisions['direct_misspecification_cost'][0];mI=decisions['independent_margin'][0];eta=F(3,10**5)
assert eta*U['BA'][1]<mJ and eta*U['AB'][1]<mI
def exp_lower(x):
s=term=F(1)
for k in range(1,200):term*=x/k;s+=term
return s
finite={}
for H in [151,152]:
lower=exp_lower(F(9,100)*H);errJ=C['BA'][1]/lower;errI=C['AB'][1]/lower
finite[str(H)]=dict(J_error_upper=errJ,I_error_upper=errI,both_preserved=errJ<mJ and errI<mI,common_added_exposure=F(29,100)*H)
assert finite['152']['both_preserved'] and not finite['151']['both_preserved']
return dict(occupation_reward=v,U=U,C=C,residual_acquisition_scale=eta,residual_errors=(eta*U['BA'][1],eta*U['AB'][1]),finite_clearing=finite,
exposure_ratio=F(29,100)*152/F(29,10),scope='Separate one-sided perturbation scenarios. Do not combine them without a new bound. Finite clearing requires autonomous existing M lineages and no immigration; H=152 is sufficient, not optimal. The baseline indefinite clearing is not drug-free and has infinite charged exposure.')
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)
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)
certificate,q,gaps=reference_certificates();dump('fresh_flow_certificates.json',certificate)
decisions=decision_consequences(q,gaps);dump('decision_consequences.json',decisions);dump('cone_invariance.json',cone_check())
print('Eight fresh outward flow enclosures prove opposite orders and the full killing-rate interval.',flush=True)
dump('scalar_feedback_floor.json',scalar_floor());dump('continuation_bounds.json',continuation_bounds(decisions))
acquisition=Acquisition(F(ACQUISITION_SCALE),tuple(map(F,ACQUISITION_PROFILE)));resistant=ResistantLineage(F(RESISTANT_DIVISION),F(RESISTANT_DEATH))
a=Phase(Action(F(ACTION_A_ERASURE),F(ACTION_A_EXTRA_DEATH)),F(PULSE_DURATION));b=Phase(Action(F(ACTION_B_ERASURE),F(ACTION_B_EXTRA_DEATH)),F(PULSE_DURATION))
models={'J':BranchingSource(ComplementarySisters(),acquisition,resistant),'I':BranchingSource(IndependentSisters(),acquisition,resistant)}
if len(INITIAL_COUNTS)!=7 or any(not isinstance(n,int) or n<0 for n in INITIAL_COUNTS):raise ValueError('Seven nonnegative integer founder counts required.')
matched={};configured={};mechanisms={}
for order,phases in [('AB',[a,b]),('BA',[b,a])]:
matched[order]=dict(phase_matrices=[[[str(v) for v in row] for row in models['J'].mean_matrix(p.action)] for p in phases],means={})
for phase in phases:assert models['J'].mean_matrix(phase.action)==models['I'].mean_matrix(phase.action)
for law,model in models.items():
terminal,slack=model.clearing_terminal();ext=model.compose(phases,terminal)
survival=1-np.prod([ext[i]**n for i,n in enumerate(INITIAL_COUNTS)])
endpoints=dict(eventual_total_survival=1-ext[5],total_survival_at_deadline=1-model.compose(phases,[0]*7)[5],
M_present_at_deadline=1-model.compose(phases,[1]*6+[0])[5],M_ever_appeared=1-model.compose(phases,[1]*6+[0],killed_acquisition=True)[5])
configured[law+order]=dict(extinction=ext,configured_population_survival=survival,AA_endpoint_diagnostics=endpoints,clearing_slack=slack)
matched[order]['means'][law]=model.means(phases,INITIAL_COUNTS)
dump('configured_comparison.json',dict(initial_configuration=INITIAL_COUNTS,acquisition_scale=acquisition.scale,acquisition_profile=acquisition.profile,
resistant_rates=vars(resistant),phases={label:dict(erasure=p.action.erasure,extra_death=p.action.extra_molecular_death,duration=p.duration) for label,p in [('A',a),('B',b)]},
results=configured,evidence='Custom numerical scenario. The fixed interval certificates are not silently transferred after editing inputs.'))
dump('matched_means.json',dict(comparisons=matched,scope='Law-to-law equality under each common deterministic schedule; AB and BA need not have equal means, and shared feedback does not preserve this argument.'))
# Mechanism and smooth boundary curves deliberately keep the paper reference.
ar=Phase(Action(F(3,10)));br=Phase(Action(F(1,100),F(120553,10**6)))
for order,phases in [('AB',[ar,br]),('BA',[br,ar])]:
result,t,states=mechanism(phases);mechanisms[order]=result
assert result['identity_residual']<1e-10
table('mechanism_'+order+'.csv',['backward_time','h_AA','epsilon2_delta_AA','finite_error_AA','response_AA','absolute_bound_AA'],zip(t,states[:,17],.01*states[:,23],states[:,11]-states[:,5],states[:,29],states[:,35]))
dump('covariance_mechanism.json',dict(results=mechanisms,evidence='Numerical sensitivities and the exact finite-amplitude response identity; no certified asymptotic remainder is claimed.'))
curve=[]
for cb in np.linspace(.12054,.12057,21):
bb=Phase(Action(F(1,100),F(str(cb))));row=dict(c=cb)
for law,cls in [('J',ComplementarySisters),('I',IndependentSisters)]:
m=BranchingSource(cls());row[law]=m.compose([ar,bb])[5]-m.compose([bb,ar])[5]
curve.append(row)
table('order_gap_curve.csv',['c_B','J_gap','I_gap'],[(r['c'],r['J'],r['I']) for r in curve])
founderrows=[]
for n in range(1,201):
row=[n]
for h,l in [(q['JBA'],q['JAB']),(q['IAB'],q['IBA'])]:row.append(float(sum(h)/2)**n-float(sum(l)/2)**n)
founderrows.append(row)
table('founder_advantages.csv',['AA_founders','J_advantage','I_advantage'],founderrows)
plot(out,curve,founderrows,mechanisms,decisions)
print(f'Certified preferred order: J={certificate["preferred"]["J"]}, I={certificate["preferred"]["I"]}; direct mistake cost {float(decisions["direct_misspecification_cost"][0]):.12g}.',flush=True)
print('Exact mean equality, founder maximum at 32, scalar feedback floor and separately priced continuation bounds checked.',flush=True)
digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest();here=Path(__file__).parent
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']},
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,curve,founders,mechanisms,decisions):
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
for law in ['J','I']:axs[0].plot([r['c'] for r in curve],[1e6*r[law] for r in curve],label=law+(' complementary' if law=='J' else ' independent'))
axs[0].axvspan(.120545,.120561,color='gray',alpha=.15,label='Certified interval');axs[0].axhline(0,color='black',lw=.7)
axs[0].set(xlabel='Additional death rate in B',ylabel='AB minus BA extinction (millionths)',title='Extinction advantage of AB over BA by\nsister law');axs[0].legend(fontsize=8)
for i,law in enumerate(['J','I']):axs[1].plot([r[0] for r in founders],[1e6*r[i+1] for r in founders],label=law)
axs[1].axvline(32,color='gray',ls=':');axs[1].set(xlabel='Identical AA founders',ylabel='Extinction advantage (millionths)',title='Pulse-order extinction advantage versus\nfounder count');axs[1].legend()
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'decisions.png',dpi=180);fig.savefig(out/'decisions.svg');plt.close(fig)
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
for order in ['AB','BA']:
data=np.genfromtxt(out/('mechanism_'+order+'.csv'),delimiter=',',names=True)
axs[0].plot(data['backward_time'],1e6*data['epsilon2_delta_AA'],label=order+' leading order')
axs[0].scatter([20],[1e6*mechanisms[order]['finite_dependence_error'][5]],marker='*',s=90)
axs[0].axvline(10,color='gray',ls=':');axs[0].set(xlabel='Backward time from the deadline',ylabel='Dependence correction (millionths)',title='Sister-dependence correction through each\nschedule');axs[0].legend(fontsize=8)
vals=[sum(decisions['independent_margin'])/2,sum(decisions['paired_correction'])/2,sum(decisions['corrected_gap'])/2]
axs[1].bar(['I order gap','Paired correction','Corrected J gap'],[1e6*float(v) for v in vals]);axs[1].axhline(0,color='black',lw=.7)
axs[1].set(ylabel='Probability difference (millionths)',title='Dependence correction and pulse-order\ndecision margin');axs[1].tick_params(axis='x',labelsize=8)
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'dependence.png',dpi=180);fig.savefig(out/'dependence.svg');plt.close(fig)
if __name__=='__main__':main()
Run output
Eight fresh outward flow enclosures prove opposite orders and the full killing-rate interval. Certified preferred order: J=BA, I=AB; direct mistake cost 1.70134081855e-06. Exact mean equality, founder maximum at 32, scalar feedback floor and separately priced continuation bounds checked.