# Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
# have to explicitly mangle our sitedir here so that "mjo.cone"
# resolves.
from os.path import abspath
from site import addsitedir
addsitedir(abspath('../../'))

from sage.all import *


def drop_dependent(vs):
    r"""
    Return the largest linearly-independent subset of ``vs``.
    """
    if len(vs) == 0:
        # ...for lazy enough definitions of linearly-independent
        return vs

    result = []
    old_V = VectorSpace(vs[0].parent().base_field(), 0)

    for v in vs:
        new_V = span(result + [v])
        if new_V.dimension() > old_V.dimension():
            result.append(v)
            old_V = new_V

    return result


def basically_the_same(K1,K2):
    r"""
    ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
    otherwise.
    """
    if K1.lattice_dim() != K2.lattice_dim():
        return False

    if K1.nrays() != K2.nrays():
        return False

    if K1.dim() != K2.dim():
        return False

    if lineality(K1) != lineality(K2):
        return False

    if K1.is_solid() != K2.is_solid():
        return False

    if K1.is_strictly_convex() != K2.is_strictly_convex():
        return False

    if len(LL(K1)) != len(LL(K2)):
        return False

    C_of_K1 = discrete_complementarity_set(K1)
    C_of_K2 = discrete_complementarity_set(K2)
    if len(C_of_K1) != len(C_of_K2):
        return False

    if len(K1.facets()) != len(K2.facets()):
        return False

    return True



def iso_space(K):
    r"""
    Construct the space `W \times W^{\perp}` isomorphic to the ambient space
    of ``K`` where `W` is equal to the span of ``K``.
    """
    V = K.lattice().vector_space()

    # Create the space W \times W^{\perp} isomorphic to V.
    # First we get an orthogonal (but not normal) basis...
    M = matrix(V.base_field(), K.rays())
    W_basis = drop_dependent(K.rays())

    W = V.subspace_with_basis(W_basis)
    W_perp = W.complement()

    return W.cartesian_product(W_perp)


def ips_iso(K):
    r"""
    Construct the IPS isomorphism and its inverse from our paper.

    Given a cone ``K``, the returned isomorphism will split its ambient
    vector space `V` into a cartesian product `W \times W^{\perp}` where
    `W` equals the span of ``K``.
    """
    V = K.lattice().vector_space()
    V_iso = iso_space(K)
    (W, W_perp) = V_iso.cartesian_factors()

    # A space equivalent to V, but using our basis.
    V_user = V.subspace_with_basis( W.basis() + W_perp.basis() )

    def phi(v):
        # Write v in terms of our custom basis, where the first dim(W)
        # coordinates are for the W-part of the basis.
        cs = V_user.coordinates(v)

        w1 = sum([ V_user.basis()[idx]*cs[idx]
                    for idx in range(0, W.dimension()) ])
        w2 = sum([ V_user.basis()[idx]*cs[idx]
                    for idx in range(W.dimension(), V.dimension()) ])

        return V_iso( (w1, w2) )


    def phi_inv( pair ):
        # Crash if the arguments are in the wrong spaces.
        V_iso(pair)

        #w = sum([ sub_w[idx]*W.basis()[idx] for idx in range(0,m) ])
        #w_prime = sum([ sub_w_prime[idx]*W_perp.basis()[idx]
        #             for idx in range(0,n-m) ])

        return sum( pair.cartesian_factors() )


    return (phi,phi_inv)


def rho(K, K2=None):
    r"""
    Restrict ``K`` into its own span, or the span of another cone.

    INPUT:

    - ``K2`` -- another cone whose lattice has the same rank as this cone.

    OUTPUT:

    A new cone in a sublattice.

    EXAMPLES::

        sage: K = Cone([(1,)])
        sage: restrict_span(K) == K
        True

        sage: K2 = Cone([(1,0)])
        sage: restrict_span(K2).rays()
        N(1)
        in 1-d lattice N
        sage: K3 = Cone([(1,0,0)])
        sage: restrict_span(K3).rays()
        N(1)
        in 1-d lattice N
        sage: restrict_span(K2) == restrict_span(K3)
        True

    TESTS:

    The projected cone should always be solid::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8)
        sage: K_S = restrict_span(K)
        sage: K_S.is_solid()
        True

    And the resulting cone should live in a space having the same
    dimension as the space we restricted it to::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8)
        sage: K_S = restrict_span(K, K.dual() )
        sage: K_S.lattice_dim() == K.dual().dim()
        True

    This function should not affect the dimension of a cone::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8)
        sage: K.dim() == restrict_span(K).dim()
        True

    Nor should it affect the lineality of a cone::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8)
        sage: lineality(K) == lineality(restrict_span(K))
        True

    No matter which space we restrict to, the lineality should not
    increase::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8)
        sage: lineality(K) >= lineality(restrict_span(K))
        True
        sage: lineality(K) >= lineality(restrict_span(K, K.dual()))
        True

    If we do this according to our paper, then the result is proper::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=False)
        sage: K_S = restrict_span(K)
        sage: P = restrict_span(K_S.dual()).dual()
        sage: P.is_proper()
        True
        sage: P = restrict_span(K_S, K_S.dual())
        sage: P.is_proper()
        True

    ::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=False)
        sage: K_S = restrict_span(K)
        sage: P = restrict_span(K_S.dual()).dual()
        sage: P.is_proper()
        True
        sage: P = restrict_span(K_S, K_S.dual())
        sage: P.is_proper()
        True

    ::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=True)
        sage: K_S = restrict_span(K)
        sage: P = restrict_span(K_S.dual()).dual()
        sage: P.is_proper()
        True
        sage: P = restrict_span(K_S, K_S.dual())
        sage: P.is_proper()
        True

    ::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=True)
        sage: K_S = restrict_span(K)
        sage: P = restrict_span(K_S.dual()).dual()
        sage: P.is_proper()
        True
        sage: P = restrict_span(K_S, K_S.dual())
        sage: P.is_proper()
        True

    Test the proposition in our paper concerning the duals, where the
    subspace `W` is the span of `K^{*}`::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=False)
        sage: K_W = restrict_span(K, K.dual())
        sage: K_star_W_star = restrict_span(K.dual()).dual()
        sage: basically_the_same(K_W, K_star_W_star)
        True

    ::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=False)
        sage: K_W = restrict_span(K, K.dual())
        sage: K_star_W_star = restrict_span(K.dual()).dual()
        sage: basically_the_same(K_W, K_star_W_star)
        True

    ::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=True)
        sage: K_W = restrict_span(K, K.dual())
        sage: K_star_W_star = restrict_span(K.dual()).dual()
        sage: basically_the_same(K_W, K_star_W_star)
        True

    ::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=True)
        sage: K_W = restrict_span(K, K.dual())
        sage: K_star_W_star = restrict_span(K.dual()).dual()
        sage: basically_the_same(K_W, K_star_W_star)
        True

    """
    if K2 is None:
        K2 = K

    phi,_ = ips_iso(K2)
    (W, W_perp) = iso_space(K2).cartesian_factors()

    ray_pairs = [ phi(r) for r in K.rays() ]

    # Shouldn't matter?
    #
    #if any([ w2 != W_perp.zero() for (_, w2) in ray_pairs ]):
    #    msg = 'Cone has nonzero components in W-perp!'
    #    raise ValueError(msg)

    # Represent the cone in terms of a basis for W, i.e. with smaller
    # vectors.
    ws = [ W.coordinate_vector(w1) for (w1, _) in ray_pairs ]

    L = ToricLattice(W.dimension())

    return Cone(ws, lattice=L)



def lineality(K):
    r"""
    Compute the lineality of this cone.

    The lineality of a cone is the dimension of the largest linear
    subspace contained in that cone.

    OUTPUT:

    A nonnegative integer; the dimension of the largest subspace
    contained within this cone.

    REFERENCES:

    .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
       University Press, Princeton, 1970.

    EXAMPLES:

    The lineality of the nonnegative orthant is zero, since it clearly
    contains no lines::

        sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
        sage: lineality(K)
        0

    However, if we add another ray so that the entire `x`-axis belongs
    to the cone, then the resulting cone will have lineality one::

        sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
        sage: lineality(K)
        1

    If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
    to the dimension of the ambient space (i.e. two)::

        sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
        sage: lineality(K)
        2

    Per the definition, the lineality of the trivial cone in a trivial
    space is zero::

        sage: K = Cone([], lattice=ToricLattice(0))
        sage: lineality(K)
        0

    TESTS:

    The lineality of a cone should be an integer between zero and the
    dimension of the ambient space, inclusive::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8)
        sage: l = lineality(K)
        sage: l in ZZ
        True
        sage: (0 <= l) and (l <= K.lattice_dim())
        True

    A strictly convex cone should have lineality zero::

        sage: set_random_seed()
        sage: K = random_cone(max_dim = 8, strictly_convex = True)
        sage: lineality(K)
        0

    """
    return K.linear_subspace().dimension()


def discrete_complementarity_set(K):
    r"""
    Compute the discrete complementarity set of this cone.

    The complementarity set of this cone is the set of all orthogonal
    pairs `(x,s)` such that `x` is in this cone, and `s` is in its
    dual. The discrete complementarity set restricts `x` and `s` to be
    generators of their respective cones.

    OUTPUT:

    A list of pairs `(x,s)` such that,

      * `x` is in this cone.
      * `x` is a generator of this cone.
      * `s` is in this cone's dual.
      * `s` is a generator of this cone's dual.
      * `x` and `s` are orthogonal.

    EXAMPLES:

    The discrete complementarity set of the nonnegative orthant consists
    of pairs of standard basis vectors::

        sage: K = Cone([(1,0),(0,1)])
        sage: discrete_complementarity_set(K)
        [((1, 0), (0, 1)), ((0, 1), (1, 0))]

    If the cone consists of a single ray, the second components of the
    discrete complementarity set should generate the orthogonal
    complement of that ray::

        sage: K = Cone([(1,0)])
        sage: discrete_complementarity_set(K)
        [((1, 0), (0, 1)), ((1, 0), (0, -1))]
        sage: K = Cone([(1,0,0)])
        sage: discrete_complementarity_set(K)
        [((1, 0, 0), (0, 1, 0)),
          ((1, 0, 0), (0, -1, 0)),
          ((1, 0, 0), (0, 0, 1)),
          ((1, 0, 0), (0, 0, -1))]

    When the cone is the entire space, its dual is the trivial cone, so
    the discrete complementarity set is empty::

        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
        sage: discrete_complementarity_set(K)
        []

    TESTS:

    The complementarity set of the dual can be obtained by switching the
    components of the complementarity set of the original cone::

        sage: set_random_seed()
        sage: K1 = random_cone(max_dim=6)
        sage: K2 = K1.dual()
        sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
        sage: actual = discrete_complementarity_set(K1)
        sage: sorted(actual) == sorted(expected)
        True

    """
    V = K.lattice().vector_space()

    # Convert the rays to vectors so that we can compute inner
    # products.
    xs = [V(x) for x in K.rays()]
    ss = [V(s) for s in K.dual().rays()]

    return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]


def LL(K):
    r"""
    Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
    on this cone.

    OUTPUT:

    A list of matrices forming a basis for the space of all
    Lyapunov-like transformations on the given cone.

    EXAMPLES:

    The trivial cone has no Lyapunov-like transformations::

        sage: L = ToricLattice(0)
        sage: K = Cone([], lattice=L)
        sage: LL(K)
        []

    The Lyapunov-like transformations on the nonnegative orthant are
    simply diagonal matrices::

        sage: K = Cone([(1,)])
        sage: LL(K)
        [[1]]

        sage: K = Cone([(1,0),(0,1)])
        sage: LL(K)
        [
        [1 0]  [0 0]
        [0 0], [0 1]
        ]

        sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
        sage: LL(K)
        [
        [1 0 0]  [0 0 0]  [0 0 0]
        [0 0 0]  [0 1 0]  [0 0 0]
        [0 0 0], [0 0 0], [0 0 1]
        ]

    Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
    `L^{3}_{\infty}` cones [Rudolf et al.]_::

        sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
        sage: LL(L31)
        [
        [1 0 0]
        [0 1 0]
        [0 0 1]
        ]

        sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
        sage: LL(L3infty)
        [
        [1 0 0]
        [0 1 0]
        [0 0 1]
        ]

    If our cone is the entire space, then every transformation on it is
    Lyapunov-like::

        sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
        sage: M = MatrixSpace(QQ,2)
        sage: M.basis() == LL(K)
        True

    TESTS:

    The inner product `\left< L\left(x\right), s \right>` is zero for
    every pair `\left( x,s \right)` in the discrete complementarity set
    of the cone::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8)
        sage: C_of_K = discrete_complementarity_set(K)
        sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
        sage: sum(map(abs, l))
        0

    The Lyapunov-like transformations on a cone and its dual are related
    by transposition, but we're not guaranteed to compute transposed
    elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
    \right)`

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8)
        sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
        sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
        sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
        sage: LL2_vecs = [ V(m.list()) for m in LL2   ]
        sage: V.span(LL1_vecs) == V.span(LL2_vecs)
        True

    """
    V = K.lattice().vector_space()

    C_of_K = discrete_complementarity_set(K)

    tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]

    # Sage doesn't think matrices are vectors, so we have to convert
    # our matrices to vectors explicitly before we can figure out how
    # many are linearly-indepenedent.
    #
    # The space W has the same base ring as V, but dimension
    # dim(V)^2. So it has the same dimension as the space of linear
    # transformations on V. In other words, it's just the right size
    # to create an isomorphism between it and our matrices.
    W = VectorSpace(V.base_ring(), V.dimension()**2)

    # Turn our matrices into long vectors...
    vectors = [ W(m.list()) for m in tensor_products ]

    # Vector space representation of Lyapunov-like matrices
    # (i.e. vec(L) where L is Luapunov-like).
    LL_vector = W.span(vectors).complement()

    # Now construct an ambient MatrixSpace in which to stick our
    # transformations.
    M = MatrixSpace(V.base_ring(), V.dimension())

    matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]

    return matrix_basis



def lyapunov_rank(K):
    r"""
    Compute the Lyapunov (or bilinearity) rank of this cone.

    The Lyapunov rank of a cone can be thought of in (mainly) two ways:

    1. The dimension of the Lie algebra of the automorphism group of the
       cone.

    2. The dimension of the linear space of all Lyapunov-like
       transformations on the cone.

    INPUT:

    A closed, convex polyhedral cone.

    OUTPUT:

    An integer representing the Lyapunov rank of the cone. If the
    dimension of the ambient vector space is `n`, then the Lyapunov rank
    will be between `1` and `n` inclusive; however a rank of `n-1` is
    not possible (see the first reference).

    .. note::

        In the references, the cones are always assumed to be proper. We
        do not impose this restriction.

    .. seealso::

        :meth:`is_proper`

    ALGORITHM:

    The codimension formula from the second reference is used. We find
    all pairs `(x,s)` in the complementarity set of `K` such that `x`
    and `s` are rays of our cone. It is known that these vectors are
    sufficient to apply the codimension formula. Once we have all such
    pairs, we "brute force" the codimension formula by finding all
    linearly-independent `xs^{T}`.

    REFERENCES:

    .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
       cone and Lyapunov-like transformations, Mathematical Programming, 147
       (2014) 155-170.

    .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
       Improper Cone. Work in-progress.

    .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
       optimality constraints for the cone of positive polynomials,
       Mathematical Programming, Series B, 129 (2011) 5-31.

    EXAMPLES:

    The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
    [Rudolf et al.]_::

        sage: positives = Cone([(1,)])
        sage: lyapunov_rank(positives)
        1
        sage: quadrant = Cone([(1,0), (0,1)])
        sage: lyapunov_rank(quadrant)
        2
       sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
        sage: lyapunov_rank(octant)
        3

    The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
    [Orlitzky/Gowda]_::

        sage: R5 = VectorSpace(QQ, 5)
        sage: gs = R5.basis() + [ -r for r in R5.basis() ]
        sage: K = Cone(gs)
        sage: lyapunov_rank(K)
        25

    The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
    [Rudolf et al.]_::

        sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
        sage: lyapunov_rank(L31)
        1

    Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::

        sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
        sage: lyapunov_rank(L3infty)
        1

    A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
    + 1` [Orlitzky/Gowda]_::

        sage: K = Cone([(1,0,0,0,0)])
        sage: lyapunov_rank(K)
        21
        sage: K.lattice_dim()**2 - K.lattice_dim() + 1
        21

    A subspace (of dimension `m`) in `n` dimensions should have a
    Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::

        sage: e1 = (1,0,0,0,0)
        sage: neg_e1 = (-1,0,0,0,0)
        sage: e2 = (0,1,0,0,0)
        sage: neg_e2 = (0,-1,0,0,0)
        sage: z = (0,0,0,0,0)
        sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
        sage: lyapunov_rank(K)
        19
        sage: K.lattice_dim()**2 - K.dim()*codim(K)
        19

    The Lyapunov rank should be additive on a product of proper cones
    [Rudolf et al.]_::

        sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
        sage: K = L31.cartesian_product(octant)
        sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
        True

    Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
    The cone ``K`` in the following example is isomorphic to the nonnegative
    octant in `\mathbb{R}^{3}`::

        sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
        sage: lyapunov_rank(K)
        3

    The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
    itself [Rudolf et al.]_::

        sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
        True

    TESTS:

    The Lyapunov rank should be additive on a product of proper cones
    [Rudolf et al.]_::

        sage: set_random_seed()
        sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
        sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
        sage: K = K1.cartesian_product(K2)
        sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
        True

    The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
    itself [Rudolf et al.]_::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8)
        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
        True

    Make sure we exercise the non-strictly-convex/non-solid case::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
        True

    Let's check the other permutations as well, just to be sure::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
        True

    ::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
        True

    ::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
        True

    The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
    be any number between `1` and `n` inclusive, excluding `n-1`
    [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
    trivial cone in a trivial space as well. However, in zero dimensions,
    the Lyapunov rank of the trivial cone will be zero::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
        sage: b = lyapunov_rank(K)
        sage: n = K.lattice_dim()
        sage: (n == 0 or 1 <= b) and b <= n
        True
        sage: b == n-1
        False

    In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
    Lyapunov rank `n-1` in `n` dimensions::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8)
        sage: b = lyapunov_rank(K)
        sage: n = K.lattice_dim()
        sage: b == n-1
        False

    The calculation of the Lyapunov rank of an improper cone can be
    reduced to that of a proper cone [Orlitzky/Gowda]_::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8)
        sage: actual = lyapunov_rank(K)
        sage: K_S = restrict_span(K)
        sage: P = restrict_span(K_S.dual()).dual()
        sage: l = lineality(K)
        sage: c = codim(K)
        sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2
        sage: actual == expected
        True

    The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
        sage: lyapunov_rank(K) == len(LL(K))
        True

    In fact the same can be said of any cone. These additional tests
    just increase our confidence that the reduction scheme works::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
        sage: lyapunov_rank(K) == len(LL(K))
        True

    ::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
        sage: lyapunov_rank(K) == len(LL(K))
        True

    ::

        sage: set_random_seed()
        sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
        sage: lyapunov_rank(K) == len(LL(K))
        True

    """
    K_orig = K
    beta = 0

    m = K.dim()
    n = K.lattice_dim()
    l = lineality(K)

    if m < n:
        # K is not solid, project onto its span.
        K = restrict_span(K)

        # Lemma 2
        beta += m*(n - m) + (n - m)**2

    if l > 0:
        # K is not pointed, project its dual onto its span.
        # Uses a proposition from our paper, i.e. this is
        # equivalent to K = restrict_span(K.dual()).dual()
        #K = restrict_span(intersect_span(K,K.dual()), K.dual())
        K = restrict_span(K, K.dual())

        # Lemma 3
        beta += m * l

    beta += len(LL(K))
    return beta