Example code

"""Treatment order is a family decision; a retained branch omits sister coupling."""
from fractions import Fraction as F

# EDITABLE INPUTS. Synthetic paper rates, scaled to founder division hazard one.
EXP_MINUS_PULSE=F(7,8)       # exact x=exp(-h); use 49/50 for the short-pulse case
FOUNDER_COVARIANCE=F(1,5)
DESCENDANT_DIVISION=F(1,100)
USE_PERSISTENT_MARGINALS=False
DESCENDANT_COVARIANCES=(F(-1,10),F(1,10))
INITIAL_COUNTS=(1,0,0)       # P,S,T; independent founders conditional on this source
MARKER_FLIP_PROBABILITY=F(1,10)
ACCEPTED_DIVISIONS=800
OBSERVED_AGREEMENTS=600      # illustrative supplied data, not an experimental claim
ERROR_PROBABILITY=F(1,20)
FOUNDER_MEAN_WAIT_HOURS=24   # illustrative scaling only
SEED=6801
MANUSCRIPT_SHA256='fdce95c103c6e2c3f7570d8993029d8e24e331f8a011a269ab6512db853bada3'

import argparse,csv,hashlib,json,platform
from pathlib import Path
from math import log,ceil
import numpy as np
from branching import SisterKernel,FamilySource,RetainedBranchExperiment,Phase,A,B
from certificates import ReferenceCertificate,boundary,pulse_extinction,constant_extinction
from assay import SymmetricReadout,PairedDecisionRule,tail,binomial_confidence_interval,project_covariance,calibrated_free_pulse_choice,covariance_sign,minimax_regret


def enclosure(value,grid=10**15):
    v=F(value);return (F(v.numerator*grid//v.denominator,grid),F(ceil(v*grid),grid))


def symbolic_checks():
    import sympy as s
    x,c=s.symbols('x c');aa=(x**4+x**3+x*x+x+1)*(3*x**5+6*x**4+9*x**3+12*x*x+8*x+4)
    qq=3*x**8+12*x**7+30*x**6+60*x**5+98*x**4+137*x**3+170*x*x+159*x+66
    def bb(r,t):return (x**(t+1)-x**(r+t))/(r-1)+(x*x-x**(t+1))/(t-1)
    def ext(r,t):return 1-x*x-bb(r[0],t[0])-bb(r[1],t[1])+(s.Rational(1,4)+c)*(bb(2*r[0],2*t[0])+bb(2*r[1],2*t[1]))+(s.Rational(1,2)-2*c)*bb(sum(r),sum(t))
    factor=x*x*(x-1)**3*(4*c*aa+(x-1)*qq)/210
    assert s.cancel(ext((2,4),(3,3))-ext((3,3),(2,4))-factor)==0
    cp=(1-x)*qq/(4*aa);P=s.Poly(s.cancel(-s.diff(cp,x)*4*aa**2),x)
    assert all(a>0 for a in P.all_coeffs())
    assert boundary(F(80120,100000))[0]>F(1,4)>boundary(F(80121,100000))[0]
    assert aa.subs(x,1)==210 and qq.subs(x,1)==735
    return dict(factorization=str(factor),positive_derivative_numerator=[int(a) for a in P.all_coeffs()],unique_boundary_x=(F(80120,100000),F(80121,100000)),
        interpretation='cp decreases with x and increases with h=-log(x). For LONGER pulses beyond the crossing AB wins for all feasible c; for shorter pulses both preferences exist.',
        source_note='The manuscript Figure 2 caption reverses its below/above-h0 description. This example follows Proposition 3 and the directly checked monotonicity.')


def main():
    parser=argparse.ArgumentParser();parser.add_argument('--output',type=Path,default=Path('outputs'));args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
    def default(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 dump(name,value):(out/name).write_text(json.dumps(value,default=default,indent=2)+'\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)
    ref=ReferenceCertificate(EXP_MINUS_PULSE,DESCENDANT_DIVISION,arbitrary_descendants=USE_PERSISTENT_MARGINALS)
    h=-log(float(EXP_MINUS_PULSE));kernel=SisterKernel(F(1,2),FOUNDER_COVARIANCE)
    descendants=(SisterKernel(F(3,5),DESCENDANT_COVARIANCES[0]),SisterKernel(F(2,5),DESCENDANT_COVARIANCES[1])) if USE_PERSISTENT_MARGINALS else None
    model=FamilySource(kernel,DESCENDANT_DIVISION,descendants);phases={'AB':[Phase(A,h),Phase(B,h)],'BA':[Phase(B,h),Phase(A,h)]}
    configured={name:dict(extinction_vector=model.extinction(seq),configuration_extinction=model.founder_configuration_extinction(seq,INITIAL_COUNTS)) for name,seq in phases.items()}
    band=ref.band('pulses');gap=band.low_minus_high_extinction_interval(FOUNDER_COVARIANCE)
    dump('configured_source.json',dict(x=EXP_MINUS_PULSE,pulse=h,epsilon=DESCENDANT_DIVISION,founder_kernel=kernel.entries,descendant_kernels=[k.entries for k in model.descendants],initial_counts=INITIAL_COUNTS,
        results=configured,certified_single_P_gap=gap,certified_choice=band.decide(FOUNDER_COVARIANCE),pulse_hours=h*FOUNDER_MEAN_WAIT_HOURS,
        scope='Exact reserve certifies one P founder under the fixed A/B response law. Configuration products and ODE values are numerical diagnostics. Independent P replication preserves the rank but changes its magnitude.'))
    dump('symbolic_checks.json',symbolic_checks())
    reference=[];reserve_outputs={}
    for x,cw,e in [(F(7,8),F(1,5),F(1,30)),(F(7,8),F(1,4),F(1,16)),(F(49,50),F(1,4),F(1))]:
        r=ReferenceCertificate(x,e,arbitrary_descendants=True);bb=r.band('pulses');left=bb.low_minus_high_extinction_interval(-cw);right=bb.low_minus_high_extinction_interval(cw)
        assert left[0]>0>right[1]
        reference.append(dict(x=x,covariance_magnitude=cw,epsilon_max=e,negative_covariance_gap=left,positive_covariance_gap=right,minimum_gap=min(left[0],-right[1]),certified_epsilon_supremum=r.reversal_range(cw)))
        reserve_outputs[str(x)]=dict(log_interval=r.h,polynomials=r.polynomials,sensitivity_intervals=r.W,occupation=r.V)
    dump('recurring_division_certificates.json',dict(cases=reference,reserves=reserve_outputs,scope='Fresh rational/log enclosures plus the manuscript comparison theorem. Arbitrary type-dependent descendant kernels satisfy the checked supersolution criterion at these pulse lengths.'))
    # Exact record-law signatures are the evidence; matched sample paths illustrate.
    low=FamilySource(SisterKernel(F(1,2),F(-1,5)),F(1,30));high=FamilySource(SisterKernel(F(1,2),F(1,5)),F(1,30))
    sig=low.observation_signature([A,B]);assert sig==high.observation_signature([A,B])
    means={a.label:low.mean_operator(a) for a in [A,B]};assert all(low.mean_operator(a)==high.mean_operator(a) for a in [A,B])
    longphases=[Phase(A,1),Phase(B,1)];records=[]
    for seed in range(SEED,SEED+40):
        a=RetainedBranchExperiment(low).sample(longphases,seed);b=RetainedBranchExperiment(high).sample(longphases,seed);assert a==b;records.append(a)
    dump('observation_equivalence.json',dict(signature=sig,mean_operators=means,scope='Common life-event clocks and retained-state kernels give exact equality of complete pruned history laws, including under shared retained-history policies. Equality of forty coupled numerical records is an illustration, not the proof. Hidden sister/family observations, survivor selection and molecular-conservation constraints are excluded.'))
    dump('retained_records.json',records)
    # Fixed finite-sample rules: distinct tasks have distinct thresholds/orientations.
    assayref=ReferenceCertificate(F(7,8),F(1,100));rules={};resolution=[]
    for name,task,n,expected in [('constant200','constant',200,(106,137)),('pulse800','pulses',800,(496,584)),('pulse1600','pulses',1600,(1009,1153))]:
        rule=PairedDecisionRule(assayref.assay_band(task),n);assert (rule.low_cut,rule.high_cut)==expected
        probabilities={str(c):{key:enclosure(v) for key,v in rule.probabilities(c).items()} for c in [F(-1,5),F(1,5)]}
        assert rule.probabilities(F(1,5))['resolution']>F(19,20)
        rules[name]=dict(n=n,low_action=rule.band.low_action,high_action=rule.band.high_action,low_cut=rule.low_cut,high_cut=rule.high_cut,p_minus=rule.p_minus,p_plus=rule.p_plus,probability_enclosures=probabilities,
            boundary_error_enclosures=(enclosure(tail(n,rule.p_minus,rule.low_cut)),enclosure(tail(n,rule.p_plus,rule.high_cut,upper=True))))
    constant=ReferenceCertificate(F(7,8),F(1,4)).band('constant');slope=constant.slope;t=constant.threshold;read=SymmetricReadout()
    separated=(tail(40,read.agreement(t-F(1,10)),25,upper=True),tail(40,read.agreement(t+F(1,10)),24))
    assert max(separated)<F(1,20) and constant.decide(t-F(1,10))=='B' and constant.decide(t+F(1,10))=='A'
    dump('paired_rules.json',dict(full_class=rules,separated40=dict(choose_A_at_least=25,error_intervals=[enclosure(v) for v in separated],epsilon_max=F(1,4),assumed_covariance_separation=F(1,10)),
        scope='Independent accepted founder divisions; known symmetric conditionally independent read errors. Full-class error is unconditional and permits unresolved output. The forty-division rule assumes separation; it cannot certify separation from the data.'))
    current=PairedDecisionRule(ref.assay_band('pulses'),ACCEPTED_DIVISIONS,SymmetricReadout(MARKER_FLIP_PROBABILITY),ERROR_PROBABILITY)
    bitinterval=binomial_confidence_interval(ACCEPTED_DIVISIONS,OBSERVED_AGREEMENTS,ERROR_PROBABILITY,bits=28)
    k=current.readout.attenuation;covinterval=project_covariance(bitinterval,(k,k));Dinterval=band.low_minus_high_extinction_interval(*covinterval);Rinterval=(-Dinterval[1],-Dinterval[0])
    dump('configured_assay.json',dict(accepted_divisions=ACCEPTED_DIVISIONS,observed_agreements=OBSERVED_AGREEMENTS,decision=current.decide(OBSERVED_AGREEMENTS),cutoffs=(current.low_cut,current.high_cut),
        bit_confidence_interval=bitinterval,covariance_interval=covinterval,confidence_set_choice=band.decide(*covinterval),extinction_gap_interval=Dinterval,survival_risk_difference_interval=Rinterval,minimax=minimax_regret(Rinterval),
        scope='Illustrative input counts. Direct tail rule uses a conservatively rounded reserve; the confidence-set route uses outward bisection. Their cutoffs need not coincide at machine-close boundaries.'))
    aliases=[(F(1,20),F(0)),(F(1,5),F(1,4))];assert all(SymmetricReadout(eta).agreement(c)==F(3,5) for c,eta in aliases)
    dump('calibration_and_cost.json',dict(same_bit_law_different_decisions=[dict(c=c,eta=eta,agreement=F(3,5),constant_choice=assayref.band('constant').decide(c)) for c,eta in aliases],
        uncertain_calibration_projection=project_covariance((F(59,100),F(61,100)),(F(1,4),F(1))),
        calibration_free_examples={str(I):calibrated_free_pulse_choice(assayref.band('pulses'),I) for I in [(F(3,10),F(49,100)),(F(81,100),F(82,100)),(F(59,100),F(61,100))]},
        correlated_error_examples=dict(positive=covariance_sign((F(2,100),F(3,100)),F(1,100)),unresolved=covariance_sign((F(2,100),F(3,100)),F(3,100))),
        paired_effort=dict(deadline_rate=1,expected_roots=2*ACCEPTED_DIVISIONS,expected_founder_time=ACCEPTED_DIVISIONS,no_timeout_roots=ACCEPTED_DIVISIONS,calibration_and_read_cost='additional'),
        scope='Calibration coverage must be charged jointly; the examples are supplied parameter sets, not estimated calibration. Correlated-error signs require a bound on mean conditional error covariance and the stated conditional-mean model. Expected collection cost is not a fixed-budget guarantee.'))
    # Matched sufficient direct-endpoint benchmark, not an optimality claim.
    g=slope/10-assayref.epsilon*assayref.W['A'][1];n=ceil(6*(F(12,1000)+2*g/3)/g**2)
    assert n==2688960
    fab=pulse_extinction(F(7,8),F(1,5),'AB');fba=pulse_extinction(F(7,8),F(1,5),'BA')
    dump('effect_and_effort.json',dict(direct_roots=2*n,per_arm=n,gap_lower=g,founders={str(n):dict(AB=fab**n,BA=fba**n,absolute_gap=fba**n-fab**n) for n in [1,2,10]},
        retained_only_worst_case_choice_error_lower=F(1,2),scope='Direct endpoint budget is a sufficient Bernstein bound for the same separated constant-action task at epsilon<=0.01. It does not establish an optimal sample-complexity advantage. One-founder survival stays near 97%; this is a ranking mechanism, not effective eradication.'))
    curve=[]
    for c in np.linspace(-.25,.25,41):
        source=FamilySource(SisterKernel(F(1,2),F(str(c))),DESCENDANT_DIVISION)
        d=source.extinction(phases['AB'])[0]-source.extinction(phases['BA'])[0];lo,hi=band.low_minus_high_extinction_interval(F(str(c)))
        curve.append((c,float(ref.reference_gap(F(str(c)))),d,float(lo),float(hi)))
    table('decision_curve.csv',['covariance','reference_gap','numerical_recurring_gap','certified_lower','certified_upper'],curve)
    scaling=[]
    for x in [F(7,8),F(49,50),F(99,100),F(999,1000),F(9999,10000)]:
        rr=ReferenceCertificate(x);cp,s=boundary(x);scaling.append((str(x),float(rr.h[1]),float(cp),float(rr.reversal_range(F(1,4))),float(rr.W['BA'][1]/s/rr.h[1])))
    table('short_pulse_scaling.csv',['x','pulse','threshold','certified_epsilon_supremum','band_halfwidth_per_epsilon_per_pulse'],scaling)
    plot(out,curve,rules,assayref)
    print(f'Configured single-P pulse choice: {band.decide(FOUNDER_COVARIANCE)}; paired-read decision: {current.decide(OBSERVED_AGREEMENTS)}.',flush=True)
    print('Exact recurring-division reversals, three full-class assay rules and the separate forty-division classification checked. Identical retained records omit the required sister dependence.',flush=True)
    here=Path(__file__).parent;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']},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,rules,reference):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    from scipy.stats import binom
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    hs=np.linspace(.005,.3,120);cp=[float(boundary(F(str(np.exp(-h))))[0]) for h in hs]
    axs[0].plot(hs,cp,label='Exact finite-duration boundary');axs[0].plot(hs,7*hs/8,ls='--',label='First-order approximation');axs[0].axhline(.25,color='gray',ls=':');axs[0].set(xlabel='Pulse length (model time)',ylabel='Sister covariance threshold',title='Sister-covariance decision threshold versus\npulse length');axs[0].legend(fontsize=8)
    data=np.array(curve);axs[1].fill_between(data[:,0],1e3*data[:,3],1e3*data[:,4],alpha=.2,label='Certified recurring-division band');axs[1].plot(data[:,0],1e3*data[:,1],label='No descendant division');axs[1].plot(data[:,0],1e3*data[:,2],ls='--',label='Recurring division: numerical');axs[1].axhline(0,color='black',lw=.7);axs[1].set(xlabel='Founder sister covariance',ylabel='AB minus BA extinction × 1000',title='Pulse-order extinction contrast by sister\ncovariance');axs[1].legend(fontsize=7)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'family_decision.png',dpi=180);fig.savefig(out/'family_decision.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    c=np.linspace(-.25,.25,301);p=.5+1.28*c
    for name,label in [('constant200','Constant: 200 pairs'),('pulse800','Pulses: 800 pairs'),('pulse1600','Pulses: 1600 pairs')]:
        r=rules[name];resolution=binom.cdf(r['low_cut'],r['n'],p)+binom.sf(r['high_cut']-1,r['n'],p);axs[0].plot(c,resolution,label=label)
    axs[0].axhline(.95,color='gray',ls=':');axs[0].set(xlabel='Sister covariance',ylabel='Probability of making a declaration',title='Probability that the paired-read rule\nresolves a decision',ylim=(0,1.02));axs[0].legend(fontsize=7)
    r=rules['pulse800'];axs[1].axvspan(0,r['low_cut'],color='green',alpha=.18,label='AB');axs[1].axvspan(r['low_cut'],r['high_cut'],color='gray',alpha=.18,label='Unresolved');axs[1].axvspan(r['high_cut'],800,color='orange',alpha=.18,label='BA');counts=np.arange(801)
    for cc in [.137064,.2]:axs[1].plot(counts,binom.pmf(counts,800,.5+1.28*cc),label=f'c={cc:.3f}')
    axs[1].set(xlabel='Agreements in 800 independent pairs',ylabel='Probability mass',xlim=(400,680),title='Paired-marker agreement counts and decision\nthresholds');axs[1].legend(fontsize=7)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'paired_assay.png',dpi=180);fig.savefig(out/'paired_assay.svg');plt.close(fig)


if __name__=='__main__':main()