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 self._e1 in K:
 323             raise ValueError('the point e1 must lie in the interior of K')
 324
 325         if not self._e2 in K:
 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         p1_value = -soln_dict['primal objective']
 437         p2_value = -soln_dict['dual objective']
 438         p1_optimal = soln_dict['x'][1:]
 439         p2_optimal = soln_dict['z'][self._K.dimension():]
 440
 441         # The "status" field contains "optimal" if everything went
 442         # according to plan. Other possible values are "primal
 443         # infeasible", "dual infeasible", "unknown", all of which mean
 444         # we didn't get a solution. The "infeasible" ones are the
 445         # worst, since they indicate that CVXOPT is convinced the
 446         # problem is infeasible (and that cannot happen).
 447         if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
 448             raise GameUnsolvableException(soln_dict)
 449         elif soln_dict['status'] == 'unknown':
 450             # When we get a status of "unknown", we may still be able
 451             # to salvage a solution out of the returned
 452             # dictionary. Often this is the result of numerical
 453             # difficulty and we can simply check that the primal/dual
 454             # objectives match (within a tolerance) and that the
 455             # primal/dual optimal solutions are within the cone (to a
 456             # tolerance as well).
 457             if (abs(p1_value - p2_value) > options.ABS_TOL):
 458                 raise GameUnsolvableException(soln_dict)
 459             if (p1_optimal not in self._K) or (p2_optimal not in self._K):
 460                 raise GameUnsolvableException(soln_dict)
 461
 462         return Solution(p1_value, p1_optimal, p2_optimal)
 463
 464
 465     def dual(self):
 466         r"""
 467         Return the dual game to this game.
 468
 469         If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
 470         then its dual is :math:`G^{*} =
 471         \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
 472         is symmetric, :math:`K^{*} = K`.
 473
 474         Examples
 475         --------
 476
 477             >>> from cones import NonnegativeOrthant
 478             >>> K = NonnegativeOrthant(3)
 479             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 480             >>> e1 = [1,1,1]
 481             >>> e2 = [1,2,3]
 482             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 483             >>> print(SLG.dual())
 484             The linear game (L, K, e1, e2) where
 485               L = [  1  -1 -12]
 486                   [ -5   2 -15]
 487                   [-15  -3   1],
 488               K = Nonnegative orthant in the real 3-space,
 489               e1 = [ 1]
 490                    [ 2]
 491                    [ 3],
 492               e2 = [ 1]
 493                    [ 1]
 494                    [ 1].
 495
 496         """
 497         # We pass ``self._L`` right back into the constructor, because
 498         # it will be transposed there. And keep in mind that ``self._K``
 499         # is its own dual.
 500         return SymmetricLinearGame(self._L,
 501                                    self._K,
 502                                    self._e2,
 503                                    self._e1)
 504
 505
 506 class SymmetricLinearGameTest(TestCase):
 507     """
 508     Tests for the SymmetricLinearGame and Solution classes.
 509     """
 510
 511     def random_orthant_params(self):
 512         """
 513         Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
 514         random game over the nonnegative orthant.
 515         """
 516         ambient_dim = randint(1, 10)
 517         K = NonnegativeOrthant(ambient_dim)
 518         e1 = [uniform(0.1, 10) for idx in range(K.dimension())]
 519         e2 = [uniform(0.1, 10) for idx in range(K.dimension())]
 520         L = [[uniform(-10, 10) for i in range(K.dimension())]
 521              for j in range(K.dimension())]
 522         return (L, K, e1, e2)
 523
 524
 525     def random_icecream_params(self):
 526         """
 527         Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
 528         random game over the ice cream cone.
 529         """
 530         # Use a minimum dimension of two to avoid divide-by-zero in
 531         # the fudge factor we make up later.
 532         ambient_dim = randint(2, 10)
 533         K = IceCream(ambient_dim)
 534         e1 = [1] # Set the "height" of e1 to one
 535         e2 = [1] # And the same for e2
 536
 537         # If we choose the rest of the components of e1,e2 randomly
 538         # between 0 and 1, then the largest the squared norm of the
 539         # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We
 540         # need to make it less than one (the height of the cone) so
 541         # that the whole thing is in the cone. The norm of the
 542         # non-height part is sqrt(dim(K) - 1), and we can divide by
 543         # twice that.
 544         fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0))
 545         e1 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
 546         e2 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
 547         L = [[uniform(-10, 10) for i in range(K.dimension())]
 548              for j in range(K.dimension())]
 549
 550         return (L, K, e1, e2)
 551
 552
 553     def assert_within_tol(self, first, second):
 554         """
 555         Test that ``first`` and ``second`` are equal within our default
 556         tolerance.
 557         """
 558         self.assertTrue(abs(first - second) < options.ABS_TOL)
 559
 560
 561     def assert_solution_exists(self, L, K, e1, e2):
 562         """
 563         Given the parameters needed to construct a SymmetricLinearGame,
 564         ensure that that game has a solution.
 565         """
 566         G = SymmetricLinearGame(L, K, e1, e2)
 567         soln = G.solution()
 568         L_matrix = matrix(L).trans()
 569         expected = inner_product(L_matrix*soln.player1_optimal(),
 570                                  soln.player2_optimal())
 571         self.assert_within_tol(soln.game_value(), expected)
 572
 573     def test_solution_exists_nonnegative_orthant(self):
 574         """
 575         Every linear game has a solution, so we should be able to solve
 576         every symmetric linear game over the NonnegativeOrthant. Pick
 577         some parameters randomly and give it a shot. The resulting
 578         optimal solutions should give us the optimal game value when we
 579         apply the payoff operator to them.
 580         """
 581         (L, K, e1, e2) = self.random_orthant_params()
 582         self.assert_solution_exists(L, K, e1, e2)
 583
 584     def test_solution_exists_ice_cream(self):
 585         """
 586         Like :meth:`test_solution_exists_nonnegative_orthant`, except
 587         over the ice cream cone.
 588         """
 589         (L, K, e1, e2) = self.random_icecream_params()
 590         self.assert_solution_exists(L, K, e1, e2)
 591
 592
 593     def test_negative_value_Z_operator(self):
 594         """
 595         Test the example given in Gowda/Ravindran of a Z-matrix with
 596         negative game value on the nonnegative orthant.
 597         """
 598         K = NonnegativeOrthant(2)
 599         e1 = [1,1]
 600         e2 = e1
 601         L = [[1,-2],[-2,1]]
 602         G = SymmetricLinearGame(L, K, e1, e2)
 603         self.assertTrue(G.solution().game_value() < -options.ABS_TOL)
 604
 605
 606     def test_nonnegative_scaling_orthant(self):
 607         """
 608         Test that scaling ``L`` by a nonnegative number scales the value
 609         of the game by the same number. Use the nonnegative orthant as
 610         our cone.
 611         """
 612         (L, K, e1, e2) = self.random_orthant_params()
 613         L = matrix(L) # So that we can scale it by alpha below.
 614         G1 = SymmetricLinearGame(L, K, e1, e2)
 615         value1 = G1.solution().game_value()
 616         alpha = uniform(0.1, 10)
 617
 618         G2 = SymmetricLinearGame(alpha*L, K, e1, e2)
 619         value2 = G2.solution().game_value()
 620         self.assert_within_tol(alpha*value1, value2)
 621
 622
 623     def test_nonnegative_scaling_icecream(self):
 624         """
 625         The same test as :meth:`test_nonnegative_scaling_orthant`,
 626         except over the ice cream cone.
 627         """
 628         (L, K, e1, e2) = self.random_icecream_params()
 629         L = matrix(L) # So that we can scale it by alpha below.
 630
 631         G1 = SymmetricLinearGame(L, K, e1, e2)
 632         value1 = G1.solution().game_value()
 633         alpha = uniform(0.1, 10)
 634
 635         G2 = SymmetricLinearGame(alpha*L, K, e1, e2)
 636         value2 = G2.solution().game_value()
 637         self.assert_within_tol(alpha*value1, value2)
 638