+def _orthogonalize_wrt(r, S):
+ r"""
+ Orthogonalize the vector ``r`` with respect to the list ``S``.
+
+ We make ``r`` orthogonal to every element of ``S``, assuming that the
+ vectors in ``S`` are already orthogonal to one another.
+
+ SETUP::
+
+ sage: from mjo.cone.decomposition import _orthogonalize_wrt
+
+ EXAMPLES::
+
+ sage: V = QQ^3
+ sage: S = [ V((1,0,0)), V((0,1,0)) ]
+ sage: r = V((1,1,1))
+ sage: _orthogonalize_wrt(r, S)
+ (0, 0, 1)
+ sage: r = V((2,-3,7))
+ sage: _orthogonalize_wrt(r, S)
+ (0, 0, 7)
+
+ """
+ from math import prod
+ t = r*prod(s.inner_product(s) for s in S)
+ t -= sum( prod( p.inner_product(p) for p in S if p != s)
+ * r.inner_product(s) * s
+ for s in S )
+ return t
+
+def _orthogonalize(S):
+ r"""
+ Orthogonalize the set ``S`` of vectors over the rationals.
+
+ Go one at a time. The first one forms an orthogonal set, then the
+ second one gets orthogonalized with respect to that, and so on.
+
+ SETUP::
+
+ sage: from mjo.cone.decomposition import _orthogonalize
+
+ EXAMPLES::
+
+ sage: A = matrix(QQ, 3, [[-4, -16, 71],
+ ....: [-2, -5, 21],
+ ....: [5, 19, -84]])
+ sage: S = A.columns()
+ sage: # compute the gram matrix of S
+ sage: matrix(QQ, 3, lambda i,j: S[i].inner_product(S[j]))
+ [ 45 169 -746]
+ [ 169 642 -2837]
+ [ -746 -2837 12538]
+ sage: O = _orthogonalize(S); O
+ [(-4, -2, 5), (-44, 113, 10), (-26325, -8100, -24300)]
+ sage: # compute the gram matrix of "O" now
+ sage: matrix(QQ, 3, lambda i,j: O[i].inner_product(O[j]))
+ [ 45 0 0]
+ [ 0 14805 0]
+ [ 0 0 1349105625]
+
+ """
+ if len(S) == 0:
+ return S
+ result = [S[0]]
+ for s in S[1:]:
+ s_ortho = _orthogonalize_wrt(s, result)
+ if not s_ortho.is_zero():
+ result.append(s_ortho)
+ return result
+
+
+def pointed_linspace_decomposition(K):
+ r"""
+ Return the pointed component of this cone along with its
+ lineality space.
+
+ Every closed convex cone can be decomposed into the direct sum of
+ a vector space (its lineality space, which is also a closed convex
+ cone) and a pointed closed convex cone that lives in the
+ orthogonal complement of the lineality space. The straightforward
+ way to obtain the generators of the pointed component is through
+ orthogonal projection, but we can't orthogonally project over the
+ rationals.
+
+ This function takes a different approach, positively scaling the
+ generators of the pointed component and then removing their lineal
+ parts. The resulting generators still generate the original cone,
+ but can (by construction) be separated into two sets corresponding
+ to the lineality space and its orthogonal complement. This is a
+ direct sum decomposition of the original cone.
+
+ SETUP::
+
+ sage: from mjo.cone.decomposition import (
+ ....: pointed_linspace_decomposition
+ ....: )
+
+ EXAMPLES:
+
+ It's difficult to see this function working because Sage (via PPL)
+ minimizes the generating sets of its cones by default. In most(?)
+ cases, this minimization process will do the same thing that this
+ method does, subtracting out the lineal portions of the generators
+ that do not live in the lineality space. Here is one example of
+ the ell-one cone living above a line::
+
+ sage: N = ToricLattice(3)
+ sage: G = [ (1,0,0), (-1,0,0),
+ ....: (1,0,1), (-1,0,1), (0,1,1), (0,-1,1) ]
+ sage: K = Cone(list(map(N,G)), check=False, normalize=False)
+ sage: P,L = pointed_linspace_decomposition(K)
+ sage: P.rays()
+ N(0, -1, 1),
+ N(0, 1, 1)
+ in 3-d lattice N
+ sage: L.rays()
+ N( 1, 0, 0),
+ N(-1, 0, 0)
+ in 3-d lattice N
+
+ These are the same rays that we get if we allow the ``Cone``
+ constructor to reduce the generating set::
+
+ sage: Cone(G).rays()
+ N( 0, 1, 1),
+ N( 0, -1, 1),
+ N( 1, 0, 0),
+ N(-1, 0, 0)
+ in 3-d lattice N
+
+ TESTS::
+
+ sage: K = random_cone(strictly_convex=False)
+ sage: lat = K.lattice()
+ sage: P,L = pointed_linspace_decomposition(K)
+ sage: V = lat.vector_space()
+ sage: all( V(p).inner_product(V(l)).is_zero()
+ ....: for p in P
+ ....: for l in L)
+ True
+ sage: J = Cone(P.rays() + L.rays(), lattice=lat)
+ sage: J.is_equivalent(K)
+ True
+
+ """
+ from sage.geometry.cone import Cone
+
+ L = K.linear_subspace()
+ lat = K.lattice()
+ if L.dimension() == 0:
+ # the lineal part is trivial
+ T = Cone([], lattice=lat)
+ return (K,T)
+
+ # We need to compute inner products and doing it in the ambient
+ # vector space (rather than in a more picky lattice) is the
+ # easiest way to do it.
+ V = lat.vector_space()
+
+ L_rays = []
+ P_rays = []
+ for r in K.rays():
+ v = V(r)
+ if v in L:
+ L_rays.append(v)
+ else:
+ P_rays.append(v)
+
+ # The set of rays that generate L can contain up to 2*dim(L)
+ # rays, but of course we need only dim(L) of them when they
+ # are orthogonal; the other dim(L) can then be obtained by
+ # negation.
+ L_ortho_basis = [ r for r in _orthogonalize(L_rays)
+ if not r.is_zero() ]
+ L_ortho_rays = [ c*b for c in [-1,1] for b in L_ortho_basis ]
+
+ P_ortho_rays = []
+ for r in P_rays:
+ p = _orthogonalize_wrt(r, L_ortho_basis)
+ P_ortho_rays.append(p)
+
+ return (Cone(P_ortho_rays, lattice=lat),
+ Cone(L_ortho_rays, lattice=lat))
+
+
def irreducible_factors(K):
r"""
Decompose a strictly convex (AKA pointed) convex cone ``K``
projects / sage.d.git / commitdiff