test/randomgen.py

   1 """
   2 Random thing generators used in the rest of the test suite.
   3 """
   4 from random import randint, uniform
   5
   6 from math import sqrt
   7 from cvxopt import matrix
   8 from dunshire.cones import NonnegativeOrthant, IceCream
   9 from dunshire.games import SymmetricLinearGame
  10 from dunshire.matrices import (append_col, append_row, identity)
  11
  12 MAX_COND = 250
  13 """
  14 The maximum condition number of a randomly-generated game.
  15 """
  16
  17 RANDOM_MAX = 10
  18 """
  19 When generating random real numbers or integers, this is used as the
  20 largest allowed magnitude. It keeps our condition numbers down and other
  21 properties within reason.
  22 """
  23
  24 def random_scalar():
  25     """
  26     Generate a random scalar in ``[-RANDOM_MAX, RANDOM_MAX]``.
  27
  28     Returns
  29     -------
  30
  31     float
  32
  33     Examples
  34     --------
  35
  36         >>> abs(random_scalar()) <= RANDOM_MAX
  37         True
  38
  39     """
  40     return uniform(-RANDOM_MAX, RANDOM_MAX)
  41
  42
  43 def random_nn_scalar():
  44     """
  45     Generate a random nonnegative scalar in ``[0, RANDOM_MAX]``.
  46
  47     Returns
  48     -------
  49
  50     float
  51
  52     Examples
  53     --------
  54
  55         >>> 0 <= random_nn_scalar() <= RANDOM_MAX
  56         True
  57
  58     """
  59     return abs(random_scalar())
  60
  61
  62 def random_natural():
  63     """
  64     Generate a random natural number between ``1 and RANDOM_MAX``
  65     inclusive.
  66
  67     Returns
  68     -------
  69
  70     int
  71
  72     Examples
  73     --------
  74
  75         >>> 1 <= random_natural() <= RANDOM_MAX
  76         True
  77
  78     """
  79     return randint(1, RANDOM_MAX)
  80
  81
  82 def random_matrix(dims):
  83     """
  84     Generate a random square matrix.
  85
  86     Parameters
  87     ----------
  88
  89     dims : int
  90         The number of rows/columns you want in the returned matrix.
  91
  92     Returns
  93     -------
  94
  95     matrix
  96         A new matrix whose entries are random floats chosen uniformly from
  97         the interval [-RANDOM_MAX, RANDOM_MAX].
  98
  99     Examples
 100     --------
 101
 102         >>> A = random_matrix(3)
 103         >>> A.size
 104         (3, 3)
 105
 106     """
 107     return matrix([[random_scalar()
 108                     for _ in range(dims)]
 109                    for _ in range(dims)])
 110
 111
 112 def random_nonnegative_matrix(dims):
 113     """
 114     Generate a random square matrix with nonnegative entries.
 115
 116     Parameters
 117     ----------
 118
 119     dims : int
 120         The number of rows/columns you want in the returned matrix.
 121
 122     Returns
 123     -------
 124
 125     matrix
 126         A new matrix whose entries are chosen by :func:`random_nn_scalar`.
 127
 128     Examples
 129     --------
 130
 131         >>> A = random_nonnegative_matrix(3)
 132         >>> A.size
 133         (3, 3)
 134         >>> all([entry >= 0 for entry in A])
 135         True
 136
 137     """
 138     return matrix([[random_nn_scalar()
 139                     for _ in range(dims)]
 140                    for _ in range(dims)])
 141
 142
 143 def random_diagonal_matrix(dims):
 144     """
 145     Generate a random square matrix with zero off-diagonal entries.
 146
 147     These matrices are Lyapunov-like on the nonnegative orthant, as is
 148     fairly easy to see.
 149
 150     Parameters
 151     ----------
 152
 153     dims : int
 154         The number of rows/columns you want in the returned matrix.
 155
 156     Returns
 157     -------
 158
 159     matrix
 160         A new matrix whose diagonal entries are random floats chosen
 161         using func:`random_scalar` and whose off-diagonal entries are
 162         zero.
 163
 164     Examples
 165     --------
 166
 167         >>> A = random_diagonal_matrix(3)
 168         >>> A.size
 169         (3, 3)
 170         >>> A[0,1] == A[0,2] == A[1,0] == A[2,0] == A[1,2] == A[2,1] == 0
 171         True
 172
 173     """
 174     return matrix([[random_scalar()*int(i == j)
 175                     for i in range(dims)]
 176                    for j in range(dims)])
 177
 178
 179 def random_skew_symmetric_matrix(dims):
 180     """
 181     Generate a random skew-symmetrix matrix.
 182
 183     Parameters
 184     ----------
 185
 186     dims : int
 187         The number of rows/columns you want in the returned matrix.
 188
 189     Returns
 190     -------
 191
 192     matrix
 193         A new skew-matrix whose strictly above-diagonal entries are
 194         random floats chosen with :func:`random_scalar`.
 195
 196     Examples
 197     --------
 198
 199         >>> A = random_skew_symmetric_matrix(3)
 200         >>> A.size
 201         (3, 3)
 202
 203         >>> from dunshire.options import ABS_TOL
 204         >>> from dunshire.matrices import norm
 205         >>> A = random_skew_symmetric_matrix(random_natural())
 206         >>> norm(A + A.trans()) < ABS_TOL
 207         True
 208
 209     """
 210     strict_ut = [[random_scalar()*int(i < j)
 211                   for i in range(dims)]
 212                  for j in range(dims)]
 213
 214     strict_ut = matrix(strict_ut, (dims, dims))
 215     return strict_ut - strict_ut.trans()
 216
 217
 218 def random_lyapunov_like_icecream(dims):
 219     r"""
 220     Generate a random matrix Lyapunov-like on the ice-cream cone.
 221
 222     The form of these matrices is cited in Gowda and Tao
 223     [GowdaTao]_. The scalar ``a`` and the vector ``b`` (using their
 224     notation) are easy to generate. The submatrix ``D`` is a little
 225     trickier, but it can be found noticing that :math:`C + C^{T} = 0`
 226     for a skew-symmetric matrix :math:`C` implying that :math:`C + C^{T}
 227     + \left(2a\right)I = \left(2a\right)I`. Thus we can stick an
 228     :math:`aI` with each of :math:`C,C^{T}` and let those be our
 229     :math:`D,D^{T}`.
 230
 231     Parameters
 232     ----------
 233
 234     dims : int
 235         The dimension of the ice-cream cone (not of the matrix you want!)
 236         on which the returned matrix should be Lyapunov-like.
 237
 238     Returns
 239     -------
 240
 241     matrix
 242         A new matrix, Lyapunov-like on the ice-cream cone in ``dims``
 243         dimensions, whose free entries are random floats chosen uniformly
 244         from the interval [-RANDOM_MAX, RANDOM_MAX].
 245
 246     References
 247     ----------
 248
 249     .. [GowdaTao] M. S. Gowda and J. Tao. On the bilinearity rank of a
 250        proper cone and Lyapunov-like transformations. Mathematical
 251        Programming, 147:155-170, 2014.
 252
 253     Examples
 254     --------
 255
 256         >>> L = random_lyapunov_like_icecream(3)
 257         >>> L.size
 258         (3, 3)
 259
 260         >>> from dunshire.options import ABS_TOL
 261         >>> from dunshire.matrices import inner_product
 262         >>> x = matrix([1,1,0])
 263         >>> s = matrix([1,-1,0])
 264         >>> abs(inner_product(L*x, s)) < ABS_TOL
 265         True
 266
 267     """
 268     a = matrix([random_scalar()], (1, 1))
 269     b = matrix([random_scalar() for _ in range(dims-1)], (dims-1, 1))
 270     D = random_skew_symmetric_matrix(dims-1) + a*identity(dims-1)
 271     row1 = append_col(a, b.trans())
 272     row2 = append_col(b, D)
 273     return append_row(row1, row2)
 274
 275
 276 def random_orthant_game():
 277     """
 278     Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
 279     random game over the nonnegative orthant, and return the
 280     corresponding :class:`SymmetricLinearGame`.
 281
 282     We keep going until we generate a game with a condition number under
 283     5000.
 284     """
 285     ambient_dim = random_natural() + 1
 286     K = NonnegativeOrthant(ambient_dim)
 287     e1 = [random_nn_scalar() for _ in range(K.dimension())]
 288     e2 = [random_nn_scalar() for _ in range(K.dimension())]
 289     L = random_matrix(K.dimension())
 290     G = SymmetricLinearGame(L, K, e1, e2)
 291
 292     if G.condition() <= MAX_COND:
 293         return G
 294     else:
 295         return random_orthant_game()
 296
 297
 298 def random_icecream_game():
 299     """
 300     Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
 301     random game over the ice-cream cone, and return the corresponding
 302     :class:`SymmetricLinearGame`.
 303     """
 304     # Use a minimum dimension of two to avoid divide-by-zero in
 305     # the fudge factor we make up later.
 306     ambient_dim = random_natural() + 1
 307     K = IceCream(ambient_dim)
 308     e1 = [1] # Set the "height" of e1 to one
 309     e2 = [1] # And the same for e2
 310
 311     # If we choose the rest of the components of e1,e2 randomly
 312     # between 0 and 1, then the largest the squared norm of the
 313     # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We
 314     # need to make it less than one (the height of the cone) so
 315     # that the whole thing is in the cone. The norm of the
 316     # non-height part is sqrt(dim(K) - 1), and we can divide by
 317     # twice that.
 318     fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0))
 319     e1 += [fudge_factor*uniform(0, 1) for _ in range(K.dimension() - 1)]
 320     e2 += [fudge_factor*uniform(0, 1) for _ in range(K.dimension() - 1)]
 321     L = random_matrix(K.dimension())
 322     G = SymmetricLinearGame(L, K, e1, e2)
 323
 324     if G.condition() <= MAX_COND:
 325         return G
 326     else:
 327         return random_icecream_game()
 328
 329
 330 def random_ll_orthant_game():
 331     """
 332     Return a random Lyapunov game over some nonnegative orthant.
 333     """
 334     G = random_orthant_game()
 335     L = random_diagonal_matrix(G._K.dimension())
 336
 337     # Replace the totally-random ``L`` with random Lyapunov-like one.
 338     G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 339
 340     while G.condition() > MAX_COND:
 341         # Try again until the condition number is satisfactory.
 342         G = random_orthant_game()
 343         L = random_diagonal_matrix(G._K.dimension())
 344         G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 345
 346     return G
 347
 348
 349 def random_ll_icecream_game():
 350     """
 351     Return a random Lyapunov game over some ice-cream cone.
 352     """
 353     G = random_icecream_game()
 354     L = random_lyapunov_like_icecream(G._K.dimension())
 355
 356     # Replace the totally-random ``L`` with random Lyapunov-like one.
 357     G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 358
 359     while G.condition() > MAX_COND:
 360         # Try again until the condition number is satisfactory.
 361         G = random_icecream_game()
 362         L = random_lyapunov_like_icecream(G._K.dimension())
 363         G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 364
 365     return G
 366
 367
 368 def random_positive_orthant_game():
 369     G = random_orthant_game()
 370     L = random_nonnegative_matrix(G._K.dimension())
 371
 372     # Replace the totally-random ``L`` with the random nonnegative one.
 373     G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 374
 375     while G.condition() > MAX_COND:
 376         # Try again until the condition number is satisfactory.
 377         G = random_orthant_game()
 378         L = random_nonnegative_matrix(G._K.dimension())
 379         G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 380
 381     return G
 382
 383
 384 def random_nn_scaling(G):
 385     alpha = random_nn_scalar()
 386     H = SymmetricLinearGame(alpha*G._L.trans(), G._K, G._e1, G._e2)
 387
 388     while H.condition() > MAX_COND:
 389         # Loop until the condition number of H doesn't suck.
 390         alpha = random_nn_scalar()
 391         H = SymmetricLinearGame(alpha*G._L.trans(), G._K, G._e1, G._e2)
 392
 393     return (alpha, H)
 394
 395 def random_translation(G):
 396     alpha = random_scalar()
 397     tensor_prod = G._e1 * G._e2.trans()
 398     M = G._L + alpha*tensor_prod
 399
 400     H = SymmetricLinearGame(M.trans(), G._K, G._e1, G._e2)
 401     while H.condition() > MAX_COND:
 402         # Loop until the condition number of H doesn't suck.
 403         alpha = random_scalar()
 404         M = G._L + alpha*tensor_prod
 405         H = SymmetricLinearGame(M.trans(), G._K, G._e1, G._e2)
 406
 407     return (alpha, H)