"""Literal seven-species, fourteen-channel finite-fuel source and operations."""
from dataclasses import dataclass
from fractions import Fraction as F
from math import prod
import numpy as np

SPECIES=('X','Y','F','PX','PY','H','W')


@dataclass(frozen=True)
class Reaction:
    name:str
    reactants:tuple
    products:tuple
    rate:F
    def __post_init__(self):
        object.__setattr__(self,'rate',F(self.rate))
        if len(self.reactants)!=7 or len(self.products)!=7 or any(not isinstance(n,int) or n<0 for n in self.reactants+self.products) or self.rate<=0:raise ValueError('Positive rate and seven nonnegative integer stoichiometries required.')
    @property
    def change(self):return tuple(b-a for a,b in zip(self.reactants,self.products))
    def propensity(self,state):
        if any(n<a for n,a in zip(state,self.reactants)):return F(0)
        return self.rate*prod(prod(range(int(n)-a+1,int(n)+1)) for n,a in zip(state,self.reactants))


def stoich(**counts):return tuple(counts.get(name,0) for name in SPECIES)


@dataclass(frozen=True)
class Rates:
    correction:F=F(2*10**7)
    reverse_correction:F=F(6,10**11)
    leakage:F=F(1,10**7)
    replication:F=F(1)
    reverse_replication:F=F(1,100)
    production:F=F(1)
    reverse_production:F=F(1,100)
    def __post_init__(self):
        for name in self.__dataclass_fields__:
            object.__setattr__(self,name,F(getattr(self,name)))
            if getattr(self,name)<=0:raise ValueError('The reversible source requires strictly positive rates.')


class ChemicalSource:
    def __init__(self,rates=Rates()):
        self.rates=rates;reactions=[]
        def pair(name,a,b,forward,reverse):
            reactions.extend([Reaction(name,a,b,forward),Reaction(name+'_reverse',b,a,reverse)])
        for label in ['X','Y']:
            pair('replicate_'+label,stoich(**{label:1,'F':1}),stoich(**{label:2}),rates.replication,rates.reverse_replication)
            pair('produce_'+label,stoich(**{label:1,'F':1}),stoich(**{label:1,'P'+label:1}),rates.production,rates.reverse_production)
        pair('leak_X_to_Y',stoich(X=1),stoich(Y=1),rates.leakage,rates.leakage)
        pair('correct_X',stoich(X=2,Y=1,H=1),stoich(X=3,W=1),rates.correction,rates.reverse_correction)
        pair('correct_Y',stoich(Y=2,X=1,H=1),stoich(Y=3,W=1),rates.correction,rates.reverse_correction)
        self.reactions=tuple(reactions)
        if any(sum(r.change[:5]) or sum(r.change[5:]) for r in self.reactions):raise ArithmeticError('Reaction table violates material invariants.')
        self.changes=np.array([r.change for r in reactions],dtype=np.int64)
    def propensities(self,state):return tuple(r.propensity(state) for r in self.reactions)
    def equilibrium_activities(self):
        # Deterministic mass-action activities use powers, not count factorials.
        r=self.rates;resident=r.replication/r.reverse_replication;product=r.production/r.reverse_production
        a=(resident,resident,F(1),product,product,F(1),r.correction/r.reverse_correction)
        flux=[rx.rate*prod(v**n for v,n in zip(a,rx.reactants)) for rx in self.reactions]
        if any(flux[j]!=flux[j+1] for j in range(0,14,2)):raise ArithmeticError('Positive activities fail detailed balance.')
        return a
    def simulate(self,initial,horizon,rng,event_budget=200000,record=True):
        if len(initial)!=7 or any(not isinstance(n,(int,np.integer)) or n<0 for n in initial) or horizon<0:raise ValueError('Seven nonnegative counts and a nonnegative horizon required.')
        state=np.array(initial,dtype=np.int64);mass=(int(state[:5].sum()),int(state[5:].sum()));t=0.;trace=[(t,*map(int,state),-1)];fired=np.zeros(14,dtype=int)
        for count in range(event_budget):
            props=np.array([float(a) for a in self.propensities(state)]);rate=props.sum()
            if rate==0:break
            dt=rng.exponential(1/rate)
            if t+dt>horizon:break
            t+=dt;j=int(rng.choice(14,p=props/rate));state+=self.changes[j];fired[j]+=1
            if min(state)<0 or (int(state[:5].sum()),int(state[5:].sum()))!=mass:raise ArithmeticError('Invalid source jump.')
            if record:trace.append((t,*map(int,state),j))
        else:return dict(status='event_budget_exhausted',state=state,time=t,trace=trace,counts=fired)
        trace.append((float(horizon),*map(int,state),-1))
        return dict(status='completed',state=state,time=float(horizon),trace=trace,counts=fired)


@dataclass(frozen=True)
class BatchOperations:
    core_inventory:int=80
    fuel_inventory:int=1
    minimum_residents:int=8
    minority_allowance:int=1
    recovery:F=F(9,10)
    partition:F=F(12,25)
    def __post_init__(self):
        for name in ['core_inventory','fuel_inventory','minimum_residents','minority_allowance']:
            value=getattr(self,name)
            if not isinstance(value,int) or value<0:raise ValueError('Nonnegative integer inventories required.')
        if self.core_inventory<2*self.minimum_residents or self.minimum_residents<=self.minority_allowance:raise ValueError('Need disjoint restart regions and room for both daughters.')
        for name in ['recovery','partition']:
            object.__setattr__(self,name,F(getattr(self,name)))
            if not 0<=getattr(self,name)<=1:raise ValueError('Operation probability outside [0,1].')
    def region(self,state,label):
        if label not in ['X','Y']:raise ValueError('Unknown program.')
        if len(state)!=7 or any(n<0 for n in state):return False
        i=0 if label=='X' else 1;j=1-i
        return bool(sum(state[:5])==self.core_inventory and state[i]>=self.minimum_residents and state[j]<=self.minority_allowance and state[3]==state[4]==state[6]==0 and state[5]==self.fuel_inventory)
    def apply(self,state,rng):
        if len(state)!=7 or any(not isinstance(n,(int,np.integer)) or n<0 for n in state):raise ValueError('Seven nonnegative integer counts required.')
        mother=np.array(state,dtype=np.int64)
        if min(mother)<0 or int(mother[:5].sum())!=self.core_inventory or int(mother[5:].sum())!=self.fuel_inventory:raise ValueError('Mother inventories do not match the operation contract.')
        product=mother[3:5].copy();credited=rng.binomial(product,float(self.recovery));mother[3:5]=0
        a=rng.binomial(mother,float(self.partition));b=mother-a;before=(a.copy(),b.copy())
        food=fuel=waste=0
        for daughter in [a,b]:
            added=self.core_inventory-int(daughter[:5].sum());food+=added;daughter[2]+=added
            waste+=int(daughter[6]);daughter[6]=0;fuel+=self.fuel_inventory-int(daughter[5]);daughter[5]=self.fuel_inventory
        if food!=self.core_inventory+int(product.sum()):raise ArithmeticError('Refill ledger mismatch.')
        return dict(credited=credited,harvested=product,unrecovered=product-credited,partition_before_refill=before,daughters=(a,b),food_supplied=food,fuel_supplied=fuel,waste_removed=waste)
    def score(self,result,label):
        if label not in ['X','Y']:raise ValueError('Unknown program.')
        i=0 if label=='X' else 1
        return bool(result['credited'][i]>=4 and result['credited'][1-i]<=1 and all(self.region(d,label) for d in result['daughters']))


class ProductRetention:
    def __init__(self,selected_product='X',gate=4,background=F(1,2)):
        if selected_product not in ['X','Y'] or not isinstance(gate,int) or gate<0 or not 0<=background<=1:raise ValueError('Invalid product retention rule.')
        self.index=0 if selected_product=='X' else 1;self.gate=gate;self.background=F(background)
    def retain(self,operation_result,rng):
        # Product only: no internal label is inspected.
        probability=1. if operation_result['credited'][self.index]>=self.gate else float(self.background)
        return [d for d in operation_result['daughters'] if rng.random()<probability]