"""Explicit reporter storage with interchangeable native/regeneration rate laws."""
from dataclasses import dataclass
from typing import Protocol
import math
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import brentq


class RateLaw(Protocol):
    def __call__(self,concentration:float)->float:...
    def slopes(self,pool:float)->tuple[float,float]:...


@dataclass(frozen=True)
class LinearRate:
    coefficient:float
    def __post_init__(self):
        if not math.isfinite(self.coefficient) or self.coefficient<=0:raise ValueError('Positive finite rate required.')
    def __call__(self,x):return self.coefficient*x
    def slopes(self,pool):return self.coefficient,self.coefficient


@dataclass(frozen=True)
class MichaelisMenten:
    capacity:float
    affinity:float
    def __post_init__(self):
        if min(self.capacity,self.affinity)<=0 or not all(map(math.isfinite,(self.capacity,self.affinity))):raise ValueError('Positive finite rate parameters required.')
    def __call__(self,x):return self.capacity*x/(self.affinity+x)
    def slopes(self,pool):return self.capacity*self.affinity/(self.affinity+pool)**2,self.capacity/self.affinity


@dataclass(frozen=True)
class Reporter:
    total:float=1.
    release:float=20.
    dissociation:float=1.
    activity:float=.08
    def __post_init__(self):
        if min(self.total,self.release,self.dissociation,self.activity)<=0 or not all(map(math.isfinite,(self.total,self.release,self.dissociation,self.activity))):raise ValueError('Positive finite reporter parameters required.')
    @property
    def clearance(self):return self.release+self.dissociation
    @property
    def association(self):return self.activity*self.clearance/(self.total*self.release)


@dataclass(frozen=True)
class AssociationSchedule:
    acquisition:float=8.
    windows:tuple=((0.,8.),)
    fade_rate:float|None=None
    def __post_init__(self):
        if self.acquisition<=0 or not math.isfinite(self.acquisition):raise ValueError('Positive acquisition required.')
        previous=0.
        for l,r in self.windows:
            if not 0<=l<r<=self.acquisition or l<previous:raise ValueError('Ordered disjoint acquisition windows required.')
            previous=r
        if self.fade_rate is not None and (not math.isfinite(self.fade_rate) or self.fade_rate<=0):raise ValueError('Positive fade rate required.')
    def multiplier(self,t):
        if t>=self.acquisition:return 0. if self.fade_rate is None else math.exp(-self.fade_rate*(t-self.acquisition))
        return float(any(l<=t<r for l,r in self.windows))
    def breaks(self,end):return sorted(set([0.,end,self.acquisition]+[t for w in self.windows for t in w]))


class CofactorReactor:
    columns=('free','complex','reporter_product','native_deficit','reference_free','integrated_complex','association_mass','native_product')
    def __init__(self,regeneration:RateLaw,native:RateLaw,reporter:Reporter=Reporter(),pool=1.,background=0.):
        if pool<=0 or background<0 or not math.isfinite(pool+background):raise ValueError('Positive pool and nonnegative background required.')
        self.regeneration,self.native,self.reporter,self.pool,self.background=regeneration,native,reporter,pool,background
    def stationary(self):
        return brentq(lambda x:self.regeneration(self.pool-x)-self.native(x)-self.background*x,0,self.pool,xtol=1e-14)
    def run(self,schedule=AssociationSchedule(),recovery=5.,offset=0.,samples=201):
        if recovery<0 or not math.isfinite(recovery) or type(samples)!=int or samples<2:raise ValueError('Valid recovery and sampling budget required.')
        r=self.reporter;C=self.pool;xs=self.stationary();initial=xs*(1+offset)
        if not math.isfinite(offset) or not -1<offset or not 0<initial<=C:raise ValueError('Initial matched preparation outside positive physical pool.')
        y=np.array([initial,0,0,0,initial,0,0,0],float);times=[];states=[];at_acquisition=None
        end=schedule.acquisition+recovery;breaks=schedule.breaks(end)
        for left,right in zip(breaks,breaks[1:]):
            if right<=left:continue
            # Keep discontinuities out of an integration segment's interior.
            active=schedule.multiplier((left+right)/2) if right<=schedule.acquisition else None
            def rhs(t,y):
                x,b,q,D,x0,I,V,PN=y
                alpha=r.association*(active if active is not None else schedule.multiplier(t))
                v=alpha*x*(r.total-b);native=self.native(x);base=self.native(x0)
                return [self.regeneration(C-x-b)-native-self.background*x-v+r.dissociation*b,
                        v-r.clearance*b,r.release*b,base-native,
                        self.regeneration(C-x0)-base-self.background*x0,b,v,native]
            sol=solve_ivp(rhs,(left,right),y,method='DOP853',rtol=2e-11,atol=2e-13,dense_output=True)
            if not sol.success:raise RuntimeError(sol.message)
            grid=np.linspace(left,right,samples);values=sol.sol(grid).T
            times.extend(grid[1:] if times else grid);states.extend(values[1:] if states else values)
            y=sol.y[:,-1]
            if right==schedule.acquisition:at_acquisition=y.copy()
        t=np.array(times);Y=np.array(states);x,b,q,D,x0,I,V,PN=Y.T;e=x0-x
        if min(x.min(),b.min(),(C-x-b).min(),(r.total-b).min()) < -1e-8:raise RuntimeError('Numerical trajectory violates conserved inventories.')
        result=dict(acquisition_signal=float(at_acquisition[2]),integrated_loss=float(D[-1]),
                    peak_relative_suppression=float(max(e/x0)),recorded_native_product=float(PN[-1]),
                    final_free_deficit=float(e[-1]),final_complex=float(b[-1]),
                    complex_ledger_residual=float(max(abs(b-V+r.clearance*I))),
                    product_ledger_residual=float(max(abs(q-r.release*I))),
                    evidence='Numerical integration of literal mass-action source; not a uniform certificate.')
        if isinstance(self.regeneration,LinearRate) and isinstance(self.native,LinearRate):
            p,k=self.regeneration.coefficient,self.native.coefficient;a=p+k+self.background
            ratio=k*(p+r.release)/(a*r.release)
            residual=a*D-k*((p+r.release)*I+b-e)
            # Association still to come is zero for ideal shutoff. For a fade,
            # bound its entire tail; do not pretend the switch has already ended.
            Vtail=0. if schedule.fade_rate is None else r.association*C*r.total*math.exp(-schedule.fade_rate*recovery)/schedule.fade_rate
            qlo=q[-1]+r.release*b[-1]/r.clearance;qhi=qlo+r.release*Vtail/r.clearance
            Rlo=k/a*(e[-1]-b[-1]+(p+r.release)*b[-1]/r.clearance)
            Rhi=Rlo+k/a*(p+r.release)*Vtail/r.clearance
            result.update(cost_coefficient=ratio,finite_account_residual=float(max(abs(residual))),
                remaining_association_upper=Vtail,complete_product_enclosure_numeric=[float(qlo),float(qhi)],
                remaining_loss_enclosure_numeric=[float(Rlo),float(Rhi)],
                integrated_plus_tail_loss_enclosure_numeric=[float(D[-1]+Rlo),float(D[-1]+Rhi)],
                independent_integral_ratio_numeric=float((D[-1]+Rlo)/qlo),
                source_identity_loss_enclosure_numeric=[float(ratio*qlo),float(ratio*qhi)])
        return t,Y,result