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 _basically_the_same(K1, K2):
  12     r"""
  13     Test whether or not ``K1`` and ``K2`` are "basically the same."
  14
  15     This is a hack to get around the fact that it's difficult to tell
  16     when two cones are linearly isomorphic. We have a proposition that
  17     equates two cones, but represented over `\mathbb{Q}`, they are
  18     merely linearly isomorphic (not equal). So rather than test for
  19     equality, we test a list of properties that should be preserved
  20     under an invertible linear transformation.
  21
  22     OUTPUT:
  23
  24     ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
  25     otherwise.
  26
  27     EXAMPLES:
  28
  29     Any proper cone with three generators in `\mathbb{R}^{3}` is
  30     basically the same as the nonnegative orthant::
  31
  32         sage: K1 = Cone([(1,0,0), (0,1,0), (0,0,1)])
  33         sage: K2 = Cone([(1,2,3), (3, 18, 4), (66, 51, 0)])
  34         sage: _basically_the_same(K1, K2)
  35         True
  36
  37     Negating a cone gives you another cone that is basically the same::
  38
  39         sage: K = Cone([(0,2,-5), (-6, 2, 4), (0, 51, 0)])
  40         sage: _basically_the_same(K, -K)
  41         True
  42
  43     TESTS:
  44
  45     Any cone is basically the same as itself::
  46
  47         sage: K = random_cone(max_ambient_dim = 8)
  48         sage: _basically_the_same(K, K)
  49         True
  50
  51     After applying an invertible matrix to the rows of a cone, the
  52     result should be basically the same as the cone we started with::
  53
  54         sage: K1 = random_cone(max_ambient_dim = 8)
  55         sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
  56         sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
  57         sage: _basically_the_same(K1, K2)
  58         True
  59
  60     """
  61     if K1.lattice_dim() != K2.lattice_dim():
  62         return False
  63
  64     if K1.nrays() != K2.nrays():
  65         return False
  66
  67     if K1.dim() != K2.dim():
  68         return False
  69
  70     if K1.lineality() != K2.lineality():
  71         return False
  72
  73     if K1.is_solid() != K2.is_solid():
  74         return False
  75
  76     if K1.is_strictly_convex() != K2.is_strictly_convex():
  77         return False
  78
  79     if len(LL(K1)) != len(LL(K2)):
  80         return False
  81
  82     C_of_K1 = K1.discrete_complementarity_set()
  83     C_of_K2 = K2.discrete_complementarity_set()
  84     if len(C_of_K1) != len(C_of_K2):
  85         return False
  86
  87     if len(K1.facets()) != len(K2.facets()):
  88         return False
  89
  90     return True
  91
  92
  93
  94 def _restrict_to_space(K, W):
  95     r"""
  96     Restrict this cone a subspace of its ambient space.
  97
  98     INPUT:
  99
 100     - ``W`` -- The subspace into which this cone will be restricted.
 101
 102     OUTPUT:
 103
 104     A new cone in a sublattice corresponding to ``W``.
 105
 106     EXAMPLES:
 107
 108     When this cone is solid, restricting it into its own span should do
 109     nothing::
 110
 111         sage: K = Cone([(1,)])
 112         sage: _restrict_to_space(K, K.span()) == K
 113         True
 114
 115     A single ray restricted into its own span gives the same output
 116     regardless of the ambient space::
 117
 118         sage: K2 = Cone([(1,0)])
 119         sage: K2_S = _restrict_to_space(K2, K2.span()).rays()
 120         sage: K2_S
 121         N(1)
 122         in 1-d lattice N
 123         sage: K3 = Cone([(1,0,0)])
 124         sage: K3_S = _restrict_to_space(K3, K3.span()).rays()
 125         sage: K3_S
 126         N(1)
 127         in 1-d lattice N
 128         sage: K2_S == K3_S
 129         True
 130
 131     TESTS:
 132
 133     The projected cone should always be solid::
 134
 135         sage: set_random_seed()
 136         sage: K = random_cone(max_ambient_dim = 8)
 137         sage: _restrict_to_space(K, K.span()).is_solid()
 138         True
 139
 140     And the resulting cone should live in a space having the same
 141     dimension as the space we restricted it to::
 142
 143         sage: set_random_seed()
 144         sage: K = random_cone(max_ambient_dim = 8)
 145         sage: K_P = _restrict_to_space(K, K.dual().span())
 146         sage: K_P.lattice_dim() == K.dual().dim()
 147         True
 148
 149     This function should not affect the dimension of a cone::
 150
 151         sage: set_random_seed()
 152         sage: K = random_cone(max_ambient_dim = 8)
 153         sage: K.dim() == _restrict_to_space(K,K.span()).dim()
 154         True
 155
 156     Nor should it affect the lineality of a cone::
 157
 158         sage: set_random_seed()
 159         sage: K = random_cone(max_ambient_dim = 8)
 160         sage: K.lineality() == _restrict_to_space(K, K.span()).lineality()
 161         True
 162
 163     No matter which space we restrict to, the lineality should not
 164     increase::
 165
 166         sage: set_random_seed()
 167         sage: K = random_cone(max_ambient_dim = 8)
 168         sage: S = K.span(); P = K.dual().span()
 169         sage: K.lineality() >= _restrict_to_space(K,S).lineality()
 170         True
 171         sage: K.lineality() >= _restrict_to_space(K,P).lineality()
 172         True
 173
 174     If we do this according to our paper, then the result is proper::
 175
 176         sage: set_random_seed()
 177         sage: K = random_cone(max_ambient_dim = 8)
 178         sage: K_S = _restrict_to_space(K, K.span())
 179         sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
 180         sage: K_SP.is_proper()
 181         True
 182         sage: K_SP = _restrict_to_space(K_S, K_S.dual().span())
 183         sage: K_SP.is_proper()
 184         True
 185
 186     Test the proposition in our paper concerning the duals and
 187     restrictions. Generate a random cone, then create a subcone of
 188     it. The operation of dual-taking should then commute with
 189     _restrict_to_space::
 190
 191         sage: set_random_seed()
 192         sage: J = random_cone(max_ambient_dim = 8)
 193         sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
 194         sage: K_W_star = _restrict_to_space(K, J.span()).dual()
 195         sage: K_star_W = _restrict_to_space(K.dual(), J.span())
 196         sage: _basically_the_same(K_W_star, K_star_W)
 197         True
 198
 199     """
 200     # First we want to intersect ``K`` with ``W``. The easiest way to
 201     # do this is via cone intersection, so we turn the subspace ``W``
 202     # into a cone.
 203     W_cone = Cone(W.basis() + [-b for b in W.basis()], lattice=K.lattice())
 204     K = K.intersection(W_cone)
 205
 206     # We've already intersected K with the span of K2, so every
 207     # generator of K should belong to W now.
 208     K_W_rays = [ W.coordinate_vector(r) for r in K.rays() ]
 209
 210     L = ToricLattice(W.dimension())
 211     return Cone(K_W_rays, lattice=L)
 212
 213
 214 def LL(K):
 215     r"""
 216     Compute a basis of Lyapunov-like transformations on this cone.
 217
 218     OUTPUT:
 219
 220     A list of matrices forming a basis for the space of all
 221     Lyapunov-like transformations on the given cone.
 222
 223     EXAMPLES:
 224
 225     The trivial cone has no Lyapunov-like transformations::
 226
 227         sage: L = ToricLattice(0)
 228         sage: K = Cone([], lattice=L)
 229         sage: LL(K)
 230         []
 231
 232     The Lyapunov-like transformations on the nonnegative orthant are
 233     simply diagonal matrices::
 234
 235         sage: K = Cone([(1,)])
 236         sage: LL(K)
 237         [[1]]
 238
 239         sage: K = Cone([(1,0),(0,1)])
 240         sage: LL(K)
 241         [
 242         [1 0]  [0 0]
 243         [0 0], [0 1]
 244         ]
 245
 246         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 247         sage: LL(K)
 248         [
 249         [1 0 0]  [0 0 0]  [0 0 0]
 250         [0 0 0]  [0 1 0]  [0 0 0]
 251         [0 0 0], [0 0 0], [0 0 1]
 252         ]
 253
 254     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 255     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 256
 257         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 258         sage: LL(L31)
 259         [
 260         [1 0 0]
 261         [0 1 0]
 262         [0 0 1]
 263         ]
 264
 265         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 266         sage: LL(L3infty)
 267         [
 268         [1 0 0]
 269         [0 1 0]
 270         [0 0 1]
 271         ]
 272
 273     If our cone is the entire space, then every transformation on it is
 274     Lyapunov-like::
 275
 276         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 277         sage: M = MatrixSpace(QQ,2)
 278         sage: M.basis() == LL(K)
 279         True
 280
 281     TESTS:
 282
 283     The inner product `\left< L\left(x\right), s \right>` is zero for
 284     every pair `\left( x,s \right)` in the discrete complementarity set
 285     of the cone::
 286
 287         sage: set_random_seed()
 288         sage: K = random_cone(max_ambient_dim=8)
 289         sage: C_of_K = K.discrete_complementarity_set()
 290         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 291         sage: sum(map(abs, l))
 292         0
 293
 294     The Lyapunov-like transformations on a cone and its dual are related
 295     by transposition, but we're not guaranteed to compute transposed
 296     elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
 297     \right)`
 298
 299         sage: set_random_seed()
 300         sage: K = random_cone(max_ambient_dim=8)
 301         sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
 302         sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
 303         sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
 304         sage: LL2_vecs = [ V(m.list()) for m in LL2   ]
 305         sage: V.span(LL1_vecs) == V.span(LL2_vecs)
 306         True
 307
 308     """
 309     V = K.lattice().vector_space()
 310
 311     C_of_K = K.discrete_complementarity_set()
 312
 313     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 314
 315     # Sage doesn't think matrices are vectors, so we have to convert
 316     # our matrices to vectors explicitly before we can figure out how
 317     # many are linearly-indepenedent.
 318     #
 319     # The space W has the same base ring as V, but dimension
 320     # dim(V)^2. So it has the same dimension as the space of linear
 321     # transformations on V. In other words, it's just the right size
 322     # to create an isomorphism between it and our matrices.
 323     W = VectorSpace(V.base_ring(), V.dimension()**2)
 324
 325     # Turn our matrices into long vectors...
 326     vectors = [ W(m.list()) for m in tensor_products ]
 327
 328     # Vector space representation of Lyapunov-like matrices
 329     # (i.e. vec(L) where L is Luapunov-like).
 330     LL_vector = W.span(vectors).complement()
 331
 332     # Now construct an ambient MatrixSpace in which to stick our
 333     # transformations.
 334     M = MatrixSpace(V.base_ring(), V.dimension())
 335
 336     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 337
 338     return matrix_basis
 339
 340
 341
 342 def lyapunov_rank(K):
 343     r"""
 344     Compute the Lyapunov rank (or bilinearity rank) of this cone.
 345
 346     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 347
 348     1. The dimension of the Lie algebra of the automorphism group of the
 349        cone.
 350
 351     2. The dimension of the linear space of all Lyapunov-like
 352        transformations on the cone.
 353
 354     INPUT:
 355
 356     A closed, convex polyhedral cone.
 357
 358     OUTPUT:
 359
 360     An integer representing the Lyapunov rank of the cone. If the
 361     dimension of the ambient vector space is `n`, then the Lyapunov rank
 362     will be between `1` and `n` inclusive; however a rank of `n-1` is
 363     not possible (see [Orlitzky/Gowda]_).
 364
 365     ALGORITHM:
 366
 367     The codimension formula from the second reference is used. We find
 368     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 369     and `s` are rays of our cone. It is known that these vectors are
 370     sufficient to apply the codimension formula. Once we have all such
 371     pairs, we "brute force" the codimension formula by finding all
 372     linearly-independent `xs^{T}`.
 373
 374     REFERENCES:
 375
 376     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 377        cone and Lyapunov-like transformations, Mathematical Programming, 147
 378        (2014) 155-170.
 379
 380     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 381        Improper Cone. Work in-progress.
 382
 383     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 384        optimality constraints for the cone of positive polynomials,
 385        Mathematical Programming, Series B, 129 (2011) 5-31.
 386
 387     EXAMPLES:
 388
 389     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 390     [Rudolf et al.]_::
 391
 392         sage: positives = Cone([(1,)])
 393         sage: lyapunov_rank(positives)
 394         1
 395         sage: quadrant = Cone([(1,0), (0,1)])
 396         sage: lyapunov_rank(quadrant)
 397         2
 398        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 399         sage: lyapunov_rank(octant)
 400         3
 401
 402     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 403     [Orlitzky/Gowda]_::
 404
 405         sage: R5 = VectorSpace(QQ, 5)
 406         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 407         sage: K = Cone(gs)
 408         sage: lyapunov_rank(K)
 409         25
 410
 411     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 412     [Rudolf et al.]_::
 413
 414         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 415         sage: lyapunov_rank(L31)
 416         1
 417
 418     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 419
 420         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 421         sage: lyapunov_rank(L3infty)
 422         1
 423
 424     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 425     + 1` [Orlitzky/Gowda]_::
 426
 427         sage: K = Cone([(1,0,0,0,0)])
 428         sage: lyapunov_rank(K)
 429         21
 430         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 431         21
 432
 433     A subspace (of dimension `m`) in `n` dimensions should have a
 434     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 435
 436         sage: e1 = (1,0,0,0,0)
 437         sage: neg_e1 = (-1,0,0,0,0)
 438         sage: e2 = (0,1,0,0,0)
 439         sage: neg_e2 = (0,-1,0,0,0)
 440         sage: z = (0,0,0,0,0)
 441         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 442         sage: lyapunov_rank(K)
 443         19
 444         sage: K.lattice_dim()**2 - K.dim()*K.codim()
 445         19
 446
 447     The Lyapunov rank should be additive on a product of proper cones
 448     [Rudolf et al.]_::
 449
 450         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 451         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 452         sage: K = L31.cartesian_product(octant)
 453         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 454         True
 455
 456     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 457     The cone ``K`` in the following example is isomorphic to the nonnegative
 458     octant in `\mathbb{R}^{3}`::
 459
 460         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 461         sage: lyapunov_rank(K)
 462         3
 463
 464     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 465     itself [Rudolf et al.]_::
 466
 467         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 468         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 469         True
 470
 471     TESTS:
 472
 473     The Lyapunov rank should be additive on a product of proper cones
 474     [Rudolf et al.]_::
 475
 476         sage: set_random_seed()
 477         sage: K1 = random_cone(max_ambient_dim=8,
 478         ....:                  strictly_convex=True,
 479         ....:                  solid=True)
 480         sage: K2 = random_cone(max_ambient_dim=8,
 481         ....:                  strictly_convex=True,
 482         ....:                  solid=True)
 483         sage: K = K1.cartesian_product(K2)
 484         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 485         True
 486
 487     The Lyapunov rank is invariant under a linear isomorphism
 488     [Orlitzky/Gowda]_::
 489
 490         sage: K1 = random_cone(max_ambient_dim = 8)
 491         sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
 492         sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
 493         sage: lyapunov_rank(K1) == lyapunov_rank(K2)
 494         True
 495
 496     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 497     itself [Rudolf et al.]_::
 498
 499         sage: set_random_seed()
 500         sage: K = random_cone(max_ambient_dim=8)
 501         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 502         True
 503
 504     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 505     be any number between `1` and `n` inclusive, excluding `n-1`
 506     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 507     trivial cone in a trivial space as well. However, in zero dimensions,
 508     the Lyapunov rank of the trivial cone will be zero::
 509
 510         sage: set_random_seed()
 511         sage: K = random_cone(max_ambient_dim=8,
 512         ....:                 strictly_convex=True,
 513         ....:                 solid=True)
 514         sage: b = lyapunov_rank(K)
 515         sage: n = K.lattice_dim()
 516         sage: (n == 0 or 1 <= b) and b <= n
 517         True
 518         sage: b == n-1
 519         False
 520
 521     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 522     Lyapunov rank `n-1` in `n` dimensions::
 523
 524         sage: set_random_seed()
 525         sage: K = random_cone(max_ambient_dim=8)
 526         sage: b = lyapunov_rank(K)
 527         sage: n = K.lattice_dim()
 528         sage: b == n-1
 529         False
 530
 531     The calculation of the Lyapunov rank of an improper cone can be
 532     reduced to that of a proper cone [Orlitzky/Gowda]_::
 533
 534         sage: set_random_seed()
 535         sage: K = random_cone(max_ambient_dim=8)
 536         sage: actual = lyapunov_rank(K)
 537         sage: K_S = _restrict_to_space(K, K.span())
 538         sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
 539         sage: l = K.lineality()
 540         sage: c = K.codim()
 541         sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
 542         sage: actual == expected
 543         True
 544
 545     The Lyapunov rank of any cone is just the dimension of ``LL(K)``::
 546
 547         sage: set_random_seed()
 548         sage: K = random_cone(max_ambient_dim=8)
 549         sage: lyapunov_rank(K) == len(LL(K))
 550         True
 551
 552     We can make an imperfect cone perfect by adding a slack variable
 553     (a Theorem in [Orlitzky/Gowda]_)::
 554
 555         sage: set_random_seed()
 556         sage: K = random_cone(max_ambient_dim=8,
 557         ....:                 strictly_convex=True,
 558         ....:                 solid=True)
 559         sage: L = ToricLattice(K.lattice_dim() + 1)
 560         sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
 561         sage: lyapunov_rank(K) >= K.lattice_dim()
 562         True
 563
 564     """
 565     beta = 0
 566
 567     m = K.dim()
 568     n = K.lattice_dim()
 569     l = K.lineality()
 570
 571     if m < n:
 572         # K is not solid, restrict to its span.
 573         K = _restrict_to_space(K, K.span())
 574
 575         # Non-solid reduction lemma.
 576         beta += (n - m)*n
 577
 578     if l > 0:
 579         # K is not pointed, restrict to the span of its dual. Uses a
 580         # proposition from our paper, i.e. this is equivalent to K =
 581         # _rho(K.dual()).dual().
 582         K = _restrict_to_space(K, K.dual().span())
 583
 584         # Non-pointed reduction lemma.
 585         beta += l * m
 586
 587     beta += len(LL(K))
 588     return beta
 589
 590
 591
 592 def is_lyapunov_like(L,K):
 593     r"""
 594     Determine whether or not ``L`` is Lyapunov-like on ``K``.
 595
 596     We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
 597     L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
 598     `\left\langle x,s \right\rangle` in the complementarity set of
 599     ``K``. It is known [Orlitzky]_ that this property need only be
 600     checked for generators of ``K`` and its dual.
 601
 602     INPUT:
 603
 604     - ``L`` -- A linear transformation or matrix.
 605
 606     - ``K`` -- A polyhedral closed convex cone.
 607
 608     OUTPUT:
 609
 610     ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
 611     and ``False`` otherwise.
 612
 613     .. WARNING::
 614
 615         If this function returns ``True``, then ``L`` is Lyapunov-like
 616         on ``K``. However, if ``False`` is returned, that could mean one
 617         of two things. The first is that ``L`` is definitely not
 618         Lyapunov-like on ``K``. The second is more of an "I don't know"
 619         answer, returned (for example) if we cannot prove that an inner
 620         product is zero.
 621
 622     REFERENCES:
 623
 624     .. [Orlitzky] M. Orlitzky. The Lyapunov rank of an
 625        improper cone (preprint).
 626
 627     EXAMPLES:
 628
 629     The identity is always Lyapunov-like in a nontrivial space::
 630
 631         sage: set_random_seed()
 632         sage: K = random_cone(min_ambient_dim = 1, max_rays = 8)
 633         sage: L = identity_matrix(K.lattice_dim())
 634         sage: is_lyapunov_like(L,K)
 635         True
 636
 637     As is the "zero" transformation::
 638
 639         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
 640         sage: R = K.lattice().vector_space().base_ring()
 641         sage: L = zero_matrix(R, K.lattice_dim())
 642         sage: is_lyapunov_like(L,K)
 643         True
 644
 645     Everything in ``LL(K)`` should be Lyapunov-like on ``K``::
 646
 647         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
 648         sage: all([is_lyapunov_like(L,K) for L in LL(K)])
 649         True
 650
 651     """
 652     return all([(L*x).inner_product(s) == 0
 653                 for (x,s) in K.discrete_complementarity_set()])
 654
 655
 656 def random_element(K):
 657     r"""
 658     Return a random element of ``K`` from its ambient vector space.
 659
 660     ALGORITHM:
 661
 662     The cone ``K`` is specified in terms of its generators, so that
 663     ``K`` is equal to the convex conic combination of those generators.
 664     To choose a random element of ``K``, we assign random nonnegative
 665     coefficients to each generator of ``K`` and construct a new vector
 666     from the scaled rays.
 667
 668     A vector, rather than a ray, is returned so that the element may
 669     have non-integer coordinates. Thus the element may have an
 670     arbitrarily small norm.
 671
 672     EXAMPLES:
 673
 674     A random element of the trivial cone is zero::
 675
 676         sage: set_random_seed()
 677         sage: K = Cone([], ToricLattice(0))
 678         sage: random_element(K)
 679         ()
 680         sage: K = Cone([(0,)])
 681         sage: random_element(K)
 682         (0)
 683         sage: K = Cone([(0,0)])
 684         sage: random_element(K)
 685         (0, 0)
 686         sage: K = Cone([(0,0,0)])
 687         sage: random_element(K)
 688         (0, 0, 0)
 689
 690     TESTS:
 691
 692     Any cone should contain an element of itself::
 693
 694         sage: set_random_seed()
 695         sage: K = random_cone(max_rays = 8)
 696         sage: K.contains(random_element(K))
 697         True
 698
 699     """
 700     V = K.lattice().vector_space()
 701     F = V.base_ring()
 702     coefficients = [ F.random_element().abs() for i in range(K.nrays()) ]
 703     vector_gens  = map(V, K.rays())
 704     scaled_gens  = [ coefficients[i]*vector_gens[i]
 705                          for i in range(len(vector_gens)) ]
 706
 707     # Make sure we return a vector. Without the coercion, we might
 708     # return ``0`` when ``K`` has no rays.
 709     v = V(sum(scaled_gens))
 710     return v
 711
 712
 713 def positive_operators(K):
 714     r"""
 715     Compute generators of the cone of positive operators on this cone.
 716
 717     OUTPUT:
 718
 719     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 720     Each matrix ``P`` in the list should have the property that ``P*x``
 721     is an element of ``K`` whenever ``x`` is an element of
 722     ``K``. Moreover, any nonnegative linear combination of these
 723     matrices shares the same property.
 724
 725     EXAMPLES:
 726
 727     The trivial cone in a trivial space has no positive operators::
 728
 729         sage: K = Cone([], ToricLattice(0))
 730         sage: positive_operators(K)
 731         []
 732
 733     Positive operators on the nonnegative orthant are nonnegative matrices::
 734
 735         sage: K = Cone([(1,)])
 736         sage: positive_operators(K)
 737         [[1]]
 738
 739         sage: K = Cone([(1,0),(0,1)])
 740         sage: positive_operators(K)
 741         [
 742         [1 0]  [0 1]  [0 0]  [0 0]
 743         [0 0], [0 0], [1 0], [0 1]
 744         ]
 745
 746     Every operator is positive on the ambient vector space::
 747
 748         sage: K = Cone([(1,),(-1,)])
 749         sage: K.is_full_space()
 750         True
 751         sage: positive_operators(K)
 752         [[1], [-1]]
 753
 754         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 755         sage: K.is_full_space()
 756         True
 757         sage: positive_operators(K)
 758         [
 759         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 760         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 761         ]
 762
 763     TESTS:
 764
 765     A positive operator on a cone should send its generators into the cone::
 766
 767         sage: K = random_cone(max_ambient_dim = 6)
 768         sage: pi_of_k = positive_operators(K)
 769         sage: all([K.contains(p*x) for p in pi_of_k for x in K.rays()])
 770         True
 771
 772     """
 773     V = K.lattice().vector_space()
 774
 775     # Sage doesn't think matrices are vectors, so we have to convert
 776     # our matrices to vectors explicitly before we can figure out how
 777     # many are linearly-indepenedent.
 778     #
 779     # The space W has the same base ring as V, but dimension
 780     # dim(V)^2. So it has the same dimension as the space of linear
 781     # transformations on V. In other words, it's just the right size
 782     # to create an isomorphism between it and our matrices.
 783     W = VectorSpace(V.base_ring(), V.dimension()**2)
 784
 785     G1 = [ V(x) for x in K.rays() ]
 786     G2 = [ V(s) for s in K.dual().rays() ]
 787
 788     tensor_products = [ s.tensor_product(x) for x in G1 for s in G2 ]
 789
 790     # Turn our matrices into long vectors...
 791     vectors = [ W(m.list()) for m in tensor_products ]
 792
 793     # Create the *dual* cone of the positive operators, expressed as
 794     # long vectors..
 795     L = ToricLattice(W.dimension())
 796     pi_dual = Cone(vectors, lattice=L)
 797
 798     # Now compute the desired cone from its dual...
 799     pi_cone = pi_dual.dual()
 800
 801     # And finally convert its rays back to matrix representations.
 802     M = MatrixSpace(V.base_ring(), V.dimension())
 803
 804     return [ M(v.list()) for v in pi_cone.rays() ]