Begin working on a two-cone pi(K1,K2).

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index 0bfc3b6a779147b0f0a30799bfe2c149cb818169..d364a01834a8e1c37e00f8468f5ed7136b7e8408 100644 (file)
--- a/mjo/cone/cone.py
+++ b/mjo/cone/cone.py
@@ -71,17 +71,26 @@ def is_lyapunov_like(L,K):
                 for (x,s) in K.discrete_complementarity_set()])
 
 
-def positive_operator_gens(K):
+def positive_operator_gens(K1, K2 = None):
     r"""
-    Compute generators of the cone of positive operators on this cone.
+    Compute generators of the cone of positive operators on this cone. A
+    linear operator on a cone is positive if the image of the cone under
+    the operator is a subset of the cone. This concept can be extended
+    to two cones, where the image of the first cone under a positive
+    operator is a subset of the second cone.
+
+    INPUT:
+
+    - ``K2`` -- (default: ``K1``) the codomain cone; the image of this
+                cone under the returned operators is a subset of ``K2``.
 
     OUTPUT:
 
-    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 list of `m`-by-``n`` matrices where ``m == K2.lattice_dim()`` and
+    ``n == K1.lattice_dim()``. Each matrix ``P`` in the list should have
+    the property that ``P*x`` is an element of ``K2`` whenever ``x`` is
+    an element of ``K1``. Moreover, any nonnegative linear combination of
+    these matrices shares the same property.
 
     REFERENCES:
 
@@ -159,50 +168,58 @@ def positive_operator_gens(K):
 
     TESTS:
 
-    Each positive operator generator should send the generators of the
-    cone into the cone::
+    Each positive operator generator should send the generators of one
+    cone into the other cone::
 
         sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=4)
-        sage: pi_of_K = positive_operator_gens(K)
-        sage: all([ K.contains(P*x) for P in pi_of_K for x in K ])
+        sage: K1 = random_cone(max_ambient_dim=4)
+        sage: K2 = random_cone(max_ambient_dim=4)
+        sage: pi_K1_K2 = positive_operator_gens(K1,K2)
+        sage: all([ K2.contains(P*x) for P in pi_K1_K2 for x in K1 ])
         True
 
-    Each positive operator generator should send a random element of the
-    cone into the cone::
+    Each positive operator generator should send a random element of one
+    cone into the other cone::
 
         sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=4)
-        sage: pi_of_K = positive_operator_gens(K)
-        sage: all([ K.contains(P*K.random_element(QQ)) for P in pi_of_K ])
+        sage: K1 = random_cone(max_ambient_dim=4)
+        sage: K2 = random_cone(max_ambient_dim=4)
+        sage: pi_K1_K2 = positive_operator_gens(K1,K2)
+        sage: all([ K2.contains(P*K1.random_element(QQ)) for P in pi_K1_K2 ])
         True
 
     A random element of the positive operator cone should send the
-    generators of the cone into the cone::
+    generators of one cone into the other cone::
 
         sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=4)
-        sage: pi_of_K = positive_operator_gens(K)
-        sage: L = ToricLattice(K.lattice_dim()**2)
-        sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
+        sage: K1 = random_cone(max_ambient_dim=4)
+        sage: K2 = random_cone(max_ambient_dim=4)
+        sage: pi_K1_K2 = positive_operator_gens(K1,K2)
+        sage: L = ToricLattice(K1.lattice_dim() * K2.lattice_dim())
+        sage: pi_cone = Cone([ g.list() for g in pi_K1_K2 ],
         ....:                lattice=L,
         ....:                check=False)
-        sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
-        sage: all([ K.contains(P*x) for x in K ])
+        sage: P = matrix(K2.lattice_dim(),
+        ....:            K1.lattice_dim(),
+        ....:            pi_cone.random_element(QQ).list())
+        sage: all([ K2.contains(P*x) for x in K1 ])
         True
 
     A random element of the positive operator cone should send a random
-    element of the cone into the cone::
+    element of one cone into the other cone::
 
         sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=4)
-        sage: pi_of_K = positive_operator_gens(K)
-        sage: L = ToricLattice(K.lattice_dim()**2)
-        sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
+        sage: K1 = random_cone(max_ambient_dim=4)
+        sage: K2 = random_cone(max_ambient_dim=4)
+        sage: pi_K1_K2 = positive_operator_gens(K1,K2)
+        sage: L = ToricLattice(K1.lattice_dim() * K2.lattice_dim())
+        sage: pi_cone = Cone([ g.list() for g in pi_K1_K2 ],
         ....:                lattice=L,
         ....:                check=False)
-        sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
-        sage: K.contains(P*K.random_element(ring=QQ))
+        sage: P = matrix(K2.lattice_dim(),
+        ....:            K1.lattice_dim(),
+        ....:            pi_cone.random_element(QQ).list())
+        sage: K2.contains(P*K1.random_element(ring=QQ))
         True
 
     The lineality space of the dual of the cone of positive operators
@@ -396,21 +413,45 @@ def positive_operator_gens(K):
         sage: L_star = W(M(L.list()).transpose().list())
         sage: pi_cone.contains(L) ==  pi_star.contains(L_star)
         True
+
+    The Lyapunov rank of the positive operator cone is the product of
+    the Lyapunov ranks of the associated cones if they're all proper::
+
+        sage: K1 = random_cone(max_ambient_dim=4,
+        ....:                  strictly_convex=True,
+        ....:                  solid=True)
+        sage: K2 = random_cone(max_ambient_dim=4,
+        ....:                  strictly_convex=True,
+        ....:                  solid=True)
+        sage: pi_K1_K2 = positive_operator_gens(K1,K2)
+        sage: L = ToricLattice(K1.lattice_dim() * K2.lattice_dim())
+        sage: pi_cone = Cone([ g.list() for g in pi_K1_K2 ],
+        ....:                lattice=L,
+        ....:                check=False)
+        sage: beta1 = K1.lyapunov_rank()
+        sage: beta2 = K2.lyapunov_rank()
+        sage: pi_cone.lyapunov_rank() == beta1*beta2
+        True
+
     """
+    if K2 is None:
+        K2 = K1
+
     # Matrices are not vectors in Sage, so we have to convert them
     # to vectors explicitly before we can find a basis. We need these
     # two values to construct the appropriate "long vector" space.
-    F = K.lattice().base_field()
-    n = K.lattice_dim()
+    F = K1.lattice().base_field()
+    n = K1.lattice_dim()
+    m = K2.lattice_dim()
 
-    tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
+    tensor_products = [ s.tensor_product(x) for x in K1 for s in K2.dual() ]
 
     # Convert those tensor products to long vectors.
-    W = VectorSpace(F, n**2)
+    W = VectorSpace(F, n*m)
     vectors = [ W(tp.list()) for tp in tensor_products ]
 
     check = True
-    if K.is_proper():
+    if K1.is_proper() and K2.is_proper():
         # All of the generators involved are extreme vectors and
         # therefore minimal [Tam]_. If this cone is neither solid nor
         # strictly convex, then the tensor product of ``s`` and ``x``
@@ -426,7 +467,7 @@ def positive_operator_gens(K):
     pi_cone = pi_dual.dual()
 
     # And finally convert its rays back to matrix representations.
-    M = MatrixSpace(F, n)
+    M = MatrixSpace(F, m, n)
     return [ M(v.list()) for v in pi_cone ]