"""Four resident species with literal extraction and retained molecular identity."""
from dataclasses import dataclass,replace
from fractions import Fraction as F
import math
import numpy as np
from scipy.optimize import brentq

E=F(1,100000);ETA=F(1,10000)
ENERGY={
 'L':np.array([[1079122,-568413,2227359,3266330],[-568413,1057785,-2194388,-3215160],
  [2227359,-2194388,6383792,9601239],[3266330,-3215160,9601239,14651126]],dtype=object),
 'H':np.array([[976196,-1500816,2954639,4411258],[-1500816,5195401,-8320868,-12539182],
  [2954639,-8320868,14137961,21329594],[4411258,-12539182,21329594,32242779]],dtype=object)}


@dataclass(frozen=True)
class Cell:
    tag:str
    size:int
    counts:tuple
    def __post_init__(self):
        if self.tag not in ('L','H') or type(self.size) is not int or self.size<1 or len(self.counts)!=4 or any(type(n) is not int or n<0 for n in self.counts):
            raise ValueError('Cell requires an immutable tag, positive size and four nonnegative integer counts')


@dataclass(frozen=True)
class Channel:
    name:str
    inputs:tuple
    outputs:tuple
    coefficient:F
    def rate(self,cell):
        a=self.coefficient*F(cell.size)**(1-sum(self.inputs))
        for n,k in zip(cell.counts,self.inputs):
            if n<k:return F(0)
            for offset in range(k):a*=n-offset
        return a
    def apply(self,cell):
        if self.rate(cell)<=0:raise ValueError('Disabled chemical event')
        return replace(cell,counts=tuple(n-a+b for n,a,b in zip(cell.counts,self.inputs,self.outputs)))


def reconstruction(z,rho=F(1,100)):
    B=60/(z+2);H=(16*z+2*z*z)/(2+ETA)
    K=(2*(1+2*ETA)*z*z-16*(1-ETA)*z)/(2+ETA)
    return (z*B+K+rho*z,B,z,H)


def residual(z,rho=F(1,100)):
    A,B,_,_=reconstruction(z,rho)
    return A+B+E*(A*A-B)-33


class ResidentChemistry:
    def __init__(self,rho='1/100'):
        self.rho=F(rho)
        if self.rho<0:raise ValueError('Nonnegative extraction coefficient required')
        pairs=[('A/Bz',(1,0,0,0),(0,1,1,0),F(1),F(1)),
            ('z/H',(0,0,1,0),(0,0,0,1),F(16),F(1)),
            ('H/2z',(0,0,0,1),(0,0,2,0),F(1),F(2)),
            ('A supply',(0,0,0,0),(1,0,0,0),F(6),F(1)),
            ('B supply',(0,0,0,0),(0,1,0,0),F(27),F(1)),
            ('B/2A',(0,1,0,0),(2,0,0,0),E,E)]
        channels=[]
        for name,a,b,kf,kr in pairs:channels += [Channel(name+'+',a,b,kf),Channel(name+'-',b,a,kr)]
        channels += [Channel('H waste',(0,0,0,1),(0,0,0,0),ETA),Channel('z extraction',(0,0,1,0),(0,0,0,0),self.rho)]
        self.channels=tuple(channels)
    def rates(self,cell):
        a,b,z,h=cell.counts;m=cell.size
        return np.array([a,b*z/m,16*z,h,h,2*z*(z-1)/m,6*m,a,27*m,b,b/100000,a*(a-1)/(100000*m),h/10000,float(self.rho)*z])
    def density(self,x,rho=None):
        a,b,z,h=x;rho=float(self.rho) if rho is None else rho
        return np.array([6-2*a+b*z+2e-5*(b-a*a),27+a-(1+z)*b-1e-5*(b-a*a),
            a-b*z-16*z-4*z*z+3*h-rho*z,16*z+2*z*z-2.0001*h])
    def centers(self):
        rho=float(self.rho);roots=[];grid=np.linspace(.01,6,1200)
        for a,b in zip(grid,grid[1:]):
            if residual(a,rho)*residual(b,rho)<0:
                z=brentq(lambda z:float(residual(z,rho)),a,b,xtol=1e-14);roots.append(np.array(reconstruction(z,rho),float))
        return roots
    def generator(self,cell,observable):
        return sum(c.rate(cell)*(observable(c.apply(cell))-observable(cell)) for c in self.channels if c.rate(cell)>0)


def rational_root(rho,tag,steps=140):
    """Strict paper bracket with rational bisection for a fixed admitted rho."""
    rho=F(rho)
    if not F(9999,10**6)<=rho<=F(1,100):raise ValueError('Outside certified extraction interval')
    lo,hi=map(F,('0.98172','0.98174') if tag=='L' else ('2.89014','2.89017'))
    fl=residual(lo,rho)
    if fl*residual(hi,rho)>=0:raise ArithmeticError('Invalid root bracket')
    for _ in range(steps):
        mid=(lo+hi)/2;fm=residual(mid,rho)
        if fm==0:return mid,mid
        if fl*fm<0:hi=mid
        else:lo=mid;fl=fm
    return lo,hi