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.matrices import norm
 204         >>> A = random_skew_symmetric_matrix(random_natural())
 205         >>> norm(A + A.trans()) < options.ABS_TOL
 206         True
 207
 208     """
 209     strict_ut = [[random_scalar()*int(i < j)
 210                   for i in range(dims)]
 211                  for j in range(dims)]
 212
 213     strict_ut = matrix(strict_ut, (dims, dims))
 214     return strict_ut - strict_ut.trans()
 215
 216
 217 def random_lyapunov_like_icecream(dims):
 218     r"""
 219     Generate a random matrix Lyapunov-like on the ice-cream cone.
 220
 221     The form of these matrices is cited in Gowda and Tao
 222     [GowdaTao]_. The scalar ``a`` and the vector ``b`` (using their
 223     notation) are easy to generate. The submatrix ``D`` is a little
 224     trickier, but it can be found noticing that :math:`C + C^{T} = 0`
 225     for a skew-symmetric matrix :math:`C` implying that :math:`C + C^{T}
 226     + \left(2a\right)I = \left(2a\right)I`. Thus we can stick an
 227     :math:`aI` with each of :math:`C,C^{T}` and let those be our
 228     :math:`D,D^{T}`.
 229
 230     Parameters
 231     ----------
 232
 233     dims : int
 234         The dimension of the ice-cream cone (not of the matrix you want!)
 235         on which the returned matrix should be Lyapunov-like.
 236
 237     Returns
 238     -------
 239
 240     matrix
 241         A new matrix, Lyapunov-like on the ice-cream cone in ``dims``
 242         dimensions, whose free entries are random floats chosen uniformly
 243         from the interval [-RANDOM_MAX, RANDOM_MAX].
 244
 245     References
 246     ----------
 247
 248     .. [GowdaTao] M. S. Gowda and J. Tao. On the bilinearity rank of a
 249        proper cone and Lyapunov-like transformations. Mathematical
 250        Programming, 147:155-170, 2014.
 251
 252     Examples
 253     --------
 254
 255         >>> L = random_lyapunov_like_icecream(3)
 256         >>> L.size
 257         (3, 3)
 258         >>> x = matrix([1,1,0])
 259         >>> s = matrix([1,-1,0])
 260         >>> abs(inner_product(L*x, s)) < options.ABS_TOL
 261         True
 262
 263     """
 264     a = matrix([random_scalar()], (1, 1))
 265     b = matrix([random_scalar() for _ in range(dims-1)], (dims-1, 1))
 266     D = random_skew_symmetric_matrix(dims-1) + a*identity(dims-1)
 267     row1 = append_col(a, b.trans())
 268     row2 = append_col(b, D)
 269     return append_row(row1, row2)
 270
 271
 272 def random_orthant_game():
 273     """
 274     Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
 275     random game over the nonnegative orthant, and return the
 276     corresponding :class:`SymmetricLinearGame`.
 277
 278     We keep going until we generate a game with a condition number under
 279     5000.
 280     """
 281     ambient_dim = random_natural() + 1
 282     K = NonnegativeOrthant(ambient_dim)
 283     e1 = [random_nn_scalar() for _ in range(K.dimension())]
 284     e2 = [random_nn_scalar() for _ in range(K.dimension())]
 285     L = random_matrix(K.dimension())
 286     G = SymmetricLinearGame(L, K, e1, e2)
 287
 288     if G.condition() <= MAX_COND:
 289         return G
 290     else:
 291         return random_orthant_game()
 292
 293
 294 def random_icecream_game():
 295     """
 296     Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
 297     random game over the ice-cream cone, and return the corresponding
 298     :class:`SymmetricLinearGame`.
 299     """
 300     # Use a minimum dimension of two to avoid divide-by-zero in
 301     # the fudge factor we make up later.
 302     ambient_dim = random_natural() + 1
 303     K = IceCream(ambient_dim)
 304     e1 = [1] # Set the "height" of e1 to one
 305     e2 = [1] # And the same for e2
 306
 307     # If we choose the rest of the components of e1,e2 randomly
 308     # between 0 and 1, then the largest the squared norm of the
 309     # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We
 310     # need to make it less than one (the height of the cone) so
 311     # that the whole thing is in the cone. The norm of the
 312     # non-height part is sqrt(dim(K) - 1), and we can divide by
 313     # twice that.
 314     fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0))
 315     e1 += [fudge_factor*uniform(0, 1) for _ in range(K.dimension() - 1)]
 316     e2 += [fudge_factor*uniform(0, 1) for _ in range(K.dimension() - 1)]
 317     L = random_matrix(K.dimension())
 318     G = SymmetricLinearGame(L, K, e1, e2)
 319
 320     if G.condition() <= MAX_COND:
 321         return G
 322     else:
 323         return random_icecream_game()
 324
 325
 326 def random_ll_orthant_game():
 327     """
 328     Return a random Lyapunov game over some nonnegative orthant.
 329     """
 330     G = random_orthant_game()
 331     L = random_diagonal_matrix(G._K.dimension())
 332
 333     # Replace the totally-random ``L`` with random Lyapunov-like one.
 334     G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 335
 336     while G.condition() > MAX_COND:
 337         # Try again until the condition number is satisfactory.
 338         G = random_orthant_game()
 339         L = random_diagonal_matrix(G._K.dimension())
 340         G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 341
 342     return G
 343
 344
 345 def random_ll_icecream_game():
 346     """
 347     Return a random Lyapunov game over some ice-cream cone.
 348     """
 349     G = random_icecream_game()
 350     L = random_lyapunov_like_icecream(G._K.dimension())
 351
 352     # Replace the totally-random ``L`` with random Lyapunov-like one.
 353     G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 354
 355     while G.condition() > MAX_COND:
 356         # Try again until the condition number is satisfactory.
 357         G = random_icecream_game()
 358         L = random_lyapunov_like_icecream(G._K.dimension())
 359         G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 360
 361     return G
 362
 363
 364 def random_positive_orthant_game():
 365     G = random_orthant_game()
 366     L = random_nonnegative_matrix(G._K.dimension())
 367
 368     # Replace the totally-random ``L`` with the random nonnegative one.
 369     G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 370
 371     while G.condition() > MAX_COND:
 372         # Try again until the condition number is satisfactory.
 373         G = random_orthant_game()
 374         L = random_nonnegative_matrix(G._K.dimension())
 375         G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 376
 377     return G
 378
 379
 380 def random_nn_scaling(G):
 381     alpha = random_nn_scalar()
 382     H = SymmetricLinearGame(alpha*G._L.trans(), G._K, G._e1, G._e2)
 383
 384     while H.condition() > MAX_COND:
 385         # Loop until the condition number of H doesn't suck.
 386         alpha = random_nn_scalar()
 387         H = SymmetricLinearGame(alpha*G._L.trans(), G._K, G._e1, G._e2)
 388
 389     return (alpha, H)
 390
 391 def random_translation(G):
 392     alpha = random_scalar()
 393     tensor_prod = G._e1 * G._e2.trans()
 394     M = G._L + alpha*tensor_prod
 395
 396     H = SymmetricLinearGame(M.trans(), G._K, G._e1, G._e2)
 397     while H.condition() > MAX_COND:
 398         # Loop until the condition number of H doesn't suck.
 399         alpha = random_scalar()
 400         M = G._L + alpha*tensor_prod
 401         H = SymmetricLinearGame(M.trans(), G._K, G._e1, G._e2)
 402
 403     return (alpha, H)