"""Literal enzyme kinetics, reciprocal observation certificates and scalar tasks."""
from dataclasses import dataclass
from fractions import Fraction as F
import math
import numpy as np
from scipy.integrate import solve_ivp, quad
from scipy.optimize import linprog, brentq
import mpmath as mp


def rational(x):return x if isinstance(x,F) else F(str(x))


@dataclass(frozen=True)
class EnzymeParameters:
    capacity: F = F(1)
    substrate_affinity: F = F(7)
    carrier_affinity: F = F(3)
    inhibition_H: F = F(56)
    inhibition_A: F = F(125)
    inhibition_B: F = F(520)

    def __post_init__(self):
        for name in self.__dataclass_fields__:
            v=rational(getattr(self,name));object.__setattr__(self,name,v)
            if v<=0:raise ValueError('All kinetic parameters must be positive.')

    def rate(self,N,S,H=0,A=0,B=0):
        if min(N,S,H,A,B)<0:raise ValueError('Nonnegative concentrations required.')
        V,Kg,Kn,Kh,Ka,Kb=self.values()
        return V*(N/Kn)*(S/Kg)/(1+(N/Kn)*(1+S/Kg)+H/Kh+A/Ka+B/Kb)

    def values(self):return tuple(getattr(self,n) for n in self.__dataclass_fields__)

    def coefficients(self):
        V,Kg,Kn,Kh,Ka,Kb=self.values()
        return (1/V,Kg/V,Kn*Kg/V,Kn*Kg/(V*Kh),Kn*Kg/(V*Ka),Kn*Kg/(V*Kb))

    @classmethod
    def from_coefficients(cls,beta):
        b=tuple(map(rational,beta))
        if len(b)!=6 or min(b)<=0:raise ValueError('Six strictly positive reciprocal coefficients required.')
        return cls(1/b[0],b[1]/b[0],b[2]/b[1],b[2]/b[3],b[2]/b[4],b[2]/b[5])


def features(N,S,H=0,A=0,B=0):
    N,S,H,A,B=map(rational,[N,S,H,A,B])
    if min(N,S)<=0 or min(H,A,B)<0:raise ValueError('Positive substrates and nonnegative inhibitors required.')
    return (F(1),1/S,1/(N*S),H/(N*S),A/(N*S),B/(N*S))


class ObservationPolyhedron:
    """Individual-realization class beta>0, rows*beta<=bounds.

    Solver output only proposes dual multipliers. Exact checking decides whether
    a proposed finite certificate works; nonemptiness needs a positive witness.
    """
    def __init__(self,rows,bounds):
        self.rows=tuple(tuple(map(rational,r)) for r in rows);self.bounds=tuple(map(rational,bounds))
        if not self.rows or len(self.rows)!=len(self.bounds) or any(len(r)!=6 for r in self.rows):raise ValueError('Six-column constraint rows required.')

    @classmethod
    def from_rate_bands(cls,observations):
        rows=[];bounds=[]
        for inputs,lo,hi in observations:
            lo,hi=map(rational,[lo,hi])
            if not 0<lo<=hi:raise ValueError('Strictly positive ordered rate band required.')
            phi=features(*inputs);rows.extend([phi,tuple(-v for v in phi)]);bounds.extend([1/lo,-1/hi])
        return cls(rows,bounds)

    def contains(self,beta):
        beta=tuple(map(rational,beta))
        return len(beta)==6 and min(beta)>0 and all(sum(a*b for a,b in zip(row,beta))<=bound for row,bound in zip(self.rows,self.bounds))

    def verify_dual(self,target,multipliers,witness=None):
        target=tuple(map(rational,target));y=tuple(map(rational,multipliers))
        if len(target)!=6 or len(y)!=len(self.rows) or min(y)<0:raise ValueError('Invalid nonnegative dual vector.')
        if any(sum(y[i]*self.rows[i][j] for i in range(len(y)))!=target[j] for j in range(6)):raise ValueError('Dual row identity does not hold exactly.')
        U=sum(a*b for a,b in zip(y,self.bounds))
        if U<=0:raise ValueError('Positive target reciprocal bound required.')
        if witness is not None and not self.contains(witness):raise ValueError('Supplied positive class witness is incompatible.')
        return dict(reciprocal_upper=U,rate_lower=1/U,multipliers=y,nonempty_verified=witness is not None)

    def propose_dual(self,target,witness=None,denominator_budget=10**6):
        target=tuple(map(rational,target))
        result=linprog(np.array(self.bounds,float),A_eq=np.array(self.rows,float).T,b_eq=np.array(target,float),bounds=(0,None),method='highs')
        if not result.success:return dict(status='not certified',reason=result.message)
        y=[F(float(v)).limit_denominator(denominator_budget) for v in result.x]
        try:return dict(status='exact certificate verified',**self.verify_dual(target,y,witness))
        except ValueError as e:return dict(status='not certified',reason=str(e))


@dataclass(frozen=True)
class ScalarTask:
    parameters: EnzymeParameters = EnzymeParameters()
    pool: F = F(56)
    substrate: F = F(7)
    inhibitor_A: F = F(0)
    inhibitor_B: F = F(0)
    load: F = F(3,10)
    initial: F = F(10)
    target: F = F(28)
    deadline: F = F(120)

    def __post_init__(self):
        for n in list(self.__dataclass_fields__)[1:]:object.__setattr__(self,n,rational(getattr(self,n)))
        if not 0<=self.initial<self.target<self.pool or min(self.substrate,self.load,self.deadline)<=0 or min(self.inhibitor_A,self.inhibitor_B)<0:raise ValueError('Invalid scalar task.')

    def collapsed(self):
        b=self.parameters.coefficients();S=self.substrate
        return b[0]+b[1]/S,b[3]/S,(b[2]+b[4]*self.inhibitor_A+b[5]*self.inhibitor_B)/S

    def rate(self,g):
        a,b,c=self.collapsed();return (self.pool-g)/(a*(self.pool-g)+b*g+c)

    def drift(self,g):return self.rate(g)-self.load

    def target_features(self):return features(self.pool-self.target,self.substrate,self.target,self.inhibitor_A,self.inhibitor_B)

    def uniform_deadline(self,reciprocal_upper):
        U=rational(reciprocal_upper)
        if U<=0:raise ValueError('Positive reciprocal upper bound required.')
        margin=1/U-self.load
        return dict(status='certified' if margin>0 else 'not certified',drift_lower=margin,
                    arrival_upper=(self.target-self.initial)/margin if margin>0 else None)

    def passage_rectangles(self,n=128):
        if type(n)!=int or not 1<=n<=65536:raise ValueError('Rectangle budget: 1..65536.')
        if self.drift(self.target)<=0:return dict(status='no finite arrival',lower=None,upper=None)
        width=(self.target-self.initial)/n
        values=[1/self.drift(self.initial+j*width) for j in range(n+1)]
        return dict(status='finite arrival',lower=width*sum(values[:-1]),upper=width*sum(values[1:]))

    def classify(self,n=128):
        if self.drift(0)<0:return dict(status='not certified',reason='Nonnegative-pool invariance is not established by the inward-field contract.')
        r=self.passage_rectangles(n)
        if r['status']=='no finite arrival':return dict(status='failure',reason='nonpositive target drift; monotonicity and uniqueness prevent finite arrival')
        r.pop('status')
        if r['upper']<=self.deadline:return dict(status='success',reason='exact upper passage bound',**r)
        if r['lower']>self.deadline:return dict(status='failure',reason='exact lower passage bound',**r)
        return dict(status='unresolved',reason='deadline overlaps passage enclosure',**r)

    def passage_numeric(self):
        if self.drift(self.target)<=0:return math.inf
        # Closed-form integral in high precision prevents near-threshold
        # cancellation. This remains a numerical value, not an enclosure.
        with mp.workdps(50):
            cv=lambda x:mp.mpf(x.numerator)/x.denominator
            a,b,c=map(cv,self.collapsed());P,x,R,q=map(cv,[self.pool,self.initial,self.target,self.load])
            A=a*P+c;B=b-a;u=P-q*A;v=1+q*B
            if v==0:result=(A*(R-x)+B*(R*R-x*x)/2)/u
            else:result=-B*(R-x)/v+(A+B*u/v)/v*mp.log((u-v*x)/(u-v*R))
            return float(result)

    def trajectory(self,points=241):
        if self.drift(0)<0:raise ValueError('Inward field at zero is required for the physical pool trajectory.')
        P,S,A,B=map(float,[self.pool,self.substrate,self.inhibitor_A,self.inhibitor_B]);q=float(self.load)
        parameters=self.parameters
        def rhs(t,state):
            # Solver stages may stray slightly outside the invariant interval.
            # The algebraic law is continued smoothly here, without clipping the state.
            g=state[0];V,Kg,Kn,Kh,Ka,Kb=map(float,parameters.values());N=P-g
            return [V*(N/Kn)*(S/Kg)/(1+(N/Kn)*(1+S/Kg)+g/Kh+A/Ka+B/Kb)-q]
        def hit(t,state):return state[0]-float(self.target)
        hit.direction=1
        sol=solve_ivp(rhs,(0,float(self.deadline)),[float(self.initial)],method='Radau',rtol=2e-10,atol=2e-11,
                      t_eval=np.linspace(0,float(self.deadline),points),events=hit)
        if not sol.success:raise ArithmeticError(sol.message)
        return sol.t,sol.y[0],float(sol.t_events[0][0]) if len(sol.t_events[0]) else None


class CapacityFamily:
    """Only V varies; fixed affinities and task define a numerical cutoff."""
    def __init__(self,task=ScalarTask()):self.task=task

    def at(self,V):
        from dataclasses import replace
        return replace(self.task,parameters=replace(self.task.parameters,capacity=rational(V)))

    def eventual_threshold(self):return self.task.load/self.at(1).rate(self.task.target)

    def numerical_cutoff(self,deadline=None):
        T=float(self.task.deadline if deadline is None else deadline)
        if not math.isfinite(T) or T<=0:raise ValueError('Positive finite deadline required.')
        lower=float(self.eventual_threshold());lo=np.nextafter(lower,math.inf);hi=max(1.,2*lo)
        while self.at(hi).passage_numeric()>T:
            hi*=2
            if hi>1e12:raise ArithmeticError('Cutoff bracket budget exceeded.')
        if self.at(lo).passage_numeric()<T:raise ArithmeticError('Cutoff below floating-point resolution; use high precision.')
        return brentq(lambda V:self.at(V).passage_numeric()-T,lo,hi,xtol=1e-14)