Example code

"""Finite-fuel chemical correction, productive readout and two-daughter reuse."""
from fractions import Fraction as F

# EDITABLE literal-source inputs. These are abstract counts/rates, not calibration.
REGIME_NAME='WO'
CORE_INVENTORY=80
FUEL_INVENTORY=1
MINIMUM_RESIDENTS=8
MINORITY_ALLOWANCE=1
INITIAL_MAJOR=8
INITIAL_MINOR=1
CORRECTION_RATE=F(2*10**7)
REVERSE_CORRECTION_RATE=F(6,10**11)
LEAKAGE_RATE=F(1,10**7)
PREFIX_DURATION=F(1,8*10**7)
CORE_DURATION=F(1)
RECOVERY=F(9,10)
PARTITION_PROBABILITY=F(12,25)
LINEAGE_ROUNDS=5
SEED=6901
EVENT_BUDGET=200000
MANUSCRIPT_SHA256='a6f9278c4667dd096f85e660e71fdd87a367643df1f8f45dd354749095c16f42'

import argparse,csv,hashlib,json,platform
from pathlib import Path
from math import log,log2
import numpy as np
from chemistry import Rates,ChemicalSource,BatchOperations,ProductRetention,SPECIES
from core import PureCore,binomial_event,cluster_partition
from correction import CorrectionWalk,REGIMES,one_minority_bound,fuel_deadline,finite_fuel_bound,selection_bound


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)
    # Recompute reference certificates; no cached manuscript numerators are inputs.
    core=PureCore();ideal=core.certify();operations=core.certify(F(12,25),F(9,10))
    assert ideal['lower']==F(2147232289,2**31) and operations['lower']==F(2147153442,2**31)
    table('core_lower_vector.csv',['residents','terminal_product','ideal_lower_numerator','imperfect_lower_numerator','denominator'],[(n,p,int(ideal['values'][i]),int(operations['values'][i]),2**31) for i,(n,p) in enumerate(core.states)])
    table('newborn_bounds.csv',['initial_residents','ideal_lower','imperfect_lower'],[(n,float(F(int(ideal['values'][core.index[(n,0)]]),2**31)),float(F(int(operations['values'][core.index[(n,0)]]),2**31))) for n in range(8,81)])
    dump('integer_core_certificates.json',dict(ideal={k:v for k,v in ideal.items() if k!='values'},imperfect={k:v for k,v in operations.items() if k!='values'},
        endpoint_reduction='Recovery is monotone; the complementary split payoff is symmetric about 1/2 and increases up to 1/2. Thus theta in [0.48,0.52], eta in [0.9,1] is covered by (0.48,0.9), not by a parameter grid.'))
    regime_results={r.name:one_minority_bound(r,operations['lower'] if r.imperfect else ideal['lower']) for r in REGIMES if r.minority==1}
    assert all(row['passes'] for row in regime_results.values())
    fuelreg=REGIMES[-1];deadline=fuel_deadline(fuelreg);fuelrows=[]
    for j in range(3):
        for r in range(8+j,81):
            law=CorrectionWalk(r).absorption(j,64);value=finite_fuel_bound(operations['lower'],law,0,fuelreg)
            fuelrows.append((r,j,law['moments'][0]['mass'],law['moments'][0]['spent'],value))
    worst=min(fuelrows,key=lambda row:row[-1]);assert worst[:2]==(10,2) and worst[-1]>F(2479,2500) and len(fuelrows)==216
    table('all_restart_fuel_bounds.csv',['resident_total','minority','correct_absorption_mass','unnormalized_spent_fuel','joint_lower'],fuelrows)
    law=CorrectionWalk(10).absorption(2,64);switch=finite_fuel_bound(operations['lower'],law,10,fuelreg);pure=80*fuelreg.epsilon_max*(1+fuelreg.prefix_duration)
    assert switch>F(37,5000) and switch/pure>935
    coarse=operations['lower']-(1-law['moments'][0]['mass'])-80*fuelreg.epsilon_max-64*80**3*fuelreg.beta_max-3216*fuelreg.prefix_duration-F(1,10**12)
    regime_results['F']=dict(lower=worst[-1],worst_state=worst[:2],all_states=216,headline=fuelreg.headline,coarse_full_inventory_charge=coarse,spent_fuel_gain=worst[-1]-coarse,deadline=deadline)
    dump('regime_certificates.json',regime_results)
    table('absorption_spent_fuel.csv',['jump','endpoint','exact_probability'],law['hits'])
    dump('variation.json',dict(moments=law['moments'],transient_mass=sum(law['transient'].values()),correct_consensus_ceiling=F(127,128),source_bound=worst[-1],gap_to_correction_ceiling=F(127,128)-worst[-1],
        mixed_opposite_pair_lower=switch,pure_opposite_pair_upper=pure,ratio_lower=switch/pure,at_least_one_switch_128_independent_lower=1-(1-F(37,5000))**128,
        scope='The mixed preparation is (8,2,70,0,0,64,0). Pure newborns have no waste: only leakage can create the first minority. The ceiling concerns the correction-only route, not every possible full-source trajectory. A label-only mutation rate cannot describe both preparations.'))
    fuelcurves=[]
    for r in [10,12,16]:
        walk=CorrectionWalk(r)
        for L in range(65):
            a=walk.absorption(2,L);fuelcurves.append((r,L,float(a['moments'][0]['mass']),float(1-walk.wrong_consensus(2))))
    table('fuel_sweep.csv',['residents','available_fuel','correct_absorption','unlimited_correct_limit'],fuelcurves)
    # Literal source uses the editable rates; the reference guarantee is gated.
    rates=Rates(CORRECTION_RATE,REVERSE_CORRECTION_RATE,LEAKAGE_RATE);source=ChemicalSource(rates)
    ops=BatchOperations(CORE_INVENTORY,FUEL_INVENTORY,MINIMUM_RESIDENTS,MINORITY_ALLOWANCE,RECOVERY,PARTITION_PROBABILITY)
    if not isinstance(INITIAL_MAJOR,int) or not isinstance(INITIAL_MINOR,int) or INITIAL_MAJOR<0 or INITIAL_MINOR<0 or INITIAL_MAJOR+INITIAL_MINOR>CORE_INVENTORY:raise ValueError('Invalid initial resident inventory.')
    initial=np.array([INITIAL_MAJOR,INITIAL_MINOR,CORE_INVENTORY-INITIAL_MAJOR-INITIAL_MINOR,0,0,FUEL_INVENTORY,0],dtype=np.int64)
    regime=next(r for r in REGIMES if r.name==REGIME_NAME)
    conditions=dict(core=CORE_INVENTORY==80,minimum=MINIMUM_RESIDENTS==8,fuel=FUEL_INVENTORY==regime.fuel,minority=MINORITY_ALLOWANCE==regime.minority,
        gamma=regime.gamma_min<=CORRECTION_RATE<=2*10**9,beta=F(1,10**13)<=REVERSE_CORRECTION_RATE<=regime.beta_max,epsilon=F(1,10**10)<=LEAKAGE_RATE<=regime.epsilon_max,
        deadline=PREFIX_DURATION==regime.prefix_duration and CORE_DURATION==1,
        operations=(F(12,25)<=PARTITION_PROBABILITY<=F(13,25) and RECOVERY>=F(9,10)) if regime.imperfect else (PARTITION_PROBABILITY==F(1,2) and RECOVERY==1),
        initial_restart=ops.region(initial,'X'))
    rng=np.random.default_rng(SEED);horizon=float(PREFIX_DURATION+CORE_DURATION)
    simulation=source.simulate(initial,horizon,rng,EVENT_BUDGET)
    if simulation['status']!='completed':raise RuntimeError('Simulation event budget exhausted; no completed batch or success score is issued.')
    operation=ops.apply(simulation['state'],rng)
    table('literal_trajectory.csv',['time',*SPECIES,'channel_index'],simulation['trace'])
    dump('literal_source.json',dict(rates=vars(rates),positive_equilibrium_activities=source.equilibrium_activities(),reaction_table=[dict(name=r.name,reactants=r.reactants,products=r.products,coefficient=r.rate,change=r.change) for r in source.reactions],
        reference_premises=conditions,reference_guarantee_applicable=all(conditions.values()),reference_lower=regime.headline if all(conditions.values()) else None,
        initial=initial,terminal=simulation['state'],reaction_counts={r.name:int(n) for r,n in zip(source.reactions,simulation['counts'])},operations=operation,joint_event=ops.score(operation,'X'),
        scope='A seeded literal Gillespie trajectory, not evidence of the high-probability theorem. Harvest uses terminal stock, not forward reaction counts. All rates run throughout the batch.'))
    # Follow daughter A chosen in advance; score both daughters each round.
    lineage=[];state=initial.copy();spent_food=spent_fuel=0;collected=np.zeros(2,dtype=int)
    for round_ in range(1,LINEAGE_ROUNDS+1):
        sim=source.simulate(state,horizon,rng,EVENT_BUDGET,record=False)
        if sim['status']!='completed':raise RuntimeError('Lineage simulation incomplete at event budget.')
        result=ops.apply(sim['state'],rng);collected+=result['credited'];spent_food+=result['food_supplied'];spent_fuel+=result['fuel_supplied']
        lineage.append(dict(round=round_,initial=state.copy(),terminal=sim['state'],operations=result,joint_event=ops.score(result,'X')));state=result['daughters'][0].copy()
    dump('lineage.json',dict(rounds=lineage,collected=collected,food_supplied=spent_food,fuel_supplied=spent_fuel,
        conditional_uniform_lower=regime.headline**LINEAGE_ROUNDS if all(conditions.values()) else None,
        scope='Preselected daughter A continues; both daughters are scored and refilled. Failed/out-of-region states are not reset or discarded. Conditional uniform bounds do not assume independent cycles and do not certify an exponential tree by q^depth.'))
    # External selection reads only product. Every daughter is refilled first.
    retention=ProductRetention();retained=[];population_supply=0;batchscores=[]
    for label in ['X','Y']:
        for k in range(10):
            state=initial.copy()
            if label=='Y':state[[0,1]]=state[[1,0]]
            sim=source.simulate(state,horizon,rng,EVENT_BUDGET,record=False)
            if sim['status']!='completed':raise RuntimeError('Selection-round simulation incomplete.')
            result=ops.apply(sim['state'],rng);population_supply+=result['food_supplied']+result['fuel_supplied'];batchscores.append(ops.score(result,label));retained.extend(retention.retain(result,rng))
    counts={label:sum(ops.region(d,label) for d in retained) for label in ['X','Y']};counts['outside']=len(retained)-counts['X']-counts['Y']
    dump('selection_and_information.json',dict(actual_retained=counts,actual_all_chemical_events=all(batchscores),actual_refill_supply=population_supply,
        fixed_reference={r.name:dict(one_round_lower=selection_bound(r.headline),ten_lineage_divisions_lower=r.headline**10,information_bits_display=1-binary_entropy(float(1-r.headline))) for r in REGIMES},
        scope='Selection certificate uses a union bound over 20 chemical failures and independent fair retention coins for Y daughters, not independence of chemical failures. Fano displays are numerical readings of the probability bound, not measured information rates.'))
    clusters={str(w):cluster_partition([w]*(32//w)) for w in [1,2,4]}
    NA=6.02214076e23;volume=1e-15;t0=NA*volume/1e6
    dump('mechanism_and_inventory.json',dict(cluster_success=clusters,elementary_first_association_upper=F(56,4*10**8),
        reference_material=dict(WO_initial=81,F_initial=144,WO_ten_divisions_supply=1701,F_ten_divisions_supply=3024,food_per_division_upper=160),
        dimensional_display=dict(assumed_volume_litres=volume,assumed_second_order_rate=1e6,time_unit_seconds=t0,WO_prefix_microseconds=t0/(8e7)*1e6,fourth_order_rate_gamma_1e8=1e8*(NA*volume)**3/t0),
        scope='Intact complexes change the partition polynomial; matching an effective correction rate cannot transfer the monomer certificate. Physical scales are illustrative and the kinetic demand is extreme. Solvent, containment, apparatus and volume-reset work are outside the modeled material bill.'))
    plot(out,core,ideal,operations,fuelcurves,regime_results,switch,pure,clusters,simulation)
    print(f'Integer core minima: {ideal["lower"]} (ideal), {operations["lower"]} (imperfect); 3240 states freshly replayed.',flush=True)
    print(f'WO joint bound {float(regime_results["WO"]["lower"]):.12g}; finite-fuel F bound {float(worst[-1]):.12g}; mixed opposite-pair lower {float(switch):.12g}.',flush=True)
    print(f'Configured certificate applicable: {all(conditions.values())}; simulated joint event: {ops.score(operation,"X")}. Simulation is illustrative.',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 binary_entropy(p):return 0. if p in [0,1] else -p*log2(p)-(1-p)*log2(1-p)


def plot(out,core,ideal,operations,fuel,regimes,switch,pure,clusters,simulation):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained');ns=list(range(8,81))
    for cert,label in [(ideal,'Ideal operations'),(operations,'90% recovery; biased division')]:axs[0].semilogy(ns,[1-int(cert['values'][core.index[(n,0)]])/cert['scale'] for n in ns],label=label)
    axs[0].set(xlabel='Initial pure residents (food = 80 − residents)',ylabel='Certified joint failure allowance',title='Joint output and two-daughter restart\nfailure bound');axs[0].legend(fontsize=8)
    for r in [10,12,16]:
        rows=[row for row in fuel if row[0]==r];axs[1].plot([row[1] for row in rows],[row[2] for row in rows],label=f'{r} residents');axs[1].axhline(rows[0][3],color='gray',ls=':',lw=.6)
    axs[1].set(xlabel='Fuel units available',ylabel='Correct-consensus probability',ylim=(.84,1.005),title='Correct-consensus probability versus\navailable fuel');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'contract.png',dpi=180);fig.savefig(out/'contract.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    axs[0].barh(['Pure newborn: upper','Mixed newborn: lower'],[float(pure),float(switch)]);axs[0].set(xscale='log',xlabel='Probability of productive opposite daughters',title='Opposite-label daughter bounds by newborn\ncomposition')
    axs[1].bar([str(w) for w in [1,2,4]],[1-float(clusters[str(w)]) for w in [1,2,4]]);axs[1].set(yscale='log',xlabel='Moieties per intact cluster (32 total)',ylabel='Failure: a daughter receives fewer than 8',title='Partition failure versus molecular cluster\nsize')
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'variation.png',dpi=180);fig.savefig(out/'variation.svg');plt.close(fig)


if __name__=='__main__':main()