Ensure that we generate min_rays generators; add more cone tests.

[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 *

# TODO: This test fails, maybe due to a bug in the existing cone code.
#     If we request enough generators to span the space, then the returned
#     cone should equal the ambient space::
#
#        sage: K = random_cone(min_dim=5, max_dim=5, min_rays=10, max_rays=10)
#        sage: K.lines().dimension() == K.lattice_dim()
#        True

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 number of generating rays is naturally limited to twice the
    dimension of the ambient space. Take for example $\mathbb{R}^{2}$.
    You could have the generators $\left\{ \pm e_{1}, \pm e_{2}
    \right\}$, with cardinality $4 = 2 \cdot 2$; however any other ray
    in the space is a nonnegative linear combination of those four.

    .. 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.
    #
    # However, since we're going to *check* whether or not we actually
    # have ``r``, we need ``r`` rays to be attainable. So we need to
    # limit ``r`` to twice the dimension of the ambient space.
    #
    r = min(r, 2*d)
    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.
    while r > K.nrays():
        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())