X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;ds=sidebyside;f=mjo%2Fcone%2Fcone.py;h=aeec0c90b5f582c8ba858d4b616fff191fb6d529;hb=c2d15c0884c8b92483f58826747887bd2bcdcdeb;hp=7e9c549eec66ede6dc0c94bd67e0e16d4d999538;hpb=455255081db7cf2fbb9a221029a19d7d17310577;p=sage.d.git

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index 7e9c549..aeec0c9 100644
--- a/mjo/cone/cone.py
+++ b/mjo/cone/cone.py
@@ -1,74 +1,436 @@
 from sage.all import *
+from sage.geometry.cone import is_Cone
 
-def is_lyapunov_like(L,K):
+def is_positive_on(L,K):
+    r"""
+    Determine whether or not ``L`` is positive on ``K``.
+
+    We say that ``L`` is positive on a closed convex cone ``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``.
+
+    To reliably check whether or not ``L`` is positive, its base ring
+    must be either exact (for example, the rationals) or ``SR``. An
+    exact ring is more reliable, but in some cases a matrix whose
+    entries contain symbolic constants like ``e`` and ``pi`` will work.
+
+    INPUT:
+
+    - ``L`` -- A matrix over either an exact ring or ``SR``.
+
+    - ``K`` -- A polyhedral closed convex cone.
+
+    OUTPUT:
+
+    If the base ring of ``L`` is exact, then ``True`` will be returned if
+    and only if ``L`` is positive on ``K``.
+
+    If the base ring of ``L`` is ``SR``, then the situation is more
+    complicated:
+
+    - ``True`` will be returned if it can be proven that ``L``
+      is positive on ``K``.
+    - ``False`` will be returned if it can be proven that ``L``
+      is not positive on ``K``.
+    - ``False`` will also be returned if we can't decide; specifically
+      if we arrive at a symbolic inequality that cannot be resolved.
+
+    .. SEEALSO::
+
+          :func:`is_cross_positive_on`,
+          :func:`is_Z_operator_on`,
+          :func:`is_lyapunov_like_on`
+
+    EXAMPLES:
+
+    Nonnegative matrices are positive operators on the nonnegative
+    orthant::
+
+        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
+
+    TESTS:
+
+    The identity operator is always positive::
+
+        sage: set_random_seed()
+        sage: K = random_cone(max_ambient_dim=8)
+        sage: L = identity_matrix(K.lattice_dim())
+        sage: is_positive_on(L,K)
+        True
+
+    The "zero" operator is always positive::
+
+        sage: K = random_cone(max_ambient_dim=8)
+        sage: R = K.lattice().vector_space().base_ring()
+        sage: L = zero_matrix(R, K.lattice_dim())
+        sage: is_positive_on(L,K)
+        True
+
+    Everything in ``K.positive_operators_gens()`` should be
+    positive on ``K``::
+
+        sage: K = random_cone(max_ambient_dim=5)
+        sage: all([ is_positive_on(L,K)                     # long time
+        ....:       for L in K.positive_operators_gens() ]) # long time
+        True
+        sage: all([ is_positive_on(L.change_ring(SR),K)     # long time
+        ....:       for L in K.positive_operators_gens() ]) # long time
+        True
+
+    Technically we could test this, but for now only closed convex cones
+    are supported as our ``K`` argument::
+
+        sage: K = [ vector([1,2,3]), vector([5,-1,7]) ]
+        sage: L = identity_matrix(3)
+        sage: is_positive_on(L,K)
+        Traceback (most recent call last):
+        ...
+        TypeError: K must be a Cone.
+
+    We can't give reliable answers over inexact rings::
+
+        sage: K = Cone([(1,2,3), (4,5,6)])
+        sage: L = identity_matrix(RR,3)
+        sage: is_positive_on(L,K)
+        Traceback (most recent call last):
+        ...
+        ValueError: The base ring of L is neither SR nor exact.
+
+    """
+
+    if not is_Cone(K):
+        raise TypeError('K must be a Cone.')
+    if not L.base_ring().is_exact() and not L.base_ring() is SR:
+        raise ValueError('The base ring of L is neither SR nor exact.')
+
+    if L.base_ring().is_exact():
+        # This should be way faster than computing the dual and
+        # checking a bunch of inequalities, but it doesn't work if
+        # ``L*x`` is symbolic. For example, ``e in Cone([(1,)])``
+        # is true, but returns ``False``.
+        return all([ L*x in K for x in K ])
+    else:
+        # 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() ])
+
+
+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 a closed convex cone``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.
+
+    To reliably check whether or not ``L`` is cross-positive, its base
+    ring must be either exact (for example, the rationals) or ``SR``. An
+    exact ring is more reliable, but in some cases a matrix whose
+    entries contain symbolic constants like ``e`` and ``pi`` will work.
+
+    INPUT:
+
+    - ``L`` -- A matrix over either an exact ring or ``SR``.
+
+    - ``K`` -- A polyhedral closed convex cone.
+
+    OUTPUT:
+
+    If the base ring of ``L`` is exact, then ``True`` will be returned if
+    and only if ``L`` is cross-positive on ``K``.
+
+    If the base ring of ``L`` is ``SR``, then the situation is more
+    complicated:
+
+    - ``True`` will be returned if it can be proven that ``L``
+      is cross-positive on ``K``.
+    - ``False`` will be returned if it can be proven that ``L``
+      is not cross-positive on ``K``.
+    - ``False`` will also be returned if we can't decide; specifically
+      if we arrive at a symbolic inequality that cannot be resolved.
+
+    .. SEEALSO::
+
+          :func:`is_positive_on`,
+          :func:`is_Z_operator_on`,
+          :func:`is_lyapunov_like_on`
+
+    EXAMPLES:
+
+    The identity operator is always cross-positive::
+
+        sage: set_random_seed()
+        sage: K = random_cone(max_ambient_dim=8)
+        sage: L = identity_matrix(K.lattice_dim())
+        sage: is_cross_positive_on(L,K)
+        True
+
+    The "zero" operator is always cross-positive::
+
+        sage: K = random_cone(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:
+
+    Everything in ``K.cross_positive_operators_gens()`` should be
+    cross-positive on ``K``::
+
+        sage: K = random_cone(max_ambient_dim=5)
+        sage: all([ is_cross_positive_on(L,K)                     # long time
+        ....:       for L in K.cross_positive_operators_gens() ]) # long time
+        True
+        sage: all([ is_cross_positive_on(L.change_ring(SR),K)     # long time
+        ....:       for L in K.cross_positive_operators_gens() ]) # long time
+        True
+
+    Technically we could test this, but for now only closed convex cones
+    are supported as our ``K`` argument::
+
+        sage: L = identity_matrix(3)
+        sage: K = [ vector([8,2,-8]), vector([5,-5,7]) ]
+        sage: is_cross_positive_on(L,K)
+        Traceback (most recent call last):
+        ...
+        TypeError: K must be a Cone.
+
+    We can't give reliable answers over inexact rings::
+
+        sage: K = Cone([(1,2,3), (4,5,6)])
+        sage: L = identity_matrix(RR,3)
+        sage: is_cross_positive_on(L,K)
+        Traceback (most recent call last):
+        ...
+        ValueError: The base ring of L is neither SR nor exact.
+
+    """
+    if not is_Cone(K):
+        raise TypeError('K must be a Cone.')
+    if not L.base_ring().is_exact() and not L.base_ring() is SR:
+        raise ValueError('The base ring of L is neither SR nor exact.')
+
+    return all([ s*(L*x) >= 0
+                 for (x,s) in K.discrete_complementarity_set() ])
+
+def is_Z_operator_on(L,K):
+    r"""
+    Determine whether or not ``L`` is a Z-operator on ``K``.
+
+    We say that ``L`` is a Z-operator on a closed convex cone``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 matrix is a Z-operator on ``K`` if and only if its negation is a
+    cross-positive operator on ``K``.
+
+    To reliably check whether or not ``L`` is a Z operator, its base
+    ring must be either exact (for example, the rationals) or ``SR``. An
+    exact ring is more reliable, but in some cases a matrix whose
+    entries contain symbolic constants like ``e`` and ``pi`` will work.
+
+    INPUT:
+
+    - ``L`` -- A matrix over either an exact ring or ``SR``.
+
+    - ``K`` -- A polyhedral closed convex cone.
+
+    OUTPUT:
+
+    If the base ring of ``L`` is exact, then ``True`` will be returned if
+    and only if ``L`` is a Z-operator on ``K``.
+
+    If the base ring of ``L`` is ``SR``, then the situation is more
+    complicated:
+
+    - ``True`` will be returned if it can be proven that ``L``
+      is a Z-operator on ``K``.
+    - ``False`` will be returned if it can be proven that ``L``
+      is not a Z-operator on ``K``.
+    - ``False`` will also be returned if we can't decide; specifically
+      if we arrive at a symbolic inequality that cannot be resolved.
+
+    .. SEEALSO::
+
+          :func:`is_positive_on`,
+          :func:`is_cross_positive_on`,
+          :func:`is_lyapunov_like_on`
+
+    EXAMPLES:
+
+    The identity operator is always a Z-operator::
+
+        sage: set_random_seed()
+        sage: K = random_cone(max_ambient_dim=8)
+        sage: L = identity_matrix(K.lattice_dim())
+        sage: is_Z_operator_on(L,K)
+        True
+
+    The "zero" operator is always a Z-operator::
+
+        sage: K = random_cone(max_ambient_dim=8)
+        sage: R = K.lattice().vector_space().base_ring()
+        sage: L = zero_matrix(R, K.lattice_dim())
+        sage: is_Z_operator_on(L,K)
+        True
+
+    TESTS:
+
+    Everything in ``K.Z_operators_gens()`` should be a Z-operator
+    on ``K``::
+
+        sage: K = random_cone(max_ambient_dim=5)
+        sage: all([ is_Z_operator_on(L,K)            # long time
+        ....:       for L in K.Z_operators_gens() ]) # long time
+        True
+        sage: all([ is_Z_operator_on(L.change_ring(SR),K) # long time
+        ....:       for L in K.Z_operators_gens() ])      # long time
+        True
+
+    Technically we could test this, but for now only closed convex cones
+    are supported as our ``K`` argument::
+
+        sage: L = identity_matrix(3)
+        sage: K = [ vector([-4,20,3]), vector([1,-5,2]) ]
+        sage: is_Z_operator_on(L,K)
+        Traceback (most recent call last):
+        ...
+        TypeError: K must be a Cone.
+
+
+    We can't give reliable answers over inexact rings::
+
+        sage: K = Cone([(1,2,3), (4,5,6)])
+        sage: L = identity_matrix(RR,3)
+        sage: is_Z_operator_on(L,K)
+        Traceback (most recent call last):
+        ...
+        ValueError: The base ring of L is neither SR nor exact.
+
+    """
+    return is_cross_positive_on(-L,K)
+
+
+def is_lyapunov_like_on(L,K):
     r"""
     Determine whether or not ``L`` is Lyapunov-like 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 Lyapunov-like on a closed convex cone ``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.
 
-    There are faster ways of checking this property. For example, we
-    could compute a `lyapunov_like_basis` of the cone, and then test
-    whether or not the given matrix is contained in the span of that
-    basis. The value of this function is that it works on symbolic
-    matrices.
+    An operator is Lyapunov-like on ``K`` if and only if both the
+    operator itself and its negation are cross-positive on ``K``.
+
+    To reliably check whether or not ``L`` is Lyapunov-like, its base
+    ring must be either exact (for example, the rationals) or ``SR``. An
+    exact ring is more reliable, but in some cases a matrix whose
+    entries contain symbolic constants like ``e`` and ``pi`` will work.
 
     INPUT:
 
-    - ``L`` -- A linear transformation or matrix.
+    - ``L`` -- A matrix over either an exact ring or ``SR``.
 
     - ``K`` -- A polyhedral closed convex cone.
 
     OUTPUT:
 
-    ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
-    and ``False`` otherwise.
+    If the base ring of ``L`` is exact, then ``True`` will be returned if
+    and only if ``L`` is Lyapunov-like on ``K``.
 
-    .. WARNING::
+    If the base ring of ``L`` is ``SR``, then the situation is more
+    complicated:
 
-        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.
+    - ``True`` will be returned if it can be proven that ``L``
+      is Lyapunov-like on ``K``.
+    - ``False`` will be returned if it can be proven that ``L``
+      is not Lyapunov-like on ``K``.
+    - ``False`` will also be returned if we can't decide; specifically
+      if we arrive at a symbolic inequality that cannot be resolved.
 
-    REFERENCES:
+    .. SEEALSO::
 
-    M. Orlitzky. The Lyapunov rank of an improper cone.
-    http://www.optimization-online.org/DB_HTML/2015/10/5135.html
+          :func:`is_positive_on`,
+          :func:`is_cross_positive_on`,
+          :func:`is_Z_operator_on`
 
     EXAMPLES:
 
-    The identity is always Lyapunov-like in a nontrivial space::
+    Diagonal matrices are Lyapunov-like operators on the nonnegative
+    orthant::
+
+        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:
+
+    The identity operator is always Lyapunov-like::
 
         sage: set_random_seed()
-        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: K = random_cone(max_ambient_dim=8)
         sage: L = identity_matrix(K.lattice_dim())
-        sage: is_lyapunov_like(L,K)
+        sage: is_lyapunov_like_on(L,K)
         True
 
-    As is the "zero" transformation::
+    The "zero" operator is always Lyapunov-like::
 
-        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: K = random_cone(max_ambient_dim=8)
         sage: R = K.lattice().vector_space().base_ring()
         sage: L = zero_matrix(R, K.lattice_dim())
-        sage: is_lyapunov_like(L,K)
+        sage: is_lyapunov_like_on(L,K)
         True
 
-        Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
-        on ``K``::
+    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(L,K) for L in K.lyapunov_like_basis() ])
+        sage: K = random_cone(max_ambient_dim=5)
+        sage: all([ is_lyapunov_like_on(L,K)            # long time
+        ....:       for L in K.lyapunov_like_basis() ]) # long time
+        True
+        sage: all([ is_lyapunov_like_on(L.change_ring(SR),K) # long time
+        ....:       for L in K.lyapunov_like_basis() ])      # long time
         True
 
+    Technically we could test this, but for now only closed convex cones
+    are supported as our ``K`` argument::
+
+        sage: L = identity_matrix(3)
+        sage: K = [ vector([2,2,-1]), vector([5,4,-3]) ]
+        sage: is_lyapunov_like_on(L,K)
+        Traceback (most recent call last):
+        ...
+        TypeError: K must be a Cone.
+
+    We can't give reliable answers over inexact rings::
+
+        sage: K = Cone([(1,2,3), (4,5,6)])
+        sage: L = identity_matrix(RR,3)
+        sage: is_lyapunov_like_on(L,K)
+        Traceback (most recent call last):
+        ...
+        ValueError: The base ring of L is neither SR nor exact.
+
     """
-    return all([(L*x).inner_product(s) == 0
-                for (x,s) in K.discrete_complementarity_set()])
+    if not is_Cone(K):
+        raise TypeError('K must be a Cone.')
+    if not L.base_ring().is_exact() and not L.base_ring() is SR:
+        raise ValueError('The base ring of L is neither SR nor exact.')
+
+    return all([ s*(L*x) == 0
+                 for (x,s) in K.discrete_complementarity_set() ])
+
 
 def LL_cone(K):
     gens = K.lyapunov_like_basis()
@@ -76,18 +438,18 @@ def LL_cone(K):
     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 
 def Sigma_cone(K):
-    gens = K.cross_positive_operator_gens()
+    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_operator_gens()
+    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_operator_gens(K2)
+    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)