"""Exact source envelopes, observation inversion, and finite recovery bounds."""
from dataclasses import dataclass
from fractions import Fraction as F


def interval(a,b):
    a,b=F(a),F(b)
    if not 0<a<=b:raise ValueError('Positive ordered interval required.')
    return a,b


def exp_bounds(x):
    # Alternating Taylor enclosure with outward decimal rounding and squaring.
    # Adapted from the separately authored manuscript-60 companion.
    x=F(x)
    if x<0:raise ValueError('Nonnegative exponent magnitude required.')
    folds=0
    while x>1:x/=2;folds+=1
    term=total=F(1)
    for n in range(1,121):term*=-x/n;total+=term
    lo,hi=max(F(0),total-term*x/121),min(F(1),total);scale=10**40
    def outward(l,u):return F((l*scale).__floor__(),scale),F((u*scale).__ceil__(),scale)
    lo,hi=outward(lo,hi)
    for _ in range(folds):lo,hi=outward(lo*lo,hi*hi)
    return lo,hi


@dataclass(frozen=True)
class DesignBox:
    low:tuple=(F(95,100),F(105,100))
    high:tuple=(F(195,100),F(205,100))
    native:tuple=(F(98,100),F(102,100))
    pool:tuple=(F(99,100),F(101,100))
    reporter:tuple=(F(99,100),F(101,100))
    activity:tuple=(F(79,1000),F(81,1000))
    release:tuple=(F(20),F(25))
    dissociation:tuple=(F(1),F(2))
    gain:tuple=(F(98,100),F(102,100))
    preparation:F=F(1,100)
    error:F=F(1,100)
    def __post_init__(self):
        for name in ['low','high','native','pool','reporter','activity','release','dissociation','gain']:interval(*getattr(self,name))
        if self.low[1]>=self.high[0] or not 0<=self.preparation<1 or self.error<0:raise ValueError('Separated classes, preparation <1 and nonnegative error required.')
    def constants(self):
        p0,p1=self.low[0],self.high[1];k0,k1=self.native;c0,c1=self.release;u0,u1=self.activity
        amin=p0+k0;hmin=self.dissociation[0]+c0;hmax=self.dissociation[1]+c1
        B=u1*(1+p1/c0)/amin;m=1-self.preparation-(1+self.preparation)*B
        xb=p1*self.pool[1]/(p1+k0)
        return dict(a=amin,h=hmin,B=B,m=m,sigma=u1*self.pool[1]/(c0*self.reporter[0]),
            suppression=(1+self.preparation)/(1-self.preparation)*B,
            e=(1+self.preparation)*xb*B,b=u1*(1+self.preparation)*xb/c0,
            A=u1*hmax*self.pool[1]/c0,d=c1-k0,ka=k1/amin,pch=(p1+c1)/hmin)

    def contract(self,T=F(8),rounded=False,allowance=F(1,20)):
        return ObservationContract(self,F(T),rounded,F(allowance))


class ObservationContract:
    def __init__(self,box,T,rounded,allowance):
        self.box,self.T,self.rounded,self.allowance=box,T,rounded,allowance;c=box.constants()
        if box.release[0]<box.native[1] or c['m']<=0 or T<1/c['h']:raise ValueError('Signal/preservation envelope assumptions fail.')
        m,sigma=c['m'],c['sigma']
        if rounded:
            if m<F(94,100) or sigma>F(5,1000):raise ValueError('Rounded manuscript constants not valid for this box.')
            m,sigma=F(94,100),F(5,1000)
        self.AL=box.gain[0]*box.activity[0]*m*box.pool[0]/(1+sigma)*(T-1/c['h'])
        self.AU=box.gain[1]*box.activity[1]*(1+box.preparation)*box.pool[1]*T
        self.U=self.AU*box.low[1]/(box.low[1]+box.native[0])
        self.L=self.AL*box.high[0]/(box.high[0]+box.native[1])
        self.threshold=(self.U+self.L)/2;self.margin=self.L-self.U-2*box.error
        self.certified=self.margin>0 and c['suppression']<=allowance

    def inverse(self,y):
        y=F(y);b=self.box;lo,hi=b.low[0],b.high[1];lower=y-b.error;upper=y+b.error
        if upper<0 or lower>=self.AU:return dict(status='incompatible',p=None,flux=None)
        if lower>0:lo=max(lo,b.native[0]*lower/(self.AU-lower))
        if upper<self.AL:hi=min(hi,b.native[1]*upper/(self.AL-upper))
        if lo>hi:return dict(status='incompatible',p=None,flux=None)
        flux=(b.native[0]*b.pool[0]*lo/(lo+b.native[0]),b.native[1]*b.pool[1]*hi/(hi+b.native[1]))
        return dict(status='outer-enclosure',p=(lo,hi),flux=flux)

    def classify(self,y,promise=False):
        outer=self.inverse(y)
        if outer['p'] is None:return 'incompatible'
        if not promise or not self.certified:return 'unresolved'
        lo,hi=outer['p'];b=self.box
        low=max(lo,b.low[0])<=min(hi,b.low[1]);high=max(lo,b.high[0])<=min(hi,b.high[1])
        if not low and not high:return 'incompatible with two-class promise'
        if low and not high:return 'low'
        if high and not low:return 'high'
        return 'unresolved'


def convolution_denominator(r,rho,nu):
    candidates=[]
    if 0<nu<=r and nu<rho:candidates.append(rho-nu)
    if 0<nu<=rho and nu<r:candidates.append(r-nu)
    if not candidates:raise ValueError('Inadmissible convolution decay rate.')
    return max(candidates)


class RecoveryCertificate:
    def __init__(self,box=DesignBox()):
        self.box=box;self.c=box.constants()
        if box.release[0]<box.native[1] or self.c['m']<=0:raise ValueError('Preservation assumptions fail.')
    def ideal(self,tau):
        tau=F(tau);c=self.c
        if c['h']<=c['a']:raise ValueError('This ideal envelope needs h_min > a_min.')
        B=c['e']+c['d']*c['b']/(c['h']-c['a'])
        return c['ka']*(B*exp_bounds(c['a']*tau)[1]+c['pch']*c['b']*exp_bounds(c['h']*tau)[1])
    def fading(self,gamma,tau,nu0,nu1):
        gamma,tau,nu0,nu1=map(F,(gamma,tau,nu0,nu1));c=self.c
        if gamma<=0 or tau<0:raise ValueError('Positive fade rate and nonnegative wait required.')
        d1=convolution_denominator(c['h'],gamma,nu1);d0=convolution_denominator(c['a'],nu1,nu0)
        beta=c['b']+c['A']/d1;B=c['e']+c['d']*beta/d0
        lower=upper=F(0)
        for factor,rate in [(c['ka']*B,nu0),(c['ka']*c['pch']*beta,nu1),(c['ka']*c['pch']*c['A']/gamma,gamma)]:
            l,u=exp_bounds(rate*tau);lower+=factor*l;upper+=factor*u
        return dict(envelope_lower=lower,loss_upper=upper,beta=beta,B=B,
                    extra_complete_loss_upper=c['ka']*(1+self.box.high[1]/self.box.release[0])*self.box.activity[1]*self.box.pool[1]/gamma)


def slope_sandwich(glo,ghi,klo,khi,c):
    glo,ghi,klo,khi,c=map(F,(glo,ghi,klo,khi,c))
    if not 0<=glo<=ghi or not 0<klo<=khi<=c:raise ValueError('Divided-slope ordering and c >= K+ required.')
    return klo*(ghi+c)/(c*(ghi+khi)),khi*(glo+c)/(c*(glo+klo))