From 138b60c0dad476335a09e81642bff241c91932ca Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 30 Dec 2015 23:23:09 -0500
Subject: [PATCH] Very rough first implementation of pointed_decomposition().

---
 mjo/cone/cone.py | 36 ++++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index a5f5f2f..68edeb4 100644
--- 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.
-- 
2.44.2