mjo/cone/cone.py

   1 from sage.all import *
   2 from sage.geometry.cone import is_Cone
   3
   4 def is_positive_on(L,K):
   5     r"""
   6     Determine whether or not ``L`` is positive on ``K``.
   7
   8     We say that ``L`` is positive on a closed convex cone ``K`` if
   9     `L\left\lparen x \right\rparen` belongs to ``K`` for all `x` in
  10     ``K``. This property need only be checked for generators of ``K``.
  11
  12     To reliably check whether or not ``L`` is positive, its base ring
  13     must be either exact (for example, the rationals) or ``SR``. An
  14     exact ring is more reliable, but in some cases a matrix whose
  15     entries contain symbolic constants like ``e`` and ``pi`` will work.
  16
  17     INPUT:
  18
  19     - ``L`` -- A matrix over either an exact ring or ``SR``.
  20
  21     - ``K`` -- A polyhedral closed convex cone.
  22
  23     OUTPUT:
  24
  25     If the base ring of ``L`` is exact, then ``True`` will be returned if
  26     and only if ``L`` is positive on ``K``.
  27
  28     If the base ring of ``L`` is ``SR``, then the situation is more
  29     complicated:
  30
  31     - ``True`` will be returned if it can be proven that ``L``
  32       is positive on ``K``.
  33     - ``False`` will be returned if it can be proven that ``L``
  34       is not positive on ``K``.
  35     - ``False`` will also be returned if we can't decide; specifically
  36       if we arrive at a symbolic inequality that cannot be resolved.
  37
  38     EXAMPLES:
  39
  40     Nonnegative matrices are positive operators on the nonnegative
  41     orthant::
  42
  43         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
  44         sage: L = random_matrix(QQ,3).apply_map(abs)
  45         sage: is_positive_on(L,K)
  46         True
  47
  48     TESTS:
  49
  50     The identity operator is always positive::
  51
  52         sage: set_random_seed()
  53         sage: K = random_cone(max_ambient_dim=8)
  54         sage: L = identity_matrix(K.lattice_dim())
  55         sage: is_positive_on(L,K)
  56         True
  57
  58     The "zero" operator is always positive::
  59
  60         sage: K = random_cone(max_ambient_dim=8)
  61         sage: R = K.lattice().vector_space().base_ring()
  62         sage: L = zero_matrix(R, K.lattice_dim())
  63         sage: is_positive_on(L,K)
  64         True
  65
  66     Everything in ``K.positive_operators_gens()`` should be
  67     positive on ``K``::
  68
  69         sage: K = random_cone(max_ambient_dim=5)
  70         sage: all([ is_positive_on(L,K)                     # long time
  71         ....:       for L in K.positive_operators_gens() ]) # long time
  72         True
  73         sage: all([ is_positive_on(L.change_ring(SR),K)     # long time
  74         ....:       for L in K.positive_operators_gens() ]) # long time
  75         True
  76
  77     Technically we could test this, but for now only closed convex cones
  78     are supported as our ``K`` argument::
  79
  80         sage: L = identity_matrix(3)
  81         sage: K = [ vector([1,2,3]), vector([5,-1,7]) ]
  82         sage: is_positive_on(L,K)
  83         Traceback (most recent call last):
  84         ...
  85         TypeError: K must be a Cone.
  86
  87     """
  88     if not is_Cone(K):
  89         raise TypeError('K must be a Cone.')
  90     if not L.base_ring().is_exact() and not L.base_ring() is SR:
  91         raise ValueError('The base ring of L is neither SR nor exact.')
  92
  93     if L.base_ring().is_exact():
  94         # This should be way faster than computing the dual and
  95         # checking a bunch of inequalities, but it doesn't work if
  96         # ``L*x`` is symbolic. For example, ``e in Cone([(1,)])``
  97         # is true, but returns ``False``.
  98         return all([ L*x in K for x in K ])
  99     else:
 100         # Fall back to inequality-checking when the entries of ``L``
 101         # might be symbolic.
 102         return all([ s*(L*x) >= 0 for x in K for s in K.dual() ])
 103
 104
 105 def is_cross_positive_on(L,K):
 106     r"""
 107     Determine whether or not ``L`` is cross-positive on ``K``.
 108
 109     We say that ``L`` is cross-positive on a closed convex cone``K`` if
 110     `\left\langle L\left\lparenx\right\rparen,s\right\rangle \ge 0` for
 111     all pairs `\left\langle x,s \right\rangle` in the complementarity
 112     set of ``K``. This property need only be checked for generators of
 113     ``K`` and its dual.
 114
 115     To reliably check whether or not ``L`` is cross-positive, its base
 116     ring must be either exact (for example, the rationals) or ``SR``. An
 117     exact ring is more reliable, but in some cases a matrix whose
 118     entries contain symbolic constants like ``e`` and ``pi`` will work.
 119
 120     INPUT:
 121
 122     - ``L`` -- A matrix over either an exact ring or ``SR``.
 123
 124     - ``K`` -- A polyhedral closed convex cone.
 125
 126     OUTPUT:
 127
 128     If the base ring of ``L`` is exact, then ``True`` will be returned if
 129     and only if ``L`` is cross-positive on ``K``.
 130
 131     If the base ring of ``L`` is ``SR``, then the situation is more
 132     complicated:
 133
 134     - ``True`` will be returned if it can be proven that ``L``
 135       is cross-positive on ``K``.
 136     - ``False`` will be returned if it can be proven that ``L``
 137       is not cross-positive on ``K``.
 138     - ``False`` will also be returned if we can't decide; specifically
 139       if we arrive at a symbolic inequality that cannot be resolved.
 140
 141     EXAMPLES:
 142
 143     The identity operator is always cross-positive::
 144
 145         sage: set_random_seed()
 146         sage: K = random_cone(max_ambient_dim=8)
 147         sage: L = identity_matrix(K.lattice_dim())
 148         sage: is_cross_positive_on(L,K)
 149         True
 150
 151     The "zero" operator is always cross-positive::
 152
 153         sage: K = random_cone(max_ambient_dim=8)
 154         sage: R = K.lattice().vector_space().base_ring()
 155         sage: L = zero_matrix(R, K.lattice_dim())
 156         sage: is_cross_positive_on(L,K)
 157         True
 158
 159     TESTS:
 160
 161     Everything in ``K.cross_positive_operators_gens()`` should be
 162     cross-positive on ``K``::
 163
 164         sage: K = random_cone(max_ambient_dim=5)
 165         sage: all([ is_cross_positive_on(L,K)                     # long time
 166         ....:       for L in K.cross_positive_operators_gens() ]) # long time
 167         True
 168         sage: all([ is_cross_positive_on(L.change_ring(SR),K)     # long time
 169         ....:       for L in K.cross_positive_operators_gens() ]) # long time
 170         True
 171
 172     Technically we could test this, but for now only closed convex cones
 173     are supported as our ``K`` argument::
 174
 175         sage: L = identity_matrix(3)
 176         sage: K = [ vector([8,2,-8]), vector([5,-5,7]) ]
 177         sage: is_cross_positive_on(L,K)
 178         Traceback (most recent call last):
 179         ...
 180         TypeError: K must be a Cone.
 181
 182     """
 183     if not is_Cone(K):
 184         raise TypeError('K must be a Cone.')
 185     if not L.base_ring().is_exact() and not L.base_ring() is SR:
 186         raise ValueError('The base ring of L is neither SR nor exact.')
 187
 188     return all([ s*(L*x) >= 0
 189                  for (x,s) in K.discrete_complementarity_set() ])
 190
 191 def is_Z_on(L,K):
 192     r"""
 193     Determine whether or not ``L`` is a Z-operator on ``K``.
 194
 195     We say that ``L`` is a Z-operator on a closed convex cone``K`` if
 196     `\left\langle L\left\lparenx\right\rparen,s\right\rangle \le 0` for
 197     all pairs `\left\langle x,s \right\rangle` in the complementarity
 198     set of ``K``. It is known that this property need only be checked
 199     for generators of ``K`` and its dual.
 200
 201     A matrix is a Z-operator on ``K`` if and only if its negation is a
 202     cross-positive operator on ``K``.
 203
 204     To reliably check whether or not ``L`` is a Z operator, its base
 205     ring must be either exact (for example, the rationals) or ``SR``. An
 206     exact ring is more reliable, but in some cases a matrix whose
 207     entries contain symbolic constants like ``e`` and ``pi`` will work.
 208
 209     INPUT:
 210
 211     - ``L`` -- A matrix over either an exact ring or ``SR``.
 212
 213     - ``K`` -- A polyhedral closed convex cone.
 214
 215     OUTPUT:
 216
 217     If the base ring of ``L`` is exact, then ``True`` will be returned if
 218     and only if ``L`` is a Z-operator on ``K``.
 219
 220     If the base ring of ``L`` is ``SR``, then the situation is more
 221     complicated:
 222
 223     - ``True`` will be returned if it can be proven that ``L``
 224       is a Z-operator on ``K``.
 225     - ``False`` will be returned if it can be proven that ``L``
 226       is not a Z-operator on ``K``.
 227     - ``False`` will also be returned if we can't decide; specifically
 228       if we arrive at a symbolic inequality that cannot be resolved.
 229
 230     EXAMPLES:
 231
 232     The identity operator is always a Z-operator::
 233
 234         sage: set_random_seed()
 235         sage: K = random_cone(max_ambient_dim=8)
 236         sage: L = identity_matrix(K.lattice_dim())
 237         sage: is_Z_on(L,K)
 238         True
 239
 240     The "zero" operator is always a Z-operator::
 241
 242         sage: K = random_cone(max_ambient_dim=8)
 243         sage: R = K.lattice().vector_space().base_ring()
 244         sage: L = zero_matrix(R, K.lattice_dim())
 245         sage: is_Z_on(L,K)
 246         True
 247
 248     TESTS:
 249
 250     Everything in ``K.Z_operators_gens()`` should be a Z-operator
 251     on ``K``::
 252
 253         sage: K = random_cone(max_ambient_dim=5)
 254         sage: all([ is_Z_on(L,K)                     # long time
 255         ....:       for L in K.Z_operators_gens() ]) # long time
 256         True
 257         sage: all([ is_Z_on(L.change_ring(SR),K)     # long time
 258         ....:       for L in K.Z_operators_gens() ]) # long time
 259         True
 260
 261     Technically we could test this, but for now only closed convex cones
 262     are supported as our ``K`` argument::
 263
 264         sage: L = identity_matrix(3)
 265         sage: K = [ vector([-4,20,3]), vector([1,-5,2]) ]
 266         sage: is_Z_on(L,K)
 267         Traceback (most recent call last):
 268         ...
 269         TypeError: K must be a Cone.
 270
 271     """
 272     return is_cross_positive_on(-L,K)
 273
 274
 275 def is_lyapunov_like_on(L,K):
 276     r"""
 277     Determine whether or not ``L`` is Lyapunov-like on ``K``.
 278
 279     We say that ``L`` is Lyapunov-like on a closed convex cone ``K`` if
 280     `\left\langle L\left\lparenx\right\rparen,s\right\rangle = 0` for
 281     all pairs `\left\langle x,s \right\rangle` in the complementarity
 282     set of ``K``. This property need only be checked for generators of
 283     ``K`` and its dual.
 284
 285     To reliably check whether or not ``L`` is Lyapunov-like, its base
 286     ring must be either exact (for example, the rationals) or ``SR``. An
 287     exact ring is more reliable, but in some cases a matrix whose
 288     entries contain symbolic constants like ``e`` and ``pi`` will work.
 289
 290     INPUT:
 291
 292     - ``L`` -- A matrix over either an exact ring or ``SR``.
 293
 294     - ``K`` -- A polyhedral closed convex cone.
 295
 296     OUTPUT:
 297
 298     If the base ring of ``L`` is exact, then ``True`` will be returned if
 299     and only if ``L`` is Lyapunov-like on ``K``.
 300
 301     If the base ring of ``L`` is ``SR``, then the situation is more
 302     complicated:
 303
 304     - ``True`` will be returned if it can be proven that ``L``
 305       is Lyapunov-like on ``K``.
 306     - ``False`` will be returned if it can be proven that ``L``
 307       is not Lyapunov-like on ``K``.
 308     - ``False`` will also be returned if we can't decide; specifically
 309       if we arrive at a symbolic inequality that cannot be resolved.
 310
 311     EXAMPLES:
 312
 313     Diagonal matrices are Lyapunov-like operators on the nonnegative
 314     orthant::
 315
 316         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 317         sage: L = diagonal_matrix(random_vector(QQ,3))
 318         sage: is_lyapunov_like_on(L,K)
 319         True
 320
 321     TESTS:
 322
 323     The identity operator is always Lyapunov-like::
 324
 325         sage: set_random_seed()
 326         sage: K = random_cone(max_ambient_dim=8)
 327         sage: L = identity_matrix(K.lattice_dim())
 328         sage: is_lyapunov_like_on(L,K)
 329         True
 330
 331     The "zero" operator is always Lyapunov-like::
 332
 333         sage: K = random_cone(max_ambient_dim=8)
 334         sage: R = K.lattice().vector_space().base_ring()
 335         sage: L = zero_matrix(R, K.lattice_dim())
 336         sage: is_lyapunov_like_on(L,K)
 337         True
 338
 339     Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
 340     on ``K``::
 341
 342         sage: K = random_cone(max_ambient_dim=5)
 343         sage: all([ is_lyapunov_like_on(L,K)            # long time
 344         ....:       for L in K.lyapunov_like_basis() ]) # long time
 345         True
 346         sage: all([ is_lyapunov_like_on(L.change_ring(SR),K) # long time
 347         ....:       for L in K.lyapunov_like_basis() ])      # long time
 348         True
 349
 350     Technically we could test this, but for now only closed convex cones
 351     are supported as our ``K`` argument::
 352
 353         sage: L = identity_matrix(3)
 354         sage: K = [ vector([2,2,-1]), vector([5,4,-3]) ]
 355         sage: is_lyapunov_like_on(L,K)
 356         Traceback (most recent call last):
 357         ...
 358         TypeError: K must be a Cone.
 359
 360     """
 361     if not is_Cone(K):
 362         raise TypeError('K must be a Cone.')
 363     if not L.base_ring().is_exact() and not L.base_ring() is SR:
 364         raise ValueError('The base ring of L is neither SR nor exact.')
 365
 366     return all([ s*(L*x) == 0
 367                  for (x,s) in K.discrete_complementarity_set() ])
 368
 369
 370 def LL_cone(K):
 371     gens = K.lyapunov_like_basis()
 372     L = ToricLattice(K.lattice_dim()**2)
 373     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 374
 375 def Sigma_cone(K):
 376     gens = K.cross_positive_operators_gens()
 377     L = ToricLattice(K.lattice_dim()**2)
 378     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 379
 380 def Z_cone(K):
 381     gens = K.Z_operators_gens()
 382     L = ToricLattice(K.lattice_dim()**2)
 383     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 384
 385 def pi_cone(K1, K2=None):
 386     if K2 is None:
 387         K2 = K1
 388     gens = K1.positive_operators_gens(K2)
 389     L = ToricLattice(K1.lattice_dim()*K2.lattice_dim())
 390     return Cone([ g.list() for g in gens ], lattice=L, check=False)