mjo/cone/irreducible_decomposition.py: get it working, roughly

diff --git a/mjo/cone/irreducible_decomposition.py b/mjo/cone/irreducible_decomposition.py
index 557f6d00b1ba6282422e6cbe345414e696cfb28f..16d2f865455a54108259f534b0f9cba74ea3dff3 100644 (file)
--- a/mjo/cone/irreducible_decomposition.py
+++ b/mjo/cone/irreducible_decomposition.py
@@ -29,11 +29,29 @@ def irreducible_decomposition(K):
     same ambient space as ``K``, rather than in its own span, to avoid
     confusing isomorphic factors that live in different subspaces.
 
-    EXAMPLES:
+    ALGORITHM:
+
+    If a proper cone ``K`` is reducible to ``K1 + K2`` in the vector
+    space ``V = V1 + V2``, then the generators (extreme rays) of ``K``
+    can be split into subsets ``G1`` and ``G2`` of ``V1`` and ``V2``
+    that generate ``K1`` and ``K2``, respectively. Following Pasechnik
+    et al., we find a basis of the ambient space consisting of extreme
+    rays of ``K``, and then express each extreme ray in terms of that
+    basis. By looking for cliques among those coordinates, we can
+    determine which, if any, generators can be split accordingly.
 
-        sage: J = cones.nonnegative_orthant(3)
-        sage: K = J.cartesian_product(J)
+    REFERENCES:
 
+    ...
+
+    EXAMPLES:
+
+        sage: from mjo.cone.irreducible_decomposition import *
+        sage: K = cones.nonnegative_orthant(3)
+        sage: irreducible_decomposition(K)
+        {1-d cone in 3-d lattice N,
+         1-d cone in 3-d lattice N,
+         1-d cone in 3-d lattice N}
 
     TESTS:
     """
@@ -41,7 +59,6 @@ def irreducible_decomposition(K):
         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
@@ -53,3 +70,18 @@ def irreducible_decomposition(K):
     # The basis
     B = M.columns(copy=False)
 
+    # The ambient space, but now with a basis of extreme rays.
+    W = V.span_of_basis(B)
+
+    # Make a graph with ray -> basis-vector edges where the
+    # coefficients are nonzero.
+    vertices = list(range(V.dimension()))
+    edges = [ (i,j)
+              for i in range(K.nrays())
+              for j in W.coordinate_vector(K.ray(i)).nonzero_positions()
+              if j != i ]
+    from sage.graphs.graph import Graph
+    G = Graph([vertices, edges], format='vertices_and_edges')
+
+    from sage.geometry.cone import Cone
+    return { Cone(K.rays(c)) for c in G.connected_components() }