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=0, max_dim=None, min_rays=0, 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. If a
  17     lower bound is left unspecified, it defaults to zero. Unspecified
  18     upper bounds will be chosen randomly.
  19
  20     INPUT:
  21
  22     - ``min_dim`` (default: zero) -- A nonnegative integer representing the
  23                                      minimum dimension of the ambient lattice.
  24
  25     - ``max_dim`` (default: random) -- A nonnegative integer representing
  26                                        the maximum dimension of the ambient
  27                                        lattice.
  28
  29     - ``min_rays`` (default: zero) -- A nonnegative integer representing the
  30                                       minimum number of generating rays of the
  31                                       cone.
  32
  33     - ``max_rays`` (default: random) -- A nonnegative integer representing the
  34                                         maximum number of generating rays of the
  35                                         cone.
  36
  37     OUTPUT:
  38
  39     A new, randomly generated cone.
  40
  41     EXAMPLES:
  42
  43     If we set the lower/upper bounds to zero, then our result is
  44     predictable::
  45
  46         sage: random_cone(0,0,0,0)
  47         0-d cone in 0-d lattice N
  48
  49     In fact, as long as we ask for zero rays, we should be able to predict
  50     the output when ``min_dim == max_dim``::
  51
  52         sage: random_cone(min_dim=4, max_dim=4, min_rays=0, max_rays=0)
  53         0-d cone in 4-d lattice N
  54
  55     TESTS:
  56
  57     It's hard to test the output of a random process, but we can at
  58     least make sure that we get a cone back::
  59
  60         sage: from sage.geometry.cone import is_Cone # long time
  61         sage: K = random_cone() # long time
  62         sage: is_Cone(K)        # long time
  63         True
  64
  65     Ensure that an exception is raised when either lower bound is greater
  66     than its respective upper bound::
  67
  68         sage: random_cone(min_dim=5, max_dim=2)
  69         Traceback (most recent call last):
  70         ...
  71         ValueError: max_dim must be greater than or equal to min_dim.
  72
  73         sage: random_cone(min_rays=5, max_rays=2)
  74         Traceback (most recent call last):
  75         ...
  76         ValueError: max_rays must be greater than or equal to min_rays.
  77
  78     """
  79
  80     # Catch obvious mistakes so that we can generate clear error
  81     # messages.
  82
  83     if min_dim < 0:
  84         raise ValueError('min_dim must be nonnegative.')
  85
  86     if min_rays < 0:
  87         raise ValueError('min_rays must be nonnegative.')
  88
  89     if max_dim is not None:
  90         if max_dim < 0:
  91             raise ValueError('max_dim must be nonnegative.')
  92         if (min_dim > max_dim):
  93             raise ValueError('max_dim must be greater than or equal to min_dim.')
  94
  95     if max_rays is not None:
  96         if max_rays < 0:
  97             raise ValueError('max_rays must be nonnegative.')
  98         if (min_rays > max_rays):
  99             raise ValueError('max_rays must be greater than or equal to min_rays.')
 100
 101
 102     def random_min_max(l,u):
 103         r"""
 104         We need to handle two cases for the upper bounds, and we need to do
 105         the same thing for max_dim/max_rays. So we consolidate the logic here.
 106         """
 107         if u is None:
 108             # The upper bound is unspecified; return a random integer
 109             # in [l,infinity).
 110             return l + ZZ.random_element().abs()
 111         else:
 112             # We have an upper bound, and it's greater than or equal
 113             # to our lower bound. So we generate a random integer in
 114             # [0,u-l], and then add it to l to get something in
 115             # [l,u]. To understand the "+1", check the
 116             # ZZ.random_element() docs.
 117             return l + ZZ.random_element(u - l + 1)
 118
 119
 120     d = random_min_max(min_dim, max_dim)
 121     r = random_min_max(min_rays, max_rays)
 122
 123     L = ToricLattice(d)
 124     rays = [L.random_element() for i in range(0,r)]
 125
 126     # The lattice parameter is required when no rays are given, so we
 127     # pass it just in case.
 128     return Cone(rays, lattice=L)
 129
 130
 131 def discrete_complementarity_set(K):
 132     r"""
 133     Compute the discrete complementarity set of this cone.
 134
 135     The complementarity set of this cone is the set of all orthogonal
 136     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
 137     dual. The discrete complementarity set restricts `x` and `s` to be
 138     generators of their respective cones.
 139
 140     OUTPUT:
 141
 142     A list of pairs `(x,s)` such that,
 143
 144       * `x` is in this cone.
 145       * `x` is a generator of this cone.
 146       * `s` is in this cone's dual.
 147       * `s` is a generator of this cone's dual.
 148       * `x` and `s` are orthogonal.
 149
 150     EXAMPLES:
 151
 152     The discrete complementarity set of the nonnegative orthant consists
 153     of pairs of standard basis vectors::
 154
 155         sage: K = Cone([(1,0),(0,1)])
 156         sage: discrete_complementarity_set(K)
 157         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
 158
 159     If the cone consists of a single ray, the second components of the
 160     discrete complementarity set should generate the orthogonal
 161     complement of that ray::
 162
 163         sage: K = Cone([(1,0)])
 164         sage: discrete_complementarity_set(K)
 165         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
 166         sage: K = Cone([(1,0,0)])
 167         sage: discrete_complementarity_set(K)
 168         [((1, 0, 0), (0, 1, 0)),
 169           ((1, 0, 0), (0, -1, 0)),
 170           ((1, 0, 0), (0, 0, 1)),
 171           ((1, 0, 0), (0, 0, -1))]
 172
 173     When the cone is the entire space, its dual is the trivial cone, so
 174     the discrete complementarity set is empty::
 175
 176         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 177         sage: discrete_complementarity_set(K)
 178         []
 179
 180     TESTS:
 181
 182     The complementarity set of the dual can be obtained by switching the
 183     components of the complementarity set of the original cone::
 184
 185         sage: K1 = random_cone(max_dim=10, max_rays=10)
 186         sage: K2 = K1.dual()
 187         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 188         sage: actual = discrete_complementarity_set(K1)
 189         sage: actual == expected
 190         True
 191
 192     """
 193     V = K.lattice().vector_space()
 194
 195     # Convert the rays to vectors so that we can compute inner
 196     # products.
 197     xs = [V(x) for x in K.rays()]
 198     ss = [V(s) for s in K.dual().rays()]
 199
 200     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 201
 202
 203 def lyapunov_rank(K):
 204     r"""
 205     Compute the Lyapunov (or bilinearity) rank of this cone.
 206
 207     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 208
 209     1. The dimension of the Lie algebra of the automorphism group of the
 210        cone.
 211
 212     2. The dimension of the linear space of all Lyapunov-like
 213        transformations on the cone.
 214
 215     INPUT:
 216
 217     A closed, convex polyhedral cone.
 218
 219     OUTPUT:
 220
 221     An integer representing the Lyapunov rank of the cone. If the
 222     dimension of the ambient vector space is `n`, then the Lyapunov rank
 223     will be between `1` and `n` inclusive; however a rank of `n-1` is
 224     not possible (see the first reference).
 225
 226     .. note::
 227
 228         In the references, the cones are always assumed to be proper. We
 229         do not impose this restriction.
 230
 231     .. seealso::
 232
 233         :meth:`is_proper`
 234
 235     ALGORITHM:
 236
 237     The codimension formula from the second reference is used. We find
 238     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 239     and `s` are rays of our cone. It is known that these vectors are
 240     sufficient to apply the codimension formula. Once we have all such
 241     pairs, we "brute force" the codimension formula by finding all
 242     linearly-independent `xs^{T}`.
 243
 244     REFERENCES:
 245
 246     1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone
 247        and Lyapunov-like transformations, Mathematical Programming, 147
 248        (2014) 155-170.
 249
 250     2. G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 251        optimality constraints for the cone of positive polynomials,
 252        Mathematical Programming, Series B, 129 (2011) 5-31.
 253
 254     EXAMPLES:
 255
 256     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`::
 257
 258         sage: positives = Cone([(1,)])
 259         sage: lyapunov_rank(positives)
 260         1
 261         sage: quadrant = Cone([(1,0), (0,1)])
 262         sage: lyapunov_rank(quadrant)
 263         2
 264         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 265         sage: lyapunov_rank(octant)
 266         3
 267
 268     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one::
 269
 270         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 271         sage: lyapunov_rank(L31)
 272         1
 273
 274     Likewise for the `L^{3}_{\infty}` cone::
 275
 276         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 277         sage: lyapunov_rank(L3infty)
 278         1
 279
 280     The Lyapunov rank should be additive on a product of cones::
 281
 282         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 283         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 284         sage: K = L31.cartesian_product(octant)
 285         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 286         True
 287
 288     Two isomorphic cones should have the same Lyapunov rank. The cone
 289     ``K`` in the following example is isomorphic to the nonnegative
 290     octant in `\mathbb{R}^{3}`::
 291
 292         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 293         sage: lyapunov_rank(K)
 294         3
 295
 296     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 297     itself::
 298
 299         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 300         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 301         True
 302
 303     TESTS:
 304
 305     The Lyapunov rank should be additive on a product of cones::
 306
 307         sage: K1 = random_cone(max_dim=10, max_rays=10)
 308         sage: K2 = random_cone(max_dim=10, max_rays=10)
 309         sage: K = K1.cartesian_product(K2)
 310         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 311         True
 312
 313     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 314     itself::
 315
 316         sage: K = random_cone(max_dim=10, max_rays=10)
 317         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 318         True
 319
 320     """
 321     V = K.lattice().vector_space()
 322
 323     C_of_K = discrete_complementarity_set(K)
 324
 325     matrices = [x.tensor_product(s) for (x,s) in C_of_K]
 326
 327     # Sage doesn't think matrices are vectors, so we have to convert
 328     # our matrices to vectors explicitly before we can figure out how
 329     # many are linearly-indepenedent.
 330     #
 331     # The space W has the same base ring as V, but dimension
 332     # dim(V)^2. So it has the same dimension as the space of linear
 333     # transformations on V. In other words, it's just the right size
 334     # to create an isomorphism between it and our matrices.
 335     W = VectorSpace(V.base_ring(), V.dimension()**2)
 336
 337     def phi(m):
 338         r"""
 339         Convert a matrix to a vector isomorphically.
 340         """
 341         return W(m.list())
 342
 343     vectors = [phi(m) for m in matrices]
 344
 345     return (W.dimension() - W.span(vectors).rank())