"""Molecular pulses, continued count histories and repeated-harvesting certificates."""
# EDITABLE INPUTS: paper-scale bounds, density illustration, and separate small-count path.
COPY_SCALE = 200_000_000_000
MISSION_CYCLES = 100
TARGET_FAILURE = '1/100'
RELEASE_SPEED = '20'
CLEAVAGE_SPEED = '3/100'
INITIAL_DENSITY = ('934/1000','935/1000','60/1000','1/1000','1/1000','1/1000')
RETENTIONS = ('1/4','3/4','2/5','3/5')
SURVIVAL = ('49/50','1','49/50','1','99/100','1')
REFILL_ERRORS = ('-1/200','1/200')
DENSITY_CYCLES = 12
SMALL_COPY_SCALE = 1000
SMALL_CYCLES = 3
SMALL_EVENT_BUDGET = 250000
RANDOM_SEED = 43092026
REFERENCE_MOLAR = '1/1000'
TIME_UNIT_SECONDS = '60'

from dataclasses import dataclass,asdict
from fractions import Fraction as Q
from pathlib import Path
from itertools import product
import argparse,copy,csv,hashlib,json,math,platform
import numpy as np
import sympy as sp
from mpmath import mp,iv
from scipy.integrate import solve_ivp

mp.dps=70;iv.dps=70
MANUSCRIPT_SHA256='235f5a2a613e48b9059f90bf38afa0b502225d32956e541096cd3141b2ea4db2'
SPECIES=('U','W','X','C1','C2','Z')
A=(1,0,1,2,2,2);B=(0,1,1,1,2,2);I=(0,0,1,1,1,2);Y=(Q(0),Q(0),Q(1),Q(9,8),Q(7,5),Q(9,5))
LEDGER=('collected_I','collected_X','food_U','food_W','gross_drive','all_wash_I','net_synthesis')
V0=200000000000
def dot(a,b):return sum(x*y for x,y in zip(a,b))
def num(ctx,x):
    x=Q(str(x));return ctx.mpf(x.numerator)/x.denominator
def ceildiv(n,d):return -((-n)//d)
def counts(values):
    values=tuple(values)
    if len(values)!=6 or any(type(n) is not int or n<0 for n in values):raise ValueError('Six nonnegative Python integer counts required.')
    return values
def restart(N,V):
    N=counts(N)
    return 159*V<=160*dot(A,N)<=161*V and 159*V<=160*dot(B,N)<=161*V and 40*dot(Y,N)>=2*V

@dataclass(frozen=True)
class Channel:
    label:int
    inputs:tuple
    outputs:tuple
    coefficient:Q
    @property
    def jump(self):return tuple(self.outputs.count(i)-self.inputs.count(i) for i in range(6))
    def propensity(self,N,V):
        value=self.coefficient*Q(V)**(1-len(self.inputs));used={}
        for i in self.inputs:value*=N[i]-used.get(i,0);used[i]=used.get(i,0)+1
        return value
    def marks(self,collect):
        wash=I[self.label-12] if 12<=self.label<=17 else 0
        return (wash if collect else 0,int(collect and self.label==14),int(self.label==10),int(self.label==11),int(self.label in (18,19)),wash,dot(I,self.jump) if self.label<10 or self.label>=18 else 0)

class Reactor:
    """Literal twenty-label source; count and density propensities share stoichiometry."""
    def __init__(self,release=RELEASE_SPEED,cleavage=CLEAVAGE_SPEED):
        self.r=Q(release);self.d=Q(cleavage)
        if self.r<=0 or self.d<=0:raise ValueError('Positive speeds required.')
        eps=Q(1,500000000);eta=Q(1,8000000000)
        pairs=[((0,1),(2,),eps,eps/10),((2,0),(3,),Q(20),Q(20)),((3,1),(4,),Q(20),Q(20)),((4,),(5,),Q(20),Q(2)),((5,),(2,2),self.r,self.r)]
        channels=[]
        for j,(left,right,kf,kr) in enumerate(pairs):channels.extend([Channel(2*j,left,right,kf),Channel(2*j+1,right,left,kr)])
        channels.extend(Channel(10+i,(),(i,),Q(1)) for i in (0,1))
        channels.extend(Channel(12+i,(i,),(),Q(1)) for i in range(6))
        channels.extend([Channel(18,(2,),(0,1),self.d),Channel(19,(0,1),(2,),self.d*eta)])
        self.channels=tuple(channels);self.jumps=np.array([c.jump for c in channels]);self.k=np.array([float(c.coefficient) for c in channels]);self.orders=np.array([len(c.inputs) for c in channels])
        self.inputs=np.full((20,2),6,int);self.offsets=np.zeros((20,2),int)
        for j,c in enumerate(channels):
            for k,i in enumerate(c.inputs):self.inputs[j,k]=i;self.offsets[j,k]=c.inputs[:k].count(i)
    def scope(self):return 19<=self.r<=21 and Q(1,50)<=self.d<=Q(1,25)
    def rates(self,state,V=None):
        factors=np.r_[state,1][self.inputs]
        if V is not None:factors=factors-self.offsets
        rates=self.k*np.prod(factors,axis=1)
        return rates if V is None else rates*np.power(float(V),1-self.orders)
    def symbolic_identities(self):
        N=sp.symbols('u w x c1 c2 z');V=sp.Symbol('V',positive=True)
        # Use the literal channels; independent concentration drift uses ordinary powers.
        rate=[c.coefficient*V**(1-len(c.inputs))*sp.prod(N[i]-c.inputs[:k].count(i) for k,i in enumerate(c.inputs)) for c in self.channels]
        density=[c.coefficient*sp.prod(N[i]/V for i in c.inputs) for c in self.channels]
        L=lambda weights:sp.expand(sum(a*dot(weights,c.jump) for a,c in zip(rate,self.channels)))
        corrections=[sp.expand(L(A)-(V-dot(A,N))),sp.expand(L(B)-(V-dot(B,N))),sp.expand(L(Y)-V*sum(a*dot(Y,c.jump) for a,c in zip(density,self.channels))-self.r*N[2]/(5*V))]
        if any(c!=0 for c in corrections):raise ArithmeticError('Generator identity failed.')
        return dict(material_A='V-A',material_B='V-B',stock_correction=str(self.r)+'*N_X/(5V)',chemical_inventory_marks=[c.marks(False)[6] for c in self.channels],identities_exact=True)

@dataclass(frozen=True)
class Intervention:
    q:object='1/4'
    survival:tuple=SURVIVAL
    refill:tuple=REFILL_ERRORS
    def __post_init__(self):
        object.__setattr__(self,'q',Q(self.q));object.__setattr__(self,'survival',tuple(map(Q,self.survival)));object.__setattr__(self,'refill',tuple(map(Q,self.refill)))
        if not Q(1,4)<=self.q<=Q(3,4) or len(self.survival)!=6 or any(not Q(49,50)<=l<=1 for l in self.survival) or len(self.refill)!=2 or any(abs(e)>Q(1,200) for e in self.refill):raise ValueError('Intervention outside the admitted box.')
    def probabilities(self,i):return self.q*self.survival[i],1-self.q,self.q*(1-self.survival[i])
    def doses(self,V):return tuple(math.floor(V*(1-self.q+e)) for e in self.refill)

@dataclass(frozen=True)
class PulseOutcome:
    retained:tuple
    withdrawn:tuple
    lost:tuple
    doses:tuple
    @property
    def start(self):return tuple(n+(self.doses[i] if i<2 else 0) for i,n in enumerate(self.retained))

class MolecularPulse:
    def sample(self,N,V,intervention,rng):
        N=counts(N)
        if type(V) is not int or V<=0:raise ValueError('Positive integer copy scale required.')
        if max(N)>=2**53:raise ValueError('Numerical multinomial sampler limited to counts below 2^53.')
        draws=[tuple(map(int,rng.multinomial(n,list(map(float,intervention.probabilities(i)))))) for i,n in enumerate(N)]
        return PulseOutcome(*(tuple(row[k] for row in draws) for k in range(3)),intervention.doses(V))
    def enumerate(self,N,V,intervention):
        """Exact rational product law for small pedagogical populations; no renormalization."""
        N=counts(N)
        if type(V) is not int or V<=0:raise ValueError('Positive integer copy scale required.')
        if sum(N)>12:raise ValueError('Exact enumeration is restricted to at most twelve molecules.')
        species=[]
        for i,n in enumerate(N):
            pr,pw,pl=intervention.probabilities(i);terms=[]
            for nr in range(n+1):
                for nw in range(n-nr+1):
                    nl=n-nr-nw;mass=math.comb(n,nr)*math.comb(n-nr,nw)*pr**nr*pw**nw*pl**nl
                    if mass:terms.append(((nr,nw,nl),mass))
            species.append(terms)
        for terms in product(*species):
            triples=[t[0] for t in terms];yield PulseOutcome(*(tuple(t[k] for t in triples) for k in range(3)),intervention.doses(V)),math.prod(t[1] for t in terms)

class CycleLedger:
    def __init__(self,original,pulse,V):
        self.original=counts(original);self.pulse=pulse;self.state=list(pulse.start);self.V=V;self.values=[0,0,*pulse.doses,0,0,0]
    def event(self,time,channel):
        if channel.propensity(self.state,self.V)<=0:raise ValueError('Unavailable reaction.')
        self.state=[n+d for n,d in zip(self.state,channel.jump)]
        self.values=[n+d for n,d in zip(self.values,channel.marks(time>=3))]
        if min(self.state)<0:raise ArithmeticError('Negative molecule count.')
    def audit(self):
        lhs=self.values[6];rhs=dot(I,self.state)-dot(I,self.original)+self.values[5]+dot(I,self.pulse.withdrawn)+dot(I,self.pulse.lost)
        if lhs!=rhs:raise ArithmeticError('Pathwise inventory mismatch.')
        return lhs
    def record(self,complete,elapsed,events):
        self.audit();v=self.values;success=restart(tuple(self.state),self.V) and v[0]>=ceildiv(self.V,56) and v[1]>=ceildiv(self.V,1080) and v[2]<=5*self.V and v[3]<=5*self.V and v[4]<=self.V//5
        return dict(start=list(self.original),pulse=asdict(self.pulse),endpoint=self.state.copy(),counters=dict(zip(LEDGER,v)),complete=complete,elapsed=elapsed,events=events,success=bool(success) if complete else None,withdrawn_I=dot(I,self.pulse.withdrawn),lost_I=dot(I,self.pulse.lost))

class StochasticCycle:
    def __init__(self,reactor,pulse=None):self.reactor=reactor;self.pulse=pulse or MolecularPulse()
    def run(self,N,V,intervention,rng,event_budget):
        ledger=CycleLedger(N,self.pulse.sample(N,V,intervention,rng),V);t=0.
        for event in range(event_budget):
            rates=self.reactor.rates(ledger.state,V);total=float(rates.sum());next_time=t+float(rng.exponential(1/total))
            if next_time>4:return ledger.record(True,4.,event)
            if next_time<=t:raise ArithmeticError('Time precision exhausted at this scale; reduce V for sampling.')
            label=min(int(np.searchsorted(np.cumsum(rates),rng.random()*total,side='right')),19)
            t=next_time;ledger.event(t,self.reactor.channels[label])
        return ledger.record(False,t,event_budget)

class HistoryMission:
    """Controllers receive all completed records, including failures; actual returns persist."""
    def __init__(self,cycle):self.cycle=cycle
    def run(self,initial,V,m,controller,seed,event_budget):
        if type(V) is not int or V<=0 or type(m) is not int or m<0 or type(event_budget) is not int or event_budget<0:raise ValueError('Integer scale, cycle count and budget required.')
        state=counts(initial);rng=np.random.default_rng(seed);history=[];used=0
        for k in range(m):
            intervention=controller(copy.deepcopy(tuple(history)));record=self.cycle.run(state,V,intervention,rng,event_budget-used);history.append(record);used+=record['events']
            if not record['complete']:break
            state=tuple(record['endpoint'])
        completed=len(history)==m and all(r['complete'] for r in history)
        return dict(history=history,complete=completed,all_success=all(r['success'] for r in history) if completed else None,total_events=used,initial_in_restart=restart(initial,V))

class ReturnedStockController:
    def __call__(self,history):
        # React to actual previous yield; all choices remain in the admitted box.
        if history and history[-1]['counters']['collected_X']==0:q='3/4'
        else:q=RETENTIONS[len(history)%len(RETENTIONS)]
        return Intervention(q)

class DensityExperiment:
    """Mean molecular pulse followed by full mass-action ODE; not a probability estimate."""
    def __init__(self,reactor):self.reactor=reactor
    def run(self,initial,cycles,method='Radau'):
        state=np.array([float(Q(x)) for x in initial]);history=[];trajectory=[]
        if state.shape!=(6,) or min(state)<0:raise ValueError('Six nonnegative concentrations required.')
        for k in range(cycles):
            intervention=Intervention(RETENTIONS[k%len(RETENTIONS)]);before=state.copy();retained=np.array([float(intervention.probabilities(i)[0]) for i in range(6)])*state
            withdrawn=state*float(1-intervention.q);lost=state*np.array([float(intervention.probabilities(i)[2]) for i in range(6)])
            doses=np.array([float(1-intervention.q+e) for e in intervention.refill]);state=retained.copy();state[:2]+=doses
            augmented=np.r_[state,[0,0,*doses,0,0,0]]
            for start,end,collect in [(0,3,False),(3,4,True)]:
                marks=np.array([c.marks(collect) for c in self.reactor.channels],float)
                def rhs(t,y):
                    rates=self.reactor.rates(y[:6]);return np.r_[rates@self.reactor.jumps,rates@marks]
                sol=solve_ivp(rhs,(start,end),augmented,method=method,rtol=2e-10,atol=1e-12,dense_output=True)
                if not sol.success:raise ArithmeticError(sol.message)
                times=np.linspace(start,end,int(40*(end-start))+1);trajectory.extend((k,float(t),*row) for t,row in zip(times,sol.sol(times).T));augmented=sol.y[:,-1]
            state=augmented[:6];v=augmented[6:];residual=v[6]-(dot(I,state)-dot(I,before)+v[5]+dot(I,withdrawn)+dot(I,lost))
            history.append(dict(cycle=k,q=str(intervention.q),endpoint=state.tolist(),counters=dict(zip(LEDGER,v.tolist())),inventory_residual=float(residual)))
        return dict(history=history,trajectory=trajectory)

class MissionCertificate:
    def __init__(self,V,reactor):
        if type(V) is not int or V<1000000 or not reactor.scope():raise ValueError('One-cycle theorem requires integer V>=1e6 and fixed paired speeds in the paper rectangle.')
        self.V=V
    def logs(self,ctx=mp):
        V=ctx.mpf(self.V);ln=ctx.log
        return {'pulse_stock':-num(ctx,'1177/1000000000')*V,'pulse_material':ln(4)-V/100000,
            'free_counter':-V/40000000,'phase_occupation':ln(100)-V/10000000000,'free_clock_early':-V,'free_clock_collection':-V/200,
            'template_counter':-V/100000,'foods':ln(2)-V/300,'gross_service':-V/2000,
            'material_initial_twice':ln(8)-V/2000,'material_leak_twice':ln(96000)+ln(V)-V/2000+num(ctx,'1/50'),
            'recovery':-3*V/5000000,'recovery_clock':-V/2,'residence_initial':-V/10000,'residence_leak':ln(12000)+ln(V)-V/10000+num(ctx,'9/500'),
            'terminal_material':ln(4)-V/320000,'terminal_clock':ln(4)-V}
    def log_error(self,ctx=mp):
        logs=self.logs(ctx);shift=logs['phase_occupation'];return shift+ctx.log(sum(ctx.exp(x-shift) for x in logs.values()))
    def evaluate(self,m):
        if type(m) is not int or m<0:raise ValueError('Nonnegative integer mission length required.')
        L=self.log_error();LI=self.log_error(iv);e=mp.exp(L);eI=iv.exp(LI)
        if m==0:lower=iv.mpf(1);ordinary=mp.mpf(1)
        elif bool(eI.b<iv.mpf(1).a):lower=(1-eI)**m;ordinary=mp.exp(m*mp.log1p(-e))
        else:lower=iv.mpf(0);ordinary=mp.mpf(0)
        digits=max(0,int(mp.floor(mp.mpf(lower.a)*10**9)))
        # Certify the printed decimal by an interval endpoint comparison.
        while digits and not bool(num(iv,Q(digits,10**9)).b<=lower.a):digits-=1
        return dict(V=self.V,m=m,error_log10=mp.nstr(L/mp.log(10),25),error_log_interval=str(LI),joint_product_lower_ordinary=mp.nstr(ordinary,25),joint_product_lower_certified=f'{digits/10**9:.9f}',linear_lower_ordinary=mp.nstr(max(0,1-m*e),25),component_log10={k:mp.nstr(v/mp.log(10),22) for k,v in self.logs().items()},scope='Repeated bound conditions on actual successful returns; it does not assume independent cycles. Initial restart and admitted interventions are required.')
    @staticmethod
    def sizing(m,failure,template_demand=0,free_demand=0):
        failure=Q(str(failure))
        if type(m) is not int or m<1 or not 0<failure<1 or any(type(n) is not int or n<0 for n in (template_demand,free_demand)):raise ValueError('Positive mission length, failure in (0,1), nonnegative integer demands required.')
        value=iv.mpf(10**10)*iv.log(101*m/num(iv,failure));V=max(V0,int(mp.ceil(mp.mpf(value.b))),ceildiv(56*template_demand,m),ceildiv(1080*free_demand,m))
        # A rational Taylor lower sum also certifies the selected logarithmic size.
        x=Q(V,10**10);target=101*m/failure;term=Q(1);total=Q(1);degree=0
        while total<target:
            degree+=1;term*=x/degree;total+=term
            if degree>10000:raise ArithmeticError('Taylor certificate budget exhausted.')
        necessary=iv.mpf(10**10)*iv.log(100*m/num(iv,failure))
        return dict(sufficient_V=V,taylor_degree=degree,taylor_sum=str(total),taylor_target=str(target),earlier_linear_V=max(V0,math.ceil(5000000*m/failure)),certificate_bottleneck_lower_ordinary=mp.nstr(10**10*mp.log(100*m/num(mp,failure)),25),certificate_bottleneck_interval=str(necessary))
    def inventory(self,m,initial=None):
        if type(m) is not int or m<1:raise ValueError('Positive mission length required for final-restart inventory bound.')
        if initial is not None and not restart(initial,self.V):raise ValueError('Initial state outside restart set.')
        V=self.V;initial_I=dot(I,initial) if initial is not None else 161*V//160
        return dict(collected_I=m*ceildiv(V,56),collected_X=m*ceildiv(V,1080),each_food=5*m*V,gross_drive=m*(V//5),final_I=ceildiv(V,28),initial_I=initial_I,net_synthesis=m*ceildiv(V,56)+ceildiv(V,28)-initial_I)

class PhaseTransport:
    def weights(self,V,n):
        if type(V) is not int or V<1 or type(n) is not int or n<0:raise ValueError('Integer V>=1 and n>=0 required.')
        q=mp.mpf(3000)*V;a=1-70/q;base=mp.exp(n*mp.log1p(-70/q))
        return tuple(base*c for c in (1,20*n/(q*a),580*n*(n-1)/(q*q*a*a),38*n/(q*a)))
    def certificates(self):
        q,n=sp.symbols('q n',positive=True);a=1-70/q
        closed=lambda n:sp.Matrix([a**n,a**n*20*n/(q*a),a**n*580*n*(n-1)/(q*q*a*a),a**n*38*n/(q*a)])
        T=sp.Matrix([[a,0,0,0],[20/q,a,0,0],[0,20/q,a,20/q],[38/q,0,0,a]])
        if any(sp.simplify(x)!=0 for x in T*closed(n)-closed(n+1)):raise ArithmeticError('Adjoint closed form failed.')
        lower=(Q(1,12),Q(20*89,36000),Q(580*89*88,108000000),Q(38*89,36000))
        if any(l<w/35 for l,w in zip(lower,Y[2:])):raise ArithmeticError('Phase-window coefficient failed.')
        # Bound exp(11/5) above by a Taylor sum and geometric remainder.
        x=Q(11,5);s=sum(x**k/math.factorial(k) for k in range(21));tail=x**21/math.factorial(21)/(1-x/22)
        if s+tail>12:raise ArithmeticError('Exponential upper certificate failed.')
        if Q(6370,2930)>Q(11,5):raise ArithmeticError('Worst-case clock-window logarithmic bound failed.')
        return dict(symbolic_recurrence=True,window_coordinate_lower=list(map(str,lower)),target=list(map(lambda w:str(w/35),Y[2:])),exp_upper=str(s+tail),phase_exponent=str(Q(1,10**9)-Q(1,700000000)+Q(12,5)*91/Q(10**12)))

def stopping_boundary_example():
    """A separate three-state counterexample to transferring a closed stopped event."""
    # States: active, boundary, outside. Row-stochastic kernels start at active.
    physical=sp.Matrix([[Q(1,2),Q(1,2),0],[0,0,1],[0,0,1]])
    stopped=sp.Matrix([[Q(1,2),Q(1,2),0],[0,1,0],[0,0,1]])
    a=(physical**2)[0,:];b=(stopped**2)[0,:]
    return dict(physical_closed=str(a[0]+a[1]),stopped_closed=str(b[0]+b[1]),physical_strict=str(a[0]),stopped_strict=str(b[0]),scope='Separate illustrative chain, not a reactor counterexample or a replay of chronological-law identification.')

def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
    reactor=Reactor();cert=MissionCertificate(COPY_SCALE,reactor);rich=tuple(int(Q(c)*COPY_SCALE) for c in INITIAL_DENSITY)
    if any(Q(c)*COPY_SCALE!=n for c,n in zip(INITIAL_DENSITY,rich)) or not restart(rich,COPY_SCALE):raise ValueError('Default initial density must map to an integer admitted restart at selected scale.')
    density=DensityExperiment(reactor).run(INITIAL_DENSITY,DENSITY_CYCLES);check=DensityExperiment(reactor).run(INITIAL_DENSITY,DENSITY_CYCLES,'BDF')
    small_initial=(0,0,SMALL_COPY_SCALE,0,0,0);sample=HistoryMission(StochasticCycle(reactor)).run(small_initial,SMALL_COPY_SCALE,SMALL_CYCLES,ReturnedStockController(),RANDOM_SEED,SMALL_EVENT_BUDGET)
    pulse=MolecularPulse().sample(rich,COPY_SCALE,Intervention(),np.random.default_rng(RANDOM_SEED))
    sizing=[dict(m=m,**cert.sizing(m,TARGET_FAILURE)) for m in (100,10000,1000000)]
    budgets=[MissionCertificate(V,reactor).evaluate(MISSION_CYCLES) for V in (V0,230359000000,300000000000)]
    inventory=cert.inventory(MISSION_CYCLES);rich_inventory=cert.inventory(MISSION_CYCLES,rich)
    NA=Q(602214076000000000000000);c=Q(REFERENCE_MOLAR);tau=Q(TIME_UNIT_SECONDS)
    if min(c,tau)<=0:raise ValueError('Positive concentration/time units required.')
    physical=dict(volume_nL=float(Q(COPY_SCALE)/NA/c*10**9),source_hours=float(4*MISSION_CYCLES*tau/3600),template_pmol=float(Q(inventory['collected_I'])/NA*10**12),free_pmol=float(Q(inventory['collected_X'])/NA*10**12),scope='Mixed effluent; free X is included in template equivalents. Time excludes preparation/handling. Supply allowances are not total energy costs.')
    result=dict(generator=reactor.symbolic_identities(),phase=PhaseTransport().certificates(),stopping_boundary=stopping_boundary_example(),mission=cert.evaluate(MISSION_CYCLES),sizing=sizing,worked_budgets=budgets,inventory_uniform=inventory,inventory_initial_specific=rich_inventory,physical=physical,
        sampled_pulse=asdict(pulse),small_count_path=sample,density_inventory_residual=max(abs(r['inventory_residual']) for r in density['history']),density_solver_difference=float(np.max(abs(np.array(density['trajectory'])-np.array(check['trajectory'])))),
        scope='Chemical/algebraic identities and arithmetic are freshly checked. Operating probability and process-identification theorems are imported; Lean is not rerun. Density and one small-count trajectory do not validate those probabilities.')
    def dump(name,data):(out/name).write_text(json.dumps(data,indent=2,default=str)+'\n')
    def table(name,header,rows):
        with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    dump('results.json',result)
    table('density_trajectory.csv',['cycle','local_time',*SPECIES,*[name+'_per_V' for name in LEDGER]],density['trajectory'])
    table('density_cycles.csv',['cycle','q',*SPECIES,*[name+'_per_V' for name in LEDGER],'inventory_residual'],[(r['cycle'],r['q'],*r['endpoint'],*r['counters'].values(),r['inventory_residual']) for r in density['history']])
    table('channels.csv',['label','reactants','products','coefficient','synthesis_mark'],[(c.label,','.join(SPECIES[i] for i in c.inputs),','.join(SPECIES[i] for i in c.outputs),str(c.coefficient),c.marks(False)[6]) for c in reactor.channels])
    sweep=[]
    for V in np.linspace(V0,4*10**11,61).astype(np.int64):
        exact=MissionCertificate(int(V),reactor).log_error();single=mp.log(101)-mp.mpf(int(V))/10**10;residual=mp.log(10**6)+mp.log(int(V)+1)-mp.mpf(int(V))/40000000
        sweep.append((int(V),float(exact/mp.log(10)),float(single/mp.log(10)),float(residual/mp.log(10))))
    table('error_budget.csv',['V','log10_full_error','log10_single_envelope','log10_residual_envelope'],sweep)
    phase=[(n/100,*map(float,PhaseTransport().weights(100,n))) for n in range(0,15001,50)]
    table('phase_weights.csv',['clock_steps_per_V','X','C1','C2','Z'],phase)
    lines=[f'Exact twenty-channel drift and phase-recurrence identities checked.',f'{MISSION_CYCLES} cycles at V={COPY_SCALE}: joint success >= {result["mission"]["joint_product_lower_certified"]}.',f'Uniform net synthesis lower bound: {inventory["net_synthesis"]}; initial-specific: {rich_inventory["net_synthesis"]}.',f'Collected template equivalents >= {inventory["collected_I"]}, including free X >= {inventory["collected_X"]}.',f'Small-count diagnostic: complete={sample["complete"]}, all_success={sample["all_success"]}, events={sample["total_events"]}.', 'The small-count path is outside the useful theorem scale. Failed cycles remain in its history; an event-budget cutoff has no success verdict.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained');data=np.array(sweep)
    axs[0].plot(data[:,0]/1e11,data[:,1],label='Full error budget');axs[0].plot(data[:,0]/1e11,data[:,2],ls='--',label='Single exponential');axs[0].set(xlabel='Copy scale V (hundreds of billions)',ylabel='log10 one-cycle failure upper bound',title='One-cycle failure bound versus\nmolecule-count scale');axs[0].legend(fontsize=8)
    p=np.array(phase)
    for i,label in enumerate(SPECIES[2:]):axs[1].plot(p[:,0],p[:,i+1]/float(Y[i+2]),label=label)
    axs[1].axhline(1/35,color='gray',ls='--');axs[1].axvspan(89,91,color='gray',alpha=.2);axs[1].set(xlabel='Uniformized clock steps / V',ylabel='Transport weight / stock weight',xlim=(70,110),ylim=(0,.25),title='Phase transport near the certified window');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'certificate.png',dpi=180);fig.savefig(out/'certificate.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained');data=np.array(density['trajectory']);t=4*data[:,0]+data[:,1]
    axs[0].plot(t,data[:,2:8]@np.array(Y,float),label='Weighted stock');axs[0].plot(t,data[:,4],label='Free template');axs[0].axhline(.05,color='gray',ls='--',label='Restart stock floor');axs[0].set(xlabel='Source-running time',ylabel='Concentration',title='Mean-pulse density illustration');axs[0].legend(fontsize=8)
    records=density['history'];x=np.arange(1,len(records)+1)
    axs[1].plot(x,[r['counters']['collected_I'] for r in records],'o-',label='Template equivalents');axs[1].plot(x,[r['counters']['collected_X'] for r in records],'o-',label='Included free X');axs[1].axhline(1/56,color='gray',ls='--');axs[1].axhline(1/1080,color='gray',ls=':');axs[1].set(xlabel='Cycle',ylabel='Collected output / V',title='Collected free template and total template\nequivalents');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'operation.png',dpi=180);fig.savefig(out/'operation.svg');plt.close(fig)
    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(),output_sha256={p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}))

if __name__=='__main__':main()