mjo/cone/schur.py

   1 """
   2 The Schur cone, as described in the "Critical angles..." papers by
   3 Iusem, Seeger, and Sossa. It defines the Schur ordering on `R^{n}`.
   4 """
   5
   6 from sage.all import *
   7
   8 def schur_cone(n, lattice=None):
   9     r"""
  10     Return the Schur cone in ``n`` dimensions that induces the
  11     majorization ordering.
  12
  13     INPUT:
  14
  15       - ``n`` -- the dimension the ambient space.
  16
  17       - ``lattice`` -- (default: ``None``) an ambient lattice of rank ``n``
  18                        to be passed to the :func:`Cone` constructor.
  19
  20     OUTPUT:
  21
  22     A rational closed convex Schur cone of dimension ``n``. Each
  23     generating ray will have the integer ring as its base ring.
  24
  25     If a ``lattice`` was specified, then the resulting cone will live in
  26     that lattice unless its rank is incompatible with the dimension
  27     ``n`` (in which case a ``ValueError`` is raised).
  28
  29     REFERENCES:
  30
  31     .. [GourionSeeger] Daniel Gourion and Alberto Seeger.
  32        Critical angles in polyhedral convex cones: numerical and
  33        statistical considerations. Mathematical Programming, 123:173-198,
  34        2010, doi:10.1007/s10107-009-0317-2.
  35
  36     .. [IusemSeegerOnPairs] Alfredo Iusem and Alberto Seeger.
  37        On pairs of vectors achieving the maximal angle of a convex cone.
  38        Mathematical Programming, 104(2-3):501-523, 2005,
  39        doi:10.1007/s10107-005-0626-z.
  40
  41     .. [SeegerSossaI] Alberto Seeger and David Sossa.
  42        Critical angles between two convex cones I. General theory.
  43        TOP, 24(1):44-65, 2016, doi:10.1007/s11750-015-0375-y.
  44
  45     SETUP::
  46
  47         sage: from mjo.cone.nonnegative_orthant import nonnegative_orthant
  48         sage: from mjo.cone.schur import schur_cone
  49
  50     EXAMPLES:
  51
  52     Verify the claim that the maximal angle between any two generators
  53     of the Schur cone and the nonnegative quintant is ``3*pi/4``::
  54
  55         sage: P = schur_cone(5)
  56         sage: Q = nonnegative_orthant(5)
  57         sage: G = ( g.change_ring(QQbar).normalized() for g in P )
  58         sage: H = ( h.change_ring(QQbar).normalized() for h in Q )
  59         sage: actual = max(arccos(u.inner_product(v)) for u in G for v in H)
  60         sage: expected = 3*pi/4
  61         sage: abs(actual - expected).n() < 1e-12
  62         True
  63
  64     The dual of the Schur cone is the "downward monotonic cone"
  65     [GourionSeeger]_, whose elements' entries are in non-increasing
  66     order::
  67
  68         sage: set_random_seed()
  69         sage: n = ZZ.random_element(10)
  70         sage: K = schur_cone(n).dual()
  71         sage: x = K.random_element()
  72         sage: all( x[i] >= x[i+1] for i in range(n-1) )
  73         True
  74
  75     TESTS:
  76
  77     We get the trivial cone when ``n`` is zero::
  78
  79         sage: schur_cone(0).is_trivial()
  80         True
  81
  82     The Schur cone induces the majorization ordering::
  83
  84         sage: set_random_seed()
  85         sage: def majorized_by(x,y):
  86         ....:     return (all(sum(x[0:i]) <= sum(y[0:i])
  87         ....:                 for i in range(x.degree()-1))
  88         ....:             and sum(x) == sum(y))
  89         sage: n = ZZ.random_element(10)
  90         sage: V = VectorSpace(QQ, n)
  91         sage: S = schur_cone(n)
  92         sage: majorized_by(V.zero(), S.random_element())
  93         True
  94         sage: x = V.random_element()
  95         sage: y = V.random_element()
  96         sage: majorized_by(x,y) == ( (y-x) in S )
  97         True
  98
  99     If a ``lattice`` was given, it is actually used::
 100
 101         sage: L = ToricLattice(3, 'M')
 102         sage: schur_cone(3, lattice=L)
 103         2-d cone in 3-d lattice M
 104
 105     Unless the rank of the lattice disagrees with ``n``::
 106
 107         sage: L = ToricLattice(1, 'M')
 108         sage: schur_cone(3, lattice=L)
 109         Traceback (most recent call last):
 110         ...
 111         ValueError: lattice rank=1 and dimension n=3 are incompatible
 112
 113     """
 114     if lattice is None:
 115         lattice = ToricLattice(n)
 116
 117     if lattice.rank() != n:
 118         raise ValueError('lattice rank=%d and dimension n=%d are incompatible'
 119                          %
 120                          (lattice.rank(), n))
 121
 122     def _f(i,j):
 123         if i == j:
 124             return 1
 125         elif j - i == 1:
 126             return -1
 127         else:
 128             return 0
 129
 130     # The "max" below catches the trivial case where n == 0.
 131     S = matrix(ZZ, max(0,n-1), n, _f)
 132
 133     return Cone(S.rows(), lattice)