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