"""Exact paper witnesses and independent numerical checks of reusable dynamics."""
import unittest
from fractions import Fraction as F
import numpy as np
from example import (Network,EquilibriumModel,RationalLift,BASE_REACTANTS,BASE_PRODUCTS,
    EQUILIBRIUM_FLUX,COMPACT_STATE,EXACT_STATE,EXACT_PADDED_ORDERS,compact_network,
    characteristic,quartic_certificate,child_certificates,base_stability_certificate,
    exact_eigenpair,determinant)


class KineticOrderTests(unittest.TestCase):
    def setUp(self): self.base=Network(BASE_REACTANTS,BASE_PRODUCTS)

    def test_exact_compact_quartic(self):
        model=EquilibriumModel(compact_network(self.base,200),COMPACT_STATE,EQUILIBRIUM_FLUX)
        self.assertEqual(characteristic(model.jacobian()),[F(1),F(87840586,423),F(325113896,5499),F(22009472,1833),F(2105344,611)])
        certificate=quartic_certificate(model.jacobian())
        self.assertEqual(F(certificate['delta2']),F(2196782093181776,178929))
        self.assertEqual(F(certificate['delta3']),-F(384290295376441702400,327976857))
        self.assertEqual(certificate['classification'],'two right-half-plane roots')
        # Independent polynomial evaluation through determinants at distinct z.
        for z in (F(-2),F(0),F(1,3),F(3),F(11)):
            value=F(0)
            for coefficient in characteristic(model.jacobian()): value=value*z+coefficient
            self.assertEqual(value,determinant(z*np.eye(4,dtype=int)-model.jacobian()))

    def test_principal_witness_and_generic_lift(self):
        orders=self.base.reactants.copy()
        for (i,j),v in zip(((2,1),(2,2),(3,0),(3,3)),EXACT_PADDED_ORDERS): orders[i,j]=v
        model=EquilibriumModel(self.base.pad_to(orders),EXACT_STATE,EQUILIBRIUM_FLUX)
        self.assertEqual(exact_eigenpair(model.jacobian())['residual'],'exactly zero')
        lift=RationalLift.build(self.base,EQUILIBRIUM_FLUX,model.reactivity())
        self.assertTrue(np.array_equal(lift.reactivity(),model.reactivity()))
        self.assertTrue(self.base.same_skeleton(lift.network))
        self.assertTrue(np.array_equal(lift.jacobian(),model.jacobian()))
        self.assertTrue(np.all(np.isfinite(model.log_rate_constants())))
        # A second network and target exercises lifting outside the paper skeleton.
        cycle=Network([[1,0],[0,1]],[[0,1],[1,0]])
        other=RationalLift.build(cycle,[2,2],[[F(1,3),0],[0,F(5,7)]])
        self.assertTrue(np.array_equal(other.reactivity(),[[F(1,3),0],[0,F(5,7)]]))

    def test_all_children_and_all_scalings_certificate(self):
        certificate=child_certificates(self.base)
        self.assertEqual(certificate['child_count'],24);self.assertEqual(certificate['covered_count'],24)
        base=base_stability_certificate(self.base,EQUILIBRIUM_FLUX)
        self.assertEqual(sorted(map(F,base['cubic_hurwitz_gap'].values())),sorted(map(F,[73728,36864,67584,58368,12288,16896,6144])))
        original={(s,r):a.tolist() for s,r,a in self.base.children()}
        padded={(s,r):a.tolist() for s,r,a in compact_network(self.base,200).children()}
        self.assertEqual(original,padded)

    def test_rate_jacobian_complex_step(self):
        model=EquilibriumModel(compact_network(self.base,200),COMPACT_STATE,EQUILIBRIUM_FLUX)
        np.testing.assert_allclose(model.rates(model.x),model.v,rtol=0,atol=0)
        np.testing.assert_allclose(model.rhs(0,model.x),np.zeros(4),rtol=0,atol=0)
        numeric=np.empty((4,4));h=1e-25
        for i in range(4):
            state=model.x.astype(complex);state[i]+=h*1j
            numeric[:,i]=model.rhs(0,state).imag/h
        np.testing.assert_allclose(numeric,np.array(model.jacobian(),dtype=float),rtol=1e-12,atol=1e-9)

    def test_nonlinear_solver_against_exact_reversible_pair(self):
        network=Network([[1,0],[0,1]],[[0,1],[1,0]])
        model=EquilibriumModel(network,[2,3],[2,2]);times=np.linspace(0,4,41)
        result=model.simulate([1,4],times,rtol=1e-10,atol=1e-12)
        exact=2-np.exp(-F(5,3).__float__()*times)
        np.testing.assert_allclose(result[0],exact,rtol=2e-10,atol=2e-10)
        np.testing.assert_allclose(result.sum(axis=0),5,rtol=0,atol=2e-9)

    def test_invalid_inputs(self):
        with self.assertRaises(ValueError): Network([[F(1,2)]],[[1]])
        with self.assertRaises(ValueError): EquilibriumModel(self.base,COMPACT_STATE,[1,1,1,1,1])
        with self.assertRaises(ValueError): compact_network(self.base,3)
        orders=self.base.reactants.copy();orders[0,1]=1
        with self.assertRaises(ValueError): self.base.pad_to(orders)
        with self.assertRaises(ValueError): RationalLift.build(self.base,EQUILIBRIUM_FLUX,[[0]*4]*5)
        model=EquilibriumModel(self.base,COMPACT_STATE,EQUILIBRIUM_FLUX)
        with self.assertRaises(ValueError): model.rates([0,1,1,1])


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