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     """
  78     if not is_Cone(K):
  79         raise TypeError('K must be a Cone')
  80     if not L.base_ring().is_exact() and not L.base_ring() is SR:
  81         raise ValueError('base ring of operator L is neither SR nor exact')
  82
  83     if L.base_ring().is_exact():
  84         # This should be way faster than computing the dual and
  85         # checking a bunch of inequalities, but it doesn't work if
  86         # ``L*x`` is symbolic. For example, ``e in Cone([(1,)])``
  87         # is true, but returns ``False``.
  88         return all([ L*x in K for x in K ])
  89     else:
  90         # Fall back to inequality-checking when the entries of ``L``
  91         # might be symbolic.
  92         return all([ s*(L*x) >= 0 for x in K for s in K.dual() ])
  93
  94
  95 def is_cross_positive_on(L,K):
  96     r"""
  97     Determine whether or not ``L`` is cross-positive on ``K``.
  98
  99     We say that ``L`` is cross-positive on a closed convex cone``K`` if
 100     `\left\langle L\left\lparenx\right\rparen,s\right\rangle \ge 0` for
 101     all pairs `\left\langle x,s \right\rangle` in the complementarity
 102     set of ``K``. This property need only be checked for generators of
 103     ``K`` and its dual.
 104
 105     To reliably check whether or not ``L`` is cross-positive, its base
 106     ring must be either exact (for example, the rationals) or ``SR``. An
 107     exact ring is more reliable, but in some cases a matrix whose
 108     entries contain symbolic constants like ``e`` and ``pi`` will work.
 109
 110     INPUT:
 111
 112     - ``L`` -- A matrix over either an exact ring or ``SR``.
 113
 114     - ``K`` -- A polyhedral closed convex cone.
 115
 116     OUTPUT:
 117
 118     If the base ring of ``L`` is exact, then ``True`` will be returned if
 119     and only if ``L`` is cross-positive on ``K``.
 120
 121     If the base ring of ``L`` is ``SR``, then the situation is more
 122     complicated:
 123
 124     - ``True`` will be returned if it can be proven that ``L``
 125       is cross-positive on ``K``.
 126     - ``False`` will be returned if it can be proven that ``L``
 127       is not cross-positive on ``K``.
 128     - ``False`` will also be returned if we can't decide; specifically
 129       if we arrive at a symbolic inequality that cannot be resolved.
 130
 131     EXAMPLES:
 132
 133     The identity operator is always cross-positive::
 134
 135         sage: set_random_seed()
 136         sage: K = random_cone(max_ambient_dim=8)
 137         sage: L = identity_matrix(K.lattice_dim())
 138         sage: is_cross_positive_on(L,K)
 139         True
 140
 141     The "zero" operator is always cross-positive::
 142
 143         sage: K = random_cone(max_ambient_dim=8)
 144         sage: R = K.lattice().vector_space().base_ring()
 145         sage: L = zero_matrix(R, K.lattice_dim())
 146         sage: is_cross_positive_on(L,K)
 147         True
 148
 149     TESTS:
 150
 151     Everything in ``K.cross_positive_operators_gens()`` should be
 152     cross-positive on ``K``::
 153
 154         sage: K = random_cone(max_ambient_dim=5)
 155         sage: all([ is_cross_positive_on(L,K)                     # long time
 156         ....:       for L in K.cross_positive_operators_gens() ]) # long time
 157         True
 158         sage: all([ is_cross_positive_on(L.change_ring(SR),K)     # long time
 159         ....:       for L in K.cross_positive_operators_gens() ]) # long time
 160         True
 161
 162     """
 163     if not is_Cone(K):
 164         raise TypeError('K must be a Cone')
 165     if not L.base_ring().is_exact() and not L.base_ring() is SR:
 166         raise ValueError('base ring of operator L is neither SR nor exact')
 167
 168     return all([ s*(L*x) >= 0
 169                  for (x,s) in K.discrete_complementarity_set() ])
 170
 171 def is_Z_on(L,K):
 172     r"""
 173     Determine whether or not ``L`` is a Z-operator on ``K``.
 174
 175     We say that ``L`` is a Z-operator on a closed convex cone``K`` if
 176     `\left\langle L\left\lparenx\right\rparen,s\right\rangle \le 0` for
 177     all pairs `\left\langle x,s \right\rangle` in the complementarity
 178     set of ``K``. It is known that this property need only be checked
 179     for generators of ``K`` and its dual.
 180
 181     A matrix is a Z-operator on ``K`` if and only if its negation is a
 182     cross-positive operator on ``K``.
 183
 184     To reliably check whether or not ``L`` is a Z operator, its base
 185     ring must be either exact (for example, the rationals) or ``SR``. An
 186     exact ring is more reliable, but in some cases a matrix whose
 187     entries contain symbolic constants like ``e`` and ``pi`` will work.
 188
 189     INPUT:
 190
 191     - ``L`` -- A matrix over either an exact ring or ``SR``.
 192
 193     - ``K`` -- A polyhedral closed convex cone.
 194
 195     OUTPUT:
 196
 197     If the base ring of ``L`` is exact, then ``True`` will be returned if
 198     and only if ``L`` is a Z-operator on ``K``.
 199
 200     If the base ring of ``L`` is ``SR``, then the situation is more
 201     complicated:
 202
 203     - ``True`` will be returned if it can be proven that ``L``
 204       is a Z-operator on ``K``.
 205     - ``False`` will be returned if it can be proven that ``L``
 206       is not a Z-operator on ``K``.
 207     - ``False`` will also be returned if we can't decide; specifically
 208       if we arrive at a symbolic inequality that cannot be resolved.
 209
 210     EXAMPLES:
 211
 212     The identity operator is always a Z-operator::
 213
 214         sage: set_random_seed()
 215         sage: K = random_cone(max_ambient_dim=8)
 216         sage: L = identity_matrix(K.lattice_dim())
 217         sage: is_Z_on(L,K)
 218         True
 219
 220     The "zero" operator is always a Z-operator::
 221
 222         sage: K = random_cone(max_ambient_dim=8)
 223         sage: R = K.lattice().vector_space().base_ring()
 224         sage: L = zero_matrix(R, K.lattice_dim())
 225         sage: is_Z_on(L,K)
 226         True
 227
 228     TESTS:
 229
 230     Everything in ``K.Z_operators_gens()`` should be a Z-operator
 231     on ``K``::
 232
 233         sage: K = random_cone(max_ambient_dim=5)
 234         sage: all([ is_Z_on(L,K)                     # long time
 235         ....:       for L in K.Z_operators_gens() ]) # long time
 236         True
 237         sage: all([ is_Z_on(L.change_ring(SR),K)     # long time
 238         ....:       for L in K.Z_operators_gens() ]) # long time
 239         True
 240
 241     """
 242     return is_cross_positive_on(-L,K)
 243
 244
 245 def is_lyapunov_like_on(L,K):
 246     r"""
 247     Determine whether or not ``L`` is Lyapunov-like on ``K``.
 248
 249     We say that ``L`` is Lyapunov-like on a closed convex cone ``K`` if
 250     `\left\langle L\left\lparenx\right\rparen,s\right\rangle = 0` for
 251     all pairs `\left\langle x,s \right\rangle` in the complementarity
 252     set of ``K``. This property need only be checked for generators of
 253     ``K`` and its dual.
 254
 255     To reliably check whether or not ``L`` is Lyapunov-like, its base
 256     ring must be either exact (for example, the rationals) or ``SR``. An
 257     exact ring is more reliable, but in some cases a matrix whose
 258     entries contain symbolic constants like ``e`` and ``pi`` will work.
 259
 260     INPUT:
 261
 262     - ``L`` -- A matrix over either an exact ring or ``SR``.
 263
 264     - ``K`` -- A polyhedral closed convex cone.
 265
 266     OUTPUT:
 267
 268     If the base ring of ``L`` is exact, then ``True`` will be returned if
 269     and only if ``L`` is Lyapunov-like on ``K``.
 270
 271     If the base ring of ``L`` is ``SR``, then the situation is more
 272     complicated:
 273
 274     - ``True`` will be returned if it can be proven that ``L``
 275       is Lyapunov-like on ``K``.
 276     - ``False`` will be returned if it can be proven that ``L``
 277       is not Lyapunov-like on ``K``.
 278     - ``False`` will also be returned if we can't decide; specifically
 279       if we arrive at a symbolic inequality that cannot be resolved.
 280
 281     EXAMPLES:
 282
 283     Diagonal matrices are Lyapunov-like operators on the nonnegative
 284     orthant::
 285
 286         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 287         sage: L = diagonal_matrix(random_vector(QQ,3))
 288         sage: is_lyapunov_like_on(L,K)
 289         True
 290
 291     TESTS:
 292
 293     The identity operator is always Lyapunov-like::
 294
 295         sage: set_random_seed()
 296         sage: K = random_cone(max_ambient_dim=8)
 297         sage: L = identity_matrix(K.lattice_dim())
 298         sage: is_lyapunov_like_on(L,K)
 299         True
 300
 301     The "zero" operator is always Lyapunov-like::
 302
 303         sage: K = random_cone(max_ambient_dim=8)
 304         sage: R = K.lattice().vector_space().base_ring()
 305         sage: L = zero_matrix(R, K.lattice_dim())
 306         sage: is_lyapunov_like_on(L,K)
 307         True
 308
 309     Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
 310     on ``K``::
 311
 312         sage: K = random_cone(max_ambient_dim=5)
 313         sage: all([ is_lyapunov_like_on(L,K)            # long time
 314         ....:       for L in K.lyapunov_like_basis() ]) # long time
 315         True
 316         sage: all([ is_lyapunov_like_on(L.change_ring(SR),K) # long time
 317         ....:       for L in K.lyapunov_like_basis() ])      # long time
 318         True
 319
 320     """
 321     if not is_Cone(K):
 322         raise TypeError('K must be a Cone')
 323     if not L.base_ring().is_exact() and not L.base_ring() is SR:
 324         raise ValueError('base ring of operator L is neither SR nor exact')
 325
 326     return all([ s*(L*x) == 0
 327                  for (x,s) in K.discrete_complementarity_set() ])
 328
 329
 330 def LL_cone(K):
 331     gens = K.lyapunov_like_basis()
 332     L = ToricLattice(K.lattice_dim()**2)
 333     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 334
 335 def Sigma_cone(K):
 336     gens = K.cross_positive_operators_gens()
 337     L = ToricLattice(K.lattice_dim()**2)
 338     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 339
 340 def Z_cone(K):
 341     gens = K.Z_operators_gens()
 342     L = ToricLattice(K.lattice_dim()**2)
 343     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 344
 345 def pi_cone(K1, K2=None):
 346     if K2 is None:
 347         K2 = K1
 348     gens = K1.positive_operators_gens(K2)
 349     L = ToricLattice(K1.lattice_dim()*K2.lattice_dim())
 350     return Cone([ g.list() for g in gens ], lattice=L, check=False)