X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fdoubly_nonnegative.py;h=03d23b4ddcbd785b48fe9dd447876cefb586a067;hb=bc4804ef6cbf4b5f12fb654193e133f719dcc623;hp=f705e5e98b083d6cf85fac39148d171cdf452d99;hpb=2d95c4e34d085c9647c73b73b2957f937cfee26b;p=sage.d.git

diff --git a/mjo/cone/doubly_nonnegative.py b/mjo/cone/doubly_nonnegative.py
index f705e5e..03d23b4 100644
--- a/mjo/cone/doubly_nonnegative.py
+++ b/mjo/cone/doubly_nonnegative.py
@@ -29,7 +29,7 @@ def is_doubly_nonnegative(A):
 
     INPUT:
 
-      - ``A`` - The matrix in question
+    - ``A`` - The matrix in question
 
     OUTPUT:
 
@@ -67,9 +67,144 @@ def is_doubly_nonnegative(A):
 
 
 
+def has_admissible_extreme_rank(A):
+    """
+    The extreme matrices of the doubly-nonnegative cone have some
+    restrictions on their ranks. This function checks to see whether or
+    not ``A`` could be extreme based on its rank.
+
+    INPUT:
+
+    - ``A`` - The matrix in question
+
+    OUTPUT:
+
+    ``False`` if the rank of ``A`` precludes it from being an extreme
+    matrix of the doubly-nonnegative cone, ``True`` otherwise.
+
+    REFERENCE:
+
+    Hamilton-Jester, Crista Lee; Li, Chi-Kwong. Extreme Vectors of
+    Doubly Nonnegative Matrices. Rocky Mountain Journal of Mathematics
+    26 (1996), no. 4, 1371--1383. doi:10.1216/rmjm/1181071993.
+    http://projecteuclid.org/euclid.rmjm/1181071993.
+
+    EXAMPLES:
+
+    The zero matrix has rank zero, which is admissible::
+
+    sage: A = zero_matrix(QQ, 5, 5)
+    sage: has_admissible_extreme_rank(A)
+    True
+
+    """
+    if not A.is_symmetric():
+        # This function is more or less internal, so blow up if passed
+        # something unexpected.
+        raise ValueError('The matrix ``A`` must be symmetric.')
+
+    r = rank(A)
+    n = ZZ(A.nrows()) # Columns would work, too, since ``A`` is symmetric.
+
+    if r == 0:
+        # Zero is in the doubly-nonnegative cone.
+        return True
+
+    # See Theorem 3.1 in the cited reference.
+    if r == 2:
+        return False
+
+    if n.mod(2) == 0:
+        # n is even
+        return r <= max(1, n-3)
+    else:
+        # n is odd
+        return r <= max(1, n-2)
+
+
+def E(matrix_space, i,j):
+    """
+    Return the ``i``,``j``th element of the standard basis in
+    ``matrix_space``.
+
+    INPUT:
+
+    - ``matrix_space`` - The underlying matrix space of whose basis
+                         the returned matrix is an element
+
+    - ``i`` - The row index of the single nonzero entry
+
+    - ``j`` - The column index of the single nonzero entry
+
+    OUTPUT:
+
+    A basis element of ``matrix_space``. It has a single \"1\" in the
+    ``i``,``j`` row,column and zeros elsewhere.
+
+    EXAMPLES::
+
+        sage: M = MatrixSpace(ZZ, 2, 2)
+        sage: E(M,0,0)
+        [1 0]
+        [0 0]
+        sage: E(M,0,1)
+        [0 1]
+        [0 0]
+        sage: E(M,1,0)
+        [0 0]
+        [1 0]
+        sage: E(M,1,1)
+        [0 0]
+        [0 1]
+        sage: E(M,2,1)
+        Traceback (most recent call last):
+        ...
+        IndexError: Index `i` is out of bounds.
+        sage: E(M,1,2)
+        Traceback (most recent call last):
+        ...
+        IndexError: Index `j` is out of bounds.
+
+    """
+    # We need to check these ourselves, see below.
+    if i >= matrix_space.nrows():
+        raise IndexError('Index `i` is out of bounds.')
+    if j >= matrix_space.ncols():
+        raise IndexError('Index `j` is out of bounds.')
+
+    # The basis here is returned as a one-dimensional list, so we need
+    # to compute the offset into it based on ``i`` and ``j``. Since we
+    # compute the index ourselves, we need to do bounds-checking
+    # manually. Otherwise for e.g. a 2x2 matrix space, the index (0,2)
+    # would be computed as offset 3 into a four-element list and we
+    # would succeed incorrectly.
+    idx = matrix_space.ncols()*i + j
+    return matrix_space.basis()[idx]
+
+
+
 def is_extreme_doubly_nonnegative(A):
     """
     Returns ``True`` if the given matrix is an extreme matrix of the
     doubly-nonnegative cone, and ``False`` otherwise.
+
+    EXAMPLES:
+
+    The zero matrix is an extreme matrix::
+
+        sage: A = zero_matrix(QQ, 5, 5)
+        sage: is_extreme_doubly_nonnegative(A)
+        True
+
     """
+
+    r = A.rank()
+
+    if r == 0:
+        # Short circuit, we know the zero matrix is extreme.
+        return True
+
+    if not is_admissible_extreme_rank(r):
+        return False
+
     raise NotImplementedError()