import unittest from dataclasses import replace import numpy as np import sympy as sp from example import (Parameters,CoinfectionModel,ReactionSource,Reaction,InvasionBlock, NormalizedLift,ParameterBox,ResidentEntropy,Interval) class ScientificChecks(unittest.TestCase): @classmethod def setUpClass(cls):cls.model=CoinfectionModel() def test_source_balance_and_boundary_residents(self): m=self.model;p=m.p;s,a,b,c=m.x expected=sp.Matrix([ p.recruitment-s*(p.mu0+p.alpha1*a+p.alpha2*b+(p.alpha3+p.beta1+p.beta2)*c), a*(p.alpha1*s-p.gamma1*b-p.eta1*c-p.mu1)+p.beta1*s*c, b*(p.alpha2*s-p.gamma2*a-p.eta2*c-p.mu2)+p.beta2*s*c, (p.gamma1+p.gamma2)*a*b+c*(p.eta1*a+p.eta2*b+p.alpha3*s-p.mu3)]) self.assertEqual((m.F-expected).applyfunc(sp.expand),sp.zeros(4,1)) self.assertEqual(sp.expand(sum(m.F)-(p.recruitment-p.mu0*s-p.mu1*a-p.mu2*b-p.mu3*c)),0) self.assertEqual(m.source.siphons(),[(),(1,3),(2,3),(1,2,3)]) self.assertFalse(m.source.is_siphon([3])) for E in m.residents():self.assertEqual(m.F.subs(dict(zip(m.x,E))),sp.zeros(4,1)) with self.assertRaises(ValueError):Parameters(beta2=0) def test_overlap_at_threshold_and_uncovered_mode(self): m=self.model for s in (0,1,4): report=m.source.overlap_certificate(m.rates,[s,0,0,0],(1,3),(2,3)) self.assertEqual(report['cross_multiplied_residual'],'0') atzero=m.source.overlap_certificate(m.rates,[1,0,0,0],(1,3),(2,3)) self.assertEqual(atzero['chi_intersection'],'z') # An additional missing compartment has an uncovered positive mode. reactions=[Reaction(r.rate_name,r.reactant+(0,),r.product+(0,)) for r in m.source.reactions] reactions.append(Reaction('extra',(0,0,0,0,1),(0,0,0,0,2))) source=ReactionSource(('s','a','b','c','d'),reactions);rates=dict(m.rates,extra=sp.Integer(1)) report=source.overlap_certificate(rates,[4,0,0,0,0],(1,3),(2,3),normal_indices=(1,2,3,4)) empty=next(r for r in report['signature_factors'] if not r['membership']) self.assertEqual(empty['factor'],'z - 1') z=sp.Symbol('z') self.assertEqual(sp.expand(sp.sympify(report['chi_full_normal'])-sp.sympify(report['chi_union'])*(z-1)),0) with self.assertRaises(ValueError):m.source.overlap_certificate(m.rates,m.residents()[1],(1,3),(2,3)) def test_exact_invasion_threshold_covectors_and_dfe_obstruction(self): low=CoinfectionModel(Parameters(beta2='1/10')) high=CoinfectionModel();zero=CoinfectionModel(Parameters(beta2='51/140')) changed=CoinfectionModel(Parameters(beta2='1/10',eta1='1/5')) self.assertEqual(low.permanence()['status'],'permanence_fails_resident_sink') self.assertEqual(high.permanence()['status'],'uniform_permanence_by_paper_theorem') self.assertEqual(zero.permanence()['status'],'threshold_not_decided') self.assertEqual(changed.permanence()['status'],'uniform_permanence_by_paper_theorem') self.assertEqual(low.J.subs(dict(zip(low.x,low.residents()[0]))),changed.J.subs(dict(zip(changed.x,changed.residents()[0])))) for m in (low,high,changed): for M in m.invasion_blocks(): cert=InvasionBlock(M).certificate();r=sp.Rational(cert['left_weight']) row=sp.Matrix([[1,r]])*M sign=1 if cert['status']=='strict_invasion' else -1 self.assertTrue(all(sign*v>0 for v in row)) missing=CoinfectionModel(Parameters(alpha1='1/10')) self.assertEqual(missing.permanence()['status'],'outside_resident_existence_hypotheses') self.assertIsNone(high.permanence()['numerical_floor']) def test_interval_certificate_entire_independent_box(self): certificate=ParameterBox(Parameters()).certificate() self.assertTrue(certificate['certified']) self.assertEqual(certificate['relative_margins'],['5601/101000','8701/25250']) self.assertEqual([row[0] for row in certificate['resident_s_intervals']],['99/202','99/101']) self.assertEqual([row[0] for row in certificate['resident_u_intervals']],['68207/19998','29003/9999']) # Additional substitutions audit the enclosure implementation, not the proof of the box. rng=np.random.default_rng(2401) for _ in range(12): center=Parameters();rates={n:getattr(center,n)*sp.Rational(int(rng.choice([99,101])),100) for n in center.__dataclass_fields__} M1,M2=CoinfectionModel(Parameters(**rates)).invasion_blocks() actual=[*(sp.Matrix([[1,2]])*M1),*(sp.Matrix([[1,1]])*M2)] actual[1]/=2 for v,bound in zip(actual,certificate['weighted_relative_column_intervals']): self.assertGreaterEqual(v,sp.Rational(bound[0]));self.assertLessEqual(v,sp.Rational(bound[1])) self.assertFalse(ParameterBox(Parameters(beta2='1/10')).certificate()['certified']) self.assertEqual((Interval(-2,3)*Interval(-5,7)).pair(),['-15','21']) with self.assertRaises(ValueError):Interval(1,2)/Interval(-1,1) def test_normalized_quotient_field_fibre_growth_and_compensation(self): m=self.model;lift=NormalizedLift(m) state=sp.Matrix([sp.Rational(3,2),sp.Rational(1,4),sp.Rational(2,3),sp.Rational(1,5)]) extended=lift.lift(state);F=m.F.subs(dict(zip(m.x,state))) s,a,b,c=m.x;U=a+lift.ru*c;V=b+lift.rv*c;J=a+b+lift.rj*c directions=sp.Matrix([a/U,b/V,a/J,b/J]) derivative=(directions.jacobian(m.x)*m.F).subs(dict(zip(m.x,state))) self.assertEqual((lift.field(extended)[4:,:]-derivative).applyfunc(sp.simplify),sp.zeros(4,1)) product=U*V*J**lift.k relative=(sp.Matrix([product]).jacobian(m.x)*m.F)[0]/product self.assertEqual(sp.simplify(relative.subs(dict(zip(m.x,state)))-lift.growth(extended)[3]),0) E0,E1,E2=m.residents() for r in (0,sp.Rational(1,3),1): self.assertGreater(lift.growth(list(E1)+[1,r,1,0])[3],0) self.assertGreater(lift.growth(list(E2)+[r,1,0,1])[3],0) difficult=NormalizedLift(CoinfectionModel(Parameters(eta1=10,eta2=10,mu3=10))) E0=difficult.model.residents()[0] hu,hv,hj,G=difficult.growth(list(E0)+[0,0,0,0]) self.assertLess(hu+hv,0);self.assertGreater(G,0);self.assertGreater(difficult.k,1) with self.assertRaises(ValueError):lift.lift([4,0,0,0]) def test_entropy_strictification_and_conditional_recovery(self): model=CoinfectionModel(Parameters(mu0='6/5',mu1='13/10')) entropy=ResidentEntropy(model);s,u=sp.symbols('s u',positive=True) sb,ub=entropy.sbar,entropy.ubar;alpha=entropy.alpha;g=entropy.g;e=entropy.epsilon v=s-sb;w=u-ub;ds=-g*v-alpha*s*w;du=alpha*u*v H=s-sb-sb*sp.log(s/sb)+u-ub-ub*sp.log(u/ub)+e*v*w exact=-(g/s-e*alpha*u)*v*v-e*g*v*w-e*alpha*s*w*w self.assertEqual(sp.simplify(sp.diff(H,s)*ds+sp.diff(H,u)*du-exact),0) self.assertLessEqual(e*alpha*entropy.R,g/(4*entropy.R)) self.assertLessEqual(e*g*entropy.R,alpha*entropy.ell/2) lift=NormalizedLift(self.model);report=lift.conditional_recovery('1/1000') self.assertTrue(report['conditional_only']);self.assertGreater(sp.Rational(report['epsilon']),0) # Recover the private-compartment forcing identity exactly, conditional on U>=mU. p=self.model.p;s,a,b,c=self.model.x;R,ell=self.model.bounds();mU=sp.Rational(report['U_floor']) Ca=p.beta1*ell*mU/lift.ru;Ka=p.mu1+(p.gamma1+p.eta1)*R+p.beta1*ell/lift.ru rhs=p.alpha1*s*a+p.gamma1*a*(R-b)+p.eta1*a*(R-c)+p.beta1*(s-ell)*c+p.beta1*ell/lift.ru*(a+lift.ru*c-mU) self.assertEqual(sp.expand(self.model.F[1]-(Ca-Ka*a)-rhs),0) def test_positive_integration_balance_and_independent_solver(self): m=self.model;initial=[.5,3.5,1e-5,1e-5];times=np.linspace(0,100,251) states=m.simulate(initial,times);other=m.simulate(initial,times,method='BDF') expected=4+(sum(initial)-4)*np.exp(-times) np.testing.assert_allclose(states.sum(axis=1),expected,atol=1e-8,rtol=1e-9) np.testing.assert_allclose(states,other,atol=2e-7,rtol=2e-6) self.assertTrue(np.all(states>0)) sink=CoinfectionModel(Parameters(beta2='1/10')).simulate(initial,times) self.assertLess(sink[-1,2]+sink[-1,3],1e-3*(initial[2]+initial[3])) if __name__=='__main__':unittest.main()