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 discrete_complementarity_set(K):
  90     r"""
  91     Compute the discrete complementarity set of this cone.
  92
  93     The complementarity set of this cone is the set of all orthogonal
  94     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
  95     dual. The discrete complementarity set restricts `x` and `s` to be
  96     generators of their respective cones.
  97
  98     OUTPUT:
  99
 100     A list of pairs `(x,s)` such that,
 101
 102       * `x` is in this cone.
 103       * `x` is a generator of this cone.
 104       * `s` is in this cone's dual.
 105       * `s` is a generator of this cone's dual.
 106       * `x` and `s` are orthogonal.
 107
 108     EXAMPLES:
 109
 110     The discrete complementarity set of the nonnegative orthant consists
 111     of pairs of standard basis vectors::
 112
 113         sage: K = Cone([(1,0),(0,1)])
 114         sage: discrete_complementarity_set(K)
 115         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
 116
 117     If the cone consists of a single ray, the second components of the
 118     discrete complementarity set should generate the orthogonal
 119     complement of that ray::
 120
 121         sage: K = Cone([(1,0)])
 122         sage: discrete_complementarity_set(K)
 123         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
 124         sage: K = Cone([(1,0,0)])
 125         sage: discrete_complementarity_set(K)
 126         [((1, 0, 0), (0, 1, 0)),
 127           ((1, 0, 0), (0, -1, 0)),
 128           ((1, 0, 0), (0, 0, 1)),
 129           ((1, 0, 0), (0, 0, -1))]
 130
 131     When the cone is the entire space, its dual is the trivial cone, so
 132     the discrete complementarity set is empty::
 133
 134         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 135         sage: discrete_complementarity_set(K)
 136         []
 137
 138     TESTS:
 139
 140     The complementarity set of the dual can be obtained by switching the
 141     components of the complementarity set of the original cone::
 142
 143         sage: K1 = random_cone(0,10,0,10)
 144         sage: K2 = K1.dual()
 145         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 146         sage: actual = discrete_complementarity_set(K1)
 147         sage: actual == expected
 148         True
 149
 150     """
 151     V = K.lattice().vector_space()
 152
 153     # Convert the rays to vectors so that we can compute inner
 154     # products.
 155     xs = [V(x) for x in K.rays()]
 156     ss = [V(s) for s in K.dual().rays()]
 157
 158     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 159
 160
 161 def lyapunov_rank(K):
 162     r"""
 163     Compute the Lyapunov (or bilinearity) rank of this cone.
 164
 165     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 166
 167     1. The dimension of the Lie algebra of the automorphism group of the
 168        cone.
 169
 170     2. The dimension of the linear space of all Lyapunov-like
 171        transformations on the cone.
 172
 173     INPUT:
 174
 175     A closed, convex polyhedral cone.
 176
 177     OUTPUT:
 178
 179     An integer representing the Lyapunov rank of the cone. If the
 180     dimension of the ambient vector space is `n`, then the Lyapunov rank
 181     will be between `1` and `n` inclusive; however a rank of `n-1` is
 182     not possible (see the first reference).
 183
 184     .. note::
 185
 186         In the references, the cones are always assumed to be proper. We
 187         do not impose this restriction.
 188
 189     .. seealso::
 190
 191         :meth:`is_proper`
 192
 193     ALGORITHM:
 194
 195     The codimension formula from the second reference is used. We find
 196     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 197     and `s` are rays of our cone. It is known that these vectors are
 198     sufficient to apply the codimension formula. Once we have all such
 199     pairs, we "brute force" the codimension formula by finding all
 200     linearly-independent `xs^{T}`.
 201
 202     REFERENCES:
 203
 204     1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone
 205        and Lyapunov-like transformations, Mathematical Programming, 147
 206        (2014) 155-170.
 207
 208     2. G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 209        optimality constraints for the cone of positive polynomials,
 210        Mathematical Programming, Series B, 129 (2011) 5-31.
 211
 212     EXAMPLES:
 213
 214     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`::
 215
 216         sage: positives = Cone([(1,)])
 217         sage: lyapunov_rank(positives)
 218         1
 219         sage: quadrant = Cone([(1,0), (0,1)])
 220         sage: lyapunov_rank(quadrant)
 221         2
 222         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 223         sage: lyapunov_rank(octant)
 224         3
 225
 226     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one::
 227
 228         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 229         sage: lyapunov_rank(L31)
 230         1
 231
 232     Likewise for the `L^{3}_{\infty}` cone::
 233
 234         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 235         sage: lyapunov_rank(L3infty)
 236         1
 237
 238     The Lyapunov rank should be additive on a product of cones::
 239
 240         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 241         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 242         sage: K = L31.cartesian_product(octant)
 243         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 244         True
 245
 246     Two isomorphic cones should have the same Lyapunov rank. The cone
 247     ``K`` in the following example is isomorphic to the nonnegative
 248     octant in `\mathbb{R}^{3}`::
 249
 250         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 251         sage: lyapunov_rank(K)
 252         3
 253
 254     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 255     itself::
 256
 257         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 258         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 259         True
 260
 261     TESTS:
 262
 263     The Lyapunov rank should be additive on a product of cones::
 264
 265         sage: K1 = random_cone(0,10,0,10)
 266         sage: K2 = random_cone(0,10,0,10)
 267         sage: K = K1.cartesian_product(K2)
 268         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 269         True
 270
 271     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 272     itself::
 273
 274         sage: K = random_cone(0,10,0,10)
 275         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 276         True
 277
 278     """
 279     V = K.lattice().vector_space()
 280
 281     C_of_K = discrete_complementarity_set(K)
 282
 283     matrices = [x.tensor_product(s) for (x,s) in C_of_K]
 284
 285     # Sage doesn't think matrices are vectors, so we have to convert
 286     # our matrices to vectors explicitly before we can figure out how
 287     # many are linearly-indepenedent.
 288     #
 289     # The space W has the same base ring as V, but dimension
 290     # dim(V)^2. So it has the same dimension as the space of linear
 291     # transformations on V. In other words, it's just the right size
 292     # to create an isomorphism between it and our matrices.
 293     W = VectorSpace(V.base_ring(), V.dimension()**2)
 294
 295     def phi(m):
 296         r"""
 297         Convert a matrix to a vector isomorphically.
 298         """
 299         return W(m.list())
 300
 301     vectors = [phi(m) for m in matrices]
 302
 303     return (W.dimension() - W.span(vectors).rank())