src/dunshire/games.py

   1 """
   2 Symmetric linear games and their solutions.
   3
   4 This module contains the main :class:`SymmetricLinearGame` class that
   5 knows how to solve a linear game.
   6 """
   7
   8 # These few are used only for tests.
   9 from math import sqrt
  10 from random import randint, uniform
  11 from unittest import TestCase
  12
  13 # These are mostly actually needed.
  14 from cvxopt import matrix, printing, solvers
  15 from cones import CartesianProduct, IceCream, NonnegativeOrthant
  16 from errors import GameUnsolvableException
  17 from matrices import append_col, append_row, identity, inner_product
  18 import options
  19
  20 printing.options['dformat'] = options.FLOAT_FORMAT
  21 solvers.options['show_progress'] = options.VERBOSE
  22
  23
  24 class Solution:
  25     """
  26     A representation of the solution of a linear game. It should contain
  27     the value of the game, and both players' strategies.
  28
  29     Examples
  30     --------
  31
  32         >>> print(Solution(10, matrix([1,2]), matrix([3,4])))
  33         Game value: 10.0000000
  34         Player 1 optimal:
  35           [ 1]
  36           [ 2]
  37         Player 2 optimal:
  38           [ 3]
  39           [ 4]
  40
  41     """
  42     def __init__(self, game_value, p1_optimal, p2_optimal):
  43         """
  44         Create a new Solution object from a game value and two optimal
  45         strategies for the players.
  46         """
  47         self._game_value = game_value
  48         self._player1_optimal = p1_optimal
  49         self._player2_optimal = p2_optimal
  50
  51     def __str__(self):
  52         """
  53         Return a string describing the solution of a linear game.
  54
  55         The three data that are described are,
  56
  57           * The value of the game.
  58           * The optimal strategy of player one.
  59           * The optimal strategy of player two.
  60
  61         The two optimal strategy vectors are indented by two spaces.
  62         """
  63         tpl = 'Game value: {:.7f}\n' \
  64               'Player 1 optimal:{:s}\n' \
  65               'Player 2 optimal:{:s}'
  66
  67         p1_str = '\n{!s}'.format(self.player1_optimal())
  68         p1_str = '\n  '.join(p1_str.splitlines())
  69         p2_str = '\n{!s}'.format(self.player2_optimal())
  70         p2_str = '\n  '.join(p2_str.splitlines())
  71
  72         return tpl.format(self.game_value(), p1_str, p2_str)
  73
  74
  75     def game_value(self):
  76         """
  77         Return the game value for this solution.
  78
  79         Examples
  80         --------
  81
  82             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  83             >>> s.game_value()
  84             10
  85
  86         """
  87         return self._game_value
  88
  89
  90     def player1_optimal(self):
  91         """
  92         Return player one's optimal strategy in this solution.
  93
  94         Examples
  95         --------
  96
  97             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  98             >>> print(s.player1_optimal())
  99             [ 1]
 100             [ 2]
 101             <BLANKLINE>
 102
 103         """
 104         return self._player1_optimal
 105
 106
 107     def player2_optimal(self):
 108         """
 109         Return player two's optimal strategy in this solution.
 110
 111         Examples
 112         --------
 113
 114             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
 115             >>> print(s.player2_optimal())
 116             [ 3]
 117             [ 4]
 118             <BLANKLINE>
 119
 120         """
 121         return self._player2_optimal
 122
 123
 124 class SymmetricLinearGame:
 125     r"""
 126     A representation of a symmetric linear game.
 127
 128     The data for a symmetric linear game are,
 129
 130       * A "payoff" operator ``L``.
 131       * A symmetric cone ``K``.
 132       * Two points ``e1`` and ``e2`` in the interior of ``K``.
 133
 134     The ambient space is assumed to be the span of ``K``.
 135
 136     With those data understood, the game is played as follows. Players
 137     one and two choose points :math:`x` and :math:`y` respectively, from
 138     their respective strategy sets,
 139
 140     .. math::
 141         \begin{aligned}
 142             \Delta_{1}
 143             &=
 144             \left\{
 145               x \in K \ \middle|\ \left\langle x, e_{2} \right\rangle = 1
 146             \right\}\\
 147             \Delta_{2}
 148             &=
 149             \left\{
 150               y \in K \ \middle|\ \left\langle y, e_{1} \right\rangle = 1
 151             \right\}.
 152         \end{aligned}
 153
 154     Afterwards, a "payout" is computed as :math:`\left\langle
 155     L\left(x\right), y \right\rangle` and is paid to player one out of
 156     player two's pocket. The game is therefore zero sum, and we suppose
 157     that player one would like to guarantee himself the largest minimum
 158     payout possible. That is, player one wishes to,
 159
 160     .. math::
 161         \begin{aligned}
 162             \text{maximize }
 163               &\underset{y \in \Delta_{2}}{\min}\left(
 164                 \left\langle L\left(x\right), y \right\rangle
 165               \right)\\
 166             \text{subject to } & x \in \Delta_{1}.
 167         \end{aligned}
 168
 169     Player two has the simultaneous goal to,
 170
 171     .. math::
 172         \begin{aligned}
 173             \text{minimize }
 174               &\underset{x \in \Delta_{1}}{\max}\left(
 175                 \left\langle L\left(x\right), y \right\rangle
 176               \right)\\
 177             \text{subject to } & y \in \Delta_{2}.
 178         \end{aligned}
 179
 180     These goals obviously conflict (the game is zero sum), but an
 181     existence theorem guarantees at least one optimal min-max solution
 182     from which neither player would like to deviate. This class is
 183     able to find such a solution.
 184
 185     Parameters
 186     ----------
 187
 188     L : list of list of float
 189         A matrix represented as a list of ROWS. This representation
 190         agrees with (for example) SageMath and NumPy, but not with CVXOPT
 191         (whose matrix constructor accepts a list of columns).
 192
 193     K : :class:`SymmetricCone`
 194         The symmetric cone instance over which the game is played.
 195
 196     e1 : iterable float
 197         The interior point of ``K`` belonging to player one; it
 198         can be of any iterable type having the correct length.
 199
 200     e2 : iterable float
 201         The interior point of ``K`` belonging to player two; it
 202         can be of any enumerable type having the correct length.
 203
 204     Raises
 205     ------
 206
 207     ValueError
 208         If either ``e1`` or ``e2`` lie outside of the cone ``K``.
 209
 210     Examples
 211     --------
 212
 213         >>> from cones import NonnegativeOrthant
 214         >>> K = NonnegativeOrthant(3)
 215         >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 216         >>> e1 = [1,1,1]
 217         >>> e2 = [1,2,3]
 218         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 219         >>> print(SLG)
 220         The linear game (L, K, e1, e2) where
 221           L = [  1  -5 -15]
 222               [ -1   2  -3]
 223               [-12 -15   1],
 224           K = Nonnegative orthant in the real 3-space,
 225           e1 = [ 1]
 226                [ 1]
 227                [ 1],
 228           e2 = [ 1]
 229                [ 2]
 230                [ 3].
 231
 232     Lists can (and probably should) be used for every argument::
 233
 234         >>> from cones import NonnegativeOrthant
 235         >>> K = NonnegativeOrthant(2)
 236         >>> L = [[1,0],[0,1]]
 237         >>> e1 = [1,1]
 238         >>> e2 = [1,1]
 239         >>> G = SymmetricLinearGame(L, K, e1, e2)
 240         >>> print(G)
 241         The linear game (L, K, e1, e2) where
 242           L = [ 1  0]
 243               [ 0  1],
 244           K = Nonnegative orthant in the real 2-space,
 245           e1 = [ 1]
 246                [ 1],
 247           e2 = [ 1]
 248                [ 1].
 249
 250     The points ``e1`` and ``e2`` can also be passed as some other
 251     enumerable type (of the correct length) without much harm, since
 252     there is no row/column ambiguity::
 253
 254         >>> import cvxopt
 255         >>> import numpy
 256         >>> from cones import NonnegativeOrthant
 257         >>> K = NonnegativeOrthant(2)
 258         >>> L = [[1,0],[0,1]]
 259         >>> e1 = cvxopt.matrix([1,1])
 260         >>> e2 = numpy.matrix([1,1])
 261         >>> G = SymmetricLinearGame(L, K, e1, e2)
 262         >>> print(G)
 263         The linear game (L, K, e1, e2) where
 264           L = [ 1  0]
 265               [ 0  1],
 266           K = Nonnegative orthant in the real 2-space,
 267           e1 = [ 1]
 268                [ 1],
 269           e2 = [ 1]
 270                [ 1].
 271
 272     However, ``L`` will always be intepreted as a list of rows, even
 273     if it is passed as a :class:`cvxopt.base.matrix` which is
 274     otherwise indexed by columns::
 275
 276         >>> import cvxopt
 277         >>> from cones import NonnegativeOrthant
 278         >>> K = NonnegativeOrthant(2)
 279         >>> L = [[1,2],[3,4]]
 280         >>> e1 = [1,1]
 281         >>> e2 = e1
 282         >>> G = SymmetricLinearGame(L, K, e1, e2)
 283         >>> print(G)
 284         The linear game (L, K, e1, e2) where
 285           L = [ 1  2]
 286               [ 3  4],
 287           K = Nonnegative orthant in the real 2-space,
 288           e1 = [ 1]
 289                [ 1],
 290           e2 = [ 1]
 291                [ 1].
 292         >>> L = cvxopt.matrix(L)
 293         >>> print(L)
 294         [ 1  3]
 295         [ 2  4]
 296         <BLANKLINE>
 297         >>> G = SymmetricLinearGame(L, K, e1, e2)
 298         >>> print(G)
 299         The linear game (L, K, e1, e2) where
 300           L = [ 1  2]
 301               [ 3  4],
 302           K = Nonnegative orthant in the real 2-space,
 303           e1 = [ 1]
 304                [ 1],
 305           e2 = [ 1]
 306                [ 1].
 307
 308     """
 309     def __init__(self, L, K, e1, e2):
 310         """
 311         Create a new SymmetricLinearGame object.
 312         """
 313         self._K = K
 314         self._e1 = matrix(e1, (K.dimension(), 1))
 315         self._e2 = matrix(e2, (K.dimension(), 1))
 316
 317         # Our input ``L`` is indexed by rows but CVXOPT matrices are
 318         # indexed by columns, so we need to transpose the input before
 319         # feeding it to CVXOPT.
 320         self._L = matrix(L, (K.dimension(), K.dimension())).trans()
 321
 322         if not K.contains_strict(self._e1):
 323             raise ValueError('the point e1 must lie in the interior of K')
 324
 325         if not K.contains_strict(self._e2):
 326             raise ValueError('the point e2 must lie in the interior of K')
 327
 328     def __str__(self):
 329         """
 330         Return a string representation of this game.
 331         """
 332         tpl = 'The linear game (L, K, e1, e2) where\n' \
 333               '  L = {:s},\n' \
 334               '  K = {!s},\n' \
 335               '  e1 = {:s},\n' \
 336               '  e2 = {:s}.'
 337         indented_L = '\n      '.join(str(self._L).splitlines())
 338         indented_e1 = '\n       '.join(str(self._e1).splitlines())
 339         indented_e2 = '\n       '.join(str(self._e2).splitlines())
 340         return tpl.format(indented_L, str(self._K), indented_e1, indented_e2)
 341
 342
 343     def solution(self):
 344         """
 345         Solve this linear game and return a :class:`Solution`.
 346
 347         Returns
 348         -------
 349
 350         :class:`Solution`
 351             A :class:`Solution` object describing the game's value and
 352             the optimal strategies of both players.
 353
 354         Raises
 355         ------
 356         GameUnsolvableException
 357             If the game could not be solved (if an optimal solution to its
 358             associated cone program was not found).
 359
 360         Examples
 361         --------
 362
 363         This example is computed in Gowda and Ravindran in the section
 364         "The value of a Z-transformation"::
 365
 366             >>> from cones import NonnegativeOrthant
 367             >>> K = NonnegativeOrthant(3)
 368             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 369             >>> e1 = [1,1,1]
 370             >>> e2 = [1,1,1]
 371             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 372             >>> print(SLG.solution())
 373             Game value: -6.1724138
 374             Player 1 optimal:
 375               [ 0.5517241]
 376               [-0.0000000]
 377               [ 0.4482759]
 378             Player 2 optimal:
 379               [0.4482759]
 380               [0.0000000]
 381               [0.5517241]
 382
 383         The value of the following game can be computed using the fact
 384         that the identity is invertible::
 385
 386             >>> from cones import NonnegativeOrthant
 387             >>> K = NonnegativeOrthant(3)
 388             >>> L = [[1,0,0],[0,1,0],[0,0,1]]
 389             >>> e1 = [1,2,3]
 390             >>> e2 = [4,5,6]
 391             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 392             >>> print(SLG.solution())
 393             Game value: 0.0312500
 394             Player 1 optimal:
 395               [0.0312500]
 396               [0.0625000]
 397               [0.0937500]
 398             Player 2 optimal:
 399               [0.1250000]
 400               [0.1562500]
 401               [0.1875000]
 402
 403         """
 404         # The cone "C" that appears in the statement of the CVXOPT
 405         # conelp program.
 406         C = CartesianProduct(self._K, self._K)
 407
 408         # The column vector "b" that appears on the right-hand side of
 409         # Ax = b in the statement of the CVXOPT conelp program.
 410         b = matrix([1], tc='d')
 411
 412         # A column of zeros that fits K.
 413         zero = matrix(0, (self._K.dimension(), 1), tc='d')
 414
 415         # The column vector "h" that appears on the right-hand side of
 416         # Gx + s = h in the statement of the CVXOPT conelp program.
 417         h = matrix([zero, zero])
 418
 419         # The column vector "c" that appears in the objective function
 420         # value <c,x> in the statement of the CVXOPT conelp program.
 421         c = matrix([-1, zero])
 422
 423         # The matrix "G" that appears on the left-hand side of Gx + s = h
 424         # in the statement of the CVXOPT conelp program.
 425         G = append_row(append_col(zero, -identity(self._K.dimension())),
 426                        append_col(self._e1, -self._L))
 427
 428         # The matrix "A" that appears on the right-hand side of Ax = b
 429         # in the statement of the CVXOPT conelp program.
 430         A = matrix([0, self._e2], (1, self._K.dimension() + 1), 'd')
 431
 432         # Actually solve the thing and obtain a dictionary describing
 433         # what happened.
 434         soln_dict = solvers.conelp(c, G, h, C.cvxopt_dims(), A, b)
 435
 436         # The "status" field contains "optimal" if everything went
 437         # according to plan. Other possible values are "primal
 438         # infeasible", "dual infeasible", "unknown", all of which
 439         # mean we didn't get a solution. That should never happen,
 440         # because by construction our game has a solution, and thus
 441         # the cone program should too.
 442         if soln_dict['status'] != 'optimal':
 443             raise GameUnsolvableException(soln_dict)
 444
 445         p1_value = soln_dict['x'][0]
 446         p1_optimal = soln_dict['x'][1:]
 447         p2_optimal = soln_dict['z'][self._K.dimension():]
 448
 449         return Solution(p1_value, p1_optimal, p2_optimal)
 450
 451     def dual(self):
 452         r"""
 453         Return the dual game to this game.
 454
 455         If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
 456         then its dual is :math:`G^{*} =
 457         \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
 458         is symmetric, :math:`K^{*} = K`.
 459
 460         Examples
 461         --------
 462
 463             >>> from cones import NonnegativeOrthant
 464             >>> K = NonnegativeOrthant(3)
 465             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 466             >>> e1 = [1,1,1]
 467             >>> e2 = [1,2,3]
 468             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 469             >>> print(SLG.dual())
 470             The linear game (L, K, e1, e2) where
 471               L = [  1  -1 -12]
 472                   [ -5   2 -15]
 473                   [-15  -3   1],
 474               K = Nonnegative orthant in the real 3-space,
 475               e1 = [ 1]
 476                    [ 2]
 477                    [ 3],
 478               e2 = [ 1]
 479                    [ 1]
 480                    [ 1].
 481
 482         """
 483         # We pass ``self._L`` right back into the constructor, because
 484         # it will be transposed there. And keep in mind that ``self._K``
 485         # is its own dual.
 486         return SymmetricLinearGame(self._L,
 487                                    self._K,
 488                                    self._e2,
 489                                    self._e1)
 490
 491
 492 class SymmetricLinearGameTest(TestCase):
 493     """
 494     Tests for the SymmetricLinearGame and Solution classes.
 495     """
 496
 497     def assert_within_tol(self, first, second):
 498         """
 499         Test that ``first`` and ``second`` are equal within our default
 500         tolerance.
 501         """
 502         self.assertTrue(abs(first - second) < options.ABS_TOL)
 503
 504
 505     def assert_solution_exists(self, L, K, e1, e2):
 506         """
 507         Given the parameters needed to construct a SymmetricLinearGame,
 508         ensure that that game has a solution.
 509         """
 510         G = SymmetricLinearGame(L, K, e1, e2)
 511         soln = G.solution()
 512         L_matrix = matrix(L).trans()
 513         expected = inner_product(L_matrix*soln.player1_optimal(),
 514                                  soln.player2_optimal())
 515         self.assert_within_tol(soln.game_value(), expected)
 516
 517     def test_solution_exists_nonnegative_orthant(self):
 518         """
 519         Every linear game has a solution, so we should be able to solve
 520         every symmetric linear game over the NonnegativeOrthant. Pick
 521         some parameters randomly and give it a shot. The resulting
 522         optimal solutions should give us the optimal game value when we
 523         apply the payoff operator to them.
 524         """
 525         ambient_dim = randint(1, 10)
 526         K = NonnegativeOrthant(ambient_dim)
 527         e1 = [uniform(0.1, 10) for idx in range(K.dimension())]
 528         e2 = [uniform(0.1, 10) for idx in range(K.dimension())]
 529         L = [[uniform(-10, 10) for i in range(K.dimension())]
 530              for j in range(K.dimension())]
 531         self.assert_solution_exists(L, K, e1, e2)
 532
 533     def test_solution_exists_ice_cream(self):
 534         """
 535         Like :meth:`test_solution_exists_nonnegative_orthant`, except
 536         over the ice cream cone.
 537         """
 538         # Use a minimum dimension of two to avoid divide-by-zero in
 539         # the fudge factor we make up later.
 540         ambient_dim = randint(2, 10)
 541         K = IceCream(ambient_dim)
 542         e1 = [1] # Set the "height" of e1 to one
 543         e2 = [1] # And the same for e2
 544
 545         # If we choose the rest of the components of e1,e2 randomly
 546         # between 0 and 1, then the largest the squared norm of the
 547         # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We
 548         # need to make it less than one (the height of the cone) so
 549         # that the whole thing is in the cone. The norm of the
 550         # non-height part is sqrt(dim(K) - 1), and we can divide by
 551         # twice that.
 552         fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0))
 553         e1 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
 554         e2 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
 555         L = [[uniform(-10, 10) for i in range(K.dimension())]
 556              for j in range(K.dimension())]
 557         self.assert_solution_exists(L, K, e1, e2)
 558