"""
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, lattice=None):
r"""
- Return the Schur cone in ``n`` dimensions.
+ Return the Schur cone in ``n`` dimensions that induces the
+ majorization ordering.
INPUT:
- - ``n`` -- the ambient dimension of the Schur cone and its ambient space.
+ - ``n`` -- the dimension the ambient space.
+
+ - ``lattice`` -- (default: ``None``) an ambient lattice of rank ``n``
+ to be passed to the :func:`Cone` constructor.
OUTPUT:
- A rational closed convex Schur cone of dimension ``n``.
+ A rational closed convex Schur cone of dimension ``n``. Each
+ generating ray will have the integer ring as its base ring.
+
+ If a ``lattice`` was specified, then the resulting cone will live in
+ that lattice unless its rank is incompatible with the dimension
+ ``n`` (in which case a ``ValueError`` is raised).
+
+ REFERENCES:
+
+ .. [GourionSeeger] Daniel Gourion and Alberto Seeger.
+ Critical angles in polyhedral convex cones: numerical and
+ statistical considerations. Mathematical Programming, 123:173-198,
+ 2010, doi:10.1007/s10107-009-0317-2.
+
+ .. [IusemSeegerOnPairs] Alfredo Iusem and Alberto Seeger.
+ On pairs of vectors achieving the maximal angle of a convex cone.
+ Mathematical Programming, 104(2-3):501-523, 2005,
+ doi:10.1007/s10107-005-0626-z.
+
+ .. [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
+
+ The dual of the Schur cone is the "downward monotonic cone"
+ [GourionSeeger]_, whose elements' entries are in non-increasing
+ order::
+
+ sage: n = ZZ.random_element(10)
+ sage: K = schur_cone(n).dual()
+ sage: x = K.random_element()
+ sage: all( x[i] >= x[i+1] for i in range(n-1) )
+ True
+
+ TESTS:
+
+ We get the trivial cone when ``n`` is zero::
+
+ sage: schur_cone(0).is_trivial()
+ True
+
+ The Schur cone induces the majorization ordering::
+
+ sage: def majorized_by(x,y):
+ ....: return (all(sum(x[0:i]) <= sum(y[0:i])
+ ....: for i in range(x.degree()-1))
+ ....: and sum(x) == sum(y))
+ sage: n = ZZ.random_element(10)
+ sage: V = VectorSpace(QQ, n)
+ sage: S = schur_cone(n)
+ sage: majorized_by(V.zero(), S.random_element())
+ True
+ sage: x = V.random_element()
+ sage: y = V.random_element()
+ sage: majorized_by(x,y) == ( (y-x) in S )
+ True
+
+ If a ``lattice`` was given, it is actually used::
+
+ sage: L = ToricLattice(3, 'M')
+ sage: schur_cone(3, lattice=L)
+ 2-d cone in 3-d lattice M
+
+ Unless the rank of the lattice disagrees with ``n``::
+
+ sage: L = ToricLattice(1, 'M')
+ sage: schur_cone(3, lattice=L)
+ Traceback (most recent call last):
+ ...
+ ValueError: lattice rank=1 and dimension n=3 are incompatible
+
"""
+ if lattice is None:
+ lattice = ToricLattice(n)
+
+ if lattice.rank() != n:
+ raise ValueError('lattice rank=%d and dimension n=%d are incompatible'
+ %
+ (lattice.rank(), n))
+
+ def _f(i,j):
+ if i == j:
+ return 1
+ elif j - i == 1:
+ return -1
+ else:
+ return 0
- 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))
+ # The "max" below catches the trivial case where n == 0.
+ S = matrix(ZZ, max(0,n-1), n, _f)
- return Cone(hs)
+ return Cone(S.rows(), lattice)
projects / sage.d.git / blobdiff