mjo/cone/decomposition.py: support non-solid cones

diff --git a/mjo/cone/decomposition.py b/mjo/cone/decomposition.py
index 4a2066783a183b5b8412556ae17ec918f02c0cff..f1a3f76625f37b8fac57dd6fcdfdb85ee111710a 100644 (file)
--- a/mjo/cone/decomposition.py
+++ b/mjo/cone/decomposition.py
@@ -1,26 +1,45 @@
 def irreducible_factors(K):
     r"""
-    Decompose ``K`` into irreducible factors.
+    Decompose a strictly convex (AKA pointed) convex cone ``K``
+    into nontrivial 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.
+    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 (subcones).
+
+    In the literature, "reducible" and "decomposable" are
+    interchangeable.
+
+    The cone being strictly convex / pointed is essential to the
+    uniqueness of the decomposition: the plane is a cone, and it can
+    be decomposed into the direct sum of the x-axis and y-axis, but
+    any other pair of (unequal) lines through the origin would work
+    just as well.
+
+    Being solid is less essential. The definition of "direct sum"
+    requires that the ambient space be fully decomposed into a set of
+    subspaces, and if the cone is not solid, we must (to satisfy the
+    definition) manufacture trivial cones to use as the factors in the
+    subspaces orthogonal to the given cone. This destroys the
+    uniqueness of the direct sum decomposition in the sense that the
+    vector subspaces corresponding to the trivial factors are
+    arbitrary. If the given cone is one-dimensional and the ambient
+    space three-dimensional, then a full decomposition would include
+    two trivial cones in any two one-dimensional subspaces of the
+    plane. As in the preceding paragraph, those one-dimensional
+    subspaces can be chosen arbitrarily.
+
+    Considering that this method does not return the ambient vector
+    space decomposition, adding trivial cones to the list of
+    irreducible factors is not helpful: any cone is equal to the sum
+    of itself with a trivial cone, or two trivial cones, etc. This is
+    all to say: this method will accept non-solid cones, but it will
+    not return any trivial factors. The result does not technically
+    correspond to a direct sum, but by omitting the trivial factors,
+    we recover uniqueness of the factors in the non-solid case.
 
     OUTPUT:
 
@@ -31,18 +50,28 @@ def irreducible_factors(K):
 
     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.
+    If a strictly convex / pointed 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 Bremner 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.
 
     REFERENCES:
 
-    ...
+    .. [HausGul2002] Raphael A. Hauser and Osman Guler.
+                     *Self-Scaled Barrier Functions on Symmetric Cones
+                     and Their Classification*.
+                     Foundations of Computational Mathematics 2(2):121-143,
+                     2002. :doi:`10.1007/s102080010022`.
+    .. [BSPRS2014] David Bremner, Mathieu Dutour Sikirić, Dmitrii V.
+                   Pasechnik, Thomas Rehn and Achill Schürmann.
+                   *Computing symmetry groups of polyhedra*.
+                   LMS Journal of Computation and Mathematics. 17(1):565-581,
+                   2014. :doi:`10.1112/S1461157014000400`.
 
     EXAMPLES:
 
@@ -53,10 +82,28 @@ def irreducible_factors(K):
          1-d cone in 3-d lattice N,
          1-d cone in 3-d lattice N}
 
+    We can decompose cones that aren't solid::
+
+        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}
+
     TESTS:
+
+    The error message looks right::
+
+        sage: K = Cone([(1,0),(-1,0)])
+        sage: irreducible_factors(K)
+        Traceback (most recent call last):
+        ...
+        ValueError: cone must be strictly convex (AKA pointed) for its
+        irreducible factors to be well-defined
+
     """
-    if not K.is_proper():
-        raise ValueError("cone must be proper for the irreducible decomposition to be well-defined")
+    if not K.is_strictly_convex():
+        raise ValueError("cone must be strictly convex (AKA pointed) for"
+                         " its irreducible factors to be well-defined")
 
     V = K.ambient_vector_space()
 
@@ -67,15 +114,17 @@ def irreducible_factors(K):
     pivots = M.echelon_form(algorithm="classical").pivots()
     M = M.matrix_from_columns(pivots)
 
-    # The basis
+    # A basis for span(K)
     B = M.columns(copy=False)
 
-    # The ambient space, but now with a basis of extreme rays.
-    W = V.span_of_basis(B)
+    # 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)
 
     # Make a graph with ray -> basis-vector edges where the
     # coefficients are nonzero.
-    vertices = list(range(V.dimension()))
+    vertices = list(range(W.dimension()))
     edges = [ (i,j)
               for i in range(K.nrays())
               for j in W.coordinate_vector(K.ray(i)).nonzero_positions()