]> gitweb.michael.orlitzky.com - dunshire.git/blobdiff - src/dunshire/games.py
Attempt to recover from "unknown" solutions.
[dunshire.git] / src / dunshire / games.py
index e2d29788f2b2612fe0d23e12be07010902c93451..e85fe63ddcc425f8938adf85392fccddc755901d 100644 (file)
@@ -122,23 +122,91 @@ class Solution:
 
 
 class SymmetricLinearGame:
-    """
+    r"""
     A representation of a symmetric linear game.
 
-    The data for a linear game are,
+    The data for a symmetric linear game are,
 
       * A "payoff" operator ``L``.
-      * A cone ``K``.
-      * A point ``e`` in the interior of ``K``.
-      * A point ``f`` in the interior of the dual of ``K``.
-
-    In a symmetric game, the cone ``K`` is be self-dual. We therefore
-    name the two interior points ``e1`` and ``e2`` to indicate that
-    they come from the same cone but are "chosen" by players one and
-    two respectively.
+      * A symmetric cone ``K``.
+      * Two points ``e1`` and ``e2`` in the interior of ``K``.
 
     The ambient space is assumed to be the span of ``K``.
 
+    With those data understood, the game is played as follows. Players
+    one and two choose points :math:`x` and :math:`y` respectively, from
+    their respective strategy sets,
+
+    .. math::
+        \begin{aligned}
+            \Delta_{1}
+            &=
+            \left\{
+              x \in K \ \middle|\ \left\langle x, e_{2} \right\rangle = 1
+            \right\}\\
+            \Delta_{2}
+            &=
+            \left\{
+              y \in K \ \middle|\ \left\langle y, e_{1} \right\rangle = 1
+            \right\}.
+        \end{aligned}
+
+    Afterwards, a "payout" is computed as :math:`\left\langle
+    L\left(x\right), y \right\rangle` and is paid to player one out of
+    player two's pocket. The game is therefore zero sum, and we suppose
+    that player one would like to guarantee himself the largest minimum
+    payout possible. That is, player one wishes to,
+
+    .. math::
+        \begin{aligned}
+            \text{maximize }
+              &\underset{y \in \Delta_{2}}{\min}\left(
+                \left\langle L\left(x\right), y \right\rangle
+              \right)\\
+            \text{subject to } & x \in \Delta_{1}.
+        \end{aligned}
+
+    Player two has the simultaneous goal to,
+
+    .. math::
+        \begin{aligned}
+            \text{minimize }
+              &\underset{x \in \Delta_{1}}{\max}\left(
+                \left\langle L\left(x\right), y \right\rangle
+              \right)\\
+            \text{subject to } & y \in \Delta_{2}.
+        \end{aligned}
+
+    These goals obviously conflict (the game is zero sum), but an
+    existence theorem guarantees at least one optimal min-max solution
+    from which neither player would like to deviate. This class is
+    able to find such a solution.
+
+    Parameters
+    ----------
+
+    L : list of list of float
+        A matrix represented as a list of ROWS. This representation
+        agrees with (for example) SageMath and NumPy, but not with CVXOPT
+        (whose matrix constructor accepts a list of columns).
+
+    K : :class:`SymmetricCone`
+        The symmetric cone instance over which the game is played.
+
+    e1 : iterable float
+        The interior point of ``K`` belonging to player one; it
+        can be of any iterable type having the correct length.
+
+    e2 : iterable float
+        The interior point of ``K`` belonging to player two; it
+        can be of any enumerable type having the correct length.
+
+    Raises
+    ------
+
+    ValueError
+        If either ``e1`` or ``e2`` lie outside of the cone ``K``.
+
     Examples
     --------
 
@@ -161,7 +229,6 @@ class SymmetricLinearGame:
                [ 2]
                [ 3].
 
-
     Lists can (and probably should) be used for every argument::
 
         >>> from cones import NonnegativeOrthant
@@ -242,23 +309,6 @@ class SymmetricLinearGame:
     def __init__(self, L, K, e1, e2):
         """
         Create a new SymmetricLinearGame object.
-
-        INPUT:
-
-          - ``L`` -- an square matrix represented as a list of lists
-             of real numbers. ``L`` itself is interpreted as a list of
-             ROWS, which agrees with (for example) SageMath and NumPy,
-             but not with CVXOPT (whose matrix constructor accepts a
-             list of columns).
-
-          - ``K`` -- a SymmetricCone instance.
-
-          - ``e1`` -- the interior point of ``K`` belonging to player one;
-            it can be of any enumerable type having the correct length.
-
-          - ``e2`` -- the interior point of ``K`` belonging to player two;
-            it can be of any enumerable type having the correct length.
-
         """
         self._K = K
         self._e1 = matrix(e1, (K.dimension(), 1))
@@ -269,10 +319,10 @@ class SymmetricLinearGame:
         # feeding it to CVXOPT.
         self._L = matrix(L, (K.dimension(), K.dimension())).trans()
 
-        if not K.contains_strict(self._e1):
+        if not self._e1 in K:
             raise ValueError('the point e1 must lie in the interior of K')
 
-        if not K.contains_strict(self._e2):
+        if not self._e2 in K:
             raise ValueError('the point e2 must lie in the interior of K')
 
     def __str__(self):
@@ -292,16 +342,20 @@ class SymmetricLinearGame:
 
     def solution(self):
         """
-        Solve this linear game and return a Solution object.
+        Solve this linear game and return a :class:`Solution`.
+
+        Returns
+        -------
 
-        OUTPUT:
+        :class:`Solution`
+            A :class:`Solution` object describing the game's value and
+            the optimal strategies of both players.
 
-        If the cone program associated with this game could be
-        successfully solved, then a Solution object containing the
-        game's value and optimal strategies is returned. If the game
-        could *not* be solved -- which should never happen -- then a
-        GameUnsolvableException is raised. It can be printed to get the
-        raw output from CVXOPT.
+        Raises
+        ------
+        GameUnsolvableException
+            If the game could not be solved (if an optimal solution to its
+            associated cone program was not found).
 
         Examples
         --------
@@ -379,25 +433,44 @@ class SymmetricLinearGame:
         # what happened.
         soln_dict = solvers.conelp(c, G, h, C.cvxopt_dims(), A, b)
 
+        p1_value = -soln_dict['primal objective']
+        p2_value = -soln_dict['dual objective']
+        p1_optimal = soln_dict['x'][1:]
+        p2_optimal = soln_dict['z'][self._K.dimension():]
+
         # The "status" field contains "optimal" if everything went
         # according to plan. Other possible values are "primal
-        # infeasible", "dual infeasible", "unknown", all of which
-        # mean we didn't get a solution. That should never happen,
-        # because by construction our game has a solution, and thus
-        # the cone program should too.
-        if soln_dict['status'] != 'optimal':
+        # infeasible", "dual infeasible", "unknown", all of which mean
+        # we didn't get a solution. The "infeasible" ones are the
+        # worst, since they indicate that CVXOPT is convinced the
+        # problem is infeasible (and that cannot happen).
+        if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
             raise GameUnsolvableException(soln_dict)
-
-        p1_value = soln_dict['x'][0]
-        p1_optimal = soln_dict['x'][1:]
-        p2_optimal = soln_dict['z'][self._K.dimension():]
+        elif soln_dict['status'] == 'unknown':
+            # When we get a status of "unknown", we may still be able
+            # to salvage a solution out of the returned
+            # dictionary. Often this is the result of numerical
+            # difficulty and we can simply check that the primal/dual
+            # objectives match (within a tolerance) and that the
+            # primal/dual optimal solutions are within the cone (to a
+            # tolerance as well).
+            if (abs(p1_value - p2_value) > options.ABS_TOL):
+                raise GameUnsolvableException(soln_dict)
+            if (p1_optimal not in self._K) or (p2_optimal not in self._K):
+                raise GameUnsolvableException(soln_dict)
 
         return Solution(p1_value, p1_optimal, p2_optimal)
 
+
     def dual(self):
-        """
+        r"""
         Return the dual game to this game.
 
+        If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
+        then its dual is :math:`G^{*} =
+        \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
+        is symmetric, :math:`K^{*} = K`.
+
         Examples
         --------
 
@@ -435,6 +508,48 @@ class SymmetricLinearGameTest(TestCase):
     Tests for the SymmetricLinearGame and Solution classes.
     """
 
+    def random_orthant_params(self):
+        """
+        Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
+        random game over the nonnegative orthant.
+        """
+        ambient_dim = randint(1, 10)
+        K = NonnegativeOrthant(ambient_dim)
+        e1 = [uniform(0.1, 10) for idx in range(K.dimension())]
+        e2 = [uniform(0.1, 10) for idx in range(K.dimension())]
+        L = [[uniform(-10, 10) for i in range(K.dimension())]
+             for j in range(K.dimension())]
+        return (L, K, e1, e2)
+
+
+    def random_icecream_params(self):
+        """
+        Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
+        random game over the ice cream cone.
+        """
+        # Use a minimum dimension of two to avoid divide-by-zero in
+        # the fudge factor we make up later.
+        ambient_dim = randint(2, 10)
+        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 idx in range(K.dimension() - 1)]
+        e2 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
+        L = [[uniform(-10, 10) for i in range(K.dimension())]
+             for j in range(K.dimension())]
+
+        return (L, K, e1, e2)
+
+
     def assert_within_tol(self, first, second):
         """
         Test that ``first`` and ``second`` are equal within our default
@@ -463,12 +578,7 @@ class SymmetricLinearGameTest(TestCase):
         optimal solutions should give us the optimal game value when we
         apply the payoff operator to them.
         """
-        ambient_dim = randint(1, 10)
-        K = NonnegativeOrthant(ambient_dim)
-        e1 = [uniform(0.1, 10) for idx in range(K.dimension())]
-        e2 = [uniform(0.1, 10) for idx in range(K.dimension())]
-        L = [[uniform(-10, 10) for i in range(K.dimension())]
-             for j in range(K.dimension())]
+        (L, K, e1, e2) = self.random_orthant_params()
         self.assert_solution_exists(L, K, e1, e2)
 
     def test_solution_exists_ice_cream(self):
@@ -476,24 +586,53 @@ class SymmetricLinearGameTest(TestCase):
         Like :meth:`test_solution_exists_nonnegative_orthant`, except
         over the ice cream cone.
         """
-        # Use a minimum dimension of two to avoid divide-by-zero in
-        # the fudge factor we make up later.
-        ambient_dim = randint(2, 10)
-        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 idx in range(K.dimension() - 1)]
-        e2 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
-        L = [[uniform(-10, 10) for i in range(K.dimension())]
-             for j in range(K.dimension())]
+        (L, K, e1, e2) = self.random_icecream_params()
         self.assert_solution_exists(L, K, e1, e2)
 
+
+    def test_negative_value_Z_operator(self):
+        """
+        Test the example given in Gowda/Ravindran of a Z-matrix with
+        negative game value on the nonnegative orthant.
+        """
+        K = NonnegativeOrthant(2)
+        e1 = [1,1]
+        e2 = e1
+        L = [[1,-2],[-2,1]]
+        G = SymmetricLinearGame(L, K, e1, e2)
+        self.assertTrue(G.solution().game_value() < -options.ABS_TOL)
+
+
+    def test_nonnegative_scaling_orthant(self):
+        """
+        Test that scaling ``L`` by a nonnegative number scales the value
+        of the game by the same number. Use the nonnegative orthant as
+        our cone.
+        """
+        (L, K, e1, e2) = self.random_orthant_params()
+        L = matrix(L) # So that we can scale it by alpha below.
+        G1 = SymmetricLinearGame(L, K, e1, e2)
+        value1 = G1.solution().game_value()
+        alpha = uniform(0.1, 10)
+
+        G2 = SymmetricLinearGame(alpha*L, K, e1, e2)
+        value2 = G2.solution().game_value()
+        self.assert_within_tol(alpha*value1, value2)
+
+
+    def test_nonnegative_scaling_icecream(self):
+        """
+        The same test as :meth:`test_nonnegative_scaling_orthant`,
+        except over the ice cream cone.
+        """
+        (L, K, e1, e2) = self.random_icecream_params()
+        L = matrix(L) # So that we can scale it by alpha below.
+
+        G1 = SymmetricLinearGame(L, K, e1, e2)
+        value1 = G1.solution().game_value()
+        alpha = uniform(0.1, 10)
+
+        G2 = SymmetricLinearGame(alpha*L, K, e1, e2)
+        value2 = G2.solution().game_value()
+        self.assert_within_tol(alpha*value1, value2)
+