]> gitweb.michael.orlitzky.com - dunshire.git/blobdiff - dunshire/games.py
Add some more docs for the player[12]_start() methods.
[dunshire.git] / dunshire / games.py
index ae1426a2c611f7e315f94fc5fea0e98f1da0905b..6ab420d1b0c26cc6de561caf85629bae28f0462a 100644 (file)
@@ -9,9 +9,10 @@ from .cones import CartesianProduct
 from .errors import GameUnsolvableException, PoorScalingException
 from .matrices import (append_col, append_row, condition_number, identity,
                        inner_product, norm, specnorm)
-from . import options
+from .options import ABS_TOL, FLOAT_FORMAT, DEBUG_FLOAT_FORMAT
+
+printing.options['dformat'] = FLOAT_FORMAT
 
-printing.options['dformat'] = options.FLOAT_FORMAT
 
 class Solution:
     """
@@ -219,8 +220,7 @@ class SymmetricLinearGame:
                [ 1],
           e2 = [ 1]
                [ 2]
-               [ 3],
-          Condition((L, K, e1, e2)) = 31.834...
+               [ 3]
 
     Lists can (and probably should) be used for every argument::
 
@@ -238,8 +238,7 @@ class SymmetricLinearGame:
           e1 = [ 1]
                [ 1],
           e2 = [ 1]
-               [ 1],
-          Condition((L, K, e1, e2)) = 1.707...
+               [ 1]
 
     The points ``e1`` and ``e2`` can also be passed as some other
     enumerable type (of the correct length) without much harm, since
@@ -261,8 +260,7 @@ class SymmetricLinearGame:
           e1 = [ 1]
                [ 1],
           e2 = [ 1]
-               [ 1],
-          Condition((L, K, e1, e2)) = 1.707...
+               [ 1]
 
     However, ``L`` will always be intepreted as a list of rows, even
     if it is passed as a :class:`cvxopt.base.matrix` which is
@@ -283,8 +281,7 @@ class SymmetricLinearGame:
           e1 = [ 1]
                [ 1],
           e2 = [ 1]
-               [ 1],
-          Condition((L, K, e1, e2)) = 6.073...
+               [ 1]
         >>> L = cvxopt.matrix(L)
         >>> print(L)
         [ 1  3]
@@ -299,8 +296,7 @@ class SymmetricLinearGame:
           e1 = [ 1]
                [ 1],
           e2 = [ 1]
-               [ 1],
-          Condition((L, K, e1, e2)) = 6.073...
+               [ 1]
 
     """
     def __init__(self, L, K, e1, e2):
@@ -334,8 +330,7 @@ class SymmetricLinearGame:
               '  L = {:s},\n' \
               '  K = {!s},\n' \
               '  e1 = {:s},\n' \
-              '  e2 = {:s},\n' \
-              '  Condition((L, K, e1, e2)) = {:f}.'
+              '  e2 = {:s}'
         indented_L = '\n      '.join(str(self.L()).splitlines())
         indented_e1 = '\n       '.join(str(self.e1()).splitlines())
         indented_e2 = '\n       '.join(str(self.e2()).splitlines())
@@ -343,8 +338,7 @@ class SymmetricLinearGame:
         return tpl.format(indented_L,
                           str(self.K()),
                           indented_e1,
-                          indented_e2,
-                          self.condition())
+                          indented_e2)
 
 
     def L(self):
@@ -619,7 +613,7 @@ class SymmetricLinearGame:
 
 
 
-    def _G(self):
+    def G(self):
         r"""
         Return the matrix ``G`` used in our CVXOPT construction.
 
@@ -646,7 +640,7 @@ class SymmetricLinearGame:
             >>> e1 = [1,2,3]
             >>> e2 = [1,1,1]
             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
-            >>> print(SLG._G())
+            >>> print(SLG.G())
             [  0.0000000  -1.0000000   0.0000000   0.0000000]
             [  0.0000000   0.0000000  -1.0000000   0.0000000]
             [  0.0000000   0.0000000   0.0000000  -1.0000000]
@@ -661,7 +655,7 @@ class SymmetricLinearGame:
                           append_col(self.e1(), -self.L()))
 
 
-    def _c(self):
+    def c(self):
         """
         Return the vector ``c`` used in our CVXOPT construction.
 
@@ -688,7 +682,7 @@ class SymmetricLinearGame:
             >>> e1 = [1,2,3]
             >>> e2 = [1,1,1]
             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
-            >>> print(SLG._c())
+            >>> print(SLG.c())
             [-1.0000000]
             [ 0.0000000]
             [ 0.0000000]
@@ -729,8 +723,8 @@ class SymmetricLinearGame:
         """
         return CartesianProduct(self._K, self._K)
 
-    def _h(self):
-        """
+    def h(self):
+        r"""
         Return the ``h`` vector used in our CVXOPT construction.
 
         The ``h`` vector appears on the right-hand side of :math:`Gx + s
@@ -756,7 +750,7 @@ class SymmetricLinearGame:
             >>> e1 = [1,2,3]
             >>> e2 = [1,1,1]
             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
-            >>> print(SLG._h())
+            >>> print(SLG.h())
             [0.0000000]
             [0.0000000]
             [0.0000000]
@@ -772,7 +766,7 @@ class SymmetricLinearGame:
 
     @staticmethod
     def b():
-        """
+        r"""
         Return the ``b`` vector used in our CVXOPT construction.
 
         The vector ``b`` appears on the right-hand side of :math:`Ax =
@@ -815,17 +809,32 @@ class SymmetricLinearGame:
         Return a feasible starting point for player one.
 
         This starting point is for the CVXOPT formulation and not for
-        the original game. The basic premise is that if you normalize
-        :meth:`e2`, then you get a point in :meth:`K` that makes a unit
-        inner product with :meth:`e2`. We then get to choose the primal
-        objective function value such that the constraint involving
-        :meth:`L` is satisfied.
+        the original game. The basic premise is that if you scale
+        :meth:`e2` by the reciprocal of its squared norm, then you get a
+        point in :meth:`K` that makes a unit inner product with
+        :meth:`e2`. We then get to choose the primal objective function
+        value such that the constraint involving :meth:`L` is satisfied.
+
+        Returns
+        -------
+
+        dict
+            A dictionary with two keys, 'x' and 's', which contain the
+            vectors of the same name in the CVXOPT primal problem
+            formulation.
+
+            The vector ``x`` consists of the primal objective function
+            value concatenated with the strategy (for player one) that
+            achieves it. The vector ``s`` is essentially a dummy
+            variable, and is computed from the equality constraing in
+            the CVXOPT primal problem.
+
         """
         p = self.e2() / (norm(self.e2()) ** 2)
         dist = self.K().ball_radius(self.e1())
         nu = - self._L_specnorm()/(dist*norm(self.e2()))
         x = matrix([nu, p], (self.dimension() + 1, 1))
-        s = - self._G()*x
+        s = - self.G()*x
 
         return {'x': x, 's': s}
 
@@ -833,6 +842,29 @@ class SymmetricLinearGame:
     def player2_start(self):
         """
         Return a feasible starting point for player two.
+
+        This starting point is for the CVXOPT formulation and not for
+        the original game. The basic premise is that if you scale
+        :meth:`e1` by the reciprocal of its squared norm, then you get a
+        point in :meth:`K` that makes a unit inner product with
+        :meth:`e1`. We then get to choose the dual objective function
+        value such that the constraint involving :meth:`L` is satisfied.
+
+        Returns
+        -------
+
+        dict
+            A dictionary with two keys, 'y' and 'z', which contain the
+            vectors of the same name in the CVXOPT dual problem
+            formulation.
+
+            The ``1``-by-``1`` vector ``y`` consists of the dual
+            objective function value. The last :meth:`dimension` entries
+            of the vector ``z`` contain the strategy (for player two)
+            that achieves it. The remaining entries of ``z`` are
+            essentially dummy variables, computed from the equality
+            constraint in the CVXOPT dual problem.
+
         """
         q = self.e1() / (norm(self.e1()) ** 2)
         dist = self.K().ball_radius(self.e2())
@@ -847,18 +879,81 @@ class SymmetricLinearGame:
 
     def _L_specnorm(self):
         """
-        Compute the spectral norm of ``L`` and cache it.
+        Compute the spectral norm of :meth:`L` and cache it.
+
+        The spectral norm of the matrix :meth:`L` is used in a few
+        places. Since it can be expensive to compute, we want to cache
+        its value. That is not possible in :func:`specnorm`, which lies
+        outside of a class, so this is the place to do it.
+
+        Returns
+        -------
+
+        float
+            A nonnegative real number; the largest singular value of
+            the matrix :meth:`L`.
+
         """
         if self._L_specnorm_value is None:
             self._L_specnorm_value = specnorm(self.L())
         return self._L_specnorm_value
 
-    def epsilon_scale(self, solution):
-        # Don't return anything smaller than 1... we can't go below
-        # out "minimum tolerance."
+
+    def tolerance_scale(self, solution):
+        r"""
+        Return a scaling factor that should be applied to ``ABS_TOL``
+        for this game.
+
+        When performing certain comparisons, the default tolernace
+        ``ABS_TOL`` may not be appropriate. For example, if we expect
+        ``x`` and ``y`` to be within ``ABS_TOL`` of each other, than the
+        inner product of ``L*x`` and ``y`` can be as far apart as the
+        spectral norm of ``L`` times the sum of the norms of ``x`` and
+        ``y``. Such a comparison is made in :meth:`solution`, and in
+        many of our unit tests.
+
+        The returned scaling factor found from the inner product mentioned
+        above is
+
+        .. math::
+
+            \left\lVert L \right\rVert_{2}
+            \left( \left\lVert \bar{x} \right\rVert
+                   + \left\lVert \bar{y} \right\rVert
+            \right),
+
+        where :math:`\bar{x}` and :math:`\bar{y}` are optimal solutions
+        for players one and two respectively. This scaling factor is not
+        formally justified, but attempting anything smaller leads to
+        test failures.
+
+        .. warning::
+
+            Optimal solutions are not unique, so the scaling factor
+            obtained from ``solution`` may not work when comparing other
+            solutions.
+
+        Parameters
+        ----------
+
+        solution : Solution
+            A solution of this game, used to obtain the norms of the
+            optimal strategies.
+
+        Returns
+        -------
+
+        float
+            A scaling factor to be multiplied by ``ABS_TOL`` when
+            making comparisons involving solutions of this game.
+
+        """
         norm_p1_opt = norm(solution.player1_optimal())
         norm_p2_opt = norm(solution.player2_optimal())
         scale = self._L_specnorm()*(norm_p1_opt + norm_p2_opt)
+
+        # Don't return anything smaller than 1... we can't go below
+        # out "minimum tolerance."
         return max(1, scale)
 
 
@@ -976,23 +1071,62 @@ class SymmetricLinearGame:
               [2.506...]
               [0.000...]
 
+        This is another one that was difficult numerically, and caused
+        trouble even after we fixed the first two::
+
+            >>> from dunshire import *
+            >>> L = [[57.22233908627052301199, 41.70631373437460354126],
+            ...      [83.04512571985074487202, 57.82581810406928468637]]
+            >>> K = NonnegativeOrthant(2)
+            >>> e1 = [7.31887017043399268346, 0.89744171905822367474]
+            >>> e2 = [0.11099824781179848388, 6.12564670639315345113]
+            >>> SLG = SymmetricLinearGame(L,K,e1,e2)
+            >>> print(SLG.solution())
+            Game value: 70.437...
+            Player 1 optimal:
+              [9.009...]
+              [0.000...]
+            Player 2 optimal:
+              [0.136...]
+              [0.000...]
+
+        And finally, here's one that returns an "optimal" solution, but
+        whose primal/dual objective function values are far apart::
+
+            >>> from dunshire import *
+            >>> L = [[ 6.49260076597376212248, -0.60528030227678542019],
+            ...      [ 2.59896077096751731972, -0.97685530240286766457]]
+            >>> K = IceCream(2)
+            >>> e1 = [1, 0.43749513972645248661]
+            >>> e2 = [1, 0.46008379832200291260]
+            >>> SLG = SymmetricLinearGame(L, K, e1, e2)
+            >>> print(SLG.solution())
+            Game value: 11.596...
+            Player 1 optimal:
+              [ 1.852...]
+              [-1.852...]
+            Player 2 optimal:
+              [ 1.777...]
+              [-1.777...]
+
         """
         try:
             opts = {'show_progress': False}
-            soln_dict = solvers.conelp(self._c(),
-                                       self._G(),
-                                       self._h(),
+            soln_dict = solvers.conelp(self.c(),
+                                       self.G(),
+                                       self.h(),
                                        self.C().cvxopt_dims(),
                                        self.A(),
                                        self.b(),
                                        primalstart=self.player1_start(),
+                                       dualstart=self.player2_start(),
                                        options=opts)
         except ValueError as error:
             if str(error) == 'math domain error':
                 # Oops, CVXOPT tried to take the square root of a
                 # negative number. Report some details about the game
                 # rather than just the underlying CVXOPT crash.
-                printing.options['dformat'] = options.DEBUG_FLOAT_FORMAT
+                printing.options['dformat'] = DEBUG_FLOAT_FORMAT
                 raise PoorScalingException(self)
             else:
                 raise error
@@ -1017,9 +1151,23 @@ class SymmetricLinearGame:
         # that CVXOPT is convinced the problem is infeasible (and that
         # cannot happen).
         if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
-            printing.options['dformat'] = options.DEBUG_FLOAT_FORMAT
+            printing.options['dformat'] = DEBUG_FLOAT_FORMAT
             raise GameUnsolvableException(self, soln_dict)
 
+        # For the game value, we could use any of:
+        #
+        #   * p1_value
+        #   * p2_value
+        #   * (p1_value + p2_value)/2
+        #   * the game payoff
+        #
+        # We want the game value to be the payoff, however, so it
+        # makes the most sense to just use that, even if it means we
+        # can't test the fact that p1_value/p2_value are close to the
+        # payoff.
+        payoff = self.payoff(p1_optimal, p2_optimal)
+        soln = Solution(payoff, p1_optimal, p2_optimal)
+
         # The "optimal" and "unknown" results, we actually treat the
         # same. Even if CVXOPT bails out due to numerical difficulty,
         # it will have some candidate points in mind. If those
@@ -1030,29 +1178,18 @@ class SymmetricLinearGame:
         # close enough (one could be low by ABS_TOL, the other high by
         # it) because otherwise CVXOPT might return "unknown" and give
         # us two points in the cone that are nowhere near optimal.
-        if abs(p1_value - p2_value) > 2*options.ABS_TOL:
-            printing.options['dformat'] = options.DEBUG_FLOAT_FORMAT
+        #
+        if abs(p1_value - p2_value) > self.tolerance_scale(soln)*ABS_TOL:
+            printing.options['dformat'] = DEBUG_FLOAT_FORMAT
             raise GameUnsolvableException(self, soln_dict)
 
         # And we also check that the points it gave us belong to the
         # cone, just in case...
         if (p1_optimal not in self._K) or (p2_optimal not in self._K):
-            printing.options['dformat'] = options.DEBUG_FLOAT_FORMAT
+            printing.options['dformat'] = DEBUG_FLOAT_FORMAT
             raise GameUnsolvableException(self, soln_dict)
 
-        # For the game value, we could use any of:
-        #
-        #   * p1_value
-        #   * p2_value
-        #   * (p1_value + p2_value)/2
-        #   * the game payoff
-        #
-        # We want the game value to be the payoff, however, so it
-        # makes the most sense to just use that, even if it means we
-        # can't test the fact that p1_value/p2_value are close to the
-        # payoff.
-        payoff = self.payoff(p1_optimal, p2_optimal)
-        return Solution(payoff, p1_optimal, p2_optimal)
+        return soln
 
 
     def condition(self):
@@ -1083,13 +1220,11 @@ class SymmetricLinearGame:
         >>> e1 = [1]
         >>> e2 = e1
         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
-        >>> actual = SLG.condition()
-        >>> expected = 1.8090169943749477
-        >>> abs(actual - expected) < options.ABS_TOL
-        True
+        >>> SLG.condition()
+        1.809...
 
         """
-        return (condition_number(self._G()) + condition_number(self.A()))/2
+        return (condition_number(self.G()) + condition_number(self.A()))/2
 
 
     def dual(self):
@@ -1121,8 +1256,7 @@ class SymmetricLinearGame:
                    [ 3],
               e2 = [ 1]
                    [ 1]
-                   [ 1],
-              Condition((L, K, e1, e2)) = 44.476...
+                   [ 1]
 
         """
         # We pass ``self.L()`` right back into the constructor, because