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     .. SEEALSO::
  39
  40           :func:`is_cross_positive_on`,
  41           :func:`is_Z_operator_on`,
  42           :func:`is_lyapunov_like_on`
  43
  44     EXAMPLES:
  45
  46     Nonnegative matrices are positive operators on the nonnegative
  47     orthant::
  48
  49         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
  50         sage: L = random_matrix(QQ,3).apply_map(abs)
  51         sage: is_positive_on(L,K)
  52         True
  53
  54     TESTS:
  55
  56     The identity operator is always positive::
  57
  58         sage: set_random_seed()
  59         sage: K = random_cone(max_ambient_dim=8)
  60         sage: L = identity_matrix(K.lattice_dim())
  61         sage: is_positive_on(L,K)
  62         True
  63
  64     The "zero" operator is always positive::
  65
  66         sage: K = random_cone(max_ambient_dim=8)
  67         sage: R = K.lattice().vector_space().base_ring()
  68         sage: L = zero_matrix(R, K.lattice_dim())
  69         sage: is_positive_on(L,K)
  70         True
  71
  72     Everything in ``K.positive_operators_gens()`` should be
  73     positive on ``K``::
  74
  75         sage: K = random_cone(max_ambient_dim=5)
  76         sage: all([ is_positive_on(L,K)                     # long time
  77         ....:       for L in K.positive_operators_gens() ]) # long time
  78         True
  79         sage: all([ is_positive_on(L.change_ring(SR),K)     # long time
  80         ....:       for L in K.positive_operators_gens() ]) # long time
  81         True
  82
  83     Technically we could test this, but for now only closed convex cones
  84     are supported as our ``K`` argument::
  85
  86         sage: K = [ vector([1,2,3]), vector([5,-1,7]) ]
  87         sage: L = identity_matrix(3)
  88         sage: is_positive_on(L,K)
  89         Traceback (most recent call last):
  90         ...
  91         TypeError: K must be a Cone.
  92
  93     We can't give reliable answers over inexact rings::
  94
  95         sage: K = Cone([(1,2,3), (4,5,6)])
  96         sage: L = identity_matrix(RR,3)
  97         sage: is_positive_on(L,K)
  98         Traceback (most recent call last):
  99         ...
 100         ValueError: The base ring of L is neither SR nor exact.
 101
 102     """
 103
 104     if not is_Cone(K):
 105         raise TypeError('K must be a Cone.')
 106     if not L.base_ring().is_exact() and not L.base_ring() is SR:
 107         raise ValueError('The base ring of L is neither SR nor exact.')
 108
 109     if L.base_ring().is_exact():
 110         # This should be way faster than computing the dual and
 111         # checking a bunch of inequalities, but it doesn't work if
 112         # ``L*x`` is symbolic. For example, ``e in Cone([(1,)])``
 113         # is true, but returns ``False``.
 114         return all([ L*x in K for x in K ])
 115     else:
 116         # Fall back to inequality-checking when the entries of ``L``
 117         # might be symbolic.
 118         return all([ s*(L*x) >= 0 for x in K for s in K.dual() ])
 119
 120
 121 def is_cross_positive_on(L,K):
 122     r"""
 123     Determine whether or not ``L`` is cross-positive on ``K``.
 124
 125     We say that ``L`` is cross-positive on a closed convex cone``K`` if
 126     `\left\langle L\left\lparenx\right\rparen,s\right\rangle \ge 0` for
 127     all pairs `\left\langle x,s \right\rangle` in the complementarity
 128     set of ``K``. This property need only be checked for generators of
 129     ``K`` and its dual.
 130
 131     To reliably check whether or not ``L`` is cross-positive, its base
 132     ring must be either exact (for example, the rationals) or ``SR``. An
 133     exact ring is more reliable, but in some cases a matrix whose
 134     entries contain symbolic constants like ``e`` and ``pi`` will work.
 135
 136     INPUT:
 137
 138     - ``L`` -- A matrix over either an exact ring or ``SR``.
 139
 140     - ``K`` -- A polyhedral closed convex cone.
 141
 142     OUTPUT:
 143
 144     If the base ring of ``L`` is exact, then ``True`` will be returned if
 145     and only if ``L`` is cross-positive on ``K``.
 146
 147     If the base ring of ``L`` is ``SR``, then the situation is more
 148     complicated:
 149
 150     - ``True`` will be returned if it can be proven that ``L``
 151       is cross-positive on ``K``.
 152     - ``False`` will be returned if it can be proven that ``L``
 153       is not cross-positive on ``K``.
 154     - ``False`` will also be returned if we can't decide; specifically
 155       if we arrive at a symbolic inequality that cannot be resolved.
 156
 157     .. SEEALSO::
 158
 159           :func:`is_positive_on`,
 160           :func:`is_Z_operator_on`,
 161           :func:`is_lyapunov_like_on`
 162
 163     EXAMPLES:
 164
 165     The identity operator is always cross-positive::
 166
 167         sage: set_random_seed()
 168         sage: K = random_cone(max_ambient_dim=8)
 169         sage: L = identity_matrix(K.lattice_dim())
 170         sage: is_cross_positive_on(L,K)
 171         True
 172
 173     The "zero" operator is always cross-positive::
 174
 175         sage: K = random_cone(max_ambient_dim=8)
 176         sage: R = K.lattice().vector_space().base_ring()
 177         sage: L = zero_matrix(R, K.lattice_dim())
 178         sage: is_cross_positive_on(L,K)
 179         True
 180
 181     TESTS:
 182
 183     Everything in ``K.cross_positive_operators_gens()`` should be
 184     cross-positive on ``K``::
 185
 186         sage: K = random_cone(max_ambient_dim=5)
 187         sage: all([ is_cross_positive_on(L,K)                     # long time
 188         ....:       for L in K.cross_positive_operators_gens() ]) # long time
 189         True
 190         sage: all([ is_cross_positive_on(L.change_ring(SR),K)     # long time
 191         ....:       for L in K.cross_positive_operators_gens() ]) # long time
 192         True
 193
 194     Technically we could test this, but for now only closed convex cones
 195     are supported as our ``K`` argument::
 196
 197         sage: L = identity_matrix(3)
 198         sage: K = [ vector([8,2,-8]), vector([5,-5,7]) ]
 199         sage: is_cross_positive_on(L,K)
 200         Traceback (most recent call last):
 201         ...
 202         TypeError: K must be a Cone.
 203
 204     We can't give reliable answers over inexact rings::
 205
 206         sage: K = Cone([(1,2,3), (4,5,6)])
 207         sage: L = identity_matrix(RR,3)
 208         sage: is_cross_positive_on(L,K)
 209         Traceback (most recent call last):
 210         ...
 211         ValueError: The base ring of L is neither SR nor exact.
 212
 213     """
 214     if not is_Cone(K):
 215         raise TypeError('K must be a Cone.')
 216     if not L.base_ring().is_exact() and not L.base_ring() is SR:
 217         raise ValueError('The base ring of L is neither SR nor exact.')
 218
 219     return all([ s*(L*x) >= 0
 220                  for (x,s) in K.discrete_complementarity_set() ])
 221
 222 def is_Z_operator_on(L,K):
 223     r"""
 224     Determine whether or not ``L`` is a Z-operator on ``K``.
 225
 226     We say that ``L`` is a Z-operator on a closed convex cone``K`` if
 227     `\left\langle L\left\lparenx\right\rparen,s\right\rangle \le 0` for
 228     all pairs `\left\langle x,s \right\rangle` in the complementarity
 229     set of ``K``. It is known that this property need only be checked
 230     for generators of ``K`` and its dual.
 231
 232     A matrix is a Z-operator on ``K`` if and only if its negation is a
 233     cross-positive operator on ``K``.
 234
 235     To reliably check whether or not ``L`` is a Z operator, its base
 236     ring must be either exact (for example, the rationals) or ``SR``. An
 237     exact ring is more reliable, but in some cases a matrix whose
 238     entries contain symbolic constants like ``e`` and ``pi`` will work.
 239
 240     INPUT:
 241
 242     - ``L`` -- A matrix over either an exact ring or ``SR``.
 243
 244     - ``K`` -- A polyhedral closed convex cone.
 245
 246     OUTPUT:
 247
 248     If the base ring of ``L`` is exact, then ``True`` will be returned if
 249     and only if ``L`` is a Z-operator on ``K``.
 250
 251     If the base ring of ``L`` is ``SR``, then the situation is more
 252     complicated:
 253
 254     - ``True`` will be returned if it can be proven that ``L``
 255       is a Z-operator on ``K``.
 256     - ``False`` will be returned if it can be proven that ``L``
 257       is not a Z-operator on ``K``.
 258     - ``False`` will also be returned if we can't decide; specifically
 259       if we arrive at a symbolic inequality that cannot be resolved.
 260
 261     .. SEEALSO::
 262
 263           :func:`is_positive_on`,
 264           :func:`is_cross_positive_on`,
 265           :func:`is_lyapunov_like_on`
 266
 267     EXAMPLES:
 268
 269     The identity operator is always a Z-operator::
 270
 271         sage: set_random_seed()
 272         sage: K = random_cone(max_ambient_dim=8)
 273         sage: L = identity_matrix(K.lattice_dim())
 274         sage: is_Z_operator_on(L,K)
 275         True
 276
 277     The "zero" operator is always a Z-operator::
 278
 279         sage: K = random_cone(max_ambient_dim=8)
 280         sage: R = K.lattice().vector_space().base_ring()
 281         sage: L = zero_matrix(R, K.lattice_dim())
 282         sage: is_Z_operator_on(L,K)
 283         True
 284
 285     TESTS:
 286
 287     Everything in ``K.Z_operators_gens()`` should be a Z-operator
 288     on ``K``::
 289
 290         sage: K = random_cone(max_ambient_dim=5)
 291         sage: all([ is_Z_operator_on(L,K)            # long time
 292         ....:       for L in K.Z_operators_gens() ]) # long time
 293         True
 294         sage: all([ is_Z_operator_on(L.change_ring(SR),K) # long time
 295         ....:       for L in K.Z_operators_gens() ])      # long time
 296         True
 297
 298     Technically we could test this, but for now only closed convex cones
 299     are supported as our ``K`` argument::
 300
 301         sage: L = identity_matrix(3)
 302         sage: K = [ vector([-4,20,3]), vector([1,-5,2]) ]
 303         sage: is_Z_operator_on(L,K)
 304         Traceback (most recent call last):
 305         ...
 306         TypeError: K must be a Cone.
 307
 308
 309     We can't give reliable answers over inexact rings::
 310
 311         sage: K = Cone([(1,2,3), (4,5,6)])
 312         sage: L = identity_matrix(RR,3)
 313         sage: is_Z_operator_on(L,K)
 314         Traceback (most recent call last):
 315         ...
 316         ValueError: The base ring of L is neither SR nor exact.
 317
 318     """
 319     return is_cross_positive_on(-L,K)
 320
 321
 322 def is_lyapunov_like_on(L,K):
 323     r"""
 324     Determine whether or not ``L`` is Lyapunov-like on ``K``.
 325
 326     We say that ``L`` is Lyapunov-like on a closed convex cone ``K`` if
 327     `\left\langle L\left\lparenx\right\rparen,s\right\rangle = 0` for
 328     all pairs `\left\langle x,s \right\rangle` in the complementarity
 329     set of ``K``. This property need only be checked for generators of
 330     ``K`` and its dual.
 331
 332     An operator is Lyapunov-like on ``K`` if and only if both the
 333     operator itself and its negation are cross-positive on ``K``.
 334
 335     To reliably check whether or not ``L`` is Lyapunov-like, its base
 336     ring must be either exact (for example, the rationals) or ``SR``. An
 337     exact ring is more reliable, but in some cases a matrix whose
 338     entries contain symbolic constants like ``e`` and ``pi`` will work.
 339
 340     INPUT:
 341
 342     - ``L`` -- A matrix over either an exact ring or ``SR``.
 343
 344     - ``K`` -- A polyhedral closed convex cone.
 345
 346     OUTPUT:
 347
 348     If the base ring of ``L`` is exact, then ``True`` will be returned if
 349     and only if ``L`` is Lyapunov-like on ``K``.
 350
 351     If the base ring of ``L`` is ``SR``, then the situation is more
 352     complicated:
 353
 354     - ``True`` will be returned if it can be proven that ``L``
 355       is Lyapunov-like on ``K``.
 356     - ``False`` will be returned if it can be proven that ``L``
 357       is not Lyapunov-like on ``K``.
 358     - ``False`` will also be returned if we can't decide; specifically
 359       if we arrive at a symbolic inequality that cannot be resolved.
 360
 361     .. SEEALSO::
 362
 363           :func:`is_positive_on`,
 364           :func:`is_cross_positive_on`,
 365           :func:`is_Z_operator_on`
 366
 367     EXAMPLES:
 368
 369     Diagonal matrices are Lyapunov-like operators on the nonnegative
 370     orthant::
 371
 372         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 373         sage: L = diagonal_matrix(random_vector(QQ,3))
 374         sage: is_lyapunov_like_on(L,K)
 375         True
 376
 377     TESTS:
 378
 379     The identity operator is always Lyapunov-like::
 380
 381         sage: set_random_seed()
 382         sage: K = random_cone(max_ambient_dim=8)
 383         sage: L = identity_matrix(K.lattice_dim())
 384         sage: is_lyapunov_like_on(L,K)
 385         True
 386
 387     The "zero" operator is always Lyapunov-like::
 388
 389         sage: K = random_cone(max_ambient_dim=8)
 390         sage: R = K.lattice().vector_space().base_ring()
 391         sage: L = zero_matrix(R, K.lattice_dim())
 392         sage: is_lyapunov_like_on(L,K)
 393         True
 394
 395     Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
 396     on ``K``::
 397
 398         sage: K = random_cone(max_ambient_dim=5)
 399         sage: all([ is_lyapunov_like_on(L,K)            # long time
 400         ....:       for L in K.lyapunov_like_basis() ]) # long time
 401         True
 402         sage: all([ is_lyapunov_like_on(L.change_ring(SR),K) # long time
 403         ....:       for L in K.lyapunov_like_basis() ])      # long time
 404         True
 405
 406     Technically we could test this, but for now only closed convex cones
 407     are supported as our ``K`` argument::
 408
 409         sage: L = identity_matrix(3)
 410         sage: K = [ vector([2,2,-1]), vector([5,4,-3]) ]
 411         sage: is_lyapunov_like_on(L,K)
 412         Traceback (most recent call last):
 413         ...
 414         TypeError: K must be a Cone.
 415
 416     We can't give reliable answers over inexact rings::
 417
 418         sage: K = Cone([(1,2,3), (4,5,6)])
 419         sage: L = identity_matrix(RR,3)
 420         sage: is_lyapunov_like_on(L,K)
 421         Traceback (most recent call last):
 422         ...
 423         ValueError: The base ring of L is neither SR nor exact.
 424
 425     An operator is Lyapunov-like on a cone if and only if both the
 426     operator and its negation are cross-positive on the cone::
 427
 428         sage: K = random_cone(max_ambient_dim=5)
 429         sage: R = K.lattice().vector_space().base_ring()
 430         sage: L = random_matrix(R, K.lattice_dim())
 431         sage: actual = is_lyapunov_like_on(L,K)          # long time
 432         sage: expected = (is_cross_positive_on(L,K) and  # long time
 433         ....:             is_cross_positive_on(-L,K))    # long time
 434         sage: actual == expected                         # long time
 435         True
 436
 437     """
 438     if not is_Cone(K):
 439         raise TypeError('K must be a Cone.')
 440     if not L.base_ring().is_exact() and not L.base_ring() is SR:
 441         raise ValueError('The base ring of L is neither SR nor exact.')
 442
 443     # Even though ``discrete_complementarity_set`` is a cached method
 444     # of cones, this is faster than calling ``is_cross_positive_on``
 445     # twice: doing so checks twice as many inequalities as the number
 446     # of equalities that we're about to check.
 447     return all([ s*(L*x) == 0
 448                  for (x,s) in K.discrete_complementarity_set() ])
 449
 450
 451 def LL_cone(K):
 452     gens = K.lyapunov_like_basis()
 453     L = ToricLattice(K.lattice_dim()**2)
 454     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 455
 456 def Sigma_cone(K):
 457     gens = K.cross_positive_operators_gens()
 458     L = ToricLattice(K.lattice_dim()**2)
 459     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 460
 461 def Z_cone(K):
 462     gens = K.Z_operators_gens()
 463     L = ToricLattice(K.lattice_dim()**2)
 464     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 465
 466 def pi_cone(K1, K2=None):
 467     if K2 is None:
 468         K2 = K1
 469     gens = K1.positive_operators_gens(K2)
 470     L = ToricLattice(K1.lattice_dim()*K2.lattice_dim())
 471     return Cone([ g.list() for g in gens ], lattice=L, check=False)