X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=test%2Frandomgen.py;fp=test%2Frandomgen.py;h=ef1b19e569825625e134eeef84be3445f2f87b3c;hb=ac39a0b32d176fa78ecd5cf4ef21676e3bd56d6c;hp=0000000000000000000000000000000000000000;hpb=f49f013a90848d4b158dabe85987581c1255164e;p=dunshire.git diff --git a/test/randomgen.py b/test/randomgen.py new file mode 100644 index 0000000..ef1b19e --- /dev/null +++ b/test/randomgen.py @@ -0,0 +1,403 @@ +""" +Random thing generators used in the rest of the test suite. +""" +from random import randint, uniform + +from math import sqrt +from cvxopt import matrix +from dunshire.cones import NonnegativeOrthant, IceCream +from dunshire.games import SymmetricLinearGame +from dunshire.matrices import (append_col, append_row, identity) + +MAX_COND = 250 +""" +The maximum condition number of a randomly-generated game. +""" + +RANDOM_MAX = 10 +""" +When generating random real numbers or integers, this is used as the +largest allowed magnitude. It keeps our condition numbers down and other +properties within reason. +""" + +def random_scalar(): + """ + Generate a random scalar in ``[-RANDOM_MAX, RANDOM_MAX]``. + + Returns + ------- + + float + + Examples + -------- + + >>> abs(random_scalar()) <= RANDOM_MAX + True + + """ + return uniform(-RANDOM_MAX, RANDOM_MAX) + + +def random_nn_scalar(): + """ + Generate a random nonnegative scalar in ``[0, RANDOM_MAX]``. + + Returns + ------- + + float + + Examples + -------- + + >>> 0 <= random_nn_scalar() <= RANDOM_MAX + True + + """ + return abs(random_scalar()) + + +def random_natural(): + """ + Generate a random natural number between ``1 and RANDOM_MAX`` + inclusive. + + Returns + ------- + + int + + Examples + -------- + + >>> 1 <= random_natural() <= RANDOM_MAX + True + + """ + return randint(1, RANDOM_MAX) + + +def random_matrix(dims): + """ + Generate a random square matrix. + + Parameters + ---------- + + dims : int + The number of rows/columns you want in the returned matrix. + + Returns + ------- + + matrix + A new matrix whose entries are random floats chosen uniformly from + the interval [-RANDOM_MAX, RANDOM_MAX]. + + Examples + -------- + + >>> A = random_matrix(3) + >>> A.size + (3, 3) + + """ + return matrix([[random_scalar() + for _ in range(dims)] + for _ in range(dims)]) + + +def random_nonnegative_matrix(dims): + """ + Generate a random square matrix with nonnegative entries. + + Parameters + ---------- + + dims : int + The number of rows/columns you want in the returned matrix. + + Returns + ------- + + matrix + A new matrix whose entries are chosen by :func:`random_nn_scalar`. + + Examples + -------- + + >>> A = random_nonnegative_matrix(3) + >>> A.size + (3, 3) + >>> all([entry >= 0 for entry in A]) + True + + """ + return matrix([[random_nn_scalar() + for _ in range(dims)] + for _ in range(dims)]) + + +def random_diagonal_matrix(dims): + """ + Generate a random square matrix with zero off-diagonal entries. + + These matrices are Lyapunov-like on the nonnegative orthant, as is + fairly easy to see. + + Parameters + ---------- + + dims : int + The number of rows/columns you want in the returned matrix. + + Returns + ------- + + matrix + A new matrix whose diagonal entries are random floats chosen + using func:`random_scalar` and whose off-diagonal entries are + zero. + + Examples + -------- + + >>> A = random_diagonal_matrix(3) + >>> A.size + (3, 3) + >>> A[0,1] == A[0,2] == A[1,0] == A[2,0] == A[1,2] == A[2,1] == 0 + True + + """ + return matrix([[random_scalar()*int(i == j) + for i in range(dims)] + for j in range(dims)]) + + +def random_skew_symmetric_matrix(dims): + """ + Generate a random skew-symmetrix matrix. + + Parameters + ---------- + + dims : int + The number of rows/columns you want in the returned matrix. + + Returns + ------- + + matrix + A new skew-matrix whose strictly above-diagonal entries are + random floats chosen with :func:`random_scalar`. + + Examples + -------- + + >>> A = random_skew_symmetric_matrix(3) + >>> A.size + (3, 3) + + >>> from dunshire.matrices import norm + >>> A = random_skew_symmetric_matrix(random_natural()) + >>> norm(A + A.trans()) < options.ABS_TOL + True + + """ + strict_ut = [[random_scalar()*int(i < j) + for i in range(dims)] + for j in range(dims)] + + strict_ut = matrix(strict_ut, (dims, dims)) + return strict_ut - strict_ut.trans() + + +def random_lyapunov_like_icecream(dims): + r""" + Generate a random matrix Lyapunov-like on the ice-cream cone. + + The form of these matrices is cited in Gowda and Tao + [GowdaTao]_. The scalar ``a`` and the vector ``b`` (using their + notation) are easy to generate. The submatrix ``D`` is a little + trickier, but it can be found noticing that :math:`C + C^{T} = 0` + for a skew-symmetric matrix :math:`C` implying that :math:`C + C^{T} + + \left(2a\right)I = \left(2a\right)I`. Thus we can stick an + :math:`aI` with each of :math:`C,C^{T}` and let those be our + :math:`D,D^{T}`. + + Parameters + ---------- + + dims : int + The dimension of the ice-cream cone (not of the matrix you want!) + on which the returned matrix should be Lyapunov-like. + + Returns + ------- + + matrix + A new matrix, Lyapunov-like on the ice-cream cone in ``dims`` + dimensions, whose free entries are random floats chosen uniformly + from the interval [-RANDOM_MAX, RANDOM_MAX]. + + References + ---------- + + .. [GowdaTao] M. S. Gowda and J. Tao. On the bilinearity rank of a + proper cone and Lyapunov-like transformations. Mathematical + Programming, 147:155-170, 2014. + + Examples + -------- + + >>> L = random_lyapunov_like_icecream(3) + >>> L.size + (3, 3) + >>> x = matrix([1,1,0]) + >>> s = matrix([1,-1,0]) + >>> abs(inner_product(L*x, s)) < options.ABS_TOL + True + + """ + a = matrix([random_scalar()], (1, 1)) + b = matrix([random_scalar() for _ in range(dims-1)], (dims-1, 1)) + D = random_skew_symmetric_matrix(dims-1) + a*identity(dims-1) + row1 = append_col(a, b.trans()) + row2 = append_col(b, D) + return append_row(row1, row2) + + +def random_orthant_game(): + """ + Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a + random game over the nonnegative orthant, and return the + corresponding :class:`SymmetricLinearGame`. + + We keep going until we generate a game with a condition number under + 5000. + """ + ambient_dim = random_natural() + 1 + K = NonnegativeOrthant(ambient_dim) + e1 = [random_nn_scalar() for _ in range(K.dimension())] + e2 = [random_nn_scalar() for _ in range(K.dimension())] + L = random_matrix(K.dimension()) + G = SymmetricLinearGame(L, K, e1, e2) + + if G.condition() <= MAX_COND: + return G + else: + return random_orthant_game() + + +def random_icecream_game(): + """ + Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a + random game over the ice-cream cone, and return the corresponding + :class:`SymmetricLinearGame`. + """ + # Use a minimum dimension of two to avoid divide-by-zero in + # the fudge factor we make up later. + ambient_dim = random_natural() + 1 + K = IceCream(ambient_dim) + e1 = [1] # Set the "height" of e1 to one + e2 = [1] # And the same for e2 + + # If we choose the rest of the components of e1,e2 randomly + # between 0 and 1, then the largest the squared norm of the + # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We + # need to make it less than one (the height of the cone) so + # that the whole thing is in the cone. The norm of the + # non-height part is sqrt(dim(K) - 1), and we can divide by + # twice that. + fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0)) + e1 += [fudge_factor*uniform(0, 1) for _ in range(K.dimension() - 1)] + e2 += [fudge_factor*uniform(0, 1) for _ in range(K.dimension() - 1)] + L = random_matrix(K.dimension()) + G = SymmetricLinearGame(L, K, e1, e2) + + if G.condition() <= MAX_COND: + return G + else: + return random_icecream_game() + + +def random_ll_orthant_game(): + """ + Return a random Lyapunov game over some nonnegative orthant. + """ + G = random_orthant_game() + L = random_diagonal_matrix(G._K.dimension()) + + # Replace the totally-random ``L`` with random Lyapunov-like one. + G = SymmetricLinearGame(L, G._K, G._e1, G._e2) + + while G.condition() > MAX_COND: + # Try again until the condition number is satisfactory. + G = random_orthant_game() + L = random_diagonal_matrix(G._K.dimension()) + G = SymmetricLinearGame(L, G._K, G._e1, G._e2) + + return G + + +def random_ll_icecream_game(): + """ + Return a random Lyapunov game over some ice-cream cone. + """ + G = random_icecream_game() + L = random_lyapunov_like_icecream(G._K.dimension()) + + # Replace the totally-random ``L`` with random Lyapunov-like one. + G = SymmetricLinearGame(L, G._K, G._e1, G._e2) + + while G.condition() > MAX_COND: + # Try again until the condition number is satisfactory. + G = random_icecream_game() + L = random_lyapunov_like_icecream(G._K.dimension()) + G = SymmetricLinearGame(L, G._K, G._e1, G._e2) + + return G + + +def random_positive_orthant_game(): + G = random_orthant_game() + L = random_nonnegative_matrix(G._K.dimension()) + + # Replace the totally-random ``L`` with the random nonnegative one. + G = SymmetricLinearGame(L, G._K, G._e1, G._e2) + + while G.condition() > MAX_COND: + # Try again until the condition number is satisfactory. + G = random_orthant_game() + L = random_nonnegative_matrix(G._K.dimension()) + G = SymmetricLinearGame(L, G._K, G._e1, G._e2) + + return G + + +def random_nn_scaling(G): + alpha = random_nn_scalar() + H = SymmetricLinearGame(alpha*G._L.trans(), G._K, G._e1, G._e2) + + while H.condition() > MAX_COND: + # Loop until the condition number of H doesn't suck. + alpha = random_nn_scalar() + H = SymmetricLinearGame(alpha*G._L.trans(), G._K, G._e1, G._e2) + + return (alpha, H) + +def random_translation(G): + alpha = random_scalar() + tensor_prod = G._e1 * G._e2.trans() + M = G._L + alpha*tensor_prod + + H = SymmetricLinearGame(M.trans(), G._K, G._e1, G._e2) + while H.condition() > MAX_COND: + # Loop until the condition number of H doesn't suck. + alpha = random_scalar() + M = G._L + alpha*tensor_prod + H = SymmetricLinearGame(M.trans(), G._K, G._e1, G._e2) + + return (alpha, H)