"""Scientific checks for the coupled reactor and the scope of its conclusions.""" import unittest from dataclasses import replace from fractions import Fraction as Q import numpy as np from example import (Parameters, CoupledReactor, StationaryResponse, ResponsePotential, NumericalSelector, FoldCertificate, Interval) class ScientificChecks(unittest.TestCase): def test_reactions_jacobian_and_activity_balance(self): for g in (Q(0),Q(1,10),Q(1)): model=CoupledReactor(Parameters(g=g,eta=Q(1,1000000))) for x in ((2,3,5,7),(Q(13,2),Q(17,3),Q(9,5),11)): literal=np.zeros(4,dtype=object) for pair in model.pairs:literal+=pair.change*pair.current(x) np.testing.assert_allclose(np.array(literal,dtype=float),model.field(0,x),atol=2e-13) f=model.field(0,x);activity=model.activity(x) self.assertAlmostEqual(activity[0]+2*activity[1]+activity[2],f[2],places=11) point=np.array(x,dtype=float);step=1e-5 numeric=np.column_stack([(model.field(0,point+np.eye(4)[i]*step)-model.field(0,point-np.eye(4)[i]*step))/(2*step) for i in range(4)]) np.testing.assert_allclose(numeric,model.jacobian(point),rtol=1e-6,atol=1e-8) self.assertEqual(model.pairs[0].forward+model.pairs[1].forward,1) self.assertEqual(model.pairs[0].reverse+model.pairs[1].reverse,1) def test_clamped_modules_preserve_all_boundary_terms(self): for g in (Q(1,10),Q(1)): for d in (Q(1,10000),Q(2)): model=CoupledReactor(Parameters(g=g,d=d)) response=StationaryResponse(model) for z in (.1,1,9): A,B=response.clamped_ab(z) np.testing.assert_allclose(model.field(0,(A,B,z,4))[:2],0,atol=1e-12) for A,B in ((1,20),(30,2),(10,15)): z,H=response.clamped_zh(A,B) np.testing.assert_allclose(model.field(0,(A,B,z,H))[2:],0,atol=2e-12) def test_exact_root_counts_positive_reconstruction_and_spectra(self): for e in (Q(1,200000),Q(1,100000),Q(1,50000)): model=CoupledReactor(Parameters(e=e));response=StationaryResponse(model) roots=response.exact_roots();self.assertEqual(len(roots),3) for i,root in enumerate(roots): self.assertEqual(root['multiplicity'],1) lo,hi=map(Q,root['z_interval']);self.assertLessEqual(hi-lo,Q(1,10**14)) x=np.array(root['state']);np.testing.assert_allclose(model.field(0,x),0,atol=1e-11) eig=np.linalg.eigvals(model.jacobian(x)) self.assertLess(max(abs(eig.imag)),1e-10) self.assertEqual(sum(eig.real>0),int(i==1)) prod=model.activity(x) self.assertLess(prod[1],-6);self.assertGreater(prod[2],4);self.assertGreater(prod[3],.0005) for g,count in [(Q(1,10),1),(Q(1999999,2000000),3)]: roots=StationaryResponse(CoupledReactor(Parameters(g=g))).exact_roots() self.assertEqual(len(roots),count) isolated=CoupledReactor(Parameters(g=0)) x=StationaryResponse(isolated).isolated_state() np.testing.assert_allclose(isolated.field(0,x),0,atol=1e-12) self.assertEqual(len(StationaryResponse(CoupledReactor(Parameters(d=1))).exact_roots()),1) def test_potential_derivative_dissipation_and_strict_selection(self): model=CoupledReactor();pot=ResponsePotential(model) roots=StationaryResponse(model).exact_roots();states=[np.array(r['state']) for r in roots] selector=NumericalSelector(pot,states[1]) self.assertEqual(selector.inspect(states[0]),'low_z_sink') self.assertEqual(selector.inspect(states[2]),'high_z_sink') self.assertEqual(selector.inspect(states[1]),'unresolved') self.assertEqual(selector.inspect((40,20,2,20)),'outside_absorbing_region') rng=np.random.default_rng(19) for _ in range(30): B=rng.uniform(2.1,30);A=rng.uniform(.1,33.9-B) x=np.array([A,B,rng.uniform(.1,11.9),rng.uniform(.1,100)]) self.assertTrue(pot.inside(x)) self.assertLessEqual(pot.derivative(x)+pot.dissipation(x),1e-10) x=np.array([12.,19.,1.2,10.]);f=model.field(0,x);dt=1e-6 finite=(pot.value(x+dt*f)-pot.value(x-dt*f))/(2*dt) self.assertAlmostEqual(finite,pot.derivative(x),delta=2e-7) def test_trajectories_select_opposite_compositions(self): model=CoupledReactor();states=[np.array(r['state']) for r in StationaryResponse(model).exact_roots()] vals,vecs=np.linalg.eig(model.jacobian(states[1]));v=vecs[:,np.argmax(vals.real)].real v=v/np.linalg.norm(v)*np.sign(v[1]);pot=ResponsePotential(model) selector=NumericalSelector(pot,states[1]) for sign,destination,label in [(1,states[0],'low_z_sink'),(-1,states[2],'high_z_sink')]: sol=model.integrate(states[1]+sign*.02*v) sampled=sol.sol(np.linspace(0,1200,121)).T self.assertLess(max(abs(sampled[-1]-destination)),.002) self.assertEqual(selector.inspect(sampled[-1]),label) energy=[pot.value(x) for x in sampled] self.assertLessEqual(max(np.diff(energy)),1e-8) def test_fixed_fold_exact_contraction(self): result=FoldCertificate.check() self.assertLess(Q(result['contraction_norm']),Q(11,10**8)) self.assertGreater(Q(result['Qg'][0]),0) self.assertLess(Q(result['Qzz'][1]),0) self.assertGreater(Q(result['hurwitz_gap'][0]),10000) self.assertEqual((Interval(-2,3)**2).as_json(),['0','9']) with self.assertRaises(ZeroDivisionError):Interval(-1,1)**-1 def test_parameter_changes_do_not_inherit_flagship_theorem(self): for change in ({'g':Q(99,100)},{'eta':Q(1,10**6)},{'d':Q(1)},{'e':Q(1,100)}): with self.assertRaises(ValueError):ResponsePotential(CoupledReactor(replace(Parameters(),**change))) with self.assertRaises(ValueError):Parameters(g=2) with self.assertRaises(ValueError):CoupledReactor().integrate((0,1,1,1)) model=CoupledReactor(Parameters(eta=Q(1,10**6))) for root in StationaryResponse(model).exact_roots(): np.testing.assert_allclose(model.field(0,root['state']),0,atol=1e-11) if __name__=='__main__':unittest.main()