# 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 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:

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

    2. 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`::

        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 `L^{3}_{1}` cone is known to have a Lyapunov rank of one::

        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::

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

    The Lyapunov rank should be additive on a product of cones::

        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. 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::

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

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

    xs = [V(x) for x in K.rays()]
    ss = [V(s) for s in K.dual().rays()]

    # WARNING: This isn't really C(K), it only contains the pairs
    # (x,s) in C(K) where x,s are extreme in their respective cones.
    C_of_K = [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]

    matrices = [x.column() * s.row() 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)

    def phi(m):
        r"""
        Convert a matrix to a vector isomorphically.
        """
        return W(m.list())

    vectors = [phi(m) for m in matrices]

    return (W.dimension() - W.span(vectors).rank())