"""Delivered shared input and exact source-to-mission certificate assembly."""
from dataclasses import dataclass
from fractions import Fraction as F
from math import factorial
import numpy as np
from scipy.integrate import solve_ivp
from population import InheritedPopulation,MolecularRates,Action


def exp_lower(x,degree):
    x=F(x)
    if x<0:raise ValueError('Positive-series lower bound requires x>=0.')
    return sum((x**k/factorial(k) for k in range(degree+1)),F(0))


@dataclass(frozen=True)
class AdministrationCourse:
    rate:F=F(29,100)
    stop:F=F(112)
    horizon:F=F(120)
    clearance_min:F=F(99,100)
    clearance_max:F=F(101,100)
    def __post_init__(self):
        for name in self.__dataclass_fields__:object.__setattr__(self,name,F(getattr(self,name)))
        if not self.rate>=0 or not 0<self.stop<=self.horizon or not 0<self.clearance_min<=self.clearance_max:raise ValueError('Invalid rectangular course or clearance box.')

    def concentration(self,time,clearance):
        k=float(clearance);t=np.asarray(time,dtype=float);stop=float(self.stop)
        if k<=0 or np.any(t<0):raise ValueError('Positive clearance and nonnegative time required.')
        return float(self.rate)/k*(-np.expm1(-k*np.minimum(t,stop)))*np.exp(-k*np.maximum(t-stop,0))

    def band_certificate(self,ramp=F(4)):
        ramp=F(ramp)
        if not 0<ramp<=self.stop:raise ValueError('Ramp must lie within administration.')
        taylor=exp_lower(self.clearance_min*ramp,15)
        low=self.rate/self.clearance_max*(1-1/taylor)
        return dict(concentration_upper=self.rate/self.clearance_min,plateau_lower=low,plateau_start=ramp,plateau_stop=self.stop,
            exponential_lower_sum=taylor,administered_amount=self.rate*self.stop,
            full_horizon_exposure_upper=self.rate*self.stop/self.clearance_min,
            washout_exposure_upper=self.rate/self.clearance_min**2,
            scope='One common input over the entire clearance interval; zero initial concentration. Washout remains in the reserve horizon.')


@dataclass(frozen=True)
class HealthyResponse:
    selectivity:F=F(2)
    extra_death:F=F(1,200)
    baseline_numerator:F=F(3,10)
    def __post_init__(self):
        for name in self.__dataclass_fields__:object.__setattr__(self,name,F(getattr(self,name)))
        if self.selectivity<=0 or self.extra_death<0 or self.baseline_numerator<0:raise ValueError('Invalid response parameters.')
    def mortality(self,concentration):return (float(self.baseline_numerator)+concentration)/float(self.selectivity)+float(self.extra_death)
    def ceiling(self,concentration_upper):return (self.baseline_numerator+concentration_upper)/self.selectivity+self.extra_death


class TargetDriftCertificate:
    def __init__(self,source=None,weights=(14,11,10,84,43,107),relative_error=F(1,100)):
        self.source=source or InheritedPopulation();self.weights=tuple(map(F,weights));self.relative_error=F(relative_error)
        if self.source.size!=len(self.weights) or min(self.weights)<=0 or not 0<=self.relative_error<1:raise ValueError('Positive matched weight and valid relative box required.')

    def rows(self,intrinsic):
        action=Action(F(intrinsic)-self.source.rates.intrinsic_erasure);Q=self.source.generator(action);A=self.source.mean_matrix(action);w=self.weights
        ratios=[sum(a*x for a,x in zip(row,w))/w[i] for i,row in enumerate(A)]
        Lw=[sum(p*(w[j]+w[k]) for j,k,p in row) for row in self.source.pairs]
        envelope=[(sum(Q[i][j]*abs(w[j]-w[i]) for j in range(self.source.size) if i!=j)+self.source.birth[i]*abs(Lw[i]-w[i])+self.source.death[i]*w[i])/w[i] for i in range(self.source.size)]
        return ratios,envelope

    def verify(self):
        # The exact paper bands are explicit; no sampling over the rate box.
        rows={str(e):self.rows(e) for e in map(F,['.01','.29','.31'])}
        nominal_growth=F(67,500);nominal_decay=F(99,1000);envelope_cap=F(27,20)
        if max(max(rows[e][0]) for e in ['1/100','31/100'])>nominal_growth:raise ValueError('Growth-band weight fails.')
        if max(max(rows[e][0]) for e in ['29/100','31/100'])>-nominal_decay:raise ValueError('Contraction-band weight fails.')
        if max(rows['31/100'][1])>envelope_cap:raise ValueError('Independent-rate perturbation envelope fails.')
        # Q off-diagonals are affine nondecreasing in e for this source.
        Qlo=self.source.generator();Qhi=self.source.generator(Action(F(3,10)))
        if any(Qhi[i][j]<Qlo[i][j] for i in range(self.source.size) for j in range(self.source.size) if i!=j):raise ValueError('The chosen endpoint does not dominate the error envelope.')
        growth=nominal_growth+self.relative_error*envelope_cap;decay=nominal_decay-self.relative_error*envelope_cap
        if decay<=0:raise ValueError('Rate uncertainty consumes the contraction margin.')
        return dict(rows=rows,envelope_cap=envelope_cap,relative_error=self.relative_error,growth=growth,decay=decay,
            scope='Independent fixed relative factors on positive chemical off-diagonals, division and deaths; reconstruct conservative chemical diagonals. Daughter law, preparation and state space fixed.')

    def certify_course(self,course,initial_counts,ramp=F(4)):
        if len(initial_counts)!=self.source.size or any(not isinstance(n,int) or n<0 for n in initial_counts):raise ValueError('Nonnegative integer founder configuration required.')
        drift=self.verify();band=course.band_certificate(ramp)
        baseline=self.source.rates.intrinsic_erasure
        if baseline!=F(1,100):raise ValueError('This fixed growth-band certificate requires the reference baseline erasure 0.01; recertify a different baseline.')
        if baseline+band['concentration_upper']>F(31,100) or baseline+band['plateau_lower']<F(29,100):raise ValueError('Delivered trajectory does not enter the certified erasure band.')
        exponent=(course.stop-ramp)*drift['decay']-ramp*drift['growth']
        if exponent<=0:raise ValueError('This growth/contraction split gives no useful decay certificate.')
        prefactor=sum(n*w for n,w in zip(initial_counts,self.weights))/min(self.weights)
        lower=exp_lower(exponent,40);risk=min(F(1),prefactor/lower)
        # Preserve the simple published bound when its sufficient conditions hold.
        published=None
        if exponent>F(43,5) and exp_lower(F(43,5),30)>5000:published=min(F(1),prefactor/5000)
        return dict(drift=drift,delivery=band,exponent=exponent,prefactor=prefactor,positive_exp_series=lower,
            exact_rational_upper=risk,published_coarsening=published,
            scope='Finite-time survival at administration stop and every later deadline without immigration. No contraction during washout is assumed.')


def numerical_target(course,clearance=1.,source=None,initial_counts=(0,0,0,0,0,4)):
    model=source or InheritedPopulation();T=float(course.horizon);b=np.array(model.birth,float);d=np.array(model.death,float)
    # Time-dependent backward field sees physical time T-s. Split at washout boundary.
    Q0=np.array(model.generator(),float);R=np.array(model.generator(Action(F(1))),float)-Q0
    def rhs(s,x):
        c=float(course.concentration(T-s,clearance));D=np.array([sum(float(p)*x[j]*x[k] for j,k,p in row) for row in model.pairs])
        return (Q0+c*R)@x+d*(1-x)+b*(D-x)
    x=np.zeros(model.size)
    breaks=[0.,T-float(course.stop),T];breaks=sorted(set(breaks))
    for left,right in zip(breaks,breaks[1:]):
        sol=solve_ivp(rhs,(left,right),x,method='DOP853',rtol=2e-11,atol=2e-13)
        if not sol.success:raise RuntimeError(sol.message)
        x=sol.y[:,-1]
    return dict(extinction_vector=x,survival=1-float(np.prod(x**np.array(initial_counts))),evidence='Nominal deterministic-course backward PGF, numerical only; not the robust drift bound.')


def assemble_mission(target_upper,reserve_upper,target_tolerance=F(1,100),reserve_tolerance=F(1,100)):
    if not 0<=target_upper<=1 or not 0<=reserve_upper<=1:raise ValueError('Marginal upper bounds must be probabilities.')
    return dict(target_upper=target_upper,reserve_loss_upper=reserve_upper,
        target_passes=target_upper<=target_tolerance,reserve_passes=reserve_upper<=reserve_tolerance,
        joint_success_lower=max(F(0),1-target_upper-reserve_upper),
        scope='Union bound on two different failure events. No independence or multiplication of success probabilities.')