"""Driven-cleavage reactor: literal counts, repeated-window ledger, and bounds."""
from __future__ import annotations

# EDITABLE INPUTS -------------------------------------------------------------
COPY_SCALE = 100_000_000
KINETIC_BIAS = '10'
DUPLEX_RELEASE = '20'
DRIVEN_CLEAVAGE = '1/20'
POST_STARTUP_DURATION = '100'
TARGET_FAILURE = '1/100'
BASAL_RATE = '1/500000000'  # fixed by the theorem; changing it exits its scope
REFERENCE_MOLAR = '1/1000000'
DILUTION_PER_HOUR = '1'
TEMPERATURE_K = '298'
SSA_EVENT_BUDGET = 20_000
SSA_SEED = 25092026
MANUSCRIPT_SHA256 = 'b03c9193c19090688e71a36a46ad31bf33dc9bde1d910c5b0cf3615ae803b05e'
# ---------------------------------------------------------------------------

from dataclasses import dataclass
from fractions import Fraction as Q
import argparse
import csv
import hashlib
import json
import math
from pathlib import Path
import platform

import numpy as np
from scipy.integrate import solve_ivp
from mpmath import mp, iv

mp.dps=70;iv.dps=70
SPECIES=('U','W','X','C1','C2','Z')
A_WEIGHTS=(1,0,1,2,2,2);B_WEIGHTS=(0,1,1,1,2,2)
Y_WEIGHTS=(Q(0),Q(0),Q(1),Q(9,8),Q(7,5),Q(9,5))
M_WEIGHTS=(2,2,4,6,8,8);S_WEIGHTS=(0,0,4,4,4,8)
LEDGER=('export','food_U','food_W','drive_F','drive_P','wash_A','wash_B','wash_mass')


def dot(a,b):return sum(x*y for x,y in zip(a,b))
def number(ctx,x):
    x=Q(str(x));return ctx.mpf(x.numerator)/x.denominator
def logsum(ctx,values):
    values=list(values);shift=values[0]
    return shift+ctx.log(sum(ctx.exp(v-shift) for v in values))


@dataclass(frozen=True)
class Parameters:
    K: object=KINETIC_BIAS
    release: object=DUPLEX_RELEASE
    delta: object=DRIVEN_CLEAVAGE
    epsilon: object=BASAL_RATE
    def __post_init__(self):
        for name in ('K','release','delta','epsilon'):
            q=Q(str(getattr(self,name)))
            if q<=0:raise ValueError('Positive reaction coefficients required.')
            object.__setattr__(self,name,q)
    @property
    def eta(self):return self.epsilon/16
    def theorem_scope(self):
        return 8<=self.K<=12 and 18<=self.release<=22 and Q(1,100)<=self.delta<=Q(1,20) and self.epsilon==Q(1,500000000)


@dataclass(frozen=True)
class Channel:
    label: int
    name: str
    inputs: tuple[int,...]
    outputs: tuple[int,...]
    coefficient: Q
    marks: tuple[int,...]=(0,)*8
    @property
    def jump(self):return tuple(self.outputs.count(i)-self.inputs.count(i) for i in range(6))
    def count_rate(self,counts,V):
        rate=self.coefficient*Q(V)**(1-len(self.inputs));seen={}
        for i in self.inputs:
            rate*=counts[i]-seen.get(i,0);seen[i]=seen.get(i,0)+1
        return rate


class DrivenReactor:
    """Preserves original labels after deleting only channels 6 and 7."""
    def __init__(self,parameters=Parameters(),enabled=True):
        self.p=parameters;self.enabled=enabled;p=parameters
        pairs=(('basal',(0,1),(2,),p.epsilon,p.epsilon/p.K),
               ('bind_U',(2,0),(3,),Q(20),Q(20)),
               ('bind_W',(3,1),(4,),Q(20),Q(20)),
               ('ligation',(4,),(5,),Q(20),Q(20)/p.K),
               ('release',(5,),(2,2),p.release,p.release))
        channels=[]
        for j,(name,left,right,kp,km) in enumerate(pairs):
            channels.extend([Channel(2*j,name+'_forward',left,right,kp),Channel(2*j+1,name+'_reverse',right,left,km)])
        for i in (0,1):
            marks=[0]*8;marks[1+i]=1
            channels.append(Channel(10+i,'feed_'+SPECIES[i],(),(i,),Q(1),tuple(marks)))
        for i in range(6):
            marks=(S_WEIGHTS[i],0,0,0,0,A_WEIGHTS[i],B_WEIGHTS[i],M_WEIGHTS[i])
            channels.append(Channel(12+i,'wash_'+SPECIES[i],(i,),(),Q(1),marks))
        channels.extend([Channel(18,'driven_cleavage',(2,),(0,1),p.delta,(0,0,0,1,0,0,0,0)),
                         Channel(19,'driven_reverse',(0,1),(2,),p.delta*p.eta,(0,0,0,0,1,0,0,0))])
        self.channels=tuple(c for c in channels if enabled or c.label not in (6,7))
        self.by_label={c.label:c for c in self.channels}
        self.jumps=np.array([c.jump for c in self.channels],int);self.marks=np.array([c.marks for c in self.channels],int)
        self.coefficients=np.array([float(c.coefficient) for c in self.channels])
        self.orders=np.array([len(c.inputs) for c in self.channels])
        self.inputs=np.full((len(self.channels),2),6,int);self.offsets=np.zeros((len(self.channels),2),int)
        for j,c in enumerate(self.channels):
            seen={}
            for k,i in enumerate(c.inputs):
                self.inputs[j,k]=i;self.offsets[j,k]=seen.get(i,0);seen[i]=seen.get(i,0)+1

    def rates(self,state,V=None):
        factors=np.r_[state,1][self.inputs]
        if V is not None:factors=np.maximum(0,factors-self.offsets)
        rate=self.coefficients*np.prod(factors,axis=1)
        return rate if V is None else rate*np.power(float(V),1-self.orders)

    def generator(self,counts,V,weights):
        rates=[c.count_rate(counts,V) for c in self.channels]
        jumps=[dot(weights,c.jump) for c in self.channels]
        return sum(a*d for a,d in zip(rates,jumps)),sum(a*d*d for a,d in zip(rates,jumps))

    def deterministic(self,H,method='Radau'):
        H=float(Q(str(H)))
        if not 1<=H<=10000:raise ValueError('ODE illustration is limited to duration 1..10000; bounds have a separate interface.')
        def rhs(t,state):
            rates=self.rates(state[:6]);return np.r_[rates@self.jumps,rates@self.marks]
        def entry(t,state):return np.array(Y_WEIGHTS,float)@state[:6]-1/2500
        entry.direction=1
        sol=solve_ivp(rhs,(0,500+H),[1,1,0,0,0,0]+[0]*8,method=method,dense_output=True,events=entry,
                       rtol=1e-10,atol=1e-16,max_step=1)
        if not sol.success:raise ArithmeticError(sol.message)
        return sol

    def stochastic(self,policy,limit=SSA_EVENT_BUDGET,seed=SSA_SEED):
        if not 1<=policy.V<10**12 or limit<1:raise ValueError('Direct SSA needs V<10^12 and a positive event budget.')
        rng=np.random.default_rng(seed);monitor=OperatingMonitor(self,policy);events=0
        while monitor.status=='active' and events<limit:
            rates=self.rates(monitor.counts,policy.V);total=float(sum(rates))
            event_time=monitor.time+rng.exponential(1/total)
            boundary=float(monitor.next_boundary())
            if event_time>boundary:
                monitor.advance(boundary);continue
            j=int(np.searchsorted(np.cumsum(rates),rng.random()*total,side='right'))
            monitor.event(event_time,self.channels[j].label);events+=1
        report=monitor.report();report.update(events=events)
        if monitor.status=='active':report.update(status='unfinished',event_verdict=None)
        return report


@dataclass(frozen=True)
class ObservationPolicy:
    V: int=COPY_SCALE
    H: object=POST_STARTUP_DURATION
    def __post_init__(self):
        object.__setattr__(self,'H',Q(str(self.H)))
        if not isinstance(self.V,int) or self.V<1 or self.H<1:raise ValueError('Positive integer copy scale and H>=1 required.')
    @property
    def end(self):return 500+self.H
    @property
    def windows(self):return self.H.numerator//self.H.denominator
    @property
    def export_cap(self):return math.ceil(Q(self.V,5000))
    @property
    def supply_cap(self):return math.ceil(2*self.V*self.end)+1


class OperatingMonitor:
    """Tracks one continuing path. Boundaries reset only the window counter.

    Events exactly at a window's right endpoint belong to that window. Entry
    at time 500 is allowed. Return at Y<=V/5000 fails (including equality).
    Natural ledgers provide a second audit of the capped measurement counters.
    """
    def __init__(self,reactor,policy):
        if policy.windows>100000:raise ValueError('Monitor window limit; use ProbabilityBudget for huge horizons.')
        self.reactor=reactor;self.policy=policy;self.counts=[policy.V,policy.V,0,0,0,0]
        self.ledger=[0]*8;self.capped_supplies=[0]*4;self.counter=0
        self.time=0.;self.phase='pre_entry';self.entry_time=None;self.status='active'
        self.started=False;self.completed=0;self.all_windows_pass=True;self.window_exports=[]

    def next_boundary(self):
        if not self.started:return Q(500)
        return Q(501+self.completed) if self.completed<self.policy.windows else self.policy.end

    def corridor(self):
        V=self.policy.V
        return all(9*V<=10*dot(w,self.counts)<=11*V for w in (A_WEIGHTS,B_WEIGHTS))

    def _boundary(self):
        if not self.started:
            self.started=True
            if self.reactor.enabled and self.phase=='pre_entry':self.status='missed_deadline'
        elif self.completed<self.policy.windows:
            self.window_exports.append(self.counter)
            self.all_windows_pass &= self.counter>=self.policy.export_cap
            self.completed+=1;self.counter=0
        if self.status=='active' and self.time>=float(self.policy.end):self.status='complete'

    def advance(self,t):
        if t<self.time or t>float(self.policy.end):raise ValueError('Invalid monitor time.')
        while self.status=='active' and float(self.next_boundary())<=t:
            self.time=float(self.next_boundary());self._boundary()
        if self.status=='active':self.time=t

    def event(self,t,label):
        if self.status!='active' or t<self.time or t>float(self.policy.end):raise ValueError('Event outside active observation interval.')
        while self.status=='active' and float(self.next_boundary())<t:
            self.time=float(self.next_boundary());self._boundary()
        if self.status!='active':return
        channel=self.reactor.by_label[label]
        if channel.count_rate(self.counts,self.policy.V)<=0:raise ValueError('Reaction is not supported by molecule counts.')
        self.time=t;self.counts=[n+d for n,d in zip(self.counts,channel.jump)]
        if min(self.counts)<0:raise AssertionError('Count transition left the orthant.')
        self.ledger=[n+d for n,d in zip(self.ledger,channel.marks)]
        self.capped_supplies=[min(self.policy.supply_cap,n+d) for n,d in zip(self.capped_supplies,channel.marks[1:5])]
        if t>500 and self.completed<self.policy.windows:
            self.counter=min(self.policy.export_cap,self.counter+channel.marks[0])
        V=self.policy.V;Y40=dot((0,0,40,45,56,72),self.counts)
        if self.reactor.enabled:
            if self.phase=='pre_entry' and 2500*Y40>=40*V:self.phase='entered';self.entry_time=t
            elif self.phase=='entered' and 5000*Y40<=40*V:self.phase='returned';self.status='return_failure'
        if not self.corridor():self.status='resource_exit'
        self.audit_ledger()
        # Apply deterministic gates after a reaction on their exact endpoint.
        if self.status=='active' and float(self.next_boundary())==t:self._boundary()

    def audit_ledger(self):
        V=self.policy.V;E,IU,IW,QF,QP,OA,OB,MW=self.ledger
        if dot(A_WEIGHTS,self.counts)+OA!=V+IU or dot(B_WEIGHTS,self.counts)+OB!=V+IW:
            raise AssertionError('Moiety ledger failed.')
        if dot(M_WEIGHTS,self.counts)+MW!=4*V+2*(IU+IW):raise AssertionError('Mass ledger failed.')
        if OA*4<E or OB*4<E:raise AssertionError('Export exceeded washed moiety accounting.')
        if self.capped_supplies!=[min(self.policy.supply_cap,n) for n in self.ledger[1:5]]:
            raise AssertionError('Capped counters no longer match natural histories.')

    def report(self):
        V=self.policy.V;T=self.policy.end;IU,IW,QF,QP=self.capped_supplies
        supplies=IU<=2*V*T and IW<=2*V*T and QF+QP<=V*T/8
        verdict=None
        if self.status!='active':
            verdict=(self.status=='complete' and self.all_windows_pass and supplies) if self.reactor.enabled else (
                self.status=='resource_exit' or (self.status=='complete' and self.all_windows_pass))
        return {'status':self.status,'time':self.time,'phase':self.phase,'entry_time':self.entry_time,
            'counts':self.counts,'natural_ledger':dict(zip(LEDGER,self.ledger)),
            'capped_supply_counts':self.capped_supplies,'current_window_counter':self.counter,
            'completed_windows':self.completed,'capped_window_exports':self.window_exports,
            'all_windows_pass':bool(self.all_windows_pass),'supplies_pass':bool(supplies),'event_verdict':verdict,
            'comparison_scope':'Enabled joint event' if self.reactor.enabled else 'Disabled output schedule OR any resource exit credited as success; no entry/residence/supply requirement.'}


class ProbabilityBudget:
    def __init__(self,V=COPY_SCALE,H=POST_STARTUP_DURATION,parameters=Parameters()):
        self.V=V;self.H=Q(str(H));self.p=parameters
        if not isinstance(V,int) or V<100000000 or self.H<1 or not parameters.theorem_scope():
            raise ValueError('Published probability scope: integer V>=1e8, H>=1, K in [8,12], r in [18,22], delta in [1/100,1/20], fixed epsilon.')

    def logs(self,ctx=mp):
        V=ctx.mpf(self.V);H=number(ctx,self.H);T=500+H;m=self.H.numerator//self.H.denominator
        return {'entry':-ctx.mpf(8047)/312500000000*V,'startup_clock':-3*V/500,
            'resource_initial':ctx.log(4)-V/1000,
            'resource_leak':ctx.log(12000)+ctx.log(V)+ctx.log(T)+number(ctx,'1/50')-V/2000,
            'return_initial':-V/62500,
            'return_leak':ctx.log(3000)+ctx.log(V)+ctx.log(T)+number(ctx,'18/125')-V/15625,
            'all_windows':ctx.log(m)+number(ctx,'1/8')-ctx.mpf(13)*V/1080000,
            'two_foods':ctx.log(2)-(2*ctx.log(2)-1)*V*T,'gross_service':-V*T/200}

    def interval_pass(self,rho):
        rho=Q(str(rho))
        if not 0<rho<1:raise ValueError('Failure target must lie between zero and one.')
        return bool(logsum(iv,self.logs(iv).values()).b<=iv.log(number(iv,rho)).a)

    def evaluate(self):
        logs=self.logs();intervals=self.logs(iv);errorlog=logsum(mp,logs.values())
        in_simple=bool(iv.log(number(iv,self.H)).b<=(iv.mpf(self.V)/1000000).a)
        result={'V':self.V,'H':str(self.H),'complete_windows':self.H.numerator//self.H.denominator,
            'component_log10':{k:mp.nstr(v/mp.log(10),25) for k,v in logs.items()},
            'failure_formula':mp.nstr(mp.exp(errorlog),35),
            'enabled_success_lower_ordinary':mp.nstr(max(mp.mpf(0),1-mp.exp(errorlog)),35),
            'log_failure_outward_interval':str(logsum(iv,intervals.values())),
            'simplified_horizon_verified':in_simple,
            'simplified_success_lower':mp.nstr(1-2*mp.exp(-mp.mpf(self.V)/50000000),35) if in_simple else None,
            'log10_exponential_horizon':mp.nstr(mp.mpf(self.V)/1000000/mp.log(10),25),
            'disabled_log10_upper':None,'scope':'Outward mpmath interval evaluation is numerical certificate checking, not Lean evidence; no independent-window assumption.'}
        if self.H.denominator==1:
            off=[-mp.mpf(self.V)*self.H.numerator/25000,logs['resource_initial'],logs['resource_leak']]
            offiv=[-iv.mpf(self.V)*self.H.numerator/25000,intervals['resource_initial'],intervals['resource_leak']]
            result.update(disabled_log10_upper=mp.nstr(logsum(mp,off)/mp.log(10),35),
                disabled_log_upper_interval=str(logsum(iv,offiv)),
                disabled_below_10_power_minus_21699=bool(logsum(iv,offiv).b<(-21699*iv.log(10)).a))
        return result

    @classmethod
    def sufficient_scale(cls,H,rho,parameters=Parameters()):
        H=Q(str(H));rho=Q(str(rho))
        if H<1 or not 0<rho<1:raise ValueError('H>=1 and failure target in (0,1) required.')
        simple=int(mp.ceil(max(100000000,1000000*mp.log(number(mp,H)),50000000*mp.log(2/number(mp,rho)))))
        for _ in range(8):
            if bool(iv.log(number(iv,H)).b<=(iv.mpf(simple)/1000000).a) and bool(iv.log(2/number(iv,rho)).b<=(iv.mpf(simple)/50000000).a):break
            simple+=1
        else:raise ArithmeticError('Sizing rounding check unresolved.')
        candidate=100000000
        for _ in range(100):
            if cls(candidate,H,parameters).interval_pass(rho):
                return {'simple_sufficient_V':simple,'individually_interval_checked_V':candidate,'minimum_claimed':False}
            candidate=(11*candidate+9)//10
        raise RuntimeError('No checked geometric candidate in 100 steps.')


class ProductionLedger:
    def __init__(self,policy,parameters=Parameters()):self.policy=policy;self.p=parameters
    def guarantees(self):
        V=self.policy.V;T=self.policy.end;m=self.policy.windows
        return {'per_window_export_mass':str(Q(V,5000)),'total_completed_export_mass':str(Q(V*m,5000)),
            'covalent_equivalents_per_window':str(Q(V,20000)),
            'each_food_budget':str(2*V*T),'gross_drive_budget':str(V*T/8),
            'resident_mass_interval':[str(Q(18,5)*V),str(Q(22,5)*V)],
            'service_per_export_upper':str(625*T/m),'recovery_ratio_lower':str(Q(m,5000)/(4+8*T)),
            'scope':'Deterministic consequences conditional on the joint event; not expected physical efficiency.'}
    def aligned_interval(self,length):
        length=Q(str(length))
        if not 0<=length<=self.policy.windows:raise ValueError('Interval must fit within the completed-window span.')
        return Q(self.policy.V,5000)*max(0,math.floor(length)-1)
    def physical(self,Cref=REFERENCE_MOLAR,dilution=DILUTION_PER_HOUR,temperature=TEMPERATURE_K):
        Cref=number(mp,Cref);dilution=number(mp,dilution);temperature=number(mp,temperature)
        if min(Cref,dilution,temperature)<=0:raise ValueError('Positive dimensional scales required.')
        NA=mp.mpf('6.02214076e23');volume=self.policy.V/(NA*Cref)
        affinity=mp.log(number(mp,self.p.K/self.p.eta))
        return {'volume_pL':mp.nstr(volume*10**12,25),'startup_hours':mp.nstr(500/dilution,25),
            'post_startup_hours':mp.nstr(number(mp,self.policy.H)/dilution,25),
            'covalent_equivalent_rate_pM_per_hour':mp.nstr(Cref*dilution/20000*10**12,25),
            'reservoir_affinity_kBT':mp.nstr(affinity,25),
            'reservoir_affinity_kJ_per_mol':mp.nstr(mp.mpf('8.314462618')*temperature*affinity/1000,25),
            'scope':'Unit translation and maintained-reservoir requirement, not calibration; collected complexes are not purified free active product.'}


def algebraic_certificate():
    """Recompute finite source identities and exact worst-box coefficient slack."""
    import sympy as sp
    counts=sp.symbols('u w x c1 c2 z');V=100000000
    reactor=DrivenReactor();off=DrivenReactor(enabled=False)
    identities={}
    for name,weights in (('A',A_WEIGHTS),('B',B_WEIGHTS)):
        drift,noise=reactor.generator(counts,V,weights)
        identities[name+'_restoring_drift']=sp.expand(drift-(V-dot(weights,counts)))==0
        # Resource quadratic rate is bounded by V+2D because washout marks are 0,1,2.
        identities[name+'_quadratic_bound']=all(Q(c.marks[5 if name=='A' else 6])**2<=2*c.marks[5 if name=='A' else 6]
            for c in reactor.channels if 12<=c.label<=17)
    p=reactor.p;u,w,x,c1,c2,z=counts
    drift,noise=reactor.generator(counts,V,Y_WEIGHTS)
    expansion=(p.epsilon+p.delta*p.eta)*u*w/V-(p.epsilon/p.K+p.delta)*x+(Q(5,2)*u/V-1)*x+(-Q(29,8)+Q(11,2)*w/V)*c1+Q(11,10)*c2+(p.release/5-Q(9,5)-8/p.K)*z-p.release/(5*V)*(x*x-x)
    identities['Y_drift_expansion']=sp.expand(drift-expansion)==0
    identities['mass_balance_each_label']=all(dot(M_WEIGHTS,c.jump)+c.marks[7]==2*(c.marks[1]+c.marks[2]) for c in reactor.channels)
    identities['internal_moiety_conservation']=all(dot(A_WEIGHTS,c.jump)==dot(B_WEIGHTS,c.jump)==0 for c in reactor.channels if c.label<10 or c.label>=18)
    identities['disabled_pathwise_covalent_creation_bound']=all(Q(dot(S_WEIGHTS,c.jump)+c.marks[0],4)<=int(c.label in (0,19)) for c in off.channels)
    offdrift,_=off.generator(counts,V,S_WEIGHTS)
    export=sum(c.count_rate(counts,V)*c.marks[0] for c in off.channels)
    identities['disabled_generator_identity']=sp.expand(offdrift+export-4*((p.epsilon+p.delta*p.eta)*u*w/V-(p.epsilon/p.K+p.delta)*x))==0
    drift_coefficients=(1-Q(1,20)-Q(22,5000)-Q(1,8000000),Q(31,40),Q(11,10),Q(4,5))
    noise_coefficients=(1+Q(1,20)+Q(11,32)+Q(22,25000)+Q(1,8000000),Q(10374,3200),Q(2669,400),Q(113,25))
    weights=Y_WEIGHTS[2:]
    drift_slack=[a-Q(39,100)*w for a,w in zip(drift_coefficients,weights)]
    noise_slack=[5*w-a for a,w in zip(noise_coefficients,weights)]
    cE=399*Q(14,25)*Q(2,25)*Q(1,500000000)-Q(1,2500*40000)
    cW=Q(1,27000)-Q(1,40000);lam=Q(121,100)*Q(321,320)*Q(1,500000000)
    disabled_slack=Q(1,21000)-1004*lam-Q(1,25000)
    if not all(identities.values()) or min(*drift_slack,*noise_slack,disabled_slack)<=0:raise AssertionError('Algebraic certificate failed.')
    return {'source_identities':identities,'drift_coefficient_slack':list(map(str,drift_slack)),
        'noise_coefficient_slack':list(map(str,noise_slack)),'entry_exponent':str(cE),'window_exponent':str(cW),
        'disabled_creation_intensity_coefficient':str(lam),'disabled_chernoff_slack':str(disabled_slack),
        'scope':'Source identities and rational coefficient checks, not a replay of the finite-kernel or Lean proof.'}


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));parser.add_argument('--simulate',action='store_true')
    args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
    parameters=Parameters();policy=ObservationPolicy();reactor=DrivenReactor(parameters);off=DrivenReactor(parameters,False)
    budget=ProbabilityBudget(policy.V,policy.H,parameters);bounds=budget.evaluate();ledger=ProductionLedger(policy,parameters)
    result={'parameters':{n:str(getattr(parameters,n)) for n in parameters.__dataclass_fields__},
        'probability':bounds,'requested_failure':TARGET_FAILURE,'requested_failure_interval_pass':budget.interval_pass(TARGET_FAILURE),
        'sufficient_scales':budget.sufficient_scale(policy.H,TARGET_FAILURE,parameters),
        'event_guarantees':ledger.guarantees(),'physical_reading':ledger.physical(),
        'startup_no_seed_probability_lower_at_time_one':mp.nstr(mp.exp(-number(mp,parameters.epsilon+parameters.delta*parameters.eta)*policy.V),25),
        'claimed_0_9238_success_interval_check':budget.interval_pass('0.0762'),
        'algebraic_certificate':algebraic_certificate()}
    result['expected_completed_export_lower_ordinary']=mp.nstr(number(mp,Q(policy.V*policy.windows,5000))*mp.mpf(bounds['enabled_success_lower_ordinary']),35)
    if args.simulate:result['limited_count_path']=reactor.stochastic(policy)
    if policy.H>10000:raise ValueError('Default illustration requires H<=10000. Use ProbabilityBudget directly for larger durations.')
    enabled=reactor.deterministic(policy.H);disabled=off.deterministic(policy.H)
    independent=reactor.deterministic(policy.H,method='BDF')
    times=np.unique(np.r_[np.linspace(0,20,201),np.linspace(20,float(policy.end),581)])
    values=enabled.sol(times).T;offvalues=disabled.sol(times).T
    result['deterministic_illustration']={'entry_time':float(enabled.t_events[0][0]) if len(enabled.t_events[0]) else None,
        'max_independent_solver_difference':float(np.max(np.abs(values-independent.sol(times).T))),
        'scope':'Density-limit ODE with ledgers per V; does not estimate finite-copy startup or event probabilities.'}
    windows=[]
    for i in range(policy.windows):
        end=501+i;start=end-1
        windows.append([i,float(enabled.sol(end)[6]-enabled.sol(start)[6]),float(disabled.sol(end)[6]-disabled.sol(start)[6]),1/5000])
    def write_csv(name,headers,rows):
        with (out/name).open('w',newline='') as f:
            writer=csv.writer(f);writer.writerow(headers);writer.writerows(rows)
    write_csv('deterministic_trajectory.csv',['time',*SPECIES,*[n+'_per_V' for n in LEDGER]],[[t,*row] for t,row in zip(times,values)])
    write_csv('window_exports.csv',['window','enabled_export_per_V','disabled_export_per_V','certified_threshold_per_V'],windows)
    sweep=[]
    for V in np.linspace(1e8,4e8,81).astype(np.int64):
        b=ProbabilityBudget(int(V),policy.H,parameters).evaluate()
        sweep.append([int(V),float(b['enabled_success_lower_ordinary']),float(b['simplified_success_lower']) if b['simplified_success_lower'] is not None else '',b['disabled_log10_upper']])
    write_csv('copy_scale_bounds.csv',['V','full_success_lower','simplified_success_lower','disabled_log10_upper'],sweep)
    write_csv('failure_components.csv',['failure_mode','log10_upper_component'],bounds['component_log10'].items())
    (out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
    lines=[f'Finite-duration joint success lower bound: {bounds["enabled_success_lower_ordinary"]}.',
        f'Simplified bound: {bounds["simplified_success_lower"]}; disabled log10 upper: {bounds["disabled_log10_upper"]}.',
        f'Sufficient scales for failure {TARGET_FAILURE}: {result["sufficient_scales"]}.',
        f'Output and supply ledger: {result["event_guarantees"]}.',
        'Entry and shared failures are paid once; window errors are added without independence.',
        'ODE trajectories are illustrative; direct stochastic paths, if requested, may remain unfinished.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained');rows=np.array([r[:3] for r in sweep],float)
    axes[0].plot(rows[:,0]/1e8,rows[:,1],label='Full finite-duration bound')
    axes[0].plot(rows[:,0]/1e8,rows[:,2],'--',label='Simplified bound')
    axes[0].set(xlabel='Copy scale V / 1e8',ylabel='Joint success probability lower bound',title='Joint operating-success lower bound');axes[0].legend(fontsize=8)
    if policy.H.denominator==1:
        axes[1].plot(rows[:,0]/1e8,[float(r[3]) for r in sweep],color='#b95024')
        axes[1].set(xlabel='Copy scale V / 1e8',ylabel='log10 disabled probability upper bound',title='Disabled-reactor output upper bound')
    else:
        axes[1].axis('off');axes[1].text(.1,.5,'Disabled finite-duration theorem shown\nonly for integer H.',transform=axes[1].transAxes)
    for ax in axes:ax.grid(alpha=.2)
    fig.savefig(out/'reliability.png',dpi=180);fig.savefig(out/'reliability.svg');plt.close(fig)
    fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
    axes[0].semilogy(times,np.maximum(values[:,:6]@np.array(Y_WEIGHTS,float),1e-14),label='Enabled Y/V')
    axes[0].semilogy(times,np.maximum(offvalues[:,:6]@np.array(Y_WEIGHTS,float),1e-14),label='Disabled Y/V')
    axes[0].axhline(1/2500,color='black',ls='--',lw=.8,label='Entry threshold')
    axes[0].set(xlabel='Dimensionless time',ylabel='Weighted catalytic concentration',title='Catalyst buildup in the deterministic\nreactor');axes[0].legend(fontsize=8)
    w=np.array(windows,float)
    axes[1].semilogy(w[:,0]+1,w[:,1],label='Enabled')
    axes[1].semilogy(w[:,0]+1,w[:,2],label='Disabled')
    axes[1].axhline(1/5000,color='black',ls='--',label='Theorem threshold')
    axes[1].set(xlabel='Completed unit window after time 500',ylabel='Deterministic covalent export / V',title='Export in consecutive operating windows');axes[1].legend(fontsize=8)
    for ax in axes:ax.grid(alpha=.2)
    fig.savefig(out/'reactor.png',dpi=180);fig.savefig(out/'reactor.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    (out/'run_metadata.json').write_text(json.dumps({'manuscript_sha256':MANUSCRIPT_SHA256,
        'source_sha256':digest(Path(__file__)),'python':platform.python_version(),
        'output_sha256':{p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}},indent=2)+'\n')


if __name__=='__main__':main()