From 21a3775c72afaab4e0b36e975f7c3498183e4386 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Tue, 1 Jul 2025 09:33:19 -0400 Subject: [PATCH] mjo/cone: improve isomorphism testing & docs, tests pass now --- mjo/cone/isomorphism.py | 123 ++++++++++++++++++++++++++++++++-------- 1 file changed, 99 insertions(+), 24 deletions(-) diff --git a/mjo/cone/isomorphism.py b/mjo/cone/isomorphism.py index f746914..61cea42 100644 --- a/mjo/cone/isomorphism.py +++ b/mjo/cone/isomorphism.py @@ -14,18 +14,19 @@ def linear_isomorphisms(K1, K2): - ``K1`` -- the first cone - - ``K2`` -- the second cone + - ``K2`` -- :class:`sage.geometry.cone.ConvexRationalPolyhedralCone`; + the cone to generate isomorphisms with OUTPUT: A generator that yields linear isomorphisms between ``K1`` and ``K2``, one at a time. The isomorphisms consist of invertible matrices ``A`` whose entries are rational numbers and satisfy - ``Cone(r*A for r in K1).is_equivalent(K2)``. In other words, they - map the cone ``K1`` onto ``K2``, acting on the rays of ``K1`` from - the right. The right-action is mainly for compatibility with the - basis matrix of a vector space, which in SageMath contains the - basis vectors as rows. + ``Cone(K1.rays()*A).is_equivalent(K2)``. In other words, they map + the cone ``K1`` onto ``K2``, acting on the rays of ``K1`` from the + right. The right-action is mainly for compatibility with the ray + matrices of cones and the basis matrices of vector spaces, which + in Sage are "lists of rows." The returned matrices are immutable so that they may be hashed, and, for example, de-duplicated using ``set()``. @@ -64,6 +65,10 @@ def linear_isomorphisms(K1, K2): to vector-space isomorphisms of the lineality spaces and othogonal complements. + SETUP:: + + sage: from mjo.cone.isomorphism import linear_isomorphisms + EXAMPLES: In this example, the lower half-space is obtained from the upper @@ -95,9 +100,9 @@ def linear_isomorphisms(K1, K2): ....: [ 0, 1, 3, 16], ....: [ 2, 6, 9, 58], ....: [ -1, -7, -16, -90]]) - sage: K2 = Cone([r*A for r in K1]) + sage: K2 = Cone(K1.rays()*A, K1.lattice()) sage: g = next(linear_isomorphisms(K1,K2)) - sage: Cone(r*g for r in K1.rays()).is_equivalent(K2) + sage: Cone(K1.rays()*g).is_equivalent(K2) True Automorphisms can be obtained by passing ``K2 == K1``. In this @@ -141,6 +146,30 @@ def linear_isomorphisms(K1, K2): sage: I in linear_isomorphisms(K,K) True + Check the properties of the generated isomorphisms. We build + ``K2`` from ``K1`` using an invertible rational matrix, so we know + that there is at least one isomorphism between them:: + + sage: K1 = random_cone(max_ambient_dim=5) + sage: R = K1.lattice().vector_space().base_ring() + sage: n = K1.lattice_dim() + sage: A = matrix.random(R, n, algorithm='unimodular') + sage: q = QQ.random_element() + sage: while q.is_zero(): + ....: q = QQ.random_element() + sage: K2 = Cone(K1.rays()*A, K1.lattice()) + sage: all( + ....: g.is_invertible() + ....: and + ....: g.is_immutable() + ....: and + ....: Cone(K1.rays()*g).is_equivalent(K2) + ....: and + ....: Cone(K2.rays()*g.inverse()).is_equivalent(K1) + ....: for g in linear_isomorphisms(K1,K2) + ....: ) + True + """ # There are no invertible maps between lattices of different # dimensions. @@ -164,8 +193,7 @@ def linear_isomorphisms(K1, K2): # Cone dimensions don't match, not isomorphic. return - # Standard linear algebra trick for orthogonal projections - # onto L and M. + # Standard trick for orthogonal projections onto L and M. A = L.basis_matrix() L_proj = A.T * (A*A.T).inverse() * A B = M.basis_matrix() @@ -176,28 +204,34 @@ def linear_isomorphisms(K1, K2): L_perps = [ r - r*L_proj for r in K1 ] M_perps = [ r - r*M_proj for r in K2 ] - # which need to be immutable to be put into sets. + # ...which need to be immutable to be put into sets. for i in range(len(L_perps)): L_perps[i].set_immutable() for j in range(len(M_perps)): M_perps[j].set_immutable() - # Our generators are assumed to be a minimal set, so none should - # be eliminated by the Cone(...) below, although some may be + # Our generators are assumed to be minimal, so none should be + # eliminated by the Cone(...) below, although some may be # reordered (which is why we still use L_perps/M_perps directly # even after constructing these cones). In any case, if we wind up # with different numbers of rays here, we aren't going to get any # isomorphisms. This is slow, but it's the most common case (cones - # not isomorphic at all), so probably worth doing before - # proceeding. + # not isomorphic), so probably worth doing before proceeding. from sage.geometry.cone import Cone K1_p = Cone(L_perps, K1.lattice()) K2_p = Cone(M_perps, K2.lattice()) if not K1_p.nrays() == K2_p.nrays(): return - # The dimension of the non-lineal part of our cones. - n = K1.dim() - K1.lineality() + # The dimension of the non-lineal part of our cones. Needs to be a + # python int so that we can pass it to itertools.permutations. + n = int(K1.dim() - K1.lineality()) + + # The base ring of the matrices we'll construct later. We pass it + # explicitly to avoid the rare situation where our rational + # matrices all contain integral entries, making sage think that we + # want an integral matrix. This becomes important when passing + # ``extend=False`` to solve_right(). R = K1.lattice().vector_space().base_ring() # QQ # If the non-lineal parts of our cones are dimension zero, both @@ -205,7 +239,7 @@ def linear_isomorphisms(K1, K2): # between them should work. This special case is needed to avoid # solving linear equations with the empty set below. from sage.matrix.constructor import matrix - if n.is_zero(): + if n == 0: X = matrix.identity(R, V.dimension()) X.set_immutable() yield X @@ -214,15 +248,15 @@ def linear_isomorphisms(K1, K2): # Compute these outside of the loop. These are extra rows that we # append the the system of equations to ensure that the lineality # and perp spaces of K1 get mapped to those of K2. - K1_extra_rows = L.basis_matrix().rows() + K1_perp.basis_matrix().rows() - K2_extra_rows = M.basis_matrix().rows() + K2_perp.basis_matrix().rows() + K1_extra_rows = A.rows() + K1_perp.basis_matrix().rows() + K2_extra_rows = B.rows() + K2_perp.basis_matrix().rows() # Now try to find isomorphisms in the first component (non-lineal # part), and extend them to the whole thing. Use lists to keep # everything in the right order. from itertools import permutations - for rs1 in permutations(L_perps, int(n)): - for rs2 in permutations(M_perps, int(n)): + for rs1 in permutations(L_perps, n): + for rs2 in permutations(M_perps, n): # Try mapping every subset of the right size in the domain # to every subset of the right size in the codomain. Extend # to the lineal part, and then check that the result is @@ -234,7 +268,7 @@ def linear_isomorphisms(K1, K2): except ValueError: continue - if Cone(r*X for r in K1.rays()).is_equivalent(K2): + if Cone(K1.rays()*X).is_equivalent(K2): X.set_immutable() yield X @@ -244,6 +278,47 @@ def is_linearly_isomorphic(K1,K2): Return ``True`` if ``K1`` and ``K2`` are linearly isomorphic over the rationals, and ``False`` otherwise. + INPUT: + + - ``K1`` -- the first cone + + - ``K2`` -- :class:`sage.geometry.cone.ConvexRationalPolyhedralCone`; + the cone to check for isomorphism with ``K1`` + + OUTPUT: + + ``True`` if there exists an invertible linear matrix with rational + entries mapping ``K1`` to ``K2``, and ``False`` otherwise. + + ALGORITHM: + + We attempt to construct an explicit isomorphism using + :meth:`linear_isomorphisms`. For more information, see the + documentation of that method. + + SETUP:: + + sage: from mjo.cone.isomorphism import is_linearly_isomorphic + + EXAMPLES: + + All simplicial cones with the same number of rays are isomorphic:: + + sage: K1 = random_cone(max_ambient_dim=5) + sage: while not K1.is_simplicial(): + ....: K1 = random_cone(max_ambient_dim=5) + sage: n = K1.lattice_dim() + sage: r = K1.nrays() + sage: K2 = random_cone(min_ambient_dim=n, + ....: max_ambient_dim=n, + ....: min_rays=r, + ....: max_rays=r, + ....: strictly_convex=True) + sage: K2.is_simplicial() + True + sage: is_linearly_isomorphic(K1,K2) + True + TESTS: Every cone is isomorphic to itself:: @@ -259,7 +334,7 @@ def is_linearly_isomorphic(K1,K2): ....: max_ambient_dim=K.lattice_dim(), ....: min_rays=K.nrays(), ....: max_rays=K.nrays()) - sage: K.is_isomorphic(J) == J.is_isomorphic(K) + sage: is_linearly_isomorphic(K,J) == is_linearly_isomorphic(K,J) True """ -- 2.49.0