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]
-
- An example over the integers::
-
- sage: V = ZZ^5
- sage: v1 = V((0, -11, 5, -1, 1))
- sage: v2 = V((-4, 1, -1, -12, -1))
- sage: v3 = V((-2, 1, 1, 38, 1))
- sage: v4 = V((-1, 0, -3, -1, 1))
- sage: v5 = V((-2, 1, 0, 1, 2))
- sage: O = _orthogonalize([v1,v2,v3,v4,v5])
- sage: all( O[i].inner_product(O[j]).is_zero()
- ....: for i in range(len(O))
- ....: for j in range(len(O))
- ....: if not i == j )
- True
-
- """
- 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):
+def nonlineal_part(K):
r"""
- Return the pointed component of this cone along with its
- lineality space.
+ Return the pointed component of this cone.
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.
+ a vector space (its :meth:`linear_subspace`) 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 (directly) 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.
+ generators that live outside of the lineality space, and then
+ removing their lineal parts. The resulting generators (along with
+ those for the lineality space) still generate the original cone
+ because the original generating set can be recovered by adding
+ back the positive multiples of the orthogonal basis for the
+ lineality space.
SETUP::
- sage: from mjo.cone.decomposition import (
- ....: pointed_linspace_decomposition
- ....: )
+ sage: from mjo.cone.decomposition import nonlineal_part
EXAMPLES:
It's difficult to see this function working because Sage (via PPL)
- minimizes the generating sets of its cones by default. In most(?)
+ minimizes the generating sets of its cones by default. In many
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
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 = nonlineal_part(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)
+ sage: K.lines()
+ N(1, 0, 0)
in 3-d lattice N
These are the same rays that we get if we allow the ``Cone``
But we get it right::
- sage: P,L = pointed_linspace_decomposition(K)
+ sage: P = nonlineal_part(K)
sage: P.rays()
N(7, 223, -146, -123)
in 4-d lattice N
- sage: L.rays()
- N( 11, -1, -1, 0),
- N(-11, 1, 1, 0)
- in 4-d lattice N
- sage: P.ray(0)*M(L.ray(0).dense_coefficient_list())
+ sage: P.ray(0)*M(K.lines()[0].dense_coefficient_list())
0
TESTS::
sage: K = random_cone(strictly_convex=False)
sage: lat = K.lattice()
- sage: P,L = pointed_linspace_decomposition(K)
+ sage: P = nonlineal_part(K)
sage: V = lat.vector_space()
sage: all( V(p).inner_product(V(l)).is_zero()
....: for p in P
- ....: for l in L)
+ ....: for l in K.lines())
True
- sage: J = Cone(P.rays() + L.rays(), lattice=lat)
+ sage: linspace_rays = [ c*l for l in K.lines() for c in [-1,1] ]
+ sage: J = Cone(list(P.rays()) + linspace_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.
+ lat = K.lattice()
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)
-
- # Both generating sets should be minimal, but PPL may be able to
- # rescale them (no guarantee we don't have a bunch of rays that
- # can all be divided by two, for example).
- return (Cone(P_ortho_rays, lattice=lat),
- Cone(L_ortho_rays, lattice=lat))
+ linspace_basis = [V(l) for l in K.lines()]
+ if not linspace_basis:
+ # K is pointed
+ return K
+ linspace = K.linear_subspace() # uses cached lines()
+
+ # WARNING: this assumes that K.lines() is an orthogonal set!
+ P_rays = [
+ _orthogonalize_wrt(v, linspace_basis)
+ for r in K.rays()
+ if not (v := V(r)) in linspace
+ ]
+
+ # The generating set should be minimal, but PPL may be able to
+ # rescale it: there's no guarantee we don't have a bunch of rays
+ # that can all be divided by two, for example.
+ from sage.geometry.cone import Cone
+ return (Cone(P_rays, lattice=lat))
def irreducible_factors(K):
projects / sage.d.git / commitdiff