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