import math import unittest from fractions import Fraction as F from collections import Counter import numpy as np import sympy as sp from polymer import Channel,PolymerCatalogue,PolymerReactor,CatalyticAssignment,RetainedReward,FOOD,SELECTED from certificate import ReactorCertificate,Bernstein,exp_negative_upper,source_mass_alpha2_n4,readout_floor,EPS from example import generator_panel class ScientificChecks(unittest.TestCase): def test_literal_catalogue_roles_and_conservation(self): cat=PolymerCatalogue(4);self.assertEqual((len(cat.words),len(cat.splits)),(30,68)) rx=PolymerReactor(cat,CatalyticAssignment(weak_food=1e-5)) for ch in rx.channels: if ch.kind=='internal':self.assertEqual(sum(len(z)*d for z,d in ch.change.items()),0) repeated=Channel('repeated',('0011',)*3,('0011','00110011'),4.) self.assertEqual(repeated.change['0011'],-2) self.assertAlmostEqual(repeated.rate({'0011':3,'00110011':0},10),.24) self.assertEqual(repeated.rate({'0011':2,'00110011':0},10),0) N=rx.initial(40);flux=rx.flux(N/40,40)*40 counts=dict(zip(cat.words,map(int,N))) np.testing.assert_allclose(flux,[ch.rate(counts,40) for ch in rx.channels],rtol=1e-14,atol=1e-16) self.assertAlmostEqual(float(rx.length@(rx.stoich@flux))/40,0,places=12) def test_exact_internal_detailed_balance_including_repeated_roles(self): # Explicit product-Poisson weight and literal falling factors, independently computed. v=F(5);theta0,theta1=F(2,3),F(4,5) pairs=[(('00','11','0011'),('0011','0011')),(('00','00','00'),('0000','00'))] for inputs,outputs in pairs: words=set(inputs)|set(outputs);nu=Counter(inputs);nv=Counter(outputs) N={z:4 for z in words};M={z:N[z]-nu[z]+nv[z] for z in words} star={z:theta0**z.count('0')*theta1**z.count('1') for z in words} def weight(state):return math.prod((v*star[z])**state[z]/math.factorial(state[z]) for z in words) def rate(state,roles):return v**(1-sum(roles.values()))*math.prod(math.prod(state[z]-i for i in range(j)) for z,j in roles.items()) self.assertEqual(weight(N)*rate(N,nu),weight(M)*rate(M,nv)) self.assertEqual(math.prod(star[z] for z in inputs),math.prod(star[z] for z in outputs)) def test_exact_scalar_certificates_and_provenance(self): a,b=sp.symbols('a b') self.assertEqual(sp.expand(4*(1-a)*(1-b)+8*(a*a+b*b)-F(16,5)-5*(a+b-F(2,5))**2-3*(a-b)**2),0) self.assertLess(F(2718281828,10**9),sum((F(1,math.factorial(j)) for j in range(40)),F(0))) self.assertGreater(1/exp_negative_upper(44),F(33,4)*10**18) ell=F(1,1358)-F(1,50000) self.assertGreater(EPS*ell*ell/10,F(1,10**17)+F(1558,9)*F(1,2*10**18)) self.assertLess(F(9,10)**10,F(1,2)) specifications=[(2*10**22,F(1,10**5),48,F(1,10**10),F(101,1000)), (2*10**22,F(1,10**4),37,F(1,10**6),F(1002,10000)), (5*10**22,F(1,10**5),40,F(1,10**18),F(135,1000)), (5*10**22,F(2,10**4),37,F(1,10**15),F(101,1000)), (10**23,F(1,10**5),40,F(1,10**37),F(142,1000)), (10**24,F(5,10**4),8,F(1,10**14),F(1,10))] for v,eta,z,error,quota in specifications: r=ReactorCertificate(v,v//10**8,eta).evaluate(F(z)) self.assertLess(r['error_upper'],error);self.assertGreater(r['quota_99'],quota) for vv in [v+1,v+3*10**18-1,3*v+7]: s=ReactorCertificate(vv,v//10**8,eta).evaluate(F(z)) self.assertGreaterEqual(s['quota_99'],r['quota_99']);self.assertLessEqual(s['error_upper'],r['error_upper']) self.assertEqual(ReactorCertificate().floor,6667) self.assertEqual(ReactorCertificate(10**20).evaluate()['status'],'UNAVAILABLE') self.assertEqual(ReactorCertificate(10**22,4,F(3,100)).evaluate()['status'],'QUOTA_NOT_CERTIFIED') def test_all_windows_inverse_and_readout(self): r=ReactorCertificate(10**24,10**15).evaluate(F(14),True) self.assertLess(r['error_upper'],F(1,10**20));self.assertGreater(r['quota_198'],F(51,100)) service=ReactorCertificate(10**24).constants() self.assertGreater(service.quota(F(796,10),F(14),F(1,40),True),F(1,10)) self.assertGreater(r['quota_198']/2400,F(1,4700)) z=service.reward.inverse(1e-6,2,99) self.assertAlmostEqual(2*math.exp(-float(service.reward.exponent(F(z),99))),1e-6,places=17) self.assertEqual(Bernstein(F(0),F(0)).upper(1,1),0) lengths={'000':3,'0011':4} self.assertEqual(readout_floor(F(135,1000),{'000':F(285,100),'0011':F(380,100)},lengths,F(2,1000))['lower'],F(12625,100000)) self.assertFalse(readout_floor(F(135,1000),{'0011':4},lengths,F(2,1000))['informative']) def test_reward_keeps_every_loss_and_food_increment(self): volume=2*10**22;counts={z:volume for z in FOOD};counts['0011']=6667 reward=RetainedReward(volume) basal=Channel('basal',('00','11'),('0011',),2e-9) selected=Channel('selected',('00','11','0011'),('0011','0011'),4.,selected_birth=True) loss=Channel('loss',('0011','0011','0011'),('00110011','0011'),4.) for ch in [basal,selected,loss]: kept,full,omitted=reward.increments(counts,ch) self.assertLessEqual(kept,full);self.assertGreaterEqual(omitted,0) self.assertLess(reward.increments(counts,basal)[0],0) # Food penalty retained even on omitted birth. self.assertEqual(reward.increments(counts,selected)[2],0) self.assertAlmostEqual(reward.increments(counts,loss)[0],math.log1p(-2/6667)) panel=generator_panel(PolymerCatalogue(4),volume) self.assertTrue(all(r['max_jump']<=r['jump_upper'] for r in panel)) self.assertGreater(panel[0]['omitted_log_variance'],panel[0]['variance']) def test_continuous_material_account_and_stochastic_no_restart(self): rx=PolymerReactor(PolymerCatalogue(4));times=np.array([0.,1.,100.,199.]);y=rx.deterministic(times);size=len(rx.catalogue.words) np.testing.assert_allclose(y[:,:size]@rx.length,10,rtol=1e-9) np.testing.assert_allclose(y[:,size+1],6*times,rtol=1e-10) np.testing.assert_allclose(y[:,:size]@rx.nonfood+y[:,size]-y[:,size+3],0,atol=1e-8) result=rx.simulate(2,199,1234,10000) self.assertTrue(result['completed']);self.assertEqual(result['ledger_residual'],0) self.assertEqual(len(result['observations']),3);self.assertFalse(result['mission']) limited=rx.simulate(2,199,1234,2) self.assertFalse(limited['completed']);self.assertIsNone(limited['mission']) def test_source_mass_and_invalid_assignment(self): r=source_mass_alpha2_n4();lo,hi=r['witness_mass'] self.assertGreater(lo,F(187365,10**8));self.assertLess(hi,F(187366,10**8)) self.assertLess((1-F(187365,10**8)*(1-F(1,10**10)))**1600,F(1,20)) cat=PolymerCatalogue(4) for assignment in [CatalyticAssignment(selected=3),CatalyticAssignment(weak_food=.1), CatalyticAssignment(background=(('0011',SELECTED,4),))]: with self.assertRaises(ValueError):PolymerReactor(cat,assignment) with self.assertRaises(ValueError):ReactorCertificate(2*10**22,4,F(1,10)) if __name__=='__main__':unittest.main()