"""A reusable mass-action assembly and exact complex-balanced realization certificates."""
# EDITABLE INPUTS: nondimensional mathematical examples, not fitted laboratory rates.
PARAMETERS = ('7/16','7/16','1','12/25','1/2','1')  # a,b,u,v,e,d
EXACT_REFERENCE = ('1/2','1/2','5/8','13/48')  # set None to solve numerically after changing parameters
INITIAL_STATES = ((.1,1.2,.2,.8),(1.4,.08,1.1,.04))
TIME_HORIZON = 35.
TIME_SAMPLES = 351
SLICE_ETAS = ('9/2','5','11/2','8')
ATTACHMENTS = (('W1',(1,0,1,0),'2','3','1/10'),('W2',(0,1,0,0),'1','2','2'))
GRID_E = ('1/20','4',81)
GRID_D = ('1/20','5/2',61)
ROOT_BISECTIONS = 100

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 scipy.integrate import solve_ivp

MANUSCRIPT_SHA256='ef9c1c93dfe470c759488cd4ad2a5352e81cd75b04412a78501fd68b13036770'
SPECIES=('A','B','z','H')
COMPLEX_NAMES=('0','A','B+z','z','H','2z','B','2A')
COMPLEXES=((0,0,0,0),(1,0,0,0),(0,1,1,0),(0,0,1,0),(0,0,0,1),(0,0,2,0),(0,1,0,0),(2,0,0,0))


@dataclass(frozen=True)
class Reaction:
    source:tuple
    target:tuple
    rate:Q


class MassActionNetwork:
    """Finite reactions, polynomial coefficients and positive-domain integration."""
    def __init__(self,species,reactions):
        self.species=tuple(species);self.reactions=tuple(reactions);n=len(species)
        if len(set(species))!=n:raise ValueError('Distinct species names required.')
        for r in self.reactions:
            if len(r.source)!=n or len(r.target)!=n or any(type(v) is not int or v<0 for v in r.source+r.target) or r.rate<0:raise ValueError('Natural exponents and nonnegative rates required.')
        self.Y=np.array([r.source for r in self.reactions],int);self.S=np.array([np.array(r.target)-r.source for r in self.reactions],float).T
        self.k=np.array([float(r.rate) for r in self.reactions]);self.n=n
    def coefficients(self):
        rows={}
        for r in self.reactions:
            row=rows.setdefault(r.source,[Q(0)]*self.n)
            for i in range(self.n):row[i]+=r.rate*(r.target[i]-r.source[i])
        return {y:tuple(c) for y,c in rows.items() if any(c)}
    def exact_field(self,state):
        state=tuple(map(Q,state))
        return tuple(sum(math.prod(x**y for x,y in zip(state,source))*c[i] for source,c in self.coefficients().items()) for i in range(self.n))
    def field(self,state):
        x=np.asarray(state,float)
        return self.S@(self.k*np.prod(x**self.Y,axis=1))
    def jacobian(self,state):
        x=np.asarray(state,float)
        if np.any(x<=0):raise ValueError('Positive state required.')
        rates=self.k*np.prod(x**self.Y,axis=1)
        return self.S@(rates[:,None]*self.Y/x)
    def balance(self,state):
        state=tuple(map(Q,state));imbalance={}
        for r in self.reactions:
            flux=r.rate*math.prod(x**y for x,y in zip(state,r.source))
            imbalance[r.source]=imbalance.get(r.source,Q(0))+flux
            imbalance[r.target]=imbalance.get(r.target,Q(0))-flux
        return imbalance
    def integrate(self,initial,times,method='Radau'):
        initial=np.asarray(initial,float);times=np.asarray(times,float)
        if initial.shape!=(self.n,) or np.any(initial<=0) or times[0]!=0 or np.any(np.diff(times)<=0):raise ValueError('Positive initial state and increasing times from zero required.')
        def fun(t,y):
            x=np.exp(y);return self.field(x)/x
        def jac(t,y):
            x=np.exp(y);return self.jacobian(x)*x[None,:]/x[:,None]-np.diag(self.field(x)/x)
        sol=solve_ivp(fun,(0,times[-1]),np.log(initial),method=method,jac=jac,t_eval=times,rtol=2e-10,atol=2e-12)
        values=np.exp(sol.y.T)
        if not sol.success or not np.all(np.isfinite(values)) or np.any(values<=0):raise ArithmeticError('Positive-domain integration failed.')
        return values


@dataclass(frozen=True)
class Parameters:
    a:Q;b:Q;u:Q;v:Q;e:Q;d:Q
    def __post_init__(self):
        for name in ('a','b','u','v','e','d'):
            object.__setattr__(self,name,Q(getattr(self,name)))
            if getattr(self,name)<=0:raise ValueError('All six parameters must be positive.')
    @classmethod
    def slice(cls,eta):
        eta=Q(eta)
        return cls((19-2*eta)/25,(eta-4)/25,1,5,eta,3)
    @classmethod
    def from_stationary_design(cls,A,B,z,u,v,e):
        """Choose rational concentrations and three rates; solve the other balances.

        Raises if this design needs a nonpositive feed or output rate.
        """
        A,B,z,u,v,e=map(Q,(A,B,z,u,v,e))
        if min(A,B,z,u,v,e)<=0:raise ValueError('Positive design inputs required.')
        K=A-B*z;J=e*(B-A*A);H=(u*z+2*v*z*z-K)/3
        if H<=0:raise ValueError('Design gives nonpositive H.')
        p=cls(2*A-B*z-2*J,(1+z)*B-A+J,u,v,e,(u*z+v*z*z)/H-2)
        return p,(A,B,z,H)
    def literal(self):
        edges=((1,2,1),(2,1,1),(3,4,self.u),(4,3,1),(5,4,self.v),(4,5,1),
               (0,1,self.a),(1,0,1),(0,6,self.b),(6,0,1),(6,7,self.e),(7,6,self.e),(4,0,self.d))
        return MassActionNetwork(SPECIES,[Reaction(COMPLEXES[i],COMPLEXES[j],Q(k)) for i,j,k in edges])
    def reduced(self,z):
        z=Q(z)
        if z<=-2:raise ValueError('Reduced functions require z > -2.')
        B=(self.a+2*self.b)/(z+2);K=(self.v*(1+2*self.d)*z*z-self.u*(1-self.d)*z)/(2+self.d)
        A=z*B+K;H=(self.u*z+self.v*z*z)/(2+self.d)
        E=self.b-(1+self.e)*B+K+self.e*A*A
        return (A,B,z,H),K,E
    def cut(self):return max(Q(0),self.u*(1-self.d)/(self.v*(1+2*self.d)))
    def locus(self):
        z0=self.cut();E0=self.reduced(z0)[2];s=self.a+self.b;t=(self.a+2*self.b)*(self.e-1)/(self.e*s*s)-2
        if E0>0:status='no admissible nonnegative-fork-current equilibrium';inside=False;margin=None
        elif self.e<=1 or t<=z0:status='disguised toric';inside=True;margin=None
        else:
            margin=s-self.reduced(t)[0][0];inside=margin>=0
            status='disguised toric' if inside else 'reservoir current budget violated'
        return dict(inside=inside,status=status,z0=z0,E_z0=E0,t=t,threshold_margin=margin)
    def admissible_root(self,steps=ROOT_BISECTIONS):
        """Exact rational enclosure; no claim about roots below the admissible cut."""
        lo=self.cut()
        if self.reduced(lo)[2]>0:raise ValueError('No root at or above the nonnegative-current cut; other stationary states are not excluded.')
        if self.reduced(lo)[2]==0:return lo,lo
        hi=lo+1+(1+self.e)*(self.a+2*self.b)/self.b
        if self.reduced(hi)[2]<=0:raise ArithmeticError('Analytic upper bracket failed.')
        for _ in range(steps):
            mid=(lo+hi)/2;residual=self.reduced(mid)[2]
            if residual==0:return mid,mid
            if residual<0:lo=mid
            else:hi=mid
        return lo,hi
    def permanence_bounds(self):
        C=(self.e+Q(3,2))*self.a+(2*self.e+1)*self.b+(self.u+2)**2/(8*self.v)
        rho=min(1/(self.e+Q(3,2)),1/(2*self.e+1),Q(1),self.d);M=2*C/rho
        ma=self.a/(2*(2+2*self.e*M));mb=self.b/(2*(1+self.e+M));mz=ma/(2*(M+self.u+2*self.v*M));mh=self.u*mz/(2*(2+self.d))
        return dict(C=C,rho=rho,eventual_upper=M,eventual_lower=min(ma,mb,mz,mh),individual_lowers=(ma,mb,mz,mh))


def currents(p,x):
    A,B,z,H=map(Q,x);K=A-B*z;J=p.e*(B-A*A)
    return dict(K=K,J=J,B_minus_J=B-J,AB_margins=(-K+2*J,K-J),zH_margins=(-p.u*z+3*H-2*p.v*z*z,p.u*z-2*H+p.v*z*z))


class SupportingCertificate:
    def __init__(self,kind):
        if kind=='K':self.h=((0,0,1,0),)+((0,1,0,0),)*5+((0,0,1,0),)*2;self.pi=(0,0,-1,0,0,0,0,0)
        elif kind=='J':self.h=((1,1,1,1),(1,1,0,1),(-1,-1,0,1),(1,1,0,1),(1,1,0,1),(1,1,-1,-1),(-1,-1,1,1),(-1,-1,1,1));self.pi=(0,-1,-1,0,-1,0,-1,0)
        else:raise ValueError('Certificate kind must be K or J.')
        self.kind=kind
    def extension(self,extra=()):
        vertices=list(COMPLEXES);h=list(self.h);pi=list(self.pi)
        for point in extra:
            if len(point)!=4:raise ValueError('Auxiliary exponent must have four coordinates.')
            values=[p-sum(a*(b-c) for a,b,c in zip(slope,point,y)) for p,slope,y in zip(self.pi,self.h,COMPLEXES)]
            active=max(range(8),key=lambda i:values[i]);vertices.append(tuple(point));h.append(self.h[active]);pi.append(values[active])
        slacks=[[sum(h[i][k]*(y[k]-x[k]) for k in range(4))+pi[j]-pi[i] for j,y in enumerate(vertices)] for i,x in enumerate(vertices)]
        if min(min(row) for row in slacks)<0:raise ArithmeticError('Supporting inequality failed.')
        return slacks
    def pairing(self,p,state):
        x=tuple(map(Q,state));coeff=p.literal().coefficients()
        return sum(math.prod(xi**yi for xi,yi in zip(x,y))*sum(hi*ci for hi,ci in zip(h,coeff[y])) for y,h in zip(COMPLEXES,self.h))


class BalancedRealization:
    """Exact witness at an exact rational stationary state; never round a failed budget."""
    def __init__(self,p,state):
        self.p=p;self.state=tuple(map(Q,state));A,B,z,H=self.state
        if min(self.state)<=0 or any(p.literal().exact_field(self.state)):raise ValueError('Exact positive stationary state required; a rounded numerical root is not an exact witness.')
        values=currents(p,self.state);K=values['K'];J=values['J']
        if K<0 or J>B:raise ValueError('Supporting-plane budget excludes every finite complex-balanced realization.')
        w=B*z;U=p.u*z;V=p.v*z*z;jp=max(J,Q(0));jm=max(-J,Q(0));tau=max(Q(0),V-Q(3,2)*H)
        self.fluxes={(0,1):p.a,(0,6):p.b,(1,2):w,(2,1):w,(1,6):K,(1,3):K,(1,0):w+jm,(1,7):jm,
            (6,7):p.e*B-jp,(6,1):2*jp,(6,0):B-jp,(7,6):p.e*A*A,(3,4):U,(3,5):tau,(3,0):tau,
            (5,4):V,(4,3):3*H-2*V+2*tau,(4,5):V-tau,(4,0):K-tau}
        if min(self.fluxes.values())<0:raise ArithmeticError('Negative realizing flux.')
        self.network=MassActionNetwork(SPECIES,[Reaction(COMPLEXES[i],COMPLEXES[j],q/math.prod(x**y for x,y in zip(self.state,COMPLEXES[i]))) for (i,j),q in self.fluxes.items() if q>0])
        if any(self.network.balance(self.state).values()) or self.network.coefficients()!=p.literal().coefficients():raise ArithmeticError('Complex balance or full polynomial equality failed.')
    def entropy(self,x):
        x=np.asarray(x,float);ref=np.array(self.state,float);return float(np.sum(x*np.log(x/ref)-x+ref))
    def entropy_derivative(self,x):return float(np.log(np.asarray(x,float)/np.array(self.state,float))@self.network.field(x))


@dataclass(frozen=True)
class CatalyticAttachment:
    name:str
    catalyst:tuple
    forward:Q
    reverse:Q
    def __post_init__(self):
        object.__setattr__(self,'forward',Q(self.forward));object.__setattr__(self,'reverse',Q(self.reverse))
        if len(self.catalyst)!=4 or any(type(i) is not int or i<0 for i in self.catalyst) or min(self.forward,self.reverse)<=0:raise ValueError('Base-only natural catalyst exponents and positive rates required.')
    @property
    def stationary(self):return self.forward/self.reverse
    @staticmethod
    def compose(base,attachments):
        attachments=tuple(attachments);m=len(attachments)
        if base.species!=SPECIES:raise ValueError('Attach to the four-species base model.')
        reactions=[Reaction(r.source+(0,)*m,r.target+(0,)*m,r.rate) for r in base.reactions]
        for j,a in enumerate(attachments):
            source=a.catalyst+(0,)*m;target=a.catalyst+tuple(int(i==j) for i in range(m))
            reactions.extend([Reaction(source,target,a.forward),Reaction(target,source,a.reverse)])
        return MassActionNetwork(SPECIES+tuple(a.name for a in attachments),reactions)


def slice_certificate():
    eta=sp.Symbol('eta',positive=True)
    A,B,z,H=sp.symbols('A B z H');a=(19-2*eta)/25;b=(eta-4)/25
    F=sp.Matrix([a-2*A+B*z+2*eta*(B-A*A),b+A-(1+z)*B-eta*(B-A*A),A-B*z-z-10*z*z+3*H,z+5*z*z-5*H])
    state={A:sp.Rational(2,5),B:sp.Rational(1,5),z:sp.Rational(1,5),H:sp.Rational(2,25)}
    if any(sp.expand(v) for v in F.subs(state)):raise ArithmeticError('Symbolic slice stationarity failed.')
    J=F.jacobian([A,B,z,H]).subs(state)
    coeff=J.charpoly().all_coeffs();_,c1,c2,c3,c4=coeff;quantities=coeff[1:]+[sp.expand(c1*c2-c3),sp.expand(c1*c2*c3-c3*c3-c1*c1*c4)]
    if not all(v>0 for expression in quantities for v in sp.Poly(expression,eta).all_coeffs()):raise ArithmeticError('All-positive Hurwitz coefficients failed.')
    D=sp.diag(sp.Rational(5,2),5,5,sp.Rational(25,2));J5=J.subs(eta,5);energy=-(D*J5+J5.T*D)/2;minors=[energy[:i,:i].det() for i in range(1,5)]
    if any(v<=0 for v in minors):raise ArithmeticError('Boundary quadratic entropy certificate failed.')
    return dict(characteristic=list(map(str,coeff)),hurwitz=list(map(str,quantities)),boundary_energy_minors=list(map(str,minors)),scope='Local stability for the entire positive slice; the paper proves a non-explicit open globally stable nontoric neighborhood, not every off-locus parameter.')


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)
    p=Parameters(*PARAMETERS);network=p.literal();classification=p.locus();root=None;witness=None
    if EXACT_REFERENCE is not None:
        reference=tuple(map(Q,EXACT_REFERENCE))
        if any(network.exact_field(reference)):raise ValueError('Update EXACT_REFERENCE or set it to None after changing PARAMETERS.')
        if classification['inside']:witness=BalancedRealization(p,reference)
    else:
        root=p.admissible_root();reference=p.reduced(sum(root)/2)[0]
    times=np.linspace(0,TIME_HORIZON,TIME_SAMPLES);trajectories=[];differences=[];realization_difference=None
    for initial in INITIAL_STATES:
        values=network.integrate(initial,times);independent=network.integrate(initial,times,'BDF');trajectories.append(values);differences.append(float(np.max(abs(values-independent))))
    if witness:realization_difference=float(np.max(abs(witness.network.integrate(INITIAL_STATES[0],times)-trajectories[0])))
    certificates={}
    for kind in ('K','J'):
        certificate=SupportingCertificate(kind);slacks=certificate.extension();expanded=certificate.extension(((9,4,7,11),(100,0,23,5),(0,0,0,12)))
        certificates[kind]=dict(slacks=slacks,off_diagonal_zeros=sum(slacks[i][j]==0 for i in range(8) for j in range(8) if i!=j),pairing=certificate.pairing(p,reference),auxiliary_pairs_checked=len(expanded)**2)
    slice_state=tuple(map(Q,('2/5','1/5','1/5','2/25')));slice_rows=[]
    for eta in SLICE_ETAS:
        q=Parameters.slice(eta);margins=currents(q,slice_state);slice_rows.append((eta,int(q.locus()['inside']),*map(str,margins['AB_margins']),str(margins['B_minus_J']),float(max(np.linalg.eigvals(q.literal().jacobian(np.array(slice_state,float))).real))))
    grid=[];elo,ehi=map(Q,GRID_E[:2]);dlo,dhi=map(Q,GRID_D[:2])
    for j in range(GRID_D[2]):
        d=dlo+(dhi-dlo)*j/(GRID_D[2]-1)
        for i in range(GRID_E[2]):
            e=elo+(ehi-elo)*i/(GRID_E[2]-1);entry=Parameters(Q(1,10),Q(1,10),5,1,e,d).locus();category=0 if entry['inside'] else (1 if entry['E_z0']>0 else 2)
            grid.append((float(e),float(d),category))
    attachments=[CatalyticAttachment(name,tuple(catalyst),a,b) for name,catalyst,a,b,initial in ATTACHMENTS]
    augmented=CatalyticAttachment.compose(network,attachments);attachment_initial=list(INITIAL_STATES[0])+[float(Q(row[4])) for row in ATTACHMENTS]
    attachment_states=augmented.integrate(attachment_initial,times)
    if witness:
        extended=CatalyticAttachment.compose(witness.network,attachments);xstar=witness.state+tuple(a.stationary for a in attachments)
        if extended.coefficients()!=augmented.coefficients() or any(extended.balance(xstar).values()):raise ArithmeticError('Attachment realization failed.')
    def write_json(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)
    table('trajectories.csv',['initial_index','time',*SPECIES,'entropy','entropy_derivative'],[(i,t,*x,witness.entropy(x) if witness else '',witness.entropy_derivative(x) if witness else '') for i,values in enumerate(trajectories) for t,x in zip(times,values)])
    table('attachments.csv',['time',*augmented.species],np.column_stack([times,attachment_states]))
    table('stable_slice.csv',['eta','disguised_toric','AB_A_margin','AB_B_margin','B_minus_J','spectral_abscissa'],slice_rows)
    table('parameter_locus.csv',['e','d','category_0_toric_1_fork_2_reservoir'],grid)
    if witness:
        table('realization.csv',['source','target','equilibrium_flux','rate','present'],[(COMPLEX_NAMES[i],COMPLEX_NAMES[j],str(q),str(q/math.prod(x**y for x,y in zip(witness.state,COMPLEXES[i]))),int(q>0)) for (i,j),q in witness.fluxes.items()])
    elif (out/'realization.csv').exists():
        (out/'realization.csv').unlink()  # This run has no exact witness; remove only its stale optional export.
    results=dict(parameters=PARAMETERS,classification=classification,reference=reference,reference_scope='exact rational stationary state' if EXACT_REFERENCE is not None else 'midpoint lift of exact root enclosure; not an exact witness',root_enclosure=root,
        currents=currents(p,reference),supporting_certificates=certificates,permanence=p.permanence_bounds(),slice=slice_certificate(),
        exact_realization_built=witness is not None,active_realization_edges=len(witness.network.reactions) if witness else None,
        numerical_checks=dict(independent_solver_differences=differences,equivalent_field_solver_difference=realization_difference,attachment_base_difference=float(np.max(abs(attachment_states[:,:4]-trajectories[0]))),minimum_concentration=float(min(x.min() for x in trajectories))),
        scope='Exact rational classification, supporting inequalities and witness reconstruction; numerical trajectories illustrate conventional permanence/global-stability theorems. Lean not rerun. Deterministic equivalence does not identify stochastic processes.')
    write_json('results.json',results)
    lines=[f'Parameter classification: {classification["status"]}.',f'Exact rational realization built: {witness is not None}; active edges: {results["active_realization_edges"]}.',
        f'Current budgets: K={results["currents"]["K"]}, B-J={results["currents"]["B_minus_J"]}; AB production margins={results["currents"]["AB_margins"]}.',
        'All 128 source-pair supporting inequalities verified; active-piece extension checked on three auxiliary complexes.',
        'Slice: realizable iff eta <= 5, strictly AB-productive iff 9/2 < eta < 9; locally stable throughout 4 < eta < 19/2.',
        f'Independent solver differences: {differences}.', 'No explicit globally stable nontoric neighborhood radius is supplied. No stochastic equivalence claimed.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    from matplotlib.colors import ListedColormap
    fig,axs=plt.subplots(1,2,figsize=(10,4.3),layout='constrained')
    g=np.array(grid);axs[0].pcolormesh(g[:,0].reshape(GRID_D[2],GRID_E[2]),g[:,1].reshape(GRID_D[2],GRID_E[2]),g[:,2].reshape(GRID_D[2],GRID_E[2]),cmap=ListedColormap(['#3e8585','#cf8952','#866b9c']),shading='nearest',vmin=0,vmax=2)
    axs[0].set(xlabel='Reverse-channel rate e',ylabel='Output rate d',title='Existence of an equivalent complex-balanced\nmodel')
    from matplotlib.patches import Patch
    axs[0].legend(handles=[Patch(color=c,label=s) for c,s in zip(['#3e8585','#cf8952','#866b9c'],['Realizable','Fork budget fails','Reservoir budget fails'])],fontsize=7)
    etas=np.linspace(4.01,9.49,200);abscissae=[max(np.linalg.eigvals(Parameters.slice(str(e)).literal().jacobian(np.array(slice_state,float))).real) for e in etas]
    axs[1].plot(etas,abscissae);axs[1].axvline(5,color='#866b9c',ls='--',label='Realizability boundary');axs[1].axhline(0,color='gray',lw=.7);axs[1].set(xlabel='Slice parameter eta',ylabel='Spectral abscissa',title='Equilibrium stability across the\nrealizability boundary');axs[1].legend(fontsize=7);axs[1].grid(alpha=.2)
    fig.savefig(out/'classification.png',dpi=180);fig.savefig(out/'classification.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4.3),layout='constrained')
    for i,name in enumerate(SPECIES):
        line,=axs[0].plot(times,trajectories[0][:,i],label=name)
        for values in trajectories[1:]:axs[0].plot(times,values[:,i],ls='--',color=line.get_color(),lw=.8)
    axs[0].set(xlabel='Time',ylabel='Concentration',title=f'{len(trajectories)} positive initial states');axs[0].legend()
    if witness:
        for i,values in enumerate(trajectories):axs[1].semilogy(times,np.maximum([witness.entropy(x) for x in values],1e-16),label=f'Initial state {i+1}')
        axs[1].set(ylabel='Horn–Jackson entropy (display floor 1e-16)',title='Horn–Jackson entropy along reactor\ntrajectories')
    else:
        for i,a in enumerate(attachments):axs[1].plot(times,attachment_states[:,4+i],label=a.name)
        axs[1].set(ylabel='Private-species concentration',title='Catalytic attachments')
    axs[1].set_xlabel('Time');axs[1].legend()
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'trajectories.png',dpi=180);fig.savefig(out/'trajectories.svg');plt.close(fig)
    digest=lambda path:hashlib.sha256(path.read_bytes()).hexdigest()
    write_json('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()