From bc4804ef6cbf4b5f12fb654193e133f719dcc623 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Sun, 2 Nov 2014 22:03:52 -0500
Subject: [PATCH] Add the E() basis matrix function.

---
 mjo/cone/doubly_nonnegative.py | 61 ++++++++++++++++++++++++++++++++++
 1 file changed, 61 insertions(+)

diff --git a/mjo/cone/doubly_nonnegative.py b/mjo/cone/doubly_nonnegative.py
index 89bacb9..03d23b4 100644
--- a/mjo/cone/doubly_nonnegative.py
+++ b/mjo/cone/doubly_nonnegative.py
@@ -122,6 +122,67 @@ def has_admissible_extreme_rank(A):
         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
-- 
2.44.2