--- /dev/null
+def rational_isomorphism_generator(K1, K2):
+ r"""
+ Generate invertible linear transformations mapping ``K1``
+ onto ``K2``.
+
+ If ``K1`` and ``K2`` are both pointed and solid in
+ ``n``-dimensional space, then they have extreme rays, and we
+ construct isomorphisms (over the rationals) by mapping size-``n``
+ subsets of ``K1``'s extreme rays to size-``n`` subsets of ``K2``'s
+ extreme rays. The resulting isomorphisms are unique only up to
+ scalar multiples; it may be possible to scale some rows of the
+ isomorphism by different positive multiples, and it is always
+ possible to scale the entire isomorphism by a single positive
+ amount. For example, the automorphism group of the nonnegative
+ orthant is allowed to scale each extreme ray by a different
+ multiple, but the automorphisms of the little-l1 cone (a square
+ raised to height one in ``RR^3``) must be scaled as a whole. In
+ general this comes down to the number of extreme rays being
+ greater than the dimension of the cone.
+
+ If ``K1`` and ``K2`` are not pointed, then in addition to the
+ procedure above, we also map a basis for the lineality space of
+ ``K1`` to a basis of the lineality space for ``K2``. Similarly,
+ when ``K1`` or ``K2`` are not solid, we map bases for their
+ orthogonal complements to one another. This procedure is not
+ unique, since there are many isomorphisms (corresponding to many
+ choices of bases) between two vector subspaces of the same
+ size. As a result you should expect a greater degree of
+ non-uniqueness in the isomorphisms when your cones are not pointed
+ or not solid. (If, say, only ``K1`` is not pointed, then we know
+ immediately that ``K2`` cannot be isomorphic to it.)
+
+ The returned matrices are immutable so that they may be hashed,
+ and, for example, de-duplicated using ``set()``.
+
+ EXAMPLES:
+
+ In this example, the lower half-space is obtained from the upper
+ half space by flipping the y-coordinate::
+
+ sage: K1 = Cone([(0,1),(1,0),(-1,0)])
+ sage: K2 = Cone([(0,-1),(1,0),(-1,0)])
+ sage: next(rational_isomorphism_generator(K1,K2))
+ [ 1 0]
+ [ 0 -1]
+
+ Two different descriptions of the plane should be isomorphic::
+
+ sage: K1 = Cone([(1,0), (-1,0), (0,1), (0,-1)])
+ sage: K2 = Cone([(1,1), (-1,1), (0,-1)])
+ sage: K1.is_equivalent(K2)
+ True
+ sage: next(rational_isomorphism_generator(K1,K2))
+ [1 0]
+ [0 1]
+
+ A randomly-generated example constructed to be isomorphic::
+
+ sage: K1 = Cone( [(17, 55, -13, 0),
+ ....: (-9, -32, 914, 0),
+ ....: ( 1, 16, -2, 7),
+ ....: (-1, -16, 2, -7)])
+ sage: A = matrix(QQ, [[ 1, 4, 7, 43],
+ ....: [ 0, 1, 3, 16],
+ ....: [ 2, 6, 9, 58],
+ ....: [ -1, -7, -16, -90]])
+ sage: K2 = Cone([r*A for r in K1])
+ sage: g = next(rational_isomorphism_generator(K1,K2))
+ sage: Cone(r*g for r in K1.rays()).is_equivalent(K2)
+ True
+
+ Automorphisms can be obtained by passing ``K2 == K1``. In this
+ case, there are many duplicates so we use ``set()`` to obtain
+ only the unique automorphisms::
+
+ sage: K1 = Cone([(1,0,1), (-1,0,1), (0,1,1), (0,-1,1)])
+ sage: set(rational_isomorphism_generator(K1,K1))
+ {[-1 0 0]
+ [ 0 -1 0]
+ [ 0 0 1],
+ [-1 0 0]
+ [ 0 1 0]
+ [ 0 0 1],
+ [ 0 -1 0]
+ [-1 0 0]
+ [ 0 0 1],
+ [ 0 -1 0]
+ [ 1 0 0]
+ [ 0 0 1],
+ [ 0 1 0]
+ [-1 0 0]
+ [ 0 0 1],
+ [0 1 0]
+ [1 0 0]
+ [0 0 1],
+ [ 1 0 0]
+ [ 0 -1 0]
+ [ 0 0 1],
+ [1 0 0]
+ [0 1 0]
+ [0 0 1]}
+
+ """
+ # There are no invertible maps between lattices of different
+ # dimensions.
+ if K1.lattice_dim() != K2.lattice_dim():
+ return
+
+ # Create two vector spaces V and W (resp.) to hold K1 and K2,
+ # and find their lineality / perp subspaces.
+ V = K1.lattice().vector_space()
+ L = K1.linear_subspace()
+ K1_perp = K1.span().vector_space().complement()
+ W = K2.lattice().vector_space()
+ M = K2.linear_subspace()
+ K2_perp = K2.span().vector_space().complement()
+
+ if L.dimension() != M.dimension():
+ # Linealities don't match, not isomorphic.
+ return
+
+ if K1_perp.dimension() != K2_perp.dimension():
+ # Cone dimensions don't match, not isomorphic.
+ return
+
+ # Standard linear algebra 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()
+ M_proj = B.T * (B*B.T).inverse() * B
+
+ # Project onto the lineality space and subtract to obtain the
+ # "pointed parts"...
+ 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.
+ 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
+ # 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.
+ 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()
+ R = K1.lattice().vector_space().base_ring() # QQ
+
+ # If the non-lineal parts of our cones are dimension zero, both
+ # cones are subspaces of the same dimension, and the identity map
+ # 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():
+ X = matrix.identity(R, V.dimension())
+ X.set_immutable()
+ yield X
+ return
+
+ # 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()
+
+ # 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)):
+ # 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
+ # actually an isomorphism.
+ A = matrix(R, list(rs1) + K1_extra_rows)
+ B = matrix(R, list(rs2) + K2_extra_rows)
+ try:
+ X = A.solve_right(B, extend=False)
+ except ValueError:
+ continue
+
+ if Cone(r*X for r in K1.rays()).is_equivalent(K2):
+ X.set_immutable()
+ yield X
projects / sage.d.git / commitdiff