"""Exact prescribed-state design, independent kinetics, and physical root audit."""
from dataclasses import dataclass
from fractions import Fraction as Q
import sympy as sp
import phos_sharp as ps

@dataclass(frozen=True)
class PrescribedStates:
    auxiliary_roots:tuple
    enzyme_total_ratio:Q
    def construct(self):
        xs=tuple(Q(str(x)) for x in self.auxiliary_roots);r=Q(str(self.enzyme_total_ratio))
        if not xs or len(xs)%2!=1 or len(set(xs))!=len(xs) or min(xs)<=1 or r<=max(map(ps.ratio,xs)):raise ValueError('Distinct odd-length roots >1, with r>max (x*x-1)/8 required.')
        rec=ps.build(sorted(xs),r)
        if not rec['positive']:raise ValueError('Coefficient positivity failed: enlarge the enzyme ratio or change the prescribed roots.')
        return rec

@dataclass(frozen=True)
class KineticFreedom:
    kinase_unbinding:tuple
    phosphatase_unbinding:tuple
    phosphatase_catalysis:tuple
    def realize(self,rec):
        n=len(rec['B']);b,be,ga=[tuple(Q(str(v)) for v in a) for a in (self.kinase_unbinding,self.phosphatase_unbinding,self.phosphatase_catalysis)]
        if any(len(a)!=n or min(a)<=0 for a in [b,be,ga]):raise ValueError('Three positive site arrays are required.')
        out=[]
        for i in range(n):
            p=rec['B'][i]/rec['A'][i];q=rec['D'][i]/rec['A'][i+1];c=ga[i]*rec['D'][i]/rec['B'][i]
            out.append(((b[i]+c)*p,b[i],c,(be[i]+ga[i])*q,be[i],ga[i]))
        return out

class EquilibriumChart:
    """General positive six-rate kinetics. Exact real-root isolation, not a grid.

    Clearing denominators creates candidates. Positivity, poles and u=r are
    checked separately. Algebraic states are retained symbolically.
    """
    def __init__(self,rates,totals):
        self.rates=[tuple(Q(str(v)) for v in row) for row in rates];self.totals=tuple(Q(str(v)) for v in totals)
        if not self.rates or any(len(row)!=6 or min(row)<=0 for row in self.rates) or len(self.totals)!=3 or min(self.totals)<=0:raise ValueError('Positive six-rate rows and positive ET,FT,ST required.')
        A=[Q(1)];B=[];D=[]
        for a,b,c,al,be,ga in self.rates:
            p=a/(b+c);q=al/(be+ga);B.append(p*A[-1]);A.append(A[-1]*c*p/(ga*q));D.append(q*A[-1])
        self.A,self.B,self.D=A,B,D;self.r=self.totals[0]/self.totals[1]
    def audit(self):
        u=sp.Symbol('u');poly=lambda a:sum(sp.Rational(v)*u**i for i,v in enumerate(a));A,B,D=map(poly,[self.A,self.B,self.D]);r,FT,ST=map(sp.Rational,[self.r,self.totals[1],self.totals[2]]);L=B-r*D;M=B-u*D
        raw=sp.Poly(sp.expand((r-u)*A*M+FT*(r-u)*(B+D)*u*L-ST*u*L*M),u)
        roots=sp.polys.polytools.intervals(raw,eps=sp.Rational(1,10**20));rows=[]
        # Evaluate signs at an exact algebraic root; SymPy's real algebraic
        # isolation decides these comparisons without rounding a near-zero value.
        for z,multiplicity in sp.polys.polytools.real_roots(raw,multiple=False):
            lv=L.subs(u,z);mv=M.subs(u,z)
            if z<=0:status='rejected_nonpositive_ratio'
            elif z==r:status='exceptional_ratio_checked_separately'
            elif lv==0 or mv==0:status='rejected_pole'
            elif (r-z)*lv<=0:status='rejected_negative_concentration'
            else:status='positive_regular_equilibrium'
            row=dict(ratio=str(z),ratio_numerical=float(z.evalf()),multiplicity=multiplicity,status=status)
            if status=='positive_regular_equilibrium':
                s=(r-z)/(z*lv);f=FT*lv/mv
                # All positive coefficients and admitted domain give positivity.
                row['state']=[str(v) for v in ps.state_from(self.A,self.B,self.D,z,s,f)]
            rows.append(row)
        exceptional=dict(status='absent',reason='B(r) differs from r D(r)')
        if L.subs(u,r)==0:
            # At u=r, enzyme totals are compatible for every s>0; substrate
            # total is strictly increasing in s and selects one positive s.
            s=sp.Symbol('s');ar=A.subs(u,r);dr=D.subs(u,r);br=B.subs(u,r)
            eq=sp.Poly(sp.expand(s*ar*(1+s*r*dr)+r*s*FT*(br+dr)-ST*(1+s*r*dr)),s)
            sr=[x for x in sp.real_roots(eq) if x>0]
            if len(sr)!=1:raise ArithmeticError('Exceptional positive root not unique.')
            sf=sr[0];f=FT/(1+sf*r*dr);state=ps.state_from(self.A,self.B,self.D,r,sf,f)
            exceptional=dict(status='positive_exceptional_equilibrium',ratio=str(r),free_substrate=str(sf),state=list(map(str,state)))
        return dict(raw_polynomial=str(raw.as_expr()),degree=raw.degree(),isolating_intervals=[dict(lower=str(a),upper=str(b),multiplicity=m) for (a,b),m in roots],candidates=rows,exceptional=exceptional,positive_count=sum(x['status']=='positive_regular_equilibrium' for x in rows)+int(exceptional['status']=='positive_exceptional_equilibrium'))

def stability(rates,state):
    mat=ps.reduced_jacobian(rates,state);cp=ps.charpoly(mat)
    try:count,signs=ps.routh_rhp(cp)
    except AssertionError:return dict(status='unresolved_regular_routh')
    return dict(status='certified_regular_routh',unstable=count,first_column_signs=signs,determinant=str(ps.det(mat)))

def substrate_window(rec,probes,half_width):
    r=rec['r'];ET,FT,ST=rec['totals'];probes=list(map(Q,probes));half_width=Q(half_width)
    if half_width<=0 or probes!=sorted(set(probes)) or not 0<probes[0]<probes[-1]<r:raise ValueError('Ordered interior probes and positive width required.')
    u=sp.Symbol('u');L=ps.add(rec['B'],ps.scale(-r,rec['D']));lp=sp.Poly(sum(sp.Rational(v)*u**i for i,v in enumerate(L)),u)
    continuous=lp.count_roots(probes[0],probes[-1])==0 and ps.val(L,probes[0])>0
    values=[ps.chart(rec['A'],rec['B'],rec['D'],r,FT,p)[0]-ST for p in probes]
    passes=continuous and all((-1)**j*v>half_width for j,v in enumerate(values))
    return dict(status='certified_at_least_'+str(len(probes)-1) if passes else 'not_certified',continuous=continuous,half_width=str(half_width),probe_residuals=list(map(str,values)),scope='Continuity plus strict alternating signs; equality with 2n-1 also invokes the published universal upper bound.')