"""Scientific regression checks of the actual source and evidence contracts."""
import unittest
import numpy as np
import sympy as sp
from reactor import EquilibriumDesign,MassActionReactor,ConstantBinding
from example import design,certify_references
import phos


class ScientificChecks(unittest.TestCase):
    def test_reaction_source_and_conservation(self):
        m=design().exact();r=MassActionReactor(m)
        self.assertEqual(m.field(),sp.zeros(12,1));self.assertEqual(m.N.rank(),9)
        self.assertEqual(tuple(m.totals()),tuple(map(sp.Rational,['119/25','1379/50','1106/25'])))
        np.testing.assert_allclose(r.conservation@r.N,0,atol=0)
        y=np.linspace(-.001,.001,9)
        np.testing.assert_allclose(r.P@r.field(0,y),r.full_field(0,r.species(y)),atol=1e-13)
        np.testing.assert_allclose(r.jacobian(0,np.zeros(9)),np.array(m.jacobian(),dtype=float),atol=1e-12)

    def test_arbitrary_site_composition_and_input_checks(self):
        d=design().append_site().append_site('1/1000');m=d.exact()
        self.assertEqual(m.n,5);self.assertEqual(m.field(),sp.zeros(18,1))
        with self.assertRaises(ValueError):design().append_site('0')
        with self.assertRaises(ValueError):design().with_reverse('-1').exact()
        with self.assertRaises(ValueError):MassActionReactor(design().exact(),relaxation=2)

    def test_static_equality_and_dynamic_disagreement(self):
        a=MassActionReactor(design().with_reverse('13/10').exact());b=MassActionReactor(design().with_reverse('3/2').exact())
        self.assertGreater(max(np.linalg.eigvals(a.jacobian(0,np.zeros(9))).real),0)
        self.assertLess(max(np.linalg.eigvals(b.jacobian(0,np.zeros(9))).real),0)
        np.testing.assert_allclose(a.static_jacobian(),b.static_jacobian(),atol=1e-12)
        for u in [np.zeros(3),np.array([.003,-.002,.001])]:
            za,fa=a.static_complexes(u);zb,fb=b.static_complexes(u)
            np.testing.assert_allclose(za,zb,atol=1e-10);np.testing.assert_allclose(fa,fb,atol=1e-10)

    def test_fresh_hopf_sign_and_reduction_certificates(self):
        records,raw=certify_references()
        self.assertLess(raw['A']['l1'].b,0);self.assertGreater(raw['H']['l1'].a,0)
        self.assertLess(raw['A']['crossing'].re.b,0);self.assertGreater(raw['A']['eps_derivative'].re.a,0)
        self.assertTrue(all(v.a>0 for v in raw['A']['routh']))
        self.assertGreater(raw['A_kinase_clamped']['l1'].a,0)
        self.assertLess(raw['A_added_site']['l1'].b,0)
        self.assertLess(sp.Rational(records['A_C2_eliminated']['l1'][1]),0)
        self.assertLess(abs(phos.mid(raw['A']['r'])-1.433228306811195),1e-13)

    def test_buffering_and_relaxation(self):
        m=design().exact(True,True);r=MassActionReactor(m)
        self.assertLess(max(np.linalg.eigvals(r.jacobian(0,np.zeros(9))).real),0)
        x=r.species(np.linspace(-.001,.001,9));self.assertEqual(r.full_field(0,x)[4],0);self.assertEqual(r.full_field(0,x)[5],0)
        a=MassActionReactor(design().exact());b=MassActionReactor(design().exact(),relaxation=.5)
        np.testing.assert_allclose(a.static_jacobian(),b.static_jacobian(),atol=1e-12)
        np.testing.assert_allclose(b.jacobian(0,np.zeros(9))[3:],2*a.jacobian(0,np.zeros(9))[3:],atol=1e-12)

    def test_reversible_currents_and_numerical_jacobian(self):
        r=MassActionReactor(design().exact(),drive=12)
        y=np.linspace(-.001,.001,9);f=r.full_field(0,r.species(y))
        np.testing.assert_allclose(r.conservation@f,0,atol=1e-13)
        fd=np.column_stack([(r.field(0,y+np.eye(9)[i]*1e-6)-r.field(0,y-np.eye(9)[i]*1e-6))/2e-6 for i in range(9)])
        np.testing.assert_allclose(r.jacobian(0,y),fd,atol=2e-7,rtol=1e-7)
        self.assertGreater(max(abs(r.field(0,np.zeros(9)))),0)
        with self.assertRaises(ValueError):r.with_rates(r.kon,r.koff,r.kcat)

    def test_finite_time_flux_ledger_and_pulse(self):
        r=MassActionReactor(design().exact());y=np.zeros(9);y[2]=.01
        sol,X=r.integrate(y,2,31,ledger=True);counts=sol.y[9:,-1]
        np.testing.assert_allclose(sol.y[:3,-1]-y[:3],counts[:3]-counts[3:],atol=3e-11)
        self.assertGreater(X.min(),0)
        D=[]
        for ratio in ['13/10','3/2']:
            p=MassActionReactor(phos.witness('A',ratio),ConstantBinding(1.5));_,z=p.integrate(np.zeros(9),.05,11);D.append(z[9,-1])
        self.assertGreater(D[1]-D[0],.0001764) # nominal result only, not the uncertain pulse theorem


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