"""Reaction identity, exact OR laws, reversible RAF peeling, and finite-size diagnostics."""
# EDITABLE STUDY INPUTS. Intensities are dimensionless, not kinetic rates.
CATALYTIC_CAP = 6
INTENSITIES = ('0.05', '0.1', '0.25', '0.5', '1')
MONTE_CARLO_TRIALS = 256
RANDOM_SEED = 36092026
CENSUS_ENUMERATION_MAX = 12
CENSUS_FORMULA_MAX = 40
TOY_PROBABILITIES = ('1/3', '1/2')
ESCAPE_CUTOFFS = (2, 3, 4)
ESCAPE_OPENNESS = '1/10'
MAX_EXACT_CATALYTIC_BITS = 16
MAX_EXPLICIT_CAP = 12

from dataclasses import dataclass
from fractions import Fraction as Q
from itertools import product
from collections import Counter
import argparse
import csv
import hashlib
import json
import math
from pathlib import Path
import platform
import numpy as np

MANUSCRIPT_SHA256='99c0aaf9ca053f4e56a901b84a96ce0c67bba99106b4616f3b5bc17474fd6b1b'


def probability(p):
    if isinstance(p,float):raise ValueError('Use an integer or rational string for a model probability.')
    p=Q(p)
    if not 0<=p<=1:raise ValueError('Probability outside [0,1]; no clamp is applied.')
    return p


def words(n):return tuple(''.join(w) for k in range(1,n+1) for w in product('AB',repeat=k))


@dataclass(frozen=True,order=True)
class Description:
    u:str
    v:str

    def __post_init__(self):
        if not self.u or not self.v:raise ValueError('Nonempty factors required.')
    @property
    def w(self):return self.u+self.v
    def canonical(self):
        return Description(self.v,self.u) if self.u+self.v==self.v+self.u and len(self.v)<len(self.u) else self
    def label(self):return self.u+' + '+self.v+' <-> '+self.w


class ReversibleNetwork:
    def __init__(self,molecules,food,channels):
        self.molecules=tuple(molecules);self.channels=tuple(channels);self.index={m:i for i,m in enumerate(self.molecules)}
        if len(self.index)!=len(self.molecules) or not set(food)<=set(self.molecules):raise ValueError('Unique molecules and a valid food set required.')
        self.food=sum(1<<self.index[x] for x in set(food));self.all_molecules=(1<<len(self.molecules))-1
        self.factors=[];self.products=[];self.endpoints=[]
        for r in self.channels:
            if any(x not in self.index for x in (r.u,r.v,r.w)):raise ValueError('Channel endpoint missing from molecule catalogue.')
            factors=(1<<self.index[r.u])|(1<<self.index[r.v]);prod=1<<self.index[r.w]
            self.factors.append(factors);self.products.append(prod);self.endpoints.append(factors|prod)
    def closure(self,active,seed=None):
        H=self.food if seed is None else seed
        if H<0 or H&~self.all_molecules:raise ValueError('Invalid seed mask.')
        active=tuple(active)
        if any(type(j) is not int or not 0<=j<len(self.channels) for j in active):raise ValueError('Unknown channel.')
        while True:
            new=H
            for j in active:
                if H&self.products[j] or H&self.factors[j]==self.factors[j]:new|=self.endpoints[j]
            if new==H:return H
            H=new
    def is_raf(self,active,marks):
        active=tuple(active);H=self.closure(active)
        return bool(active) and all(self.endpoints[j]&H==self.endpoints[j] and marks[j]&H for j in active)
    def gateways(self,productive=False):
        return tuple(j for j in range(len(self.channels)) if (self.products[j]&self.food or self.factors[j]&self.food==self.factors[j]) and (not productive or not self.products[j]&self.food))
    def generated_words(self,mask):return [w for i,w in enumerate(self.molecules) if mask>>i&1]


class ChannelQuotient:
    def __init__(self,split):
        self.split=split;canonical=sorted(set(r.canonical() for r in split.channels));lookup={r:j for j,r in enumerate(canonical)}
        self.mapping=tuple(lookup[r.canonical()] for r in split.channels)
        fibres=[[] for _ in canonical]
        for i,j in enumerate(self.mapping):fibres[j].append(i)
        self.fibres=tuple(tuple(f) for f in fibres)
        self.quotient=ReversibleNetwork(split.molecules,split.generated_words(split.food),canonical)
        if any(len(f)>2 for f in self.fibres):raise ValueError('Input repeats identical descriptions; this binary-word quotient assumes each ordered pair appears once.')
        if len(set(split.channels))!=len(split.channels):raise ValueError('Duplicate records require their own dependence model.')
    @classmethod
    def binary(cls,n):
        if type(n) is not int or not 2<=n<=MAX_EXPLICIT_CAP:raise ValueError('Explicit catalogue cap outside implementation budget.')
        molecules=words(n);return cls(ReversibleNetwork(molecules,words(2),[Description(w[:i],w[i:]) for w in molecules for i in range(1,len(w))]))
    def project_marks(self,marks):
        if len(marks)!=len(self.split.channels):raise ValueError('Wrong split matrix size.')
        result=[]
        for fibre in self.fibres:
            value=0
            for i in fibre:value|=marks[i]
            result.append(value)
        return tuple(result)
    def or_parameters(self,p):
        p=probability(p);return tuple(1-(1-p)**len(f) for f in self.fibres)
    def lift_witness(self,selected,split_marks):
        """Choose representatives using catalytic witnesses, not arbitrary labels."""
        projected=self.project_marks(split_marks);selected=tuple(selected)
        if not self.quotient.is_raf(selected,projected):raise ValueError('Input is not a quotient RAF for the projected marks.')
        H=self.quotient.closure(selected)
        lifted=tuple(next(i for i in self.fibres[j] if split_marks[i]&H) for j in selected)
        if not self.split.is_raf(lifted,split_marks) or self.split.closure(lifted)!=H:raise ArithmeticError('Witness lift failed.')
        return lifted


class CatalyticPeeler:
    def __init__(self,network):self.network=network
    def run(self,marks,prefix=False):
        q=self.network
        if len(marks)!=len(q.channels) or any(type(m) is not int or m<0 or m&~q.all_molecules for m in marks):raise ValueError('One valid catalyst mask per reversible channel required.')
        active=frozenset(j for j,m in enumerate(marks) if m);history=[tuple(sorted(active))]
        for _ in range(len(q.channels)+1):
            H=q.closure(active);pool=(1<<H.bit_count())-1 if prefix else H
            next_active=frozenset(j for j,m in enumerate(marks) if m&pool)
            if not next_active<=active:raise ArithmeticError('Peeling history is not decreasing.')
            if next_active==active:
                usable=tuple(j for j in sorted(active) if q.endpoints[j]&H==q.endpoints[j])
                if not prefix and usable and not q.is_raf(usable,marks):raise ArithmeticError('Terminal RAF failed literal check.')
                return {'history':tuple(history),'closure':H,'usable':usable,'has_raf':bool(usable)}
            active=next_active;history.append(tuple(sorted(active)))
        raise ArithmeticError('Peeling did not terminate within its finite bound.')


class IndependentCatalysis:
    def __init__(self,rng):self.rng=rng
    def sample(self,network,p):
        """Binomial count plus uniform subset = independent Bernoulli coordinates."""
        p=float(probability(p));M=len(network.molecules);counts=self.rng.binomial(M,p,len(network.channels));marks=[]
        for count in counts:
            mask=0
            for i in self.rng.choice(M,size=int(count),replace=False):mask|=1<<int(i)
            marks.append(mask)
        return tuple(marks)


class ExactCatalyticLaw:
    def __init__(self,network):
        self.network=network;self.M=len(network.molecules);self.bits=self.M*len(network.channels)
        if self.bits>MAX_EXACT_CATALYTIC_BITS:raise ValueError('Exact catalytic enumeration budget exceeded.')
    def configurations(self):
        for mask in range(1<<self.bits):yield tuple((mask>>(j*self.M))&((1<<self.M)-1) for j in range(len(self.network.channels)))
    def probability(self,parameters,event=None):
        parameters=tuple(probability(p) for p in parameters)
        if len(parameters)!=len(self.network.channels):raise ValueError('One parameter per channel required.')
        peeler=CatalyticPeeler(self.network);event=event or (lambda marks:peeler.run(marks)['has_raf']);total=Q(0)
        for marks in self.configurations():
            if not event(marks):continue
            weight=Q(1)
            for p,m in zip(parameters,marks):k=m.bit_count();weight*=p**k*(1-p)**(self.M-k)
            total+=weight
        return total
    def history_polynomials(self,prefix=False):
        """Counts by history and number of marked coordinates, valid at every p."""
        histories={};peeler=CatalyticPeeler(self.network)
        for marks in self.configurations():
            history=peeler.run(marks,prefix)['history'];histories.setdefault(history,Counter())[sum(m.bit_count() for m in marks)]+=1
        return histories
    def history_probability(self,history,p):
        """Exact column-deletion law, including the final stability step."""
        p=probability(p);history=tuple(frozenset(A) for A in history)
        if not history or any(not A<=set(range(len(self.network.channels))) for A in history):raise ValueError('Invalid history.')
        if any(not B<A for A,B in zip(history,history[1:])):return Q(0)
        extended=history+(history[-1],)
        sizes=[self.M]+[self.network.closure(A).bit_count() for A in history]
        s=1-p;weight=Q(1)
        for j in range(len(self.network.channels)):
            drop=next((t for t,A in enumerate(extended) if j not in A),None)
            weight*=1-s**sizes[-1] if drop is None else s**sizes[0] if drop==0 else s**sizes[drop]-s**sizes[drop-1]
        return weight


class ChannelCensus:
    @staticmethod
    def formula(n):
        if type(n) is not int or not 2<=n<=10000:raise ValueError('Count cap must be an integer from 2 to 10000.')
        primitive=[0]*(n//3+1)
        for d in range(1,len(primitive)):primitive[d]=2**d-sum(primitive[e] for e in range(1,d) if d%e==0)
        deleted=sum(primitive[d]*((n//d-1)**2//4) for d in range(1,len(primitive)))
        M=2**(n+1)-2;R=(n-2)*2**(n+1)+4
        return {'n':n,'M':M,'split':R,'quotient':R-deleted,'deleted':deleted,'elementary_loss_bound':n*(2**((n-1)//2+1)-2)}
    @staticmethod
    def openness(p,k):
        p=probability(p)
        if type(k) is not int or k<0:raise ValueError('Nonnegative pool size required.')
        if not k or not p:return 0.
        if p==1:return 1.
        return -math.expm1(k*math.log1p(-float(p)))


def seed_cutoff(a0,m):
    """Exact integer verification of the sufficient contour cutoff; no enumeration."""
    a0=probability(a0)
    if not a0 or type(m) is not int or m<2:raise ValueError('Positive openness and integer m >= 2 required.')
    b=a0/2;r=1-b*b;C=2*7**64*81;D=m*(m+1)+1
    decay=-math.log1p(-float(b*b))
    if decay<=0 or math.log(C*D)/decay>1000000:raise ValueError('Exact seed-cutoff verification exceeds arithmetic budget.')
    k=max(1,math.ceil(math.log(C*D)/decay))
    holds=lambda j:C*D*pow(r.numerator,j)<=pow(r.denominator,j)
    while not holds(k):k+=1
    while k>1 and holds(k-1):k-=1
    return {'a0':str(a0),'m':m,'k':k,'L':10*k,'exact_integer_inequality':True,'scope':'Sufficient supplied-seed cutoff; not a convergence rate or practical survival estimate.'}


def wilson(successes,trials):
    if trials==0:return None
    z=1.959963984540054;p=successes/trials;den=1+z*z/trials;center=(p+z*z/(2*trials))/den
    radius=z*math.sqrt(p*(1-p)/trials+z*z/(4*trials*trials))/den
    return [0. if successes==0 else max(0.,center-radius),1. if successes==trials else min(1.,center+radius)]


def toy():
    return ChannelQuotient(ReversibleNetwork(('A','AA','AAA','AAAA'),('A','AA'),(Description('A','AA'),Description('AA','A'),Description('A','AAA'))))


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)
    q=toy();split_law=ExactCatalyticLaw(q.split);quotient_law=ExactCatalyticLaw(q.quotient);toyrows=[]
    for pstr in TOY_PROBABILITIES:
        p=probability(pstr);split=split_law.probability([p]*3);projected=quotient_law.probability(q.or_parameters(p));homogeneous=quotient_law.probability([p]*2)
        if split!=projected or not split>homogeneous:raise ArithmeticError('Toy OR-law comparison failed (use interior probabilities for this comparison).')
        toyrows.append({'p':str(p),'split_raf':str(split),'projected_or_raf':str(projected),'homogeneous_quotient_raf':str(homogeneous)})
    original=quotient_law.history_polynomials();prefix=quotient_law.history_polynomials(True)
    if original!=prefix:raise ArithmeticError('Complete-history law failed.')
    for history,counts in original.items():
        for pstr in TOY_PROBABILITIES:
            p=probability(pstr)
            enumerated=sum(count*p**k*(1-p)**(quotient_law.bits-k) for k,count in counts.items())
            if enumerated!=quotient_law.history_probability(history,p):raise ArithmeticError('Column-deletion history weights disagree with enumeration.')
    lifts=0
    for marks in split_law.configurations():
        a=CatalyticPeeler(q.split).run(marks);b=CatalyticPeeler(q.quotient).run(q.project_marks(marks))
        if a['has_raf']!=b['has_raf']:raise ArithmeticError('Configuration-wise RAF projection failed.')
        if b['has_raf']:q.lift_witness(b['usable'],marks);lifts+=1
    census=[]
    for n in range(2,CENSUS_FORMULA_MAX+1):
        row=ChannelCensus.formula(n)
        if n<=CENSUS_ENUMERATION_MAX:
            catalog=ChannelQuotient.binary(n);net=catalog.quotient
            if len(catalog.split.channels)!=row['split'] or len(net.channels)!=row['quotient']:raise ArithmeticError('Count formula disagrees with enumeration.')
            row.update(gateways=len(net.gateways()),productive_gateways=len(net.gateways(True)))
            if n>=4 and (row['gateways'],row['productive_gateways'])!=(34,30):raise ArithmeticError('Gateway census mismatch.')
        census.append(row)
    normalization=[]
    for row in census:
        for lam in INTENSITIES:
            p=Q(lam)*row['n']/row['quotient']
            if not 0<=p<=1:continue
            normalization.append({'n':row['n'],'lambda':lam,'p':str(p),'normalization_ratio':str(Q(row['n']*row['M'],row['quotient'])),
                'full_pool_openness':ChannelCensus.openness(p,row['M']),'limiting_openness':-math.expm1(-float(Q(lam)))})
    rng=np.random.default_rng(RANDOM_SEED);sampler=IndependentCatalysis(rng);net=ChannelQuotient.binary(CATALYTIC_CAP).quotient;finite=[]
    if type(MONTE_CARLO_TRIALS) is not int or not 1<=MONTE_CARLO_TRIALS<=100000:raise ValueError('Trial count outside budget.')
    for lam in INTENSITIES:
        p=probability(Q(lam)*CATALYTIC_CAP/len(net.channels));hits=gateways=food_only=0
        for _ in range(MONTE_CARLO_TRIALS):
            marks=sampler.sample(net,p);result=CatalyticPeeler(net).run(marks);gateway=any(marks[j] for j in net.gateways())
            if result['has_raf'] and not gateway:raise ArithmeticError('RAF without catalytic gateway.')
            hits+=result['has_raf'];gateways+=gateway;food_only+=result['has_raf'] and result['closure']==net.food
        finite.append({'n':CATALYTIC_CAP,'lambda':lam,'p':str(p),'trials':MONTE_CARLO_TRIALS,'rafs':hits,'gateway_samples':gateways,'food_only_rafs':food_only,
            'raf_estimate':hits/MONTE_CARLO_TRIALS,'raf_wilson95':wilson(hits,MONTE_CARLO_TRIALS),
            'gateway_probability_exact_formula_float':ChannelCensus.openness(p,len(net.gateways())*len(net.molecules)),
            'conditional_estimate':hits/gateways if gateways else None,'conditional_wilson95':wilson(hits,gateways)})
    openness=probability(ESCAPE_OPENNESS);escape=[]
    for K in ESCAPE_CUTOFFS:
        if type(K) is not int or K<2:raise ValueError('Escape cutoff must be an integer at least two.')
        network=ChannelQuotient.binary(2*K).quotient;hits=0
        for _ in range(MONTE_CARLO_TRIALS):
            active=tuple(int(j) for j in np.flatnonzero(rng.random(len(network.channels))<float(openness)));H=network.closure(active)
            hits+=any(len(w)>K for w in network.generated_words(H))
        escape.append({'K':K,'cap':2*K,'openness':str(openness),'hits':hits,'trials':MONTE_CARLO_TRIALS,'estimate':hits/MONTE_CARLO_TRIALS,'wilson95':wilson(hits,MONTE_CARLO_TRIALS),
            'exact_E2':str(1-(1-openness)**30) if K==2 else None,
            'scope':'Estimate of finite escape probability e_K, an upper bound on infinite survival; interval is sampling uncertainty, not a certified survival bracket.'})
    result={'toy':toyrows,'projection':{'configurations_checked':1<<split_law.bits,'successful_witness_lifts':lifts},
        'history_identity':{'configurations_checked':1<<quotient_law.bits,'distinct_histories':len(original),'exact_count_polynomials_equal':True},
        'seed_cutoffs':[seed_cutoff('1/10',10),seed_cutoff('1/2',10)],'finite_catalytic':finite,'finite_escape':escape,
        'scope':'The theorem identifies the limit as S_qt(1-exp(-lambda)); these computations do not evaluate that unknown infinite profile or prove a finite-size convergence rate. Lean not rerun.'}
    def write_json(name,data):(out/name).write_text(json.dumps(data,indent=2)+'\n')
    write_json('results.json',result);write_json('census.json',census);write_json('history_polynomials.json',[{'history':[list(A) for A in history],'counts_by_mark_number':dict(sorted(counts.items()))} for history,counts in sorted(original.items())])
    def table(name,header,rows):
        with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    table('normalization.csv',list(normalization[0]),[list(r.values()) for r in normalization])
    table('finite_catalytic.csv',['n','lambda','trials','rafs','food_only_rafs','estimate','wilson_low','wilson_high','conditional_estimate'],[(r['n'],r['lambda'],r['trials'],r['rafs'],r['food_only_rafs'],r['raf_estimate'],*r['raf_wilson95'],r['conditional_estimate']) for r in finite])
    table('toy_probabilities.csv',list(toyrows[0]),[list(r.values()) for r in toyrows])
    lines=[f'Toy p={r["p"]}: split={r["split_raf"]}, OR projection={r["projected_or_raf"]}, homogeneous quotient={r["homogeneous_quotient_raf"]}.' for r in toyrows]
    lines += [f'All {1<<split_law.bits} split configurations preserve RAF existence; {lifts} witnesses lifted.',f'Complete peeling histories have identical count polynomials across {1<<quotient_law.bits} configurations.',
        f'Finite catalytic estimates use n={CATALYTIC_CAP}, {MONTE_CARLO_TRIALS} trials per intensity; they are not the limiting survival curve.',
        '34 gateways bound finite RAF existence; 30 productive gateways bound infinite survival. Food-only RAFs are retained.']
    (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');positions=np.arange(len(toyrows));width=.25
    for i,(key,label) in enumerate((('split_raf','Split independent marks'),('projected_or_raf','Quotient with OR law'),('homogeneous_quotient_raf','Homogeneous quotient'))):axs[0].bar(positions+(i-1)*width,[float(Q(r[key])) for r in toyrows],width,label=label)
    axs[0].set(xticks=positions,xticklabels=[r['p'] for r in toyrows],xlabel='Catalytic probability p',ylabel='Exact toy RAF probability',ylim=(0,1.08),title='RAF probability under three\nchannel-assignment rules');axs[0].legend(fontsize=8,loc='lower right')
    for lam in INTENSITIES:
        points=[r for r in normalization if r['lambda']==lam];line,=axs[1].plot([r['n'] for r in points],[r['full_pool_openness'] for r in points],label='lambda='+lam)
        axs[1].axhline(-math.expm1(-float(Q(lam))),color=line.get_color(),ls=':',lw=.8)
    axs[1].set(xlabel='Maximum word length n',ylabel='Full-pool channel openness (not RAF probability)',title='Channel openness versus maximum polymer\nlength');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(axis='y',alpha=.2)
    fig.savefig(out/'identity_and_normalization.png',dpi=180);fig.savefig(out/'identity_and_normalization.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained');x=[float(Q(r['lambda'])) for r in finite];y=[r['raf_estimate'] for r in finite]
    axs[0].errorbar(x,y,yerr=[[v-r['raf_wilson95'][0] for v,r in zip(y,finite)],[r['raf_wilson95'][1]-v for v,r in zip(y,finite)]],fmt='o-',label='Finite RAF estimate, pointwise Wilson 95%')
    axs[0].plot(x,[r['gateway_probability_exact_formula_float'] for r in finite],'--',label='Finite catalytic gateway upper bound')
    axs[0].set(xlabel='Dimensionless intensity lambda',ylabel='Probability',ylim=(0,1.05),title=f'Actual quotient catalysis at n={CATALYTIC_CAP}');axs[0].legend(fontsize=7)
    y=[r['estimate'] for r in escape]
    axs[1].errorbar([r['K'] for r in escape],y,yerr=[[v-r['wilson95'][0] for v,r in zip(y,escape)],[r['wilson95'][1]-v for v,r in zip(y,escape)]],fmt='o',label='Finite escape estimate, pointwise Wilson 95%')
    axs[1].axhline(float(1-(1-openness)**30),ls='--',color='gray',label='Exact productive-gateway escape E2')
    axs[1].set(xlabel='Escape cutoff K (computed with cap 2K)',ylabel='Finite escape probability',ylim=(0,1.05),title=f'Static open channels at a={ESCAPE_OPENNESS}');axs[1].legend(fontsize=7)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'finite_diagnostics.png',dpi=180);fig.savefig(out/'finite_diagnostics.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    write_json('run_metadata.json',{'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':digest(Path(__file__)),'python':platform.python_version(),'numpy':np.__version__,'random_seed':RANDOM_SEED,'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()