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     """
 126     A representation of a symmetric linear game.
 127
 128     The data for a linear game are,
 129
 130       * A "payoff" operator ``L``.
 131       * A cone ``K``.
 132       * A point ``e`` in the interior of ``K``.
 133       * A point ``f`` in the interior of the dual of ``K``.
 134
 135     In a symmetric game, the cone ``K`` is be self-dual. We therefore
 136     name the two interior points ``e1`` and ``e2`` to indicate that
 137     they come from the same cone but are "chosen" by players one and
 138     two respectively.
 139
 140     The ambient space is assumed to be the span of ``K``.
 141
 142     Examples
 143     --------
 144
 145         >>> from cones import NonnegativeOrthant
 146         >>> K = NonnegativeOrthant(3)
 147         >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 148         >>> e1 = [1,1,1]
 149         >>> e2 = [1,2,3]
 150         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 151         >>> print(SLG)
 152         The linear game (L, K, e1, e2) where
 153           L = [  1  -5 -15]
 154               [ -1   2  -3]
 155               [-12 -15   1],
 156           K = Nonnegative orthant in the real 3-space,
 157           e1 = [ 1]
 158                [ 1]
 159                [ 1],
 160           e2 = [ 1]
 161                [ 2]
 162                [ 3].
 163
 164
 165     Lists can (and probably should) be used for every argument::
 166
 167         >>> from cones import NonnegativeOrthant
 168         >>> K = NonnegativeOrthant(2)
 169         >>> L = [[1,0],[0,1]]
 170         >>> e1 = [1,1]
 171         >>> e2 = [1,1]
 172         >>> G = SymmetricLinearGame(L, K, e1, e2)
 173         >>> print(G)
 174         The linear game (L, K, e1, e2) where
 175           L = [ 1  0]
 176               [ 0  1],
 177           K = Nonnegative orthant in the real 2-space,
 178           e1 = [ 1]
 179                [ 1],
 180           e2 = [ 1]
 181                [ 1].
 182
 183     The points ``e1`` and ``e2`` can also be passed as some other
 184     enumerable type (of the correct length) without much harm, since
 185     there is no row/column ambiguity::
 186
 187         >>> import cvxopt
 188         >>> import numpy
 189         >>> from cones import NonnegativeOrthant
 190         >>> K = NonnegativeOrthant(2)
 191         >>> L = [[1,0],[0,1]]
 192         >>> e1 = cvxopt.matrix([1,1])
 193         >>> e2 = numpy.matrix([1,1])
 194         >>> G = SymmetricLinearGame(L, K, e1, e2)
 195         >>> print(G)
 196         The linear game (L, K, e1, e2) where
 197           L = [ 1  0]
 198               [ 0  1],
 199           K = Nonnegative orthant in the real 2-space,
 200           e1 = [ 1]
 201                [ 1],
 202           e2 = [ 1]
 203                [ 1].
 204
 205     However, ``L`` will always be intepreted as a list of rows, even
 206     if it is passed as a :class:`cvxopt.base.matrix` which is
 207     otherwise indexed by columns::
 208
 209         >>> import cvxopt
 210         >>> from cones import NonnegativeOrthant
 211         >>> K = NonnegativeOrthant(2)
 212         >>> L = [[1,2],[3,4]]
 213         >>> e1 = [1,1]
 214         >>> e2 = e1
 215         >>> G = SymmetricLinearGame(L, K, e1, e2)
 216         >>> print(G)
 217         The linear game (L, K, e1, e2) where
 218           L = [ 1  2]
 219               [ 3  4],
 220           K = Nonnegative orthant in the real 2-space,
 221           e1 = [ 1]
 222                [ 1],
 223           e2 = [ 1]
 224                [ 1].
 225         >>> L = cvxopt.matrix(L)
 226         >>> print(L)
 227         [ 1  3]
 228         [ 2  4]
 229         <BLANKLINE>
 230         >>> G = SymmetricLinearGame(L, K, e1, e2)
 231         >>> print(G)
 232         The linear game (L, K, e1, e2) where
 233           L = [ 1  2]
 234               [ 3  4],
 235           K = Nonnegative orthant in the real 2-space,
 236           e1 = [ 1]
 237                [ 1],
 238           e2 = [ 1]
 239                [ 1].
 240
 241     """
 242     def __init__(self, L, K, e1, e2):
 243         """
 244         Create a new SymmetricLinearGame object.
 245
 246         INPUT:
 247
 248           - ``L`` -- an square matrix represented as a list of lists
 249              of real numbers. ``L`` itself is interpreted as a list of
 250              ROWS, which agrees with (for example) SageMath and NumPy,
 251              but not with CVXOPT (whose matrix constructor accepts a
 252              list of columns).
 253
 254           - ``K`` -- a SymmetricCone instance.
 255
 256           - ``e1`` -- the interior point of ``K`` belonging to player one;
 257             it can be of any enumerable type having the correct length.
 258
 259           - ``e2`` -- the interior point of ``K`` belonging to player two;
 260             it can be of any enumerable type having the correct length.
 261
 262         """
 263         self._K = K
 264         self._e1 = matrix(e1, (K.dimension(), 1))
 265         self._e2 = matrix(e2, (K.dimension(), 1))
 266
 267         # Our input ``L`` is indexed by rows but CVXOPT matrices are
 268         # indexed by columns, so we need to transpose the input before
 269         # feeding it to CVXOPT.
 270         self._L = matrix(L, (K.dimension(), K.dimension())).trans()
 271
 272         if not K.contains_strict(self._e1):
 273             raise ValueError('the point e1 must lie in the interior of K')
 274
 275         if not K.contains_strict(self._e2):
 276             raise ValueError('the point e2 must lie in the interior of K')
 277
 278     def __str__(self):
 279         """
 280         Return a string representation of this game.
 281         """
 282         tpl = 'The linear game (L, K, e1, e2) where\n' \
 283               '  L = {:s},\n' \
 284               '  K = {!s},\n' \
 285               '  e1 = {:s},\n' \
 286               '  e2 = {:s}.'
 287         indented_L = '\n      '.join(str(self._L).splitlines())
 288         indented_e1 = '\n       '.join(str(self._e1).splitlines())
 289         indented_e2 = '\n       '.join(str(self._e2).splitlines())
 290         return tpl.format(indented_L, str(self._K), indented_e1, indented_e2)
 291
 292
 293     def solution(self):
 294         """
 295         Solve this linear game and return a Solution object.
 296
 297         OUTPUT:
 298
 299         If the cone program associated with this game could be
 300         successfully solved, then a Solution object containing the
 301         game's value and optimal strategies is returned. If the game
 302         could *not* be solved -- which should never happen -- then a
 303         GameUnsolvableException is raised. It can be printed to get the
 304         raw output from CVXOPT.
 305
 306         Examples
 307         --------
 308
 309         This example is computed in Gowda and Ravindran in the section
 310         "The value of a Z-transformation"::
 311
 312             >>> from cones import NonnegativeOrthant
 313             >>> K = NonnegativeOrthant(3)
 314             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 315             >>> e1 = [1,1,1]
 316             >>> e2 = [1,1,1]
 317             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 318             >>> print(SLG.solution())
 319             Game value: -6.1724138
 320             Player 1 optimal:
 321               [ 0.5517241]
 322               [-0.0000000]
 323               [ 0.4482759]
 324             Player 2 optimal:
 325               [0.4482759]
 326               [0.0000000]
 327               [0.5517241]
 328
 329         The value of the following game can be computed using the fact
 330         that the identity is invertible::
 331
 332             >>> from cones import NonnegativeOrthant
 333             >>> K = NonnegativeOrthant(3)
 334             >>> L = [[1,0,0],[0,1,0],[0,0,1]]
 335             >>> e1 = [1,2,3]
 336             >>> e2 = [4,5,6]
 337             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 338             >>> print(SLG.solution())
 339             Game value: 0.0312500
 340             Player 1 optimal:
 341               [0.0312500]
 342               [0.0625000]
 343               [0.0937500]
 344             Player 2 optimal:
 345               [0.1250000]
 346               [0.1562500]
 347               [0.1875000]
 348
 349         """
 350         # The cone "C" that appears in the statement of the CVXOPT
 351         # conelp program.
 352         C = CartesianProduct(self._K, self._K)
 353
 354         # The column vector "b" that appears on the right-hand side of
 355         # Ax = b in the statement of the CVXOPT conelp program.
 356         b = matrix([1], tc='d')
 357
 358         # A column of zeros that fits K.
 359         zero = matrix(0, (self._K.dimension(), 1), tc='d')
 360
 361         # The column vector "h" that appears on the right-hand side of
 362         # Gx + s = h in the statement of the CVXOPT conelp program.
 363         h = matrix([zero, zero])
 364
 365         # The column vector "c" that appears in the objective function
 366         # value <c,x> in the statement of the CVXOPT conelp program.
 367         c = matrix([-1, zero])
 368
 369         # The matrix "G" that appears on the left-hand side of Gx + s = h
 370         # in the statement of the CVXOPT conelp program.
 371         G = append_row(append_col(zero, -identity(self._K.dimension())),
 372                        append_col(self._e1, -self._L))
 373
 374         # The matrix "A" that appears on the right-hand side of Ax = b
 375         # in the statement of the CVXOPT conelp program.
 376         A = matrix([0, self._e2], (1, self._K.dimension() + 1), 'd')
 377
 378         # Actually solve the thing and obtain a dictionary describing
 379         # what happened.
 380         soln_dict = solvers.conelp(c, G, h, C.cvxopt_dims(), A, b)
 381
 382         # The "status" field contains "optimal" if everything went
 383         # according to plan. Other possible values are "primal
 384         # infeasible", "dual infeasible", "unknown", all of which
 385         # mean we didn't get a solution. That should never happen,
 386         # because by construction our game has a solution, and thus
 387         # the cone program should too.
 388         if soln_dict['status'] != 'optimal':
 389             raise GameUnsolvableException(soln_dict)
 390
 391         p1_value = soln_dict['x'][0]
 392         p1_optimal = soln_dict['x'][1:]
 393         p2_optimal = soln_dict['z'][self._K.dimension():]
 394
 395         return Solution(p1_value, p1_optimal, p2_optimal)
 396
 397     def dual(self):
 398         """
 399         Return the dual game to this game.
 400
 401         Examples
 402         --------
 403
 404             >>> from cones import NonnegativeOrthant
 405             >>> K = NonnegativeOrthant(3)
 406             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 407             >>> e1 = [1,1,1]
 408             >>> e2 = [1,2,3]
 409             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 410             >>> print(SLG.dual())
 411             The linear game (L, K, e1, e2) where
 412               L = [  1  -1 -12]
 413                   [ -5   2 -15]
 414                   [-15  -3   1],
 415               K = Nonnegative orthant in the real 3-space,
 416               e1 = [ 1]
 417                    [ 2]
 418                    [ 3],
 419               e2 = [ 1]
 420                    [ 1]
 421                    [ 1].
 422
 423         """
 424         # We pass ``self._L`` right back into the constructor, because
 425         # it will be transposed there. And keep in mind that ``self._K``
 426         # is its own dual.
 427         return SymmetricLinearGame(self._L,
 428                                    self._K,
 429                                    self._e2,
 430                                    self._e1)
 431
 432
 433 class SymmetricLinearGameTest(TestCase):
 434     """
 435     Tests for the SymmetricLinearGame and Solution classes.
 436     """
 437
 438     def assert_within_tol(self, first, second):
 439         """
 440         Test that ``first`` and ``second`` are equal within our default
 441         tolerance.
 442         """
 443         self.assertTrue(abs(first - second) < options.ABS_TOL)
 444
 445
 446     def assert_solution_exists(self, L, K, e1, e2):
 447         """
 448         Given the parameters needed to construct a SymmetricLinearGame,
 449         ensure that that game has a solution.
 450         """
 451         G = SymmetricLinearGame(L, K, e1, e2)
 452         soln = G.solution()
 453         L_matrix = matrix(L).trans()
 454         expected = inner_product(L_matrix*soln.player1_optimal(),
 455                                  soln.player2_optimal())
 456         self.assert_within_tol(soln.game_value(), expected)
 457
 458     def test_solution_exists_nonnegative_orthant(self):
 459         """
 460         Every linear game has a solution, so we should be able to solve
 461         every symmetric linear game over the NonnegativeOrthant. Pick
 462         some parameters randomly and give it a shot. The resulting
 463         optimal solutions should give us the optimal game value when we
 464         apply the payoff operator to them.
 465         """
 466         ambient_dim = randint(1, 10)
 467         K = NonnegativeOrthant(ambient_dim)
 468         e1 = [uniform(0.1, 10) for idx in range(K.dimension())]
 469         e2 = [uniform(0.1, 10) for idx in range(K.dimension())]
 470         L = [[uniform(-10, 10) for i in range(K.dimension())]
 471              for j in range(K.dimension())]
 472         self.assert_solution_exists(L, K, e1, e2)
 473
 474     def test_solution_exists_ice_cream(self):
 475         """
 476         Like :meth:`test_solution_exists_nonnegative_orthant`, except
 477         over the ice cream cone.
 478         """
 479         # Use a minimum dimension of two to avoid divide-by-zero in
 480         # the fudge factor we make up later.
 481         ambient_dim = randint(2, 10)
 482         K = IceCream(ambient_dim)
 483         e1 = [1] # Set the "height" of e1 to one
 484         e2 = [1] # And the same for e2
 485
 486         # If we choose the rest of the components of e1,e2 randomly
 487         # between 0 and 1, then the largest the squared norm of the
 488         # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We
 489         # need to make it less than one (the height of the cone) so
 490         # that the whole thing is in the cone. The norm of the
 491         # non-height part is sqrt(dim(K) - 1), and we can divide by
 492         # twice that.
 493         fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0))
 494         e1 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
 495         e2 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
 496         L = [[uniform(-10, 10) for i in range(K.dimension())]
 497              for j in range(K.dimension())]
 498         self.assert_solution_exists(L, K, e1, e2)
 499