"""Chemical completion, marked reactor, cycle thermodynamics and certified sizing."""
# EDITABLE INPUTS: schematic dimensionless chemistry and illustrative physical units.
COPY_SCALE = 100_000_000
RELEASE_SPEED = '20'
CLEAVAGE_SPEED = '3/100'
FUEL_ACTIVITY = '1'
WASTE_ACTIVITY = '1'
POST_STARTUP_DURATION = '100'
TARGET_FAILURE = '1/100'
REFERENCE_MOLAR = '1/1000'
TIME_UNIT_SECONDS = '60'
TEMPERATURE_K = '298.15'
BATH_TOLERANCE = '1/100'
SSA_EVENT_BUDGET = 20000
SSA_SEED = 41092026

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

mp.dps=70;iv.dps=70
MANUSCRIPT_SHA256='2b3d026fac78d73fdb8733abd0fb3b79a45b4a8f50f454159589c97a298eae6c'
SPECIES=('U','W','X','C1','C2','Z','F','P')
COMPOSITION=((1,0,0),(0,1,0),(1,1,0),(2,1,0),(2,2,0),(2,2,0),(0,0,1),(0,0,1))
BOLTZMANN=(Q(1),Q(1),Q(10),Q(10),Q(10),Q(100),Q(1,80000000000),Q(1))  # exp(-g/RT)
EPSILON=Q(1,500000000);ETA=Q(1,8000000000)


def number(ctx,value):
    value=Q(str(value));return ctx.mpf(value.numerator)/value.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 ChemicalPair:
    name:str
    left:tuple
    right:tuple
    forward:Q
    reverse:Q
    labels:tuple
    @property
    def jump(self):return tuple(self.right.count(i)-self.left.count(i) for i in range(8))
    def thermodynamic_ratio(self):return math.prod(b**v for b,v in zip(BOLTZMANN,self.jump))
    def projected(self,activities):
        result=[]
        for inputs,outputs,k,label in [(self.left,self.right,self.forward,self.labels[0]),(self.right,self.left,self.reverse,self.labels[1])]:
            marks=[0]*8
            if label==18:marks[3]=1
            if label==19:marks[4]=1
            effective=k*math.prod(activities[i-6] for i in inputs if i>=6)
            result.append(reference.Channel(label,self.name+('_forward' if label==self.labels[0] else '_reverse'),tuple(i for i in inputs if i<6),tuple(i for i in outputs if i<6),effective,tuple(marks)))
        return result


class CompletedExporter:
    def __init__(self,release=RELEASE_SPEED,delta=CLEAVAGE_SPEED,fuel=FUEL_ACTIVITY,waste=WASTE_ACTIVITY,enabled=True):
        self.release,self.delta,self.fuel,self.waste=map(Q,(release,delta,fuel,waste));self.enabled=enabled
        if min(self.release,self.delta,self.fuel,self.waste)<=0:raise ValueError('Positive paired speeds and reservoir activities required.')
        self.p=reference.Parameters(K=10,release=self.release,delta=self.delta)
        self.pairs=(ChemicalPair('basal',(0,1),(2,),EPSILON,EPSILON/10,(0,1)),ChemicalPair('bind_U',(2,0),(3,),Q(20),Q(20),(2,3)),
            ChemicalPair('bind_W',(3,1),(4,),Q(20),Q(20),(4,5)),ChemicalPair('ligation',(4,),(5,),Q(20),Q(2),(6,7)),
            ChemicalPair('release',(5,),(2,2),self.release,self.release,(8,9)),ChemicalPair('driven',(2,6),(0,1,7),self.delta,self.delta*ETA,(18,19)))
        channels=[]
        for j,pair in enumerate(self.pairs):
            if enabled or j!=3:channels.extend(pair.projected((self.fuel,self.waste)))
        for i in (0,1):
            marks=[0]*8;marks[1+i]=1;channels.append(reference.Channel(10+i,'feed_'+SPECIES[i],(),(i,),Q(1),tuple(marks)))
        for i in range(6):
            marks=(reference.S_WEIGHTS[i],0,0,0,0,reference.A_WEIGHTS[i],reference.B_WEIGHTS[i],reference.M_WEIGHTS[i])
            channels.append(reference.Channel(12+i,'wash_'+SPECIES[i],(i,),(),Q(1),marks))
        self.channels=tuple(sorted(channels,key=lambda c:c.label));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):
            used={}
            for k,i in enumerate(c.inputs):self.inputs[j,k]=i;self.offsets[j,k]=used.get(i,0);used[i]=used.get(i,0)+1
    def certificate_scope(self):return 19<=self.release<=21 and Q(1,50)<=self.delta<=Q(1,25) and self.fuel==self.waste==1
    def rates(self,state,V=None):return reference.DrivenReactor.rates(self,state,V)
    def deterministic(self,H,method='Radau'):
        H=float(Q(str(H)))
        if not 1<=H<=10000:raise ValueError('Density illustration requires 1<=H<=10000.')
        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(reference.Y_WEIGHTS,float)@state[:6]-1/2500
        entry.direction=1
        solution=solve_ivp(rhs,(0,500+H),[1,1,0,0,0,0]+[0]*8,method=method,dense_output=True,events=entry,rtol=1e-11,atol=1e-18,max_step=.5)
        if not solution.success:raise ArithmeticError(solution.message)
        return solution
    def stochastic(self,policy,limit=SSA_EVENT_BUDGET,seed=SSA_SEED):return reference.DrivenReactor.stochastic(self,policy,limit,seed)
    def certificates(self):
        S=sp.Matrix([p.jump for p in self.pairs]).T;E=sp.Matrix(COMPOSITION).T;cp=sp.Matrix([-1,1,1,1,1,0]);cd=sp.Matrix([1,0,0,0,0,1])
        if E*S!=sp.zeros(3,6) or S[:6,:].rank()!=4 or S.rank()!=5 or S*cp!=sp.zeros(8,1) or list(S*cd)!=[0,0,0,0,0,0,-1,1]:raise ArithmeticError('Material balance or cycle classification failed.')
        ratios=[pair.forward/pair.reverse for pair in self.pairs]
        if any(r!=p.thermodynamic_ratio() for r,p in zip(ratios,self.pairs)):raise ArithmeticError('No common standard potentials for the supplied coefficients.')
        if math.prod(r**int(c) for r,c in zip(ratios,cp))!=1 or math.prod(r**int(c) for r,c in zip(ratios,cd))!=80000000000:raise ArithmeticError('Cycle affinity failed.')
        expected=reference.DrivenReactor(self.p,self.enabled)
        projection=self.fuel==self.waste==1 and all((c.inputs,c.outputs,c.coefficient,c.marks)==(expected.by_label[c.label].inputs,expected.by_label[c.label].outputs,expected.by_label[c.label].coefficient,expected.by_label[c.label].marks) for c in self.channels)
        if self.fuel==self.waste==1 and not projection:raise ArithmeticError('Label-wise source identity failed.')
        return dict(composition=COMPOSITION,full_stoichiometry=S.tolist(),internal_rank=4,full_rank=5,passive_cycle=list(cp),drive_cycle=list(cd),rate_ratios=list(map(str,ratios)),boltzmann_weights=list(map(str,BOLTZMANN)),drive_affinity='log(80000000000)',unit_activity_marked_source_identity=projection)
    def neighboring_ratio(self,pair_index,surplus,V):
        if type(V) is not int or V<=0 or len(surplus)!=6 or any(type(n) is not int or n<0 for n in surplus):raise ValueError('Natural surplus counts and positive integer copy scale required.')
        pair=self.pairs[pair_index];forward,reverse=pair.projected((self.fuel,self.waste))
        before=tuple(n+pair.left.count(i) for i,n in enumerate(surplus));after=tuple(n+pair.right.count(i) for i,n in enumerate(surplus))
        ratio=forward.count_rate(before,V)/reverse.count_rate(after,V)
        # exp[-(G(after)-G(before))], evaluated exactly using rational Boltzmann weights.
        state_ratio=math.prod((BOLTZMANN[i]*V)**(after[i]-before[i])*Q(math.factorial(before[i]),math.factorial(after[i])) for i in range(6))
        reservoir=math.prod(BOLTZMANN[i]**pair.jump[i]*activity**(-pair.jump[i]) for i,activity in ((6,self.fuel),(7,self.waste)))
        if ratio!=state_ratio*reservoir:raise ArithmeticError('Neighboring-count detailed balance failed.')
        return dict(pair=pair.name,before=before,after=after,propensity_ratio=str(ratio),count_potential_factor=str(state_ratio),reservoir_factor=str(reservoir))


class OperatingCertificate:
    def __init__(self,V,H,model):
        self.V=V;self.H=Q(str(H));self.model=model
        if type(V) is not int or V<=0 or self.H<0 or not model.certificate_scope() or not model.enabled:raise ValueError('Enabled certificate requires integer V>0, H>=0, paired rectangle and unit maintained activities.')
    def logs(self,ctx=mp):
        V=ctx.mpf(self.V);H=number(ctx,self.H);T=500+H;m=math.floor(self.H)
        result={'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,'residence_initial':-V/62500,
            'residence_leak':ctx.log(3000)+ctx.log(V)+ctx.log(T)+number(ctx,'18/125')-V/15625,
            'two_foods':ctx.log(2)-(2*ctx.log(2)-1)*V*T,'gross_service':-V*T/200}
        if m:result['windows']=ctx.log(m)+number(ctx,'1/8')-ctx.mpf(13)*V/1080000
        return result
    def passes(self,failure):
        failure=Q(str(failure))
        if not 0<failure<1:raise ValueError('Failure target in (0,1) required.')
        return bool(logsum(iv,self.logs(iv).values()).b<=iv.log(number(iv,failure)).a)
    def evaluate(self):
        logs=self.logs();L=logsum(mp,logs.values());interval=logsum(iv,self.logs(iv).values());success=max(mp.mpf(0),1-mp.exp(L))
        # Downward six-decimal claim is accepted only after an outward interval comparison.
        digits=int(mp.floor(success*1000000));claim=Q(digits,1000000)
        while claim>0 and not self.passes(1-claim):digits-=1;claim=Q(digits,1000000)
        off=None
        if self.H.denominator==1 and self.H>=1:
            parts=[-mp.mpf(self.V)*int(self.H)/25000,logs['resource_initial'],logs['resource_leak']]
            exactlogs=self.logs(iv);parts_iv=[-iv.mpf(self.V)*int(self.H)/25000,exactlogs['resource_initial'],exactlogs['resource_leak']]
            logoff=logsum(mp,parts);Loff=logsum(iv,parts_iv);exponent=max(0,int(mp.floor(-logoff/mp.log(10))))
            while exponent and not bool(Loff.b<(-exponent*iv.log(10)).a):exponent-=1
            off=dict(log10_ordinary=mp.nstr(logoff/mp.log(10),25),certified_less_than_10_power_minus=exponent if bool(Loff.b<iv.mpf(0).a) else None,log_outward_interval=str(Loff))
        return dict(V=self.V,H=str(self.H),success_lower_ordinary=mp.nstr(success,30),success_lower_certified_six_decimals=f'{digits/1000000:.6f}',failure_log_interval=str(interval),failure_log10_ordinary=mp.nstr(L/mp.log(10),25),component_log10={k:mp.nstr(v/mp.log(10),25) for k,v in logs.items()},disabled=off)
    @staticmethod
    def necessary_scale(failure,delta=Q(1,25)):
        failure=Q(str(failure));delta=Q(str(delta))
        if not 0<failure<1 or not Q(1,50)<=delta<=Q(1,25):raise ValueError('Target and cleavage speed outside initiation theorem scope.')
        value=iv.log(1/number(iv,failure))/(500*number(iv,EPSILON+delta*ETA));lo=mp.mpf(value.a);hi=mp.mpf(value.b)
        if mp.ceil(lo)!=mp.ceil(hi):raise ArithmeticError('Necessary integer rounding unresolved.')
        return int(mp.ceil(lo))
    @classmethod
    def sufficient_scales(cls,H,failure,model):
        H=Q(str(H));failure=Q(str(failure))
        if H<1 or not 0<failure<1:raise ValueError('Sizing requires H>=1 and target failure in (0,1).')
        required=max(iv.mpf(100000000).b,(1000000*iv.log(number(iv,H))).b,(50000000*iv.log(2/number(iv,failure))).b)
        simple=int(mp.ceil(mp.mpf(required)));lower=100000000;upper=max(simple,lower)
        if not cls(upper,H,model).passes(failure):raise ArithmeticError('Simple sufficient size failed full-budget check.')
        # Each error component decreases with V on V>=1e8; search only this domain.
        while lower<upper:
            midpoint=(lower+upper)//2
            if cls(midpoint,H,model).passes(failure):upper=midpoint
            else:lower=midpoint+1
        if not cls(lower,H,model).passes(failure):raise ArithmeticError('Final sufficient size unresolved.')
        return dict(simple_sufficient_V=simple,full_budget_interval_checked_V=lower,scope='Smallest accepted integer in the V>=1e8 monotone search domain; not the smallest physically successful reactor.')


class PhysicalAccounting:
    def __init__(self,V,H,concentration=REFERENCE_MOLAR,time_unit=TIME_UNIT_SECONDS,temperature=TEMPERATURE_K):
        self.V=V;self.H=Q(str(H));self.c=Q(str(concentration));self.tau=Q(str(time_unit));self.temperature=Q(str(temperature))
        if V<=0 or self.H<1 or min(self.c,self.tau,self.temperature)<=0:raise ValueError('Positive scale, H>=1 and physical units required.')
    def evaluate(self,model,rho=BATH_TOLERANCE):
        rho=Q(str(rho))
        if not 0<rho<1 or not model.certificate_scope() or not model.enabled:raise ValueError('Conditional accounting requires the enabled unit-activity certificate and bath tolerance in (0,1).')
        V=self.V;T=500+self.H;m=math.floor(self.H);gross=V*T/8;stock=math.ceil(gross/rho);NA=mp.mpf('6.02214076e23');R=mp.mpf('8.31446261815324');force=R*number(mp,self.temperature)*mp.log(80000000000)
        amount=lambda n:mp.nstr(number(mp,n)/NA,25)
        rates=[]
        for pair in model.pairs:
            rates.append(dict(pair=pair.name,forward_order=len(pair.left),reverse_order=len(pair.right),forward_physical=str(pair.forward/(self.tau*self.c**(len(pair.left)-1))),reverse_physical=str(pair.reverse/(self.tau*self.c**(len(pair.right)-1)))))
        return dict(volume_pL=mp.nstr(V/(NA*number(mp,self.c))*10**12,25),startup_hours=str(500*self.tau/3600),total_hours=str(T*self.tau/3600),
            template_equivalents_per_window=str(Q(V,20000)),completed_template_equivalents=str(Q(V*m,20000)),completed_export_moles=amount(Q(V*m,20000)),
            each_food_allowance=str(2*V*T),initial_each_food=V,gross_driven_allowance=str(gross),gross_driven_moles=amount(gross),service_per_equivalent_upper=str(2500*T/m),
            reservoir_force_kJ_per_mole=mp.nstr(force/1000,25),chemical_work_bound_nJ=mp.nstr(force*number(mp,gross)/NA*10**9,25),full_complex_dimensional_coefficients=rates,
            conditional_bath_each_stock=stock,bath_volume_nL=mp.nstr(number(mp,stock)/(NA*number(mp,self.c))*10**9,25),
            activity_ratio_corridor=[str(ETA*(1-rho)/(1+rho)),str(ETA*(1+rho)/(1-rho))],force_variation_bound_RT=mp.nstr(mp.log(number(mp,(1+rho)/(1-rho))),25),
            scope='Allowances and output guarantees are conditional on the maintained-model event. Bath sizing is conditional accounting only; no operating probability transfers to a moving bath. Export includes complexes/duplexes, not purified free X. Work excludes feeds, pumping, maintenance and purification.')


class FiniteBath:
    """Explicit exchange bookkeeping; activities change the generator."""
    def __init__(self,each_stock):
        if type(each_stock) is not int or each_stock<=0:raise ValueError('Positive integer bath stock required.')
        self.initial=each_stock;self.fuel=each_stock;self.waste=each_stock;self.gross=0
    def exchange(self,label):
        if label not in (18,19):return
        if label==18:
            if self.fuel==0:raise ValueError('Fuel exhausted.')
            self.fuel-=1;self.waste+=1
        else:
            if self.waste==0:raise ValueError('Waste exhausted.')
            self.fuel+=1;self.waste-=1
        self.gross+=1
    @property
    def activities(self):return Q(self.fuel,self.initial),Q(self.waste,self.initial)


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)
    model=CompletedExporter();disabled=CompletedExporter(enabled=False);chemistry=model.certificates();disabled.certificates()
    neighbors=[model.neighboring_ratio(i,(3,5,2,4,1,7),31) for i in range(6)]
    certificate=OperatingCertificate(COPY_SCALE,POST_STARTUP_DURATION,model);bounds=certificate.evaluate();sizes=certificate.sufficient_scales(POST_STARTUP_DURATION,TARGET_FAILURE,model)
    necessary=certificate.necessary_scale(TARGET_FAILURE);physical=PhysicalAccounting(COPY_SCALE,POST_STARTUP_DURATION).evaluate(model)
    H=Q(POST_STARTUP_DURATION)
    if not 1<=H<=10000:raise ValueError('Default trajectory illustration requires 1<=H<=10000; the bound interface supports H>=0 independently.')
    sol=model.deterministic(H);independent=model.deterministic(H,'BDF');off=disabled.deterministic(H)
    times=np.unique(np.r_[np.linspace(0,20,201),np.linspace(20,float(500+H),581)]);states=sol.sol(times).T;other=off.sol(times).T
    rows=[OperatingCertificate(V,H,model).evaluate() for V in (100000000,200000000,300000000)]
    scale_sweep=[]
    for V in np.linspace(100000000,300000000,51).astype(np.int64):
        cert=OperatingCertificate(int(V),H,model);entry=cert.evaluate();floor=-500*number(mp,EPSILON+Q(1,25)*ETA)*int(V)/mp.log(10)
        scale_sweep.append((int(V),float(floor),float(entry['failure_log10_ordinary'])))
    result=dict(chemistry=chemistry,neighboring_count_checks=neighbors,operating_bound=bounds,worked_scales=rows,target_failure=TARGET_FAILURE,target_certified_at_selected_scale=certificate.passes(TARGET_FAILURE),
        necessary_uniform_copy_scale=necessary,sufficient_sizes=sizes,physical_accounting=physical,deterministic_solver_difference=float(np.max(abs(states-independent.sol(times).T))),
        preparation='Food only: (V,V,0,0,0,0). Density trajectories illustrate dynamics, not finite-copy probability.',
        manuscript_arithmetic_note='The sentence after Corollary 6.2 says the full budget certifies 1e8 for 99% success; direct evaluation and Table 5 give only 0.923849 there. This package uses the evaluated bound and computes the target-specific sufficient scale.',
        scope='Chemical/source identities are freshly checked. Probability bounds are inherited symbolic theorems evaluated with outward intervals. Lean and infinite-state process identification are not rerun.')
    if args.simulate:result['limited_stochastic_path']=model.stochastic(reference.ObservationPolicy(COPY_SCALE,H))
    def write_json(name,data):(out/name).write_text(json.dumps(data,indent=2,default=str)+'\n')
    def table(name,header,data):
        with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(data)
    table('trajectory.csv',['time',*SPECIES[:6],*[name+'_per_V' for name in reference.LEDGER]],np.column_stack([times,states]))
    table('scale_bounds.csv',['V','necessary_failure_floor_log10','sufficient_failure_ceiling_log10'],scale_sweep)
    table('channels.csv',['label','name','internal_inputs','internal_outputs','effective_coefficient','export_mark','feed_U','feed_W','drive_F','drive_P'],[(c.label,c.name,','.join(SPECIES[i] for i in c.inputs),','.join(SPECIES[i] for i in c.outputs),str(c.coefficient),*c.marks[:5]) for c in model.channels])
    table('window_exports.csv',['window','enabled_marks_per_V','disabled_marks_per_V','guaranteed_marks_per_V'],[(i,float(sol.sol(501+i)[6]-sol.sol(500+i)[6]),float(off.sol(501+i)[6]-off.sol(500+i)[6]),1/5000) for i in range(math.floor(H))])
    write_json('results.json',result)
    lines=[f'Balanced six-pair chemistry; internal rank 4, full rank 5; exact label/mark projection verified.',
        f'Joint operating success >= {bounds["success_lower_certified_six_decimals"]}; disabled bound < 10^-{bounds["disabled"]["certified_less_than_10_power_minus"]} for the configured integer horizon.' if bounds['disabled'] and bounds['disabled']['certified_less_than_10_power_minus'] is not None else f'Joint operating success >= {bounds["success_lower_certified_six_decimals"]}; no informative disabled power bound reported.',
        f'For failure {TARGET_FAILURE}: necessary uniform V >= {necessary}; sufficient sizes {sizes}.',f'Physical volume {physical["volume_pL"]} pL; startup {physical["startup_hours"]} hours; specified reservoir work allowance {physical["chemical_work_bound_nJ"]} nJ.',
        'The 1e8 example does not certify 99% success. Mixed exported species are not purified product. Moving-bath probabilities need a new theorem.']
    (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(scale_sweep)
    axs[0].plot(data[:,0]/1e6,data[:,1],label='Necessary failure floor');axs[0].plot(data[:,0]/1e6,data[:,2],label='Sufficient failure ceiling');axs[0].set(xlabel='Copy scale V (millions)',ylabel='log10 failure probability',title='Bounds on the same joint operating event');axs[0].legend(fontsize=8)
    Y=np.array(reference.Y_WEIGHTS,float);axs[1].semilogy(times,np.maximum(states[:,:6]@Y,1e-14),label='Enabled density model');axs[1].semilogy(times,np.maximum(other[:,:6]@Y,1e-14),label='Ligation pair deleted');axs[1].axhline(1/2500,color='gray',ls='--',label='Entry threshold');axs[1].set(xlabel='Dimensionless time',ylabel='Weighted catalyst concentration Y/V',title='Food-only density trajectories (not probabilities)');axs[1].legend(fontsize=7)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'operation.png',dpi=180);fig.savefig(out/'operation.svg');plt.close(fig)
    fig,ax=plt.subplots(figsize=(8,4.2),layout='constrained');labels=['U feed allowance (upper)','W feed allowance (upper)','Gross drive allowance (upper)','Each initial food stock','Template equivalents (lower)'];values=[float(2*COPY_SCALE*(500+H))]*2+[float(COPY_SCALE*(500+H)/8),COPY_SCALE,float(COPY_SCALE*math.floor(H)/20000)]
    ax.barh(labels,values,color=['#417d8c']*3+['#777777','#bd5a24']);ax.set_xscale('log');ax.set(xlabel='Whole-run counts conditional on the operating event',title='Supply allowances include the startup period');ax.invert_yaxis();ax.grid(axis='x',alpha=.2)
    fig.savefig(out/'accounting.png',dpi=180);fig.savefig(out/'accounting.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    write_json('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,source_sha256=digest(Path(__file__)),reference_source_sha256=digest(Path(__file__).with_name('reference_source.py')),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()