mjo/cone: begin working on the irreducible decomposition

diff --git a/mjo/cone/irreducible_decomposition.py b/mjo/cone/irreducible_decomposition.py
new file mode 100644 (file)
index 0000000..557f6d0
--- /dev/null
+++ b/mjo/cone/irreducible_decomposition.py
@@ -0,0 +1,55 @@
+def irreducible_decomposition(K):
+    r"""
+    Decompose ``K`` into irreducible factors.
+
+    A convex cone is said to be *reducible* if it can be expressed as
+    the direct sum of two subcones, each of which is nontrivial.  An
+    *irreducible* cone, then, is one that is not reducible. Every
+    proper cone can be uniquely decomposed (up to isomorphism and the
+    order of the factors) in to a direct sum of irreducible
+    cones. Furthermore, any convex cone (even an improper one) can be
+    decomposed into a sum of its lineality space, its orthogonal
+    complement, and a unique proper cone, each of which is itself a
+    closed convex cone. In this manner, every closed convex cone can
+    be decomposed into,
+
+    * A linear subspace contained in the cone (its lineality space)
+    * A subspace of the ambient space that isorthogonal to the cone
+    * A collection of irreducible proper cones
+
+    The lineality space and orthogonal complement are unique, but can
+    not be decomposed uniquely into irreducible factors, since any
+    choice of basis for either of these spaces would lead to a
+    different decomposition as a sum of rays.
+
+    OUTPUT:
+
+    A set of ``sage.geometry.cone.ConvexRationalPolyhedralCone``, each
+    of which is nontrivial and irreducible. Each factor lives in the
+    same ambient space as ``K``, rather than in its own span, to avoid
+    confusing isomorphic factors that live in different subspaces.
+
+    EXAMPLES:
+
+        sage: J = cones.nonnegative_orthant(3)
+        sage: K = J.cartesian_product(J)
+
+
+    TESTS:
+    """
+    if not K.is_proper():
+        raise ValueError("cone must be proper for the irreducible decomposition to be well-defined")
+
+    V = K.ambient_vector_space()
+    n = V.dimension()
+
+    # Create a column matrix from the generators of K, and then
+    # perform Gaussian elimination so that the resulting pivots
+    # specify a linearly-independent subset of the original columns.
+    M = K.rays().column_matrix().change_ring(V.base_ring())
+    pivots = M.echelon_form(algorithm="classical").pivots()
+    M = M.matrix_from_columns(pivots)
+
+    # The basis
+    B = M.columns(copy=False)
+