| import unittest |
| from itertools import combinations |
| import numpy as np |
| import sympy as s |
| from ahti.core import (HodgeModel, closest_vector, ResourceLimit, integer_affine_solve, |
| constrained_quantize, simplicial_complex, smith_verified) |
| from ahti.moments import encode,decode,signed_permutations,match_face |
| from ahti.geometry import tetrahedral_operators,jacobian,no_flip_certificate,monotonicity_certificate |
|
|
| Q=s.Rational |
|
|
| class ExactAlgebra(unittest.TestCase): |
| def test_three_way_decomposition(self): |
| B=s.Matrix([2,0,0,0]); C=s.Matrix([[0,2,1,1]]) |
| V=s.Matrix([[1,2,0,0],[0,1,1,0],[0,0,1,1],[0,0,0,1]]) |
| model=HodgeModel.build(B,C,V.T*V+s.eye(4)) |
| a=s.Matrix([Q(1,2),Q(2,3),Q(-3,5),Q(7,11)]) |
| e,h,l=model.split(a) |
| self.assertEqual(e+h+l,a) |
| self.assertEqual(model.norm2(a),sum(model.norm2(x) for x in (e,h,l))) |
| self.assertEqual(C*h,s.zeros(1,1)) |
| self.assertEqual(B.T*model.W*h,s.zeros(1,1)) |
| def test_reject_noncomplex(self): |
| with self.assertRaises(ValueError): HodgeModel.build(s.eye(2),s.eye(2)) |
| def test_reject_bad_metric_and_float(self): |
| with self.assertRaises(ValueError): HodgeModel.build(s.zeros(2,0),s.zeros(0,2),s.diag(1,-1)) |
| with self.assertRaises(ValueError): HodgeModel.build(s.Matrix([[1.0]]),s.zeros(0,1)) |
| def test_saturated_integer_kernel(self): |
| model=HodgeModel.build(s.zeros(3,0),s.Matrix([[2,1,1]])) |
| z=s.Matrix([0,1,-1]) |
| result=model.repair(z,'integral') |
| self.assertEqual(result['cost'],0) |
| self.assertEqual(model.integral_lattice()['free_rank'],2) |
| def test_torsion_is_not_lost(self): |
| model=HodgeModel.build(s.diag(-2,-2,0),s.zeros(0,3)) |
| lattice=model.integral_lattice() |
| self.assertEqual(lattice['torsion'],[2,2]) |
| self.assertEqual(lattice['free_rank'],1) |
| result=model.repair(s.Matrix([Q(1,3),Q(1,4),Q(7,5)]),'integral') |
| self.assertEqual(result['cost'],Q(4,25)) |
| self.assertEqual(result['repaired'],s.Matrix([Q(1,3),Q(1,4),1])) |
| def test_exact_projection_ignores_torsion_in_real_target(self): |
| model=HodgeModel.build(s.Matrix([[-2]]),s.zeros(0,1)) |
| result=model.repair(s.Matrix([Q(1,3)]),'exact') |
| self.assertEqual(result['cost'],0) |
| self.assertEqual(model.integral_lattice()['torsion'],[2]) |
| def test_local_and_global_costs_add(self): |
| model=HodgeModel.build(s.zeros(3,0),s.Matrix([[1,1,0]])) |
| a=s.Matrix([Q(2,5),Q(1,5),Q(7,5)]) |
| result=model.repair(a,'integral') |
| lattice=model.integral_lattice()['harmonic_basis'] |
| e,h,l=model.split(a) |
| cvp=closest_vector(h,lattice,model.W) |
| self.assertEqual(result['cost'],model.norm2(l)+cvp['cost']) |
| def test_cvp_ties(self): |
| result=closest_vector(s.Matrix([Q(1,2)]),s.eye(1)) |
| self.assertEqual(result['ties'],[(0,),(1,)]) |
| self.assertEqual(result['cost'],Q(1,4)) |
| def test_oblique_cvp_against_larger_enumeration(self): |
| L=s.Matrix([[1,Q(3,5)],[0,Q(4,5)]]) |
| for i in range(-5,6): |
| h=s.Matrix([Q(i,7),Q(2*i+1,9)]) |
| result=closest_vector(h,L) |
| brute=min(((h-L*s.Matrix([k,j])).dot(h-L*s.Matrix([k,j]))) |
| for k in range(-6,7) for j in range(-6,7)) |
| self.assertEqual(result['cost'],brute) |
| def test_resource_limit_is_not_a_certificate(self): |
| with self.assertRaises(ResourceLimit): |
| closest_vector(s.Matrix([Q(1,2),Q(1,2)]),s.eye(2),max_candidates=1) |
| def test_integer_affine_solutions(self): |
| A=s.Matrix([[2,4,0],[0,0,3]]); c=s.Matrix([6,9]) |
| k,V=integer_affine_solve(A,c) |
| self.assertEqual(A*k,c);self.assertEqual(A*V,s.zeros(A.rows,V.cols)) |
| with self.assertRaises(ValueError):integer_affine_solve(s.Matrix([[2]]),s.Matrix([1])) |
| def test_smith_transformations(self): |
| for A in (s.Matrix([[2,1,1]]),s.diag(2,6,0),s.zeros(0,3),s.zeros(3,0)): |
| d,u,v=smith_verified(A) |
| self.assertEqual(d,u*A*v) |
| self.assertEqual(abs(u.det()),1);self.assertEqual(abs(v.det()),1) |
| def test_boundary_quantization(self): |
| B=s.Matrix([[-1,1,0],[0,-1,1],[1,0,-1]]) |
| model=HodgeModel.build(B,s.zeros(0,3)) |
| a=s.Matrix([Q(1,5),Q(2,5),Q(2,5)]) |
| M=s.Matrix([[1,0,0],[0,1,0]]);N=s.zeros(2,1);d=s.Matrix([0,0]) |
| out=constrained_quantize(model,a,M,N,d) |
| self.assertEqual(M*out['u']+N*out['k'],d) |
| Z=model.integral_lattice()['integer_representatives'] |
| self.assertEqual(B*out['u']+Z*out['k'],out['repaired']) |
| brute=[] |
| for k in range(-5,6): |
| x=B[:,2];y=a-Z*s.Matrix([k]) |
| t=(x.T*y)[0]/(x.T*x)[0] |
| brute.append((y-x*t).dot(y-x*t)) |
| self.assertEqual(out['cost'],min(brute)) |
| def test_boundary_impossible_integer_constraint(self): |
| model=HodgeModel.build(s.zeros(1,1),s.zeros(0,1)) |
| with self.assertRaises(ValueError): |
| constrained_quantize(model,s.Matrix([0]),s.zeros(1,1),s.Matrix([[2]]),s.Matrix([1])) |
| def test_flat_simplicial_connection(self): |
| R=s.Matrix([[0,-1,0],[1,0,0],[0,0,1]]) |
| gauges=[s.eye(3),R,R*R] |
| edges=[(0,1),(0,2),(1,2)] |
| transports={(i,j):gauges[i]*gauges[j].T for i,j in edges} |
| B,C=simplicial_complex(3,edges,[(0,1,2)],transports) |
| self.assertEqual(C*B,s.zeros(3,9)) |
| model=HodgeModel.build(B,C) |
| self.assertEqual(model.PH.rank(),0) |
| def test_reject_curved_triangle(self): |
| R=s.diag(-1,-1,1) |
| with self.assertRaises(ValueError): |
| simplicial_complex(3,[(0,1),(0,2),(1,2)],[(0,1,2)],{(0,1):R}) |
|
|
| class Geometry(unittest.TestCase): |
| @classmethod |
| def setUpClass(cls): |
| cls.X=s.Matrix([[0,0,0],[1,0,0],[0,1,0],[0,0,1],[0,0,-1]]) |
| cls.ops=tetrahedral_operators(cls.X,[(0,1,2,3),(0,1,2,4)]) |
| cls.model=HodgeModel.build(cls.ops['B'],cls.ops['C'],cls.ops['W']) |
| cls.u=s.Matrix([v for row in cls.X.tolist() for v in row]) |
| cls.a=cls.ops['B']*cls.u |
| def test_reconstruction(self): |
| for E in self.ops['reconstruction']: |
| self.assertEqual(jacobian(E,self.a),s.eye(3)) |
| def test_geometric_weight_identity(self): |
| rhs=sum(mu*sum(x*x for x in jacobian(E,self.a)) |
| for mu,E in zip(self.ops['volumes'],self.ops['reconstruction'])) |
| self.assertEqual(self.model.norm2(self.a),rhs) |
| def test_noisy_repair_and_orientation(self): |
| noisy=self.a.copy();noisy[0]+=Q(1,100) |
| out=self.model.repair(noisy,'exact') |
| self.assertGreater(out['cost'],0) |
| self.assertEqual(self.ops['C']*out['repaired'],s.zeros(self.ops['C'].rows,1)) |
| for E in self.ops['reconstruction']: |
| cert=no_flip_certificate(s.eye(3),jacobian(E,out['repaired'])) |
| self.assertTrue(cert['verified']) |
| def test_certificate_rejects_flip(self): |
| self.assertFalse(no_flip_certificate(s.eye(3),s.diag(-1,1,1))['verified']) |
| def test_certificate_is_sufficient_not_necessary(self): |
| self.assertFalse(no_flip_certificate(s.eye(3),2*s.eye(3))['verified']) |
| self.assertTrue(monotonicity_certificate(s.eye(3),2*s.eye(3))['verified']) |
| def test_accretive_shear(self): |
| J=s.Matrix([[1,2,0],[-2,1,0],[0,0,1]]) |
| self.assertTrue(monotonicity_certificate(s.eye(3),J,Q(1,2))['verified']) |
| self.assertGreater(J.det(),0) |
| def test_degenerate_tetra_rejected(self): |
| with self.assertRaises(ValueError): |
| tetrahedral_operators([[0,0,0],[1,0,0],[0,1,0],[1,1,0]],[(0,1,2,3)]) |
|
|
| class Moments(unittest.TestCase): |
| def test_group_sizes(self): |
| self.assertEqual(len(signed_permutations()),48) |
| self.assertEqual(len(signed_permutations(True)),24) |
| def test_decoder_80_frames(self): |
| rng=np.random.default_rng(20261005) |
| for _ in range(80): |
| U,_=np.linalg.qr(rng.normal(size=(3,3))) |
| V,_=np.linalg.qr(rng.normal(size=(3,3))) |
| J=U@np.diag(np.exp(rng.uniform(-.6,.6,3)))@V.T |
| moments=encode(J) |
| out=decode(**moments) |
| error=min(np.linalg.norm(out['J']-R@J) for R in signed_permutations()) |
| self.assertLess(error,1e-8) |
| self.assertEqual(out['contractions'],13) |
| def test_inconsistent_sixth_moment_rejected(self): |
| m=encode(np.eye(3));m['L6'][0,0,0,0,0,0]+=.1 |
| with self.assertRaises(ValueError):decode(**m) |
| def test_singular_frame_rejected(self): |
| with self.assertRaises(ValueError):decode(**encode(np.diag([1.,1.,0.]))) |
| def test_invalid_fourth_moment_rejected(self): |
| m=encode(np.eye(3));m['L4'][0,0,0,0]+=.1 |
| with self.assertRaises(ValueError):decode(m['L2'],m['L4']) |
| def test_unique_proper_face_match(self): |
| E=np.array([[1.,0.],[0.,1.],[0.,0.]]) |
| R=np.array([[0,-1,0],[1,0,0],[0,0,1]]) |
| out=match_face(np.eye(3),R,E) |
| self.assertEqual(out['matching_count_within_tolerance'],1) |
| self.assertTrue(np.array_equal(out['R'],R)) |
| both=match_face(np.eye(3),R,E,proper=False) |
| self.assertEqual(both['matching_count_within_tolerance'],2) |
|
|
|
|
| class ExactTensorAdapter(unittest.TestCase): |
| def test_rational_tensor_witness(self): |
| from ahti.moments import exact_moments,decode_rational |
| J=s.Matrix([[1,Q(1,3),0],[0,2,Q(1,5)],[Q(1,7),0,1]]) |
| m=exact_moments(J) |
| out=decode_rational(**m) |
| self.assertEqual(exact_moments(out['J']),m) |
| self.assertEqual(out['status'],'exact_rational_moment_witness') |
| def test_frame_to_coordinate_adapter(self): |
| from ahti.geometry import compile_ordinary_frames |
| X=[[0,0,0],[1,0,0],[0,1,0],[0,0,1],[0,0,-1]] |
| R=s.Matrix([[0,-1,0],[1,0,0],[0,0,1]]) |
| out=compile_ordinary_frames(X,[(0,1,2,3),(0,1,2,4)],[s.eye(3),R]) |
| self.assertEqual(out['period_obstruction'],s.zeros(out['a'].rows,1)) |
| self.assertEqual(out['model'].repair(out['a'])['cost'],0) |
| def test_incompatible_face_rejected(self): |
| from ahti.geometry import compile_ordinary_frames |
| X=[[0,0,0],[1,0,0],[0,1,0],[0,0,1],[0,0,-1]] |
| with self.assertRaises(ValueError): |
| compile_ordinary_frames(X,[(0,1,2,3),(0,1,2,4)],[s.eye(3),s.diag(2,1,1)]) |
|
|
| class IntrinsicMetric(unittest.TestCase): |
| def test_exact_pullback_metric(self): |
| from ahti.moments import exact_moment_metric |
| out=exact_moment_metric(s.Matrix([[1,Q(1,3),0],[0,2,Q(1,5)],[0,0,1]])) |
| D,G=out['differential'],out['ambient_semidefinite_metric'] |
| v=s.Matrix([Q(i-4,7) for i in range(9)]) |
| self.assertEqual((D*v).T*G*(D*v),v.T*v) |
| def test_quotient_distance_scale_equivariance(self): |
| from ahti.moments import exact_quotient_distance2 |
| J=s.eye(3);K=s.Matrix([[1,1,0],[0,1,0],[0,0,2]]) |
| a=exact_quotient_distance2(J,K)['squared_distance'] |
| b=exact_quotient_distance2(3*J,3*K)['squared_distance'] |
| self.assertEqual(b,9*a) |
| def test_quotient_distance_zero_on_same_orbit(self): |
| from ahti.moments import exact_quotient_distance2 |
| R=s.Matrix([[0,-1,0],[1,0,0],[0,0,1]]) |
| self.assertEqual(exact_quotient_distance2(s.eye(3),R)['squared_distance'],0) |
|
|
| if __name__=='__main__':unittest.main(verbosity=2) |
|
|
