"""Resident chemistry, intact-cell transfer and two-cycle chemical-state selection."""
# EDITABLE INPUTS: theorem witness and separate, deliberately uncertified small diagnostics.
NEWBORN_SIZE = 65536 * 10**18
RETAINED_CELLS = 4 * 10**9
GROWTH_COUPLING = '1/100000000000'
SERVICE_FAILURE_ALLOWANCE = '1/1000'
DIAGNOSTIC_TRANSFER_TOLERANCE = '1/50'
DIAGNOSTIC_SIZE = 10
DIAGNOSTIC_CELLS = 4
DIAGNOSTIC_GROWTH = '2'             # Accelerated exploration; outside theorem scope.
DIAGNOSTIC_RECOVERY = '.05'        # The theorem instead requires 5376.
EVENT_BUDGET = 100000
RANDOM_SEED = 44092026

from dataclasses import dataclass,replace,asdict
from fractions import Fraction as Q
from itertools import combinations
from pathlib import Path
import argparse,csv,hashlib,json,math,platform
import numpy as np
import sympy as sp
from mpmath import mp,iv
from scipy.integrate import solve_ivp
from resident import Cell,ResidentChemistry,ChemicalRegions,Interval,MATRICES,A,reconstruction,residual

mp.dps=80;iv.dps=80
MANUSCRIPT_SHA256='d4d09105a0cfd39230861a7911081421a89f536afc2caf28d9297d00cab282ca'
ALPHA=Q(1,10**12);MIN_SIZE=14*10**19
SPECIES=('A','B','z','H');RESERVOIRS=('F','G','RA','RB','WH')
# Gross exchange for each directed resident label, including supplied and collected molecules.
EXCHANGE=(None,None,0,0,None,None,2,2,3,3,1,1,4)
def number(ctx,x):
    x=Q(str(x));return ctx.mpf(x.numerator)/x.denominator
def dot(a,b):return sum(x*y for x,y in zip(a,b))

class ChemicalModel(ResidentChemistry):
    def __init__(self):
        super().__init__();self.jumps=np.array([np.subtract(c.outputs,c.inputs) for c in self.channels]);self.coefficients=np.array([float(c.coefficient) for c in self.channels])
    def rates(self,cell):
        a,b,z,h=cell.counts;m=cell.size
        return np.array([a,b*z/m,16*z,h,h,2*z*(z-1)/m,6*m,a,27*m,b,b/100000,a*(a-1)/(100000*m),h/10000],float)
    def density(self,x):
        a,b,z,h=x
        return np.array([6-2*a+b*z+2e-5*(b-a*a),27+a-(1+z)*b-1e-5*(b-a*a),a-b*z-16*z-4*z*z+3*h,16*z+2*z*z-2.0001*h])
    def certificates(self):
        # Full species A,B,z,H,F,G,RA,RB,WH, with maintained reservoirs projected out.
        left=((0,),(2,4),(3,),(6,),(7,),(1,5));right=((1,2),(3,),(2,2),(0,),(1,),(0,0))
        element=(2,1,1,2,1,3,2,1,2);boltzmann=(Q(1),Q(1),Q(1),Q(2),Q(1,8),Q(1),Q(1,6),Q(1,27),Q(1))
        for j,(l,r) in enumerate(zip(left,right)):
            if sum(element[i] for i in l)!=sum(element[i] for i in r):raise ArithmeticError('Full material balance failed.')
            ratio=math.prod(boltzmann[i] for i in r)/math.prod(boltzmann[i] for i in l)
            if ratio!=self.channels[2*j].coefficient/self.channels[2*j+1].coefficient:raise ArithmeticError('Common chemical potentials failed.')
            for channel,reactants,products in [(self.channels[2*j],l,r),(self.channels[2*j+1],r,l)]:
                if channel.inputs!=tuple(reactants.count(i) for i in range(4)) or channel.outputs!=tuple(products.count(i) for i in range(4)):raise ArithmeticError('Reservoir projection failed.')
        if element[3]!=element[8]:raise ArithmeticError('Sink material imbalance.')
        matrices={}
        for tag,P in MATRICES.items():
            P=sp.Matrix(P)/10**6
            if any((P-sp.eye(4)/200)[:k,:k].det()<=0 or (42*sp.eye(4)-P)[:k,:k].det()<=0 for k in range(1,5)):raise ArithmeticError('Energy sandwich failed.')
            matrices[tag]=P.tolist()
        z=sp.Symbol('z');coords=reconstruction(z);a,b,zz,h=coords;eps=sp.Rational(1,100000)
        f=(6-2*a+b*zz+2*eps*(b-a*a),27+a-(1+zz)*b-eps*(b-a*a),a-b*zz-16*zz-4*zz*zz+3*h,16*zz+2*zz*zz-sp.Rational(20001,10000)*h)
        r=residual(z)
        if any(sp.simplify(x)!=0 for x in (f[0]+2*r,f[1]-r,f[2],f[3])):raise ArithmeticError('Stationary reconstruction failed.')
        return dict(element=element,boltzmann_weights=list(map(str,boltzmann)),energy_matrices=matrices,stationary_curve_identity=True,scope='Exact chemical, source projection and energy-matrix checks. The local stochastic exponential-generator inequality is imported from the paper.')
    def recovery_density(self,cell,duration=5376,method='Radau'):
        initial=np.array(cell.counts,float)/cell.size
        sol=solve_ivp(lambda t,x:self.density(x),(0,duration),initial,method=method,rtol=2e-12,atol=1e-14,dense_output=True)
        if not sol.success:raise ArithmeticError(sol.message)
        return sol

@dataclass(frozen=True)
class Population:
    cells:tuple
    precursor:int
    omega:int
    endpoint:int
    divisions:int=0
    @property
    def size(self):return sum(c.size for c in self.cells)

class ComplementaryDivision:
    def split(self,counts,rng):
        if max(counts)>=2**53:raise ValueError('Numerical division requires molecule counts below 2^53.')
        first=tuple(int(rng.binomial(n,.5)) for n in counts)
        return first,tuple(n-a for n,a in zip(counts,first))

class PopulationSource:
    """Physical source only: no type-dependent rates, reset, or analytical failure censoring."""
    def __init__(self,N,gamma,chemistry=None,division=None):
        self.N=N;self.gamma=Q(str(gamma));self.chemistry=chemistry or ChemicalModel();self.division=division or ComplementaryDivision()
        if type(N) is not int or N<1 or self.gamma<=0:raise ValueError('Positive integer N and growth coefficient required.')
    def refill(self,cells):
        cells=tuple(cells)
        if not cells or any(not self.N<=c.size<2*self.N for c in cells):raise ValueError('Cells must lie in the live-size interval.')
        W=sum(c.size for c in cells);return Population(cells,4*W,4*W,W)
    def growth_rate(self,pop,cell):return self.gamma*Q(pop.precursor,pop.omega)*cell.counts[2] if pop.precursor else Q(0)
    def step(self,pop,index,label,rng,allocation=None):
        if type(index) is not int or not 0<=index<len(pop.cells) or type(label) is not int or not 0<=label<=13:raise ValueError('Unknown cell or event label.')
        cell=pop.cells[index];Qnext=pop.precursor;divisions=pop.divisions
        if label<13:children=(self.chemistry.channels[label].apply(cell),)
        elif label==13:
            if self.growth_rate(pop,cell)<=0:raise ValueError('Disabled growth event.')
            n=list(cell.counts);n[2]-=1;grown=Cell(cell.tag,cell.size+1,tuple(n));Qnext-=1
            if grown.size==2*self.N:
                if allocation is None:first,second=self.division.split(grown.counts,rng)
                else:
                    first=tuple(allocation);second=tuple(n-a for n,a in zip(grown.counts,first))
                    if len(first)!=4 or min(first+second)<0:raise ValueError('Invalid complementary allocation.')
                children=(Cell(cell.tag,self.N,first),Cell(cell.tag,self.N,second));divisions+=1
            elif grown.size<2*self.N:children=(grown,)
            else:raise ValueError('Division threshold exceeded.')
        else:raise ValueError('Unknown event label.')
        result=replace(pop,cells=pop.cells[:index]+children+pop.cells[index+1:],precursor=Qnext,divisions=divisions)
        if result.size+result.precursor!=pop.size+pop.precursor:raise ArithmeticError('Precursor/size conservation failed.')
        return result
    def run(self,pop,duration,rng,event_budget,collect_endpoint=True,quota=None,hard_limit=None):
        if duration<0 or type(event_budget) is not int or event_budget<0:raise ValueError('Nonnegative duration/event budget required.')
        if quota is not None and (type(quota) is not int or quota<1):raise ValueError('Positive integer service quota required.')
        if hard_limit is not None and (type(hard_limit) is not int or hard_limit<0):raise ValueError('Nonnegative integer hard budget required.')
        time=0.;gross=[0]*5;steps=0;status='event_budget';counter=0
        if duration==0:return pop,dict(status='duration',time=0.,events=0,gross_exchange=gross,service_counter=0,saturated=False,complete=True)
        for _ in range(event_budget):
            if collect_endpoint and pop.precursor==pop.endpoint:status='endpoint';break
            arrays=[np.r_[self.chemistry.rates(c),float(self.growth_rate(pop,c))] for c in pop.cells];rates=np.concatenate(arrays);total=float(rates.sum())
            proposed=time+float(rng.exponential(1/total))
            if proposed>duration:time=float(duration);status='duration';break
            if proposed<=time:raise ArithmeticError('Time resolution exhausted; use a smaller diagnostic.')
            chosen=min(int(np.searchsorted(np.cumsum(rates),rng.random()*total,side='right')),len(rates)-1);index,label=divmod(chosen,14)
            reservoir=EXCHANGE[label] if label<13 else None
            if hard_limit is not None and reservoir is not None and gross[reservoir]>=hard_limit:status='hard_service_cutoff';break
            pop=self.step(pop,index,label,rng);steps+=1;time=proposed
            if reservoir is not None:gross[reservoir]+=1
            counter=counter+1 if quota is None else min(quota,counter+1)
            if collect_endpoint and pop.precursor==pop.endpoint:status='endpoint';break
        # Saturating the counter never changes a physical transition or terminates the source.
        return pop,dict(status=status,time=time,events=steps,gross_exchange=gross,service_counter=counter,saturated=quota is not None and counter==quota,complete=status in ('endpoint','duration'))

class UniformTransfer:
    def sample(self,cells,M,rng):
        cells=tuple(cells)
        if type(M) is not int or not 1<=M<=len(cells):raise ValueError('Sample size must fit the endpoint population.')
        chosen=set(map(int,rng.choice(len(cells),M,replace=False)))
        return tuple(c for i,c in enumerate(cells) if i in chosen),tuple(c for i,c in enumerate(cells) if i not in chosen)
    def enumerate(self,cells,M,tolerance):
        cells=tuple(cells);eps=Q(tolerance);n=len(cells)
        if not 0<eps<1 or not 1<=M<=n or math.comb(n,M)>100000:raise ValueError('Invalid transfer or exact enumeration exceeds 100,000 subsets.')
        means={tag:Q(M,n)*sum(c.size for c in cells if c.tag==tag) for tag in ('H','L')};values=[];bad=0
        for indices in combinations(range(n),M):
            sizes={tag:sum(cells[i].size for i in indices if cells[i].tag==tag) for tag in ('H','L')}
            failure=any(abs(sizes[tag]-means[tag])>eps*means[tag] for tag in means);bad+=failure;values.append(sizes)
        moments={}
        for tag in means:
            mean=sum(Q(row[tag]) for row in values)/len(values);variance=sum((row[tag]-mean)**2 for row in values)/len(values)
            moments[tag]=dict(mean=str(mean),variance=str(variance))
        return dict(subsets=len(values),bad=bad,failure=str(Q(bad,len(values))),means={k:str(v) for k,v in means.items()},moments=moments)
    @staticmethod
    def weighted_moments(weights,M):
        weights=tuple(map(Q,weights));n=len(weights)
        if not 1<=M<=n or any(not 0<=w<=1 for w in weights):raise ValueError('Weights in [0,1] and valid sample size required.')
        p=Q(M,n);joint=Q(M*(M-1),n*(n-1)) if n>1 else Q(0)
        mean=p*sum(weights);variance=p*(1-p)*sum(w*w for w in weights)+2*(joint-p*p)*sum(weights[i]*weights[j] for i in range(n) for j in range(i+1,n))
        if not 0<=variance<=mean:raise ArithmeticError('Without-replacement variance bound failed.')
        return dict(mean=str(mean),variance=str(variance),pair_covariance=str(joint-p*p) if n>1 else None)

class SerialProtocol:
    def __init__(self,source,M,recovery_duration=5376):
        if type(M) is not int or M<1 or recovery_duration<0:raise ValueError('Positive retained count and nonnegative recovery duration required.')
        self.source=source;self.M=M;self.recovery_duration=recovery_duration;self.transfer=UniformTransfer()
    def run(self,cells,cycles,seed,event_budget):
        if type(cycles) is not int or cycles<0 or type(event_budget) is not int or event_budget<0:raise ValueError('Nonnegative integer cycle count and event budget required.')
        rng=np.random.default_rng(seed);pop=self.source.refill(cells);history=[];used=0
        for cycle in range(cycles):
            W0=pop.size;pop,batch=self.source.run(pop,float(8/self.source.gamma),rng,event_budget-used);used+=batch['events']
            if batch['status']!='endpoint':return dict(complete=False,stage='batch',history=history,partial=batch,population=asdict(pop),events=used)
            if len(pop.cells)<self.M:return dict(complete=False,stage='insufficient_cells',history=history,events=used)
            selected,discarded=self.transfer.sample(pop.cells,self.M,rng);before=selected;residual=pop.precursor
            recovering=replace(pop,cells=selected,precursor=0)
            recovered,recovery=self.source.run(recovering,self.recovery_duration,rng,event_budget-used,collect_endpoint=False);used+=recovery['events']
            record=dict(cycle=cycle,batch=batch,recovery=recovery,start_size=W0,consumed_precursor=3*W0,residual_discarded=residual,discarded_size=sum(c.size for c in discarded),selected_size=sum(c.size for c in selected),selected=[asdict(c) for c in before],recovered=[asdict(c) for c in recovered.cells])
            history.append(record)
            if recovery['status']!='duration':return dict(complete=False,stage='recovery',history=history,population=asdict(recovered),events=used)
            pop=self.source.refill(recovered.cells)
        return dict(complete=True,history=history,population=asdict(pop),events=used,scope='Physical diagnostic completion, not a theorem-success classification; analytical energy and odds marks are not applied.')

class TwoCycleCertificate:
    def __init__(self,N=NEWBORN_SIZE,M=RETAINED_CELLS,gamma=GROWTH_COUPLING,service=SERVICE_FAILURE_ALLOWANCE):
        self.N=N;self.M=M;self.gamma=Q(str(gamma));self.service=Q(str(service));self.u=Q(N)*ALPHA*A;self.T=8/self.gamma
        if type(N) is not int or N<MIN_SIZE or type(M) is not int or M<2 or M%2 or not 0<self.gamma<=Q(1,10**11) or not 0<self.service<1:raise ValueError('Requires admitted N, even M>=2, 0<gamma<=1e-11 and service allowance in (0,1).')
        self.qR=34000*N;self.qB=280000*N*M;self.JB=math.ceil(self.T*self.qB/self.service)+1;self.JR=math.ceil(M*5376*self.qR/self.service)+1
    def logs(self,p,ctx=mp):
        p=Q(str(p))
        if not 0<p<=1:raise ValueError('Positive count floor <=1 required.')
        N=ctx.mpf(self.N);M=ctx.mpf(self.M);u=number(ctx,self.u);T=number(ctx,self.T);ln=ctx.log
        return dict(outer_initial=ln(M)-7*u,daughter_initial=ln(14*M)-4*u,outer_drift=ln(M*T*u/42)-15*u/2,division_initial=ln(M)-u,division_drift=ln(M*T*u/42)-3*u/2,partition=ln(56*M)-N/(35*10**12),deadline=-N/2500,size_odds=-19*N/500000,
            batch_service=ln(T*self.qB/self.JB),transfer=ln(number(ctx,32/(Q(1,50)**2*p*self.M))),recovery_initial=ln(M)-u,recovery_return=ln(2*M)-u/2,recovery_exit=ln(M)-8*u,recovery_exit_drift=ln(16*M*u)-31*u/2,recovery_service=ln(M*5376*self.qR/self.JR))
    def evaluate(self):
        terms=[self.logs(p) for p in ('1/2','1/100')];intervals=[self.logs(p,iv) for p in ('1/2','1/100')]
        error=sum(mp.exp(x) for d in terms for x in d.values());error_interval=sum(iv.exp(x) for d in intervals for x in d.values());lower=1-error_interval
        digits=max(0,int(mp.floor(mp.mpf(lower.a)*10**9)))
        while digits and not bool(number(iv,Q(digits,10**9)).b<=lower.a):digits-=1
        return dict(N=self.N,M=self.M,u=str(self.u),success_lower_ordinary=mp.nstr(max(0,1-error),30),success_lower_certified=f'{digits/10**9:.9f}',error_interval=str(error_interval),cycle_component_log10=[{k:mp.nstr(v/mp.log(10),25) for k,v in d.items()} for d in terms],quotas=dict(qB=self.qB,qR=self.qR,JB=self.JB,JR=self.JR),fresh_precursor_bound=12*self.N*self.M,gross_each_reservoir_bound=2*(self.JB+self.M*self.JR),elapsed_model_time=str(2*(self.T+5376)))
    def rational_witness(self):
        if self.gamma!=Q(1,10**11) or self.u not in (128,256):raise ValueError('Paper rational envelopes implemented for u=128 or 256 at gamma=1e-11.')
        env=Q(1,10**38);half=Q(1,2**64) if self.u==128 else env;N=self.N;M=self.M;u=self.u;T=self.T
        if min(7*u,4*u,15*u/2,u,3*u/2,Q(N,35*10**12),Q(N,2500),Q(19*N,500000),8*u,31*u/2)<128:raise ArithmeticError('Envelope scope failed.')
        chemical=(M+14*M+M*T*u/42+M+M*T*u/42+56*M+2)*env
        recovery=M*((2+16*u)*env+2*half);transfer=Q(8000000,M)
        bound=2*(chemical+recovery+transfer+2*self.service)
        return dict(nonservice_uniform_bound=str(chemical+recovery+transfer),two_cycle_failure_bound=str(bound),success_lower=str(max(0,1-bound)),scope='Exact rational envelope uses e>=2, e^-128<=1e-38 and e^-64<=2^-64, with uniform p=1/100 in both cycles.')

class OddsAccounting:
    @staticmethod
    def coordinates(cells):
        counts={t:sum(c.tag==t for c in cells) for t in ('H','L')};sizes={t:sum(c.size for c in cells if c.tag==t) for t in ('H','L')}
        if not all(counts.values()):return dict(counts=counts,sizes=sizes,log_count=None,log_size=None,phase=None)
        C=mp.log(mp.mpf(counts['H'])/counts['L']);S=mp.log(mp.mpf(sizes['H'])/sizes['L']);phase=mp.log(mp.mpf(sizes['H'])*counts['L']/(sizes['L']*counts['H']))
        return dict(counts=counts,sizes=sizes,log_count=mp.nstr(C,30),log_size=mp.nstr(S,30),phase=mp.nstr(phase,30))
    @staticmethod
    def theorem_thresholds():
        def values(ctx):
            g=number(ctx,'3/5')*ctx.log(4)-number(ctx,'19/500')-ctx.log(number(ctx,'51/49'));one=g-ctx.log(2);two=2*g-ctx.log(2)
            return g,one,two,1/(1+ctx.exp(-two))
        g,one,two,f=values(mp);enclosures=values(iv)
        if not bool(enclosures[3].a>number(iv,'0.6930453').b):raise ArithmeticError('Decimal fraction threshold unresolved.')
        return dict(cycle_size_gain=mp.nstr(g,30),one_count_gain=mp.nstr(one,30),two_count_gain=mp.nstr(two,30),two_high_fraction=mp.nstr(f,30),two_fraction_certified_lower='0.6930453',two_count_interval=str(enclosures[2]),two_fraction_interval=str(enclosures[3]),eligibility_floors=['1/2','49/1600','2401/1280000'],exact_two_gain='(7/5)log(2)-19/250-2log(51/49)',scope='Two cycles from balanced ready newborns only. Intermediate division phases telescope; no separate log(2) penalty per cycle.')

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=ChemicalModel();regions=ChemicalRegions();certificate=TwoCycleCertificate();founders={t:regions.newborn(t,NEWBORN_SIZE) for t in ('L','H')}
    transfer=UniformTransfer();rows=[]
    for name,sizes in [('equal',[10]*8),('high_larger',[19]*4+[10]*4),('low_larger',[10]*4+[19]*4),('mixed',[19,10,19,10,10,19,10,19])]:
        cells=tuple(Cell('H' if i<4 else 'L',size,(20,20,30 if i<4 else 10,40)) for i,size in enumerate(sizes))
        for eps in (DIAGNOSTIC_TRANSFER_TOLERANCE,'1/2'):rows.append(dict(endpoint=name,epsilon=eps,**transfer.enumerate(cells,4,eps)))
    # Actual intact selected states, at large size, are followed by a density recovery illustration.
    endpoint=tuple(replace(founders['H' if i<4 else 'L'],counts=tuple(n+(NEWBORN_SIZE//50000 if j==2 else 0) for j,n in enumerate(founders['H' if i<4 else 'L'].counts))) for i in range(8))
    selected,discarded=transfer.sample(endpoint,4,np.random.default_rng(RANDOM_SEED));recovery=[];difference=0.
    for index,cell in enumerate(selected):
        sol=model.recovery_density(cell);other=model.recovery_density(cell,method='BDF');times=np.unique(np.r_[0,np.geomspace(.001,5376,180)])
        values=sol.sol(times).T;difference=max(difference,float(np.max(abs(values-other.sol(times).T))))
        recovery.extend((index,cell.tag,float(t),*x) for t,x in zip(times,values))
    result=dict(chemistry=model.certificates(),founders={t:dict(cell=asdict(c),energy=regions.energy(c).json(),readout=c.readout,multiplicity=RETAINED_CELLS//2) for t,c in founders.items()},configured=certificate.evaluate(),paper_witnesses=[dict(evaluated=c.evaluate(),rational=c.rational_witness()) for c in (TwoCycleCertificate(65536*10**18,4*10**9,'1/100000000000','1/1000'),TwoCycleCertificate(131072*10**18,10**15,'1/100000000000','1/1000'))],thresholds=OddsAccounting.theorem_thresholds(),exact_transfer=rows,
        selected_recovery=dict(selected=[asdict(c) for c in selected],discarded=[asdict(c) for c in discarded],selected_energy_enclosures=[regions.energy(c).json() for c in selected],solver_difference=difference,scope='Constructed admitted endpoint population, then intact sampling and full resident density ODE on the actual selected count/size inputs; fixed size and no precursor. Not a batch simulation or the finite-count recovery probability.'),scope='Bound evaluation imports the paper probability theorems; no Lean replay. Material, rate ratios, stationary identities, energy matrices and tiny transfer laws are freshly checked.')
    if args.simulate:
        cells=[]
        for i in range(DIAGNOSTIC_CELLS):
            tag='H' if i<DIAGNOSTIC_CELLS//2 else 'L';coords=regions.centers[tag];cells.append(Cell(tag,DIAGNOSTIC_SIZE,tuple(math.floor(DIAGNOSTIC_SIZE*(v.lo+v.hi)/2) for v in coords)))
        result['uncertified_small_mission']=SerialProtocol(PopulationSource(DIAGNOSTIC_SIZE,DIAGNOSTIC_GROWTH),DIAGNOSTIC_CELLS,float(Q(DIAGNOSTIC_RECOVERY))).run(cells,2,RANDOM_SEED,EVENT_BUDGET)
    def dump(name,data):(out/name).write_text(json.dumps(data,indent=2,default=str)+'\n')
    def table(name,header,rows):
        with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    dump('results.json',result);table('recovery.csv',['selected_index','ancestry','time',*SPECIES],recovery)
    table('transfer.csv',['endpoint','epsilon','subsets','bad','failure'],[(r['endpoint'],r['epsilon'],r['subsets'],r['bad'],r['failure']) for r in rows])
    sweep=[]
    for N in [int(x)*10**18 for x in np.linspace(140,70000,71)]:
        c=TwoCycleCertificate(N,RETAINED_CELLS);d=c.logs('1/100');groups={k:mp.log(sum(mp.exp(v) for name,v in d.items() if name in names))/mp.log(10) for k,names in {'chemical':('outer_initial','daughter_initial','outer_drift','division_initial','division_drift','partition','deadline','size_odds'),'recovery':('recovery_initial','recovery_return','recovery_exit','recovery_exit_drift'),'transfer':('transfer',)}.items()};sweep.append((N,*map(float,groups.values())))
    table('error_sweep.csv',['N','chemical_log10','recovery_log10','transfer_log10'],sweep)
    lines=[f'Thirteen resident channels, balanced completion, stationary curves and exact energy sandwich checked.',f'Configured two-cycle success >= {result["configured"]["success_lower_certified"]}; measured high fraction strictly above {result["thresholds"]["two_fraction_certified_lower"]} on success.',f'Exact eight-cell transfer failures: {[r["failure"] for r in rows]}.',f'Actual selected states retained for precursor-free density recovery; independent solver difference {difference:.3g}.','The certificate covers two cycles from ready newborn founders, not indefinite selection. Small accelerated simulations are separately labeled and outside its hypotheses.']
    (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(sweep)
    for j,label in enumerate(['Batch chemistry','Recovery','Transfer']):axs[0].plot(data[:,0]/1e22,data[:,j+1],label=label)
    axs[0].set(xlabel='Newborn size N / 10^22',ylabel='log10 per-cycle error upper bound',ylim=(-50,20),title='Batch, transfer and recovery error bounds');axs[0].legend(fontsize=8)
    x=np.arange(4);axs[1].bar(x-.18,[float(Q(rows[2*i]['failure'])) for i in range(4)],.36,label=f'Tolerance {DIAGNOSTIC_TRANSFER_TOLERANCE}');axs[1].bar(x+.18,[float(Q(rows[2*i+1]['failure'])) for i in range(4)],.36,label='Tolerance 1/2');axs[1].set_xticks(x,['Equal','High larger','Low larger','Mixed']);axs[1].set(ylabel='Exact transfer failure probability',title='Four of eight intact cells, all 70 subsets');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(axis='y',alpha=.2)
    fig.savefig(out/'selection.png',dpi=180);fig.savefig(out/'selection.svg');plt.close(fig)
    fig,ax=plt.subplots(figsize=(8,4.2),layout='constrained')
    for tag,color in [('L','#417d8c'),('H','#bd5a24')]:
        for i,c in enumerate(selected):
            if c.tag!=tag:continue
            data=np.array([r[2:] for r in recovery if r[0]==i],float);center=float((regions.centers[tag][2].lo+regions.centers[tag][2].hi)/2)
            ax.semilogx(np.maximum(data[:,0],1e-3),data[:,3]-center,color=color,label=f'Selected {tag}, index {i}')
    ax.set(xlabel='Precursor-free recovery time',ylabel='z concentration minus stationary z',title='Chemical recovery after intact-cell\ntransfer');ax.legend(fontsize=8);ax.grid(alpha=.2);ax.ticklabel_format(axis='y',style='sci',scilimits=(0,0))
    fig.savefig(out/'recovery.png',dpi=180);fig.savefig(out/'recovery.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,source_sha256=digest(Path(__file__)),resident_sha256=digest(Path(__file__).with_name('resident.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()