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