Remove unused code and implement the improved Lyapunov rank algorithm.

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index 60f9c34ec8bc271d65812859f51ca77636c8cbbc..421fb3c8f04abb23fd420271ed10e0a7e65815d4 100644 (file)
--- a/mjo/cone/cone.py
+++ b/mjo/cone/cone.py
@@ -7,6 +7,7 @@ addsitedir(abspath('../../'))
 
 from sage.all import *
 
+
 def project_span(K):
     r"""
     Project ``K`` into its own span.
@@ -32,89 +33,30 @@ def project_span(K):
 
     The projected cone should always be solid::
 
-        sage: K = random_cone()
+        sage: K = random_cone(max_dim = 10)
         sage: K_S = project_span(K)
         sage: K_S.is_solid()
         True
 
     If we do this according to our paper, then the result is proper::
 
-        sage: K = random_cone()
+        sage: K = random_cone(max_dim = 10)
         sage: K_S = project_span(K)
         sage: P = project_span(K_S.dual()).dual()
         sage: P.is_proper()
         True
 
     """
-    F = K.lattice().base_field()
-    Q = K.lattice().quotient(K.sublattice_complement())
+    L = K.lattice()
+    F = L.base_field()
+    Q = L.quotient(K.sublattice_complement())
     vecs = [ vector(F, reversed(list(Q(r)))) for r in K.rays() ]
 
-    L = None
+    newL = None
     if len(vecs) == 0:
-        L = ToricLattice(0)
-
-    return Cone(vecs, lattice=L)
-
-
-def rename_lattice(L,s):
-    r"""
-    Change all names of the given lattice to ``s``.
-    """
-    L._name = s
-    L._dual_name = s
-    L._latex_name = s
-    L._latex_dual_name = s
-
-def span_iso(K):
-    r"""
-    Return an isomorphism (and its inverse) that will send ``K`` into a
-    lower-dimensional space isomorphic to its span (and back).
-
-    EXAMPLES:
-
-    The inverse composed with the isomorphism should be the identity::
-
-        sage: K = random_cone(max_dim=10)
-        sage: (phi, phi_inv) = span_iso(K)
-        sage: phi_inv(phi(K)) == K
-        True
-
-    The image of ``K`` under the isomorphism should have full dimension::
-
-        sage: K = random_cone(max_dim=10)
-        sage: (phi, phi_inv) = span_iso(K)
-        sage: phi(K).dim() == phi(K).lattice_dim()
-        True
-
-    """
-    phi_domain = K.sublattice().vector_space()
-    phi_codo = VectorSpace(phi_domain.base_field(), phi_domain.dimension())
-
-    # S goes from the new space to the cone space.
-    S = linear_transformation(phi_codo, phi_domain, phi_domain.basis())
-
-    # phi goes from the cone space to the new space.
-    def phi(J_orig):
-        r"""
-        Takes a cone ``J`` and sends it into the new space.
-        """
-        newrays = map(S.inverse(), J_orig.rays())
-        L = None
-        if len(newrays) == 0:
-            L = ToricLattice(0)
+        newL = ToricLattice(0)
 
-        return Cone(newrays, lattice=L)
-
-    def phi_inverse(J_sub):
-        r"""
-        The inverse to phi which goes from the new space to the cone space.
-        """
-        newrays = map(S, J_sub.rays())
-        return Cone(newrays, lattice=K.lattice())
-
-
-    return (phi, phi_inverse)
+    return Cone(vecs, lattice=newL)
 
 
 
@@ -365,6 +307,15 @@ def lyapunov_rank(K):
         sage: lyapunov_rank(octant)
         3
 
+    The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
+    [Orlitzky/Gowda]_::
+
+        sage: R5 = VectorSpace(QQ, 5)
+        sage: gens = R5.basis() + [ -r for r in R5.basis() ]
+        sage: K = Cone(gens)
+        sage: lyapunov_rank(K)
+        25
+
     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
     [Rudolf et al.]_::
 
@@ -378,7 +329,30 @@ def lyapunov_rank(K):
         sage: lyapunov_rank(L3infty)
         1
 
-    The Lyapunov rank should be additive on a product of cones
+    A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
+    + 1` [Orlitzky/Gowda]_::
+
+        sage: K = Cone([(1,0,0,0,0)])
+        sage: lyapunov_rank(K)
+        21
+        sage: K.lattice_dim()**2 - K.lattice_dim() + 1
+        21
+
+    A subspace (of dimension `m`) in `n` dimensions should have a
+    Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
+
+        sage: e1 = (1,0,0,0,0)
+        sage: neg_e1 = (-1,0,0,0,0)
+        sage: e2 = (0,1,0,0,0)
+        sage: neg_e2 = (0,-1,0,0,0)
+        sage: zero = (0,0,0,0,0)
+        sage: K = Cone([e1, neg_e1, e2, neg_e2, zero, zero, zero])
+        sage: lyapunov_rank(K)
+        19
+        sage: K.lattice_dim()**2 - K.dim()*(K.lattice_dim() - K.dim())
+        19
+
+    The Lyapunov rank should be additive on a product of proper cones
     [Rudolf et al.]_::
 
         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
@@ -404,11 +378,11 @@ def lyapunov_rank(K):
 
     TESTS:
 
-    The Lyapunov rank should be additive on a product of cones
+    The Lyapunov rank should be additive on a product of proper cones
     [Rudolf et al.]_::
 
-        sage: K1 = random_cone(max_dim=10, max_rays=10)
-        sage: K2 = random_cone(max_dim=10, max_rays=10)
+        sage: K1 = random_cone(max_dim=10, strictly_convex=True, solid=True)
+        sage: K2 = random_cone(max_dim=10, strictly_convex=True, solid=True)
         sage: K = K1.cartesian_product(K2)
         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
         True
@@ -446,79 +420,42 @@ def lyapunov_rank(K):
     The calculation of the Lyapunov rank of an improper cone can be
     reduced to that of a proper cone [Orlitzky/Gowda]_::
 
-        sage: K = random_cone(max_dim=15, solid=False, strictly_convex=False)
+        sage: K = random_cone(max_dim=10)
         sage: actual = lyapunov_rank(K)
-        sage: (phi1, _) = span_iso(K)
-        sage: K_S = phi1(K)
-        sage: (phi2, _) = span_iso(K_S.dual())
-        sage: J_T = phi2(K_S.dual()).dual()
+        sage: K_S = project_span(K)
+        sage: P = project_span(K_S.dual()).dual()
         sage: l = K.linear_subspace().dimension()
         sage: codim = K.lattice_dim() - K.dim()
-        sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
+        sage: expected = lyapunov_rank(P) + K.dim()*(l + codim) + codim**2
         sage: actual == expected
         True
 
-    Repeat the previous test with different ``random_cone()`` params::
+    The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
 
-        sage: K = random_cone(max_dim=15, solid=False, strictly_convex=True)
-        sage: actual = lyapunov_rank(K)
-        sage: (phi1, _) = span_iso(K)
-        sage: K_S = phi1(K)
-        sage: (phi2, _) = span_iso(K_S.dual())
-        sage: J_T = phi2(K_S.dual()).dual()
-        sage: l = K.linear_subspace().dimension()
-        sage: codim = K.lattice_dim() - K.dim()
-        sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
-        sage: actual == expected
+        sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
+        sage: lyapunov_rank(K) == len(LL(K))
         True
 
-        sage: K = random_cone(max_dim=15, solid=True, strictly_convex=False)
-        sage: actual = lyapunov_rank(K)
-        sage: (phi1, _) = span_iso(K)
-        sage: K_S = phi1(K)
-        sage: (phi2, _) = span_iso(K_S.dual())
-        sage: J_T = phi2(K_S.dual()).dual()
-        sage: l = K.linear_subspace().dimension()
-        sage: codim = K.lattice_dim() - K.dim()
-        sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
-        sage: actual == expected
-        True
+    """
+    beta = 0
 
-        sage: K = random_cone(max_dim=15, solid=True, strictly_convex=True)
-        sage: actual = lyapunov_rank(K)
-        sage: (phi1, _) = span_iso(K)
-        sage: K_S = phi1(K)
-        sage: (phi2, _) = span_iso(K_S.dual())
-        sage: J_T = phi2(K_S.dual()).dual()
-        sage: l = K.linear_subspace().dimension()
-        sage: codim = K.lattice_dim() - K.dim()
-        sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
-        sage: actual == expected
-        True
+    m = K.dim()
+    n = K.lattice_dim()
+    l = K.linear_subspace().dimension()
 
-        sage: K = random_cone(max_dim=15)
-        sage: actual = lyapunov_rank(K)
-        sage: (phi1, _) = span_iso(K)
-        sage: K_S = phi1(K)
-        sage: (phi2, _) = span_iso(K_S.dual())
-        sage: J_T = phi2(K_S.dual()).dual()
-        sage: l = K.linear_subspace().dimension()
-        sage: codim = K.lattice_dim() - K.dim()
-        sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
-        sage: actual == expected
-        True
+    if m < n:
+        # K is not solid, project onto its span.
+        K = project_span(K)
 
-    And test with the project_span function::
+        # Lemma 2
+        beta += m*(n - m) + (n - m)**2
 
-        sage: K = random_cone(max_dim=15)
-        sage: actual = lyapunov_rank(K)
-        sage: K_S = project_span(K)
-        sage: P = project_span(K_S.dual()).dual()
-        sage: l = K.linear_subspace().dimension()
-        sage: codim = K.lattice_dim() - K.dim()
-        sage: expected = lyapunov_rank(P) + K.dim()*(l + codim) + codim**2
-        sage: actual == expected
-        True
+    if l > 0:
+        # K is not pointed, project its dual onto its span.
+        K = project_span(K.dual()).dual()
 
-    """
-    return len(LL(K))
+        # Lemma 3
+        beta += m * l
+
+    beta += len(LL(K))
+    return beta