"""Independent algebra, model composition and conservation checks."""
import unittest
from collections import Counter
import numpy as np
import sympy as sp
from example import (ConversionArm,ConversionSystem,negative_child,positive_child,
    negative_certificates,FlowerSpectrum,SaturatingKinetics,exact_routh)

class ScienceTests(unittest.TestCase):
    def test_general_conversion_roles_and_reduction(self):
        arms=(ConversionArm('a','b','enzyme','bound','x'),ConversionArm('b','b','enzyme','other','y'))
        network=ConversionSystem(['a','b'],['enzyme'],iter(arms))
        self.assertTrue(np.all(network.conservation()@network.S==0))
        for child in network.selections():
            self.assertEqual(child.reduction(False),sp.Matrix(child.matrix).det())
        with self.assertRaises(ValueError):ConversionSystem(['a'],['a'],[])
        with self.assertRaises(ValueError):network.child(['a'],['xv'])
        with self.assertRaises(ValueError):list(network.selections(limit=2))

    def test_census_and_independent_determinants(self):
        counts=Counter()
        for i,child in enumerate(ConversionSystem.futile(2).selections()):
            determinant=child.reduction(False);counts[determinant]+=1
            if i%23==0:self.assertEqual(determinant,sp.Matrix(child.matrix).det())
        self.assertEqual(counts,{-1:543,0:3377,1:535})

    def test_ordinary_and_scaled_negative_cores(self):
        result=negative_certificates(ConversionSystem.futile(3))
        self.assertEqual(result['A_minus']['routh']['right_half_plane_roots'],2)
        self.assertEqual(result['B']['routh']['right_half_plane_roots'],0)
        self.assertEqual(result['B_scaled']['routh']['right_half_plane_roots'],2)
        self.assertEqual(result['B']['lyapunov_leading_minors'],['1443','1042019','242679562','26482784643','3275722759311'])
        for n in (4,6):self.assertTrue(np.array_equal(negative_child(ConversionSystem.futile(n)).matrix,negative_child(ConversionSystem.futile(3)).matrix))

    def test_positive_embedding_and_extra_unstable_roots(self):
        z=sp.Symbol('z')
        for n in range(2,7):
            child=positive_child(ConversionSystem.futile(n),n);T=np.diag([-1]+[1]*(2*n-1)+[-1])
            self.assertTrue(np.array_equal(T@child.matrix@T,FlowerSpectrum(n).matrix()))
            self.assertEqual(child.reduction(False),1)
            self.assertEqual(sp.expand(sp.Matrix(child.matrix).charpoly(z).as_expr()-((z+1)**(2*n+1)-(z+1)**(2*n-1)-(z+1))),0)
        self.assertEqual(exact_routh(sp.Poly((z+1)**21-(z+1)**19-(z+1),z).all_coeffs())['right_half_plane_roots'],3)

    def test_rate_loss_threshold_and_asymptotics(self):
        for hub_loss,sign in [('1/2',1),('1',0),('2',-1)]:
            f=FlowerSpectrum(3,losses=[hub_loss]+[0]*6)
            self.assertEqual(sp.sign(f.gain_at_zero()-1),sign)
            self.assertAlmostEqual(f.abscissa(),max(np.linalg.eigvals(f.matrix()).real),places=11)
        f=FlowerSpectrum(3,[1,2,3,1,2,1,2]);a=f.abscissa()
        self.assertGreater(a,FlowerSpectrum.unit_alpha(3));self.assertLess(a,3*FlowerSpectrum.unit_alpha(3))
        self.assertAlmostEqual(FlowerSpectrum(3,f.rates,['1/5']*7).abscissa(),a-.2,places=12)
        self.assertTrue(np.allclose(f.matrix()@f.eigenvector(a),a*f.eigenvector(a),atol=1e-12))
        for n in (2,3,10,100,10000):
            lower,upper,_=FlowerSpectrum.lambert_bracket(n);self.assertLess(lower,FlowerSpectrum.unit_alpha(n));self.assertLess(FlowerSpectrum.unit_alpha(n),upper)

    def test_full_exact_kinetic_network(self):
        network=ConversionSystem.futile(3);model=SaturatingKinetics.paper(network)
        self.assertEqual(sp.Matrix(network.S).rank(),9)
        coefficients=(100*model.exact_jacobian).charpoly().all_coeffs()
        self.assertEqual(coefficients,[1,618,139432,13778512,543179240,14097062574,1032202684766,3925174058599,3726175329872,-160160507095,0,0,0])
        rates,field,jac=model.rates_and_jacobian(np.zeros(12))
        self.assertTrue(np.array_equal(field,np.zeros(12)));self.assertTrue(np.allclose(jac,np.array(model.exact_jacobian,float)))
        # A nonunit reference uses the same prescribed reactivity, not an implicit rescaling.
        custom=SaturatingKinetics(network,[sp.Rational(3,2)]*12,model.flux,model.R)
        self.assertTrue(np.allclose(custom.rates_and_jacobian(np.zeros(12))[2],jac))

    def test_trajectory_conservation_and_nonlinear_jacobian(self):
        network=ConversionSystem.futile(3);model=SaturatingKinetics.paper(network);L=network.conservation()
        initial=np.array(network.S[:,0],float)*1e-4;times=np.linspace(0,12,61)
        a=model.integrate(initial,times);b=model.integrate(initial,times,'BDF')
        self.assertLess(np.max(abs(a-b)),2e-9);self.assertLess(np.max(abs(a@L.T-(1+initial)@L.T)),1e-11)
        eta=np.linspace(-.05,.04,12);jac=model.rates_and_jacobian(eta)[2];step=1e-6
        numeric=np.column_stack([(model.rates_and_jacobian(eta+step*np.eye(12)[i])[1]-model.rates_and_jacobian(eta-step*np.eye(12)[i])[1])/(2*step) for i in range(12)])
        self.assertLess(np.max(abs(jac-numeric)),1e-9)

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