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