import unittest
from fractions import Fraction as Q
import numpy as np
import phos_sharp as ps
from construction import PrescribedStates,KineticFreedom,EquilibriumChart,stability,substrate_window
from reactor import Reactor

class ScientificChecks(unittest.TestCase):
    def setUp(self):self.rec=PrescribedStates(('2','3','4','5','6'),Q(5)).construct()
    def test_literal_mass_action_and_prescribed_ratios(self):
        for s in self.rec['states']:
            self.assertEqual(s['z'][4]/s['z'][5],s['u']);self.assertEqual(ps.totals(s['z']),(10,2,12));self.assertFalse(any(ps.vector_field(self.rec['rates'],s['z'])))
        self.assertEqual([stability(self.rec['rates'],s['z'])['unstable'] for s in self.rec['states']],[0,1,0,1,0])
    def test_positivity_is_not_only_r_greater_than_roots(self):
        roots=('1059/1000','1098/1000','1003/1000','1029/1000','1038/1000')
        with self.assertRaises(ValueError):PrescribedStates(roots,Q(64,100)).construct()
        self.assertTrue(PrescribedStates(roots,Q(66,100)).construct()['positive'])
        with self.assertRaises(ValueError):PrescribedStates(('2','2','3'),5).construct()
    def test_equilibrium_count_and_extraneous_roots(self):
        a=EquilibriumChart(self.rec['rates'],self.rec['totals']).audit();self.assertEqual(a['positive_count'],5)
        self.assertTrue(any(row['status']!='positive_regular_equilibrium' for row in a['candidates']))
    def test_exceptional_chart_not_lost(self):
        a=EquilibriumChart([(2,1,1,2,1,1)],(1,1,2)).audit();self.assertEqual(a['positive_count'],1);self.assertEqual(a['exceptional']['status'],'positive_exceptional_equilibrium')
    def test_substrate_window_and_failed_widening(self):
        probes=list(map(Q,['1/20','69/128','4327/3200','7917/3200','6373/1600','399/80']))
        self.assertEqual(substrate_window(self.rec,probes,Q(1,500))['status'],'certified_at_least_5');self.assertEqual(substrate_window(self.rec,probes,Q(1,100))['status'],'not_certified')
    def test_independent_kinetic_freedom(self):
        rates=KineticFreedom((1,2,3),(2,3,4),(3,4,5)).realize(self.rec)
        for s in self.rec['states']:self.assertFalse(any(ps.vector_field(rates,s['z'])))
        self.assertNotEqual(ps.reduced_jacobian(rates,self.rec['states'][0]['z']),ps.reduced_jacobian(self.rec['rates'],self.rec['states'][0]['z']))
    def test_reactor_jacobian_and_conserved_totals(self):
        r=Reactor(self.rec['rates'],self.rec['totals']);z=np.array(self.rec['states'][0]['z'],float);q=z[r.indices]
        np.testing.assert_allclose(r.jacobian(0,q),np.array(ps.reduced_jacobian(self.rec['rates'],self.rec['states'][0]['z']),float),atol=1e-12)
        q[2]*=.99;run=r.integrate(r.species(q),100,samples=11);self.assertLess(run['conservation_drift'],1e-12)
        z[0]+=.1
        with self.assertRaises(ValueError):r.integrate(z,100)

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