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