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 = 100
  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     Examples
 318     --------
 319
 320         >>> random_orthant_game()
 321         <dunshire.games.SymmetricLinearGame object at 0x...>
 322
 323     """
 324     ambient_dim = random_natural() + 1
 325     K = NonnegativeOrthant(ambient_dim)
 326     e1 = [0.1 + random_nn_scalar() for _ in range(K.dimension())]
 327     e2 = [0.1 + random_nn_scalar() for _ in range(K.dimension())]
 328     L = random_matrix(K.dimension())
 329     G = SymmetricLinearGame(L, K, e1, e2)
 330
 331     if G.condition() <= MAX_COND:
 332         return G
 333     else:
 334         return random_orthant_game()
 335
 336
 337 def random_icecream_game():
 338     """
 339     Generate a random game over the ice-cream cone.
 340
 341     We generate each of ``L``, ``K``, ``e1``, and ``e2`` randomly within
 342     the constraints of the ice-cream cone, and then construct a game
 343     from them. The process is repeated until we generate a game with a
 344     condition number under ``MAX_COND``.
 345
 346     Returns
 347     -------
 348
 349     SymmetricLinearGame
 350         A random game over some ice-cream cone.
 351
 352     Examples
 353     --------
 354
 355         >>> random_icecream_game()
 356         <dunshire.games.SymmetricLinearGame object at 0x...>
 357
 358     """
 359     # Use a minimum dimension of two to avoid divide-by-zero in
 360     # the fudge factor we make up later.
 361     ambient_dim = random_natural() + 1
 362     K = IceCream(ambient_dim)
 363     e1 = [1] # Set the "height" of e1 to one
 364     e2 = [1] # And the same for e2
 365
 366     # If we choose the rest of the components of e1,e2 randomly
 367     # between 0 and 1, then the largest the squared norm of the
 368     # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We
 369     # need to make it less than one (the height of the cone) so
 370     # that the whole thing is in the cone. The norm of the
 371     # non-height part is sqrt(dim(K) - 1), and we can divide by
 372     # twice that.
 373     fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0))
 374     e1 += [fudge_factor*uniform(0, 1) for _ in range(K.dimension() - 1)]
 375     e2 += [fudge_factor*uniform(0, 1) for _ in range(K.dimension() - 1)]
 376     L = random_matrix(K.dimension())
 377     G = SymmetricLinearGame(L, K, e1, e2)
 378
 379     if G.condition() <= MAX_COND:
 380         return G
 381     else:
 382         return random_icecream_game()
 383
 384
 385 def random_ll_orthant_game():
 386     """
 387     Return a random Lyapunov game over some nonnegative orthant.
 388
 389     We first construct a :func:`random_orthant_game` and then modify it
 390     to have a :func:`random_diagonal_matrix` as its operator. Such
 391     things are Lyapunov-like on the nonnegative orthant. That process is
 392     repeated until the condition number of the resulting game is within
 393     ``MAX_COND``.
 394
 395     Returns
 396     -------
 397
 398     SymmetricLinearGame
 399         A random game over some nonnegative orthant whose ``payoff`` method
 400         is based on a Lyapunov-like ``L`` operator.
 401
 402     Examples
 403     --------
 404
 405         >>> random_ll_orthant_game()
 406         <dunshire.games.SymmetricLinearGame object at 0x...>
 407
 408     """
 409     G = random_orthant_game()
 410     L = random_diagonal_matrix(G.dimension())
 411
 412     # Replace the totally-random ``L`` with random Lyapunov-like one.
 413     G = SymmetricLinearGame(L, G.K(), G.e1(), G.e2())
 414
 415     while G.condition() > MAX_COND:
 416         # Try again until the condition number is satisfactory.
 417         G = random_orthant_game()
 418         L = random_diagonal_matrix(G.dimension())
 419         G = SymmetricLinearGame(L, G.K(), G.e1(), G.e2())
 420
 421     return G
 422
 423
 424 def random_ll_icecream_game():
 425     """
 426     Return a random Lyapunov game over some ice-cream cone.
 427
 428     We first construct a :func:`random_icecream_game` and then modify it
 429     to have a :func:`random_lyapunov_like_icecream` operator. That
 430     process is repeated until the condition number of the resulting game
 431     is within ``MAX_COND``.
 432
 433     Returns
 434     -------
 435
 436     SymmetricLinearGame
 437         A random game over some ice-cream cone whose ``payoff`` method
 438         is based on a Lyapunov-like ``L`` operator.
 439
 440     Examples
 441     --------
 442
 443         >>> random_ll_icecream_game()
 444         <dunshire.games.SymmetricLinearGame object at 0x...>
 445
 446     """
 447     G = random_icecream_game()
 448     L = random_lyapunov_like_icecream(G.dimension())
 449
 450     # Replace the totally-random ``L`` with random Lyapunov-like one.
 451     G = SymmetricLinearGame(L, G.K(), G.e1(), G.e2())
 452
 453     while G.condition() > MAX_COND:
 454         # Try again until the condition number is satisfactory.
 455         G = random_icecream_game()
 456         L = random_lyapunov_like_icecream(G.dimension())
 457         G = SymmetricLinearGame(L, G.K(), G.e1(), G.e2())
 458
 459     return G
 460
 461
 462 def random_positive_orthant_game():
 463     """
 464     Return a random game over the nonnegative orthant with a positive
 465     operator.
 466
 467     We first construct a :func:`random_orthant_game` and then modify it
 468     to have a :func:`random_nonnegative_matrix` as its operator. That
 469     process is repeated until the condition number of the resulting game
 470     is within ``MAX_COND``.
 471
 472     Returns
 473     -------
 474
 475     SymmetricLinearGame
 476         A random game over some nonnegative orthant whose ``payoff`` method
 477         is based on a positive ``L`` operator.
 478
 479     Examples
 480     --------
 481
 482         >>> random_positive_orthant_game()
 483         <dunshire.games.SymmetricLinearGame object at 0x...>
 484
 485     """
 486
 487     G = random_orthant_game()
 488     L = random_nonnegative_matrix(G.dimension())
 489
 490     # Replace the totally-random ``L`` with the random nonnegative one.
 491     G = SymmetricLinearGame(L, G.K(), G.e1(), G.e2())
 492
 493     while G.condition() > MAX_COND:
 494         # Try again until the condition number is satisfactory.
 495         G = random_orthant_game()
 496         L = random_nonnegative_matrix(G.dimension())
 497         G = SymmetricLinearGame(L, G.K(), G.e1(), G.e2())
 498
 499     return G
 500
 501
 502 def random_nn_scaling(G):
 503     """
 504     Scale the given game by a random nonnegative amount.
 505
 506     We re-attempt the scaling with a new random number until the
 507     resulting scaled game has an acceptable condition number.
 508
 509     Parameters
 510     ----------
 511
 512     G : SymmetricLinearGame
 513         The game that you would like to scale.
 514
 515     Returns
 516     -------
 517     (float, SymmetricLinearGame)
 518         A pair containing the both the scaling factor and the new scaled game.
 519
 520     Examples
 521     --------
 522
 523         >>> from dunshire.matrices import norm
 524         >>> from dunshire.options import ABS_TOL
 525         >>> G = random_orthant_game()
 526         >>> (alpha, H) = random_nn_scaling(G)
 527         >>> alpha >= 0
 528         True
 529         >>> G.K() == H.K()
 530         True
 531         >>> norm(G.e1() - H.e1()) < ABS_TOL
 532         True
 533         >>> norm(G.e2() - H.e2()) < ABS_TOL
 534         True
 535
 536     """
 537     alpha = random_nn_scalar()
 538     H = SymmetricLinearGame(alpha*G.L().trans(), G.K(), G.e1(), G.e2())
 539
 540     while H.condition() > MAX_COND:
 541         # Loop until the condition number of H doesn't suck.
 542         alpha = random_nn_scalar()
 543         H = SymmetricLinearGame(alpha*G.L().trans(), G.K(), G.e1(), G.e2())
 544
 545     return (alpha, H)
 546
 547
 548 def random_translation(G):
 549     """
 550     Translate the given game by a random amount.
 551
 552     We re-attempt the translation with new random scalars until the
 553     resulting translated game has an acceptable condition number.
 554
 555     Parameters
 556     ----------
 557
 558     G : SymmetricLinearGame
 559         The game that you would like to translate.
 560
 561     Returns
 562     -------
 563     (float, SymmetricLinearGame)
 564         A pair containing the both the translation distance and the new
 565         scaled game.
 566
 567     Examples
 568     --------
 569
 570         >>> from dunshire.matrices import norm
 571         >>> from dunshire.options import ABS_TOL
 572         >>> G = random_orthant_game()
 573         >>> (alpha, H) = random_translation(G)
 574         >>> G.K() == H.K()
 575         True
 576         >>> norm(G.e1() - H.e1()) < ABS_TOL
 577         True
 578         >>> norm(G.e2() - H.e2()) < ABS_TOL
 579         True
 580
 581     """
 582     alpha = random_scalar()
 583     tensor_prod = G.e1() * G.e2().trans()
 584     M = G.L() + alpha*tensor_prod
 585
 586     H = SymmetricLinearGame(M.trans(), G.K(), G.e1(), G.e2())
 587     while H.condition() > MAX_COND:
 588         # Loop until the condition number of H doesn't suck.
 589         alpha = random_scalar()
 590         M = G.L() + alpha*tensor_prod
 591         H = SymmetricLinearGame(M.trans(), G.K(), G.e1(), G.e2())
 592
 593     return (alpha, H)