Add the new mjo.cone.faces module with the face_generated_by() function.

diff --git a/mjo/cone/all.py b/mjo/cone/all.py
index c8957a0ee58f704d71ffb6cc4283893c78073d82..00ceb9ead486d54b561ded31582afb065810d012 100644 (file)
--- a/mjo/cone/all.py
+++ b/mjo/cone/all.py
@@ -11,6 +11,7 @@ addsitedir(abspath('../../'))
 from mjo.cone.cone import *
 from mjo.cone.completely_positive import *
 from mjo.cone.doubly_nonnegative import *
+from mjo.cone.faces import *
 from mjo.cone.permutation_invariant import *
 from mjo.cone.rearrangement import *
 from mjo.cone.symmetric_psd import *

diff --git a/mjo/cone/faces.py b/mjo/cone/faces.py
new file mode 100644 (file)
index 0000000..baa0708
--- /dev/null
+++ b/mjo/cone/faces.py
@@ -0,0 +1,91 @@
+from sage.all import *
+
+def face_generated_by(K,S):
+    r"""
+    Return the intersection of all faces of ``K`` that contain ``S``.
+
+    REFERENCES:
+
+    .. [Tam] Bit-Shun Tam. On the duality operator of a convex cone. Linear
+       Algebra and its Applications, 64:33-56, 1985, doi:10.1016/0024-3795(85)
+       90265-4.
+
+    SETUP::
+
+        sage: from mjo.cone.faces import face_generated_by
+
+    EXAMPLES:
+
+    The face generated by a standard basis vector in the nonnegative
+    orthant should be the cone consisting of all nonnegative multiples
+    of that standard basis vector::
+
+        sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
+        sage: F = Cone([(1,0,0)])
+        sage: face_generated_by(K,F) == F
+        True
+
+    If we take a nontrivial combination of standard basis vectors in the
+    nonnegative orthant, then the face that the combination generates
+    should be the cone generated by those two basis vectors::
+
+        sage: e1 = vector(QQ, [1,0,0])
+        sage: e2 = vector(QQ, [0,1,0])
+        sage: e3 = vector(QQ, [0,0,1])
+        sage: K = Cone([e1,e2,e2])
+        sage: F = [e1/2 + e2/2]
+        sage: face_generated_by(K,F).is_equivalent(Cone([e1,e2]))
+        True
+
+     An error is raised if ``S`` is not contained in ``K``::
+
+        sage: K = Cone([(1,0),(0,1)])
+        sage: S = [(1,1), (-1,3)]
+        sage: face_generated_by(K,S)
+        Traceback (most recent call last):
+        ...
+        ValueError: S is not a subset of the cone
+
+    TESTS:
+
+    The face generated by a set should be a face::
+
+        sage: set_random_seed()
+        sage: K = random_cone(max_ambient_dim=8, max_rays=10)
+        sage: S = [K.random_element() for idx in range(0,5)]
+        sage: F = face_generated_by(K, S)
+        sage: F.is_face_of(K)
+        True
+
+    The face generated a set should always contain that set::
+
+        sage: set_random_seed()
+        sage: K = random_cone(max_ambient_dim=8, max_rays=10)
+        sage: S = [K.random_element() for idx in range(0,5)]
+        sage: F = face_generated_by(K, S)
+        sage: all([F.contains(x) for x in S])
+        True
+
+    The generators of a proper cone are all extreme vectors of the cone,
+    and therefore generate their own faces::
+
+        sage: set_random_seed()
+        sage: K = random_cone(max_ambient_dim=8,
+        ....:                 max_rays=10,
+        ....:                 strictly_convex=True,
+        ....:                 solid=True)
+        sage: all([face_generated_by(K, [r]) == Cone([r]) for r in K])
+        True
+
+    """
+    face_lattice = K.face_lattice()
+
+    # A list comprehension doesn't work and don't ask me why.
+    candidates = [F for F in face_lattice if all([F.contains(x) for x in S])]
+
+    # K itself is a face of K, so unless we were given a set S that
+    # isn't a subset of K, the candidates list will be nonempty.
+    if len(candidates) == 0:
+        raise ValueError('S is not a subset of the cone')
+    else:
+        return face_lattice.sorted(candidates)[0]