Example code
Product collected from a microbial preparation may come from fresh conversion, pre-existing product or stored precursors. This example follows mobile product and convertible reserves through two collection windows and a recovery step. It computes the exact minimum fresh input compatible with the declared material and measurement bounds, then constructs a history that attains it.
For the paper's synthetic 10 mL example, the washed record certifies 1.69 µmol of fresh credited inventory, versus 1.52 µmol for the matched unwashed record. At a target of 1.6 µmol, the first is “at or above” and the second remains unresolved. The code retains every active expression and each material transfer, so users can see which premise supports the result.
The small wash advantage is not free removal of background: known mobile product disappears from both the observed output and its explanation. Here the gain comes from attenuating first-window measurement uncertainty. Uncertainty in stored reserves and loss of productive activity can consume all of that margin.
Editable inputs control stock, reserve or uptake premises, retention, recovery inputs, measured amounts and error bounds. Reusable classes validate finite histories, reproduce zero-fresh alternatives, propagate amount errors, and return at-or-above, below, unresolved or incompatible outcomes. Exact witnesses also show why low collected output supplies no upper bound on fresh formation.
The certificate concerns entry into a declared product-equivalent inventory during the completed windows. It does not establish substrate-specific atom attribution, viability or future function. All numerical inputs are synthetic; the package includes exact parameter sweeps, 3,316 replayed histories and seven scientific test groups, without rerunning Lean.
Python source
"""Sharp fresh-conversion certificate with material histories and design sweeps."""
from fractions import Fraction as F
from pathlib import Path
from dataclasses import asdict
import argparse,csv,hashlib,itertools,json,platform,random
from material import Window,Recovery,History,attaining_history,amount_error
from certificate import DirectReserve,UptakeBound,UnboundedReserve,MaterialCertificate,Observation,report
# EDITABLE INPUTS: micromol product equivalents per ORIGINAL 10 mL aliquot.
ORIGINAL_VOLUME_ML='10'
INITIAL_STOCK_BOUND='10'
PREWASH_RESERVE_BOUND='2'
INITIAL_RESERVE_BOUND='2' # distinct uptake-premise input
CUMULATIVE_UPTAKE_BOUND='0.1' # gross uptake, not net disappearance
MOBILE_RETENTION='0.05'
RESERVE_RETENTION='0.9'
RECOVERY_INPUT_BOUND='0.2'
REPORTED_COLLECTIONS=('6','4')
AMOUNT_ERROR_BOUNDS=('0.2','0.2')
FRESH_TARGET='1.6' # paper's principal target, not legacy 1.5
FRESH_UPPER=None # separately calibrated; never inferred from low output
ADDITIONAL_ACTIVITY_RETAINED='1' # relative to already-sampled unwashed comparator
RESERVE_SLACK='0'
WINDOW_DURATIONS_HOURS=('1','1') # illustrative only: average, not steady-state rate
SEED=62092026
RANDOM_REPLAYS=400
MANUSCRIPT_SHA256='e6e69d291896c7e93c96d47920557c1a1eb16bd3f6c15c83c5697207f75f45df'
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)
def dump(name,data):
(out/name).write_text(json.dumps(data,indent=2,default=lambda x:{'exact':str(x),'decimal':float(x)} if isinstance(x,F) else str(x))+'\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)
B,J,R,T,e,s,H=map(F,(INITIAL_STOCK_BOUND,PREWASH_RESERVE_BOUND,INITIAL_RESERVE_BOUND,CUMULATIVE_UPTAKE_BOUND,MOBILE_RETENTION,RESERVE_RETENTION,RECOVERY_INPUT_BOUND))
observation=Observation(*map(F,REPORTED_COLLECTIONS),*map(F,AMOUNT_ERROR_BOUNDS));threshold=F(FRESH_TARGET)
cert=lambda premise,ee=e:MaterialCertificate(B,ee,s,H,premise)
direct=cert(DirectReserve(J));uptake=cert(UptakeBound(R,T));unbounded=cert(UnboundedReserve())
decisions={name:report(observation,[c],threshold,FRESH_UPPER) for name,c in [('direct',direct),('uptake',uptake),('no-reserve-information',unbounded)]}
decisions['combined']=report(observation,[direct,uptake],threshold,FRESH_UPPER)
witnesses={};q=observation.lower_collections()
if q is not None:
for name,c in [('direct',direct),('uptake',uptake),('no-reserve-information',unbounded)]:
h=c.witness(*q);witnesses[name]=dict(history=asdict(h),fresh=h.fresh,valid_for_premise=c.valid_for(h),events=h.ledger())
# Attainment from any split applies to the DIRECT reserve class.
splits=[direct.witness(*q,v).fresh for v in [F(0),B/2,B]]
unbounded_fresh=direct.witness(*q).extra_unobserved_formation(100)
witnesses['same-collections-extra-fresh-loss']=dict(history=asdict(unbounded_fresh),fresh=unbounded_fresh.fresh,
valid_for_premise=direct.valid_for(unbounded_fresh),note='Low collection does not bound fresh input above: the extra material may be lost.')
else:splits=[]
# Paper reference histories are fixed explicit examples, separate from edits.
first=Window(F(8),F(2),collection=F(6));washed=History(first,Window(F(1,10),F(9,5),fresh=F(21,10),release=F(39,10),collection=F(4)),Recovery())
unwashed=History(first,Window(F(9,5),F(9,5),fresh=F(21,10),release=F(39,10),collection=F(57,10)),Recovery(F(9,10),F(9,10)))
stored=History(Window(F(8),F(2),uptake=F(2),collection=F(6)),Window(F(0),F(18,5),release=F(18,5),collection=F(18,5)),Recovery())
ambiguity=History(Window(F(8),F(2),uptake=F(2),collection=F(29,5)),Window(F(1,100),F(19,5),release=F(19,5),collection=F(381,100)),Recovery(reserve_input=F(1,5)))
paper_cert=lambda j,ee=F(1,20),hh=F(1,5):MaterialCertificate(F(10),ee,F(9,10),hh,DirectReserve(j))
paper_obs=Observation(F(6),F(4),F(1,5),F(1,5))
paper_w=report(paper_obs,[paper_cert(F(2))],F(8,5));paper_u=report(Observation(F(6),F(57,10),F(1,5),F(1,5)),[paper_cert(F(2),F(9,10))],F(8,5))
counter=dict(storage=asdict(stored),true_fresh=stored.fresh,
wrong_initial_substitution=paper_cert(F(2),hh=F(0)).evaluate(6,F(18,5))['lower'],
wrong_premise_rejected=not paper_cert(F(2),hh=F(0)).valid_for(stored),
correct_direct=paper_cert(F(4),hh=F(0)).evaluate(6,F(18,5))['lower'],
correct_uptake=MaterialCertificate(F(10),F(1,20),F(9,10),F(0),UptakeBound(F(2),F(2))).evaluate(6,F(18,5))['lower'],
same_record_zero_fresh_history=asdict(ambiguity),within_output_envelopes=abs(ambiguity.first.collection-6)<=F(1,5) and abs(ambiguity.second.collection-4)<=F(1,5))
reserve_rows=[];activity_rows=[]
for n in range(101):
j=F(n,20);r=paper_cert(j).evaluate(F(29,5),F(19,5))
reserve_rows.append([float(j),float(r['lower']),';'.join(r['active'])])
for delta in [F(0),F(1,20),F(1,10),F(1,5)]:
for n in range(101):
a=F(n,100);second=F(17,10)+F(21,10)*a
r=paper_cert(2+delta).evaluate(F(29,5),second)
assert r['lower']==max(F(0),F(21,10)*a-F(41,100)-F(17,20)*delta)
activity_rows.append([float(a),float(delta),float(r['lower']),r['lower']>F(38,25),r['lower']>=F(8,5)])
a,delta=map(F,(ADDITIONAL_ACTIVITY_RETAINED,RESERVE_SLACK))
if not 0<=a<=1 or delta<0:raise ValueError('Activity must be in [0,1]; reserve slack must be nonnegative.')
design=dict(activity=a,slack=delta,certificate=paper_cert(2+delta).evaluate(F(29,5),F(17,10)+F(21,10)*a)['lower'],
strict_improvement_activity_boundary=F(193,210)+F(17,42)*delta,
target_activity_boundary=F(67,70)+F(17,42)*delta,
margin_budget=F(17,100)-F(17,20)*delta-F(21,10)*(1-a),
reserve_boundaries=dict(positive_strict=F(339,85),target_inclusive=F(179,85),improvement_strict=F(11,5)),
uptake_boundaries=dict(target_inclusive=F(9,85),improvement_strict=F(1,5)))
# Matched source: second collection changes with mobile retention.
matched=[]
for ee in [F(0),F(1,20),F(1,2),F(9,10)]:
Q2=ee*2+F(9,10)*2+F(21,10)
nominal=paper_cert(F(2),ee).evaluate(6,Q2)['lower']
robust=paper_cert(F(2),ee).evaluate(F(29,5),Q2-F(1,5))['lower']
matched.append(dict(e=ee,second_collection=Q2,nominal=nominal,robust=robust))
# Exact endpoint and random schedule replay, not a replacement for the proof.
replay_count=0
for q1,q2,stock,reserve,input_ in itertools.product([F(0),F(1),F(3)],repeat=5):
for ee,ss in [(F(0),F(0)),(F(0),F(1)),(F(1,2),F(1,2)),(F(1),F(1))]:
c=MaterialCertificate(stock,ee,ss,input_,DirectReserve(reserve))
for split in [F(0),stock/2,stock]:
h=c.witness(q1,q2,split);h.ledger();assert c.valid_for(h);replay_count+=1
rng=random.Random(SEED)
for _ in range(RANDOM_REPLAYS):
draw=lambda:F(rng.randrange(101),10)
stock,q1,q2,rb,tb,inp=draw(),draw(),draw(),draw(),draw(),draw();ee=F(rng.randrange(11),10);ss=F(rng.randrange(int(ee*10),11),10)
c=MaterialCertificate(stock,ee,ss,inp,UptakeBound(rb,tb));h=c.witness(q1,q2);h.ledger();assert c.valid_for(h);replay_count+=1
durations=tuple(map(F,WINDOW_DURATIONS_HOURS))
if min(durations)<=0:raise ValueError('Positive recorded window durations required.')
summary=dict(configured_decisions=decisions,all_initial_splits_attain=splits,original_volume_ml=F(ORIGINAL_VOLUME_ML),
total_window_average_lower=None if decisions['direct']['lower'] is None else decisions['direct']['lower']/sum(durations),
paper_reference=dict(washed=paper_w,unwashed=paper_u,washed_history=asdict(washed),unwashed_history=asdict(unwashed),realized_fresh=washed.fresh),
exact_replay_count=replay_count,matched_retention=matched,
measurement_example=dict(volume_ml=F(2),concentration_umol_per_ml=F(3),volume_error=F(1,20),concentration_error=F(1,10),
amount_error=amount_error(2,3,F(1,20),F(1,10))),
interpretation=['Credits fresh entry into declared inventory, which may include intermediates; not substrate-specific atom attribution.',
'No inference of viability, growth, abundance, future function or a steady-state rate.',
'Initial stock must include all old convertible reserves. Challenge substrate must be assigned once, outside before fresh entry or inside old stock.',
'Collection physically leaves the accounting boundary once; wash waste is separate.',
'All inputs are synthetic. Lean is not rerun; exact rational computation checks the source algebra and constructed histories.'])
dump('summary.json',summary);dump('attaining_histories.json',witnesses);dump('counterexamples.json',counter);dump('design.json',design)
table('reserve_sweep.csv',['reserve_ceiling','lower','active_expressions'],reserve_rows)
table('activity_sweep.csv',['activity','reserve_slack','lower','strictly_beats_unwashed','meets_target'],activity_rows)
plot(out,reserve_rows,activity_rows,washed,ambiguity)
print(f'Configured direct-reserve report: {decisions["direct"]["verdict"]}; lower = {decisions["direct"]["lower"]}')
print(f'{replay_count} exact endpoint/random histories replayed. Paper reference: washed 1.69, unwashed 1.52.',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={n:digest(here/n) for n in ['material.py','certificate.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,reserve,activity,positive,zero):
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
axs[0].plot([r[0] for r in reserve],[r[1] for r in reserve]);axs[0].axhline(1.6,color='black',ls=':',label='Fresh target 1.6');axs[0].axhline(1.52,color='gray',ls='--',label='Unwashed certificate 1.52')
axs[0].set(xlabel='Pre-wash reserve ceiling (micromol equivalents)',ylabel='Certified fresh lower amount',title='Fresh-conversion bound versus initial\nreserve allowance');axs[0].legend(fontsize=8)
for d in [0,.05,.1,.2]:
rows=[r for r in activity if r[1]==d];axs[1].plot([r[0] for r in rows],[r[2] for r in rows],label=f'Reserve slack {d}')
axs[1].axhline(1.52,color='gray',ls='--');axs[1].axhline(1.6,color='black',ls=':');axs[1].set(xlabel='Additional productive activity retained',ylabel='Washed fresh lower amount',title='Fresh-conversion bound versus retained\nactivity');axs[1].legend(fontsize=8)
for ax in axs:ax.grid(alpha=.2)
fig.savefig(out/'design.png',dpi=180);fig.savefig(out/'design.svg');plt.close(fig)
fig,ax=plt.subplots(figsize=(10,4),layout='constrained');x=np.arange(4)
for shift,h,label in [(-.18,positive,'Positive reference history'),(.18,zero,'Zero-fresh alternative')]:
values=[h.first.collection,h.second.collection,h.first.final_reserve,h.fresh]
ax.bar(x+shift,list(map(float,values)),width=.36,label=label)
for i,v in [(0,6),(1,4)]:ax.errorbar(i,v,yerr=.2,fmt='none',color='black',capsize=12,lw=2)
ax.set(xticks=x,xticklabels=['First collection','Second collection','Pre-wash reserve','Fresh total'],ylabel='Micromol equivalents per original aliquot',title='Two material histories compatible with the\ncollection data');ax.legend();ax.grid(axis='y',alpha=.2)
fig.savefig(out/'ambiguity.png',dpi=180);fig.savefig(out/'ambiguity.svg');plt.close(fig)
if __name__=='__main__':main()
Run output
Configured direct-reserve report: at or above; lower = 169/100 3316 exact endpoint/random histories replayed. Paper reference: washed 1.69, unwashed 1.52.