]> gitweb.michael.orlitzky.com - dunshire.git/blobdiff - src/dunshire/games.py
Add tests for Lyapunov games over the ice-cream cone.
[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 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
@@ -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.
     """
-    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():
@@ -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())]
-    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():
@@ -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)]
-    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):
@@ -581,14 +632,13 @@ class SymmetricLinearGameTest(TestCase):
         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.
-        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)
 
@@ -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.
         """
-        # 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()
 
@@ -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.
         """
-        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()
 
-        # 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)
-        L = matrix(L).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.
@@ -719,16 +762,11 @@ class SymmetricLinearGameTest(TestCase):
         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()
 
@@ -770,17 +808,14 @@ class SymmetricLinearGameTest(TestCase):
         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()
 
-        # 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)
 
@@ -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.
         """
-        (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)
+
+
+    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)