From 11ebb337a712d805112d43110b1d23c1a4252160 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Sun, 25 Sep 2016 16:17:07 -0400
Subject: [PATCH] Add an is_cross_positive() function and implement
 is_lyapunov_like() using it.

---
 mjo/cone/cone.py | 78 ++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 76 insertions(+), 2 deletions(-)

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index 7e9c549..be05f5e 100644
--- a/mjo/cone/cone.py
+++ b/mjo/cone/cone.py
@@ -1,5 +1,72 @@
 from sage.all import *
 
+def is_cross_positive(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 >= 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.
+
+    INPUT:
+
+    - ``L`` -- A linear transformation or matrix.
+
+    - ``K`` -- A polyhedral closed convex cone.
+
+    OUTPUT:
+
+    ``True`` if it can be proven that ``L`` is cross-positive on ``K``,
+    and ``False`` otherwise.
+
+    .. WARNING::
+
+        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:
+
+    The identity is always cross-positive in a nontrivial 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_cross_positive(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(L,K)
+        True
+
+        Everything in ``K.cross_positive_operator_gens()`` should be
+        cross-positive on ``K``::
+
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
+        sage: all([ is_cross_positive(L,K)
+        ....:       for L in K.cross_positive_operator_gens() ])
+        True
+
+    """
+    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')
+
+
 def is_lyapunov_like(L,K):
     r"""
     Determine whether or not ``L`` is Lyapunov-like on ``K``.
@@ -67,8 +134,15 @@ def is_lyapunov_like(L,K):
         True
 
     """
-    return all([(L*x).inner_product(s) == 0
-                for (x,s) in K.discrete_complementarity_set()])
+    if L.base_ring().is_exact() or L.base_ring() is SR:
+        V = VectorSpace(K.lattice().base_field(), K.lattice_dim()**2)
+        LL_of_K = V.span([ V(m.list()) for m in K.lyapunov_like_basis() ])
+        return V(L.list()) in LL_of_K
+    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()
-- 
2.43.2