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