import unittest,json
from pathlib import Path
from fractions import Fraction as Q
import numpy as np
import phos_sharp as ps
import phos_capacity as pc
from model import Geometry,Kinetics,Reactor,ordered_design,stability
from recovery import LocalRecovery

class ScientificChecks(unittest.TestCase):
    @classmethod
    def setUpClass(cls):
        root=Path(__file__).parent;cls.ex=json.loads((root/'worked_example.json').read_text());cls.cert=json.loads((root/'operational_certificate.json').read_text());cls.rates=[[Q(v) for v in row] for row in cls.cert['rates']];cls.states=[list(map(Q,p['state'])) for p in cls.ex['record']['profiles']]
    def test_exact_five_equilibria_and_full_stability(self):
        counts=[]
        for z in self.states:
            self.assertTrue(all(v==0 for v in ps.vector_field(self.rates,z)));self.assertEqual(ps.totals(z),(18,2,20));counts.append(stability(self.rates,z)['unstable'])
        self.assertEqual(counts,[0,1,0,1,0]);readout=[z[3]+z[11] for z in self.states];self.assertEqual(readout,sorted(readout))
    def test_kinetic_retuning_keeps_geometry_and_totals(self):
        g=Geometry(('2','7/2','5','13/2','8'),Q(9));rec=g.construct()
        for eps in (Q(1),Q(1,10)):
            rates=Kinetics((1,2,3),(2,3,4),(3,4,5),eps).rates(g)
            for state in rec['states']:self.assertTrue(all(v==0 for v in ps.vector_field(rates,state['z'])));self.assertEqual(ps.totals(state['z']),(18,2,20))
        with self.assertRaises(ValueError):Kinetics((1,),(1,),(1,)).rates(g)
    def test_full_reactor_matches_independent_exact_jacobian(self):
        reactor=Reactor(self.rates,(18,2,20));z=np.array(self.states[2],float);q=z[reactor.indices];J=reactor.jacobian(0,q)
        np.testing.assert_allclose(J,np.array(ps.reduced_jacobian(self.rates,self.states[2]),float),atol=1e-14)
        h=1e-6;numeric=np.column_stack([(reactor.field(0,q+h*np.eye(9)[i])-reactor.field(0,q-h*np.eye(9)[i]))/(2*h) for i in range(9)])
        np.testing.assert_allclose(J,numeric,atol=1e-7,rtol=1e-7)
        q[2]*=1.001;run=reactor.integrate(reactor.species(q),10,samples=21);self.assertLess(run['conservation_drift'],1e-12)
    def test_fresh_recovery_and_out_of_scope_radius(self):
        cert=LocalRecovery(self.cert,self.states);good=cert.evaluate();self.assertTrue(good['all_certified']);self.assertTrue(good['readout_intervals_disjoint']);self.assertLessEqual(max(s['tenfold_seconds_upper'] for s in good['sinks']),32918000)
        self.assertFalse(cert.evaluate('1/1000')['all_certified']);self.assertFalse(cert.evaluate(measurement_error='1')['readout_intervals_disjoint'])
    def test_ordered_design_and_budget_status(self):
        design=ordered_design(3,4);self.assertEqual(design['status'],'certified_design');self.assertEqual(design['r'],4);self.assertEqual(design['delta'],Q(1,2));self.assertEqual(design['unstable'],[0,1,0,1,0])
        self.assertEqual(ordered_design(3,0)['status'],'unknown_budget')
    def test_loaded_matrix_and_diagonal_stability(self):
        n=4;state=pc.coalesced(n,Q(100))['states'][0]['z'];M,h,K=pc.loaded(n,state,Q(200),Q(2));self.assertTrue(all(v>0 for v in pc.ldl_pivots(M)));self.assertEqual(ps.det(K),0)
        lim=pc.limit_data(n);v=(3*n-2)*Q(2,9)**(n-1);self.assertEqual(1-lim['G'],v*(50*n-45)/(9*(14-8*v)))
        self.assertTrue(all(sum(a*b for a,b in zip(row,lim['z']))==-(1-lim['G']) for row in lim['Bstar']))
    def test_base_case_and_routh_uncertainty(self):
        rates=[(Q(1),)*6];z=[Q(v) for v in (2,2,1,1,1,1)];self.assertEqual(stability(rates,z)['unstable'],0)
        co=pc.coalesced(3,Q(4));self.assertEqual(stability(co['rates'],co['states'][0]['z'])['status'],'unresolved_regular_routh')
        with self.assertRaises(ValueError):Geometry(('1','2','3'),Q(9))
        with self.assertRaises(ValueError):Reactor(rates,(0,2,6))

if __name__=='__main__':unittest.main()