mjo/cone/decomposition.py: finish nonlineal_part() updates

This was remarkably close to working for something untested; only a
few minor tweaks and typos remained.

diff --git a/mjo/cone/decomposition.py b/mjo/cone/decomposition.py
index 61c3994f5ef1b0f04060ce203cdad2e29d334fd9..3d04c8dba929d40ea9807bdae4aae19a178767c5 100644 (file)
--- a/mjo/cone/decomposition.py
+++ b/mjo/cone/decomposition.py
@@ -205,6 +205,7 @@ def nonlineal_part(K, orthogonal=False):
 
         sage: K = random_cone(strictly_convex=False)
         sage: nonlineal_part(K).is_strictly_convex()
+        True
         sage: nonlineal_part(K, orthogonal=True).is_strictly_convex()
         True
         sage: V = K.lattice().vector_space()
@@ -275,19 +276,26 @@ def nonlineal_part(K, orthogonal=False):
         Vector space of degree 4 and dimension 1 over Rational Field
         Basis matrix:
         [0 0 0 1]
-        sage: P = nonlinear_part(J)
+        sage: P = nonlineal_part(J)
         sage: list(P.rays()) == [g1,g2,g3,g4]
         True
-        sage: L.intersection(span(J.rays())).dimension()
+        sage: L.intersection(V.span(J.rays())).dimension()
         1
 
     """
+    if len(K.lines()) == 0:
+        # K is pointed. We uses len(lines) instead of is_pointed()
+        # because lines() caches the line set that we may need later.
+        return K
+
+    from sage.geometry.cone import Cone
     lat = K.lattice()
 
     if not orthogonal:
         # Open question: is this generating set guaranteed to be
         # minimal?
-        return Cone([r for r in K if not r in K.lines()], lattice=lat)
+        pm_lines = [c*l for c in [-1,1] for l in K.lines()]
+        return Cone([r for r in K if not r in pm_lines], lattice=lat)
 
     # We need to compute inner products and doing it in the ambient
     # vector space (rather than in a more picky lattice) is the
@@ -309,5 +317,4 @@ def nonlineal_part(K, orthogonal=False):
     # The generating set should be minimal, but PPL may be able to
     # rescale it: there's no guarantee we don't have a bunch of rays
     # that can all be divided by two, for example.
-    from sage.geometry.cone import Cone
     return (Cone(P_rays, lattice=lat))