"""Independent source, reduction, boundary and numerical regression checks.""" import itertools import unittest from dataclasses import replace import numpy as np import sympy as sp from example import (Q, UnitPath, separated_type_ii, coincident_type_ii, TypeVTopology, StationaryCertificate, vector) class ScientificChecks(unittest.TestCase): def test_source_and_physical_jacobian(self): source=separated_type_ii(((2,3),(1,),(2,1))) n=len(source.names); model=source.realize([1]*n,[Q(i%3,10) for i in range(n)]) symbols=sp.symbols(f'x0:{n}',positive=True) independent=model.field(symbols).jacobian(symbols)*sp.diag(*symbols) self.assertEqual(sp.simplify(independent-model.scaled_jacobian(symbols)),sp.zeros(n)) self.assertEqual(source.P[source.names.index('G0_1'),source.names.index('F0')],2) self.assertEqual(source.P[source.names.index('T1'),source.names.index('G0_1')],3) for i in range(3): self.assertEqual(source.P[source.names.index(f'F{i}'),source.names.index(f'T{i}')],1) with self.assertRaises(ValueError): separated_type_ii(((1,), (1,))) with self.assertRaises(ValueError): separated_type_ii(((1,), (), (1,))) def test_passive_path_recurrence_and_internal_balance(self): for losses in [(0,), (1,), (0,2,0), (1,0,3)]: m=len(losses); path=UnitPath(tuple(range(1,m+2)),tuple(range(2,m+3)),losses) x=path.reconstruct(Q(7,3),Q(5,2)); c,b,L,R=path.compress() currents=[path.forward[k]*x[k]-path.reverse[k]*x[k+1] for k in range(m+1)] self.assertEqual(currents[0],(c+L)*x[0]-b*x[-1]) self.assertEqual(currents[-1],c*x[0]-(b+R)*x[-1]) for k in range(m): self.assertEqual(currents[k]-currents[k+1],Q(losses[k])*x[k+1]) self.assertTrue(all(v>0 for v in x)); self.assertGreater(c,0);self.assertGreater(b,0) self.assertEqual(UnitPath((1,1),(1,1),(1,)).compress(),(Q(1,3),)*4) self.assertEqual(UnitPath((1,1),(1,1),(0,)).compress(),(Q(1,2),Q(1,2),0,0)) with self.assertRaises(ValueError): UnitPath((1,1),(1,0),(0,)) def test_expanded_type_ii_stationary_schur_reduction(self): source=separated_type_ii(); n=len(source.names); model=source.realize([1]*n,[0]*n) path=UnitPath((2,1,3),(1,4,2),(0,Q(1,3))) expanded,compressed=model.expand_unit_edge('T1','F1',path) retained=vector([Q(i+2,3) for i in range(n)]) route=path.reconstruct(retained[source.names.index('T1')],retained[source.names.index('F1')]) full=sp.Matrix([*retained,*route[1:-1]]) self.assertEqual(expanded.field(full)[:n],list(compressed.field(retained))) self.assertEqual(expanded.field(full)[n:], [0,0]) with self.assertRaises(ValueError): model.expand_unit_edge('F0','T1',path) def test_all_eight_type_v_patterns_and_ratio_tangent(self): for linked in itertools.product((0,2),repeat=3): top=TypeVTopology(linked); source=top.source(); n=len(source.names) # Singleton loss supports include paths, bases and fork sources. for support in range(n): x=vector([Q(i+2,3) for i in range(n)]); d=[Q(1,5) if i==support else Q(0) for i in range(n)] model=source.realize(x,d); reduced=top.reduce(model); u=vector([x[source.names.index(f'U{i}')] for i in range(3)]) self.assertEqual(reduced.residual(u),sp.zeros(3,1)); self.assertEqual(reduced.reconstruct(u),x) self.assertNotEqual(reduced.tangent(u).det(),0) # Off-root equation transformation agrees with literal pair balances. trial=vector([Q(2,3),Q(5,4),Q(7,5)]); full=reduced.reconstruct(trial); f=model.field(full) E=vector([sum(f[source.names.index(s)] for s in top.names_for_pair(i)) for i in range(3)]) self.assertEqual(reduced.residual(trial),2*(sp.ones(3)-2*sp.eye(3)).inv()*E) for ratios in itertools.product((Q(1,2),Q(1),Q(3,2)),repeat=3): witness=reduced.ratio_witness(ratios,[1,1,1]) self.assertEqual(witness['status'],'identical' if ratios==(1,1,1) else 'cannot_both_be_stationary') def test_stationary_factorization_and_multiaffine_bound(self): source=separated_type_ii(); n=len(source.names) for d in ([0]*n,[Q(1,3) if i==0 else 0 for i in range(n)],[Q(i%2,5) for i in range(n)]): x=vector([Q(i+2,2) for i in range(n)]); model=source.realize(x,d,Q(2,3)); cert=StationaryCertificate(model,x) self.assertEqual(list(cert.t),[Q(2,3)]*n) audit=cert.principal_audit(); self.assertEqual(audit['status'],'all_principal_minors_nonnegative') self.assertGreaterEqual(Q(audit['det_minus_H']),Q(audit['lower_bound_from_nonnegative_minors'])) self.assertEqual(cert.principal_audit(1)['status'],'unresolved_budget') with self.assertRaises(ValueError): StationaryCertificate(model,[2]*n) # Admissibility must handle a strictly negative current, not clip it. other=coincident_type_ii((2,3,4)); m=other.realize([1]*3,[1,0,0]); c=StationaryCertificate(m,[1]*3) self.assertLess(c.J[0],0); self.assertEqual(c.shift[0],-c.J[0]); self.assertGreater(c.p[0],0) def test_coincident_branch_inverse_and_zero_loss_stability(self): for weights in ((1,1,1),(2,3,4)): source=coincident_type_ii(weights); a,b,c=weights self.assertEqual(source.N.det(),a+b+c+a*b*c) self.assertTrue(all(source.N.inv()[i,j]>0 for i in range(3) for j in range(3) if i!=j)) for source in (source,separated_type_ii(),TypeVTopology((3,0,2)).source()): n=len(source.names); model=source.realize([Q(i+2,3) for i in range(n)],[Q(i%2,10) for i in range(n)]) zero=replace(model,loss=(0,)*n); x=zero.detailed_balance(); f,H=zero.numeric(x) self.assertLess(np.max(np.abs(f/x)),1e-10) N=np.array(source.N,float); q=np.array(zero.reverse,float)*np.exp(np.array(source.P,float).T@np.log(x)) symmetric=-np.diag(1/np.sqrt(x))@N@np.diag(q)@N.T@np.diag(1/np.sqrt(x)) self.assertLess(max(np.linalg.eigvalsh(symmetric)),0) np.testing.assert_allclose(np.sort(np.linalg.eigvals(H/x[None,:]).real),np.linalg.eigvalsh(symmetric),rtol=1e-10,atol=1e-10) def test_fixed_rate_continuation_and_independent_roots(self): from scipy.optimize import root source=separated_type_ii(); n=len(source.names); model=source.realize([1]*n,[Q(i%2,5) for i in range(n)]) rows=model.continuation([0,.01,.1,.5,1]) self.assertEqual(len(rows),5); self.assertTrue(all(r['status']=='numerical_root' for r in rows)) np.testing.assert_allclose(rows[-1]['x'],np.ones(n),atol=1e-8) check=root(lambda x:model.numeric(x)[0],np.full(n,1.01),jac=lambda x:model.numeric(x)[1]/x[None,:]) self.assertTrue(check.success); np.testing.assert_allclose(check.x,rows[-1]['x'],atol=1e-9) self.assertLess(rows[0]['spectral_abscissa'],0) if __name__=='__main__': unittest.main()