"""Common conservative transport, pulsed recovery, and honest production accounts.

Run python example.py --output outputs. Normalized deterministic concentrations;
the network theorem is imported, while its algebra and constants are checked here.
"""
from __future__ import annotations
import argparse
import csv
from dataclasses import dataclass
from fractions import Fraction as Q
import hashlib
import json
import math
from pathlib import Path
import platform
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import expm
import sympy as sp

# EDITABLE INPUTS -----------------------------------------------------------
TOPOLOGY = 'path'                       # uncoupled, path, ring, star (four nodes)
EXCHANGE_SCALE = 1.0                    # common to ALL species; no theorem upper cap
RELEASE = (19., 21., 19., 21.)          # fixed node rates, each in [19,21]
CLEAVAGE = (.04, .02, .02, .04)         # each in [.02,.04]
INITIAL_STATE = (
    (151/160,151/160,1/20,0,0,0),
    (1303/1440,1367/1440,0,2/45,0,0),
    (1033/1120,1033/1120,0,0,1/28,0),
    (1351/1440,1351/1440,0,0,0,1/36),
)
ROUTINE_CYCLES = 32                     # conditioning pulse + 12, then pulse + 4
BASE_RETENTION = (.25, .75, .4, .6)
SURVIVAL = (.98, 1., .98, 1., .99, 1.)
FOOD_ERROR = (-.005, -.005)
CONCENTRATION_MOLAR = .001              # illustrative, not fitted chemistry
NODE_VOLUME_LITERS = .001
TIME_SECONDS = 60.
DEMAND_I, DEMAND_X = '3/2', '1/5'       # normalized network amounts; overlapping outputs
AVAILABLE_TIME = 64
FOOD_U_BUDGET, FOOD_W_BUDGET, SERVICE_BUDGET = 300, 280, 12
RTOL, ATOL = 1e-9, 1e-12
MANUSCRIPT_SHA256 = 'c8b16eaf2abf49162495bbfb0ea0e69908a3a28d5e34968d60dcd399daf6518f'
# --------------------------------------------------------------------------
SPECIES = ('U','W','X','C1','C2','Z')
A = np.array([1,0,1,2,2,2.]); B = np.array([0,1,1,1,2,2.])
Y = np.array([0,0,1,9/8,7/5,9/5]); I = np.array([0,0,1,1,1,2.])
WEIGHTS = np.array([A,B,Y,I])
STOICH = np.array([[-1,-1,0,0,0,1],[-1,0,-1,0,0,1],
    [1,-1,0,0,2,-1],[0,1,-1,0,0,0],[0,0,1,-1,0,0],[0,0,0,1,-1,0]])
EPS, ETA = Q(1,500000000), Q(1,8000000000)


class LocalChemistry:
    """Six reversible pairs at each node; vectorized mass action, unit flow."""
    def __init__(self, release=RELEASE, cleavage=CLEAVAGE):
        self.r=np.asarray(release,float).copy(); self.d=np.asarray(cleavage,float).copy()
        if self.r.ndim!=1 or self.r.shape!=self.d.shape or not len(self.r) or not np.all(np.isfinite([self.r,self.d])) or np.any(self.r<0) or np.any(self.d<0):
            raise ValueError('Equal nonempty vectors of finite nonnegative node rates required.')

    @property
    def certified(self): return bool(np.all((self.r>=19)&(self.r<=21)&(self.d>=.02)&(self.d<=.04)))

    def fluxes(self,c):
        u,w,x,c1,c2,z=c.T
        return np.array([float(EPS)*(u*w-x/10),20*(x*u-c1),20*(c1*w-c2),
                         20*c2-2*z,self.r*(z-x*x),self.d*(x-float(ETA)*u*w)]).T

    def field(self,c): return np.array([1,1,0,0,0,0])-c+self.fluxes(c)@STOICH.T

    def service(self,c): return self.d*(c[:,2]+float(ETA)*c[:,0]*c[:,1])


class CommonExchange:
    """Equal-volume symmetric exchange; arbitrary finite graph, same D for all species."""
    def __init__(self, weights):
        k=np.asarray(weights,float).copy()
        if k.ndim!=2 or not len(k) or k.shape[0]!=k.shape[1] or not np.all(np.isfinite(k)) or np.any(k<0) or np.any(np.diag(k)!=0) or not np.array_equal(k,k.T):
            raise ValueError('Finite symmetric nonnegative zero-diagonal matrix required.')
        self.k=k; self.D=k-np.diag(k.sum(axis=1)); self.n=len(k)

    def field(self,c): return self.D@c
    def gross_handling(self,c): return float(self.k.sum(axis=1)@c.sum(axis=1))
    def handling_bound(self,duration): return 11/5*duration*self.k.sum()
    def directed_labels(self):
        return [(i,j,s,float(self.k[i,j])) for i in range(self.n) for j in range(self.n)
                if self.k[i,j]>0 for s in SPECIES]

    @classmethod
    def four_nodes(cls,name,scale=1.):
        if not np.isfinite(scale) or scale<0: raise ValueError('Nonnegative finite scale required.')
        edges={'uncoupled':[], 'path':[(0,1),(1,2),(2,3)],
               'ring':[(0,1),(1,2),(2,3),(3,0)],'star':[(0,1),(0,2),(0,3)]}[name]
        k=np.zeros((4,4))
        for i,j in edges:k[i,j]=k[j,i]=scale*(1/3 if name=='star' else .5)
        return cls(k)


@dataclass(frozen=True)
class Pulse:
    retained: tuple
    survival: tuple = SURVIVAL
    food_error: tuple = FOOD_ERROR

    def arrays(self,n):
        q=np.asarray(self.retained,float); ell=np.asarray(self.survival,float)
        e=np.asarray(self.food_error,float)
        if q.shape!=(n,) or ell.shape not in ((6,),(n,6)) or e.shape not in ((2,),(n,2)):
            raise ValueError('One retention per node; survival/refill may be common or node-specific.')
        if not all(np.all(np.isfinite(v)) for v in (q,ell,e)) or np.any((q<.25)|(q>.75)) or np.any((ell<.98)|(ell>1)) or np.any(abs(e)>.005):
            raise ValueError('Pulse outside paper bounds.')
        return q[:,None],np.broadcast_to(ell,(n,6)),np.broadcast_to(e,(n,2))

    def apply(self,c):
        q,ell,e=self.arrays(len(c)); after=q*ell*c; food=1-q+e
        after[:,:2]+=food
        return after,(1-q)*c,q*(1-ell)*c,food


class HistoryPolicy:
    """Example admissible feedback: preserve more after a relatively small last output.

    Receives snapshots of completed cycles and current state; never changes kinetics.
    Replace choose() to investigate another synchronized policy.
    """
    def choose(self,c,history):
        q=np.roll(np.asarray(BASE_RETENTION),len(history)).copy()
        if history:
            last=np.array(history[-1]['collection_I']); q[last<np.median(last)]=.75
        return Pulse(tuple(q))


class ReactorNetwork:
    def __init__(self,chemistry,exchange):
        if len(chemistry.r)!=exchange.n: raise ValueError('Rate and graph node counts differ.')
        self.chemistry=chemistry; self.exchange=exchange; self.n=exchange.n

    def validate_state(self,c):
        c=np.asarray(c,float).copy()
        if c.shape!=(self.n,6) or not np.all(np.isfinite(c)) or np.any(c<0): raise ValueError('Nonnegative finite N by 6 state required.')
        return c

    def in_region(self,c,ready=False):
        obs=np.asarray(c)@WEIGHTS.T; lo,hi=(159/160,161/160) if ready else (.9,1.1)
        return bool(np.all((obs[:,:2]>=lo-1e-14)&(obs[:,:2]<=hi+1e-14)) and np.all(obs[:,2]>=(.05 if ready else .0002)-1e-14))

    def rhs(self,t,v):
        c=v[:6*self.n].reshape(self.n,6); j=self.chemistry.fluxes(c)
        counters=np.column_stack((c@I,c[:,2],self.chemistry.service(c),j[:,0]+j[:,3]-j[:,5]))
        return np.r_[(self.chemistry.field(c)+self.exchange.field(c)).ravel(),counters.ravel(),self.exchange.gross_handling(c)]

    def evolve(self,initial,duration,method='Radau'):
        c=self.validate_state(initial)
        if not np.isfinite(duration) or duration<=0: raise ValueError('Positive duration required.')
        sol=solve_ivp(self.rhs,(0,duration),np.r_[c.ravel(),np.zeros(4*self.n+1)],
                      method=method,rtol=RTOL,atol=ATOL,dense_output=True)
        if not sol.success: raise RuntimeError(sol.message)
        time=np.linspace(0,duration,round(duration*50)+1)
        samples=sol.sol(time)[:6*self.n].T.reshape(-1,self.n,6)
        if samples.min()<-1e-9: raise RuntimeError('Negative numerical concentration; no clipping applied.')
        final=samples[-1]; totals=sol.y[6*self.n:-1,-1].reshape(self.n,4)
        collection=(sol.sol(duration)-sol.sol(duration-1))[6*self.n:-1].reshape(self.n,4) if duration>=1 else None
        material=np.stack([1+expm((self.exchange.D-np.eye(self.n))*t)@(c@WEIGHTS[:2].T-1) for t in time])
        residual=float(np.max(abs(samples@WEIGHTS[:2].T-material)))
        inventory=float(abs((final@I).sum()-(c@I).sum()+totals[:,0].sum()-totals[:,3].sum()))
        return dict(end=final,time=time,samples=samples,totals=totals,collection=collection,
                    handling=float(sol.y[-1,-1]),material_error=residual,inventory_error=inventory)


class HarvestMission:
    """Actual endpoint feeds the next pulse; all effluent and losses stay in the ledger."""
    def __init__(self,network,policy):self.network=network; self.policy=policy

    def run(self,initial,cycles,conditioning=True,method='Radau'):
        if not isinstance(cycles,int) or isinstance(cycles,bool) or cycles<0:raise ValueError('Nonnegative integer cycle count required.')
        c=self.network.validate_state(initial)
        if not self.network.chemistry.certified or not self.network.in_region(c,ready=not conditioning):
            raise ValueError('Mission certificate requires paper rates and admitted/ready initial state.')
        first=float((c@I).sum()); history=[]; traces=[]; all_totals=np.zeros((self.network.n,4))
        food=np.zeros(2); withdraw=loss=handling=0.; elapsed=0.; max_material=max_inventory=0.
        import copy
        for stage in range(-int(conditioning),cycles):
            pulse=self.policy.choose(c.copy(),copy.deepcopy(history)); after,removed,waste,added=pulse.apply(c)
            withdraw+=float((removed@I).sum()); loss+=float((waste@I).sum()); food+=added.sum(axis=0)
            duration=12 if stage<0 else 4
            part=self.network.evolve(after,duration,method); c=part['end']; food+=duration*self.network.n
            all_totals+=part['totals']; handling+=part['handling']
            max_material=max(max_material,part['material_error']); max_inventory=max(max_inventory,part['inventory_error'])
            traces.extend((elapsed+t,stage,i,*row) for t,cs in zip(part['time'],part['samples']) for i,row in enumerate(cs))
            elapsed+=duration
            if stage>=0:history.append(dict(cycle=stage+1,retained=list(pulse.retained),collection_I=part['collection'][:,0].tolist(),
                collection_X=part['collection'][:,1].tolist(),endpoint_Y=(c@Y).tolist()))
        net=float(all_totals[:,3].sum()); wash=float(all_totals[:,0].sum()); final=float((c@I).sum())
        return dict(history=history,trace=traces,end=c.tolist(),food=food.tolist(),gross_service=float(all_totals[:,2].sum()),
            gross_handling=handling,initial_I=first,final_I=final,all_effluent_I=wash,withdrawn_I=withdraw,lost_I=loss,
            net_synthesis=net,inventory_residual=abs(net-(final-first+wash+withdraw+loss)),
            material_error=max_material,flow_inventory_error=max_inventory,duration=elapsed)


class MissionCertificate:
    """Exact arithmetic on the imported deterministic theorem, not an ODE proof."""
    def __init__(self,n):
        if not isinstance(n,int) or isinstance(n,bool) or n<1:raise ValueError('Positive integer node count required.')
        self.n=n

    def bounds(self,m,conditioning=True):
        if not isinstance(m,int) or isinstance(m,bool) or m<0:raise ValueError('Nonnegative integer cycle count required.')
        n=self.n; start=Q(11,10) if conditioning else Q(161,160)
        return dict(cycles=m,time=12*int(conditioning)+4*m,collection_I=Q(n*m,28),collection_X=Q(n*m,160),
            each_food=Q(n*(2551*int(conditioning)+951*m),200),gross_service=Q(n*(27*int(conditioning)+9*m),50),
            net_synthesis= n*(Q(m+1,28)-start))

    def size(self,demand_I,demand_X,time,food_U,food_W,service):
        vals=list(map(Q,(demand_I,demand_X,time,food_U,food_W,service)))
        if min(vals)<0:raise ValueError('Nonnegative demands and budgets required.')
        di,dx,t,bu,bw,bg=vals;n=self.n
        available=min((t-12)/4,(200*bu/n-2551)/951,(200*bw/n-2551)/951,(50*bg/n-27)/9)
        if available<0:raise ValueError('Conditioning is not funded; zero routine cycles is not a feasible mission.')
        required=math.ceil(max(28*di/n,160*dx/n)); funded=math.floor(available)
        return dict(required_cycles=required,funded_cycles=funded,feasible=required<=funded,required_bounds=self.bounds(required))


class PhaseCertificate:
    """Lower-system weights convert total catalytic stock into free X locally.

    This Metzler lower matrix is not the reactor Jacobian. Nonnegative truncated
    exponential rows provide exact reusable weight certificates at rational lags.
    """
    M=sp.Matrix([[-48,20,0,38],[0,-43,20,0],[0,0,-41,2],[0,0,20,-24]])
    weights=sp.Matrix([[1,sp.Rational(9,8),sp.Rational(7,5),sp.Rational(9,5)]])

    def polynomial(self,degree=4):
        if not isinstance(degree,int) or degree<0:raise ValueError('Nonnegative degree required.')
        h=sp.Symbol('h',nonnegative=True); K=self.M+48*sp.eye(4); e=sp.Matrix([[1,0,0,0]])
        p=sum((h**k/sp.factorial(k)*e*K**k for k in range(degree+1)),sp.zeros(1,4))
        return h,p,sp.simplify(p*K-p.diff(h))

    def exact_floor(self):
        h,p,residual=self.polynomial(); ratios=[Q(v)/Q(w) for v,w in zip(p.subs(h,sp.Rational(1,28)),self.weights)]
        z=Q(6,7); upper=sum(z**k/math.factorial(k) for k in range(8))+z**8*Q(9,8)/math.factorial(8)
        floor=min(ratios)/20/upper**2
        assert floor>Q(1,160) and all(sp.Poly(v,h).is_zero or all(c>=0 for c in sp.Poly(v,h).all_coeffs()) for v in residual)
        return dict(polynomial=[str(v) for v in p],residual=[str(v) for v in residual],lag='1/28',
            ratios=list(map(str,ratios)),exp_6_over_7_upper=str(upper),certified_floor=str(floor),reported_floor='1/160')


def algebra_checks():
    """Fresh symbolic identities and rational comparisons; no stored verdicts loaded."""
    u,w,x,c1,c2,z,r,d=sp.symbols('u w x c1 c2 z r d'); c=sp.Matrix([u,w,x,c1,c2,z])
    eps,eta=map(sp.Rational,(EPS,ETA)); j=sp.Matrix([eps*(u*w-x/10),20*(x*u-c1),20*(c1*w-c2),20*c2-2*z,r*(z-x*x),d*(x-eta*u*w)])
    f=sp.Matrix([1,1,0,0,0,0])-c+sp.Matrix(STOICH)*j
    a,b,yy,ii=[sp.Matrix(list(map(lambda v:sp.Rational(str(v)),row))) for row in [A,B,[0,0,1,Q(9,8),Q(7,5),Q(9,5)],I]]
    assert sp.expand(a.dot(f)-1+a.dot(c))==sp.expand(b.dot(f)-1+b.dot(c))==0
    assert sp.expand(ii.dot(f)+ii.dot(c)-j[0]-j[3]+j[5])==0
    # Guard Y<=1/20, A,B>=.9: nonnegative slack products prove local growth.
    drift=sp.expand(yy.dot(f)); lower=(eps+d*eta)*u*w+(sp.Rational(1,9)-eps/10)*x+sp.Rational(51,280)*c1+c2/6+(r-19)*z/5
    aa,bb,ee,rho,sigma=sp.symbols('a b e rho sigma',nonnegative=True)
    gap=sp.expand(drift-sp.Rational(2,3)*yy.dot(c)-lower).subs({u:sp.Rational(73,90)+aa,w:sp.Rational(29,35)+bb,d:sp.Rational(1,25)-ee})
    gap=sp.expand(gap+r*x*x/5-sp.Rational(21,5)*(x/20-sigma)+rho*x*x/5)
    assert all(v>=0 for v in sp.Poly(gap,aa,bb,ee,rho,sigma,x,c1,c2,z).coeffs())
    # Shared thermochemistry includes the driven F/P pair; transport changes location only.
    full=sp.Matrix(STOICH).col_join(sp.Matrix([[0,0,0,0,0,-1],[0,0,0,0,0,1]]))
    composition=sp.Matrix([[1,0,1,2,2,2,0,0],[0,1,1,1,2,2,0,0],[0,0,0,0,0,0,1,1]])
    boltz=[sp.Integer(v) for v in [1,1,10,10,10,100]]+[sp.Rational(1,80000000000),sp.Integer(1)]
    assert composition*full==sp.zeros(3,6)
    ratios=[sp.prod(boltz[s]**full[s,jj] for s in range(8)) for jj in range(6)]
    assert ratios==[10,1,1,10,1,8000000000]
    # Exact boundary witness: only X moves; both nodes start at the SAME minimum Y.
    local=drift.subs({u:sp.Rational(151,160),w:sp.Rational(151,160),x:sp.Rational(1,20),c1:0,c2:0,z:0,r:21,d:sp.Rational(1,25)})
    selective=local-sp.Rational(1,10)
    assert selective==sp.Rational(-227999990907999,5120000000000000)<0<local
    taylor=lambda t,n:sum(Q(t)**k/math.factorial(k) for k in range(n+1))
    assert taylor(3,4)>=16 and taylor(Q(36,5),10)>=Q(50000,49) and taylor(Q(3,2),4)>=Q(200,49)
    lo=min(q*Q(49,50)*Q(9,10)+1-q-Q(1,200) for q in [Q(1,4),Q(3,4)])
    hi=max(q*Q(11,10)+1-q+Q(1,200) for q in [Q(1,4),Q(3,4)])
    assert (lo,hi)==(Q(1813,2000),Q(27,25))
    return dict(material_and_inventory='exact symbolic identities',guard_slack=str(gap),
        thermochemical_rate_ratios=list(map(str,ratios)),pulse_material_interval=[str(lo),str(hi)],
        selective=dict(local_Y_drift=str(local),X_only_exchange='-1/10',coupled_Y_drift=str(selective),
                       scope='Failure of the common-minimum argument; not proof of extinction.'),phase=PhaseCertificate().exact_floor())


def pilot(name,scale=1.,method='Radau'):
    network=ReactorNetwork(LocalChemistry(),CommonExchange.four_nodes(name,scale))
    c=network.validate_state(INITIAL_STATE); post=Pulse(BASE_RETENTION).apply(c)[0]
    p=network.evolve(post,4,method)
    result=dict(topology=name,scale=scale,collection_I_min=float(p['collection'][:,0].min()),collection_X_min=float(p['collection'][:,1].min()),
        final_Y_min=float((p['end']@Y).min()),final_material_deviation=float(abs(p['end']@WEIGHTS[:2].T-1).max()),
        handling=p['handling'],handling_bound=network.exchange.handling_bound(4),material_error=p['material_error'],inventory_error=p['inventory_error'])
    return result,p


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)
    algebra=algebra_checks(); net=ReactorNetwork(LocalChemistry(),CommonExchange.four_nodes(TOPOLOGY,EXCHANGE_SCALE))
    mission=HarvestMission(net,HistoryPolicy()).run(INITIAL_STATE,ROUTINE_CYCLES)
    cert=MissionCertificate(net.n); bounds=cert.bounds(ROUTINE_CYCLES)
    pilots=[]; curves=[]
    for name in ('uncoupled','path','ring','star'):
        result,part=pilot(name);pilots.append(result)
        curves.extend((name,t,i,*row) for t,cs in zip(part['time'],part['samples']) for i,row in enumerate(cs))
    _,independent=pilot('path',method='BDF');_,reference=pilot('path')
    disagreement=float(np.max(abs(independent['samples']-reference['samples'])))
    assert disagreement<1e-7 and mission['inventory_residual']<1e-7 and mission['material_error']<1e-7
    for h in mission['history']:
        assert min(h['collection_I'])>=1/28 and min(h['collection_X'])>=1/160
    assert max(mission['food'])<=float(bounds['each_food']) and mission['gross_service']<=float(bounds['gross_service'])
    assert mission['net_synthesis']>=float(bounds['net_synthesis'])
    traces=mission.pop('trace');scale=CONCENTRATION_MOLAR*NODE_VOLUME_LITERS*1e6
    physical=dict(amount_per_normalized_unit_umol=scale,seconds_per_time_unit=TIME_SECONDS,
        node_volume_mL=1000*NODE_VOLUME_LITERS,eleven_cycles={k:float(v)*scale for k,v in cert.bounds(11).items() if k not in ('cycles','time')},
        eleven_cycles_minutes=56*TIME_SECONDS/60,scope='Illustrative unit conversion; no measured chemistry, pumping energy, or finite reservoir depletion model.')
    result=dict(manuscript_sha256=MANUSCRIPT_SHA256,algebra=algebra,pilots=pilots,independent_solver_max_difference=disagreement,
        mission=mission,certificate=bounds,mission_sizing=cert.size(DEMAND_I,DEMAND_X,AVAILABLE_TIME,FOOD_U_BUDGET,FOOD_W_BUDGET,SERVICE_BUDGET),
        physical=physical,exchange_labels=net.exchange.directed_labels(),scope='Deterministic synchronized equal-volume network. The global theorem is imported; exact algebra and numerical diagnostics are freshly run. No Lean replay or stochastic claim.')
    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('mission.csv',['time','cycle_index','node',*SPECIES],traces);table('pilots.csv',['topology','time','node',*SPECIES],curves)
    sweep=[]
    for strength in (0.,.05,.2,1.,5.,20.):
        p,_=pilot('path',strength);sweep.append([strength,p['collection_I_min'],p['collection_X_min'],p['handling'],p['handling_bound']])
    table('exchange_sweep.csv',['scale','min_collection_I','min_collection_X','handling','handling_bound'],sweep)
    lines=[f'Exact material, inventory, guarded-growth, thermochemical and backward-weight identities passed.',
        f'Each node, each routine cycle: I output >= 1/28 and free X >= 1/160; these outputs overlap.',
        f'{ROUTINE_CYCLES} cycles plus conditioning: actual net synthesis {mission["net_synthesis"]:.6g}; theorem lower bound {bounds["net_synthesis"]}.',
        f'Global inventory residual {mission["inventory_residual"]:.3g}; independent solver difference {disagreement:.3g}.',
        f'Worked mission requires {result["mission_sizing"]["required_cycles"]} cycles; budgets fund {result["mission_sizing"]["funded_cycles"]}.',
        'Selective-X witness has negative Y drift at a shared minimum; no extinction conclusion follows.']
    (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,3,figsize=(12,3.8),layout='constrained')
    for name in ('uncoupled','path','ring','star'):
        data=np.array([r[1:] for r in curves if r[0]==name]);states=data[:,2:].reshape(-1,4,6);t=data[::4,0]
        axs[0].plot(t,(states@Y).min(axis=1),label=name.title());axs[1].plot(t,states[:,:,2].min(axis=1),label=name.title())
    axs[0].axhline(.05,color='k',ls=':',label='Ready stock 1/20');axs[1].hlines(1/160,3,4,color='k',linestyle=':',label='Collection floor 1/160')
    axs[0].set(xlabel='Time after routine pulse',ylabel='Minimum node Y',title='Minimum catalytic stock after withdrawal');axs[1].set(xlabel='Time after routine pulse',ylabel='Minimum node free X',title='Minimum free-template concentration during\nrecovery');axs[0].legend(fontsize=7);axs[1].legend(fontsize=7)
    xs=np.arange(4);axs[2].bar(xs-.18,[p['handling'] for p in pilots],.36,label='Integrated movement');axs[2].bar(xs+.18,[p['handling_bound'] for p in pilots],.36,label='Theorem allowance');axs[2].set_xticks(xs,['None','Path','Ring','Star']);axs[2].set(ylabel='Gross inter-node movement',title='Transport amounts by network structure');axs[2].legend(fontsize=7)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'networks.png',dpi=180);fig.savefig(out/'networks.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,3.8),layout='constrained');h=mission['history'];cy=[v['cycle'] for v in h]
    for i in range(net.n):axs[0].plot(cy,[v['collection_X'][i] for v in h],'.-',label=f'Node {i+1}')
    axs[0].axhline(1/160,color='k',ls=':',label='Guaranteed 1/160');axs[0].set(xlabel='Routine cycle',ylabel='Credited free-X effluent',title='Free-template collection over successive\ncycles');axs[0].legend(fontsize=7)
    m=np.arange(0,max(33,ROUTINE_CYCLES+1));axs[1].plot(m,[float(cert.bounds(int(k))['net_synthesis']) for k in m],label='Admitted start + conditioning');axs[1].plot(m,[float(cert.bounds(int(k),False)['net_synthesis']) for k in m],label='Ready start');axs[1].axhline(0,color='k',lw=.7);axs[1].set(xlabel='Routine cycles',ylabel='Net synthesis lower bound',title='Net synthesis bounds after initial-stock\nsubtraction');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'mission.png',dpi=180);fig.savefig(out/'mission.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()