"""Driven-cleavage reactor: literal counts, repeated-window ledger, and bounds."""
from __future__ import annotations

# EDITABLE INPUTS -------------------------------------------------------------
COPY_SCALE = 100_000_000
KINETIC_BIAS = '10'
DUPLEX_RELEASE = '20'
DRIVEN_CLEAVAGE = '1/20'
POST_STARTUP_DURATION = '100'
TARGET_FAILURE = '1/100'
BASAL_RATE = '1/500000000'  # fixed by the theorem; changing it exits its scope
REFERENCE_MOLAR = '1/1000000'
DILUTION_PER_HOUR = '1'
TEMPERATURE_K = '298'
SSA_EVENT_BUDGET = 20_000
SSA_SEED = 25092026
MANUSCRIPT_SHA256 = 'b03c9193c19090688e71a36a46ad31bf33dc9bde1d910c5b0cf3615ae803b05e'
# ---------------------------------------------------------------------------

from dataclasses import dataclass
from fractions import Fraction as Q
import argparse
import csv
import hashlib
import json
import math
from pathlib import Path
import platform

import numpy as np
from scipy.integrate import solve_ivp
from mpmath import mp, iv

mp.dps=70;iv.dps=70
SPECIES=('U','W','X','C1','C2','Z')
A_WEIGHTS=(1,0,1,2,2,2);B_WEIGHTS=(0,1,1,1,2,2)
Y_WEIGHTS=(Q(0),Q(0),Q(1),Q(9,8),Q(7,5),Q(9,5))
M_WEIGHTS=(2,2,4,6,8,8);S_WEIGHTS=(0,0,4,4,4,8)
LEDGER=('export','food_U','food_W','drive_F','drive_P','wash_A','wash_B','wash_mass')


def dot(a,b):return sum(x*y for x,y in zip(a,b))
def number(ctx,x):
    x=Q(str(x));return ctx.mpf(x.numerator)/x.denominator
def logsum(ctx,values):
    values=list(values);shift=values[0]
    return shift+ctx.log(sum(ctx.exp(v-shift) for v in values))


@dataclass(frozen=True)
class Parameters:
    K: object=KINETIC_BIAS
    release: object=DUPLEX_RELEASE
    delta: object=DRIVEN_CLEAVAGE
    epsilon: object=BASAL_RATE
    def __post_init__(self):
        for name in ('K','release','delta','epsilon'):
            q=Q(str(getattr(self,name)))
            if q<=0:raise ValueError('Positive reaction coefficients required.')
            object.__setattr__(self,name,q)
    @property
    def eta(self):return self.epsilon/16
    def theorem_scope(self):
        return 8<=self.K<=12 and 18<=self.release<=22 and Q(1,100)<=self.delta<=Q(1,20) and self.epsilon==Q(1,500000000)


@dataclass(frozen=True)
class Channel:
    label: int
    name: str
    inputs: tuple[int,...]
    outputs: tuple[int,...]
    coefficient: Q
    marks: tuple[int,...]=(0,)*8
    @property
    def jump(self):return tuple(self.outputs.count(i)-self.inputs.count(i) for i in range(6))
    def count_rate(self,counts,V):
        rate=self.coefficient*Q(V)**(1-len(self.inputs));seen={}
        for i in self.inputs:
            rate*=counts[i]-seen.get(i,0);seen[i]=seen.get(i,0)+1
        return rate


class DrivenReactor:
    """Preserves original labels after deleting only channels 6 and 7."""
    def __init__(self,parameters=Parameters(),enabled=True):
        self.p=parameters;self.enabled=enabled;p=parameters
        pairs=(('basal',(0,1),(2,),p.epsilon,p.epsilon/p.K),
               ('bind_U',(2,0),(3,),Q(20),Q(20)),
               ('bind_W',(3,1),(4,),Q(20),Q(20)),
               ('ligation',(4,),(5,),Q(20),Q(20)/p.K),
               ('release',(5,),(2,2),p.release,p.release))
        channels=[]
        for j,(name,left,right,kp,km) in enumerate(pairs):
            channels.extend([Channel(2*j,name+'_forward',left,right,kp),Channel(2*j+1,name+'_reverse',right,left,km)])
        for i in (0,1):
            marks=[0]*8;marks[1+i]=1
            channels.append(Channel(10+i,'feed_'+SPECIES[i],(),(i,),Q(1),tuple(marks)))
        for i in range(6):
            marks=(S_WEIGHTS[i],0,0,0,0,A_WEIGHTS[i],B_WEIGHTS[i],M_WEIGHTS[i])
            channels.append(Channel(12+i,'wash_'+SPECIES[i],(i,),(),Q(1),marks))
        channels.extend([Channel(18,'driven_cleavage',(2,),(0,1),p.delta,(0,0,0,1,0,0,0,0)),
                         Channel(19,'driven_reverse',(0,1),(2,),p.delta*p.eta,(0,0,0,0,1,0,0,0))])
        self.channels=tuple(c for c in channels if enabled or c.label not in (6,7))
        self.by_label={c.label:c for c in self.channels}
        self.jumps=np.array([c.jump for c in self.channels],int);self.marks=np.array([c.marks for c in self.channels],int)
        self.coefficients=np.array([float(c.coefficient) for c in self.channels])
        self.orders=np.array([len(c.inputs) for c in self.channels])
        self.inputs=np.full((len(self.channels),2),6,int);self.offsets=np.zeros((len(self.channels),2),int)
        for j,c in enumerate(self.channels):
            seen={}
            for k,i in enumerate(c.inputs):
                self.inputs[j,k]=i;self.offsets[j,k]=seen.get(i,0);seen[i]=seen.get(i,0)+1

    def rates(self,state,V=None):
        factors=np.r_[state,1][self.inputs]
        if V is not None:factors=np.maximum(0,factors-self.offsets)
        rate=self.coefficients*np.prod(factors,axis=1)
        return rate if V is None else rate*np.power(float(V),1-self.orders)

    def generator(self,counts,V,weights):
        rates=[c.count_rate(counts,V) for c in self.channels]
        jumps=[dot(weights,c.jump) for c in self.channels]
        return sum(a*d for a,d in zip(rates,jumps)),sum(a*d*d for a,d in zip(rates,jumps))

    def deterministic(self,H,method='Radau'):
        H=float(Q(str(H)))
        if not 1<=H<=10000:raise ValueError('ODE illustration is limited to duration 1..10000; bounds have a separate interface.')
        def rhs(t,state):
            rates=self.rates(state[:6]);return np.r_[rates@self.jumps,rates@self.marks]
        def entry(t,state):return np.array(Y_WEIGHTS,float)@state[:6]-1/2500
        entry.direction=1
        sol=solve_ivp(rhs,(0,500+H),[1,1,0,0,0,0]+[0]*8,method=method,dense_output=True,events=entry,
                       rtol=1e-10,atol=1e-16,max_step=1)
        if not sol.success:raise ArithmeticError(sol.message)
        return sol

    def stochastic(self,policy,limit=SSA_EVENT_BUDGET,seed=SSA_SEED):
        if not 1<=policy.V<10**12 or limit<1:raise ValueError('Direct SSA needs V<10^12 and a positive event budget.')
        rng=np.random.default_rng(seed);monitor=OperatingMonitor(self,policy);events=0
        while monitor.status=='active' and events<limit:
            rates=self.rates(monitor.counts,policy.V);total=float(sum(rates))
            event_time=monitor.time+rng.exponential(1/total)
            boundary=float(monitor.next_boundary())
            if event_time>boundary:
                monitor.advance(boundary);continue
            j=int(np.searchsorted(np.cumsum(rates),rng.random()*total,side='right'))
            monitor.event(event_time,self.channels[j].label);events+=1
        report=monitor.report();report.update(events=events)
        if monitor.status=='active':report.update(status='unfinished',event_verdict=None)
        return report


@dataclass(frozen=True)
class ObservationPolicy:
    V: int=COPY_SCALE
    H: object=POST_STARTUP_DURATION
    def __post_init__(self):
        object.__setattr__(self,'H',Q(str(self.H)))
        if not isinstance(self.V,int) or self.V<1 or self.H<1:raise ValueError('Positive integer copy scale and H>=1 required.')
    @property
    def end(self):return 500+self.H
    @property
    def windows(self):return self.H.numerator//self.H.denominator
    @property
    def export_cap(self):return math.ceil(Q(self.V,5000))
    @property
    def supply_cap(self):return math.ceil(2*self.V*self.end)+1


class OperatingMonitor:
    """Tracks one continuing path. Boundaries reset only the window counter.

    Events exactly at a window's right endpoint belong to that window. Entry
    at time 500 is allowed. Return at Y<=V/5000 fails (including equality).
    Natural ledgers provide a second audit of the capped measurement counters.
    """
    def __init__(self,reactor,policy):
        if policy.windows>100000:raise ValueError('Monitor window limit; use ProbabilityBudget for huge horizons.')
        self.reactor=reactor;self.policy=policy;self.counts=[policy.V,policy.V,0,0,0,0]
        self.ledger=[0]*8;self.capped_supplies=[0]*4;self.counter=0
        self.time=0.;self.phase='pre_entry';self.entry_time=None;self.status='active'
        self.started=False;self.completed=0;self.all_windows_pass=True;self.window_exports=[]

    def next_boundary(self):
        if not self.started:return Q(500)
        return Q(501+self.completed) if self.completed<self.policy.windows else self.policy.end

    def corridor(self):
        V=self.policy.V
        return all(9*V<=10*dot(w,self.counts)<=11*V for w in (A_WEIGHTS,B_WEIGHTS))

    def _boundary(self):
        if not self.started:
            self.started=True
            if self.reactor.enabled and self.phase=='pre_entry':self.status='missed_deadline'
        elif self.completed<self.policy.windows:
            self.window_exports.append(self.counter)
            self.all_windows_pass &= self.counter>=self.policy.export_cap
            self.completed+=1;self.counter=0
        if self.status=='active' and self.time>=float(self.policy.end):self.status='complete'

    def advance(self,t):
        if t<self.time or t>float(self.policy.end):raise ValueError('Invalid monitor time.')
        while self.status=='active' and float(self.next_boundary())<=t:
            self.time=float(self.next_boundary());self._boundary()
        if self.status=='active':self.time=t

    def event(self,t,label):
        if self.status!='active' or t<self.time or t>float(self.policy.end):raise ValueError('Event outside active observation interval.')
        while self.status=='active' and float(self.next_boundary())<t:
            self.time=float(self.next_boundary());self._boundary()
        if self.status!='active':return
        channel=self.reactor.by_label[label]
        if channel.count_rate(self.counts,self.policy.V)<=0:raise ValueError('Reaction is not supported by molecule counts.')
        self.time=t;self.counts=[n+d for n,d in zip(self.counts,channel.jump)]
        if min(self.counts)<0:raise AssertionError('Count transition left the orthant.')
        self.ledger=[n+d for n,d in zip(self.ledger,channel.marks)]
        self.capped_supplies=[min(self.policy.supply_cap,n+d) for n,d in zip(self.capped_supplies,channel.marks[1:5])]
        if t>500 and self.completed<self.policy.windows:
            self.counter=min(self.policy.export_cap,self.counter+channel.marks[0])
        V=self.policy.V;Y40=dot((0,0,40,45,56,72),self.counts)
        if self.reactor.enabled:
            if self.phase=='pre_entry' and 2500*Y40>=40*V:self.phase='entered';self.entry_time=t
            elif self.phase=='entered' and 5000*Y40<=40*V:self.phase='returned';self.status='return_failure'
        if not self.corridor():self.status='resource_exit'
        self.audit_ledger()
        # Apply deterministic gates after a reaction on their exact endpoint.
        if self.status=='active' and float(self.next_boundary())==t:self._boundary()

    def audit_ledger(self):
        V=self.policy.V;E,IU,IW,QF,QP,OA,OB,MW=self.ledger
        if dot(A_WEIGHTS,self.counts)+OA!=V+IU or dot(B_WEIGHTS,self.counts)+OB!=V+IW:
            raise AssertionError('Moiety ledger failed.')
        if dot(M_WEIGHTS,self.counts)+MW!=4*V+2*(IU+IW):raise AssertionError('Mass ledger failed.')
        if OA*4<E or OB*4<E:raise AssertionError('Export exceeded washed moiety accounting.')
        if self.capped_supplies!=[min(self.policy.supply_cap,n) for n in self.ledger[1:5]]:
            raise AssertionError('Capped counters no longer match natural histories.')

    def report(self):
        V=self.policy.V;T=self.policy.end;IU,IW,QF,QP=self.capped_supplies
        supplies=IU<=2*V*T and IW<=2*V*T and QF+QP<=V*T/8
        verdict=None
        if self.status!='active':
            verdict=(self.status=='complete' and self.all_windows_pass and supplies) if self.reactor.enabled else (
                self.status=='resource_exit' or (self.status=='complete' and self.all_windows_pass))
        return {'status':self.status,'time':self.time,'phase':self.phase,'entry_time':self.entry_time,
            'counts':self.counts,'natural_ledger':dict(zip(LEDGER,self.ledger)),
            'capped_supply_counts':self.capped_supplies,'current_window_counter':self.counter,
            'completed_windows':self.completed,'capped_window_exports':self.window_exports,
            'all_windows_pass':bool(self.all_windows_pass),'supplies_pass':bool(supplies),'event_verdict':verdict,
            'comparison_scope':'Enabled joint event' if self.reactor.enabled else 'Disabled output schedule OR any resource exit credited as success; no entry/residence/supply requirement.'}