mjo/cone/decomposition.py: add tests, fix the edges of the edges

A miscellaneous selection:

  * The vertex numbering was off because pivots() renumbers them.
  * Updated the docs to fix some incorrect claims.
  * Added lots of new tests.
  * Special case for irreducible cones to return K itself.
  * Use subspace_with_basis() instead of span() to preserve
    our user basis.

Tests now pass, and the documentation is converging upon the truth.

diff --git a/mjo/cone/decomposition.py b/mjo/cone/decomposition.py
index f1a3f76625f37b8fac57dd6fcdfdb85ee111710a..5375f15562764b8be285dbe14a5240262a006e48 100644 (file)
--- a/mjo/cone/decomposition.py
+++ b/mjo/cone/decomposition.py
@@ -7,7 +7,7 @@ def irreducible_factors(K):
     the direct sum of two subcones. An *irreducible* cone is a cone
     that is not reducible. Every cone that :meth:`is_proper` can be
     uniquely decomposed -- up to isomorphism and the order of the
-    factors -- as a direct sum of irreducible, proper, nontrivial
+    factors -- as a direct sum of irreducible, pointed, nontrivial
     factors (subcones).
 
     In the literature, "reducible" and "decomposable" are
@@ -44,9 +44,11 @@ def irreducible_factors(K):
     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.
+    of which is nontrivial, strictly convex (pointed), 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. As a consequence, factors of
+    solid cones will not in general be solid.
 
     ALGORITHM:
 
@@ -54,11 +56,11 @@ def irreducible_factors(K):
     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 Bremner et al., we find a basis of the
-    ambient space consisting of extreme rays of ``K``, and then
+    respectively. Following Bremner et al., we find a basis for the
+    span of ``K`` 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.
+    groups of rays that require a common subset of basis elements, we
+    determine which, if any, generators can be split accordingly.
 
     REFERENCES:
 
@@ -73,21 +75,79 @@ def irreducible_factors(K):
                    LMS Journal of Computation and Mathematics. 17(1):565-581,
                    2014. :doi:`10.1112/S1461157014000400`.
 
+    SETUP:
+
+        sage: from mjo.cone.decomposition import irreducible_factors
+
     EXAMPLES:
 
-        sage: from mjo.cone.decomposition import *
+    The nonnegative orthant is a direct sum of its extreme rays::
+
         sage: K = cones.nonnegative_orthant(3)
-        sage: irreducible_factors(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}
+        sage: expected = {Cone([r]) for r in K.rays()}
+        sage: irreducible_factors(K) == expected
+        True
 
-    We can decompose cones that aren't solid::
+    An irreducible example (the l1-norm cone)::
 
-        sage: K = Cone([(1,0,0),(0,1,0)])
-        sage: irreducible_factors(K)
-        {1-d cone in 3-d lattice N,
-         1-d cone in 3-d lattice N}
+        sage: K = Cone([(1,0,1), (0,1,1), (-1,0,1), (0,-1,1)])
+        sage: irreducible_factors(K) == {K}
+        True
+
+    We can decompose reducible cones that aren't solid::
+
+        sage: K = Cone([(1,0,0), (0,1,0)])
+        sage: expected = {Cone([r]) for r in K.rays()}
+        sage: irreducible_factors(K) == expected
+        True
+
+    And nothing bad happens with irreducible ones::
+
+        sage: K = Cone([(1,0,1,0), (0,1,1,0), (-1,0,1,0), (0,-1,1,0)])
+        sage: irreducible_factors(K) == {K}
+        True
+
+    The rays of the Schur cone are linearly-independent::
+
+        sage: K = cones.schur(5)
+        sage: expected = {Cone([r]) for r in K.rays()}
+        sage: irreducible_factors(K) == expected
+        True
+
+    The Cartesian product can be used to construct larger
+    examples. Here, two of the factors are irreducible, but the
+    nonnegative orthant should reduce to two rays::
+
+        sage: J1 = Cone([(1,0,1,0), (0,1,1,0), (-1,0,1,0), (0,-1,1,0)])
+        sage: J2 = cones.nonnegative_orthant(2)
+        sage: J3 = cones.rearrangement(2, 5)
+        sage: J1
+        3-d cone in 4-d lattice N
+        sage: J3
+        5-d cone in 5-d lattice N
+        sage: K = J1.cartesian_product(J2).cartesian_product(J3)
+        sage: sorted(irreducible_factors(K))
+        [5-d cone in 11-d lattice N+N+N,
+         1-d cone in 11-d lattice N+N+N,
+         1-d cone in 11-d lattice N+N+N,
+         3-d cone in 11-d lattice N+N+N]
+
+    All proper cones in two dimensions are isomorphic to the
+    nonnegative orthant and should decompose::
+
+        sage: K = Cone([(1,2), (-3,4)])
+        sage: sorted(irreducible_factors(K))
+        [1-d cone in 2-d lattice N, 1-d cone in 2-d lattice N]
+
+    In three dimensions, we can hit the nonnegative orthant with an
+    invertible map, and it will still decompose::
+
+        sage: A = matrix(QQ, [[1, -1/2, 0], [1, 1, -2], [-2, 0, -1]])
+        sage: K = cones.nonnegative_orthant(3)
+        sage: AK = Cone([ r*A for r in K.rays() ])
+        sage: expected = {Cone([r]) for r in AK.rays()}
+        sage: irreducible_factors(AK) == expected
+        True
 
     TESTS:
 
@@ -120,17 +180,22 @@ def irreducible_factors(K):
     # The span of K, but now with a basis of extreme rays. If K is not
     # solid, B will not be a basis of the entire space V, only of a
     # subspace.
-    W = V.span(B)
+    W = V.subspace_with_basis(B)
 
-    # Make a graph with ray -> basis-vector edges where the
-    # coefficients are nonzero.
+    # Make a graph with ray -> ray edges where the coefficients in the
+    # W-representation nonzero. Beware: while they ultimately represent
+    # rays of K, the W-coordinates have been renumbered by pivots().
     vertices = list(range(W.dimension()))
-    edges = [ (i,j)
+    edges = [ (i,pivots[j])
               for i in range(K.nrays())
               for j in W.coordinate_vector(K.ray(i)).nonzero_positions()
-              if j != i ]
+              if pivots[j] != i ]
     from sage.graphs.graph import Graph
     G = Graph([vertices, edges], format='vertices_and_edges')
 
     from sage.geometry.cone import Cone
+    if G.connected_components_number() == 1:
+        # Special case where we don't want to pointlessly return an
+        # equivalent but unequal copy of the input cone.
+        return {K}
     return { Cone(K.rays(c)) for c in G.connected_components() }