"""Scientific checks independent of the reaction-list implementation."""
import unittest
import numpy as np
from example import (A, B, I, Y, Reactor, Rates, Intervention, Integrator,
                     Schedule, INITIAL_STATE, operate)


class ReactorTests(unittest.TestCase):
    def test_literal_field_and_balances(self):
        model = Reactor()
        for c in [np.array(INITIAL_STATE), np.array([.3,.2,.1,.02,.03,.04]), np.zeros(6)]:
            u,w,x,c1,c2,z = c
            j0=u*w/500000000-x/5000000000
            j1=20*x*u-20*c1; j2=20*c1*w-20*c2
            j3=20*c2-2*z; j4=19*(z-x*x); j5=.04*x-.04*u*w/8000000000
            expected=np.array([1-u-j0-j1+j5,1-w-j0-j2+j5,-x+j0-j1+2*j4-j5,
                               -c1+j1-j2,-c2+j2-j3,-z+j3-j4])
            np.testing.assert_allclose(model.derivative(c),expected,atol=1e-14)
            self.assertAlmostEqual(A@model.derivative(c),1-A@c)
            self.assertAlmostEqual(B@model.derivative(c),1-B@c)
            self.assertAlmostEqual(I@model.derivative(c),j0+j3-j5-I@c)
        np.testing.assert_array_equal(A@model.stoichiometry,np.zeros(6))
        np.testing.assert_array_equal(B@model.stoichiometry,np.zeros(6))

    def test_pulse_accounting_and_invalid_inputs(self):
        c=np.array(INITIAL_STATE)
        pulse=Intervention(.4,(.98,1,.98,1,.99,1),(-.005,.005))
        after,withdrawn,loss=pulse.apply(c)
        np.testing.assert_allclose(after+withdrawn+loss,c+np.r_[pulse.refill,np.zeros(4)])
        self.assertAlmostEqual(I@after+I@withdrawn+I@loss,I@c)
        for build in [lambda: Rates(release=-1),lambda: Rates(cleavage=float('nan')),
                      lambda: Intervention(1.1),lambda: Intervention(survival=(1,)),
                      lambda: Schedule(cycles=1.5),lambda: pulse.apply([-1]*6)]:
            with self.assertRaises(ValueError): build()
        self.assertFalse(Rates(release=22).in_paper_box)
        self.assertFalse(Intervention(.1).in_paper_box)

    def test_growth_guard_is_essential(self):
        # Paper's counterexample to globally positive catalytic drift.
        c=np.array([0,0,1,0,0,0.])
        self.assertLess(Y@Reactor().derivative(c),0)
        # Pure catalytic-phase states on the low-stock boundary.
        for phase in range(2,6):
            c=np.zeros(6); c[phase]=.05/Y[phase]
            c[0]=.9-A@c; c[1]=.9-B@c
            self.assertGreaterEqual(Y@Reactor().derivative(c)-2/3*(Y@c),-1e-14)

    def test_weak_preparation_and_continued_cycles(self):
        # Slowest pure phase from the paper, with minimum admitted stock.
        c=np.zeros(6); c[4]=(1/5000)/Y[4]
        c[0]=.9-A@c; c[1]=.9-B@c
        result=operate(Integrator(Reactor()),initial=c,schedule=Schedule(cycles=3))
        self.assertTrue(result.summary['theorem_inputs_satisfied'])
        self.assertTrue(result.summary['all_cycle_endpoints_in_operating_region'])
        self.assertLess(result.summary['maximum_material_error'],1e-8)
        self.assertLess(abs(result.summary['final_inventory_balance_residual']),1e-8)
        for row in result.cycles:
            self.assertGreaterEqual(row['template_collected'],1/28)
            self.assertGreaterEqual(row['free_X_collected'],1/540)
            self.assertLessEqual(row['gross_service'],9/50)
            self.assertLessEqual(row['food_U'],951/200)
            self.assertGreaterEqual(row['cumulative_synthesis'],row['paper_synthesis_lower_bound'])
        for n in range(1,4):
            end=[r for r in result.trajectories if r['cycle']==n-1][-1]
            for name in ('U','W','X','C1','C2','Z'):
                self.assertEqual(result.pulses[n][f'before_{name}'],end[name])

    def test_exploration_has_no_theorem_bound(self):
        result=operate(Integrator(Reactor(Rates(release=22))),schedule=Schedule(cycles=1))
        self.assertFalse(result.summary['theorem_inputs_satisfied'])
        self.assertIsNone(result.cycles[0]['paper_synthesis_lower_bound'])


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