import unittest
from fractions import Fraction as F
import numpy as np
import sympy as sp
from example import *
from certificates import local_certificate, donor_certificate, Interval


class ScientificChecks(unittest.TestCase):
    def test_literal_kernel_and_independent_rates(self):
        for dynamic in (False,True):
            m=Reactor(Resident(),Consumers.equal(3),Reservoir.reference(F(1,20)) if dynamic else None)
            rates=tuple(F(i+1,23) for i in range(len(m.reactions)));m=m.with_rates(rates)
            x=[F(i+1,7) for i in range(len(m.names))];A,B,z,H=x[:4];R=x[4] if dynamic else 1;X=x[m.offset:]
            a,b,p,q,al,be,ef,er,u,v,h1,h2,dh=rates[:13];c=m.consumers;load=R*z*sum(k*t for k,t in zip(c.copying,X))
            expected=[a-(p+al)*A+q*B*z+2*ef*B-2*er*A*A,b+p*A-(q*z+be+ef)*B+er*A*A,p*A-q*B*z-u*z+(h1+2*h2)*H-2*v*z*z-load,u*z+v*z*z-(h1+h2+dh)*H]
            if dynamic:expected+=[rates[-2]-rates[-1]*R-load]
            expected += [t*(k*R*z-mu-rho*t) for t,k,mu,rho in zip(X,c.copying,c.loss,c.limitation)]
            self.assertEqual(m.field(x),expected)
            symbols=sp.symbols('x:'+str(len(x)));derivative=sp.Matrix(m.field(symbols)).jacobian(symbols).subs(dict(zip(symbols,x)))
            self.assertEqual(sp.Matrix(m.jacobian(x)),derivative)
            # Any missing consumer stays absent; global permanence requires positive initial data.
            x[-1]=0;self.assertEqual(m.field(x)[-1],0)
        self.assertEqual(len(Reactor(Resident(),Consumers.equal(2),Reservoir.reference(1)).reactions),21)

    def test_composition_and_weighted_variance(self):
        x=sp.symbols('x0:3',positive=True);g=sp.Symbol('g');rho=[F(2),F(3),F(6)];c=Consumers((1,)*3,(F(1,2),)*3,tuple(rho));k,q=c.target_data();self.assertEqual(k,1)
        S=sum(x);H=sum(r*t*t for r,t in zip(rho,x))
        self.assertEqual(sp.expand(H-k*S*S-sum(r*(t-qq*S)**2 for r,t,qq in zip(rho,x,q))),0)
        for i in range(3):
            derivative=sp.expand((x[i]*(g-3*x[i])*S-x[i]*(g*S-3*sum(t*t for t in x)))/S**2)
            self.assertEqual(sp.simplify(derivative-3*S*(x[i]/S)*(sum(t*t for t in x)/S**2-x[i]/S)),0)
        a=Consumers.target(('1/4','3/4'));self.assertEqual(a.limitation,(F(4),F(4,3)))
        with self.assertRaises(ValueError):Consumers.target((.2,.7))
        with self.assertRaises(ValueError):Consumers((1,2),(.5,.5),(2,2)).target_data()

    def test_supply_and_negative_control_identities(self):
        m=Reactor(Resident(),Consumers.equal(2),Reservoir.reference(F(1,20)));x=list(map(F,[3,10,1,8,1,2,3]));f=m.field(x);S=sum(x[5:]);Q=sum(v*v for v in x[5:])
        self.assertEqual(f[4]+sum(f[5:]),F(1,20)*(1-x[4])-S/2-2*Q)
        z,v,Z=sp.symbols('z v Z');self.assertEqual(sp.expand(F(1,4)+z*v+2*Z**2*v**2-(2*(z*v+F(1,4))**2+F(1,8)+2*(Z**2-z*z)*v*v)),0)
        c=Reactor(Resident(),Consumers((1,1),(F(1,2),F(3,4)),(0,0)));x=list(map(F,[3,10,1,8,2,3]));f=c.field(x)
        self.assertEqual((f[-1]*x[-2]-x[-1]*f[-2])/x[-2]**2,-F(1,4)*x[-1]/x[-2])
        for d in (0,1e-10,.05,10):
            b=OperatingAnalysis.ceiling(d);self.assertAlmostEqual(b*b+b/2,d)

    def test_exact_certificates_reject_larger_unproved_region(self):
        c=local_certificate();self.assertTrue(c['accepted']);self.assertEqual(c['consumer_floor'],'203270509477/10000000000000');self.assertGreater(F(c['boundary_margin']),F('5.10e-8'))
        self.assertFalse(local_certificate(radius=F(1,1000))['accepted'])
        self.assertTrue(donor_certificate()['accepted'])
        i=Interval(-2,3);self.assertLessEqual((i*i).lo,0);self.assertGreaterEqual((i*i).hi,9)

    def test_global_floor_representation(self):
        a=PermanenceBounds.aggregate(2,'A');b=PermanenceBounds.aggregate(2,'B',F(1,20))
        self.assertEqual(a.exponent,300000000000);self.assertEqual(a.prefactor,F(1,8*(60000002)))
        self.assertEqual(b.exponent,300000002880);self.assertGreater(b.prefactor,0)
        self.assertIn('aggregate_floor/1000',PermanenceBounds.report(2)['B2']['allowed_rate_radius'])
        with self.assertRaises(ValueError):PermanenceBounds.aggregate(2,'B',0)

    def test_numerical_recovery_balances_and_spectrum(self):
        m=Reactor(Resident.reference(),Consumers.equal(2),Reservoir.reference(.05));times=np.linspace(0,1200,601)
        run=Simulator(m).run(m.initial(),times);other=Simulator(m).run(m.initial(),times,method='BDF')
        self.assertLess(run.balance_error(),2e-8);self.assertLess(np.max(abs(run.state-other.state)),3e-6)
        self.assertLess(abs(run.consumers[-1,0]/run.total[-1]-.5),1e-8)
        eq,error=OperatingAnalysis.equilibrium(m,run.state[-1]);self.assertLess(error,2e-9);self.assertAlmostEqual(sum(eq[5:]),.0406741018954,places=9)
        split=OperatingAnalysis.spectral_split(m,eq);self.assertLess(split['composition_block_error'],1e-10);self.assertLess(split['coupling_error'],1e-10)
        # The full and reduced target manifolds share operating equations at fixed kappa.
        unequal=Reactor(m.resident,Consumers.target((.25,.75)),m.reservoir);u=eq.copy();u[5:]=sum(eq[5:])*np.array([.25,.75]);self.assertLess(max(abs(np.array(unequal.field(u),float))),2e-9)

    def test_finite_inventory_and_model_guards(self):
        m=Reactor(Resident.reference(),Consumers.equal(2),Reservoir.reference(0));initial=m.initial();run=Simulator(m).run(initial,np.linspace(0,100,201))
        self.assertLess(run.ledgers[-1,1],initial[4]+1e-9);self.assertLess(run.ledgers[-1,0],2*(initial[4]+sum(initial[5:]))+1e-9);self.assertLess(run.total[-1],1e-15);self.assertLess(run.balance_error(),2e-8)
        with self.assertRaises(ValueError):m.initial(consumers=(0,1))
        with self.assertRaises(ValueError):m.with_rates([1,2])
        scales=UnitScales(2,3);self.assertEqual(scales.rate_constant(12,3),1);self.assertEqual(scales.rate_constant(12,0),8);self.assertEqual(scales.flux(6),4)


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