Factor out the is_full_space() function.

[sage.d.git] / mjo / cone / cone.py

# 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 is_full_space(K):
    r"""
    Return whether or not this cone is equal to its ambient vector space.

    OUTPUT:

    ``True`` if this cone is the entire vector space and ``False``
    otherwise.

    EXAMPLES:

    A ray in two dimensions is not equal to the entire space::

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

    Neither is the nonnegative orthant::

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

    The right half-space contains a vector subspace, but it is still not
    equal to the entire plane::

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

    But if we include nonnegative sums from both axes, then the resulting
    cone is the entire two-dimensional space::

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

    """
    return K.linear_subspace() == K.lattice().vector_space()


def random_cone(min_dim=0, max_dim=None, min_rays=0, max_rays=None):
    r"""
    Generate a random rational convex polyhedral cone.

    Lower and upper bounds may be provided for both the dimension of the
    ambient space and the number of generating rays of the cone. If a
    lower bound is left unspecified, it defaults to zero. Unspecified
    upper bounds will be chosen randomly.

    The lower bound on the number of rays is limited to twice the
    maximum dimension of the ambient vector space. To see why, consider
    the space $\mathbb{R}^{2}$, and suppose we have generated four rays,
    $\left\{ \pm e_{1}, \pm e_{2} \right\}$. Clearly any other ray in
    the space is a nonnegative linear combination of those four,
    so it is hopeless to generate more. It is therefore an error
    to request more in the form of ``min_rays``.

    .. NOTE:

        If you do not explicitly request more than ``2 * max_dim`` rays,
        a larger number may still be randomly generated. In that case,
        the returned cone will simply be equal to the entire space.

    INPUT:

    - ``min_dim`` (default: zero) -- A nonnegative integer representing the
                                     minimum dimension of the ambient lattice.

    - ``max_dim`` (default: random) -- A nonnegative integer representing
                                       the maximum dimension of the ambient
                                       lattice.

    - ``min_rays`` (default: zero) -- A nonnegative integer representing the
                                      minimum number of generating rays of the
                                      cone.

    - ``max_rays`` (default: random) -- A nonnegative integer representing the
                                        maximum number of generating rays of
                                        the cone.

    OUTPUT:

    A new, randomly generated cone.

    A ``ValueError` will be thrown under the following conditions:

      * Any of ``min_dim``, ``max_dim``, ``min_rays``, or ``max_rays``
        are negative.

      * ``max_dim`` is less than ``min_dim``.

      * ``max_rays`` is less than ``min_rays``.

      * ``min_rays`` is greater than twice ``max_dim``.

    EXAMPLES:

    If we set the lower/upper bounds to zero, then our result is
    predictable::

        sage: random_cone(0,0,0,0)
        0-d cone in 0-d lattice N

    We can predict the dimension when ``min_dim == max_dim``::

        sage: random_cone(min_dim=4, max_dim=4, min_rays=0, max_rays=0)
        0-d cone in 4-d lattice N

    Likewise for the number of rays when ``min_rays == max_rays``::

        sage: random_cone(min_dim=10, max_dim=10, min_rays=10, max_rays=10)
        10-d cone in 10-d lattice N

    TESTS:

    It's hard to test the output of a random process, but we can at
    least make sure that we get a cone back::

        sage: from sage.geometry.cone import is_Cone # long time
        sage: K = random_cone() # long time
        sage: is_Cone(K)        # long time
        True

    The upper/lower bounds are respected::

        sage: K = random_cone(min_dim=5, max_dim=10, min_rays=3, max_rays=4)
        sage: 5 <= K.lattice_dim() and K.lattice_dim() <= 10
        True
        sage: 3 <= K.nrays() and K.nrays() <= 4
        True

    Ensure that an exception is raised when either lower bound is greater
    than its respective upper bound::

        sage: random_cone(min_dim=5, max_dim=2)
        Traceback (most recent call last):
        ...
        ValueError: max_dim cannot be less than min_dim.

        sage: random_cone(min_rays=5, max_rays=2)
        Traceback (most recent call last):
        ...
        ValueError: max_rays cannot be less than min_rays.

    And if we request too many rays::

        sage: random_cone(min_rays=5, max_dim=1)
        Traceback (most recent call last):
        ...
        ValueError: min_rays cannot be larger than twice max_dim.

    """

    # Catch obvious mistakes so that we can generate clear error
    # messages.

    if min_dim < 0:
        raise ValueError('min_dim must be nonnegative.')

    if min_rays < 0:
        raise ValueError('min_rays must be nonnegative.')

    if max_dim is not None:
        if max_dim < 0:
            raise ValueError('max_dim must be nonnegative.')
        if (max_dim < min_dim):
            raise ValueError('max_dim cannot be less than min_dim.')
        if min_rays > 2*max_dim:
            raise ValueError('min_rays cannot be larger than twice max_dim.')

    if max_rays is not None:
        if max_rays < 0:
            raise ValueError('max_rays must be nonnegative.')
        if (max_rays < min_rays):
            raise ValueError('max_rays cannot be less than min_rays.')


    def random_min_max(l,u):
        r"""
        We need to handle two cases for the upper bounds, and we need to do
        the same thing for max_dim/max_rays. So we consolidate the logic here.
        """
        if u is None:
            # The upper bound is unspecified; return a random integer
            # in [l,infinity).
            return l + ZZ.random_element().abs()
        else:
            # We have an upper bound, and it's greater than or equal
            # to our lower bound. So we generate a random integer in
            # [0,u-l], and then add it to l to get something in
            # [l,u]. To understand the "+1", check the
            # ZZ.random_element() docs.
            return l + ZZ.random_element(u - l + 1)

    d = random_min_max(min_dim, max_dim)
    r = random_min_max(min_rays, max_rays)

    L = ToricLattice(d)

    # The rays are trickier to generate, since we could generate v and
    # 2*v as our "two rays." In that case, the resuting cone would
    # have one generating ray. To avoid such a situation, we start by
    # generating ``r`` rays where ``r`` is the number we want to end
    # up with...
    rays = [L.random_element() for i in range(0, r)]

    # (The lattice parameter is required when no rays are given, so we
    # pass it just in case ``r == 0``).
    K = Cone(rays, lattice=L)

    # Now if we generated two of the "same" rays, we'll have fewer
    # generating rays than ``r``. In that case, we keep making up new
    # rays and recreating the cone until we get the right number of
    # independent generators. We can obviously stop if ``K`` is the
    # entire ambient vector space.
    while r > K.nrays() and not is_full_space(K):
        rays.append(L.random_element())
        K = Cone(rays)

    return K


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: K1 = random_cone(max_dim=10, max_rays=10)
        sage: K2 = K1.dual()
        sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
        sage: actual = discrete_complementarity_set(K1)
        sage: actual == 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 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

    TESTS:

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

        sage: K1 = random_cone(max_dim=10, max_rays=10)
        sage: K2 = random_cone(max_dim=10, max_rays=10)
        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::

        sage: K = random_cone(max_dim=10, max_rays=10)
        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
        True

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

    C_of_K = discrete_complementarity_set(K)

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