"""Exact finite RAF bounds, productive histories, and a budgeted profile evaluator.

Three distinct observables: finite catalytic RAF probability, static escape,
and infinite survival. No Monte Carlo curve is presented as the limiting profile.
"""
from __future__ import annotations
import argparse
from collections import Counter
import csv
from dataclasses import dataclass
from fractions import Fraction as Q
from functools import lru_cache
from itertools import product
import hashlib
import json
import math
from pathlib import Path
import platform

# EDITABLE INPUTS -----------------------------------------------------------
WORD_CAP = 16
INTENSITY = Q(1,10000)                # lambda; catalytic p=lambda*n/R_n
STATIC_OPENNESS = Q(1,10000)          # a; distinct from catalytic p and 1-exp(-lambda)
PROFILE_TOLERANCE = Q(1,10**12)
HISTORY_DEPTH = 2                    # exhaustive tiny census, not full cap-16 fields
MAX_HISTORY_BOUND = 80
MAX_STATIC_DECISIONS = 20000          # resource exhaustion raises; never means zero
MAX_SEED_INDEX = 10000
MAX_EXACT_BITS = 1000000
SEED_COST_OPENNESS = Q(1,2)
SEED_COST_M = 100
MANUSCRIPT_SHA256 = '1d561295f36a0415f2aa64ea12e4a310e078a3fbbc3be17eb3b8122d5e388c45'
# --------------------------------------------------------------------------
FOOD=frozenset(''.join(v) for k in (1,2) for v in product('01',repeat=k))
C_STAR=2*7**64*81


class ResourceLimit(RuntimeError):
    """Incomplete computation; carries no probability verdict."""


@dataclass(frozen=True)
class Interval:
    lower: Q
    upper: Q
    def __post_init__(self):
        if self.lower>self.upper:raise ValueError('Reversed interval.')
    @property
    def width(self):return self.upper-self.lower
    def json(self):return dict(lower=str(self.lower),upper=str(self.upper),width=str(self.width))


@dataclass(frozen=True,order=True)
class Channel:
    left: str
    right: str
    product: str
    @property
    def endpoints(self):return frozenset((self.left,self.right,self.product))
    def enabled(self,available):return self.product in available or self.left in available and self.right in available
    def productive(self,available):return self.enabled(available) and not self.endpoints<=available


class PolymerCatalogue:
    """Ordered split identity or commuting-factor quotient; product never changes."""
    def __init__(self,mode):
        if mode not in ('sp','qt'):raise ValueError('Use sp (split) or qt (quotient).')
        self.mode=mode
    def channel(self,u,v):
        if not u or not v or set(u+v)-set('01'):raise ValueError('Nonempty binary factors required.')
        if self.mode=='qt' and u+v==v+u and (len(u),u)>(len(v),v):u,v=v,u
        return Channel(u,v,u+v)
    @staticmethod
    def words(n):
        return [''.join(v) for k in range(1,n+1) for v in product('01',repeat=k)]
    def channels(self,n):
        return tuple(sorted({self.channel(w[:j],w[j:]) for w in self.words(n) for j in range(1,len(w))}))
    @staticmethod
    def molecule_count(n):return 2**(n+1)-2
    @staticmethod
    @lru_cache(None)
    def primitive_count(n):return 2**n-sum(PolymerCatalogue.primitive_count(d) for d in range(1,n) if n%d==0)
    def channel_count(self,n):
        if not isinstance(n,int) or n<2:raise ValueError('Integer cap at least two required.')
        split=(n-2)*2**(n+1)+4
        if self.mode=='sp':return split
        loss=0
        for d in range(1,n//3+1):
            limit=n//d
            pairs=sum(max(0,limit-2*i) for i in range(1,limit//2+1))
            loss+=self.primitive_count(d)*pairs
        return split-loss
    @property
    def growth_count(self):return 32 if self.mode=='sp' else 30
    @property
    def gateway_count(self):return self.growth_count+4
    @property
    def singleton_count(self):return 24+7*self.growth_count
    def productive(self,available):
        candidates={self.channel(u,v) for u in available for v in available if u+v not in available}
        candidates.update(self.channel(w[:j],w[j:]) for w in available for j in range(1,len(w))
                          if w[:j] not in available or w[j:] not in available)
        return tuple(sorted(candidates))


class ReversibleClosure:
    @staticmethod
    def compute(channels,food=FOOD):
        available=set(food)
        while True:
            before=len(available)
            for ch in channels:
                if ch.enabled(available):available.update(ch.endpoints)
            if len(available)==before:return frozenset(available)

    @staticmethod
    def maximal_raf(channels,catalysis):
        """Deletion algorithm with one catalyst mark serving both directions.

        Product self-catalysis is allowed; no catalyzed startup order is imposed.
        catalysis maps a Channel to its set of catalyst word types.
        """
        remaining=set(channels)
        while remaining:
            available=ReversibleClosure.compute(remaining)
            good={ch for ch in remaining if ch.endpoints<=available and set(catalysis.get(ch,()))&available}
            if good==remaining:return frozenset(good)
            remaining=good
        return frozenset()


class ProductiveHistories:
    def __init__(self,catalogue):self.catalogue=catalogue
    @staticmethod
    def budget(j):return (6+2*j)**2+4*2**j-j
    @lru_cache(None)
    def coefficient(self,r):
        if not isinstance(r,int) or r<0:raise ValueError('Nonnegative integer depth required.')
        return 1 if r==0 else self.catalogue.growth_count*math.prod(self.budget(j) for j in range(1,r))
    def bound(self,a,max_r=MAX_HISTORY_BOUND):
        a=probability(a)
        if max_r<0:raise ValueError('Nonnegative maximum depth required.')
        values=[Q(self.coefficient(r))*a**r for r in range(max_r+1)]
        r=min(range(len(values)),key=values.__getitem__)
        return r,values[r]
    def census(self,depth=HISTORY_DEPTH):
        """Deduplicate states, but retain ordered-path multiplicities and witness support."""
        states={FOOD:(1,frozenset())};rows=[]
        for j in range(depth+1):
            nxt={};transitions=0;maxchoices=maxT=maxQ=0
            for available,(count,used) in sorted(states.items(),key=lambda v:tuple(sorted(v[0]))):
                assert ReversibleClosure.compute(used)==available and all(ch.endpoints<=available for ch in used)
                T=sum(map(len,available));slots=T-len(available)
                assert len(available)<=6+2*j and max(map(len,available))<=2**(j+1)
                assert T<=2**(j+2)+6 and slots<=4*2**j-j
                choices=self.catalogue.productive(available);maxchoices=max(maxchoices,len(choices));maxT=max(maxT,T);maxQ=max(maxQ,slots)
                assert len(choices)<=(self.catalogue.growth_count if j==0 else self.budget(j))
                for ch in choices:
                    target=available|ch.endpoints;support=used|{ch}
                    assert ch not in used and ReversibleClosure.compute(support)==target
                    transitions+=count
                    if j<depth:
                        old=nxt.get(target,(0,support));nxt[target]=(old[0]+count,old[1])
            rows.append(dict(depth=j,states=len(states),ordered_histories=sum(v[0] for v in states.values()),
                next_histories=transitions,max_choices=maxchoices,max_total_length=maxT,max_cleavage_slots=maxQ))
            assert transitions<=self.coefficient(j+1)
            states=nxt
        return rows


def probability(v):
    v=Q(v)
    if not 0<=v<=1:raise ValueError('Probability outside [0,1].')
    return v


class FinitePrediction:
    def __init__(self,catalogue,n):
        if not isinstance(n,int) or n<4:raise ValueError('Finite prediction requires n>=4.')
        self.catalogue=catalogue;self.n=n;self.M=catalogue.molecule_count(n);self.R=catalogue.channel_count(n)
    def canonical_probability(self,intensity):return probability(Q(intensity)*self.n/self.R)
    def singleton(self,p):return 1-(1-probability(p))**self.catalogue.singleton_count
    def bounds(self,p):
        p=probability(p);a=self.catalogue.singleton_count;g=self.catalogue.gateway_count
        lower=self.singleton(p);bonf=max(Q(0),a*p-math.comb(a,2)*p*p)
        # All RAFs require a gateway. This exact bound covers caps too small for r>=1.
        candidates=[dict(kind='singleton_remainder',r=None,upper=min(Q(1),lower+math.comb(self.M*self.R,2)*p*p))]
        # Bernoulli union bound avoids raising a huge denominator to M*G.
        candidates.append(dict(kind='gateway_union',r=None,upper=min(Q(1),g*self.M*p)))
        r=1;history=ProductiveHistories(self.catalogue)
        while 2**(r+2)<=self.n:
            B=2**(r+1);hist=history.coefficient(r)*(self.M*p)**r;short=g*self.catalogue.molecule_count(B)*p
            candidates.append(dict(kind='history_plus_short_gateway',r=r,B=B,history=hist,short=short,upper=min(Q(1),hist+short)))
            r+=1
        best=min(candidates,key=lambda v:v['upper'])
        assert lower<=best['upper']
        return dict(n=self.n,molecules=self.M,channels=self.R,p=p,singleton=lower,bonferroni_lower=bonf,
                    interval=Interval(lower,best['upper']).json(),best=best,candidates=candidates)


class StaticField:
    """Exact independent-open finite field. Adaptive queries reveal only fresh marks.

    Every previously open channel is settled, so current closure and closed marks
    suffice for recursion. Memoization shares subproblems; budget exhaustion raises.
    """
    def __init__(self,channels):self.channels=tuple(sorted(set(channels)))
    def escape(self,a,B,max_decisions=MAX_STATIC_DECISIONS):
        a=probability(a)
        if B<2:raise ValueError('Escape threshold must be at least food length two.')
        if a==0:return Q(0)
        if a==1:return Q(int(max(map(len,ReversibleClosure.compute(self.channels)))>B))
        decisions=0
        @lru_cache(None)
        def visit(available,closed):
            nonlocal decisions
            if max(map(len,available))>B:return Q(1)
            enabled=next((i for i,ch in enumerate(self.channels) if i not in closed and ch.productive(available)),None)
            if enabled is None:return Q(0)
            decisions+=1
            if decisions>max_decisions:raise ResourceLimit('Static decision budget exhausted; escape probability not computed.')
            ch=self.channels[enabled]
            return a*visit(available|ch.endpoints,closed)+(1-a)*visit(available,closed|{enabled})
        return visit(FOOD,frozenset())

    @classmethod
    def full_escape(cls,catalogue,a,B,max_decisions=MAX_STATIC_DECISIONS,max_channels=200):
        a=probability(a)
        if B==2:return 1-(1-a)**catalogue.growth_count
        if catalogue.channel_count(2*B)>max_channels:
            raise ResourceLimit('Full witness-cap catalogue exceeds channel budget; no escape verdict.')
        return cls(catalogue.channels(2*B)).escape(a,B,max_decisions)


class RepairConstruction:
    def __init__(self,catalogue):self.catalogue=catalogue
    def support(self,carrier,target):
        # Append letters, then split off the whole target. Repeated identities merge.
        support=set();word=carrier
        for letter in target:
            ch=self.catalogue.channel(word,letter);support.add(ch);word+=letter
        support.add(self.catalogue.channel(carrier,target))
        assert len(support)<=len(target)+1 and all(len(ch.product)>len(carrier) for ch in support)
        assert target in ReversibleClosure.compute(support,FOOD|{carrier})
        return tuple(sorted(support))


class RationalConversion:
    @staticmethod
    def openness(z,tolerance):
        z=Q(z);tol=Q(tolerance)
        if z<0 or tol<=0:raise ValueError('Nonnegative intensity and positive tolerance required.')
        if z==0:return Interval(Q(0),Q(0))
        total=Q(1);term=Q(1);n=0
        while True:
            if n+2>=2*z:
                tail=2*term*z/(n+1);answer=Interval(1-1/total,1-1/(total+tail))
                if answer.width<=tol:return answer
            n+=1;term*=z/n;total+=term
    @staticmethod
    def intensity(a,tolerance):
        a=probability(a);tol=Q(tolerance)
        if a==1 or tol<=0:raise ValueError('a<1 and positive tolerance required.')
        total=Q(0);power=a;n=1
        while True:
            total+=power/n;remainder=power*a/((n+1)*(1-a))
            if remainder<=tol:return Interval(total,total+remainder)
            n+=1;power*=a


class EffectiveApproximation:
    """Resource-bounded implementation of the proof, with useful bound shortcuts.

    Failure to fit the resource budget is explicit. It is not a numerical estimate
    of the survival profile and does not contradict mathematical computability.
    """
    def __init__(self,catalogue):self.catalogue=catalogue
    def seed(self,q,m,max_index=MAX_SEED_INDEX,max_bits=MAX_EXACT_BITS):
        q=probability(q)
        if not q or not isinstance(m,int) or m<2:raise ValueError('q>0, integer m>=2 required.')
        b=q if self.catalogue.mode=='sp' else q/2;base=1-b*b;factor=C_STAR*(m*(m+1)+1)
        def passes(k):
            if k>max_index or k*max(base.numerator.bit_length(),base.denominator.bit_length())>max_bits:
                raise ResourceLimit('Exact seed search exceeds index/bit budget; seed not computed.')
            return factor*base**k<=1
        hi=1
        while not passes(hi):hi*=2
        lo=0
        while hi-lo>1:
            mid=(lo+hi)//2
            if passes(mid):hi=mid
            else:lo=mid
        return 10*hi
    def record_upper(self,q,L,tolerance):
        q=probability(q);tol=Q(tolerance)
        if not q or L<1 or tol<=0:raise ValueError('q>0, L>=1 and positive tolerance required.')
        if (L+1)*max(q.numerator.bit_length(),q.denominator.bit_length())>MAX_EXACT_BITS:raise ResourceLimit('Record construction exceeds rational bit budget.')
        # (1-t)^r <= 1/(1+r*t) by Bernoulli, with t=q^(L+1).
        # This computable index is conservative, not the paper's least index.
        return max(1,math.ceil((Q(2**(L+1))/tol-1)/q**(L+1)))
    def evaluate(self,a,tolerance,max_cap_bits=20):
        a=probability(a);tol=Q(tolerance)
        if tol<=0:raise ValueError('Positive tolerance required.')
        if a in (0,1):return dict(method='exact_endpoint',interval=Interval(a,a).json())
        r,upper=ProductiveHistories(self.catalogue).bound(a)
        if upper<tol:return dict(method='history_upper_bound',history_depth=r,interval=Interval(Q(0),upper).json())
        m=math.ceil(4/tol)+2;L=self.seed(a,m);record=self.record_upper(a,L,tol/4)
        if record>max_cap_bits:raise ResourceLimit(f'Profile cap cannot be materialized: seed L={L}, conservative record-index bit length={record.bit_length()}; no interior evaluation.')
        B=2**record*(L+2)-L
        e=StaticField.full_escape(self.catalogue,a,B)
        return dict(method='finite_escape_with_uniform_error',seed=L,record=record,
                    interval=Interval(max(Q(0),e-tol/2),min(Q(1),e)).json())


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)
    models={};bounds_rows=[];finite_rows=[];history_rows=[]
    for mode in ('sp','qt'):
        cat=PolymerCatalogue(mode);hist=ProductiveHistories(cat);pred=FinitePrediction(cat,WORD_CAP)
        gateways=[ch for ch in cat.channels(4) if ch.enabled(FOOD)]
        witnesses=[(x,ch) for ch in gateways for x in ReversibleClosure.compute([ch])]
        assert len(witnesses)==cat.singleton_count and len(gateways)==cat.gateway_count
        for n in range(2,10):assert len(cat.channels(n))==cat.channel_count(n)
        census=hist.census();history_rows.extend((mode,*row.values()) for row in census)
        finite=pred.bounds(pred.canonical_probability(INTENSITY));certificates=[]
        for power,r,target in [(3,8,5),(4,12,15),(6,18,45),(20,65,629)]:
            bound=Q(hist.coefficient(r),10**(power*r));assert bound<Q(1,10**target)
            certificates.append(dict(a=f'1e-{power}',r=r,bound=str(bound),strictly_below=f'1e-{target}'))
        evaluator=EffectiveApproximation(cat);L=evaluator.seed(SEED_COST_OPENNESS,SEED_COST_M)
        record=evaluator.record_upper(SEED_COST_OPENNESS,L,Q(1,100))
        # This estimate describes the LEAST record scale; the exact conservative index above is larger.
        least_log2=(L+1)+math.log2((L+1)*math.log(2)+math.log(100))
        try:interior=evaluator.evaluate(Q(1,2),Q(1,100))
        except ResourceLimit as exc:interior=dict(status='resource_limit',reason=str(exc))
        repair=RepairConstruction(cat).support('000','101')
        toy=tuple(sorted(set(repair)|{cat.channel('0','00'),cat.channel('00','0')}))
        toy_escape=StaticField(toy).escape(Q(1,2),4)
        models[mode]=dict(finite=finite,gateway_count=len(gateways),singleton_coordinates=len(witnesses),history_census=census,
            low_openness_certificates=certificates,configured_profile=evaluator.evaluate(STATIC_OPENNESS,PROFILE_TOLERANCE),
            seed_cost=dict(q=str(SEED_COST_OPENNESS),m=SEED_COST_M,L=L,conservative_record_index_bit_length=record.bit_length(),
                least_record_log2_approx=least_log2,note='Cost of the proof, not an estimate of the profile. Conservative index uses a rational Bernoulli bound.'),
            interior_attempt=interior,repair=[(c.left,c.right,c.product) for c in repair],restricted_field_escape=dict(a='1/2',B=4,value=str(toy_escape),
                channels=[(c.left,c.right,c.product) for c in toy],scope='Only this explicit small subnetwork is random; omitted channels closed. Not full polymer escape or survival.'),
            full_escape_B2_at_half=str(StaticField.full_escape(cat,Q(1,2),2)))
        for exponent in range(2,21):
            a=Q(1,10**exponent);r,v=hist.bound(a);bounds_rows.append((mode,exponent,r,str(v),math.log10(v.numerator)-math.log10(v.denominator)))
        for exponent in range(2,10):
            lam=Q(1,10**exponent);f=pred.bounds(pred.canonical_probability(lam));_,limit=hist.bound(lam)
            finite_rows.append((mode,exponent,str(f['singleton']),f['interval']['upper'],str(limit),float(f['singleton']),float(Q(f['interval']['upper'])),float(limit)))
    conversion=RationalConversion.openness(INTENSITY,Q(1,10**16));inverse=RationalConversion.intensity(Q(1,2),Q(1,10**10))
    assert conversion.upper<INTENSITY
    result=dict(manuscript_sha256=MANUSCRIPT_SHA256,models=models,intensity_to_openness=conversion.json(),half_openness_to_intensity=inverse.json(),
        dependence_counterexample=dict(model='One shared Bernoulli(a) mark opens all channels or none.',survival='a exactly',scope='Same marginals as the independent model, different survival. History powers require independence.'),
        scope='Exact finite algebra and enumeration plus imported universal theorems. No Lean replay; no interior survival plot or kinetic interpretation.')
    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('history_bounds.csv',['mode','log10_inverse_a','r','upper_exact','log10_upper_display'],bounds_rows)
    table('finite_vs_limit.csv',['mode','log10_inverse_lambda','finite_lower_exact','finite_upper_exact','limit_upper_exact','finite_lower_display','finite_upper_display','limit_upper_display'],finite_rows)
    table('history_census.csv',['mode',*models['sp']['history_census'][0].keys()],history_rows)
    lines=['Split/quotient catalogues, gateway and singleton coordinates checked by literal enumeration.',
        f'Configured n={WORD_CAP}, lambda={INTENSITY}:']
    for mode,data in models.items():
        f=data['finite'];lines.append(f'  {mode}: singleton lower {float(f["singleton"]):.9g}; RAF upper {float(Q(f["interval"]["upper"])):.9g}; p={f["p"]}.')
    lines+=['Both profiles: exact S(1e-20) < 1e-629. These are upper bounds, not the profile values.',
            'Interior evaluation at a=1/2 stops at an explicit resource limit; it returns no probability estimate.',
            'One shared channel mark gives survival a, showing why marginal probabilities alone do not justify history powers.']
    (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),layout='constrained')
    for mode,label in [('sp','Split'),('qt','Quotient')]:
        data=[r for r in bounds_rows if r[0]==mode];axs[0].plot([r[1] for r in data],[r[4] for r in data],'.-',label=label)
    # Earlier sparse-support envelope: exact integer binomials, float only at rendering.
    old=[math.log10(math.comb(PolymerCatalogue('sp').channel_count(2*(2**(r+1)+2)),r)) for r in range(1,8)]
    xs=range(2,21);axs[0].plot(xs,[min([0]+[v-r*x for r,v in enumerate(old,1)]) for x in xs],'--',color='gray',label='Earlier sparse-support bound')
    axs[0].set(xlabel='log10(1 / static openness a)',ylabel='log10 upper bound on survival',title='Upper bounds on infinite polymer survival');axs[0].legend(fontsize=8)
    for mode,color in [('sp','#246d91'),('qt','#bd5a24')]:
        d=[r for r in finite_rows if r[0]==mode];x=[r[1] for r in d]
        axs[1].fill_between(x,[math.log10(r[5]) for r in d],[math.log10(r[6]) for r in d],color=color,alpha=.2,label=f'{mode}: finite RAF interval')
        axs[1].plot(x,[math.log10(r[7]) for r in d],'--',color=color,label=f'{mode}: limit upper bound')
    axs[1].set(xlabel='log10(1 / intensity lambda)',ylabel='log10 probability / bound',title=f'Finite cap {WORD_CAP} and infinite-survival\nbounds');axs[1].legend(fontsize=7)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'bounds.png',dpi=180);fig.savefig(out/'bounds.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,3.6),layout='constrained')
    for mode,label in [('sp','Split'),('qt','Quotient')]:
        rows=models[mode]['history_census'];axs[0].plot([r['depth']+1 for r in rows],[r['next_histories'] for r in rows],'.-',label=f'{label}: counted')
        h=ProductiveHistories(PolymerCatalogue(mode));axs[0].plot([r['depth']+1 for r in rows],[h.coefficient(r['depth']+1) for r in rows],':',label=f'{label}: bound')
    axs[0].set(yscale='log',xlabel='Productive steps',ylabel='Ordered histories',title='Counts and upper bounds for productive\nhistories');axs[0].legend(fontsize=7)
    x=[0,1];axs[1].bar(x,[models[m]['seed_cost']['L'] for m in ('sp','qt')],color=['#246d91','#bd5a24']);axs[1].set_xticks(x,['Split','Quotient']);axs[1].set(ylabel='Required seed length L',title=f'Proof cost at q={SEED_COST_OPENNESS}, m={SEED_COST_M}')
    for i,m in enumerate(('sp','qt')):axs[1].text(i,models[m]['seed_cost']['L'],str(models[m]['seed_cost']['L']),ha='center',va='bottom')
    axs[1].margins(y=.15)
    for ax in axs:ax.grid(axis='y',alpha=.2)
    fig.savefig(out/'construction.png',dpi=180);fig.savefig(out/'construction.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__)),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()