Example code

"""Executable assay-design casebook: competing explanations and scoped decisions.

The source models remain independent; only their reporting records are compared.
"""
from fractions import Fraction as F
from pathlib import Path
from dataclasses import asdict
import argparse,csv,hashlib,json,platform
import numpy as np
from enzyme import CapacityFamily
from population import Population,TwoBandClass,ReadoutCalibration,row_span_certificate
from sandwich import BindingSite,SandwichSource,ObservationBudget,paper_certificates
from material import Window,Recovery,History
from material_certificate import MaterialCertificate,DirectReserve,Observation,report
from reporter import CofactorReactor,LinearRate,Reporter,AssociationSchedule
from amplification import AmplificationSource,IdentityChannel,exp_negative
from timing_certificate import timing_dual_certificate
from decisions import Claim,DeterministicEnclosure,IssuanceGuarantee,PredictionRegion,specimen_bound,path_floor,negative_recovery_bound,founder_cdf,thinned_cdf_two

# EDITABLE INPUTS: independent synthetic models; units do not transfer between cases.
RECOVERY_FRACTION_TARGET='2/3'
READOUT_SENSITIVITY=('0.9','1')
READOUT_FALSE_POSITIVE=('0','0.02')
READOUT_POSITIVE_FRACTION=('0.68','0.72')
PER_READING_ABSOLUTE_ERROR='0.009'  # simultaneous bound; strict requirement <.01
DILUTION_FACTOR=10.
SPECIMEN_FRACTION='0.45'
SPECIMEN_COUNT_THRESHOLD=6
NATIVE_REFERENCE_PAIRS=9
FALSE_ISSUANCE_ALLOWANCE='0.05'
PATH_STAGE_COVERAGE='0.8'
PATH_MISMATCH='0.2'
UNCOVERED_TYPE_FRACTION='0.05'
RECORDING_PROBABILITY='0.9'
RECOVERABLE_FRACTION_THRESHOLD='0.1'
CALIBRATION_FAILURE='0.01'
ELIGIBLE_UNITS=400
FRESH_INVENTORY_TARGET='1.6'
MICROBIAL_RESERVE_CEILING='2'
REPORTER_RELEASE='20'
SOURCE_MODELS_ARE_SYNTHETIC=True
MANUSCRIPT_SHA256='648302e05466fa7eb109e8e5299bacb868589480ddbdee61d888cef9754c11c8'


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'))
    out=parser.parse_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)
    family=CapacityFamily();bands=TwoBandClass();bulk=bands.bounds(family)
    if bulk['interval'] is None:raise RuntimeError('Kinetic band classification must precede the population claim.')
    readout=ReadoutCalibration(tuple(map(F,READOUT_SENSITIVITY)),tuple(map(F,READOUT_FALSE_POSITIVE)),tuple(map(F,READOUT_POSITIVE_FRACTION)))
    repaired=readout.interval(bulk['interval'])
    uniform=Population((F(1),),(F(1),));heterogeneous=Population((F(9,14),F(11,8)),(F(21,41),F(20,41)))
    mean_rows=[]
    for N,S,H in [(F(28),F(7),F(28)),(F(10),F(3),F(46)),(F(50),F(12),F(6))]:
        a,b=[p.pooled_rate(family,(N,S,H)) for p in (uniform,heterogeneous)]
        assert a==b;mean_rows.append(dict(inputs=(N,S,H),uniform=a,heterogeneous=b))
    claim=Claim('Fraction reaching carrier state 28 from 10','equally weighted units','fixed scalar source, deadline 120, separated capacity bands',F(RECOVERY_FRACTION_TARGET))
    population=dict(claim=asdict(claim),kinetic_contract=bands.classify_bands(family),same_pooled_records=mean_rows,
        populations=[dict(capacities=p.capacities,weights=p.weights,mean=p.mean,recovery=p.recovery_enclosure(family)) for p in [uniform,heterogeneous]],
        pooled_interval=bulk['interval'],pooled_decision=DeterministicEnclosure(*bulk['interval'],True).decide(claim),
        repaired_interval=repaired,repaired_decision=DeterministicEnclosure(*(repaired or (None,None)),feasible_witness=repaired is not None,proved_empty=repaired is None).decide(claim),
        repaired_witnesses=[] if repaired is None else [dict(population=asdict(bands.witness(p)),readout=readout.witness(p)) for p in repaired],
        feature_identifiability=row_span_certificate([[1,1,1],[F(9,14),1,F(11,8)]],[0,1,1]),
        evidence='Fresh exact kinetic band checks, population moments and endpoint witnesses; no whole-cell viability conclusion.')
    source=SandwichSource(BindingSite(1,1),BindingSite(1,1));cert=paper_certificates()['main'];proof=cert.check()
    budget=ObservationBudget(F(PER_READING_ABSOLUTE_ERROR),F(PER_READING_ABSOLUTE_ERROR))
    records=[]
    for x,rho in [(.08,1.),(50.,1.),(.1,1.),(100.,.001)]:
        neat,diluted=source.paired(x,rho,DILUTION_FACTOR);z=F(str(float(diluted-neat)))
        # The masking counterexample is outside the positive native-availability contract.
        calibrated=(.8<=rho<=1 and 9<=DILUTION_FACTOR<=11)
        records.append(dict(native=x,availability=rho,neat=float(neat),diluted=float(diluted),contrast=float(z),
            within_native_availability_contract=calibrated,without_promise=budget.report(z,cert,False,calibrated),with_promise=budget.report(z,cert,True,calibrated)))
    masked=[]
    for spike in [0.,.5,2.]:
        a=source.signal((.1+spike)/DILUTION_FACTOR);b=source.signal((100*.001+spike)/DILUTION_FACTOR)
        masked.append(dict(fully_accessible_spike=spike,low_record=float(a),masked_high_record=float(b)))
    sandwich=dict(continuum_certificate=proof,margin=cert.margin,error_budget=asdict(budget),contrast_ranges=budget.bounds(cert.margin),
        strict_separation=2*(budget.neat_error+budget.diluted_error)<budget.gain_min*cert.margin,
        records=records,spike_ambiguity=masked,
        scope='Fresh-reagent pre-dilution; independent sites, native availability [.8,1], declared high range [25,100]. Outside a two-class promise, only class exclusions follow.')
    g,z=specimen_bound(F(SPECIMEN_FRACTION),SPECIMEN_COUNT_THRESHOLD,NATIVE_REFERENCE_PAIRS)
    rule=IssuanceGuarantee(g,F(FALSE_ISSUANCE_ALLOWANCE),'Independent native-reference pairs, each input containing exactly one target, followed by an independent future specimen with shared conditional recovery state.')
    # Sharp adversary: (X,Y)=(0,0) with 1-z, (1,0) with z.
    H=(1-F(SPECIMEN_FRACTION))**SPECIMEN_COUNT_THRESHOLD;conditional_miss=1-(1-H)*z
    specimen=dict(bound=g,gate_probability_adversary=z**NATIVE_REFERENCE_PAIRS,negative_given_gate=conditional_miss,
        joint_at_adversary=z**NATIVE_REFERENCE_PAIRS*conditional_miss,blank_availability_synthetic=F(24,25)**NATIVE_REFERENCE_PAIRS,
        complete_record=rule.apply(True,True),missing_record=rule.apply(True,True,False),reference_inputs=2*NATIVE_REFERENCE_PAIRS,
        target_preparations=2,claim=f'fewer than {SPECIMEN_COUNT_THRESHOLD} original native targets')
    floor=path_floor(F(PATH_STAGE_COVERAGE),F(PATH_MISMATCH));gpath=(1-F(UNCOVERED_TYPE_FRACTION))*F(RECORDING_PROBABILITY)*floor['floor']
    paths=dict(**floor,recorded_floor=gpath,all_negative_bound=negative_recovery_bound(ELIGIBLE_UNITS,F(RECOVERABLE_FRACTION_THRESHOLD),gpath,F(CALIBRATION_FAILURE)),
        stage_mismatch_zero_path=dict(stages=(1,0,0,1),complete_paths=(0,0)),
        full_range=[dict(coverage=c,**path_floor(c,2)) for c in [F(0),F(1,2),F(1),F(3,2),F(2)]],
        scope='Type-level stage and representation assumptions; conditional calibration validity; independent eligible units and equal allocation.')
    predictions=[PredictionRegion(k,founder_cdf(k,shared),F(19,20),'shared' if shared else 'independent').evaluate() for shared in [False,True] for k in [5,6,7]]
    amplification=AmplificationSource();b4=amplification.certify(4);b5=amplification.certify(5);b7=amplification.certify(7)
    dual=timing_dual_certificate();intercept=b5.miss[0]+5*b4.blank[0]
    assert intercept>F(1055,10000)
    amplification_case=dict(timing_dual=dual,dual_intercept_lower=intercept,miss_lower_at_blank_001=intercept-F(5,100),
        identity_at_seven=IdentityChannel().guaranteed_errors(b7),loading_floor_enclosure=tuple(v*F(99,100) for v in exp_negative(4)),
        scope='Finite h=capacity=5 source. Sign cells and scalar enclosures are recomputed; the source likelihood-ratio argument is imported from the companion. Poisson-empty loaded reactions remain in the miss denominator.')
    mc=MaterialCertificate(F(10),F(1,20),F(9,10),F(1,5),DirectReserve(F(MICROBIAL_RESERVE_CEILING)))
    mo=Observation(F(6),F(4),F(1,5),F(1,5));washed=report(mo,[mc],F(FRESH_INVENTORY_TARGET))
    uc=MaterialCertificate(F(10),F(9,10),F(9,10),F(1,5),DirectReserve(F(MICROBIAL_RESERVE_CEILING)))
    unwashed=report(Observation(F(6),F(57,10),F(1,5),F(1,5)),[uc],F(FRESH_INVENTORY_TARGET))
    # This perspective's Appendix A.6 uses the allowed recovery input .2 and
    # collections at the UPPER edges 4.2/5.9; unlike the technical paper's 4/5.7.
    first=Window(F(8),F(2),collection=F(6));histories=[]
    for e in [F(1,20),F(9,10)]:
        h=History(first,Window(2*e,F(2),fresh=F(21,10),release=F(41,10),collection=2*e+F(41,10)),Recovery(e,F(9,10),reserve_input=F(1,5)))
        histories.append(dict(history=asdict(h),fresh=h.fresh,ledger=h.ledger()))
    microbial=dict(washed=washed,unwashed=unwashed,attaining=asdict(mc.witness(*mo.lower_collections())),matched_upper_edge_histories=histories)
    r=CofactorReactor(LinearRate(1),LinearRate(1),Reporter(release=float(F(REPORTER_RELEASE))))
    t,Y,account=r.run(AssociationSchedule(),30.)
    reporters=dict(configured=account,storage_factors={c:1+1/F(c) for c in ['.2','1','5','20']},
        cost_coefficients={c:F(1,2)*(1+1/F(c)) for c in ['.2','1','5','20']},
        scope='Explicit bound complex, matched preparation, finite association exposure; native-specific calibration excludes hidden background use.')
    table('reporter_trajectory.csv',['time',*r.columns],zip(t,*Y.T))
    decisions_demo=[]
    for name,enclosure in [('supported',DeterministicEnclosure(F(7,10),F(4,5),True)),('excluded',DeterministicEnclosure(F(0),F(1,2),True)),
        ('unresolved',DeterministicEnclosure(F(1,2),F(1),True)),('missing',DeterministicEnclosure(None,None,evaluable=False)),
        ('incompatible',DeterministicEnclosure(F(4,5),F(1,2))),('nonempty-outer-without-witness',DeterministicEnclosure(F(7,10),F(4,5)))]:
        decisions_demo.append(dict(record=name,result=enclosure.decide(claim)))
    thinning=[dict(detection=z,independent=thinned_cdf_two(z),shared=thinned_cdf_two(z,True),difference=thinned_cdf_two(z)-thinned_cdf_two(z,True)) for z in [F(1),F(3,4),F(1,2),F(1,4)]]
    cases=dict(enzyme=population,sandwich=sandwich,specimen=specimen,recovery_paths=paths,small_population=dict(predictions=predictions,recorded_count_at_most_two=thinning,scope='Independent thinning changes the recorded endpoint. Equality of one statistic at detection .5 does not establish equality of full laws.'),amplification=amplification_case,microbial=microbial,reporter=reporters)
    for name,data in cases.items():dump(name+'.json',data)
    workflow=dict(steps=['Write the intended claim.','Fix the original unit, denominator, task and deadline.','Trace material and observations through preparation and measurement.',
        'Construct another explanation that fits the record but changes the decision.','Choose an observation or intervention and budget its uncertainty, material, time and disturbance.',
        'Validate the complete prespecified rule in its intended setting, retaining unresolved and unevaluable outcomes.'],
        future_records=decisions_demo,denominator=len(decisions_demo),
        evidence_scope='This casebook executes mathematical source contracts; it performs no biological validation. Its different guarantee types cannot be substituted for one another.')
    dump('workflow.json',workflow);plot(out,population,source,microbial)
    print('Eight independent source cases evaluated; all reporting outcomes retained.')
    print(f'Population: {population["pooled_decision"]} -> {population["repaired_decision"]}; microbial: {unwashed["verdict"]} -> {washed["verdict"]}.',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,population,source,microbial):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    for y,iv,label in [(1,population['pooled_interval'],'Pooled activity'),(0,population['repaired_interval'],'Functional readout')]:
        if iv:axs[0].plot(list(map(float,iv)),[y,y],lw=7,label=label)
    axs[0].axvline(float(F(RECOVERY_FRACTION_TARGET)),color='black',ls=':');axs[0].set(yticks=[0,1],yticklabels=['With readout','Pooled only'],xlabel='Compatible task-completion fraction',ylim=(-.4,1.4),title='Recovery-fraction bounds with and without a\nfunctional readout')
    u=np.geomspace(.01,1000,400);axs[1].semilogx(u,source.signal(u),label='Neat');axs[1].semilogx(u,source.signal(u/DILUTION_FACTOR),label='Diluted with fresh reagents');axs[1].scatter([.08,50],source.signal(np.array([.08,50])),color='black')
    axs[1].set(xlabel='Accessible analyte concentration',ylabel='Sandwich signal',title='Dilution and native availability');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'observations.png',dpi=180);fig.savefig(out/'observations.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained');x=np.arange(2)
    axs[0].bar(x-.18,[4,5.7],.36,label='Second reading');axs[0].bar(x+.18,[float(microbial['washed']['lower']),float(microbial['unwashed']['lower'])],.36,label='Fresh-entry lower bound')
    axs[0].axhline(float(F(FRESH_INVENTORY_TARGET)),color='black',ls=':');axs[0].set(xticks=x,xticklabels=['Washed','Unwashed'],ylabel='Equivalents per original aliquot',title='Collected output and fresh-conversion\nbounds after washing');axs[0].legend(fontsize=8)
    cs=np.geomspace(.1,30,200);axs[1].semilogx(cs,.5*(1+1/cs),label='Explicit reporter storage');axs[1].axhline(.5,color='gray',ls='--',label='Instantaneous turnover')
    axs[1].set(xlabel='Reporter catalytic release rate',ylabel='Native loss / reporter product',title='Native-output loss per reporter product');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'interventions.png',dpi=180);fig.savefig(out/'interventions.svg');plt.close(fig)


if __name__=='__main__':main()