mjo/cone/cone.py

   1 # Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
   2 # have to explicitly mangle our sitedir here so that "mjo.cone"
   3 # resolves.
   4 from os.path import abspath
   5 from site import addsitedir
   6 addsitedir(abspath('../../'))
   7
   8 from sage.all import *
   9
  10
  11 def random_cone(min_dim=None, max_dim=None, min_rays=None, max_rays=None):
  12     r"""
  13     Generate a random rational convex polyhedral cone.
  14
  15     Lower and upper bounds may be provided for both the dimension of the
  16     ambient space and the number of generating rays of the cone. Any
  17     parameters left unspecified will be chosen randomly.
  18
  19     INPUT:
  20
  21     - ``min_dim`` (default: random) -- The minimum dimension of the ambient
  22                                        lattice.
  23
  24     - ``max_dim`` (default: random) -- The maximum dimension of the ambient
  25                                        lattice.
  26
  27     - ``min_rays`` (default: random) -- The minimum number of generating rays
  28                                         of the cone.
  29
  30     - ``max_rays`` (default: random) -- The maximum number of generating rays
  31                                         of the cone.
  32
  33     OUTPUT:
  34
  35     A new, randomly generated cone.
  36
  37     TESTS:
  38
  39     It's hard to test the output of a random process, but we can at
  40     least make sure that we get a cone back::
  41
  42         sage: from sage.geometry.cone import is_Cone
  43         sage: K = random_cone()
  44         sage: is_Cone(K) # long time
  45         True
  46
  47     """
  48
  49     def random_min_max(l,u):
  50         r"""
  51         We need to handle four cases to prevent us from doing
  52         something stupid like having an upper bound that's lower than
  53         our lower bound. And we would need to repeat all of that logic
  54         for the dimension/rays, so we consolidate it here.
  55         """
  56         if l is None and u is None:
  57             # They're both random, just return a random nonnegative
  58             # integer.
  59             return ZZ.random_element().abs()
  60
  61         if l is not None and u is not None:
  62             # Both were specified. Again, just make up a number and
  63             # return it. If the user wants to give us u < l then he
  64             # can have an exception.
  65             return ZZ.random_element(l,u)
  66
  67         if l is not None and u is None:
  68             # In this case, we're generating the upper bound randomly
  69             # GIVEN A LOWER BOUND. So we add a random nonnegative
  70             # integer to the given lower bound.
  71             u = l + ZZ.random_element().abs()
  72             return ZZ.random_element(l,u)
  73
  74         # Here we must be in the only remaining case, where we are
  75         # given an upper bound but no lower bound. We might as well
  76         # use zero.
  77         return ZZ.random_element(0,u)
  78
  79     d = random_min_max(min_dim, max_dim)
  80     r = random_min_max(min_rays, max_rays)
  81
  82     L = ToricLattice(d)
  83     rays = [L.random_element() for i in range(0,r)]
  84
  85     # We pass the lattice in case there are no rays.
  86     return Cone(rays, lattice=L)
  87
  88
  89 def lyapunov_rank(K):
  90     r"""
  91     Compute the Lyapunov (or bilinearity) rank of this cone.
  92
  93     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
  94
  95     1. The dimension of the Lie algebra of the automorphism group of the
  96        cone.
  97
  98     2. The dimension of the linear space of all Lyapunov-like
  99        transformations on the cone.
 100
 101     INPUT:
 102
 103     A closed, convex polyhedral cone.
 104
 105     OUTPUT:
 106
 107     An integer representing the Lyapunov rank of the cone. If the
 108     dimension of the ambient vector space is `n`, then the Lyapunov rank
 109     will be between `1` and `n` inclusive; however a rank of `n-1` is
 110     not possible (see the first reference).
 111
 112     .. note::
 113
 114         In the references, the cones are always assumed to be proper. We
 115         do not impose this restriction.
 116
 117     .. seealso::
 118
 119         :meth:`is_proper`
 120
 121     ALGORITHM:
 122
 123     The codimension formula from the second reference is used. We find
 124     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 125     and `s` are rays of our cone. It is known that these vectors are
 126     sufficient to apply the codimension formula. Once we have all such
 127     pairs, we "brute force" the codimension formula by finding all
 128     linearly-independent `xs^{T}`.
 129
 130     REFERENCES:
 131
 132     1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone
 133        and Lyapunov-like transformations, Mathematical Programming, 147
 134        (2014) 155-170.
 135
 136     2. G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 137        optimality constraints for the cone of positive polynomials,
 138        Mathematical Programming, Series B, 129 (2011) 5-31.
 139
 140     EXAMPLES:
 141
 142     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`::
 143
 144         sage: positives = Cone([(1,)])
 145         sage: lyapunov_rank(positives)
 146         1
 147         sage: quadrant = Cone([(1,0), (0,1)])
 148         sage: lyapunov_rank(quadrant)
 149         2
 150         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 151         sage: lyapunov_rank(octant)
 152         3
 153
 154     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one::
 155
 156         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 157         sage: lyapunov_rank(L31)
 158         1
 159
 160     Likewise for the `L^{3}_{\infty}` cone::
 161
 162         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 163         sage: lyapunov_rank(L3infty)
 164         1
 165
 166     The Lyapunov rank should be additive on a product of cones::
 167
 168         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 169         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 170         sage: K = L31.cartesian_product(octant)
 171         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 172         True
 173
 174     Two isomorphic cones should have the same Lyapunov rank. The cone
 175     ``K`` in the following example is isomorphic to the nonnegative
 176     octant in `\mathbb{R}^{3}`::
 177
 178         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 179         sage: lyapunov_rank(K)
 180         3
 181
 182     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 183     itself::
 184
 185         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 186         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 187         True
 188
 189     TESTS:
 190
 191     The Lyapunov rank should be additive on a product of cones::
 192
 193         sage: K1 = random_cone(0,10,0,10)
 194         sage: K2 = random_cone(0,10,0,10)
 195         sage: K = K1.cartesian_product(K2)
 196         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 197         True
 198
 199     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 200     itself::
 201
 202         sage: K = random_cone(0,10,0,10)
 203         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 204         True
 205
 206     """
 207     V = K.lattice().vector_space()
 208
 209     xs = [V(x) for x in K.rays()]
 210     ss = [V(s) for s in K.dual().rays()]
 211
 212     # WARNING: This isn't really C(K), it only contains the pairs
 213     # (x,s) in C(K) where x,s are extreme in their respective cones.
 214     C_of_K = [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 215
 216     matrices = [x.column() * s.row() for (x,s) in C_of_K]
 217
 218     # Sage doesn't think matrices are vectors, so we have to convert
 219     # our matrices to vectors explicitly before we can figure out how
 220     # many are linearly-indepenedent.
 221     #
 222     # The space W has the same base ring as V, but dimension
 223     # dim(V)^2. So it has the same dimension as the space of linear
 224     # transformations on V. In other words, it's just the right size
 225     # to create an isomorphism between it and our matrices.
 226     W = VectorSpace(V.base_ring(), V.dimension()**2)
 227
 228     def phi(m):
 229         r"""
 230         Convert a matrix to a vector isomorphically.
 231         """
 232         return W(m.list())
 233
 234     vectors = [phi(m) for m in matrices]
 235
 236     return (W.dimension() - W.span(vectors).rank())