"""Exact cusp, equivalence and kinetic checks independent of saved outputs."""
import copy
import unittest
import numpy as np
import sympy as sp
from example import (Catalogue,PlanarNetwork,AlgebraicTemplate,RationalCuspReactor,
    rational_identities,REACTIONS,x,y,t)

class ScienceTests(unittest.TestCase):
    @classmethod
    def setUpClass(cls):cls.catalogue=Catalogue()
    def test_literal_universe_and_full_coverage(self):
        result=self.catalogue.census()
        self.assertEqual((result['literal_sets'],result['eligible_sets'],result['species_orbits'],result['mechanism_classes']),(142506,60036,30051,9999))
        self.assertEqual(result['layers'],{'determinant':9063,'fold':744,'cubic':140,'cusp':52})
        with self.assertRaises(ValueError):self.catalogue.census(limit=10)

    def test_labelled_transport_and_ray_collision(self):
        source=PlanarNetwork((0,5,12,18,22));target=PlanarNetwork((22,18,12,5,2))
        match=target.transport_from(source);rates=list(map(sp.Rational,['1/11','3/11','3/11','3/11','1/11']))
        transported=[rates[j]*f for j,f in zip(match['source_indices'],match['rate_factors'])]
        self.assertEqual((source.symbolic_field(rates)-target.symbolic_field(transported)).applyfunc(sp.expand),sp.zeros(2,1))
        self.assertEqual(source.key(),target.key());self.assertTrue(self.catalogue.classify(target)['admits_cusp'])
        swapped=PlanarNetwork(source.swapped());mapping=swapped.transport_from(source)
        self.assertTrue(mapping['species_exchange'])
        a=PlanarNetwork((0,5,12,18,20));b=PlanarNetwork((0,5,12,20,24))
        self.assertEqual(sorted(a.vectors),sorted(b.vectors));self.assertNotEqual(a.key(),b.key())
        self.assertTrue(self.catalogue.classify(a)['admits_cusp']);self.assertEqual(self.catalogue.classify(b)['layer'],'cubic')

    def test_rational_and_algebraic_cusp_jets(self):
        for ordinal in (0,2,7):
            result=AlgebraicTemplate(self.catalogue.templates[ordinal]).verify()
            self.assertLess(sp.Rational(result['sign_enclosures']['tau'][1]),0)
            self.assertLess(sp.Rational(result['sign_enclosures']['c'][1]),0)
            if ordinal==0:self.assertNotEqual(result['raw_minor'],result['corrected_minor'])
        row=AlgebraicTemplate(self.catalogue.templates[2]).verify()
        self.assertEqual(row['c'],'-75/17303');self.assertEqual(row['corrected_minor'],'-297/845')
        self.assertEqual(row['raw_minor'],row['corrected_minor'])  # This particular pair has zero determinant correction.
        bad=copy.deepcopy(self.catalogue.templates[0]);bad['U']='0';bad['V']='0'
        with self.assertRaises(ArithmeticError):AlgebraicTemplate(bad).verify()

    def test_exact_cubic_and_three_equilibria(self):
        identities=rational_identities();self.assertEqual(sp.expand(sp.sympify(identities['scalar'])-(sp.Symbol('mu')+sp.Symbol('nu')*sp.Symbol('z')-sp.Symbol('z')**3)),0)
        model=RationalCuspReactor.bistable('1/10');rows=model.equilibria()
        self.assertEqual([r['stability'] for r in rows],['stable','saddle','stable'])
        self.assertEqual([r['determinant_midpoint'] for r in rows],['3/6050','-3/12100','3/6050'])
        self.assertEqual([r['trace_midpoint'] for r in rows],['-61/55','-1299/1100','-69/55'])
        # At the cusp the positive root has multiplicity three; at a fold, two.
        self.assertEqual(RationalCuspReactor.unfolding(0,0).equilibria()[0]['multiplicity'],3)
        self.assertEqual(sorted(r['multiplicity'] for r in RationalCuspReactor.unfolding('-1/500','3/100').equilibria()),[1,2])
        self.assertEqual(len(RationalCuspReactor.unfolding('1/1000','1/100').equilibria()),1)

    def test_general_rates_and_scalar_scope(self):
        model=RationalCuspReactor(['1/8','2/7','4/9','1/3','1/12']);F=model.network.symbolic_field(model.rates)
        self.assertEqual(sp.expand(F[0].subs(y,model.rates[3]/model.rates[2]*x*x)-model.equilibrium_polynomial().as_expr()),0)
        self.assertEqual(sp.expand(F[1].subs(y,model.rates[3]/model.rates[2]*x*x)),0)
        # y=x^2 is not an invariant curve: y'=0 there while 2x*x' generally is not.
        witness=RationalCuspReactor.bistable('1/10');field=witness.field((.8,.64))
        self.assertAlmostEqual(field[1],0);self.assertNotAlmostEqual(2*.8*field[0],0)
        with self.assertRaises(ValueError):RationalCuspReactor.bistable(1)

    def test_template_instantiation(self):
        template=AlgebraicTemplate(self.catalogue.templates[0]);model=template.numerical_reactor()
        self.assertLess(np.max(abs(model.field((1,1)))),1e-14)
        target=PlanarNetwork((2,5,12,18,20));transported=template.numerical_reactor(target)
        self.assertLess(np.max(abs(model.field((.7,1.3))-transported.field((.7,1.3)))),1e-14)

    def test_full_trajectory_bistability(self):
        model=RationalCuspReactor.bistable('1/10');times=np.linspace(0,30000,101)
        for initial,expected in [((.85,.7),(.9,.81)),((1.15,1.4),(1.1,1.21))]:
            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[-1]-expected)),2e-6)
            self.assertGreater(a.min(),0)

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