"""Exact stationary reductions and reusable mass-action models for mixed degradation.

Run: python example.py --output outputs
Inputs are dimensionless constructive examples, not fitted chemical constants.
"""
from __future__ import annotations

# EDITABLE INPUTS: zero is a valid loss; all reversible rates must stay positive.
TYPE_II_ACTIVE_WEIGHTS = ((2,), (1,), (3,))  # one or more active edges per gap
TYPE_II_LOSS_BY_NAME = {"F0": "1/5", "T1": "0", "F1": "0", "T2": "1/3", "F2": "0", "T0": "0"}
TARGET_CONCENTRATION = "1"  # construct rates with this exact positive state
MIN_ONE_WAY_FLOW = "1"     # positive margin above the admissibility boundary
LOSS_SWEEP = tuple(i / 40 for i in range(41))  # fixed rates, d = scale * target d
TYPE_V_PATH_EDGES = (3, 0, 2)  # 0 identifies U_i=S_i; positive values give unit paths
TYPE_V_INTERNAL_LOSS = "1/4"
TYPE_V_BASE_LOSS = ("0", "1/5", "0")
TYPE_V_FORK_LOSS = ("0", "0", "1/3")
MAX_PRINCIPAL_MINORS = 1024
MANUSCRIPT_SHA256 = "6d260d5d6596b099430287d32563d91ebf3f22ccbee5c5b9a889b54470ccf3ed"

import argparse
import csv
from dataclasses import dataclass, replace
import hashlib
import itertools
import json
from pathlib import Path
import platform

import numpy as np
import sympy as sp
from scipy.optimize import least_squares

Q = sp.Rational


def vector(values):
    return sp.Matrix([Q(v) for v in values])


def positive(values, name, allow_zero=False):
    if any(v < 0 if allow_zero else v <= 0 for v in values):
        raise ValueError(f"{name} must be {'nonnegative' if allow_zero else 'positive'}.")


@dataclass(frozen=True)
class UnitPath:
    """Exact stationary two-port; only one-to-one reactions may be put here."""
    forward: tuple
    reverse: tuple
    internal_loss: tuple

    def __post_init__(self):
        for field in ('forward', 'reverse', 'internal_loss'):
            object.__setattr__(self, field, tuple(Q(v) for v in getattr(self, field)))
        if not self.forward or len(self.forward) != len(self.reverse) or len(self.internal_loss) != len(self.forward)-1:
            raise ValueError('A path needs m+1 reversible edges and m internal losses.')
        positive(self.forward + self.reverse, 'Path rates')
        positive(self.internal_loss, 'Internal losses', True)

    def compress(self):
        c, b, left, right = self.forward[0], self.reverse[0], Q(0), Q(0)
        for a, beta, d in zip(self.forward[1:], self.reverse[1:], self.internal_loss):
            pivot = b + right + a + d
            c, b, left, right = c*a/pivot, b*beta/pivot, left+c*(right+d)/pivot, beta*(right+d)/pivot
        return c, b, left, right

    def reconstruct(self, X, Y):
        """Solve the independent tridiagonal balance, not the compression recursion."""
        X, Y = Q(X), Q(Y)
        positive((X,Y), 'Endpoints')
        m = len(self.internal_loss)
        if not m:
            return sp.Matrix([X,Y])
        matrix, rhs = sp.zeros(m), sp.zeros(m,1)
        for k in range(m):
            matrix[k,k] = self.reverse[k] + self.forward[k+1] + self.internal_loss[k]
            if k: matrix[k,k-1] = -self.forward[k]
            else: rhs[k] += self.forward[0]*X
            if k < m-1: matrix[k,k+1] = -self.reverse[k+1]
            else: rhs[k] += self.reverse[-1]*Y
        return sp.Matrix([X, *matrix.inv()*rhs, Y])


@dataclass(frozen=True)
class SourceNetwork:
    """Column r consumes one X_r and produces column r of P."""
    names: tuple
    P: sp.ImmutableMatrix
    family: str = 'unclassified'

    def __post_init__(self):
        object.__setattr__(self, 'P', sp.ImmutableMatrix(self.P))
        if len(set(self.names)) != len(self.names) or self.P.shape != (len(self.names),)*2:
            raise ValueError('Unique species and square product matrix required.')
        if any(v.is_integer is not True or v < 0 for v in self.P):
            raise ValueError('Product multiplicities must be nonnegative integers.')
        if self.N.det() == 0: raise ValueError('This stationary-flow interface requires invertible N.')

    @property
    def N(self): return self.P-sp.eye(len(self.names))

    def monomials(self, x):
        return sp.Matrix([sp.prod(x[i]**self.P[i,r] for i in range(len(x))) for r in range(len(x))])

    def realize(self, x, d, margin=MIN_ONE_WAY_FLOW):
        """Construct exact rates at a chosen root; permits backward net currents."""
        x, d = vector(x), vector(d)
        if len(x) != len(self.names) or len(d) != len(x): raise ValueError('Wrong vector length.')
        positive(x, 'Concentrations'); positive(d, 'Losses', True); positive((Q(margin),), 'Flow margin')
        currents = self.N.inv()*sp.diag(*d)*x
        q = vector([max(0,-j)+Q(margin) for j in currents]); p = q+currents
        return MassActionModel(self, tuple(p[i]/x[i] for i in range(len(x))),
            tuple(q[i]/m for i,m in enumerate(self.monomials(x))), tuple(d))


def separated_type_ii(weights=TYPE_II_ACTIVE_WEIGHTS):
    """Arbitrary nonempty active gaps; each terminal return is a unit edge."""
    if len(weights) < 3 or any(not gap or any(int(m)!=m or m<1 for m in gap) for gap in weights):
        raise ValueError('At least three forks, nonempty gaps, positive integer active weights.')
    names=[]; sectors=[]
    for j,gap in enumerate(weights):
        nodes=[f'F{j}', *[f'G{j}_{k}' for k in range(1,len(gap))], f'T{(j+1)%len(weights)}']
        names.extend(nodes); sectors.append(nodes)
    index={s:i for i,s in enumerate(names)}; P=sp.zeros(len(names))
    for j,(gap,nodes) in enumerate(zip(weights,sectors)):
        for k,m in enumerate(gap): P[index[nodes[k+1]],index[nodes[k]]] = m
        P[index[f'T{j}'],index[f'F{j}']] += 1
        P[index[f'F{(j+1)%len(weights)}'],index[nodes[-1]]] = 1
    return SourceNetwork(tuple(names), P, 'separated_type_ii_skeleton')


def coincident_type_ii(weights=(2,1,3)):
    if len(weights)!=3 or any(int(m)!=m or m<1 for m in weights): raise ValueError('Three positive integer weights required.')
    P=sp.zeros(3)
    for i,m in enumerate(weights): P[(i+1)%3,i]=m; P[(i+2)%3,i]=1
    return SourceNetwork(('X0','X1','X2'), P, 'coincident_type_ii')


@dataclass(frozen=True)
class MassActionModel:
    source: SourceNetwork
    forward: tuple
    reverse: tuple
    loss: tuple

    def __post_init__(self):
        for name in ('forward','reverse','loss'):
            object.__setattr__(self,name,tuple(Q(v) for v in getattr(self,name)))
            if len(getattr(self,name)) != len(self.source.names): raise ValueError('Wrong rate vector length.')
        positive(self.forward+self.reverse, 'Reversible rates'); positive(self.loss,'Losses',True)

    def flows(self,x):
        return sp.diag(*self.forward)*sp.Matrix(x), sp.diag(*self.reverse)*self.source.monomials(x), sp.diag(*self.loss)*sp.Matrix(x)

    def field(self,x):
        p,q,e=self.flows(x)
        return self.source.N*(p-q)-e

    def scaled_jacobian(self,x):
        p,q,e=self.flows(x); N=self.source.N
        return N*sp.diag(*p)-N*sp.diag(*q)*self.source.P.T-sp.diag(*e)

    def detailed_balance(self):
        """Positive zero-loss root for these reversible constants; not the mixed-loss root."""
        N=np.array(self.source.N,float)
        return np.exp(np.linalg.solve(N.T,np.log(np.array(self.forward,float)/np.array(self.reverse,float))))

    def numeric(self,x):
        x=np.asarray(x,float)
        if np.any(x<=0) or not np.all(np.isfinite(x)): raise ValueError('Numerical state must be finite and positive.')
        P=np.array(self.source.P,float); N=P-np.eye(len(x))
        p=np.array(self.forward,float)*x; q=np.array(self.reverse,float)*np.exp(P.T@np.log(x)); e=np.array(self.loss,float)*x
        H=N@np.diag(p)-N@np.diag(q)@P.T-np.diag(e)
        return N@(p-q)-e, H

    def continuation(self,scales,initial=None):
        """Log-coordinate local solves. Failure is unresolved, never nonexistence."""
        guess=self.detailed_balance() if initial is None else np.array(initial,float)
        rows=[]
        for scale in scales:
            if scale<0: raise ValueError('Physical loss scales must be nonnegative.')
            model=replace(self,loss=tuple(Q(str(scale))*d for d in self.loss))
            # Dividing f by x avoids accepting the origin as a tiny-residual root.
            def residual(z):
                x=np.exp(z); f,_=model.numeric(x); return f/x
            def jac(z):
                x=np.exp(z); f,H=model.numeric(x); return H/x[:,None]-np.diag(f/x)
            try:
                fit=least_squares(residual,np.log(guess),jac=jac,bounds=(-25,25),xtol=1e-13,ftol=1e-13,gtol=1e-13,max_nfev=1000)
            except (ValueError, FloatingPointError, np.linalg.LinAlgError) as exc:
                rows.append({'scale':float(scale),'status':'unresolved','relative_residual':None,'reason':str(exc)}); break
            x=np.exp(fit.x); error=float(np.max(np.abs(residual(fit.x))))
            if not fit.success or error>1e-8:
                rows.append({'scale':float(scale),'status':'unresolved','relative_residual':error}); break
            f,H=model.numeric(x); physical=H/x[None,:]
            rows.append({'scale':float(scale),'status':'numerical_root','x':x.tolist(),
                'relative_residual':error,'spectral_abscissa':float(max(np.linalg.eigvals(physical).real)),
                'det_minus_H':float(np.linalg.det(-H))})
            guess=x
        return rows

    def expand_unit_edge(self,start,end,path):
        """Replace a literal unit reaction; active m>1 or fork edges are rejected."""
        i=self.source.names.index(start); j=self.source.names.index(end)
        if self.source.P[:,i] != sp.eye(len(self.loss))[:,j]: raise ValueError('Only a standalone unit edge may be expanded.')
        extra=tuple(f'{start}_path{k}' for k in range(1,len(path.forward)))
        if set(extra)&set(self.source.names): raise ValueError('Path names collide.')
        names=self.source.names+extra; n=len(names); P=sp.zeros(n)
        P[:len(self.loss),:len(self.loss)]=self.source.P; P[:,i]=sp.zeros(n,1)
        route=[i,*range(len(self.loss),n),j]; a=list(self.forward)+[Q(1)]*len(extra); b=list(self.reverse)+[Q(1)]*len(extra)
        for k,(u,v) in enumerate(zip(route,route[1:])):
            P[v,u]=1; a[u]=path.forward[k]; b[u]=path.reverse[k]
        model=MassActionModel(SourceNetwork(names,P,'expanded_unit_path'),tuple(a),tuple(b),self.loss+path.internal_loss)
        c,beta,left,right=path.compress(); ca=list(self.forward); cb=list(self.reverse); cd=list(self.loss)
        ca[i]=c; cb[i]=beta; cd[i]+=left; cd[j]+=right
        return model,replace(self,forward=tuple(ca),reverse=tuple(cb),loss=tuple(cd))


class StationaryCertificate:
    """Exact identities on a stationary flow triple; no division by a loss."""
    def __init__(self,model,x):
        self.model=model; self.x=vector(x); positive(self.x,'State')
        if model.field(self.x)!=sp.zeros(len(x),1): raise ValueError('Exact stationarity required; an approximate root is not a certificate.')
        self.p,self.q,self.e=model.flows(self.x); self.N=model.source.N; V=self.N.inv(); self.J=V*self.e
        self.shift=vector([max(0,-j) for j in self.J]); self.t=self.q-self.shift
        self.B=(sp.diag(*self.J)-V*sp.diag(*self.e))*V.T
        self.A=sp.diag(*self.shift)-self.B; self.H=model.scaled_jacobian(self.x)
        if -self.H != self.N*(sp.diag(*self.q)-self.B)*self.N.T: raise AssertionError('Factorization failed.')

    def principal_audit(self,budget=MAX_PRINCIPAL_MINORS):
        n=len(self.x)
        if 2**n>budget: return {'status':'unresolved_budget','required_minors':2**n}
        minors=[]; expansion=Q(0)
        for mask in range(2**n):
            idx=[i for i in range(n) if mask>>i&1]
            value=self.A.extract(idx,idx).det() if idx else Q(1)
            minors.append(value); expansion+=value*sp.prod(self.t[i] for i in range(n) if i not in idx)
        if expansion!=(self.A+sp.diag(*self.t)).det(): raise AssertionError('Multiaffine expansion failed.')
        lower=self.N.det()**2*sp.prod(self.t); actual=(-self.H).det()
        return {'status':'all_principal_minors_nonnegative' if min(minors)>=0 else 'negative_principal_minor',
            'count':len(minors),'zero_minors':sum(m==0 for m in minors),'minimum':str(min(minors)),
            'det_minus_H':str(actual),'lower_bound_from_nonnegative_minors':str(lower) if min(minors)>=0 else None,
            'ratio_to_lower_bound':str(actual/lower) if min(minors)>=0 else None,
            'scope':'Exact certificate for this triple. The paper proves the universal bound for separated Type II skeletons only.'}


@dataclass(frozen=True)
class TypeVTopology:
    """All eight identifications and arbitrary finite linked paths."""
    edges: tuple

    def __post_init__(self):
        if len(self.edges)!=3 or any(int(n)!=n or n<0 for n in self.edges): raise ValueError('Three nonnegative integer edge counts required.')

    def names_for_pair(self,i):
        return (f'U{i}',) if self.edges[i]==0 else (f'U{i}',*[f'Z{i}_{k}' for k in range(1,self.edges[i])],f'S{i}')

    def source(self):
        names=tuple(s for i in range(3) for s in self.names_for_pair(i)); index={s:i for i,s in enumerate(names)}; P=sp.zeros(len(names))
        for i in range(3):
            route=self.names_for_pair(i)
            for a,b in zip(route,route[1:]): P[index[b],index[a]]=1
            for j in range(3):
                if j!=i: P[index[f'U{j}'],index[route[-1]]]=1
        return SourceNetwork(names,P,'type_v')

    def reduce(self,model):
        if model.source.names!=self.source().names or model.source.P!=self.source().P: raise ValueError('Model does not match this Type V topology.')
        return ReducedTypeV(self,model)


class ReducedTypeV:
    def __init__(self,topology,model):
        self.topology=topology; self.model=model; self.sectors=[]
        self.alpha=[]; self.beta=[]; self.lam=[]; self.mu=[]
        for i in range(3):
            route=topology.names_for_pair(i); ids=[model.source.names.index(s) for s in route]; last=ids[-1]
            f,g=model.forward[last],model.reverse[last]; d=model.loss[ids[0]]
            if len(ids)==1:
                self.sectors.append(None); alpha,beta,lam,mu=f,g,d,Q(0)
            else:
                path=UnitPath(tuple(model.forward[k] for k in ids[:-1]),tuple(model.reverse[k] for k in ids[:-1]),tuple(model.loss[k] for k in ids[1:-1]))
                a,b,left,right=path.compress(); d+=left; h=model.loss[last]+right; D=b+f+h
                self.sectors.append((path,a,g,D)); alpha,beta,lam,mu=f*a/D,g*(b+h)/D,d+h*a/D,h*g/D
            self.alpha.append(alpha); self.beta.append(beta); self.lam.append(lam); self.mu.append(mu)

    def G(self,u):
        return sp.Matrix([2*self.beta[i]*u[(i+1)%3]*u[(i+2)%3]/u[i]+self.lam[(i+1)%3]*u[(i+1)%3]/u[i]+self.lam[(i+2)%3]*u[(i+2)%3]/u[i]+self.mu[(i+1)%3]*u[(i+2)%3]+self.mu[(i+2)%3]*u[(i+1)%3] for i in range(3)])

    def residual(self,u): return sp.diag(*u)*(2*sp.Matrix(self.alpha)-self.G(u))

    def reconstruct(self,u):
        u=vector(u); positive(u,'Base concentrations'); values=[]
        for i,sector in enumerate(self.sectors):
            if sector is None: values.append(u[i])
            else:
                path,a,g,D=sector; S=(a*u[i]+g*u[(i+1)%3]*u[(i+2)%3])/D
                values.extend(path.reconstruct(u[i],S))
        return sp.Matrix(values)

    def tangent(self,u):
        z=sp.symbols('u0:3',positive=True)
        return self.G(z).jacobian(z).subs(dict(zip(z,u)))*sp.diag(*u)

    def ratio_witness(self,u,v):
        """A strict row defect shows unequal positive u,v cannot both be roots."""
        u,v=vector(u),vector(v); positive(u,'u');positive(v,'v'); r=[u[i]/v[i] for i in range(3)]
        if all(a==1 for a in r): return {'status':'identical'}
        upward=sorted(r)[1]>=1; i=(min if upward else max)(range(3),key=lambda k:r[k]); j,k=(i+1)%3,(i+2)%3
        terms=[2*self.beta[i]*v[j]*v[k]*(r[j]*r[k]-r[i]),self.lam[j]*v[j]*(r[j]-r[i]),self.lam[k]*v[k]*(r[k]-r[i]),self.mu[j]*v[i]*v[k]*r[i]*(r[k]-1),self.mu[k]*v[i]*v[j]*r[i]*(r[j]-1)]
        direction=1 if upward else -1
        if direction*terms[0]<=0 or any(direction*t<0 for t in terms): raise AssertionError('Ratio sign structure failed.')
        return {'status':'cannot_both_be_stationary','row':i,'terms':list(map(str,terms)),'strict_defect':str(sum(terms))}


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)
    source=separated_type_ii(); x=vector([TARGET_CONCENTRATION]*len(source.names)); d=vector([TYPE_II_LOSS_BY_NAME.get(n,'0') for n in source.names])
    model=source.realize(x,d); cert=StationaryCertificate(model,x); audit=cert.principal_audit(); sweep=model.continuation(LOSS_SWEEP)
    topology=TypeVTopology(TYPE_V_PATH_EDGES); vs=topology.source(); vd=[]
    for i in range(3):
        route=topology.names_for_pair(i)
        if len(route)==1: vd.append(Q(TYPE_V_BASE_LOSS[i])+Q(TYPE_V_FORK_LOSS[i]))
        else: vd.extend([Q(TYPE_V_BASE_LOSS[i]),*[Q(TYPE_V_INTERNAL_LOSS)]*(len(route)-2),Q(TYPE_V_FORK_LOSS[i])])
    vx=vector([TARGET_CONCENTRATION]*len(vs.names)); vm=vs.realize(vx,vd); reduced=topology.reduce(vm); u=vector([TARGET_CONCENTRATION]*3)
    if reduced.reconstruct(u)!=vx or reduced.residual(u)!=sp.zeros(3,1): raise AssertionError('Type V exact reconstruction failed.')
    # Internal-only loss survives as endpoint leakage; no loss division is used.
    path=UnitPath((1,1),(1,1),(1,)); path_values=path.reconstruct(2,3)
    coincident=coincident_type_ii(); cm=coincident.realize([1,1,1],[0,Q(1,3),0])
    result={'scope':'Dimensionless exact stationary examples and numerical local continuation. Universal uniqueness is the paper theorem, not a solver inference; Lean is not rerun.',
        'type_ii':{'species':source.names,'product_matrix':[[int(v) for v in row] for row in source.P.tolist()],
            'forward':list(map(str,model.forward)),'reverse':list(map(str,model.reverse)),'loss':list(map(str,d)),
            'exact_state':list(map(str,x)),'currents':list(map(str,cert.J)),'admissibility_shift':list(map(str,cert.shift)),
            'min_one_way_flows':list(map(str,cert.t)),'certificate':audit,'numerical_continuation':sweep},
        'type_v':{'edges':TYPE_V_PATH_EDGES,'species':vs.names,'exact_state':list(map(str,vx)),
            'forward':list(map(str,vm.forward)),'reverse':list(map(str,vm.reverse)),'loss':list(map(str,vd)),
            'alpha':list(map(str,reduced.alpha)),'beta':list(map(str,reduced.beta)),'lambda':list(map(str,reduced.lam)),'mu':list(map(str,reduced.mu)),
            'tangent_determinant':str(reduced.tangent(u).det()),'ratio_witness':reduced.ratio_witness([Q(3,2),1,Q(1,2)],u)},
        'internal_only_loss':{'compression':list(map(str,path.compress())),'endpoints_2_3_reconstruction':list(map(str,path_values)),
            'left_current':str(2-path_values[1]),'right_current':str(path_values[1]-3),'leakage':str(path_values[1])},
        'coincident_branch':{'det_N':str(coincident.N.det()),'inverse':[[str(a) for a in row] for row in coincident.N.inv().tolist()],
            'exact_root_verified':cm.field([1,1,1])==sp.zeros(3,1),'separated_determinant_bound_claimed':False}}
    def csv_file(name,headers,rows):
        with (out/name).open('w',newline='') as f:
            w=csv.writer(f);w.writerow(headers);w.writerows(rows)
    csv_file('loss_continuation.csv',['loss_scale','status',*source.names,'relative_residual','spectral_abscissa','det_minus_H'],
        [[r['scale'],r['status'],*r.get('x',['']*len(x)),r['relative_residual'],r.get('spectral_abscissa',''),r.get('det_minus_H','')] for r in sweep])
    csv_file('source_reactions.csv',['model','source','products','forward','reverse','degradation'],
        [[name,s,' + '.join(f'{network.P[i,r]} {network.names[i]}' for i in range(len(network.names)) if network.P[i,r]),str(m.forward[r]),str(m.reverse[r]),str(m.loss[r])]
        for name,m in [('Type II',model),('Type V',vm)] for network in [m.source] for r,s in enumerate(network.names)])
    path_rows=[]
    for loss in [Q(k,10) for k in range(31)]:
        p=UnitPath((1,1),(1,1),(loss,)); c,b,l,r=p.compress(); vals=p.reconstruct(2,3)
        path_rows.append([str(loss),str(c),str(b),str(l),str(r),str(vals[1]),str((c+l)*2-b*3),str(c*2-(b+r)*3)])
    csv_file('passive_path.csv',['internal_loss','c','beta','lambda_left','lambda_right','internal_concentration','left_current','right_current'],path_rows)
    (out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
    lines=[f'Type II: exact positive stationary state; {audit["count"]} principal minors checked.',
        f'det(-H)={audit["det_minus_H"]}; lower bound={audit["lower_bound_from_nonnegative_minors"]}.',
        f'Type V: {len(vs.names)} literal species reduce to three base variables with exact reconstruction.',
        'An internal-only loss gives two endpoint leaks of 1/3, even when both direct endpoint losses vanish.',
        f'Fixed-rate continuation: {sum(r["status"]=="numerical_root" for r in sweep)}/{len(LOSS_SWEEP)} numerical roots.',
        'At most one positive root is the theorem; existence, stability at arbitrary loss and global continuation are not asserted.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n'); print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained'); roots=[r for r in sweep if r['status']=='numerical_root']
    for i,n in enumerate(source.names): axes[0].plot([r['scale'] for r in roots],[r['x'][i] for r in roots],label=n)
    axes[0].set(xlabel='Scale of selective degradation (fixed reversible rates)',ylabel='Positive stationary concentration',title='Numerical local continuation'); axes[0].legend(ncol=2,fontsize=8)
    axes[1].plot([r['scale'] for r in roots],[r['spectral_abscissa'] for r in roots],color='#b95024'); axes[1].axhline(0,color='black',ls='--',lw=.8)
    axes[1].set(xlabel='Scale of selective degradation',ylabel='Largest real part of physical Jacobian eigenvalues',title='Stability diagnostic for these computed states')
    for ax in axes: ax.grid(alpha=.2)
    fig.savefig(out/'continuation.png',dpi=180); fig.savefig(out/'continuation.svg'); plt.close(fig)
    rows=np.array([[float(Q(v)) for v in row] for row in path_rows]); fig,axes=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
    axes[0].plot(rows[:,0],rows[:,1],label='Transmission c = beta'); axes[0].plot(rows[:,0],rows[:,3],label='Leakage lambda_left = lambda_right')
    axes[0].set(xlabel='Loss at internal species Z',ylabel='Effective coefficient',title='X <-> Z <-> Y: all four edge rates = 1'); axes[0].legend(fontsize=8)
    axes[1].plot(rows[:,0],rows[:,6],label='Current leaving X'); axes[1].plot(rows[:,0],rows[:,7],label='Current entering Y')
    axes[1].plot(rows[:,0],rows[:,0]*rows[:,5],'--',label='Their difference: internal loss')
    axes[1].set(xlabel='Loss at internal species Z',ylabel='Stationary current (X=2, Y=3)',title='Endpoint currents after internal-species\nelimination'); axes[1].legend(fontsize=8)
    for ax in axes: ax.grid(alpha=.2)
    fig.savefig(out/'path_reduction.png',dpi=180); fig.savefig(out/'path_reduction.svg'); plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    (out/'run_metadata.json').write_text(json.dumps({'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':digest(Path(__file__)),
        'python':platform.python_version(),'output_sha256':{p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}},indent=2)+'\n')


if __name__=='__main__': main()