Example code

"""A shared delivered course, inherited target and replenishing reserve.

Run: python example.py --output outputs
All rates and times are synthetic paper units, not calibrated clinical inputs.
"""
from fractions import Fraction as F

# EDITABLE SCENARIO. Fractions preserve decimal inputs exactly.
CAPACITY=400
SAFE_THRESHOLD=200                  # failure means strictly below this count
INITIAL_RESERVE=278
RETURN_ANCHOR=278
RETAINED_RENEWAL=F(1)               # post-intervention effective renewal
MORTALITY_CEILING=F(61,200)
SELECTIVITY_MIN=F(2)
ADDITIONAL_HEALTHY_DEATH_MAX=F(1,200)
ADMINISTRATION_RATE=F(29,100)
ADMINISTRATION_STOP=F(112)
DEADLINE=F(120)
CLEARANCE_RANGE=(F(99,100),F(101,100))
TARGET_INITIAL_COUNTS=(0,0,0,0,0,4)  # source order 00,0R,RR,A0,AR,AA
TARGET_RELATIVE_RATE_ERROR=F(1,100)
TARGET_TOLERANCE=F(1,100)
RESERVE_TOLERANCE=F(1,100)
MANUSCRIPT_SHA256='0328a5cb97c00c10c60e89758afed595db34421ce43e2e59af4190ad82e052ee'

import argparse,csv,hashlib,json,platform
from pathlib import Path
from math import ceil,floor,exp,log
import numpy as np
from reserve import Reserve,ReserveChain,latency_floor,fluid_limit_diagnostic
from mission import AdministrationCourse,HealthyResponse,TargetDriftCertificate,assemble_mission,numerical_target


def serialize(value):
    if isinstance(value,F):return str(value)
    if isinstance(value,np.ndarray):return value.tolist()
    if isinstance(value,np.generic):return value.item()
    raise TypeError(type(value).__name__)


def capacity_scan(maximum=300):
    rows=[]
    for K in range(1,maximum+1):
        h=ceil(F(K,2));M=floor(F(139*K,200))
        if M<h:continue
        result=Reserve(K,h,F(1),F(61,200)).anchored(M,120).evaluate(M)
        rows.append((K,h,M,result['sharp_bound'],result['simple_product_bound']))
    return rows


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 dump(name,value):(out/name).write_text(json.dumps(value,default=serialize,indent=2)+'\n',encoding='utf-8')
    def table(name,headers,rows):
        with (out/name).open('w',newline='',encoding='utf-8') as f:
            writer=csv.writer(f);writer.writerow(headers);writer.writerows(rows)
    course=AdministrationCourse(ADMINISTRATION_RATE,ADMINISTRATION_STOP,DEADLINE,*CLEARANCE_RANGE)
    healthy=HealthyResponse(SELECTIVITY_MIN,ADDITIONAL_HEALTHY_DEATH_MAX)
    reserve=Reserve(CAPACITY,SAFE_THRESHOLD,RETAINED_RENEWAL,MORTALITY_CEILING)
    target=TargetDriftCertificate(relative_error=TARGET_RELATIVE_RATE_ERROR).certify_course(course,TARGET_INITIAL_COUNTS)
    band=target['delivery'];mortality=healthy.ceiling(band['concentration_upper'])
    if mortality>MORTALITY_CEILING:raise ValueError('Shared delivered input exceeds the declared reserve mortality ceiling.')
    constraints=dict(administration=course.rate<=F(3,10),concentration=band['concentration_upper']<=F(3,10),amount=band['administered_amount']<=33,exposure=band['full_horizon_exposure_upper']<=33)
    certificate=reserve.anchored(RETURN_ANCHOR,DEADLINE);safety=certificate.evaluate(INITIAL_RESERVE)
    mission=assemble_mission(target['exact_rational_upper'],safety['sharp_bound'],TARGET_TOLERANCE,RESERVE_TOLERANCE)
    mission['resource_constraints']=constraints;mission['all_declared_requirements_pass']=all(constraints.values()) and mission['target_passes'] and mission['reserve_passes']
    dump('configured_mission.json',dict(target=target,reserve=safety,mortality_bound=mortality,mission=mission,
        comparison_assumptions='Unit births/deaths; actual birth at count h is at least retained_renewal*h*(1-h/K), death at most m*h, no immigration into target, no delayed births or catastrophe jumps. These are model premises, not fitted conclusions.'))
    # Fixed reference casebook, explicitly independent of edited scenario.
    ref=Reserve(400,200,F(1),F(61,200));cert=ref.anchored(278,120)
    preparation={str(h):cert.evaluate(h) for h in [199,200,208,209,221,240,278,400]}
    mixture={199:F(1,100),279:F(99,100)}
    dump('preparation.json',dict(reference=preparation,random_preparation=mixture,random_bound=cert.random_preparation(mixture),
        published_joint_bound=assemble_mission(F(107,12500),cert.evaluate(278)['simple_product_bound']),
        scope='An initial count below threshold is already a failure. A high mean preparation does not eliminate its low-count tail. Published target coarsening is 0.00856; the fresh exponential series gives a stronger bound.'))
    table('return_products.csv',['count_j','exact_product','decimal_product'],[(j,p,float(p)) for j,p in cert.products.items()])
    rows=capacity_scan();minimum=lambda col:next(K for K,h,M,sharp,simple in rows if (sharp,simple)[col]<=F(1,100))
    assert (minimum(0),minimum(1))==(232,252)
    table('capacity_scan.csv',['capacity','threshold','anchor','sharp_bound','simple_bound'],[(K,h,M,float(s),float(p)) for K,h,M,s,p in rows])
    dump('capacity_law.json',dict(minimum_sharp=minimum(0),minimum_simple=minimum(1),scan_maximum=300,reference=ref.capacity_envelope(F(1,2),120),
        scope='Minimum within this specified rounded design family and sufficient certificate, not a necessary biological capacity. The action law requires 1-m/r > threshold fraction.'))
    # Actual time-varying comparison uses the same delivered concentration.
    death=lambda t:healthy.mortality(float(course.concentration(t,float(course.clearance_min))))
    chain=ReserveChain(reserve,death);times,loss,diagnostic=chain.solve(INITIAL_RESERVE,DEADLINE)
    _,terminal,terminal_diag=chain.solve(INITIAL_RESERVE,DEADLINE,killed=False)
    _,constant_loss,constant_diag=ReserveChain(reserve).solve(INITIAL_RESERVE,DEADLINE)
    concentration=course.concentration(times,float(course.clearance_min))
    table('shared_course.csv',['time','administered_rate','concentration_slow_clearance','healthy_mortality','any_time_loss','terminal_shortfall','constant_ceiling_loss'],
        zip(times,np.where(times<float(course.stop),float(course.rate),0),concentration,healthy.mortality(concentration),loss,terminal,constant_loss))
    dump('numerical_diagnostics.json',dict(target=numerical_target(course,initial_counts=TARGET_INITIAL_COUNTS),reserve=diagnostic,terminal_chain=terminal_diag,constant_chain=constant_diag,
        any_time_loss=float(loss[-1]),terminal_shortfall=float(terminal[-1]),constant_ceiling_loss=float(constant_loss[-1]),scope='Floating-point solutions are diagnostics, not the exact rational path-risk certificate.'))
    renewal=[]
    for K,h in [(4,3),(4,2),(5,2)]:
        for r in [4,8,16,32,64,128]:
            R=Reserve(K,h,F(r),F(1,2));p=ReserveChain(R).high_precision_constant_risk(K,F(2));d=K-h;coefficient=R.high_renewal_coefficient(2)
            renewal.append((K,h,d,r,p,float(p)*r**d,float(coefficient),float(coefficient/F(r)**d)))
    table('high_renewal.csv',['capacity','threshold','buffer','renewal','risk_60_digit_arithmetic','scaled_risk','limiting_coefficient','sufficient_bound'],renewal)
    latency=latency_floor(20,10,F(1,2));assert latency>F(41,100)
    dump('boundary_cases.json',dict(latency_floor=latency,latency_description='No replacement before latency, empty initial pipeline, 20 initially alive units each surviving with probability at most 1/2; first loss cannot be undone.',
        zero_buffer=dict(K=20,threshold=20,mortality=F(1,2),horizon=F(2),risk_formula='1-exp(-K*m*T)',numerical_risk=1-exp(-20)),
        low_equilibrium=fluid_limit_diagnostic(1,F(3,5),F(1,2),8,400),
        baseline_unsafe_upper_terminal_safety=2*exp(-17.88),
        scope='The baseline-unsafe example and constant-source fluid obstruction are not a treatment selectivity frontier. Arbitrarily large renewal cannot repair zero buffer or an already missed latency requirement.'))
    plot(out,rows,preparation,renewal,times,concentration,loss,terminal,constant_loss,safety)
    print(f'Configured mission passes: {mission["all_declared_requirements_pass"]}; target risk <= {float(mission["target_upper"]):.10g}, reserve path loss <= {float(mission["reserve_loss_upper"]):.10g}.',flush=True)
    print(f'Joint success >= {float(mission["joint_success_lower"]):.10g}; no independence assumed. Rounded capacity certificate first passes at 232 (sharp) and 252 (simple).',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,rows,preparation,renewal,t,c,loss,terminal,constant,safety):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    axs[0].plot(t,c,label='Delivered concentration');axs[0].axvspan(4,112,color='green',alpha=.1,label='Certified contraction band');axs[0].axvline(112,color='gray',ls=':');axs[0].set(xlabel='Paper time units',ylabel='Delivered input',title='Delivered input during treatment and\nwashout');axs[0].legend(fontsize=8)
    for y,label in [(loss,'Any-time loss, shared course'),(terminal,'Terminal shortfall, shared course'),(constant,'Any-time loss, constant ceiling')]:axs[1].semilogy(t[1:],np.maximum(y[1:],1e-30),label=label)
    axs[1].axhline(float(safety['sharp_bound']),color='black',ls='--',label='Exact comparison certificate');axs[1].set(xlabel='Paper time units',ylabel='Reserve failure probability',ylim=(1e-14,1e-4),title='Any-time reserve failure and terminal\nshortfall');axs[1].legend(fontsize=7)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'mission.png',dpi=180);fig.savefig(out/'mission.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    axs[0].semilogy([r[0] for r in rows],[float(r[3]) for r in rows],label='Return-normalized certificate');axs[0].semilogy([r[0] for r in rows],[float(r[4]) for r in rows],label='Simple product certificate');axs[0].axhline(.01,color='black',ls=':');axs[0].set(xlabel='Capacity K (threshold ceil(K/2))',ylabel='All-time risk upper bound',title='Reserve-failure upper bounds versus\ncapacity');axs[0].legend(fontsize=8)
    for K,h in [(4,3),(4,2),(5,2)]:
        rr=[r for r in renewal if r[0]==K and r[1]==h];axs[1].semilogx([r[3] for r in rr],[r[5]/r[6] for r in rr],marker='o',label=f'K={K}, threshold={h}, buffer={K-h}')
    axs[1].axhline(1,color='black',ls=':');axs[1].set(xlabel='Retained renewal rate',ylabel='Risk × renewal^buffer / coefficient',title='Scaled reserve risk versus renewal rate');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'capacity.png',dpi=180);fig.savefig(out/'capacity.svg');plt.close(fig)


if __name__=='__main__':main()