]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/cone/cone.py
Ensure that we generate min_rays generators; add more cone tests.
[sage.d.git] / mjo / cone / cone.py
index 507b6cea5c24d0a4eb745d9245c85a1fbf4ddb90..132c6d9e8e63ebbc5e0065becd813f341e993dd0 100644 (file)
@@ -7,83 +7,198 @@ addsitedir(abspath('../../'))
 
 from sage.all import *
 
-
-def random_cone(min_dim=None, max_dim=None, min_rays=None, max_rays=None):
+# TODO: This test fails, maybe due to a bug in the existing cone code.
+#     If we request enough generators to span the space, then the returned
+#     cone should equal the ambient space::
+#
+#        sage: K = random_cone(min_dim=5, max_dim=5, min_rays=10, max_rays=10)
+#        sage: K.lines().dimension() == K.lattice_dim()
+#        True
+
+def random_cone(min_dim=0, max_dim=None, min_rays=0, 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.
+    ambient space and the number of generating rays of the cone. If a
+    lower bound is left unspecified, it defaults to zero. Unspecified
+    upper bounds will be chosen randomly.
+
+    The number of generating rays is naturally limited to twice the
+    dimension of the ambient space. Take for example $\mathbb{R}^{2}$.
+    You could have the generators $\left\{ \pm e_{1}, \pm e_{2}
+    \right\}$, with cardinality $4 = 2 \cdot 2$; however any other ray
+    in the space is a nonnegative linear combination of those four.
+
+    .. NOTE:
+
+        If you do not explicitly request more than ``2 * max_dim`` rays,
+        a larger number may still be randomly generated. In that case,
+        the returned cone will simply be equal to the entire space.
 
     INPUT:
 
-    - ``min_dim`` (default: random) -- The minimum dimension of the ambient
-                                       lattice.
+    - ``min_dim`` (default: zero) -- A nonnegative integer representing the
+                                     minimum dimension of the ambient lattice.
 
-    - ``max_dim`` (default: random) -- The maximum dimension of the ambient
+    - ``max_dim`` (default: random) -- A nonnegative integer representing
+                                       the maximum dimension of the ambient
                                        lattice.
 
-    - ``min_rays`` (default: random) -- The minimum number of generating rays
-                                        of the cone.
+    - ``min_rays`` (default: zero) -- A nonnegative integer representing the
+                                      minimum number of generating rays of the
+                                      cone.
 
-    - ``max_rays`` (default: random) -- The maximum number of generating rays
-                                        of the cone.
+    - ``max_rays`` (default: random) -- A nonnegative integer representing the
+                                        maximum number of generating rays of
+                                        the cone.
 
     OUTPUT:
 
     A new, randomly generated cone.
 
+    A ``ValueError` will be thrown under the following conditions:
+
+      * Any of ``min_dim``, ``max_dim``, ``min_rays``, or ``max_rays``
+        are negative.
+
+      * ``max_dim`` is less than ``min_dim``.
+
+      * ``max_rays`` is less than ``min_rays``.
+
+      * ``min_rays`` is greater than twice ``max_dim``.
+
+    EXAMPLES:
+
+    If we set the lower/upper bounds to zero, then our result is
+    predictable::
+
+        sage: random_cone(0,0,0,0)
+        0-d cone in 0-d lattice N
+
+    We can predict the dimension when ``min_dim == max_dim``::
+
+        sage: random_cone(min_dim=4, max_dim=4, min_rays=0, max_rays=0)
+        0-d cone in 4-d lattice N
+
+    Likewise for the number of rays when ``min_rays == max_rays``::
+
+        sage: random_cone(min_dim=10, max_dim=10, min_rays=10, max_rays=10)
+        10-d cone in 10-d lattice N
+
     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
+        sage: from sage.geometry.cone import is_Cone # long time
+        sage: K = random_cone() # long time
+        sage: is_Cone(K)        # long time
         True
 
+    The upper/lower bounds are respected::
+
+        sage: K = random_cone(min_dim=5, max_dim=10, min_rays=3, max_rays=4)
+        sage: 5 <= K.lattice_dim() and K.lattice_dim() <= 10
+        True
+        sage: 3 <= K.nrays() and K.nrays() <= 4
+        True
+
+    Ensure that an exception is raised when either lower bound is greater
+    than its respective upper bound::
+
+        sage: random_cone(min_dim=5, max_dim=2)
+        Traceback (most recent call last):
+        ...
+        ValueError: max_dim cannot be less than min_dim.
+
+        sage: random_cone(min_rays=5, max_rays=2)
+        Traceback (most recent call last):
+        ...
+        ValueError: max_rays cannot be less than min_rays.
+
+    And if we request too many rays::
+
+        sage: random_cone(min_rays=5, max_dim=1)
+        Traceback (most recent call last):
+        ...
+        ValueError: min_rays cannot be larger than twice max_dim.
+
     """
 
+    # Catch obvious mistakes so that we can generate clear error
+    # messages.
+
+    if min_dim < 0:
+        raise ValueError('min_dim must be nonnegative.')
+
+    if min_rays < 0:
+        raise ValueError('min_rays must be nonnegative.')
+
+    if max_dim is not None:
+        if max_dim < 0:
+            raise ValueError('max_dim must be nonnegative.')
+        if (max_dim < min_dim):
+            raise ValueError('max_dim cannot be less than min_dim.')
+        if min_rays > 2*max_dim:
+            raise ValueError('min_rays cannot be larger than twice max_dim.')
+
+    if max_rays is not None:
+        if max_rays < 0:
+            raise ValueError('max_rays must be nonnegative.')
+        if (max_rays < min_rays):
+            raise ValueError('max_rays cannot be less than min_rays.')
+
+
     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.
+        We need to handle two cases for the upper bounds, and we need to do
+        the same thing for max_dim/max_rays. So we consolidate the logic 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)
+        if u is None:
+            # The upper bound is unspecified; return a random integer
+            # in [l,infinity).
+            return l + ZZ.random_element().abs()
+        else:
+            # We have an upper bound, and it's greater than or equal
+            # to our lower bound. So we generate a random integer in
+            # [0,u-l], and then add it to l to get something in
+            # [l,u]. To understand the "+1", check the
+            # ZZ.random_element() docs.
+            return l + ZZ.random_element(u - l + 1)
+
 
     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)
+    # The rays are trickier to generate, since we could generate v and
+    # 2*v as our "two rays." In that case, the resuting cone would
+    # have one generating ray. To avoid such a situation, we start by
+    # generating ``r`` rays where ``r`` is the number we want to end
+    # up with.
+    #
+    # However, since we're going to *check* whether or not we actually
+    # have ``r``, we need ``r`` rays to be attainable. So we need to
+    # limit ``r`` to twice the dimension of the ambient space.
+    #
+    r = min(r, 2*d)
+    rays = [L.random_element() for i in range(0, r)]
+
+    # (The lattice parameter is required when no rays are given, so we
+    # pass it just in case ``r == 0``).
+    K = Cone(rays, lattice=L)
+
+    # Now if we generated two of the "same" rays, we'll have fewer
+    # generating rays than ``r``. In that case, we keep making up new
+    # rays and recreating the cone until we get the right number of
+    # independent generators.
+    while r > K.nrays():
+        rays.append(L.random_element())
+        K = Cone(rays)
+
+    return K
 
 
 def discrete_complementarity_set(K):
@@ -140,7 +255,7 @@ def discrete_complementarity_set(K):
     The complementarity set of the dual can be obtained by switching the
     components of the complementarity set of the original cone::
 
-        sage: K1 = random_cone(0,10,0,10)
+        sage: K1 = random_cone(max_dim=10, max_rays=10)
         sage: K2 = K1.dual()
         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
         sage: actual = discrete_complementarity_set(K1)
@@ -262,8 +377,8 @@ def lyapunov_rank(K):
 
     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: K1 = random_cone(max_dim=10, max_rays=10)
+        sage: K2 = random_cone(max_dim=10, max_rays=10)
         sage: K = K1.cartesian_product(K2)
         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
         True
@@ -271,7 +386,7 @@ def lyapunov_rank(K):
     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
     itself::
 
-        sage: K = random_cone(0,10,0,10)
+        sage: K = random_cone(max_dim=10, max_rays=10)
         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
         True