"""Chemical copying as a complete cycle: return regions, allocation, recovery.

Three manuscript models remain distinct. Numerical trajectories do not replace
the paper's large generator certificates; scalar consequences are checked here.
"""
from __future__ import annotations

from dataclasses import dataclass
from fractions import Fraction as Q
from pathlib import Path
from functools import lru_cache
import argparse
import csv
import hashlib
import itertools
import json
import math
import platform
import mpmath as mp
import numpy as np
from scipy.integrate import solve_ivp

# EDITABLE INPUTS: constructed source (counts, seconds, prescribed volume).
BIRTH_MARKER = 53
CROSS_CATALYSIS = 0.05
REVERSE_GROWTH = 1e-7
DIVISION_DEADLINE = 20.0
PHYSICAL_EVENT_QUOTA = 100_000_000
COMPUTE_EVENT_LIMIT = 500_000   # unfinished is not a scientific failure
COPY_SEED = 17092026
PREPARATIONS = ((53,0),(530,53))
# Redundancy theorem: different four-species chemistry, different copy number.
MODULE_COUNT = 1
GROWTH_RATE = Q(1,10**11)
TARGET_ERROR = Q(1,100)
HORIZONS = tuple(range(1,41))
# Nominal eight-species chemistry: molar concentrations, seconds.
SPECIES = ('S','C','P','E','D','U','V','I')
FEED = (0.05,0,0,0,0.1,0,0,0.00231)
DILUTION = 0.002
EXCHANGE = 0.65
LIGATION = 0.41
HYDROLYSIS = 9.26e-6
INHIBITION = 150.0
RECOVERY_TIME = 500.0
OMEGA = 24*10**18             # molecules per molar; NOT directly simulated
PERTURBATION = 0.001           # nominal-region diagnostic, not uncertainty certification
MANUSCRIPT_SHA256 = '50450e40c8f43be86391b39a88460038b1d3d931ce7d2bba0be788916d1af07e'
REGION_SHA256 = 'e214ee720185c3e078d4880be7bfa247a473c4084165182477be5f7fae379e5a'


@dataclass(frozen=True)
class CountBox:
    lower: tuple[int, ...]
    upper: tuple[int, ...]

    def __post_init__(self):
        if len(self.lower) != len(self.upper) or any(not 0 <= lo <= hi for lo,hi in zip(self.lower,self.upper)):
            raise ValueError('Nonempty nonnegative coordinate intervals required.')

    def contains(self, counts):
        return len(counts) == len(self.lower) and all(lo <= n <= hi for n,lo,hi in zip(counts,self.lower,self.upper))


class ComplementaryAllocation:
    @staticmethod
    @lru_cache(maxsize=4096)
    def coordinate_return(total, lower, upper):
        """Exact P(lower <= D,total-D <= upper), D~Bin(total,1/2)."""
        if not isinstance(total,int) or not 0 <= total <= 20000 or lower > upper:
            raise ValueError('Exact allocation guard: integer total in [0,20000].')
        lo,hi = max(0,lower,total-upper),min(total,upper,total-lower)
        return Q(sum(math.comb(total,i) for i in range(lo,hi+1)),2**total)

    @classmethod
    def joint_return(cls, parent, box):
        if len(parent) != len(box.lower):
            raise ValueError('Allocation dimension differs from admission box.')
        return math.prod(cls.coordinate_return(int(n),lo,hi) for n,lo,hi in zip(parent,box.lower,box.upper))

    @staticmethod
    def draw(parent,rng):
        first = rng.binomial(np.asarray(parent,dtype=np.int64),.5)
        return first,np.asarray(parent,dtype=np.int64)-first


class ConstructedSource:
    boxes = (CountBox((6,0),(106,8)),CountBox((371,6),(742,132)))

    def __init__(self,epsilon=CROSS_CATALYSIS,rho=REVERSE_GROWTH,birth=BIRTH_MARKER):
        if epsilon < 0 or rho < 0 or birth < 1:
            raise ValueError('Nonnegative interactions and positive marker required.')
        self.epsilon,self.rho,self.birth = epsilon,rho,birth
        jumps = []
        kinds = []
        resident = ((1,0),(-1,0),(2,-1),(-2,1),(-1,1),(1,-1),(-1,0),(1,0),(-1,0),(1,0))
        for i in range(2):
            for dx,dy in resident:
                jump = [0]*5; jump[2*i:2*i+2] = [dx,dy]
                jumps.append(jump);kinds.append('resident')
            for dx,dm,kind in ((-1,1,'growth'),(1,-1,'reverse_growth')):
                jump = [0]*5;jump[2*i]=dx;jump[4]=dm
                jumps.append(jump);kinds.append(kind)
            for dx,dy in ((-1,1),(1,-1)):
                jump = [0]*5;jump[2*i:2*i+2]=[dx,dy]
                jumps.append(jump);kinds.append('cross_catalysis')
        jumps.append([0]*5);kinds.append('self_event')
        self.jumps,self.kinds = np.array(jumps,dtype=np.int64),tuple(kinds)

    def rates(self,state,fluid=False):
        m = state[4]
        if m <= 0:
            raise ValueError('Positive marker required.')
        rates = []
        for i in range(2):
            x,y = map(float,state[2*i:2*i+2]);other_y=float(state[2*(1-i)+1])
            xx = x*x if fluid else x*max(0,x-1)
            rates.extend((193*m/5,193*x/750000,1555*y,311*xx/(30000*m),
                          15*xx/m,15*x*y/m,100*x*y/m,y/1000,105*x/2,21*m/40000,
                          x/10,self.rho*m,self.epsilon*x*other_y/m,self.epsilon*y*other_y/m))
        return np.array([*rates,1.0])

    def admitted(self,state,word):
        return state[4] == self.birth and all(self.boxes[label].contains(state[2*i:2*i+2]) for i,label in enumerate(word))

    def run(self,word,seed=COPY_SEED,initial=None,limit=COMPUTE_EVENT_LIMIT,
            deadline=DIVISION_DEADLINE,quota=PHYSICAL_EVENT_QUOTA):
        if len(word) != 2 or any(label not in (0,1) for label in word) or limit < 1 or quota < 1 or deadline <= 0:
            raise ValueError('Two binary labels and positive limits required.')
        state = np.array((*PREPARATIONS[word[0]],*PREPARATIONS[word[1]],self.birth) if initial is None else initial,dtype=np.int64)
        if not self.admitted(state,word):
            raise ValueError('Initial parent must be admitted.')
        rng,time = np.random.default_rng(seed),0.0
        trace = [(time,*state.tolist())]
        status = 'unfinished'
        for step in range(limit):
            rates = self.rates(state)
            total = rates.sum()
            time_next = time+rng.exponential(1/total)
            if time_next > deadline:
                time=deadline;status='deadline';break
            time=time_next
            event=int(np.searchsorted(np.cumsum(rates),rng.random()*total,side='right'))
            state += self.jumps[event]
            assert np.min(state) >= 0
            if step+1 >= quota:
                status='event_quota';break
            if self.kinds[event] == 'reverse_growth':
                status='reverse_growth_monitor';break
            if any(state[2*i] > 16*state[4] or state[2*i+1] > 4*state[4] for i in range(2)):
                status='domain_exit';break
            if state[4] == 2*self.birth:
                status='division';break
            if step % 2000 == 0:
                trace.append((time,*state.tolist()))
        trace.append((time,*state.tolist()))
        result={'word':list(word),'status':status,'time':time,'events':step+1,
                'parent_at_stop':state.tolist(),'joint_partition_return':None,'both_return':None,'trace':trace}
        if status == 'division':
            payoff = math.prod(ComplementaryAllocation.joint_return(state[2*i:2*i+2],self.boxes[label]) for i,label in enumerate(word))
            A,B = ComplementaryAllocation.draw(state[:4],rng)
            result.update(joint_partition_return=str(payoff),
                          both_return=all(self.boxes[label].contains(d[2*i:2*i+2]) for i,label in enumerate(word) for d in (A,B)),
                          daughters=[[*A.tolist(),self.birth],[*B.tolist(),self.birth]])
        elif status != 'unfinished':
            result['both_return']=False
        return result

    @staticmethod
    def paper_budget():
        """Scalar assembly conditional on the paper's checked generator rows.

        This function does NOT regenerate/check the 22-million-row certificate.
        """
        premonitor=Q(9911,10000)-20*Q(2,10**6)-Q(1,50000)/Q(5,8)
        quota=Q(11,4)**80/Q(27,10)**100
        assert Q(2718,1000)**25 > 250000**2
        assert Q(101,100)**100 >= Q(27,10) and quota < Q(2,10**7)
        final=Q(991,1000)-80*Q(2,10**7)*53-Q(2,10**7)
        return {'premonitor':str(premonitor),'quota_upper_bound':str(quota),
                'assembled_lower_bound':str(final),'reported_lower_bound':'9901/10000',
                'generator_table_rechecked':False,'capacity_per_pool_molecules':4000*53+1001*10**8,
                'maximum_newborn_residents_per_module':874}


class OffspringKernel:
    """Finite joint daughter kernel; missing row mass is failed copying."""
    def __init__(self,states,rows):
        self.states,self.rows=tuple(states),rows
        if set(rows) != set(states):
            raise ValueError('One row per admitted restart state required.')
        for row in rows.values():
            if any(a not in states or b not in states or p < 0 for (a,b),p in row.items()) or sum(row.values()) > 1:
                raise ValueError('Invalid substochastic joint row.')

    def success(self,generations,family=False):
        if generations < 0:
            raise ValueError('Nonnegative horizon required.')
        values={s:Q(1) for s in self.states}
        for _ in range(generations):
            values={s:sum((p*values[a]*(values[b] if family else 1)
                           for (a,b),p in self.rows[s].items()),Q(0)) for s in self.states}
        return values

    @property
    def uniform_one_cycle(self):
        return min(sum(row.values()) for row in self.rows.values())


class RedundancyBounds:
    def __init__(self,k=MODULE_COUNT,gamma=GROWTH_RATE):
        if not isinstance(k,int) or k < 1 or not 0 < gamma <= Q(1,10**11):
            raise ValueError('The four-species bound requires k >= 1 and 0 < gamma <= 1e-11.')
        self.k,self.gamma=k,gamma

    def required(self,depth,error=TARGET_ERROR,family=False):
        if not isinstance(depth,int) or depth < 1 or not 0 < error < 1:
            raise ValueError('Positive horizon and error in (0,1) required.')
        with mp.workdps(80):
            divisions=(1 << depth)-1 if family else depth
            eta=mp.mpf(error.numerator)/error.denominator
            C=2703+mp.mpf(11)*self.gamma.denominator/(5*self.gamma.numerator)
            sufficient=max(14*10**19,int(mp.ceil(1024*10**18*mp.log(self.k*divisions*C/eta))))
            necessary=mp.log(self.k*divisions/(-mp.log1p(-eta)))/(280*mp.log(2))
            return {'depth':depth,'divisions':divisions,'sufficient_N':sufficient,
                    'necessary_real_lower_bound':mp.nstr(max(0,necessary),30),
                    'criterion':'immediate_return','family':family}


class FourSpeciesModule:
    """Reference chemistry for the redundancy theorem: A, B, z, H."""
    def __init__(self,e=1e-5,d0=1e-4):
        self.e,self.d0=e,d0
        # Consumption, production, reference coefficient; source/sink channels
        # represent the theorem's open chemistry, not the controlled pools above.
        pairs=(((0,),(1,2),1.,1.),((2,),(3,),16.,1.),
               ((2,2),(3,),2.,1.),((),(0,),6.,1.),((),(1,),27.,1.),
               ((1,),(0,0),e,e))
        self.reactions=tuple(r for a,b,k,l in pairs for r in ((a,b,k),(b,a,l)))+(((3,),(),d0),)

    def fluid(self,u):
        A,B,z,H=u;e,d0=self.e,self.d0
        return np.array((6-2*A+z*B+2*e*(B-A*A),27+A-(1+z)*B-e*(B-A*A),
                         A-B*z-16*z-4*z*z+3*H,16*z+2*z*z-(2+d0)*H))

    @staticmethod
    def stationary_reduction(z):
        e,d0=Q(1,100000),Q(1,10000)
        B=60/(z+2);H=(16*z+2*z*z)/(2+d0)
        K=(2*(1+2*d0)*z*z-16*(1-d0)*z)/(2+d0)
        A=z*B+K
        return (A,B,z,H),27-(1+e)*B+K+e*A*A


class CoupledModules:
    """Literal finite-count generator with shared consuming growth and exchange.

    Physical coefficients scale with k exactly as in this theorem. Building a
    generator does not assert that arbitrary parameters satisfy local quality.
    """
    def __init__(self,k,weights,gamma=float(GROWTH_RATE),module=None):
        weights=np.asarray(weights,dtype=float)
        if (k < 1 or weights.shape != (k,k) or np.any(weights < 0)
                or not np.array_equal(weights,weights.T) or np.any(np.diag(weights) != 0) or gamma <= 0):
            raise ValueError('Symmetric nonnegative exchange without self-edges required.')
        self.k,self.weights,self.gamma,self.module=k,weights,gamma,module or FourSpeciesModule()

    def rates_and_jumps(self,counts,marker):
        counts=np.asarray(counts,dtype=np.int64)
        if counts.shape != (self.k,4) or np.min(counts) < 0 or marker <= 0:
            raise ValueError('Nonnegative module counts and positive shared marker required.')
        rates,jumps=[],[]
        for i in range(self.k):
            for inputs,outputs,coefficient in self.module.reactions:
                order=len(inputs)
                rate=coefficient*self.k**(order-2)*float(marker)**(1-order)
                used={};jump=np.zeros(4*self.k+1,dtype=np.int64)
                for x in inputs:
                    rate*=max(0,int(counts[i,x])-used.get(x,0));used[x]=used.get(x,0)+1;jump[4*i+x]-=1
                for x in outputs:jump[4*i+x]+=1
                rates.append(rate);jumps.append(jump)
            jump=np.zeros(4*self.k+1,dtype=np.int64);jump[4*i+2]=-1;jump[-1]=1
            rates.append(self.gamma*counts[i,2]/self.k);jumps.append(jump)
            for j in range(self.k):
                if self.weights[i,j]:
                    jump=np.zeros(4*self.k+1,dtype=np.int64);jump[4*i+2]=-1;jump[4*j+2]=1
                    rates.append(self.weights[i,j]*counts[i,2]/self.k);jumps.append(jump)
        return np.array(rates),np.array(jumps)

    def interaction_within_paper_allowance(self):
        return bool(np.max(self.weights.sum(axis=1)) <= 1e-11 and self.gamma <= 1e-11)


@dataclass(frozen=True)
class ReturnEllipsoid:
    label: str
    center: tuple[Q, ...]
    metric: tuple[tuple[Q, ...], ...]
    eta: Q
    omega: int

    @classmethod
    def load(cls,path=Path(__file__).with_name('terminal_regions.json')):
        data=json.loads(path.read_text(encoding='utf-8'))
        return tuple(cls(r['label'],tuple(map(Q,r['center_M'])),
                         tuple(tuple(map(Q,row)) for row in r['terminal_metric']),Q(r['eta']),r['volume_scale'])
                     for r in data['regions'])

    def exact_count_energy(self,counts):
        if len(counts) != len(self.center):
            raise ValueError('Wrong region dimension.')
        y=tuple(Q(int(n),self.omega)-z for n,z in zip(counts,self.center))
        if len(y) != len(self.center):
            raise ValueError('Wrong region dimension.')
        return sum(y[i]*self.metric[i][j]*y[j] for i in range(8) for j in range(8))/self.eta

    def energy(self,concentration):
        y=np.asarray(concentration)-np.array(self.center,dtype=float)
        return float(y@np.array(self.metric,dtype=float)@y/float(self.eta))

    def cube_tolerance(self):
        return math.sqrt(float(self.eta/sum(abs(v) for row in self.metric for v in row)))


@dataclass(frozen=True)
class SemenovModel:
    exchange: float=EXCHANGE
    ligation: float=LIGATION
    hydrolysis: float=HYDROLYSIS
    inhibition: float=INHIBITION
    dilution: float=DILUTION
    feed: tuple[float,...]=FEED

    def __post_init__(self):
        if (min(self.exchange,self.ligation,self.hydrolysis,self.inhibition,self.dilution) <= 0
                or len(self.feed) != 8 or min(self.feed) < 0):
            raise ValueError('Positive kinetic coefficients and eight nonnegative feeds required.')

    def chemistry(self):
        # Input/output species indices; inactive products remain outside tracking.
        pairs=(((4,3),(1,5)),((4,2),(1,6)),((6,3),(2,5)))
        reactions=[]
        for left,right in pairs:
            reactions.extend(((left,right,self.exchange),(right,left,self.exchange)))
        reactions += [((0,1),(2,3),self.ligation),((0,),(3,),self.hydrolysis)]
        reactions += [((x,7),(),self.inhibition) for x in (1,2,3)]
        return tuple(reactions)

    def rhs(self,time,u):
        derivative=self.dilution*(np.array(self.feed)-u)
        for inputs,outputs,k in self.chemistry():
            rate=k*math.prod(u[i] for i in inputs)
            for i in inputs: derivative[i]-=rate
            for i in outputs: derivative[i]+=rate
        return derivative

    def recovery(self,region,rtol=2e-10,atol=2e-14):
        initial=(np.array(region.center,dtype=float)+np.array(self.feed))/2
        times=np.linspace(0,RECOVERY_TIME,251)
        answer=solve_ivp(self.rhs,(0,RECOVERY_TIME),initial,t_eval=times,method='BDF',rtol=rtol,atol=atol)
        if not answer.success: raise RuntimeError(answer.message)
        return times,answer.y.T

    def jump_rates(self,counts,omega):
        """Literal 11 chemical + 8 feed + 8 outflow rates and jumps.

        This reusable kernel is not a practicable SSA at the certified volume.
        """
        rates,jumps=[],[]
        for inputs,outputs,k in self.chemistry():
            used={};rate=k*float(omega)**(1-len(inputs));jump=[0]*8
            for i in inputs:
                rate*=max(0,int(counts[i])-used.get(i,0));used[i]=used.get(i,0)+1;jump[i]-=1
            for i in outputs:jump[i]+=1
            rates.append(rate);jumps.append(jump)
        for i in range(8):
            jump=[0]*8;jump[i]=1;rates.append(float(omega)*self.dilution*self.feed[i]);jumps.append(jump)
        for i in range(8):
            jump=[0]*8;jump[i]=-1;rates.append(self.dilution*int(counts[i]));jumps.append(jump)
        return np.array(rates),np.array(jumps)


def semenov_scalar_budget(label):
    """Exact endpoint budget, conditional on the supplied paper's tube bounds."""
    eta,L0,vs,r=((Q(6,10**12),1601,Q(115,1000),Q(9,100)),
                 (Q(9,10**14),3818,Q(114,1000),Q(26,100)))[label]
    omega,K=24*10**18,12*10**18
    m0=Q(omega)*Q(15231,100000)/2
    lam=Q(omega)*Q(15231,100000)*Q(2,1000)
    m1=L0*vs/(omega*eta)
    m2=Q(L0**2)/eta**2*(24*vs**2/omega**2+8*vs/omega**3)
    energy=r**6
    terms={'inherited_parent':2*energy-energy**2,'second_moment':(4+48*r)*m1,
           'fourth_moment':12*m2,'refill_inventory':64*m0/K**2,
           'recovery_and_inventory_drift':1000*(Q(1,10**9)+4*lam/K**2)}
    total=sum(terms.values())
    assert total <= Q(1,250)
    return {'terms':{k:str(v) for k,v in terms.items()},'joint_failure_upper_bound':str(2*total),
            'polynomial_tube_certificates_rechecked':False}


def write_csv(path,headers,rows):
    with path.open('w',newline='') as stream:
        writer=csv.writer(stream);writer.writerow(headers);writer.writerows(rows)


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=ConstructedSource()
    copies=[]
    for index,word in enumerate(itertools.product((0,1),repeat=2)):
        result=source.run(word,COPY_SEED+index)
        write_csv(out/f'copy_{word[0]}{word[1]}.csv',['time','X0','Y0','X1','Y1','marker'],result.pop('trace'))
        copies.append(result)
    regions=ReturnEllipsoid.load()
    if any(region.omega != OMEGA for region in regions):
        raise ValueError('Changing OMEGA requires compatible return-region data; nominal certification does not transfer automatically.')
    model=SemenovModel()
    from dataclasses import replace
    recoveries=[]
    trajectories=[]
    sensitivity=[]
    for label,region in enumerate(regions):
        times,trajectory=model.recovery(region)
        _,fine=model.recovery(region,rtol=2e-12,atol=2e-16)
        preparation=[int(region.omega*z) for z in region.center]
        recoveries.append({'label':region.label,'endpoint_energy':region.energy(fine[-1]),
            'thiol_endpoint_molar':float(sum(fine[-1,1:4])),
            'tolerance_endpoint_difference':float(np.max(np.abs(fine[-1]-trajectory[-1]))),
            'cube_tolerance_molar':region.cube_tolerance(),'integer_preparation':preparation,
            'exact_preparation_energy':str(region.exact_count_energy(preparation)),
            'scalar_budget':semenov_scalar_budget(label)})
        trajectories.append(fine)
        write_csv(out/f'recovery_{region.label}.csv',['time',*SPECIES],zip(times,*fine.T))
        for group in ('exchange','ligation','hydrolysis','inhibition','dilution'):
            for sign in (-1,1):
                changed=replace(model,**{group:getattr(model,group)*(1+sign*PERTURBATION)})
                _,run=changed.recovery(region)
                sensitivity.append({'label':region.label,'group':group,'relative_change':sign*PERTURBATION,
                    'terminal_energy_in_nominal_region':region.energy(run[-1]),
                    'decoded_high':bool(sum(run[-1,1:4]) > .005)})
    redundancy=[RedundancyBounds().required(G,family=family) for family in (False,True) for G in HORIZONS]
    kernel=OffspringKernel(('interior','edge'),{
        'interior':{('interior','interior'):Q(9,10),('edge','edge'):Q(9,100)},
        'edge':{('interior','edge'):Q(1,2),('edge','interior'):Q(12,25)}})
    kernel_result={'role':'illustrative state-dependent kernel, not fitted chemical data',
                  'uniform_one_cycle':str(kernel.uniform_one_cycle),
                  'lineage_5':{s:str(p) for s,p in kernel.success(5).items()},
                  'family_depth_5':{s:str(p) for s,p in kernel.success(5,True).items()}}
    root_brackets=[]
    for left,right in (('0.99579401232','0.99579401233'),('2.97636724376','2.97636724377')):
        _,fl=FourSpeciesModule.stationary_reduction(Q(left))
        _,fr=FourSpeciesModule.stationary_reduction(Q(right))
        assert fl*fr < 0
        root_brackets.append({'z_interval':[left,right],'residual_signs':[1 if fl > 0 else -1,1 if fr > 0 else -1]})
    result={'constructed_runs':copies,'constructed_scalar_budget':source.paper_budget(),
            'four_species_stationary_brackets':root_brackets,
            'nominal_recovery':recoveries,'sensitivity':sensitivity,'redundancy':redundancy,'joint_kernel_example':kernel_result}
    (out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
    lines=[f'Constructed word {r["word"]}: {r["status"]} at {r["time"]:.6g} s; both admitted={r["both_return"]}.' for r in copies]
    lines += [f'Nominal {r["label"]} recovery: terminal energy {r["endpoint_energy"]:.8g}; cube tolerance {r["cube_tolerance_molar"]*1e9:.6g} nM.' for r in recoveries]
    lines += ['Scalar probability budgets reproduced; the large generator/polynomial certificate tables are not rechecked.',
              'Four SSA paths and deterministic recovery curves do not estimate the certified uniform copying probability.']
    (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=(11,4.4),layout='constrained')
    for region,trajectory in zip(regions,trajectories):
        axes[0].plot(times,trajectory[:,1:4].sum(axis=1)*1000,label=region.label)
    axes[0].axhline(5,color='black',ls='--',label='Thiol decoder threshold')
    axes[0].set(xlabel='Recovery time (s)',ylabel='Free thiol (mM)',title='Nominal split, refill and recovery')
    for family,label in ((False,'Designated lineage'),(True,'Complete binary family')):
        rows=[r for r in redundancy if r['family']==family]
        axes[1].plot([r['depth'] for r in rows],[r['sufficient_N']/1e22 for r in rows],label=label)
    axes[1].set(xlabel='Lineage divisions / family depth',ylabel='Sufficient N / 10^22',title='Different model: four-species redundancy bound')
    for ax in axes:ax.legend(fontsize=8);ax.grid(alpha=.2)
    fig.savefig(out/'copying.png',dpi=180);fig.savefig(out/'copying.svg');plt.close(fig)
    fig,axes=plt.subplots(1,2,figsize=(11,4.4),layout='constrained')
    totals=range(241)
    joint=[float(ComplementaryAllocation.coordinate_return(n,6,106)) for n in totals]
    marginal=[sum(math.comb(n,d) for d in range(6,min(106,n)+1))/2**n for n in totals]
    axes[0].plot(totals,joint,label='Both complementary daughters admitted')
    axes[0].plot(totals,np.square(marginal),'--',label='Product of marginals (incorrect for siblings)')
    axes[0].set(xlabel='Parent X count at division',ylabel='Probability for the X coordinate',
                title='One coordinate: both daughters need 6 to 106')
    horizons=np.arange(1,11)
    axes[1].plot(horizons,.9901**horizons,label='Designated lineage: p^G')
    axes[1].plot(horizons,.9901**(2**horizons-1),label='Complete family: p^(2^G - 1)')
    axes[1].set(xlabel='Lineage divisions / family depth',ylabel='Guaranteed lower bound',
                yscale='log',title='Uniform one-cycle bound p = 0.9901')
    for ax in axes:ax.legend(fontsize=8);ax.grid(alpha=.2)
    fig.savefig(out/'allocation.png',dpi=180);fig.savefig(out/'allocation.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__)),'region_source_sha256':digest(Path(__file__).with_name('terminal_regions.json')),
        '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()