Very rough first implementation of pointed_decomposition().

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index a5f5f2f4fcf1ac603a9d8d48a75a0924003bd8cf..68edeb4a012b1bca7c4f9664970d2d89371b4926 100644 (file)
--- a/mjo/cone/cone.py
+++ b/mjo/cone/cone.py
@@ -141,6 +141,42 @@ def random_element(K):
     return V(sum(scaled_gens))
 
 
+def pointed_decomposition(K):
+    """
+    Every convex cone is the direct sum of a pointed cone and a linear
+    subspace. Return a pair ``(P,S)`` of cones such that ``P`` is
+    pointed, ``S`` is a subspace, and ``K`` is the direct sum of ``P``
+    and ``S``.
+
+    OUTPUT:
+
+    An ordered pair ``(P,S)`` of closed convex polyhedral cones where
+    ``P`` is pointed, ``S`` is a subspace, and ``K`` is the direct sum
+    of ``P`` and ``S``.
+
+    TESTS:
+
+        sage: set_random_seed()
+        sage: K = random_cone(max_ambient_dim=8)
+        sage: (P,S) = pointed_decomposition(K)
+        sage: x = random_element(K)
+        sage: P.contains(x) or S.contains(x)
+        True
+        sage: x.is_zero() or (P.contains(x) != S.contains(x))
+        True
+    """
+    linspace_gens  = [ copy(b) for b in K.linear_subspace().basis() ]
+    linspace_gens += [ -b for b in linspace_gens ]
+
+    S = Cone(linspace_gens, K.lattice())
+
+    # Since ``S`` is a subspace, its dual is its orthogonal complement
+    # (albeit in the wrong lattice).
+    S_perp = Cone(S.dual(), K.lattice())
+    P = K.intersection(S_perp)
+
+    return (P,S)
+
 def positive_operator_gens(K):
     r"""
     Compute generators of the cone of positive operators on this cone.