From 3b659b1d0440daf3ff7bd8cf3cf53f90523a1609 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Thu, 14 May 2015 16:18:17 -0400
Subject: [PATCH] Add a random_cone() function and use it in two tests.

---
 mjo/cone/cone.py | 95 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 95 insertions(+)

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index 2304338..2296e3f 100644
--- a/mjo/cone/cone.py
+++ b/mjo/cone/cone.py
@@ -8,6 +8,84 @@ addsitedir(abspath('../../'))
 from sage.all import *
 
 
+def random_cone(min_dim=None, max_dim=None, min_rays=None, max_rays=None):
+    r"""
+    Generate a random rational convex polyhedral cone.
+
+    Lower and upper bounds may be provided for both the dimension of the
+    ambient space and the number of generating rays of the cone. Any
+    parameters left unspecified will be chosen randomly.
+
+    INPUT:
+
+    - ``min_dim`` (default: random) -- The minimum dimension of the ambient
+                                       lattice.
+
+    - ``max_dim`` (default: random) -- The maximum dimension of the ambient
+                                       lattice.
+
+    - ``min_rays`` (default: random) -- The minimum number of generating rays
+                                        of the cone.
+
+    - ``max_rays`` (default: random) -- The maximum number of generating rays
+                                        of the cone.
+
+    OUTPUT:
+
+    A new, randomly generated cone.
+
+    TESTS:
+
+    It's hard to test the output of a random process, but we can at
+    least make sure that we get a cone back::
+
+        sage: from sage.geometry.cone import is_Cone
+        sage: K = random_cone()
+        sage: is_Cone(K) # long time
+        True
+
+    """
+
+    def random_min_max(l,u):
+        r"""
+        We need to handle four cases to prevent us from doing
+        something stupid like having an upper bound that's lower than
+        our lower bound. And we would need to repeat all of that logic
+        for the dimension/rays, so we consolidate it here.
+        """
+        if l is None and u is None:
+            # They're both random, just return a random nonnegative
+            # integer.
+            return ZZ.random_element().abs()
+
+        if l is not None and u is not None:
+            # Both were specified. Again, just make up a number and
+            # return it. If the user wants to give us u < l then he
+            # can have an exception.
+            return ZZ.random_element(l,u)
+
+        if l is not None and u is None:
+            # In this case, we're generating the upper bound randomly
+            # GIVEN A LOWER BOUND. So we add a random nonnegative
+            # integer to the given lower bound.
+            u = l + ZZ.random_element().abs()
+            return ZZ.random_element(l,u)
+
+        # Here we must be in the only remaining case, where we are
+        # given an upper bound but no lower bound. We might as well
+        # use zero.
+        return ZZ.random_element(0,u)
+
+    d = random_min_max(min_dim, max_dim)
+    r = random_min_max(min_rays, max_rays)
+
+    L = ToricLattice(d)
+    rays = [L.random_element() for i in range(0,r)]
+
+    # We pass the lattice in case there are no rays.
+    return Cone(rays, lattice=L)
+
+
 def lyapunov_rank(K):
     r"""
     Compute the Lyapunov (or bilinearity) rank of this cone.
@@ -108,6 +186,23 @@ def lyapunov_rank(K):
         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
         True
 
+    TESTS:
+
+    The Lyapunov rank should be additive on a product of cones::
+
+        sage: K1 = random_cone(0,10,0,10)
+        sage: K2 = random_cone(0,10,0,10)
+        sage: K = K1.cartesian_product(K2)
+        sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
+        True
+
+    The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
+    itself::
+
+        sage: K = random_cone(0,10,0,10)
+        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
+        True
+
     """
     V = K.lattice().vector_space()
 
-- 
2.44.2