"""Fresh local robustness and recovery calculations for the rational benchmark."""
from fractions import Fraction as Q
from math import ceil
import phos_sharp as ps
import phos_capacity as pc
from paper_checks import structured_bounds

class LocalRecovery:
    def __init__(self,certificate,states):
        self.data=certificate;self.states=states;self.rates=[[Q(v) for v in row] for row in certificate['rates']];self.totals=list(map(Q,certificate['totals']))
        if len(self.rates)!=3:raise ValueError('This benchmark certificate has exactly three sites.')
    def evaluate(self,eta=None,measurement_error='1/1000',concentration='1/10',time_unit='10'):
        eta=Q(self.data['common_independent_relative_rate_total_radius'] if eta is None else eta);measurement=Q(measurement_error);C0=Q(concentration);T0=Q(time_unit)
        if not 0<=eta<1 or measurement<0 or min(C0,T0)<=0:raise ValueError('Admissible relative radius, nonnegative error and positive units required.')
        results=[]
        for stored in self.data['sinks']:
            j=stored['index'];state=list(map(Q,self.states[j]));q=state[1:4]+state[6:];scales=list(map(Q,stored['coordinate_scales']))
            if q!=scales or ps.totals(state)!=tuple(self.totals) or any(ps.vector_field(self.rates,state)):raise ArithmeticError('Nominal source does not match the certificate.')
            J=ps.reduced_jacobian(self.rates,state);J0=[[J[i][k]*scales[k]/scales[i] for k in range(9)] for i in range(9)]
            P=[[Q(v) for v in row] for row in stored['P']];JtP=pc.matmul(pc.transpose(J0),P);Qm=[[-JtP[i][k]-JtP[k][i] for k in range(9)] for i in range(9)]
            if P!=pc.transpose(P) or not all(v>0 for v in pc.ldl_pivots(P)+pc.ldl_pivots(Qm)):raise ArithmeticError('Nominal Lyapunov definiteness failed.')
            rho=Q(stored['rho']);b=structured_bounds(q,self.totals,self.rates,scales,P,pc.inverse(J0),rho,eta)
            M=pc.norm_inf(P);pmin=1/pc.norm_inf(pc.inverse(P));qlow=1/pc.norm_inf(pc.inverse(Qm));c=qlow-b['lyapunov_loss_bound'];rin=Q(stored['initial_scaled_euclidean_radius'])
            valid=b['contraction']<=Q(1,2) and b['equilibrium_shift_bound']<=rho/4 and b['positive_species_lower']>0 and c>0 and rin*rin*M<=pmin*(rho/2)**2
            if not valid:
                results.append(dict(index=j,status='not_certified',bounds={k:str(v) for k,v in b.items()},scope='Failure of this sufficient certificate does not establish loss of the physical state.'));continue
            readout=state[3]+state[11];radius=readout*b['equilibrium_shift_bound'];interval=[C0*(readout-radius)-measurement,C0*(readout+radius)+measurement]
            results.append(dict(index=j,status='certified',parameter_radius=str(eta),box_radius=str(rho),initial_euclidean_radius=str(rin),M=str(M),pmin=str(pmin),c=str(c),bounds={k:str(v) for k,v in b.items()},tenfold_seconds_upper=ceil(6*T0*M/c),decay_exponent_per_second=str(c/(2*T0*M)),readout_interval_micromolar=list(map(str,interval)),scope='Joint independent perturbations of all eighteen rates and three totals. Initial error is measured from the perturbed equilibrium, with nominal coordinate scales. Return concerns an envelope inside a local invariant ellipsoid.'))
        certified=all(r['status']=='certified' for r in results)
        disjoint=certified and all(Q(a['readout_interval_micromolar'][1])<Q(b['readout_interval_micromolar'][0]) for a,b in zip(results,results[1:]))
        return dict(sinks=results,all_certified=certified,readout_intervals_disjoint=disjoint)