From a6b746e90c7710c76de9a559258f4769f0111279 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Fri, 2 Nov 2018 10:45:37 -0400
Subject: [PATCH] cone/schur.py: clean up and add some tests.

---
 mjo/cone/all.py   |  1 +
 mjo/cone/schur.py | 75 ++++++++++++++++++++++++++++++-----------------
 2 files changed, 49 insertions(+), 27 deletions(-)

diff --git a/mjo/cone/all.py b/mjo/cone/all.py
index e7d3bff..80d5b0b 100644
--- a/mjo/cone/all.py
+++ b/mjo/cone/all.py
@@ -9,4 +9,5 @@ from mjo.cone.faces import *
 from mjo.cone.nonnegative_orthant import *
 from mjo.cone.permutation_invariant import *
 from mjo.cone.rearrangement import *
+from mjo.cone.schur import *
 from mjo.cone.symmetric_psd import *
diff --git a/mjo/cone/schur.py b/mjo/cone/schur.py
index bea76b6..fd806d9 100644
--- a/mjo/cone/schur.py
+++ b/mjo/cone/schur.py
@@ -1,45 +1,66 @@
 """
 The Schur cone, as described in the "Critical angles..." papers by
-Seeger and Sossa.
+Iusem, Seeger, and Sossa. It defines the Schur ordering on `R^{n}`.
 """
 
 from sage.all import *
 
-def is_pointed(P,Q):
-    newP = Cone([ list(p) for p in P ])
-    newQ = Cone([ list(q) for q in Q ])
-    return newP.intersection(-newQ).is_trivial()
-
-def ext_angles(P,Q):
-    angles = []
-    for p in P:
-        for q in Q:
-            p = vector(SR, p)
-            q = vector(SR, q)
-            p = p/p.norm()
-            q = q/q.norm()
-            angles.append(arccos(p.inner_product(q)))
-
-    return angles
-
-def schur(n):
+def schur_cone(n):
     r"""
     Return the Schur cone in ``n`` dimensions.
 
     INPUT:
 
-    - ``n`` -- the ambient dimension of the Schur cone and its ambient space.
+    - ``n`` -- the dimension the ambient space.
 
     OUTPUT:
 
     A rational closed convex Schur cone of dimension ``n``.
+
+    REFERENCES:
+
+    .. [SeegerSossaI] Alberto Seeger and David Sossa.
+       Critical angles between two convex cones I. General theory.
+       TOP, 24(1):44-65, 2016, doi:10.1007/s11750-015-0375-y.
+
+    SETUP::
+
+        sage: from mjo.cone.nonnegative_orthant import nonnegative_orthant
+        sage: from mjo.cone.schur import schur_cone
+
+    EXAMPLES:
+
+    Verify the claim that the maximal angle between any two generators
+    of the Schur cone and the nonnegative quintant is ``3*pi/4``::
+
+        sage: P = schur_cone(5)
+        sage: Q = nonnegative_orthant(5)
+        sage: G = [ g.change_ring(QQbar).normalized() for g in P ]
+        sage: H = [ h.change_ring(QQbar).normalized() for h in Q ]
+        sage: actual = max([arccos(u.inner_product(v)) for u in G for v in H])
+        sage: expected = 3*pi/4
+        sage: abs(actual - expected).n() < 1e-12
+        True
+
+    TESTS:
+
+    We get the trivial cone when ``n`` is zero::
+
+        sage: schur_cone(0).is_trivial()
+        True
+
     """
 
-    hs = []
-    for i in range(1,n):
-        h_i = [0]*n
-        h_i[i-1] = QQ(1)
-        h_i[i] = -QQ(1)
-        hs.append(vector(QQ,n,h_i))
+    def _f(i,j):
+        if i == j:
+            return 1
+        elif j - i == 1:
+            return -1
+        else:
+            return 0
+
+    # The "max" below catches the trivial case where n == 0.
+    S = matrix(ZZ, max(0,n-1), n, _f)
 
-    return Cone(hs)
+    # Likewise, when n == 0, we need to specify the lattice.
+    return Cone(S.rows(), ToricLattice(n))
-- 
2.43.2