mjo/cone/cone.py

   1 from sage.all import *
   2
   3 def is_positive_on(L,K):
   4     r"""
   5     Determine whether or not ``L`` is positive on ``K``.
   6
   7     We say that ``L`` is positive on ``K`` if `L\left\lparen x
   8     \right\rparen` belongs to ``K`` for all `x` in ``K``. This
   9     property need only be checked for generators of ``K``.
  10
  11     INPUT:
  12
  13     - ``L`` -- A linear transformation or matrix.
  14
  15     - ``K`` -- A polyhedral closed convex cone.
  16
  17     OUTPUT:
  18
  19     ``True`` if it can be proven that ``L`` is positive on ``K``,
  20     and ``False`` otherwise.
  21
  22     .. WARNING::
  23
  24         If this function returns ``True``, then ``L`` is positive
  25         on ``K``. However, if ``False`` is returned, that could mean one
  26         of two things. The first is that ``L`` is definitely not
  27         positive on ``K``. The second is more of an "I don't know"
  28         answer, returned (for example) if we cannot prove that an inner
  29         product is nonnegative.
  30
  31     EXAMPLES:
  32
  33     Positive operators on the nonnegative orthant are nonnegative
  34     matrices::
  35
  36         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
  37         sage: L = random_matrix(QQ,3).apply_map(abs)
  38         sage: is_positive_on(L,K)
  39         True
  40
  41     TESTS:
  42
  43     The identity is always positive in a nontrivial space::
  44
  45         sage: set_random_seed()
  46         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
  47         sage: L = identity_matrix(K.lattice_dim())
  48         sage: is_positive_on(L,K)
  49         True
  50
  51     As is the "zero" transformation::
  52
  53         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
  54         sage: R = K.lattice().vector_space().base_ring()
  55         sage: L = zero_matrix(R, K.lattice_dim())
  56         sage: is_positive_on(L,K)
  57         True
  58
  59     Everything in ``K.positive_operators_gens()`` should be
  60     positive on ``K``::
  61
  62         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
  63         sage: all([ is_positive_on(L,K)
  64         ....:       for L in K.positive_operators_gens() ])
  65         True
  66         sage: all([ is_positive_on(L.change_ring(SR),K)
  67         ....:       for L in K.positive_operators_gens() ])
  68         True
  69
  70     """
  71     if L.base_ring().is_exact():
  72         # This could potentially be extended to other types of ``K``...
  73         return all([ L*x in K for x in K ])
  74     elif L.base_ring() is SR:
  75         # Fall back to inequality-checking when the entries of ``L``
  76         # might be symbolic.
  77         return all([ s*(L*x) >= 0 for x in K for s in K.dual() ])
  78     else:
  79         # The only inexact ring that we're willing to work with is SR,
  80         # since it can still be exact when working with symbolic
  81         # constants like pi and e.
  82         raise ValueError('base ring of operator L is neither SR nor exact')
  83
  84
  85 def is_cross_positive_on(L,K):
  86     r"""
  87     Determine whether or not ``L`` is cross-positive on ``K``.
  88
  89     We say that ``L`` is cross-positive on ``K`` if `\left\langle
  90     L\left\lparenx\right\rparen,s\right\rangle \ge 0` for all pairs
  91     `\left\langle x,s \right\rangle` in the complementarity set of
  92     ``K``. This property need only be checked for generators of
  93     ``K`` and its dual.
  94
  95     INPUT:
  96
  97     - ``L`` -- A linear transformation or matrix.
  98
  99     - ``K`` -- A polyhedral closed convex cone.
 100
 101     OUTPUT:
 102
 103     ``True`` if it can be proven that ``L`` is cross-positive on ``K``,
 104     and ``False`` otherwise.
 105
 106     .. WARNING::
 107
 108         If this function returns ``True``, then ``L`` is cross-positive
 109         on ``K``. However, if ``False`` is returned, that could mean one
 110         of two things. The first is that ``L`` is definitely not
 111         cross-positive on ``K``. The second is more of an "I don't know"
 112         answer, returned (for example) if we cannot prove that an inner
 113         product is nonnegative.
 114
 115     EXAMPLES:
 116
 117     The identity is always cross-positive in a nontrivial space::
 118
 119         sage: set_random_seed()
 120         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
 121         sage: L = identity_matrix(K.lattice_dim())
 122         sage: is_cross_positive_on(L,K)
 123         True
 124
 125     As is the "zero" transformation::
 126
 127         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
 128         sage: R = K.lattice().vector_space().base_ring()
 129         sage: L = zero_matrix(R, K.lattice_dim())
 130         sage: is_cross_positive_on(L,K)
 131         True
 132
 133     TESTS:
 134
 135     Everything in ``K.cross_positive_operators_gens()`` should be
 136     cross-positive on ``K``::
 137
 138         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
 139         sage: all([ is_cross_positive_on(L,K)
 140         ....:       for L in K.cross_positive_operators_gens() ])
 141         True
 142         sage: all([ is_cross_positive_on(L.change_ring(SR),K)
 143         ....:       for L in K.cross_positive_operators_gens() ])
 144         True
 145
 146     """
 147     if L.base_ring().is_exact() or L.base_ring() is SR:
 148         return all([ s*(L*x) >= 0
 149                      for (x,s) in K.discrete_complementarity_set() ])
 150     else:
 151         # The only inexact ring that we're willing to work with is SR,
 152         # since it can still be exact when working with symbolic
 153         # constants like pi and e.
 154         raise ValueError('base ring of operator L is neither SR nor exact')
 155
 156
 157 def is_Z_on(L,K):
 158     r"""
 159     Determine whether or not ``L`` is a Z-operator on ``K``.
 160
 161     We say that ``L`` is a Z-operator on ``K`` if `\left\langle
 162     L\left\lparenx\right\rparen,s\right\rangle \le 0` for all pairs
 163     `\left\langle x,s \right\rangle` in the complementarity set of
 164     ``K``. It is known that this property need only be
 165     checked for generators of ``K`` and its dual.
 166
 167     A matrix is a Z-operator on ``K`` if and only if its negation is a
 168     cross-positive operator on ``K``.
 169
 170     INPUT:
 171
 172     - ``L`` -- A linear transformation or matrix.
 173
 174     - ``K`` -- A polyhedral closed convex cone.
 175
 176     OUTPUT:
 177
 178     ``True`` if it can be proven that ``L`` is a Z-operator on ``K``,
 179     and ``False`` otherwise.
 180
 181     .. WARNING::
 182
 183         If this function returns ``True``, then ``L`` is a Z-operator
 184         on ``K``. However, if ``False`` is returned, that could mean one
 185         of two things. The first is that ``L`` is definitely not
 186         a Z-operator on ``K``. The second is more of an "I don't know"
 187         answer, returned (for example) if we cannot prove that an inner
 188         product is nonnegative.
 189
 190     EXAMPLES:
 191
 192     The identity is always a Z-operator in a nontrivial space::
 193
 194         sage: set_random_seed()
 195         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
 196         sage: L = identity_matrix(K.lattice_dim())
 197         sage: is_Z_on(L,K)
 198         True
 199
 200     As is the "zero" transformation::
 201
 202         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
 203         sage: R = K.lattice().vector_space().base_ring()
 204         sage: L = zero_matrix(R, K.lattice_dim())
 205         sage: is_Z_on(L,K)
 206         True
 207
 208     TESTS:
 209
 210     Everything in ``K.Z_operators_gens()`` should be a Z-operator
 211     on ``K``::
 212
 213         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
 214         sage: all([ is_Z_on(L,K)
 215         ....:       for L in K.Z_operators_gens() ])
 216         True
 217         sage: all([ is_Z_on(L.change_ring(SR),K)
 218         ....:       for L in K.Z_operators_gens() ])
 219         True
 220
 221     """
 222     return is_cross_positive_on(-L,K)
 223
 224
 225 def is_lyapunov_like_on(L,K):
 226     r"""
 227     Determine whether or not ``L`` is Lyapunov-like on ``K``.
 228
 229     We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
 230     L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
 231     `\left\langle x,s \right\rangle` in the complementarity set of
 232     ``K``. This property need only be checked for generators of
 233     ``K`` and its dual.
 234
 235     INPUT:
 236
 237     - ``L`` -- A linear transformation or matrix.
 238
 239     - ``K`` -- A polyhedral closed convex cone.
 240
 241     OUTPUT:
 242
 243     ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
 244     and ``False`` otherwise.
 245
 246     .. WARNING::
 247
 248         If this function returns ``True``, then ``L`` is Lyapunov-like
 249         on ``K``. However, if ``False`` is returned, that could mean one
 250         of two things. The first is that ``L`` is definitely not
 251         Lyapunov-like on ``K``. The second is more of an "I don't know"
 252         answer, returned (for example) if we cannot prove that an inner
 253         product is zero.
 254
 255     EXAMPLES:
 256
 257     Lyapunov-like operators on the nonnegative orthant are diagonal
 258     matrices::
 259
 260         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 261         sage: L = diagonal_matrix(random_vector(QQ,3))
 262         sage: is_lyapunov_like_on(L,K)
 263         True
 264
 265     TESTS:
 266
 267     The identity is always Lyapunov-like in a nontrivial space::
 268
 269         sage: set_random_seed()
 270         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
 271         sage: L = identity_matrix(K.lattice_dim())
 272         sage: is_lyapunov_like_on(L,K)
 273         True
 274
 275     As is the "zero" transformation::
 276
 277         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
 278         sage: R = K.lattice().vector_space().base_ring()
 279         sage: L = zero_matrix(R, K.lattice_dim())
 280         sage: is_lyapunov_like_on(L,K)
 281         True
 282
 283     Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
 284     on ``K``::
 285
 286         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
 287         sage: all([ is_lyapunov_like_on(L,K)
 288         ....:       for L in K.lyapunov_like_basis() ])
 289         True
 290         sage: all([ is_lyapunov_like_on(L.change_ring(SR),K)
 291         ....:       for L in K.lyapunov_like_basis() ])
 292         True
 293
 294     """
 295     if L.base_ring().is_exact() or L.base_ring() is SR:
 296         # The "fast method" of creating a vector space based on a
 297         # ``lyapunov_like_basis`` is actually slower than this.
 298         return all([ s*(L*x) == 0
 299                      for (x,s) in K.discrete_complementarity_set() ])
 300     else:
 301         # The only inexact ring that we're willing to work with is SR,
 302         # since it can still be exact when working with symbolic
 303         # constants like pi and e.
 304         raise ValueError('base ring of operator L is neither SR nor exact')
 305
 306 def LL_cone(K):
 307     gens = K.lyapunov_like_basis()
 308     L = ToricLattice(K.lattice_dim()**2)
 309     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 310
 311 def Sigma_cone(K):
 312     gens = K.cross_positive_operators_gens()
 313     L = ToricLattice(K.lattice_dim()**2)
 314     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 315
 316 def Z_cone(K):
 317     gens = K.Z_operators_gens()
 318     L = ToricLattice(K.lattice_dim()**2)
 319     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 320
 321 def pi_cone(K1, K2=None):
 322     if K2 is None:
 323         K2 = K1
 324     gens = K1.positive_operators_gens(K2)
 325     L = ToricLattice(K1.lattice_dim()*K2.lattice_dim())
 326     return Cone([ g.list() for g in gens ], lattice=L, check=False)