"""Finite-batch enrichment under one shared precursor pool.

Default: evaluate the manuscript bounds and inspect a literal count-model kernel.
Optional --simulate runs a budget-limited direct jump simulation, not a proof.
"""
from fractions import Fraction as Q
# EDITABLE INPUTS: manuscript units, not laboratory-calibrated concentrations.
NEWBORN_SIZE = 111 * 10**20
HIGH_FOUNDERS = 1
LOW_FOUNDERS = 1
GROWTH_COUPLING = Q(1,10**11)
LOG_ODDS_GAIN = Q(1,10)
SCALED_DEADLINE = Q(8)  # gamma*T; physical T is this divided by gamma
TARGET_ERROR = Q(1,100)
COPY_NUMBER_SWEEP = (14*10**19,111*10**20,308*10**20)
RANDOM_SEED = 2101
EVENT_BUDGET = 1000
MANUSCRIPT_SHA256 = 'b93345aac9b134b314f9e6115f46585fd3e5fc60c2ece42f5da6eacab32199ed'

import argparse
import csv
import hashlib
import json
import math
import platform
from dataclasses import dataclass, replace
from pathlib import Path
import mpmath as mp
import numpy as np
import sympy as sp

A = Q(1,512000000)
B = 16*A
ALPHA = Q(1,10**12)
LAMBDA = Q(1,10**15)
DELTA = Q(1,10**6)
C = Q(1,1024*10**18)
THETA = Q(1,1000)
N_MIN = 14*10**19
GAMMA_MAX = Q(1,10**11)
SPECIES = ('A','B','z','H')
MATRICES = {
    'L': ((1115801,-613708,2346047,3444885),(-613708,1111499,-2339313,-3434384),
          (2346047,-2339313,6767757,10178806),(3444885,-3434384,10178806,15517433)),
    'H': ((846084,-1159490,2352853,3510058),(-1159490,4333988,-6781643,-10233028),
          (2352853,-6781643,11398763,17222276),(3510058,-10233028,17222276,26082111)),
}


@dataclass(frozen=True)
class Interval:
    lo: Q
    hi: Q | None = None
    def __post_init__(self):
        object.__setattr__(self,'lo',Q(self.lo))
        object.__setattr__(self,'hi',Q(self.lo if self.hi is None else self.hi))
        if self.lo>self.hi:raise ValueError('Inverted interval')
    @staticmethod
    def cast(v):return v if isinstance(v,Interval) else Interval(v)
    def __add__(self,v):v=self.cast(v);return Interval(self.lo+v.lo,self.hi+v.hi)
    __radd__=__add__
    def __neg__(self):return Interval(-self.hi,-self.lo)
    def __sub__(self,v):return self+-self.cast(v)
    def __rsub__(self,v):return self.cast(v)+-self
    def __mul__(self,v):
        v=self.cast(v);values=[x*y for x in (self.lo,self.hi) for y in (v.lo,v.hi)]
        return Interval(min(values),max(values))
    __rmul__=__mul__
    def __truediv__(self,v):
        v=self.cast(v)
        if v.lo<=0<=v.hi:raise ZeroDivisionError('Interval denominator includes zero')
        return self*Interval(1/v.hi,1/v.lo)
    def __rtruediv__(self,v):return self.cast(v)/self
    def json(self):return [str(self.lo),str(self.hi)]


def reconstruction(z):
    b=60/(z+2);k=(20004*z*z-159984*z)/20001
    return (z*b+k,b,z,(16*z+2*z*z)/Q(20001,10000))


def residual(z):
    aa,bb,_,_=reconstruction(z)
    return Q(1,100000)*aa*aa+aa+Q(99999,100000)*bb-33


class ChemicalRegions:
    """Rational energy enclosures around tightly isolated stationary roots."""
    def __init__(self,refinements=150):
        self.roots={};self.centers={}
        for tag,left,right in [('L',Q('0.99579401232'),Q('0.99579401233')),
                               ('H',Q('2.97636724376'),Q('2.97636724377'))]:
            fl=residual(left);assert fl*residual(right)<0
            for _ in range(refinements):
                mid=(left+right)/2;fm=residual(mid)
                if fm==0:left=right=mid;break
                if fl*fm<0:right=mid
                else:left=mid;fl=fm
            self.roots[tag]=Interval(left,right)
            self.centers[tag]=reconstruction(self.roots[tag])

    def energy(self,cell):
        delta=[Q(n,cell.size)-center for n,center in zip(cell.counts,self.centers[cell.tag])]
        value=sum((Q(MATRICES[cell.tag][i][j],10**6)*delta[i]*delta[j]
                   for i in range(4) for j in range(4)),Interval(0))
        return Interval(max(Q(0),value.lo),max(Q(0),value.hi))

    def admits(self,cell,threshold,closed=False):
        energy=self.energy(cell)
        if (energy.hi<=threshold if closed else energy.hi<threshold):return True
        if (energy.lo>threshold if closed else energy.lo>=threshold):return False
        raise ArithmeticError('Energy comparison unresolved; refine the stationary-root brackets.')

    def newborn(self,tag,N):
        centers=self.centers[tag]
        def nearest(v):return (2*v.numerator+v.denominator)//(2*v.denominator)
        cell=Cell(tag,N,tuple(nearest(N*(v.lo+v.hi)/2) for v in centers))
        if not self.admits(cell,4*A,closed=True):
            raise ValueError('Rounded preparation is outside the closed newborn region at this N.')
        return cell


@dataclass(frozen=True)
class Cell:
    tag: str
    size: int
    counts: tuple[int,int,int,int]
    def __post_init__(self):
        if self.tag not in ('H','L') or type(self.size) is not int or self.size<1:
            raise ValueError('A cell needs a valid tag and positive integer size.')
        if len(self.counts)!=4 or any(type(n) is not int or n<0 for n in self.counts):
            raise ValueError('Resident counts must be four nonnegative Python integers.')
    @property
    def readout(self):
        value=10*self.counts[2]-self.counts[1]
        return 'H' if value>0 else 'L' if value<0 else 'ambiguous'


@dataclass(frozen=True)
class Channel:
    name: str
    inputs: tuple[int,...]
    outputs: tuple[int,...]
    coefficient: Q
    def rate(self,cell):
        value=self.coefficient*Q(cell.size)**(1-sum(self.inputs))
        for count,order in zip(cell.counts,self.inputs):
            if count<order:return Q(0)
            for offset in range(order):value*=count-offset
        return value
    def apply(self,cell):
        if self.rate(cell)<=0:raise ValueError('A zero-propensity channel cannot fire.')
        return replace(cell,counts=tuple(n-out+new for n,out,new in zip(cell.counts,self.inputs,self.outputs)))


class ResidentChemistry:
    def __init__(self):
        self.channels=(
            Channel('A to B+z',(1,0,0,0),(0,1,1,0),Q(1)),
            Channel('B+z to A',(0,1,1,0),(1,0,0,0),Q(1)),
            Channel('z to H',(0,0,1,0),(0,0,0,1),Q(16)),
            Channel('H to z',(0,0,0,1),(0,0,1,0),Q(1)),
            Channel('H to 2z',(0,0,0,1),(0,0,2,0),Q(1)),
            Channel('2z to H',(0,0,2,0),(0,0,0,1),Q(2)),
            Channel('feed A',(0,0,0,0),(1,0,0,0),Q(6)),
            Channel('remove A',(1,0,0,0),(0,0,0,0),Q(1)),
            Channel('feed B',(0,0,0,0),(0,1,0,0),Q(27)),
            Channel('remove B',(0,1,0,0),(0,0,0,0),Q(1)),
            Channel('B to 2A',(0,1,0,0),(2,0,0,0),Q(1,100000)),
            Channel('2A to B',(2,0,0,0),(0,1,0,0),Q(1,100000)),
            Channel('remove H',(0,0,0,1),(0,0,0,0),Q(1,10000)),
        )

    def generator(self,cell,observable):
        base=observable(cell)
        return sum((ch.rate(cell)*(observable(ch.apply(cell))-base)
                    for ch in self.channels if ch.rate(cell)),Q(0))


@dataclass(frozen=True)
class BatchParameters:
    N: int = NEWBORN_SIZE
    h: int = HIGH_FOUNDERS
    low: int = LOW_FOUNDERS
    gamma: Q = GROWTH_COUPLING
    gain: Q = LOG_ODDS_GAIN
    scaled_deadline: Q = SCALED_DEADLINE
    def __post_init__(self):
        if any(type(v) is not int or v<1 for v in (self.N,self.h,self.low)):
            raise ValueError('N and both founder numbers must be positive integers.')
        for key in ('gamma','gain','scaled_deadline'):
            object.__setattr__(self,key,Q(str(getattr(self,key))))
        if self.gamma<=0 or self.scaled_deadline<0:raise ValueError('Invalid growth rate or deadline.')
    @property
    def M(self):return self.h+self.low
    @property
    def horizon(self):return self.scaled_deadline/self.gamma
    @property
    def theorem_scope(self):return self.N>=N_MIN and self.gamma<=GAMMA_MAX


@dataclass(frozen=True)
class Population:
    precursor: int
    cells: tuple[Cell,...]
    divisions: int = 0
    status: str = 'active'
    @property
    def total_size(self):return sum(c.size for c in self.cells)
    def tag_size(self,tag):return sum(c.size for c in self.cells if c.tag==tag)
    def tag_count(self,tag):return sum(c.tag==tag for c in self.cells)


class ComplementaryPartition:
    def draw(self,counts,rng):
        if max(counts)>2**53:
            raise ArithmeticError('Binomial sampling exceeds the configured exact-integer input range.')
        first=tuple(int(rng.binomial(n,.5)) for n in counts)
        return first,tuple(n-d for n,d in zip(counts,first))

    @staticmethod
    def allocation_probability(counts,first):
        if max(counts)>20000:raise ValueError('Exact allocation probability size guard.')
        if any(d<0 or d>n for n,d in zip(counts,first)):return Q(0)
        return Q(math.prod(math.comb(n,d) for n,d in zip(counts,first)),2**sum(counts))


class PopulationKernel:
    """Physical event outcomes retained even on a simultaneous terminal flag."""
    def __init__(self,parameters=BatchParameters(),regions=None,partition=None):
        self.p=parameters;self.regions=ChemicalRegions() if regions is None else regions
        self.partition=ComplementaryPartition() if partition is None else partition
        self.chemistry=ResidentChemistry()

    def initial(self):
        cells=tuple(self.regions.newborn(tag,self.p.N) for tag,count in [('H',self.p.h),('L',self.p.low)] for _ in range(count))
        return Population(4*self.p.N*self.p.M,cells)

    def growth_rate(self,pop,cell):
        return self.p.gamma*Q(pop.precursor,4*self.p.N*self.p.M)*cell.counts[2] if pop.precursor else Q(0)

    def events(self,pop):
        if pop.status!='active':return []
        events=[]
        for i,cell in enumerate(pop.cells):
            for k,ch in enumerate(self.chemistry.channels):
                rate=ch.rate(cell)
                if rate:events.append((i,k,rate))
            rate=self.growth_rate(pop,cell)
            if rate:events.append((i,13,rate))
        return events

    def step(self,pop,index,channel,rng,allocation=None):
        if pop.status!='active':raise ValueError('Terminal records are absorbing.')
        if not 0<=index<len(pop.cells) or not 0<=channel<=13:raise ValueError('Unknown event.')
        cell=pop.cells[index];Qnext=pop.precursor;divisions=pop.divisions
        if channel<13:
            updated=self.chemistry.channels[channel].apply(cell);newcells=(updated,)
            status='active' if self.regions.admits(updated,B) else 'outer_failure'
        else:
            if not self.growth_rate(pop,cell):raise ValueError('Disabled growth event.')
            n=list(cell.counts);n[2]-=1
            updated=Cell(cell.tag,cell.size+1,tuple(n));Qnext-=1
            division=updated.size==2*self.p.N
            if updated.size>2*self.p.N:raise ValueError('Cell exceeded the division threshold.')
            if division:
                first,second=self.partition.draw(updated.counts,rng) if allocation is None else (tuple(allocation),tuple(n-d for n,d in zip(updated.counts,allocation)))
                newcells=(Cell(cell.tag,self.p.N,first),Cell(cell.tag,self.p.N,second));divisions+=1
            else:newcells=(updated,)
            # Outer flags, parent division energy, daughter partition, endpoint.
            if not self.regions.admits(updated,B):status='outer_failure'
            elif division and not self.regions.admits(updated,2*A,closed=True):status='division_energy_failure'
            elif division and any(not self.regions.admits(c,4*A) for c in newcells):status='partition_failure'
            elif Qnext==self.p.N*self.p.M:status='nutrient_endpoint'
            else:status='active'
        outcome=Population(Qnext,pop.cells[:index]+newcells+pop.cells[index+1:],divisions,status)
        assert outcome.precursor+outcome.total_size==pop.precursor+pop.total_size
        assert len(outcome.cells)-len(pop.cells)==divisions-pop.divisions
        assert all(self.p.N<=c.size<2*self.p.N for c in outcome.cells)
        return outcome

    def simulate(self,budget=EVENT_BUDGET,seed=RANDOM_SEED):
        if type(budget) is not int or budget<1:raise ValueError('Use a positive integer event budget.')
        rng=np.random.default_rng(seed);pop=self.initial();time=0.;trace=[]
        status='event_budget';steps=0
        for steps in range(budget):
            events=self.events(pop);rates=np.array([float(e[2]) for e in events]);total=rates.sum()
            wait=float(rng.exponential(1/total))
            if time+wait>float(self.p.horizon):time=float(self.p.horizon);status='deadline';break
            event=events[int(rng.choice(len(events),p=rates/total))]
            try:following=self.step(pop,event[0],event[1],rng)
            except ArithmeticError as exc:status='computational_limit: '+str(exc);break
            pop=following;time+=wait
            if steps%100==0 or pop.status!='active':
                trace.append([time,pop.precursor,pop.total_size,pop.tag_count('H'),pop.tag_count('L'),pop.status])
            if pop.status!='active':status=pop.status;break
        determined=status=='deadline' or pop.status!='active'
        success=self.success(pop) if determined else None
        final=[time,pop.precursor,pop.total_size,pop.tag_count('H'),pop.tag_count('L'),pop.status]
        if not trace or trace[-1]!=final:trace.append(final)
        return {'time':time,'steps_attempted':steps+1,'status':status,'success':success,
                'theorem_parameter_scope':self.p.theorem_scope,'trace':trace}

    def success(self,pop):
        if pop.status!='nutrient_endpoint':return False
        if any(c.readout!=c.tag for c in pop.cells):return False
        H=pop.tag_count('H');L=pop.tag_count('L')
        if not H or not L:return False
        with mp.workdps(80):return bool(mp.log(mp.mpf(H)*self.p.low/(L*self.p.h))>=mpq(self.p.gain))


def mpq(v):
    v=Q(v);return mp.mpf(v.numerator)/v.denominator


class ErrorBudget:
    """High-precision evaluation of proved expressions, not empirical failure rates."""
    def __init__(self,parameters=BatchParameters()):self.p=parameters

    def logs(self):
        p=self.p;N=mp.mpf(p.N);M=mp.mpf(p.M);T=mpq(p.horizon);u=mpq(ALPHA*A)*N
        logT=mp.log(T) if T else -mp.inf
        branches={
            'outer_initial':mp.log(7*M)-12*u,
            'outer_drift':mp.log(M/84)+logT+mp.log(u)-mp.mpf('15.5')*u,
            'division_initial':mp.log(M)+2*u-mpq(LAMBDA)*N,
            'division_drift':mp.log(M/84)+logT+mp.log(u)-mp.mpf('1.5')*u,
            'partition':mp.log(24*M)-N*mpq(DELTA**2/35),
        }
        raw=mp.log(sum(mp.exp(v) for v in branches.values()))
        selection=-mpq(THETA)*N*(mp.log(2)/5-mpq(p.gain))
        deadline=-mpq(THETA)*N*(mpq(Q(9,40)*p.scaled_deadline)-mp.log(4))
        simple=mp.log(M*(32+T/42))-mpq(C)*N
        sharp=mp.log(M*(32+T*u/42))-3*u/2
        double=mp.log(M*(32+T/21))-u
        return branches,{'raw':raw,'simplified':simple,'retained':sharp,'double_exponent':double},selection,deadline

    def evaluate(self):
        with mp.workdps(80):
            branches,chem,sel,deadline=self.logs()
            totals={name:mp.log(sum(mp.exp(v) for v in (value,sel,deadline))) for name,value in chem.items()}
            format=lambda v:mp.nstr(v,35)
            return {'theorem_parameter_scope':self.p.theorem_scope,
                'branch_log10':{k:format(v/mp.log(10)) for k,v in branches.items()},
                'chemical_error':{k:format(mp.exp(v)) for k,v in chem.items()},
                'total_error_log10':{k:format(v/mp.log(10)) for k,v in totals.items()},
                'success_lower':{k:format(max(mp.mpf(0),-mp.expm1(v))) if self.p.theorem_scope else None for k,v in totals.items()},
                'selection_log10':format(sel/mp.log(10)),'deadline_log10':format(deadline/mp.log(10)),
                'selection_margin':format(mp.log(2)/5-mpq(self.p.gain)),
                'deadline_margin':format(mpq(Q(9,40)*self.p.scaled_deadline)-mp.log(4))}

    def sufficient_size(self,error=TARGET_ERROR,double_exponent=False):
        p=self.p
        if not 0<Q(error)<1 or p.gamma>GAMMA_MAX:raise ValueError('Use error in (0,1) and certified gamma.')
        if p.gain!=Q(1,10) or p.scaled_deadline!=8:
            raise ValueError('This closed-form sufficient size is for s=0.1 and gamma*T=8.')
        with mp.workdps(80):
            coefficient=p.M*(32+(8 if double_exponent else 4)/(21*mpq(p.gamma)))
            exponent=mpq(C)*(2 if double_exponent else 1)
            candidates=[mp.mpf(N_MIN),mp.log(3*coefficient/mpq(error))/exponent,
                        mp.mpf(500000)/19*mp.log(3/mpq(error)),2500*mp.log(3/mpq(error))]
            # One extra integer avoids understating a numerically evaluated ceiling.
            return int(mp.ceil(max(candidates)))+1


class AncestralObservables:
    def __init__(self,parameters):self.p=parameters
    def log_phi(self,pop):
        p=self.p
        with mp.workdps(80):
            H=mp.mpf(pop.tag_size('H'));L=mp.mpf(pop.tag_size('L'));W=H+L
            return mpq(THETA)*p.N*(-mp.log(H/(p.N*p.h))+mp.log(L/(p.N*p.low))+mp.mpf(3)/5*mp.log(W/(p.N*p.M)))
    def size_to_count_ratio(self,pop):
        H=pop.tag_count('H');L=pop.tag_count('L')
        if not H or not L:raise ValueError('Both tags must be present.')
        return {'size_odds':Q(pop.tag_size('H'),pop.tag_size('L')),'count_odds':Q(H,L),
                'phase_factor':Q(pop.tag_size('H')*L,pop.tag_size('L')*H)}

    def relative_odds_drift(self,pop,high_z,low_z):
        """L Phi / (beta Phi), useful for checking the source-safe rate sandwich."""
        with mp.workdps(80):
            H=mp.mpf(pop.tag_size('H'));L=mp.mpf(pop.tag_size('L'));W=H+L
            thetaN=mpq(THETA)*self.p.N
            vH=thetaN*(-mp.log1p(1/H)+mp.mpf(3)/5*mp.log1p(1/W))
            vL=thetaN*(mp.log1p(1/L)+mp.mpf(3)/5*mp.log1p(1/W))
            return mpq(high_z)*H*mp.expm1(vH)+mpq(low_z)*L*mp.expm1(vL)

    def log_position_exponential(self,cell,regions):
        return LAMBDA*(cell.size-2*self.p.N)+ALPHA*self.p.N*regions.energy(cell)


def scalar_certificates():
    assert ALPHA*A==2*C
    assert LAMBDA>=Q(7,2)*ALPHA*A and DELTA**2/35>=Q(3,2)*ALPHA*A
    assert 32*GAMMA_MAX*LAMBDA<=ALPHA*A/1344
    assert (4*ALPHA*A-LAMBDA)*N_MIN<=Q(-4479965,32)<-1
    assert 42*16*DELTA**2<A
    assert 200000000+8000000000*N_MIN*GAMMA_MAX**2<=N_MIN*A/672
    assert -Q(3,20000)+Q(6,390625)<0
    # Exact positive-definite matrix checks replace numerical eigenvalue checks.
    for tag,upper in [('L',24),('H',42)]:
        P=sp.Matrix(MATRICES[tag])/10**6
        for matrix in (P-sp.eye(4)/200,upper*sp.eye(4)-P):
            assert all(matrix[:k,:k].det()>0 for k in range(1,5))
    candidates=[Q(h,h+l) for h in range(1,9) for l in range(1,9) if h>l and h+l<=8]
    assert min(candidates)==Q(4,7)
    return {'exact_scalar_margins':True,'exact_matrix_sandwiches':True,'two_founder_frequency':'4/7',
            'scope':'does not rerun full nonlinear finite-count generator or Lean certificates'}


def frozen_composition_diagnostic(p,regions):
    """Deterministic phase illustration, explicitly NOT the resident count model.

    Chemistry is clamped at each stationary z. Exposure q satisfies dE/dt=gamma*Q/Omega.
    Total ancestral size/N is h exp(zH E)+low exp(zL E). Division conserves it,
    and equal daughters give count_i=founders_i*2**floor(z_i E/log(2)).
    """
    zH=float((regions.roots['H'].lo+regions.roots['H'].hi)/2)
    zL=float((regions.roots['L'].lo+regions.roots['L'].hi)/2)
    left=0.;right=10.
    for _ in range(70):
        mid=(left+right)/2
        if p.h*math.exp(zH*mid)+p.low*math.exp(zL*mid)>4*p.M:right=mid
        else:left=mid
    end=(left+right)/2;rows=[]
    for exposure in np.linspace(0,end,201):
        BH=p.h*math.exp(zH*exposure);BL=p.low*math.exp(zL*exposure)
        CH=p.h*2**math.floor(zH*exposure/math.log(2));CL=p.low*2**math.floor(zL*exposure/math.log(2))
        rows.append([exposure,BH,BL,CH,CL,BH/BL,CH/CL,(5*p.M-BH-BL)/(4*p.M)])
    return rows


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'))
    parser.add_argument('--simulate',action='store_true',help='Attempt a budget-limited literal count path.')
    args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
    p=BatchParameters();regions=ChemicalRegions();kernel=PopulationKernel(p,regions)
    initial=kernel.initial();assert initial.precursor+initial.total_size==5*p.N*p.M
    budget=ErrorBudget(p);evaluation=budget.evaluate();checks=scalar_certificates()
    table=[{'N':N,**ErrorBudget(replace(p,N=N)).evaluate()} for N in COPY_NUMBER_SWEEP]
    phase=frozen_composition_diagnostic(p,regions)
    result={'parameters':{'N':p.N,'h':p.h,'low':p.low,'gamma':str(p.gamma),'gain':str(p.gain),'horizon':str(p.horizon)},
            'budget':evaluation,'copy_number_table':table,'scalar_checks':checks,
            'founders':[{'tag':c.tag,'size':c.size,'counts':c.counts,'energy_interval':regions.energy(c).json(),
                         'readout':c.readout} for c in initial.cells],
            'stationary_root_intervals':{tag:v.json() for tag,v in regions.roots.items()},
            'initial_total_event_rate':str(sum(e[2] for e in kernel.events(initial))),
            'sufficient_N_original':budget.sufficient_size() if p.gain==Q(1,10) and p.scaled_deadline==8 and p.gamma<=GAMMA_MAX else None,
            'sufficient_N_double_exponent':budget.sufficient_size(double_exponent=True) if p.gain==Q(1,10) and p.scaled_deadline==8 and p.gamma<=GAMMA_MAX else None,
            'simulation':kernel.simulate() if args.simulate else None,
            'diagnostic_scope':'frozen-composition deterministic phase illustration, not a count-model simulation'}
    (out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
    def csv_file(name,headers,rows):
        with (out/name).open('w',newline='') as stream:
            w=csv.writer(stream);w.writerow(headers);w.writerows(rows)
    csv_file('bound_table.csv',['N','raw_chemical_error','simplified_chemical_error','raw_success_lower','simplified_success_lower'],
             [(r['N'],r['chemical_error']['raw'],r['chemical_error']['simplified'],r['success_lower']['raw'],r['success_lower']['simplified']) for r in table])
    csv_file('phase_diagnostic.csv',['exposure','high_size_per_N','low_size_per_N','high_count','low_count','size_odds','count_odds','precursor_fraction'],phase)
    sweep=[]
    for N in np.geomspace(float(N_MIN),4e22,180):
        d=ErrorBudget(replace(p,N=int(N))).evaluate()
        sweep.append([int(N),d['total_error_log10']['raw'],d['total_error_log10']['simplified'],d['success_lower']['raw'],d['success_lower']['simplified']])
    csv_file('copy_number_sweep.csv',['N','raw_total_log10','simplified_total_log10','raw_success_lower','simplified_success_lower'],sweep)
    lines=[f'Parameter scope for the theorem: {p.theorem_scope}; N={p.N}; M={p.M}.',
           f'Raw-formula success lower bound: {evaluation["success_lower"]["raw"]}.',
           f'Simplified-formula success lower bound: {evaluation["success_lower"]["simplified"]}.',
           f'Sufficient N for error {TARGET_ERROR} (original): {result["sufficient_N_original"]}.',
           f'Sufficient N using double exponent: {result["sufficient_N_double_exponent"]}.',
           'Two-founder success implies high-state frequency at least 4/7.',
           'The phase diagnostic clamps chemistry; it does not estimate chemical fidelity.',
           'The count kernel is available; default outputs evaluate bounds rather than simulating an astronomical batch.']
    (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),layout='constrained');data=np.array(sweep,dtype=float)
    for col,color,label in [(1,'#287a99','Raw chemical sum + population tails'),(2,'#777777','Simplified chemical envelope + tails')]:
        axes[0].plot(data[:,0]/1e22,data[:,col],color=color,label=label)
        axes[1].plot(data[:,0]/1e22,data[:,col+2],color=color,label=label)
    axes[0].axhline(0,color='black',lw=.7);axes[0].set(ylabel='log10 of total error bound',ylim=(-35,12))
    axes[1].set(ylabel='Success probability lower bound',ylim=(-.02,1.02));axes[1].legend(fontsize=7)
    for ax in axes:ax.set_xlabel('Newborn size N / 10^22');ax.grid(alpha=.2)
    fig.suptitle('Success bounds versus newborn molecule\ncount')
    fig.savefig(out/'bounds.png',dpi=180);fig.savefig(out/'bounds.svg');plt.close(fig)
    fig,axes=plt.subplots(1,2,figsize=(10,4),layout='constrained');data=np.array(phase)
    axes[0].plot(data[:,0],data[:,5],label='Ancestral size odds');axes[0].step(data[:,0],data[:,6],where='post',label='Cell-count odds')
    axes[0].set(ylabel='High / low',xlabel='Integrated precursor exposure');axes[0].legend(fontsize=8)
    axes[1].plot(data[:,0],data[:,7]);axes[1].axhline(.25,color='black',ls='--',label='Batch endpoint')
    axes[1].set(ylabel='Remaining precursor / initial precursor',xlabel='Integrated precursor exposure');axes[1].legend(fontsize=8)
    for ax in axes:ax.grid(alpha=.2)
    fig.suptitle('Lineage size and cell count with fixed\ncompositions')
    fig.savefig(out/'division_phase.png',dpi=180);fig.savefig(out/'division_phase.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()