]> gitweb.michael.orlitzky.com - dunshire.git/blobdiff - src/dunshire/games.py
Add a .pylintrc file and fix a few warnings.
[dunshire.git] / src / dunshire / games.py
index 692a0ae550305c3ae56833d500c263498f070889..3d4b09ad8c0c12ff1bba0294a1c5c87ae8ff72bf 100644 (file)
@@ -14,7 +14,8 @@ from unittest import TestCase
 from cvxopt import matrix, printing, solvers
 from cones import CartesianProduct, IceCream, NonnegativeOrthant
 from errors import GameUnsolvableException
 from cvxopt import matrix, printing, solvers
 from cones import CartesianProduct, IceCream, NonnegativeOrthant
 from errors import GameUnsolvableException
-from matrices import append_col, append_row, identity, inner_product, norm
+from matrices import (append_col, append_row, eigenvalues_re, identity,
+                      inner_product, norm)
 import options
 
 printing.options['dformat'] = options.FLOAT_FORMAT
 import options
 
 printing.options['dformat'] = options.FLOAT_FORMAT
@@ -504,13 +505,63 @@ class SymmetricLinearGame:
 
 
 
 
 
 
-def _random_square_matrix(dims):
+def _random_matrix(dims):
     """
     """
-    Generate a random square (``dims``-by-``dims``) matrix,
-    represented as a list of rows. This is used only by the
+    Generate a random square (``dims``-by-``dims``) matrix. This is used
+    only by the :class:`SymmetricLinearGameTest` class.
+    """
+    return matrix([[uniform(-10, 10) for i in range(dims)]
+                   for j in range(dims)])
+
+def _random_nonnegative_matrix(dims):
+    """
+    Generate a random square (``dims``-by-``dims``) matrix with
+    nonnegative entries. This is used only by the
+    :class:`SymmetricLinearGameTest` class.
+    """
+    L = _random_matrix(dims)
+    return matrix([abs(entry) for entry in L], (dims, dims))
+
+def _random_diagonal_matrix(dims):
+    """
+    Generate a random square (``dims``-by-``dims``) matrix with nonzero
+    entries only on the diagonal. This is used only by the
     :class:`SymmetricLinearGameTest` class.
     """
     :class:`SymmetricLinearGameTest` class.
     """
-    return [[uniform(-10, 10) for i in range(dims)] for j in range(dims)]
+    return matrix([[uniform(-10, 10)*int(i == j) for i in range(dims)]
+                   for j in range(dims)])
+
+
+def _random_skew_symmetric_matrix(dims):
+    """
+    Generate a random skew-symmetrix (``dims``-by-``dims``) matrix.
+
+    Examples
+    --------
+
+       >>> A = _random_skew_symmetric_matrix(randint(1, 10))
+       >>> norm(A + A.trans()) < options.ABS_TOL
+       True
+
+    """
+    strict_ut = [[uniform(-10, 10)*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):
+    """
+    Generate a random Lyapunov-like matrix over the ice-cream cone in
+    ``dims`` dimensions.
+    """
+    a = matrix([uniform(-10,10)], (1,1))
+    b = matrix([uniform(-10,10) for idx 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_params():
 
 
 def _random_orthant_params():
@@ -523,8 +574,8 @@ def _random_orthant_params():
     K = NonnegativeOrthant(ambient_dim)
     e1 = [uniform(0.5, 10) for idx in range(K.dimension())]
     e2 = [uniform(0.5, 10) for idx in range(K.dimension())]
     K = NonnegativeOrthant(ambient_dim)
     e1 = [uniform(0.5, 10) for idx in range(K.dimension())]
     e2 = [uniform(0.5, 10) for idx in range(K.dimension())]
-    L = _random_square_matrix(K.dimension())
-    return (L, K, e1, e2)
+    L = _random_matrix(K.dimension())
+    return (L, K, matrix(e1), matrix(e2))
 
 
 def _random_icecream_params():
 
 
 def _random_icecream_params():
@@ -550,9 +601,9 @@ def _random_icecream_params():
     fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0))
     e1 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
     e2 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
     fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0))
     e1 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
     e2 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
-    L = _random_square_matrix(K.dimension())
+    L = _random_matrix(K.dimension())
 
 
-    return (L, K, e1, e2)
+    return (L, K, matrix(e1), matrix(e2))
 
 
 class SymmetricLinearGameTest(TestCase):
 
 
 class SymmetricLinearGameTest(TestCase):
@@ -581,14 +632,13 @@ class SymmetricLinearGameTest(TestCase):
         Given the parameters needed to construct a SymmetricLinearGame,
         ensure that that game has a solution.
         """
         Given the parameters needed to construct a SymmetricLinearGame,
         ensure that that game has a solution.
         """
-        G = SymmetricLinearGame(L, K, e1, e2)
-        soln = G.solution()
-
         # The matrix() constructor assumes that ``L`` is a list of
         # columns, so we transpose it to agree with what
         # SymmetricLinearGame() thinks.
         # The matrix() constructor assumes that ``L`` is a list of
         # columns, so we transpose it to agree with what
         # SymmetricLinearGame() thinks.
-        L_matrix = matrix(L).trans()
-        expected = inner_product(L_matrix*soln.player1_optimal(),
+        G = SymmetricLinearGame(L.trans(), K, e1, e2)
+        soln = G.solution()
+
+        expected = inner_product(L*soln.player1_optimal(),
                                  soln.player2_optimal())
         self.assert_within_tol(soln.game_value(), expected)
 
                                  soln.player2_optimal())
         self.assert_within_tol(soln.game_value(), expected)
 
@@ -632,9 +682,6 @@ class SymmetricLinearGameTest(TestCase):
         Test that scaling ``L`` by a nonnegative number scales the value
         of the game by the same number.
         """
         Test that scaling ``L`` by a nonnegative number scales the value
         of the game by the same number.
         """
-        # Make ``L`` a matrix so that we can scale it by alpha. Its
-        # random, so who cares if it gets transposed.
-        L = matrix(L)
         game1 = SymmetricLinearGame(L, K, e1, e2)
         value1 = game1.solution().game_value()
 
         game1 = SymmetricLinearGame(L, K, e1, e2)
         value1 = game1.solution().game_value()
 
@@ -667,23 +714,19 @@ class SymmetricLinearGameTest(TestCase):
         Check that translating ``L`` by alpha*(e1*e2.trans()) increases
         the value of the associated game by alpha.
         """
         Check that translating ``L`` by alpha*(e1*e2.trans()) increases
         the value of the associated game by alpha.
         """
-        e1 = matrix(e1, (K.dimension(), 1))
-        e2 = matrix(e2, (K.dimension(), 1))
-        game1 = SymmetricLinearGame(L, K, e1, e2)
+        # We need to use ``L`` later, so make sure we transpose it
+        # before passing it in as a column-indexed matrix.
+        game1 = SymmetricLinearGame(L.trans(), K, e1, e2)
         soln1 = game1.solution()
         value1 = soln1.game_value()
         x_bar = soln1.player1_optimal()
         y_bar = soln1.player2_optimal()
 
         soln1 = game1.solution()
         value1 = soln1.game_value()
         x_bar = soln1.player1_optimal()
         y_bar = soln1.player2_optimal()
 
-        # Make ``L`` a CVXOPT matrix so that we can do math with
-        # it. Note that this gives us the "correct" representation of
-        # ``L`` (in agreement with what G has), but COLUMN indexed.
         alpha = uniform(-10, 10)
         alpha = uniform(-10, 10)
-        L = matrix(L).trans()
         tensor_prod = e1*e2.trans()
 
         tensor_prod = e1*e2.trans()
 
-        # Likewise, this is the "correct" representation of ``M``, but
-        # COLUMN indexed...
+        # This is the "correct" representation of ``M``, but COLUMN
+        # indexed...
         M = L + alpha*tensor_prod
 
         # so we have to transpose it when we feed it to the constructor.
         M = L + alpha*tensor_prod
 
         # so we have to transpose it when we feed it to the constructor.
@@ -719,16 +762,11 @@ class SymmetricLinearGameTest(TestCase):
         value that is the negation of the original game. Comes from
         some corollary.
         """
         value that is the negation of the original game. Comes from
         some corollary.
         """
-        e1 = matrix(e1, (K.dimension(), 1))
-        e2 = matrix(e2, (K.dimension(), 1))
-        game1 = SymmetricLinearGame(L, K, e1, e2)
-
-        # Make ``L`` a CVXOPT matrix so that we can do math with
-        # it. Note that this gives us the "correct" representation of
-        # ``L`` (in agreement with what G has), but COLUMN indexed.
-        L = matrix(L).trans()
+        # We need to use ``L`` later, so make sure we transpose it
+        # before passing it in as a column-indexed matrix.
+        game1 = SymmetricLinearGame(L.trans(), K, e1, e2)
 
 
-        # Likewise, this is the "correct" representation of ``M``, but
+        # This is the "correct" representation of ``M``, but
         # COLUMN indexed...
         M = -L.trans()
 
         # COLUMN indexed...
         M = -L.trans()
 
@@ -770,17 +808,14 @@ class SymmetricLinearGameTest(TestCase):
         Two orthogonality relations hold at an optimal solution, and we
         check them here.
         """
         Two orthogonality relations hold at an optimal solution, and we
         check them here.
         """
-        game = SymmetricLinearGame(L, K, e1, e2)
+        # We need to use ``L`` later, so make sure we transpose it
+        # before passing it in as a column-indexed matrix.
+        game = SymmetricLinearGame(L.trans(), K, e1, e2)
         soln = game.solution()
         x_bar = soln.player1_optimal()
         y_bar = soln.player2_optimal()
         value = soln.game_value()
 
         soln = game.solution()
         x_bar = soln.player1_optimal()
         y_bar = soln.player2_optimal()
         value = soln.game_value()
 
-        # Make these matrices so that we can compute with them.
-        L = matrix(L).trans()
-        e1 = matrix(e1, (K.dimension(), 1))
-        e2 = matrix(e2, (K.dimension(), 1))
-
         ip1 = inner_product(y_bar, L*x_bar - value*e1)
         self.assert_within_tol(ip1, 0)
 
         ip1 = inner_product(y_bar, L*x_bar - value*e1)
         self.assert_within_tol(ip1, 0)
 
@@ -814,11 +849,60 @@ class SymmetricLinearGameTest(TestCase):
         This test theoretically applies to the ice-cream cone as well,
         but we don't know how to make positive operators on that cone.
         """
         This test theoretically applies to the ice-cream cone as well,
         but we don't know how to make positive operators on that cone.
         """
-        (L, K, e1, e2) = _random_orthant_params()
+        (_, K, e1, e2) = _random_orthant_params()
 
 
-        # Make the entries of ``L`` nonnegative... this makes it a
-        # positive operator on ``K``.
-        L = [[abs(entry) for entry in row] for row in L]
+        # Ignore that L, we need a nonnegative one.
+        L = _random_nonnegative_matrix(K.dimension())
 
         game = SymmetricLinearGame(L, K, e1, e2)
         self.assertTrue(game.solution().game_value() >= -options.ABS_TOL)
 
         game = SymmetricLinearGame(L, K, e1, e2)
         self.assertTrue(game.solution().game_value() >= -options.ABS_TOL)
+
+
+    def assert_lyapunov_works(self, L, K, e1, e2):
+        """
+        Check that Lyapunov games act the way we expect.
+        """
+        game = SymmetricLinearGame(L, K, e1, e2)
+        soln = game.solution()
+
+        # We only check for positive/negative stability if the game
+        # value is not basically zero. If the value is that close to
+        # zero, we just won't check any assertions.
+        if soln.game_value() > options.ABS_TOL:
+            # L should be positive stable
+            ps = all([eig > -options.ABS_TOL for eig in  eigenvalues_re(L)])
+            self.assertTrue(ps)
+        elif soln.game_value() < -options.ABS_TOL:
+            # L should be negative stable
+            ns = all([eig < options.ABS_TOL for eig in  eigenvalues_re(L)])
+            self.assertTrue(ns)
+
+        # The dual game's value should always equal the primal's.
+        dualsoln = game.dual().solution()
+        self.assert_within_tol(dualsoln.game_value(), soln.game_value())
+
+
+    def test_lyapunov_orthant(self):
+        """
+        Test that a Lyapunov game on the nonnegative orthant works.
+        """
+        (L, K, e1, e2) = _random_orthant_params()
+
+        # Ignore that L, we need a diagonal (Lyapunov-like) one.
+        # (And we don't need to transpose those.)
+        L = _random_diagonal_matrix(K.dimension())
+
+        self.assert_lyapunov_works(L, K, e1, e2)
+
+
+    def test_lyapunov_icecream(self):
+        """
+        Test that a Lyapunov game on the ice-cream cone works.
+        """
+        (L, K, e1, e2) = _random_icecream_params()
+
+        # Ignore that L, we need a diagonal (Lyapunov-like) one.
+        # (And we don't need to transpose those.)
+        L = _random_lyapunov_like_icecream(K.dimension())
+
+        self.assert_lyapunov_works(L, K, e1, e2)