Fix the is_positive_on test and give better examples.

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

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index 49df3e9bc73b5071b4120f3ebfd41acb69921d01..01ef5f9d5cc56b087344237b9b5c915423ac0921 100644 (file)
--- a/mjo/cone/cone.py
+++ b/mjo/cone/cone.py
@@ -1,14 +1,12 @@
 from sage.all import *
 
 from sage.all import *
 

-def is_lyapunov_like(L,K):
+def is_positive_on(L,K):

     r"""
     r"""

-    Determine whether or not ``L`` is Lyapunov-like on ``K``.
+    Determine whether or not ``L`` is positive on ``K``.

 
 

-    We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
-    L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
-    `\left\langle x,s \right\rangle` in the complementarity set of
-    ``K``. It is known [Orlitzky]_ that this property need only be
-    checked for generators of ``K`` and its dual.
+    We say that ``L`` is positive on ``K`` if `L\left\lparen x
+    \right\rparen` belongs to ``K`` for all `x` in ``K``. This
+    property need only be checked for generators of ``K``.

 
     INPUT:
 
 
     INPUT:
 


@@ -18,296 +16,311 @@ def is_lyapunov_like(L,K):
 
     OUTPUT:
 
 
     OUTPUT:
 

-    ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
+    ``True`` if it can be proven that ``L`` is positive on ``K``,

     and ``False`` otherwise.
 
     .. WARNING::
 
     and ``False`` otherwise.
 
     .. WARNING::
 

-        If this function returns ``True``, then ``L`` is Lyapunov-like
+        If this function returns ``True``, then ``L`` is positive

         on ``K``. However, if ``False`` is returned, that could mean one
         of two things. The first is that ``L`` is definitely not
         on ``K``. However, if ``False`` is returned, that could mean one
         of two things. The first is that ``L`` is definitely not

-        Lyapunov-like on ``K``. The second is more of an "I don't know"
+        positive on ``K``. The second is more of an "I don't know"

         answer, returned (for example) if we cannot prove that an inner
         answer, returned (for example) if we cannot prove that an inner

-        product is zero.
+        product is nonnegative.

 
 

-    REFERENCES:
+    EXAMPLES:

 
 

-    M. Orlitzky. The Lyapunov rank of an improper cone.
-    http://www.optimization-online.org/DB_HTML/2015/10/5135.html
+    Positive operators on the nonnegative orthant are nonnegative
+    matrices::

 
 

-    EXAMPLES:
+        sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
+        sage: L = random_matrix(QQ,3).apply_map(abs)
+        sage: is_positive_on(L,K)
+        True

 
 

-    The identity is always Lyapunov-like in a nontrivial space::
+    TESTS:
+
+    The identity is always positive in a nontrivial space::

 
         sage: set_random_seed()
 
         sage: set_random_seed()

-        sage: K = random_cone(min_ambient_dim = 1, max_rays = 8)
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)

         sage: L = identity_matrix(K.lattice_dim())
         sage: L = identity_matrix(K.lattice_dim())

-        sage: is_lyapunov_like(L,K)
+        sage: is_positive_on(L,K)

         True
 
     As is the "zero" transformation::
 
         True
 
     As is the "zero" transformation::
 

-        sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)

         sage: R = K.lattice().vector_space().base_ring()
         sage: L = zero_matrix(R, K.lattice_dim())
         sage: R = K.lattice().vector_space().base_ring()
         sage: L = zero_matrix(R, K.lattice_dim())

-        sage: is_lyapunov_like(L,K)
+        sage: is_positive_on(L,K)

         True
 
         True
 

-        Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
-        on ``K``::
+    Everything in ``K.positive_operators_gens()`` should be
+    positive on ``K``::

 
 

-        sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
-        sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
+        sage: all([ is_positive_on(L,K)
+        ....:       for L in K.positive_operators_gens() ])
+        True
+        sage: all([ is_positive_on(L.change_ring(SR),K)
+        ....:       for L in K.positive_operators_gens() ])

         True
 
     """
         True
 
     """

-    return all([(L*x).inner_product(s) == 0
-                for (x,s) in K.discrete_complementarity_set()])
+    if L.base_ring().is_exact():
+        # This could potentially be extended to other types of ``K``...
+        return all([ L*x in K for x in K ])
+    elif L.base_ring() is SR:
+        # Fall back to inequality-checking when the entries of ``L``
+        # might be symbolic.
+        return all([ s*(L*x) >= 0 for x in K for s in K.dual() ])
+    else:
+        # The only inexact ring that we're willing to work with is SR,
+        # since it can still be exact when working with symbolic
+        # constants like pi and e.
+        raise ValueError('base ring of operator L is neither SR nor exact')
+
+
+def is_cross_positive_on(L,K):
+    r"""
+    Determine whether or not ``L`` is cross-positive on ``K``.

 
 

+    We say that ``L`` is cross-positive on ``K`` if `\left\langle
+    L\left\lparenx\right\rparen,s\right\rangle \ge 0` for all pairs
+    `\left\langle x,s \right\rangle` in the complementarity set of
+    ``K``. This property need only be checked for generators of
+    ``K`` and its dual.

 
 

-def random_element(K):
-    r"""
-    Return a random element of ``K`` from its ambient vector space.
+    INPUT:
+
+    - ``L`` -- A linear transformation or matrix.
+
+    - ``K`` -- A polyhedral closed convex cone.
+
+    OUTPUT:

 
 

-    ALGORITHM:
+    ``True`` if it can be proven that ``L`` is cross-positive on ``K``,
+    and ``False`` otherwise.

 
 

-    The cone ``K`` is specified in terms of its generators, so that
-    ``K`` is equal to the convex conic combination of those generators.
-    To choose a random element of ``K``, we assign random nonnegative
-    coefficients to each generator of ``K`` and construct a new vector
-    from the scaled rays.
+    .. WARNING::

 
 

-    A vector, rather than a ray, is returned so that the element may
-    have non-integer coordinates. Thus the element may have an
-    arbitrarily small norm.
+        If this function returns ``True``, then ``L`` is cross-positive
+        on ``K``. However, if ``False`` is returned, that could mean one
+        of two things. The first is that ``L`` is definitely not
+        cross-positive on ``K``. The second is more of an "I don't know"
+        answer, returned (for example) if we cannot prove that an inner
+        product is nonnegative.

 
     EXAMPLES:
 
 
     EXAMPLES:
 

-    A random element of the trivial cone is zero::
+    The identity is always cross-positive in a nontrivial space::

 
         sage: set_random_seed()
 
         sage: set_random_seed()

-        sage: K = Cone([], ToricLattice(0))
-        sage: random_element(K)
-        ()
-        sage: K = Cone([(0,)])
-        sage: random_element(K)
-        (0)
-        sage: K = Cone([(0,0)])
-        sage: random_element(K)
-        (0, 0)
-        sage: K = Cone([(0,0,0)])
-        sage: random_element(K)
-        (0, 0, 0)
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: L = identity_matrix(K.lattice_dim())
+        sage: is_cross_positive_on(L,K)
+        True
+
+    As is the "zero" transformation::
+
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: R = K.lattice().vector_space().base_ring()
+        sage: L = zero_matrix(R, K.lattice_dim())
+        sage: is_cross_positive_on(L,K)
+        True

 
     TESTS:
 
 
     TESTS:
 

-    Any cone should contain an element of itself::
+    Everything in ``K.cross_positive_operators_gens()`` should be
+    cross-positive on ``K``::

 
 

-        sage: set_random_seed()
-        sage: K = random_cone(max_rays = 8)
-        sage: K.contains(random_element(K))
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
+        sage: all([ is_cross_positive_on(L,K)
+        ....:       for L in K.cross_positive_operators_gens() ])
+        True
+        sage: all([ is_cross_positive_on(L.change_ring(SR),K)
+        ....:       for L in K.cross_positive_operators_gens() ])

         True
 
     """
         True
 
     """

-    V = K.lattice().vector_space()
-    F = V.base_ring()
-    coefficients = [ F.random_element().abs() for i in range(K.nrays()) ]
-    vector_gens  = map(V, K.rays())
-    scaled_gens  = [ coefficients[i]*vector_gens[i]
-                         for i in range(len(vector_gens)) ]
+    if L.base_ring().is_exact() or L.base_ring() is SR:
+        return all([ s*(L*x) >= 0
+                     for (x,s) in K.discrete_complementarity_set() ])
+    else:
+        # The only inexact ring that we're willing to work with is SR,
+        # since it can still be exact when working with symbolic
+        # constants like pi and e.
+        raise ValueError('base ring of operator L is neither SR nor exact')

 
 

-    # Make sure we return a vector. Without the coercion, we might
-    # return ``0`` when ``K`` has no rays.
-    v = V(sum(scaled_gens))
-    return v

 
 

-
-def positive_operators(K):
+def is_Z_on(L,K):

     r"""
     r"""

-    Compute generators of the cone of positive operators on this cone.
+    Determine whether or not ``L`` is a Z-operator on ``K``.

 
 

-    OUTPUT:
+    We say that ``L`` is a Z-operator on ``K`` if `\left\langle
+    L\left\lparenx\right\rparen,s\right\rangle \le 0` for all pairs
+    `\left\langle x,s \right\rangle` in the complementarity set of
+    ``K``. It is known that this property need only be
+    checked for generators of ``K`` and its dual.

 
 

-    A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
-    Each matrix ``P`` in the list should have the property that ``P*x``
-    is an element of ``K`` whenever ``x`` is an element of
-    ``K``. Moreover, any nonnegative linear combination of these
-    matrices shares the same property.
+    A matrix is a Z-operator on ``K`` if and only if its negation is a
+    cross-positive operator on ``K``.

 
 

-    EXAMPLES:
+    INPUT:
+
+    - ``L`` -- A linear transformation or matrix.
+
+    - ``K`` -- A polyhedral closed convex cone.

 
 

-    The trivial cone in a trivial space has no positive operators::
+    OUTPUT:

 
 

-        sage: K = Cone([], ToricLattice(0))
-        sage: positive_operators(K)
-        []
+    ``True`` if it can be proven that ``L`` is a Z-operator on ``K``,
+    and ``False`` otherwise.

 
 

-    Positive operators on the nonnegative orthant are nonnegative matrices::
+    .. WARNING::

 
 

-        sage: K = Cone([(1,)])
-        sage: positive_operators(K)
-        [[1]]
+        If this function returns ``True``, then ``L`` is a Z-operator
+        on ``K``. However, if ``False`` is returned, that could mean one
+        of two things. The first is that ``L`` is definitely not
+        a Z-operator on ``K``. The second is more of an "I don't know"
+        answer, returned (for example) if we cannot prove that an inner
+        product is nonnegative.

 
 

-        sage: K = Cone([(1,0),(0,1)])
-        sage: positive_operators(K)
-        [
-        [1 0]  [0 1]  [0 0]  [0 0]
-        [0 0], [0 0], [1 0], [0 1]
-        ]
+    EXAMPLES:

 
 

-    Every operator is positive on the ambient vector space::
+    The identity is always a Z-operator in a nontrivial space::

 
 

-        sage: K = Cone([(1,),(-1,)])
-        sage: K.is_full_space()
+        sage: set_random_seed()
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: L = identity_matrix(K.lattice_dim())
+        sage: is_Z_on(L,K)

         True
         True

-        sage: positive_operators(K)
-        [[1], [-1]]

 
 

-        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
-        sage: K.is_full_space()
+    As is the "zero" transformation::
+
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: R = K.lattice().vector_space().base_ring()
+        sage: L = zero_matrix(R, K.lattice_dim())
+        sage: is_Z_on(L,K)

         True
         True

-        sage: positive_operators(K)
-        [
-        [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
-        [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
-        ]

 
     TESTS:
 
 
     TESTS:
 

-    A positive operator on a cone should send its generators into the cone::
+    Everything in ``K.Z_operators_gens()`` should be a Z-operator
+    on ``K``::

 
 

-        sage: K = random_cone(max_ambient_dim = 6)
-        sage: pi_of_K = positive_operators(K)
-        sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
+        sage: all([ is_Z_on(L,K)
+        ....:       for L in K.Z_operators_gens() ])
+        True
+        sage: all([ is_Z_on(L.change_ring(SR),K)
+        ....:       for L in K.Z_operators_gens() ])

         True
 
     """
         True
 
     """

-    # 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.
-    V = K.lattice().vector_space()
-    W = VectorSpace(V.base_ring(), V.dimension()**2)
-
-    tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
+    return is_cross_positive_on(-L,K)

 
 

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

 
 

-    # Create the *dual* cone of the positive operators, expressed as
-    # long vectors..
-    L = ToricLattice(W.dimension())
-    pi_dual = Cone(vectors, lattice=L)
-
-    # Now compute the desired cone from its dual...
-    pi_cone = pi_dual.dual()
+def is_lyapunov_like_on(L,K):
+    r"""
+    Determine whether or not ``L`` is Lyapunov-like on ``K``.

 
 

-    # And finally convert its rays back to matrix representations.
-    M = MatrixSpace(V.base_ring(), V.dimension())
+    We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
+    L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
+    `\left\langle x,s \right\rangle` in the complementarity set of
+    ``K``. This property need only be checked for generators of
+    ``K`` and its dual.

 
 

-    return [ M(v.list()) for v in pi_cone.rays() ]
+    INPUT:

 
 

+    - ``L`` -- A linear transformation or matrix.

 
 

-def Z_transformations(K):
-    r"""
-    Compute generators of the cone of Z-transformations on this cone.
+    - ``K`` -- A polyhedral closed convex cone.

 
     OUTPUT:
 
 
     OUTPUT:
 

-    A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
-    Each matrix ``L`` in the list should have the property that
-    ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
-    discrete complementarity set of ``K``. Moreover, any nonnegative
-    linear combination of these matrices shares the same property.
-
-    EXAMPLES:
+    ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
+    and ``False`` otherwise.

 
 

-    Z-transformations on the nonnegative orthant are just Z-matrices.
-    That is, matrices whose off-diagonal elements are nonnegative::
-
-        sage: K = Cone([(1,0),(0,1)])
-        sage: Z_transformations(K)
-        [
-        [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
-        [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
-        ]
-        sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
-        sage: all([ z[i][j] <= 0 for z in Z_transformations(K)
-        ....:                    for i in range(z.nrows())
-        ....:                    for j in range(z.ncols())
-        ....:                    if i != j ])
-        True
+    .. WARNING::

 
 

-    The trivial cone in a trivial space has no Z-transformations::
+        If this function returns ``True``, then ``L`` is Lyapunov-like
+        on ``K``. However, if ``False`` is returned, that could mean one
+        of two things. The first is that ``L`` is definitely not
+        Lyapunov-like on ``K``. The second is more of an "I don't know"
+        answer, returned (for example) if we cannot prove that an inner
+        product is zero.

 
 

-        sage: K = Cone([], ToricLattice(0))
-        sage: Z_transformations(K)
-        []
+    EXAMPLES:

 
 

-    Z-transformations on a subspace are Lyapunov-like and vice-versa::
+    Lyapunov-like operators on the nonnegative orthant are diagonal
+    matrices::

 
 

-        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
-        sage: K.is_full_space()
-        True
-        sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
-        sage: zs  = span([ vector(z.list()) for z in Z_transformations(K) ])
-        sage: zs == lls
+        sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
+        sage: L = diagonal_matrix(random_vector(QQ,3))
+        sage: is_lyapunov_like_on(L,K)

         True
 
     TESTS:
 
         True
 
     TESTS:
 

-    The Z-property is possessed by every Z-transformation::
+    The identity is always Lyapunov-like in a nontrivial space::

 
         sage: set_random_seed()
 
         sage: set_random_seed()

-        sage: K = random_cone(max_ambient_dim = 6)
-        sage: Z_of_K = Z_transformations(K)
-        sage: dcs = K.discrete_complementarity_set()
-        sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
-        ....:                                  for (x,s) in dcs])
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: L = identity_matrix(K.lattice_dim())
+        sage: is_lyapunov_like_on(L,K)

         True
 
         True
 

-    The lineality space of Z is LL::
+    As is the "zero" transformation::

 
 

-        sage: set_random_seed()
-        sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6)
-        sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
-        sage: z_cone  = Cone([ z.list() for z in Z_transformations(K) ])
-        sage: z_cone.linear_subspace() == lls
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: R = K.lattice().vector_space().base_ring()
+        sage: L = zero_matrix(R, K.lattice_dim())
+        sage: is_lyapunov_like_on(L,K)
+        True
+
+    Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
+    on ``K``::
+
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
+        sage: all([ is_lyapunov_like_on(L,K)
+        ....:       for L in K.lyapunov_like_basis() ])
+        True
+        sage: all([ is_lyapunov_like_on(L.change_ring(SR),K)
+        ....:       for L in K.lyapunov_like_basis() ])

         True
 
     """
         True
 
     """

-    # 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.
-    V = K.lattice().vector_space()
-    W = VectorSpace(V.base_ring(), V.dimension()**2)
-
-    C_of_K = K.discrete_complementarity_set()
-    tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
-
-    # Turn our matrices into long vectors...
-    vectors = [ W(m.list()) for m in tensor_products ]
-
-    # Create the *dual* cone of the cross-positive operators,
-    # expressed as long vectors..
-    L = ToricLattice(W.dimension())
-    Sigma_dual = Cone(vectors, lattice=L)
-
-    # Now compute the desired cone from its dual...
-    Sigma_cone = Sigma_dual.dual()
-
-    # And finally convert its rays back to matrix representations.
-    # But first, make them negative, so we get Z-transformations and
-    # not cross-positive ones.
-    M = MatrixSpace(V.base_ring(), V.dimension())
-
-    return [ -M(v.list()) for v in Sigma_cone.rays() ]
+    if L.base_ring().is_exact() or L.base_ring() is SR:
+        # The "fast method" of creating a vector space based on a
+        # ``lyapunov_like_basis`` is actually slower than this.
+        return all([ s*(L*x) == 0
+                     for (x,s) in K.discrete_complementarity_set() ])
+    else:
+        # The only inexact ring that we're willing to work with is SR,
+        # since it can still be exact when working with symbolic
+        # constants like pi and e.
+        raise ValueError('base ring of operator L is neither SR nor exact')
+
+def LL_cone(K):
+    gens = K.lyapunov_like_basis()
+    L = ToricLattice(K.lattice_dim()**2)
+    return Cone([ g.list() for g in gens ], lattice=L, check=False)
+
+def Sigma_cone(K):
+    gens = K.cross_positive_operators_gens()
+    L = ToricLattice(K.lattice_dim()**2)
+    return Cone([ g.list() for g in gens ], lattice=L, check=False)
+
+def Z_cone(K):
+    gens = K.Z_operators_gens()
+    L = ToricLattice(K.lattice_dim()**2)
+    return Cone([ g.list() for g in gens ], lattice=L, check=False)
+
+def pi_cone(K1, K2=None):
+    if K2 is None:
+        K2 = K1
+    gens = K1.positive_operators_gens(K2)
+    L = ToricLattice(K1.lattice_dim()*K2.lattice_dim())
+    return Cone([ g.list() for g in gens ], lattice=L, check=False)