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(row_count, column_count=None):
  83     """
  84     Generate a random matrix.
  85
  86     Parameters
  87     ----------
  88
  89     row_count : int
  90         The number of rows you want in the returned matrix.
  91
  92     column_count: int
  93         The number of columns you want in the returned matrix (default:
  94         the same as ``row_count``).
  95
  96     Returns
  97     -------
  98
  99     matrix
 100         A new matrix whose entries are random floats chosen uniformly from
 101         the interval [-RANDOM_MAX, RANDOM_MAX].
 102
 103     Examples
 104     --------
 105
 106         >>> A = random_matrix(3)
 107         >>> A.size
 108         (3, 3)
 109
 110         >>> A = random_matrix(3,2)
 111         >>> A.size
 112         (3, 2)
 113
 114     """
 115     if column_count is None:
 116         column_count = row_count
 117
 118     entries = [random_scalar() for _ in range(row_count*column_count)]
 119     return matrix(entries, (row_count, column_count))
 120
 121
 122 def random_nonnegative_matrix(row_count, column_count=None):
 123     """
 124     Generate a random matrix with nonnegative entries.
 125
 126     Parameters
 127     ----------
 128
 129     row_count : int
 130         The number of rows you want in the returned matrix.
 131
 132     column_count : int
 133         The number of columns you want in the returned matrix (default:
 134         the same as ``row_count``).
 135
 136     Returns
 137     -------
 138
 139     matrix
 140         A new matrix whose entries are chosen by :func:`random_nn_scalar`.
 141
 142     Examples
 143     --------
 144
 145         >>> A = random_nonnegative_matrix(3)
 146         >>> A.size
 147         (3, 3)
 148         >>> all([entry >= 0 for entry in A])
 149         True
 150
 151         >>> A = random_nonnegative_matrix(3,2)
 152         >>> A.size
 153         (3, 2)
 154         >>> all([entry >= 0 for entry in A])
 155         True
 156
 157     """
 158     if column_count is None:
 159         column_count = row_count
 160
 161     entries = [random_nn_scalar() for _ in range(row_count*column_count)]
 162     return matrix(entries, (row_count, column_count))
 163
 164
 165 def random_diagonal_matrix(dims):
 166     """
 167     Generate a random square matrix with zero off-diagonal entries.
 168
 169     These matrices are Lyapunov-like on the nonnegative orthant, as is
 170     fairly easy to see.
 171
 172     Parameters
 173     ----------
 174
 175     dims : int
 176         The number of rows/columns you want in the returned matrix.
 177
 178     Returns
 179     -------
 180
 181     matrix
 182         A new matrix whose diagonal entries are random floats chosen
 183         using func:`random_scalar` and whose off-diagonal entries are
 184         zero.
 185
 186     Examples
 187     --------
 188
 189         >>> A = random_diagonal_matrix(3)
 190         >>> A.size
 191         (3, 3)
 192         >>> A[0,1] == A[0,2] == A[1,0] == A[2,0] == A[1,2] == A[2,1] == 0
 193         True
 194
 195     """
 196     return matrix([[random_scalar()*int(i == j)
 197                     for i in range(dims)]
 198                    for j in range(dims)])
 199
 200
 201 def random_skew_symmetric_matrix(dims):
 202     """
 203     Generate a random skew-symmetrix matrix.
 204
 205     Parameters
 206     ----------
 207
 208     dims : int
 209         The number of rows/columns you want in the returned matrix.
 210
 211     Returns
 212     -------
 213
 214     matrix
 215         A new skew-matrix whose strictly above-diagonal entries are
 216         random floats chosen with :func:`random_scalar`.
 217
 218     Examples
 219     --------
 220
 221         >>> A = random_skew_symmetric_matrix(3)
 222         >>> A.size
 223         (3, 3)
 224
 225         >>> from dunshire.options import ABS_TOL
 226         >>> from dunshire.matrices import norm
 227         >>> A = random_skew_symmetric_matrix(random_natural())
 228         >>> norm(A + A.trans()) < ABS_TOL
 229         True
 230
 231     """
 232     strict_ut = [[random_scalar()*int(i < j)
 233                   for i in range(dims)]
 234                  for j in range(dims)]
 235
 236     strict_ut = matrix(strict_ut, (dims, dims))
 237     return strict_ut - strict_ut.trans()
 238
 239
 240 def random_lyapunov_like_icecream(dims):
 241     r"""
 242     Generate a random matrix Lyapunov-like on the ice-cream cone.
 243
 244     The form of these matrices is cited in Gowda and Tao
 245     [GowdaTao]_. The scalar ``a`` and the vector ``b`` (using their
 246     notation) are easy to generate. The submatrix ``D`` is a little
 247     trickier, but it can be found noticing that :math:`C + C^{T} = 0`
 248     for a skew-symmetric matrix :math:`C` implying that :math:`C + C^{T}
 249     + \left(2a\right)I = \left(2a\right)I`. Thus we can stick an
 250     :math:`aI` with each of :math:`C,C^{T}` and let those be our
 251     :math:`D,D^{T}`.
 252
 253     Parameters
 254     ----------
 255
 256     dims : int
 257         The dimension of the ice-cream cone (not of the matrix you want!)
 258         on which the returned matrix should be Lyapunov-like.
 259
 260     Returns
 261     -------
 262
 263     matrix
 264         A new matrix, Lyapunov-like on the ice-cream cone in ``dims``
 265         dimensions, whose free entries are random floats chosen uniformly
 266         from the interval [-RANDOM_MAX, RANDOM_MAX].
 267
 268     References
 269     ----------
 270
 271     .. [GowdaTao] M. S. Gowda and J. Tao. On the bilinearity rank of a
 272        proper cone and Lyapunov-like transformations. Mathematical
 273        Programming, 147:155-170, 2014.
 274
 275     Examples
 276     --------
 277
 278         >>> L = random_lyapunov_like_icecream(3)
 279         >>> L.size
 280         (3, 3)
 281
 282         >>> from dunshire.options import ABS_TOL
 283         >>> from dunshire.matrices import inner_product
 284         >>> x = matrix([1,1,0])
 285         >>> s = matrix([1,-1,0])
 286         >>> abs(inner_product(L*x, s)) < ABS_TOL
 287         True
 288
 289     """
 290     a = matrix([random_scalar()], (1, 1))
 291     b = matrix([random_scalar() for _ in range(dims-1)], (dims-1, 1))
 292     D = random_skew_symmetric_matrix(dims-1) + a*identity(dims-1)
 293     row1 = append_col(a, b.trans())
 294     row2 = append_col(b, D)
 295     return append_row(row1, row2)
 296
 297
 298 def random_orthant_game():
 299     """
 300     Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
 301     random game over the nonnegative orthant, and return the
 302     corresponding :class:`SymmetricLinearGame`.
 303
 304     We keep going until we generate a game with a condition number under
 305     MAX_COND.
 306     """
 307     ambient_dim = random_natural() + 1
 308     K = NonnegativeOrthant(ambient_dim)
 309     e1 = [random_nn_scalar() for _ in range(K.dimension())]
 310     e2 = [random_nn_scalar() for _ in range(K.dimension())]
 311     L = random_matrix(K.dimension())
 312     G = SymmetricLinearGame(L, K, e1, e2)
 313
 314     if G.condition() <= MAX_COND:
 315         return G
 316     else:
 317         return random_orthant_game()
 318
 319
 320 def random_icecream_game():
 321     """
 322     Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
 323     random game over the ice-cream cone, and return the corresponding
 324     :class:`SymmetricLinearGame`.
 325
 326     We keep going until we generate a game with a condition number under
 327     MAX_COND.
 328     """
 329     # Use a minimum dimension of two to avoid divide-by-zero in
 330     # the fudge factor we make up later.
 331     ambient_dim = random_natural() + 1
 332     K = IceCream(ambient_dim)
 333     e1 = [1] # Set the "height" of e1 to one
 334     e2 = [1] # And the same for e2
 335
 336     # If we choose the rest of the components of e1,e2 randomly
 337     # between 0 and 1, then the largest the squared norm of the
 338     # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We
 339     # need to make it less than one (the height of the cone) so
 340     # that the whole thing is in the cone. The norm of the
 341     # non-height part is sqrt(dim(K) - 1), and we can divide by
 342     # twice that.
 343     fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0))
 344     e1 += [fudge_factor*uniform(0, 1) for _ in range(K.dimension() - 1)]
 345     e2 += [fudge_factor*uniform(0, 1) for _ in range(K.dimension() - 1)]
 346     L = random_matrix(K.dimension())
 347     G = SymmetricLinearGame(L, K, e1, e2)
 348
 349     if G.condition() <= MAX_COND:
 350         return G
 351     else:
 352         return random_icecream_game()
 353
 354
 355 def random_ll_orthant_game():
 356     """
 357     Return a random Lyapunov game over some nonnegative orthant.
 358
 359     We first construct a :func:`random_orthant_game` and then modify it
 360     to have a :func:`random_diagonal_matrix` as its operator. Such
 361     things are Lyapunov-like on the nonnegative orthant. That process is
 362     repeated until the condition number of the resulting game is within
 363     ``MAX_COND``.
 364     """
 365     G = random_orthant_game()
 366     L = random_diagonal_matrix(G._K.dimension())
 367
 368     # Replace the totally-random ``L`` with random Lyapunov-like 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_diagonal_matrix(G._K.dimension())
 375         G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 376
 377     return G
 378
 379
 380 def random_ll_icecream_game():
 381     """
 382     Return a random Lyapunov game over some ice-cream cone.
 383
 384     We first construct a :func:`random_icecream_game` and then modify it
 385     to have a :func:`random_lyapunov_like_icecream` operator. That
 386     process is repeated until the condition number of the resulting game
 387     is within ``MAX_COND``.
 388     """
 389     G = random_icecream_game()
 390     L = random_lyapunov_like_icecream(G._K.dimension())
 391
 392     # Replace the totally-random ``L`` with random Lyapunov-like one.
 393     G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 394
 395     while G.condition() > MAX_COND:
 396         # Try again until the condition number is satisfactory.
 397         G = random_icecream_game()
 398         L = random_lyapunov_like_icecream(G._K.dimension())
 399         G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 400
 401     return G
 402
 403
 404 def random_positive_orthant_game():
 405     """
 406     Return a random game over the nonnegative orthant with a positive
 407     operator.
 408
 409     We first construct a :func:`random_orthant_game` and then modify it
 410     to have a :func:`random_nonnegative_matrix` as its operator. That
 411     process is repeated until the condition number of the resulting game
 412     is within ``MAX_COND``.
 413     """
 414
 415     G = random_orthant_game()
 416     L = random_nonnegative_matrix(G._K.dimension())
 417
 418     # Replace the totally-random ``L`` with the random nonnegative one.
 419     G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 420
 421     while G.condition() > MAX_COND:
 422         # Try again until the condition number is satisfactory.
 423         G = random_orthant_game()
 424         L = random_nonnegative_matrix(G._K.dimension())
 425         G = SymmetricLinearGame(L, G._K, G._e1, G._e2)
 426
 427     return G
 428
 429
 430 def random_nn_scaling(G):
 431     """
 432     Scale the given game by a random nonnegative amount.
 433
 434     Parameters
 435     ----------
 436
 437     G : :class:`SymmetricLinearGame`
 438         The game that you would like to scale.
 439
 440     Returns
 441     -------
 442     (float, :class:`SymmetricLinearGame`)
 443         A pair containing the both the scaling factor and the new scaled game.
 444
 445     """
 446     alpha = random_nn_scalar()
 447     H = SymmetricLinearGame(alpha*G._L.trans(), G._K, G._e1, G._e2)
 448
 449     while H.condition() > MAX_COND:
 450         # Loop until the condition number of H doesn't suck.
 451         alpha = random_nn_scalar()
 452         H = SymmetricLinearGame(alpha*G._L.trans(), G._K, G._e1, G._e2)
 453
 454     return (alpha, H)
 455
 456
 457 def random_translation(G):
 458     """
 459     Translate the given game by a random amount.
 460
 461     Parameters
 462     ----------
 463
 464     G : :class:`SymmetricLinearGame`
 465         The game that you would like to translate.
 466
 467     Returns
 468     -------
 469     (float, :class:`SymmetricLinearGame`)
 470         A pair containing the both the translation distance and the new
 471         scaled game.
 472
 473     """
 474     alpha = random_scalar()
 475     tensor_prod = G._e1 * G._e2.trans()
 476     M = G._L + alpha*tensor_prod
 477
 478     H = SymmetricLinearGame(M.trans(), G._K, G._e1, G._e2)
 479     while H.condition() > MAX_COND:
 480         # Loop until the condition number of H doesn't suck.
 481         alpha = random_scalar()
 482         M = G._L + alpha*tensor_prod
 483         H = SymmetricLinearGame(M.trans(), G._K, G._e1, G._e2)
 484
 485     return (alpha, H)