Example code
If no sampled unit shows recovery, how many recoverable units could still have been missed? The answer depends on how reliably each recovery condition activates a unit, produces progeny and records them. This example turns the paper's two-stage, two-condition model into a reusable source and planning tool, while keeping non-detection separate from irreversible death.
The model treats activation, subsequent progeny production and recording as successive probabilistic steps. The example retains failed and unrecorded paths, heterogeneous recoverable types, and an allowed fraction of undetectable units. A stage may work under one condition and another stage under the other, yet neither condition completes the full detection sequence.
The full coverage range is 0–2, because coverage bounds a sum across conditions. This removes the apparent 1/4 ceiling from the restricted model. The code computes the sharp detection floor and an explicit worst-case source for every allocation weight. Equal allocation is optimal, but a nonbinding mismatch constraint can produce a genuine plateau of equally good designs.
For the paper's synthetic robust assumptions, exact rational powers confirm 2,300 independent units, split equally, as the minimum balanced design for excluding a recoverable fraction of 1% at 5% total error. Strengthening coverage to 1.5 reduces the balanced requirement to 750 units. More sampling cannot eliminate the declared calibration-failure floor.
Composable classes handle latent types, uniform or mean-square mismatch, source paths, fixed versus randomized allocation, adaptive negative histories, exact integer planning, and record evaluation. Positive or unevaluable records do not become an all-negative certificate. The default output keeps the calculation conditional because the synthetic inputs do not establish an actual assay's calibration.
The package includes editable inputs, source, scientific checks, planning tables, allocation witnesses, and a finite-follow-up example where different recoverable fractions produce the same observed event law. The resulting exclusion is conditional on the stated population, reference task, sampling and calibration contract. Lean is not rerun, and no biological death claim follows from the demonstration.
Python source
"""Paired recovery: conditional coverage, literal source paths, and finite decisions."""
# EDITABLE SYNTHETIC INPUTS: these are premises, not empirically fitted assay values.
COVERAGE='19/20' # sum across two conditions; natural range [0,2]
MISMATCH='1/2' # uniform or covered-population RMS bound
RECORDING_FLOOR='9/10'
EXCEPTIONAL_MASS='1/20' # among reference-recoverable types
CALIBRATION_FAILURE='1/100' # simultaneous conditional-environment allowance
TARGET_FRACTION='1/100'
TOTAL_ERROR='1/20'
OBSERVED_PER_CONDITION=(1150,1150)
POSITIVE_RECORDS=0
UNEVALUABLE_RECORDS=0
REFERENCE_TASK='Prespecified reference recovery task (synthetic placeholder)'
ELIGIBLE_POPULATION='Independently sampled eligible units after exposure (synthetic placeholder)'
EXPOSURE_DURATION=1.
OBSERVATION_DEADLINE=24.
ACTIVATION_DEADLINE=8.
EXACT_POWER_BIT_BUDGET=8_000_000 # refuse huge rational powers; no silent approximate certificate
import argparse
from dataclasses import dataclass,asdict
from fractions import Fraction as F
import csv
import hashlib
import json
import math
from pathlib import Path
import platform
from typing import Protocol
import numpy as np
from mpmath import mp,iv
from scipy.optimize import minimize
mp.dps=70;iv.dps=70
MANUSCRIPT_SHA256='6b957282b3d264241c9e0c6a5636a45e73229a71721ce5fecd9942089b48302b'
def probability(x,name='probability'):
q=F(x)
if not 0<=q<=1:raise ValueError(name+' must lie in [0,1]')
return q
@dataclass(frozen=True)
class TwoStageType:
x:F; y:F; a:F; b:F
recording1:F=F(1)
recording2:F=F(1)
covered:bool=True
def __post_init__(self):
for name in ['x','y','a','b','recording1','recording2']:object.__setattr__(self,name,probability(getattr(self,name),name))
@property
def responses(self):return (self.x*self.y*self.recording1,self.a*self.b*self.recording2)
def terminals(self,p,condition):
"""Non-target, unactivated, no progeny, unrecorded path, recorded positive."""
p=probability(p)
if condition not in (0,1):raise ValueError('Condition must be 0 or 1')
x,y,r=(self.x,self.y,self.recording1) if condition==0 else (self.a,self.b,self.recording2)
values=(1-p,p*(1-x),p*x*(1-y),p*x*y*(1-r),p*x*y*r)
if sum(values)!=1:raise ArithmeticError('Terminal law is not normalized')
return values
@dataclass(frozen=True)
class CoverageContract:
c:F=F(COVERAGE)
mismatch:F=F(MISMATCH)
kappa:F=F(RECORDING_FLOOR)
eta:F=F(EXCEPTIONAL_MASS)
def __post_init__(self):
object.__setattr__(self,'c',F(self.c));object.__setattr__(self,'mismatch',F(self.mismatch))
if not 0<=self.c<=2 or self.mismatch<0:raise ValueError('Require c in [0,2] and nonnegative mismatch')
for name in ['kappa','eta']:object.__setattr__(self,name,probability(getattr(self,name),name))
@property
def h(self):return min(self.c,2-self.c)
@property
def t(self):return min(self.mismatch,self.h)
@property
def floor(self):return (self.c*self.c-self.t*self.t)/4
@property
def recorded_floor(self):return (1-self.eta)*self.kappa*self.floor
def profile(self,w):
"""Uniform per-type mismatch allocation profile; not a fixed-count likelihood."""
w=probability(w);a=abs(2*w-1);c,h,t=self.c,self.h,self.t
if a*c<=h-t:return (c*c-t*t-a*a*c*c)/4
if a*c<=h:return (c*c+h*h-2*h*t-2*a*c*(h-t))/4
return (c*c+h*h-2*a*c*h)/4
def witness(self,w=F(1,2)):
w=probability(w);a=abs(2*w-1);c,h,t=self.c,self.h,self.t
if a*c<=h-t:v,tau=a*c,t
elif a*c<=h:v,tau=h-t,t
else:v,tau=h,F(0)
sigma=-v if w>=F(1,2) else v
x=(c+sigma+tau)/2;y=(c+sigma-tau)/2
witness=TwoStageType(x,y,c-x,c-y)
if w*x*y+(1-w)*(c-x)*(c-y)!=self.profile(w):raise ArithmeticError('Allocation witness mismatch')
return witness
def worst_population(self):
w=self.witness();covered=TwoStageType(w.x,w.y,w.a,w.b,self.kappa,self.kappa)
invisible=TwoStageType(0,0,0,0,0,0,False)
return TypePopulation(((1-self.eta,covered),(self.eta,invisible)))
@dataclass(frozen=True)
class TypePopulation:
types:tuple
def __post_init__(self):
pairs=tuple((probability(w,'type weight'),t) for w,t in self.types)
if not pairs or sum(w for w,t in pairs)!=1:raise ValueError('Recoverable type weights must sum to one')
object.__setattr__(self,'types',pairs)
@property
def responses(self):return tuple(sum(w*t.responses[j] for w,t in self.types) for j in range(2))
def verify(self,contract,mode='rms'):
if mode not in ('rms','uniform'):raise ValueError('Mismatch mode must be rms or uniform')
mass=sum(w for w,t in self.types if t.covered)
if mass<1-contract.eta:return dict(admitted=False,reason='Too much exceptional recoverable mass')
if mass==0:return dict(admitted=contract.eta==1,covered_mass=mass,mean_square=None)
square=F(0)
for w,t in self.types:
if not t.covered or w==0:continue
if t.x+t.a<contract.c or t.y+t.b<contract.c or min(t.recording1,t.recording2)<contract.kappa:
return dict(admitted=False,reason='Covered type violates coverage/recording premise')
if mode=='uniform' and abs(t.x-t.y)>contract.mismatch:return dict(admitted=False,reason='Uniform mismatch exceeded')
square+=w*(t.x-t.y)**2/mass
if square>contract.mismatch**2:return dict(admitted=False,reason='Conditional covered-population RMS budget exceeded')
actual=sum(self.responses)/2
if actual<contract.recorded_floor:raise ArithmeticError('Coverage inequality failed')
return dict(admitted=True,covered_mass=mass,mean_square=square,actual_response=actual,guaranteed_response=contract.recorded_floor)
def pooled_condition1(self):
x=sum(w*t.x for w,t in self.types);y=sum(w*t.y for w,t in self.types);xy=sum(w*t.x*t.y for w,t in self.types)
return dict(mean_activation=x,mean_second_stage=y,complete_path=xy,conditional_progeny=xy/x if x else None)
def terminals(self,p,condition):return tuple(sum(w*t.terminals(p,condition)[j] for w,t in self.types) for j in range(5))
@dataclass(frozen=True)
class Allocation:
n1:int
n2:int
def __post_init__(self):
if any(type(n) is not int or n<0 for n in (self.n1,self.n2)):raise ValueError('Nonnegative integer group sizes required')
@property
def total(self):return self.n1+self.n2
def negative_mass(self,p,population):
p=probability(p);u,v=population.responses
return (1-p*u)**self.n1*(1-p*v)**self.n2
def random_assignment_mass(self,p,population,w=F(1,2)):
p=probability(p);w=probability(w);u,v=population.responses
return (1-p*(w*u+(1-w)*v))**self.total
class AdaptivePolicy(Protocol):
def weight(self,history): ...
class AlternatingPolicy:
def weight(self,history):return F(len(history)%2)
class NegativeHistoryLaw:
"""Exact finite tree over the four negative terminal events, with feedback."""
def __init__(self,population,p,policy):self.population=population;self.p=probability(p);self.policy=policy
def mass(self,n,max_leaves=100000):
if type(n) is not int or n<0 or 8**n>max_leaves:raise ValueError('Exact history enumeration exceeds its explicit budget')
def recurse(history):
if len(history)==n:return F(1)
w=probability(self.policy.weight(history),'policy weight');result=F(0)
for cond,choice in [(0,w),(1,1-w)]:
if not choice:continue
masses=self.population.terminals(self.p,cond)
for terminal,mass in enumerate(masses[:4]):
if mass:result+=choice*mass*recurse(history+((cond,terminal),))
return result
return recurse(())
class ExactBudgetExceeded(RuntimeError):pass
@dataclass(frozen=True)
class NegativePlanner:
g:F
theta:F=F(TARGET_FRACTION)
delta:F=F(CALIBRATION_FAILURE)
alpha:F=F(TOTAL_ERROR)
bit_budget:int=EXACT_POWER_BIT_BUDGET
def __post_init__(self):
for name in ['g','theta','delta','alpha']:object.__setattr__(self,name,probability(getattr(self,name),name))
if not 0<self.alpha<1:raise ValueError('Total error must be strictly between zero and one')
if type(self.bit_budget) is not int or self.bit_budget<1:raise ValueError('Positive exact-arithmetic budget required')
@property
def z(self):return self.theta*self.g
@property
def feasible(self):return (self.z>0 and self.delta<self.alpha) or (self.z==1 and self.delta<=self.alpha)
def total_error(self,n):
if type(n) is not int or n<1:raise ValueError('Positive integer sample size required')
b=1-self.z
if n*max(b.numerator.bit_length(),b.denominator.bit_length())>self.bit_budget:
raise ExactBudgetExceeded('Exact power too large; no exact integer certificate returned')
return self.delta+(1-self.delta)*b**n
def plan(self):
if not self.feasible:return dict(feasible=False,reason='Zero detection floor or calibration allowance blocks finite certification')
if self.z==1:return dict(feasible=True,random_minimum=1,balanced_minimum=2,total_error=float(self.delta),integer_neighbors_exact=True)
B=(self.alpha-self.delta)/(1-self.delta)
n=max(1,int(mp.ceil(mp.log(mp.mpf(B.numerator)/B.denominator)/mp.log1p(-mp.mpf(self.z.numerator)/self.z.denominator))))
while self.total_error(n)>self.alpha:n+=1
while n>1 and self.total_error(n-1)<=self.alpha:n-=1
even=n+n%2
if even>2 and self.total_error(even-2)<=self.alpha:raise ArithmeticError('Balanced design not minimal')
return dict(feasible=True,random_minimum=n,balanced_minimum=even,total_error=float(self.total_error(even)),
error_at_random_predecessor=float(self.total_error(n-1)) if n>1 else None,
error_at_balanced_predecessor=float(self.total_error(even-2)) if even>2 else None,
integer_neighbors_exact=True,all_negative_endpoint=self.endpoint(even))
def endpoint(self,n):
"""Outward interval evaluation; unavailable is different from an endpoint of one."""
if type(n) is not int or n<=0 or self.g==0 or self.delta>=self.alpha:return None
B=(self.alpha-self.delta)/(1-self.delta)
q=lambda v:iv.mpf(v.numerator)/v.denominator
value=-iv.expm1(iv.log(q(B))/n)/q(self.g)
return [max(0.,min(1.,math.nextafter(float(value.a),-math.inf))),max(0.,min(1.,math.nextafter(float(value.b),math.inf)))]
def record(self,allocation,positives=0,unevaluable=0,contract_validated=False):
if any(type(n) is not int or n<0 for n in (positives,unevaluable)) or positives+unevaluable>allocation.total:raise ValueError('Invalid record counts')
if positives or unevaluable:return dict(status='no certificate',confidence_endpoint=1,reason='Positive or unevaluable records are not all-confirmed-negative')
if not contract_validated:return dict(status='conditional calculation only',confidence_endpoint=None,reason='The source/calibration/sampling contract has not been established')
if allocation.n1!=allocation.n2 or allocation.total<2:return dict(status='unavailable',confidence_endpoint=None,reason='This certificate requires a fixed balanced design')
if not self.feasible:return dict(status='unavailable',confidence_endpoint=None,reason='No feasible finite guarantee')
error=self.total_error(allocation.total)
return dict(status='excludes target fraction' if error<=self.alpha else 'insufficient sample',
error_upper=float(error),confidence_endpoint=self.endpoint(allocation.total),
interpretation='Repeated-sampling conditional error guarantee; neither sterility nor a posterior death probability')
def independent_minimum(contract,w):
"""Numerical check directly over the four-variable source constraints."""
c,d,w=float(contract.c),float(contract.mismatch),float(w)
constraints=[{'type':'ineq','fun':lambda q:q[0]+q[2]-c},
{'type':'ineq','fun':lambda q:q[1]+q[3]-c},
{'type':'ineq','fun':lambda q:d-q[0]+q[1]},
{'type':'ineq','fun':lambda q:d+q[0]-q[1]}]
results=[];lo=max(0,c-1);hi=min(1,c);width=min(d,hi-lo)
starts=[[c/2]*4,[1,1,1,1],[lo,lo,1,1]]
for x,y in [(lo,lo+width),(lo+width,lo),(hi,hi-width),(hi-width,hi)]:
starts.append([x,y,max(0,c-x),max(0,c-y)])
for start in starts:
sol=minimize(lambda q:w*q[0]*q[1]+(1-w)*q[2]*q[3],start,bounds=[(0,1)]*4,
constraints=constraints,method='SLSQP',options={'ftol':1e-13,'maxiter':300})
if sol.success and min(f['fun'](sol.x) for f in constraints)>-1e-8:results.append(float(sol.fun))
if not results:raise RuntimeError('Independent optimization failed')
return min(results)
def main():
parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',default='outputs');args=parser.parse_args()
out=Path(args.output);out.mkdir(parents=True,exist_ok=True)
if not all(math.isfinite(v) for v in (EXPOSURE_DURATION,ACTIVATION_DEADLINE,OBSERVATION_DEADLINE)) or EXPOSURE_DURATION<0 or not 0<ACTIVATION_DEADLINE<OBSERVATION_DEADLINE:
raise ValueError('Nonnegative exposure and 0 < activation deadline < observation deadline required')
contract=CoverageContract();population=contract.worst_population();planner=NegativePlanner(contract.recorded_floor)
audit=population.verify(contract);plan=planner.plan()
rows=[]
for c,d,k,e,delta in [('1','1/2','1','0','0'),('19/20','1/2','1','0','0'),
('19/20','1/2','9/10','1/20','0'),('19/20','1/2','9/10','1/20','1/100'),('3/2','1/2','9/10','1/20','1/100')]:
con=CoverageContract(c,d,k,e);p=NegativePlanner(con.recorded_floor,delta=delta).plan()
rows.append(dict(c=c,d=d,kappa=k,eta=e,delta=delta,floor=con.floor,g=con.recorded_floor,**p))
assert [r['balanced_minimum'] for r in rows]==[1598,1836,2148,2300,750]
profile_rows=[];maxgap=0.
for c,d in [('19/20','1/2'),('3/2','1/2'),('3/2','1/10')]:
con=CoverageContract(c,d)
for n in range(41):
w=F(n,40);v=con.profile(w);witness=con.witness(w)
numeric=independent_minimum(con,w) if n%5==0 else None
if numeric is not None:maxgap=max(maxgap,abs(numeric-float(v)))
profile_rows.append([c,d,float(w),float(v),numeric,*map(float,(witness.x,witness.y,witness.a,witness.b))])
if maxgap>1e-8:raise ArithmeticError('Independent optimization disagrees with exact profile')
obstruction=TypePopulation(((F(1,2),TwoStageType(1,0,0,1)),(F(1,2),TwoStageType(0,1,1,0))))
class FeedbackPolicy:
def weight(self,history):return F(3,4) if history and history[-1][1]==1 else F(1,3)
adaptive=NegativeHistoryLaw(population,F(TARGET_FRACTION),FeedbackPolicy()).mass(4)
assert adaptive==(1-F(TARGET_FRACTION)*contract.recorded_floor)**4
# Identified-set witness: two recoverable fractions generate the same finite event CDF.
time=np.linspace(0,OBSERVATION_DEADLINE,121);event=.1*(1-np.exp(-time/6))
selected=event/event[-1]
identified=dict(observable_event_at_deadline=float(event[-1]),compatible_recoverable_interval=[float(event[-1]),1],
example_recoverable_fractions=[.2,.8],conditional_cdf_at_deadline=[float(event[-1]/.2),float(event[-1]/.8)],
scope='Illustrative finite-follow-up law; unseen detection mass is placed after the deadline')
allocation=Allocation(*OBSERVED_PER_CONDITION)
record=planner.record(allocation,POSITIVE_RECORDS,UNEVALUABLE_RECORDS,contract_validated=False)
conditional=planner.record(allocation,POSITIVE_RECORDS,UNEVALUABLE_RECORDS,contract_validated=True)
result=dict(inputs=dict(contract=asdict(contract),theta=planner.theta,delta=planner.delta,alpha=planner.alpha,
reference_task=REFERENCE_TASK,eligible_population=ELIGIBLE_POPULATION,exposure=EXPOSURE_DURATION,activation_deadline=ACTIVATION_DEADLINE,observation_deadline=OBSERVATION_DEADLINE),
configured_plan=plan,sharpness_population=audit,record=record,if_contract_is_established=conditional,
planning_rows=rows,allocation_optimization_gap=maxgap,pooled_trap=obstruction.pooled_condition1(),
adaptive_negative_mass=str(adaptive),identified_set_example=identified,
scope='Exact source algebra and integer decisions; no empirical calibration or Lean rerun')
dump=lambda name,obj:(out/name).write_text(json.dumps(obj,indent=2,default=lambda v:str(v) if isinstance(v,F) else float(v))+'\n',encoding='utf-8')
dump('results.json',result)
def table(name,header,values):
with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(values)
table('allocation.csv',['c','d','weight','exact_profile','independent_numeric_minimum','x','y','a','b'],profile_rows)
table('planning.csv',['c','d','kappa','eta','delta','g','random_minimum','balanced_minimum','total_error'],[[r[k] for k in ['c','d','kappa','eta','delta','g','random_minimum','balanced_minimum','total_error']] for r in rows])
table('finite_followup.csv',['time','observable_event_probability','conditional_CDF_p_0.2','conditional_CDF_p_0.8','selected_positive_timing_CDF'],zip(time,event,event/.2,event/.8,selected))
lines=[f'Configured recorded floor g={contract.recorded_floor}; minimum balanced sample {plan.get("balanced_minimum","unavailable")}.',
'Paper balanced planning rows: '+', '.join(str(r['balanced_minimum']) for r in rows)+'.',
f'Independent four-variable allocation optimization maximum gap: {maxgap:.3g}.',
f'Observed-record status: {record["status"]}; under the explicit contract: {conditional["status"]}.',
'All inputs are synthetic. No negative record is interpreted as irreversible death.']
(out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
plot(out,profile_rows,rows)
digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,python=platform.python_version(),source_sha256=digest(Path(__file__)),
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,profiles,rows):
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
cs=np.linspace(0,2,241)
for d in [0,.25,.5,1]:axs[0].plot(cs,[(c*c-min(d,c,2-c)**2)/4 for c in cs],label=f'Mismatch {d}')
axs[0].axvline(1,color='k',ls=':');axs[0].set(title='Detection floor across the allowed coverage\nrange',xlabel='Coverage c',ylabel='Complete-path floor');axs[0].legend(fontsize=7)
for c,d in [('19/20','1/2'),('3/2','1/2'),('3/2','1/10')]:
data=[r for r in profiles if r[:2]==[c,d]];axs[1].plot([r[2] for r in data],[r[3] for r in data],label=f'c={float(F(c))}, d={float(F(d))}')
checked=[r for r in data if r[4] is not None];axs[1].scatter([r[2] for r in checked],[r[4] for r in checked],s=8,color='k')
axs[1].set(title='Detection floor versus allocation between\nconditions',xlabel='Condition-1 allocation weight',ylabel='Worst-case complete-path response');axs[1].legend(fontsize=7)
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'coverage_allocation.png',dpi=180);fig.savefig(out/'coverage_allocation.svg');plt.close(fig)
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
axs[0].bar(np.arange(5),[r['balanced_minimum'] for r in rows]);axs[0].set_xticks(np.arange(5),['Ideal','Coverage','Recording','Calibration','Full range'],rotation=20);axs[0].set(title='Required sample sizes under five\ncalibration assumptions',ylabel='Independent units, balanced design')
n=np.arange(2,4001,2)
for row,label in [(rows[3],'c=0.95'),(rows[4],'c=1.5')]:
g=float(row['g']);axs[1].plot(n,.01+.99*np.exp(n*np.log1p(-.01*g)),label=label)
axs[1].axhline(.05,color='k',ls=':',label='Total error target');axs[1].axhline(.01,color='grey',ls='--',label='Calibration floor');axs[1].set(yscale='log',ylim=(.008,1),title='False-exclusion bound versus sample size',xlabel='Independent units',ylabel='Worst-case false-certificate probability');axs[1].legend(fontsize=7)
for ax in axs:ax.grid(axis='y',alpha=.2)
fig.savefig(out/'planning.png',dpi=180);fig.savefig(out/'planning.svg');plt.close(fig)
if __name__=='__main__':main()
Run output
Configured recorded floor g=44631/320000; minimum balanced sample 2300. Paper balanced planning rows: 1598, 1836, 2148, 2300, 750. Independent four-variable allocation optimization maximum gap: 5e-16. Observed-record status: conditional calculation only; under the explicit contract: excludes target fraction. All inputs are synthetic. No negative record is interpreted as irreversible death.