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 drop_dependent(vs):
  12     r"""
  13     Return the largest linearly-independent subset of ``vs``.
  14     """
  15     if len(vs) == 0:
  16         # ...for lazy enough definitions of linearly-independent
  17         return vs
  18
  19     result = []
  20     old_V = VectorSpace(vs[0].parent().base_field(), 0)
  21
  22     for v in vs:
  23         new_V = span(result + [v])
  24         if new_V.dimension() > old_V.dimension():
  25             result.append(v)
  26             old_V = new_V
  27
  28     return result
  29
  30
  31 def basically_the_same(K1,K2):
  32     r"""
  33     ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
  34     otherwise.
  35     """
  36     if K1.lattice_dim() != K2.lattice_dim():
  37         return False
  38
  39     if K1.nrays() != K2.nrays():
  40         return False
  41
  42     if K1.dim() != K2.dim():
  43         return False
  44
  45     if lineality(K1) != lineality(K2):
  46         return False
  47
  48     if K1.is_solid() != K2.is_solid():
  49         return False
  50
  51     if K1.is_strictly_convex() != K2.is_strictly_convex():
  52         return False
  53
  54     if len(LL(K1)) != len(LL(K2)):
  55         return False
  56
  57     C_of_K1 = discrete_complementarity_set(K1)
  58     C_of_K2 = discrete_complementarity_set(K2)
  59     if len(C_of_K1) != len(C_of_K2):
  60         return False
  61
  62     if len(K1.facets()) != len(K2.facets()):
  63         return False
  64
  65     return True
  66
  67
  68
  69 def iso_space(K):
  70     r"""
  71     Construct the space `W \times W^{\perp}` isomorphic to the ambient space
  72     of ``K`` where `W` is equal to the span of ``K``.
  73     """
  74     V = K.lattice().vector_space()
  75
  76     # Create the space W \times W^{\perp} isomorphic to V.
  77     # First we get an orthogonal (but not normal) basis...
  78     M = matrix(V.base_field(), K.rays())
  79     W_basis = drop_dependent(K.rays())
  80
  81     W = V.subspace_with_basis(W_basis)
  82     W_perp = W.complement()
  83
  84     return W.cartesian_product(W_perp)
  85
  86
  87 def ips_iso(K):
  88     r"""
  89     Construct the IPS isomorphism and its inverse from our paper.
  90
  91     Given a cone ``K``, the returned isomorphism will split its ambient
  92     vector space `V` into a cartesian product `W \times W^{\perp}` where
  93     `W` equals the span of ``K``.
  94     """
  95     V = K.lattice().vector_space()
  96     V_iso = iso_space(K)
  97     (W, W_perp) = V_iso.cartesian_factors()
  98
  99     # A space equivalent to V, but using our basis.
 100     V_user = V.subspace_with_basis( W.basis() + W_perp.basis() )
 101
 102     def phi(v):
 103         # Write v in terms of our custom basis, where the first dim(W)
 104         # coordinates are for the W-part of the basis.
 105         cs = V_user.coordinates(v)
 106
 107         w1 = sum([ V_user.basis()[idx]*cs[idx]
 108                     for idx in range(0, W.dimension()) ])
 109         w2 = sum([ V_user.basis()[idx]*cs[idx]
 110                     for idx in range(W.dimension(), V.dimension()) ])
 111
 112         return V_iso( (w1, w2) )
 113
 114
 115     def phi_inv( pair ):
 116         # Crash if the arguments are in the wrong spaces.
 117         V_iso(pair)
 118
 119         #w = sum([ sub_w[idx]*W.basis()[idx] for idx in range(0,m) ])
 120         #w_prime = sum([ sub_w_prime[idx]*W_perp.basis()[idx]
 121         #             for idx in range(0,n-m) ])
 122
 123         return sum( pair.cartesian_factors() )
 124
 125
 126     return (phi,phi_inv)
 127
 128
 129 def rho(K, K2=None):
 130     r"""
 131     Restrict ``K`` into its own span, or the span of another cone.
 132
 133     INPUT:
 134
 135     - ``K2`` -- another cone whose lattice has the same rank as this cone.
 136
 137     OUTPUT:
 138
 139     A new cone in a sublattice.
 140
 141     EXAMPLES::
 142
 143         sage: K = Cone([(1,)])
 144         sage: restrict_span(K) == K
 145         True
 146
 147         sage: K2 = Cone([(1,0)])
 148         sage: restrict_span(K2).rays()
 149         N(1)
 150         in 1-d lattice N
 151         sage: K3 = Cone([(1,0,0)])
 152         sage: restrict_span(K3).rays()
 153         N(1)
 154         in 1-d lattice N
 155         sage: restrict_span(K2) == restrict_span(K3)
 156         True
 157
 158     TESTS:
 159
 160     The projected cone should always be solid::
 161
 162         sage: set_random_seed()
 163         sage: K = random_cone(max_dim = 8)
 164         sage: K_S = restrict_span(K)
 165         sage: K_S.is_solid()
 166         True
 167
 168     And the resulting cone should live in a space having the same
 169     dimension as the space we restricted it to::
 170
 171         sage: set_random_seed()
 172         sage: K = random_cone(max_dim = 8)
 173         sage: K_S = restrict_span(K, K.dual() )
 174         sage: K_S.lattice_dim() == K.dual().dim()
 175         True
 176
 177     This function should not affect the dimension of a cone::
 178
 179         sage: set_random_seed()
 180         sage: K = random_cone(max_dim = 8)
 181         sage: K.dim() == restrict_span(K).dim()
 182         True
 183
 184     Nor should it affect the lineality of a cone::
 185
 186         sage: set_random_seed()
 187         sage: K = random_cone(max_dim = 8)
 188         sage: lineality(K) == lineality(restrict_span(K))
 189         True
 190
 191     No matter which space we restrict to, the lineality should not
 192     increase::
 193
 194         sage: set_random_seed()
 195         sage: K = random_cone(max_dim = 8)
 196         sage: lineality(K) >= lineality(restrict_span(K))
 197         True
 198         sage: lineality(K) >= lineality(restrict_span(K, K.dual()))
 199         True
 200
 201     If we do this according to our paper, then the result is proper::
 202
 203         sage: set_random_seed()
 204         sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=False)
 205         sage: K_S = restrict_span(K)
 206         sage: P = restrict_span(K_S.dual()).dual()
 207         sage: P.is_proper()
 208         True
 209         sage: P = restrict_span(K_S, K_S.dual())
 210         sage: P.is_proper()
 211         True
 212
 213     ::
 214
 215         sage: set_random_seed()
 216         sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=False)
 217         sage: K_S = restrict_span(K)
 218         sage: P = restrict_span(K_S.dual()).dual()
 219         sage: P.is_proper()
 220         True
 221         sage: P = restrict_span(K_S, K_S.dual())
 222         sage: P.is_proper()
 223         True
 224
 225     ::
 226
 227         sage: set_random_seed()
 228         sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=True)
 229         sage: K_S = restrict_span(K)
 230         sage: P = restrict_span(K_S.dual()).dual()
 231         sage: P.is_proper()
 232         True
 233         sage: P = restrict_span(K_S, K_S.dual())
 234         sage: P.is_proper()
 235         True
 236
 237     ::
 238
 239         sage: set_random_seed()
 240         sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=True)
 241         sage: K_S = restrict_span(K)
 242         sage: P = restrict_span(K_S.dual()).dual()
 243         sage: P.is_proper()
 244         True
 245         sage: P = restrict_span(K_S, K_S.dual())
 246         sage: P.is_proper()
 247         True
 248
 249     Test the proposition in our paper concerning the duals, where the
 250     subspace `W` is the span of `K^{*}`::
 251
 252         sage: set_random_seed()
 253         sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=False)
 254         sage: K_W = restrict_span(K, K.dual())
 255         sage: K_star_W_star = restrict_span(K.dual()).dual()
 256         sage: basically_the_same(K_W, K_star_W_star)
 257         True
 258
 259     ::
 260
 261         sage: set_random_seed()
 262         sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=False)
 263         sage: K_W = restrict_span(K, K.dual())
 264         sage: K_star_W_star = restrict_span(K.dual()).dual()
 265         sage: basically_the_same(K_W, K_star_W_star)
 266         True
 267
 268     ::
 269
 270         sage: set_random_seed()
 271         sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=True)
 272         sage: K_W = restrict_span(K, K.dual())
 273         sage: K_star_W_star = restrict_span(K.dual()).dual()
 274         sage: basically_the_same(K_W, K_star_W_star)
 275         True
 276
 277     ::
 278
 279         sage: set_random_seed()
 280         sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=True)
 281         sage: K_W = restrict_span(K, K.dual())
 282         sage: K_star_W_star = restrict_span(K.dual()).dual()
 283         sage: basically_the_same(K_W, K_star_W_star)
 284         True
 285
 286     """
 287     if K2 is None:
 288         K2 = K
 289
 290     phi,_ = ips_iso(K2)
 291     (W, W_perp) = iso_space(K2).cartesian_factors()
 292
 293     ray_pairs = [ phi(r) for r in K.rays() ]
 294
 295     # Shouldn't matter?
 296     #
 297     #if any([ w2 != W_perp.zero() for (_, w2) in ray_pairs ]):
 298     #    msg = 'Cone has nonzero components in W-perp!'
 299     #    raise ValueError(msg)
 300
 301     # Represent the cone in terms of a basis for W, i.e. with smaller
 302     # vectors.
 303     ws = [ W.coordinate_vector(w1) for (w1, _) in ray_pairs ]
 304
 305     L = ToricLattice(W.dimension())
 306
 307     return Cone(ws, lattice=L)
 308
 309
 310
 311 def lineality(K):
 312     r"""
 313     Compute the lineality of this cone.
 314
 315     The lineality of a cone is the dimension of the largest linear
 316     subspace contained in that cone.
 317
 318     OUTPUT:
 319
 320     A nonnegative integer; the dimension of the largest subspace
 321     contained within this cone.
 322
 323     REFERENCES:
 324
 325     .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
 326        University Press, Princeton, 1970.
 327
 328     EXAMPLES:
 329
 330     The lineality of the nonnegative orthant is zero, since it clearly
 331     contains no lines::
 332
 333         sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
 334         sage: lineality(K)
 335         0
 336
 337     However, if we add another ray so that the entire `x`-axis belongs
 338     to the cone, then the resulting cone will have lineality one::
 339
 340         sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
 341         sage: lineality(K)
 342         1
 343
 344     If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
 345     to the dimension of the ambient space (i.e. two)::
 346
 347         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 348         sage: lineality(K)
 349         2
 350
 351     Per the definition, the lineality of the trivial cone in a trivial
 352     space is zero::
 353
 354         sage: K = Cone([], lattice=ToricLattice(0))
 355         sage: lineality(K)
 356         0
 357
 358     TESTS:
 359
 360     The lineality of a cone should be an integer between zero and the
 361     dimension of the ambient space, inclusive::
 362
 363         sage: set_random_seed()
 364         sage: K = random_cone(max_dim = 8)
 365         sage: l = lineality(K)
 366         sage: l in ZZ
 367         True
 368         sage: (0 <= l) and (l <= K.lattice_dim())
 369         True
 370
 371     A strictly convex cone should have lineality zero::
 372
 373         sage: set_random_seed()
 374         sage: K = random_cone(max_dim = 8, strictly_convex = True)
 375         sage: lineality(K)
 376         0
 377
 378     """
 379     return K.linear_subspace().dimension()
 380
 381
 382 def discrete_complementarity_set(K):
 383     r"""
 384     Compute the discrete complementarity set of this cone.
 385
 386     The complementarity set of this cone is the set of all orthogonal
 387     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
 388     dual. The discrete complementarity set restricts `x` and `s` to be
 389     generators of their respective cones.
 390
 391     OUTPUT:
 392
 393     A list of pairs `(x,s)` such that,
 394
 395       * `x` is in this cone.
 396       * `x` is a generator of this cone.
 397       * `s` is in this cone's dual.
 398       * `s` is a generator of this cone's dual.
 399       * `x` and `s` are orthogonal.
 400
 401     EXAMPLES:
 402
 403     The discrete complementarity set of the nonnegative orthant consists
 404     of pairs of standard basis vectors::
 405
 406         sage: K = Cone([(1,0),(0,1)])
 407         sage: discrete_complementarity_set(K)
 408         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
 409
 410     If the cone consists of a single ray, the second components of the
 411     discrete complementarity set should generate the orthogonal
 412     complement of that ray::
 413
 414         sage: K = Cone([(1,0)])
 415         sage: discrete_complementarity_set(K)
 416         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
 417         sage: K = Cone([(1,0,0)])
 418         sage: discrete_complementarity_set(K)
 419         [((1, 0, 0), (0, 1, 0)),
 420           ((1, 0, 0), (0, -1, 0)),
 421           ((1, 0, 0), (0, 0, 1)),
 422           ((1, 0, 0), (0, 0, -1))]
 423
 424     When the cone is the entire space, its dual is the trivial cone, so
 425     the discrete complementarity set is empty::
 426
 427         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 428         sage: discrete_complementarity_set(K)
 429         []
 430
 431     TESTS:
 432
 433     The complementarity set of the dual can be obtained by switching the
 434     components of the complementarity set of the original cone::
 435
 436         sage: set_random_seed()
 437         sage: K1 = random_cone(max_dim=6)
 438         sage: K2 = K1.dual()
 439         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 440         sage: actual = discrete_complementarity_set(K1)
 441         sage: sorted(actual) == sorted(expected)
 442         True
 443
 444     """
 445     V = K.lattice().vector_space()
 446
 447     # Convert the rays to vectors so that we can compute inner
 448     # products.
 449     xs = [V(x) for x in K.rays()]
 450     ss = [V(s) for s in K.dual().rays()]
 451
 452     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 453
 454
 455 def LL(K):
 456     r"""
 457     Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
 458     on this cone.
 459
 460     OUTPUT:
 461
 462     A list of matrices forming a basis for the space of all
 463     Lyapunov-like transformations on the given cone.
 464
 465     EXAMPLES:
 466
 467     The trivial cone has no Lyapunov-like transformations::
 468
 469         sage: L = ToricLattice(0)
 470         sage: K = Cone([], lattice=L)
 471         sage: LL(K)
 472         []
 473
 474     The Lyapunov-like transformations on the nonnegative orthant are
 475     simply diagonal matrices::
 476
 477         sage: K = Cone([(1,)])
 478         sage: LL(K)
 479         [[1]]
 480
 481         sage: K = Cone([(1,0),(0,1)])
 482         sage: LL(K)
 483         [
 484         [1 0]  [0 0]
 485         [0 0], [0 1]
 486         ]
 487
 488         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 489         sage: LL(K)
 490         [
 491         [1 0 0]  [0 0 0]  [0 0 0]
 492         [0 0 0]  [0 1 0]  [0 0 0]
 493         [0 0 0], [0 0 0], [0 0 1]
 494         ]
 495
 496     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 497     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 498
 499         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 500         sage: LL(L31)
 501         [
 502         [1 0 0]
 503         [0 1 0]
 504         [0 0 1]
 505         ]
 506
 507         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 508         sage: LL(L3infty)
 509         [
 510         [1 0 0]
 511         [0 1 0]
 512         [0 0 1]
 513         ]
 514
 515     If our cone is the entire space, then every transformation on it is
 516     Lyapunov-like::
 517
 518         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 519         sage: M = MatrixSpace(QQ,2)
 520         sage: M.basis() == LL(K)
 521         True
 522
 523     TESTS:
 524
 525     The inner product `\left< L\left(x\right), s \right>` is zero for
 526     every pair `\left( x,s \right)` in the discrete complementarity set
 527     of the cone::
 528
 529         sage: set_random_seed()
 530         sage: K = random_cone(max_dim=8)
 531         sage: C_of_K = discrete_complementarity_set(K)
 532         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 533         sage: sum(map(abs, l))
 534         0
 535
 536     The Lyapunov-like transformations on a cone and its dual are related
 537     by transposition, but we're not guaranteed to compute transposed
 538     elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
 539     \right)`
 540
 541         sage: set_random_seed()
 542         sage: K = random_cone(max_dim=8)
 543         sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
 544         sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
 545         sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
 546         sage: LL2_vecs = [ V(m.list()) for m in LL2   ]
 547         sage: V.span(LL1_vecs) == V.span(LL2_vecs)
 548         True
 549
 550     """
 551     V = K.lattice().vector_space()
 552
 553     C_of_K = discrete_complementarity_set(K)
 554
 555     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 556
 557     # Sage doesn't think matrices are vectors, so we have to convert
 558     # our matrices to vectors explicitly before we can figure out how
 559     # many are linearly-indepenedent.
 560     #
 561     # The space W has the same base ring as V, but dimension
 562     # dim(V)^2. So it has the same dimension as the space of linear
 563     # transformations on V. In other words, it's just the right size
 564     # to create an isomorphism between it and our matrices.
 565     W = VectorSpace(V.base_ring(), V.dimension()**2)
 566
 567     # Turn our matrices into long vectors...
 568     vectors = [ W(m.list()) for m in tensor_products ]
 569
 570     # Vector space representation of Lyapunov-like matrices
 571     # (i.e. vec(L) where L is Luapunov-like).
 572     LL_vector = W.span(vectors).complement()
 573
 574     # Now construct an ambient MatrixSpace in which to stick our
 575     # transformations.
 576     M = MatrixSpace(V.base_ring(), V.dimension())
 577
 578     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 579
 580     return matrix_basis
 581
 582
 583
 584 def lyapunov_rank(K):
 585     r"""
 586     Compute the Lyapunov (or bilinearity) rank of this cone.
 587
 588     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 589
 590     1. The dimension of the Lie algebra of the automorphism group of the
 591        cone.
 592
 593     2. The dimension of the linear space of all Lyapunov-like
 594        transformations on the cone.
 595
 596     INPUT:
 597
 598     A closed, convex polyhedral cone.
 599
 600     OUTPUT:
 601
 602     An integer representing the Lyapunov rank of the cone. If the
 603     dimension of the ambient vector space is `n`, then the Lyapunov rank
 604     will be between `1` and `n` inclusive; however a rank of `n-1` is
 605     not possible (see the first reference).
 606
 607     .. note::
 608
 609         In the references, the cones are always assumed to be proper. We
 610         do not impose this restriction.
 611
 612     .. seealso::
 613
 614         :meth:`is_proper`
 615
 616     ALGORITHM:
 617
 618     The codimension formula from the second reference is used. We find
 619     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 620     and `s` are rays of our cone. It is known that these vectors are
 621     sufficient to apply the codimension formula. Once we have all such
 622     pairs, we "brute force" the codimension formula by finding all
 623     linearly-independent `xs^{T}`.
 624
 625     REFERENCES:
 626
 627     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 628        cone and Lyapunov-like transformations, Mathematical Programming, 147
 629        (2014) 155-170.
 630
 631     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 632        Improper Cone. Work in-progress.
 633
 634     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 635        optimality constraints for the cone of positive polynomials,
 636        Mathematical Programming, Series B, 129 (2011) 5-31.
 637
 638     EXAMPLES:
 639
 640     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 641     [Rudolf et al.]_::
 642
 643         sage: positives = Cone([(1,)])
 644         sage: lyapunov_rank(positives)
 645         1
 646         sage: quadrant = Cone([(1,0), (0,1)])
 647         sage: lyapunov_rank(quadrant)
 648         2
 649        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 650         sage: lyapunov_rank(octant)
 651         3
 652
 653     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 654     [Orlitzky/Gowda]_::
 655
 656         sage: R5 = VectorSpace(QQ, 5)
 657         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 658         sage: K = Cone(gs)
 659         sage: lyapunov_rank(K)
 660         25
 661
 662     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 663     [Rudolf et al.]_::
 664
 665         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 666         sage: lyapunov_rank(L31)
 667         1
 668
 669     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 670
 671         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 672         sage: lyapunov_rank(L3infty)
 673         1
 674
 675     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 676     + 1` [Orlitzky/Gowda]_::
 677
 678         sage: K = Cone([(1,0,0,0,0)])
 679         sage: lyapunov_rank(K)
 680         21
 681         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 682         21
 683
 684     A subspace (of dimension `m`) in `n` dimensions should have a
 685     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 686
 687         sage: e1 = (1,0,0,0,0)
 688         sage: neg_e1 = (-1,0,0,0,0)
 689         sage: e2 = (0,1,0,0,0)
 690         sage: neg_e2 = (0,-1,0,0,0)
 691         sage: z = (0,0,0,0,0)
 692         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 693         sage: lyapunov_rank(K)
 694         19
 695         sage: K.lattice_dim()**2 - K.dim()*codim(K)
 696         19
 697
 698     The Lyapunov rank should be additive on a product of proper cones
 699     [Rudolf et al.]_::
 700
 701         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 702         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 703         sage: K = L31.cartesian_product(octant)
 704         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 705         True
 706
 707     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 708     The cone ``K`` in the following example is isomorphic to the nonnegative
 709     octant in `\mathbb{R}^{3}`::
 710
 711         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 712         sage: lyapunov_rank(K)
 713         3
 714
 715     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 716     itself [Rudolf et al.]_::
 717
 718         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 719         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 720         True
 721
 722     TESTS:
 723
 724     The Lyapunov rank should be additive on a product of proper cones
 725     [Rudolf et al.]_::
 726
 727         sage: set_random_seed()
 728         sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
 729         sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
 730         sage: K = K1.cartesian_product(K2)
 731         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 732         True
 733
 734     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 735     itself [Rudolf et al.]_::
 736
 737         sage: set_random_seed()
 738         sage: K = random_cone(max_dim=8)
 739         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 740         True
 741
 742     Make sure we exercise the non-strictly-convex/non-solid case::
 743
 744         sage: set_random_seed()
 745         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
 746         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 747         True
 748
 749     Let's check the other permutations as well, just to be sure::
 750
 751         sage: set_random_seed()
 752         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
 753         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 754         True
 755
 756     ::
 757
 758         sage: set_random_seed()
 759         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
 760         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 761         True
 762
 763     ::
 764
 765         sage: set_random_seed()
 766         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 767         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 768         True
 769
 770     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 771     be any number between `1` and `n` inclusive, excluding `n-1`
 772     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 773     trivial cone in a trivial space as well. However, in zero dimensions,
 774     the Lyapunov rank of the trivial cone will be zero::
 775
 776         sage: set_random_seed()
 777         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 778         sage: b = lyapunov_rank(K)
 779         sage: n = K.lattice_dim()
 780         sage: (n == 0 or 1 <= b) and b <= n
 781         True
 782         sage: b == n-1
 783         False
 784
 785     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 786     Lyapunov rank `n-1` in `n` dimensions::
 787
 788         sage: set_random_seed()
 789         sage: K = random_cone(max_dim=8)
 790         sage: b = lyapunov_rank(K)
 791         sage: n = K.lattice_dim()
 792         sage: b == n-1
 793         False
 794
 795     The calculation of the Lyapunov rank of an improper cone can be
 796     reduced to that of a proper cone [Orlitzky/Gowda]_::
 797
 798         sage: set_random_seed()
 799         sage: K = random_cone(max_dim=8)
 800         sage: actual = lyapunov_rank(K)
 801         sage: K_S = restrict_span(K)
 802         sage: P = restrict_span(K_S.dual()).dual()
 803         sage: l = lineality(K)
 804         sage: c = codim(K)
 805         sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2
 806         sage: actual == expected
 807         True
 808
 809     The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
 810
 811         sage: set_random_seed()
 812         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 813         sage: lyapunov_rank(K) == len(LL(K))
 814         True
 815
 816     In fact the same can be said of any cone. These additional tests
 817     just increase our confidence that the reduction scheme works::
 818
 819         sage: set_random_seed()
 820         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
 821         sage: lyapunov_rank(K) == len(LL(K))
 822         True
 823
 824     ::
 825
 826         sage: set_random_seed()
 827         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
 828         sage: lyapunov_rank(K) == len(LL(K))
 829         True
 830
 831     ::
 832
 833         sage: set_random_seed()
 834         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
 835         sage: lyapunov_rank(K) == len(LL(K))
 836         True
 837
 838     """
 839     K_orig = K
 840     beta = 0
 841
 842     m = K.dim()
 843     n = K.lattice_dim()
 844     l = lineality(K)
 845
 846     if m < n:
 847         # K is not solid, project onto its span.
 848         K = restrict_span(K)
 849
 850         # Lemma 2
 851         beta += m*(n - m) + (n - m)**2
 852
 853     if l > 0:
 854         # K is not pointed, project its dual onto its span.
 855         # Uses a proposition from our paper, i.e. this is
 856         # equivalent to K = restrict_span(K.dual()).dual()
 857         #K = restrict_span(intersect_span(K,K.dual()), K.dual())
 858         K = restrict_span(K, K.dual())
 859
 860         # Lemma 3
 861         beta += m * l
 862
 863     beta += len(LL(K))
 864     return beta