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 from cvxopt import matrix, printing, solvers
   8 from .cones import CartesianProduct
   9 from .errors import GameUnsolvableException, PoorScalingException
  10 from .matrices import (append_col, append_row, condition_number, identity,
  11                        inner_product, norm, specnorm)
  12 from .options import ABS_TOL, FLOAT_FORMAT, DEBUG_FLOAT_FORMAT
  13
  14 printing.options['dformat'] = FLOAT_FORMAT
  15
  16
  17 class Solution:
  18     """
  19     A representation of the solution of a linear game. It should contain
  20     the value of the game, and both players' strategies.
  21
  22     Examples
  23     --------
  24
  25         >>> print(Solution(10, matrix([1,2]), matrix([3,4])))
  26         Game value: 10.000...
  27         Player 1 optimal:
  28           [ 1]
  29           [ 2]
  30         Player 2 optimal:
  31           [ 3]
  32           [ 4]
  33
  34     """
  35     def __init__(self, game_value, p1_optimal, p2_optimal):
  36         """
  37         Create a new Solution object from a game value and two optimal
  38         strategies for the players.
  39         """
  40         self._game_value = game_value
  41         self._player1_optimal = p1_optimal
  42         self._player2_optimal = p2_optimal
  43
  44     def __str__(self):
  45         """
  46         Return a string describing the solution of a linear game.
  47
  48         The three data that are described are,
  49
  50           * The value of the game.
  51           * The optimal strategy of player one.
  52           * The optimal strategy of player two.
  53
  54         The two optimal strategy vectors are indented by two spaces.
  55         """
  56         tpl = 'Game value: {:.7f}\n' \
  57               'Player 1 optimal:{:s}\n' \
  58               'Player 2 optimal:{:s}'
  59
  60         p1_str = '\n{!s}'.format(self.player1_optimal())
  61         p1_str = '\n  '.join(p1_str.splitlines())
  62         p2_str = '\n{!s}'.format(self.player2_optimal())
  63         p2_str = '\n  '.join(p2_str.splitlines())
  64
  65         return tpl.format(self.game_value(), p1_str, p2_str)
  66
  67
  68     def game_value(self):
  69         """
  70         Return the game value for this solution.
  71
  72         Examples
  73         --------
  74
  75             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  76             >>> s.game_value()
  77             10
  78
  79         """
  80         return self._game_value
  81
  82
  83     def player1_optimal(self):
  84         """
  85         Return player one's optimal strategy in this solution.
  86
  87         Examples
  88         --------
  89
  90             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  91             >>> print(s.player1_optimal())
  92             [ 1]
  93             [ 2]
  94             <BLANKLINE>
  95
  96         """
  97         return self._player1_optimal
  98
  99
 100     def player2_optimal(self):
 101         """
 102         Return player two's optimal strategy in this solution.
 103
 104         Examples
 105         --------
 106
 107             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
 108             >>> print(s.player2_optimal())
 109             [ 3]
 110             [ 4]
 111             <BLANKLINE>
 112
 113         """
 114         return self._player2_optimal
 115
 116
 117 class SymmetricLinearGame:
 118     r"""
 119     A representation of a symmetric linear game.
 120
 121     The data for a symmetric linear game are,
 122
 123       * A "payoff" operator ``L``.
 124       * A symmetric cone ``K``.
 125       * Two points ``e1`` and ``e2`` in the interior of ``K``.
 126
 127     The ambient space is assumed to be the span of ``K``.
 128
 129     With those data understood, the game is played as follows. Players
 130     one and two choose points :math:`x` and :math:`y` respectively, from
 131     their respective strategy sets,
 132
 133     .. math::
 134         \begin{aligned}
 135             \Delta_{1}
 136             &=
 137             \left\{
 138               x \in K \ \middle|\ \left\langle x, e_{2} \right\rangle = 1
 139             \right\}\\
 140             \Delta_{2}
 141             &=
 142             \left\{
 143               y \in K \ \middle|\ \left\langle y, e_{1} \right\rangle = 1
 144             \right\}.
 145         \end{aligned}
 146
 147     Afterwards, a "payout" is computed as :math:`\left\langle
 148     L\left(x\right), y \right\rangle` and is paid to player one out of
 149     player two's pocket. The game is therefore zero sum, and we suppose
 150     that player one would like to guarantee himself the largest minimum
 151     payout possible. That is, player one wishes to,
 152
 153     .. math::
 154         \begin{aligned}
 155             \text{maximize }
 156               &\underset{y \in \Delta_{2}}{\min}\left(
 157                 \left\langle L\left(x\right), y \right\rangle
 158               \right)\\
 159             \text{subject to } & x \in \Delta_{1}.
 160         \end{aligned}
 161
 162     Player two has the simultaneous goal to,
 163
 164     .. math::
 165         \begin{aligned}
 166             \text{minimize }
 167               &\underset{x \in \Delta_{1}}{\max}\left(
 168                 \left\langle L\left(x\right), y \right\rangle
 169               \right)\\
 170             \text{subject to } & y \in \Delta_{2}.
 171         \end{aligned}
 172
 173     These goals obviously conflict (the game is zero sum), but an
 174     existence theorem guarantees at least one optimal min-max solution
 175     from which neither player would like to deviate. This class is
 176     able to find such a solution.
 177
 178     Parameters
 179     ----------
 180
 181     L : list of list of float
 182         A matrix represented as a list of **rows**. This representation
 183         agrees with (for example) `SageMath <http://www.sagemath.org/>`_
 184         and `NumPy <http://www.numpy.org/>`_, but not with CVXOPT (whose
 185         matrix constructor accepts a list of columns). In reality, ``L``
 186         can be any iterable type of the correct length; however, you
 187         should be extremely wary of the way we interpret anything other
 188         than a list of rows.
 189
 190     K : dunshire.cones.SymmetricCone
 191         The symmetric cone instance over which the game is played.
 192
 193     e1 : iterable float
 194         The interior point of ``K`` belonging to player one; it
 195         can be of any iterable type having the correct length.
 196
 197     e2 : iterable float
 198         The interior point of ``K`` belonging to player two; it
 199         can be of any enumerable type having the correct length.
 200
 201     Raises
 202     ------
 203
 204     ValueError
 205         If either ``e1`` or ``e2`` lie outside of the cone ``K``.
 206
 207     Examples
 208     --------
 209
 210         >>> from dunshire import *
 211         >>> K = NonnegativeOrthant(3)
 212         >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 213         >>> e1 = [1,1,1]
 214         >>> e2 = [1,2,3]
 215         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 216         >>> print(SLG)
 217         The linear game (L, K, e1, e2) where
 218           L = [  1  -5 -15]
 219               [ -1   2  -3]
 220               [-12 -15   1],
 221           K = Nonnegative orthant in the real 3-space,
 222           e1 = [ 1]
 223                [ 1]
 224                [ 1],
 225           e2 = [ 1]
 226                [ 2]
 227                [ 3]
 228
 229     Lists can (and probably should) be used for every argument::
 230
 231         >>> from dunshire import *
 232         >>> K = NonnegativeOrthant(2)
 233         >>> L = [[1,0],[0,1]]
 234         >>> e1 = [1,1]
 235         >>> e2 = [1,1]
 236         >>> G = SymmetricLinearGame(L, K, e1, e2)
 237         >>> print(G)
 238         The linear game (L, K, e1, e2) where
 239           L = [ 1  0]
 240               [ 0  1],
 241           K = Nonnegative orthant in the real 2-space,
 242           e1 = [ 1]
 243                [ 1],
 244           e2 = [ 1]
 245                [ 1]
 246
 247     The points ``e1`` and ``e2`` can also be passed as some other
 248     enumerable type (of the correct length) without much harm, since
 249     there is no row/column ambiguity::
 250
 251         >>> import cvxopt
 252         >>> import numpy
 253         >>> from dunshire import *
 254         >>> K = NonnegativeOrthant(2)
 255         >>> L = [[1,0],[0,1]]
 256         >>> e1 = cvxopt.matrix([1,1])
 257         >>> e2 = numpy.matrix([1,1])
 258         >>> G = SymmetricLinearGame(L, K, e1, e2)
 259         >>> print(G)
 260         The linear game (L, K, e1, e2) where
 261           L = [ 1  0]
 262               [ 0  1],
 263           K = Nonnegative orthant in the real 2-space,
 264           e1 = [ 1]
 265                [ 1],
 266           e2 = [ 1]
 267                [ 1]
 268
 269     However, ``L`` will always be intepreted as a list of rows, even
 270     if it is passed as a :class:`cvxopt.base.matrix` which is
 271     otherwise indexed by columns::
 272
 273         >>> import cvxopt
 274         >>> from dunshire import *
 275         >>> K = NonnegativeOrthant(2)
 276         >>> L = [[1,2],[3,4]]
 277         >>> e1 = [1,1]
 278         >>> e2 = e1
 279         >>> G = SymmetricLinearGame(L, K, e1, e2)
 280         >>> print(G)
 281         The linear game (L, K, e1, e2) where
 282           L = [ 1  2]
 283               [ 3  4],
 284           K = Nonnegative orthant in the real 2-space,
 285           e1 = [ 1]
 286                [ 1],
 287           e2 = [ 1]
 288                [ 1]
 289         >>> L = cvxopt.matrix(L)
 290         >>> print(L)
 291         [ 1  3]
 292         [ 2  4]
 293         <BLANKLINE>
 294         >>> G = SymmetricLinearGame(L, K, e1, e2)
 295         >>> print(G)
 296         The linear game (L, K, e1, e2) where
 297           L = [ 1  2]
 298               [ 3  4],
 299           K = Nonnegative orthant in the real 2-space,
 300           e1 = [ 1]
 301                [ 1],
 302           e2 = [ 1]
 303                [ 1]
 304
 305     """
 306     def __init__(self, L, K, e1, e2):
 307         """
 308         Create a new SymmetricLinearGame object.
 309         """
 310         self._K = K
 311         self._e1 = matrix(e1, (K.dimension(), 1))
 312         self._e2 = matrix(e2, (K.dimension(), 1))
 313
 314         # Our input ``L`` is indexed by rows but CVXOPT matrices are
 315         # indexed by columns, so we need to transpose the input before
 316         # feeding it to CVXOPT.
 317         self._L = matrix(L, (K.dimension(), K.dimension())).trans()
 318
 319         if not self._e1 in K:
 320             raise ValueError('the point e1 must lie in the interior of K')
 321
 322         if not self._e2 in K:
 323             raise ValueError('the point e2 must lie in the interior of K')
 324
 325         # Initial value of cached method.
 326         self._L_specnorm_value = None
 327
 328
 329     def __str__(self):
 330         """
 331         Return a string representation of this game.
 332         """
 333         tpl = 'The linear game (L, K, e1, e2) where\n' \
 334               '  L = {:s},\n' \
 335               '  K = {!s},\n' \
 336               '  e1 = {:s},\n' \
 337               '  e2 = {:s}'
 338         indented_L = '\n      '.join(str(self.L()).splitlines())
 339         indented_e1 = '\n       '.join(str(self.e1()).splitlines())
 340         indented_e2 = '\n       '.join(str(self.e2()).splitlines())
 341
 342         return tpl.format(indented_L,
 343                           str(self.K()),
 344                           indented_e1,
 345                           indented_e2)
 346
 347
 348     def L(self):
 349         """
 350         Return the matrix ``L`` passed to the constructor.
 351
 352         Returns
 353         -------
 354
 355         matrix
 356             The matrix that defines this game's :meth:`payoff` operator.
 357
 358         Examples
 359         --------
 360
 361             >>> from dunshire import *
 362             >>> K = NonnegativeOrthant(3)
 363             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 364             >>> e1 = [1,1,1]
 365             >>> e2 = [1,2,3]
 366             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 367             >>> print(SLG.L())
 368             [  1  -5 -15]
 369             [ -1   2  -3]
 370             [-12 -15   1]
 371             <BLANKLINE>
 372
 373         """
 374         return self._L
 375
 376
 377     def K(self):
 378         """
 379         Return the cone over which this game is played.
 380
 381         Returns
 382         -------
 383
 384         SymmetricCone
 385             The :class:`SymmetricCone` over which this game is played.
 386
 387         Examples
 388         --------
 389
 390             >>> from dunshire import *
 391             >>> K = NonnegativeOrthant(3)
 392             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 393             >>> e1 = [1,1,1]
 394             >>> e2 = [1,2,3]
 395             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 396             >>> print(SLG.K())
 397             Nonnegative orthant in the real 3-space
 398
 399         """
 400         return self._K
 401
 402
 403     def e1(self):
 404         """
 405         Return player one's interior point.
 406
 407         Returns
 408         -------
 409
 410         matrix
 411             The point interior to :meth:`K` affiliated with player one.
 412
 413         Examples
 414         --------
 415
 416             >>> from dunshire import *
 417             >>> K = NonnegativeOrthant(3)
 418             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 419             >>> e1 = [1,1,1]
 420             >>> e2 = [1,2,3]
 421             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 422             >>> print(SLG.e1())
 423             [ 1]
 424             [ 1]
 425             [ 1]
 426             <BLANKLINE>
 427
 428         """
 429         return self._e1
 430
 431
 432     def e2(self):
 433         """
 434         Return player two's interior point.
 435
 436         Returns
 437         -------
 438
 439         matrix
 440             The point interior to :meth:`K` affiliated with player one.
 441
 442         Examples
 443         --------
 444
 445             >>> from dunshire import *
 446             >>> K = NonnegativeOrthant(3)
 447             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 448             >>> e1 = [1,1,1]
 449             >>> e2 = [1,2,3]
 450             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 451             >>> print(SLG.e2())
 452             [ 1]
 453             [ 2]
 454             [ 3]
 455             <BLANKLINE>
 456
 457         """
 458         return self._e2
 459
 460
 461     def payoff(self, strategy1, strategy2):
 462         r"""
 463         Return the payoff associated with ``strategy1`` and ``strategy2``.
 464
 465         The payoff operator takes pairs of strategies to a real
 466         number. For example, if player one's strategy is :math:`x` and
 467         player two's strategy is :math:`y`, then the associated payoff
 468         is :math:`\left\langle L\left(x\right),y \right\rangle \in
 469         \mathbb{R}`. Here, :math:`L` denotes the same linear operator as
 470         :meth:`L`. This method computes the payoff given the two
 471         players' strategies.
 472
 473         Parameters
 474         ----------
 475
 476         strategy1 : matrix
 477             Player one's strategy.
 478
 479         strategy2 : matrix
 480             Player two's strategy.
 481
 482         Returns
 483         -------
 484
 485         float
 486             The payoff for the game when player one plays ``strategy1``
 487             and player two plays ``strategy2``.
 488
 489         Examples
 490         --------
 491
 492         The value of the game should be the payoff at the optimal
 493         strategies::
 494
 495             >>> from dunshire import *
 496             >>> K = NonnegativeOrthant(3)
 497             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 498             >>> e1 = [1,1,1]
 499             >>> e2 = [1,1,1]
 500             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 501             >>> soln = SLG.solution()
 502             >>> x_bar = soln.player1_optimal()
 503             >>> y_bar = soln.player2_optimal()
 504             >>> SLG.payoff(x_bar, y_bar) == soln.game_value()
 505             True
 506
 507         """
 508         return inner_product(self.L()*strategy1, strategy2)
 509
 510
 511     def dimension(self):
 512         """
 513         Return the dimension of this game.
 514
 515         The dimension of a game is not needed for the theory, but it is
 516         useful for the implementation. We define the dimension of a game
 517         to be the dimension of its underlying cone. Or what is the same,
 518         the dimension of the space from which the strategies are chosen.
 519
 520         Returns
 521         -------
 522
 523         int
 524             The dimension of the cone :meth:`K`, or of the space where
 525             this game is played.
 526
 527         Examples
 528         --------
 529
 530         The dimension of a game over the nonnegative quadrant in the
 531         plane should be two (the dimension of the plane)::
 532
 533             >>> from dunshire import *
 534             >>> K = NonnegativeOrthant(2)
 535             >>> L = [[1,-5],[-1,2]]
 536             >>> e1 = [1,1]
 537             >>> e2 = [1,4]
 538             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 539             >>> SLG.dimension()
 540             2
 541
 542         """
 543         return self.K().dimension()
 544
 545
 546     def _zero(self):
 547         """
 548         Return a column of zeros that fits ``K``.
 549
 550         This is used in our CVXOPT construction.
 551
 552         .. warning::
 553
 554             It is not safe to cache any of the matrices passed to
 555             CVXOPT, because it can clobber them.
 556
 557         Returns
 558         -------
 559
 560         matrix
 561             A ``self.dimension()``-by-``1`` column vector of zeros.
 562
 563         Examples
 564         --------
 565
 566             >>> from dunshire import *
 567             >>> K = NonnegativeOrthant(3)
 568             >>> L = identity(3)
 569             >>> e1 = [1,1,1]
 570             >>> e2 = e1
 571             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 572             >>> print(SLG._zero())
 573             [0.0000000]
 574             [0.0000000]
 575             [0.0000000]
 576             <BLANKLINE>
 577
 578         """
 579         return matrix(0, (self.dimension(), 1), tc='d')
 580
 581
 582     def A(self):
 583         r"""
 584         Return the matrix ``A`` used in our CVXOPT construction.
 585
 586         This matrix :math:`A` appears on the right-hand side of
 587         :math:`Ax = b` in the `statement of the CVXOPT conelp program
 588         <http://cvxopt.org/userguide/coneprog.html#linear-cone-programs>`_.
 589
 590         .. warning::
 591
 592             It is not safe to cache any of the matrices passed to
 593             CVXOPT, because it can clobber them.
 594
 595         Returns
 596         -------
 597
 598         matrix
 599             A ``1``-by-``(1 + self.dimension())`` row vector. Its first
 600             entry is zero, and the rest are the entries of :meth:`e2`.
 601
 602         Examples
 603         --------
 604
 605             >>> from dunshire import *
 606             >>> K = NonnegativeOrthant(3)
 607             >>> L = [[1,1,1],[1,1,1],[1,1,1]]
 608             >>> e1 = [1,1,1]
 609             >>> e2 = [1,2,3]
 610             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 611             >>> print(SLG.A())
 612             [0.0000000 1.0000000 2.0000000 3.0000000]
 613             <BLANKLINE>
 614
 615         """
 616         return matrix([0, self.e2()], (1, self.dimension() + 1), 'd')
 617
 618
 619
 620     def G(self):
 621         r"""
 622         Return the matrix ``G`` used in our CVXOPT construction.
 623
 624         Thus matrix :math:`G` appears on the left-hand side of :math:`Gx
 625         + s = h` in the `statement of the CVXOPT conelp program
 626         <http://cvxopt.org/userguide/coneprog.html#linear-cone-programs>`_.
 627
 628         .. warning::
 629
 630             It is not safe to cache any of the matrices passed to
 631             CVXOPT, because it can clobber them.
 632
 633         Returns
 634         -------
 635
 636         matrix
 637             A ``2*self.dimension()``-by-``(1 + self.dimension())`` matrix.
 638
 639         Examples
 640         --------
 641
 642             >>> from dunshire import *
 643             >>> K = NonnegativeOrthant(3)
 644             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 645             >>> e1 = [1,2,3]
 646             >>> e2 = [1,1,1]
 647             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 648             >>> print(SLG.G())
 649             [  0.0000000  -1.0000000   0.0000000   0.0000000]
 650             [  0.0000000   0.0000000  -1.0000000   0.0000000]
 651             [  0.0000000   0.0000000   0.0000000  -1.0000000]
 652             [  1.0000000  -4.0000000  -5.0000000  -6.0000000]
 653             [  2.0000000  -7.0000000  -8.0000000  -9.0000000]
 654             [  3.0000000 -10.0000000 -11.0000000 -12.0000000]
 655             <BLANKLINE>
 656
 657         """
 658         identity_matrix = identity(self.dimension())
 659         return append_row(append_col(self._zero(), -identity_matrix),
 660                           append_col(self.e1(), -self.L()))
 661
 662
 663     def c(self):
 664         r"""
 665         Return the vector ``c`` used in our CVXOPT construction.
 666
 667         The column vector :math:`c` appears in the objective function
 668         value :math:`\left\langle c,x \right\rangle` in the `statement
 669         of the CVXOPT conelp program
 670         <http://cvxopt.org/userguide/coneprog.html#linear-cone-programs>`_.
 671
 672         .. warning::
 673
 674             It is not safe to cache any of the matrices passed to
 675             CVXOPT, because it can clobber them.
 676
 677         Returns
 678         -------
 679
 680         matrix
 681             A :meth:`dimension`-by-``1`` column vector.
 682
 683         Examples
 684         --------
 685
 686             >>> from dunshire import *
 687             >>> K = NonnegativeOrthant(3)
 688             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 689             >>> e1 = [1,2,3]
 690             >>> e2 = [1,1,1]
 691             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 692             >>> print(SLG.c())
 693             [-1.0000000]
 694             [ 0.0000000]
 695             [ 0.0000000]
 696             [ 0.0000000]
 697             <BLANKLINE>
 698
 699         """
 700         return matrix([-1, self._zero()])
 701
 702
 703     def C(self):
 704         """
 705         Return the cone ``C`` used in our CVXOPT construction.
 706
 707         This is the cone over which the `CVXOPT conelp program
 708         <http://cvxopt.org/userguide/coneprog.html#linear-cone-programs>`_
 709         takes place.
 710
 711         Returns
 712         -------
 713
 714         CartesianProduct
 715             The cartesian product of ``K`` with itself.
 716
 717         Examples
 718         --------
 719
 720             >>> from dunshire import *
 721             >>> K = NonnegativeOrthant(3)
 722             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 723             >>> e1 = [1,2,3]
 724             >>> e2 = [1,1,1]
 725             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 726             >>> print(SLG.C())
 727             Cartesian product of dimension 6 with 2 factors:
 728               * Nonnegative orthant in the real 3-space
 729               * Nonnegative orthant in the real 3-space
 730
 731         """
 732         return CartesianProduct(self._K, self._K)
 733
 734     def h(self):
 735         r"""
 736         Return the ``h`` vector used in our CVXOPT construction.
 737
 738         The :math:`h` vector appears on the right-hand side of :math:`Gx
 739         + s = h` in the `statement of the CVXOPT conelp program
 740         <http://cvxopt.org/userguide/coneprog.html#linear-cone-programs>`_.
 741
 742         .. warning::
 743
 744             It is not safe to cache any of the matrices passed to
 745             CVXOPT, because it can clobber them.
 746
 747         Returns
 748         -------
 749
 750         matrix
 751             A ``2*self.dimension()``-by-``1`` column vector of zeros.
 752
 753         Examples
 754         --------
 755
 756             >>> from dunshire import *
 757             >>> K = NonnegativeOrthant(3)
 758             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 759             >>> e1 = [1,2,3]
 760             >>> e2 = [1,1,1]
 761             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 762             >>> print(SLG.h())
 763             [0.0000000]
 764             [0.0000000]
 765             [0.0000000]
 766             [0.0000000]
 767             [0.0000000]
 768             [0.0000000]
 769             <BLANKLINE>
 770
 771         """
 772
 773         return matrix([self._zero(), self._zero()])
 774
 775
 776     @staticmethod
 777     def b():
 778         r"""
 779         Return the ``b`` vector used in our CVXOPT construction.
 780
 781         The vector :math:`b` appears on the right-hand side of :math:`Ax
 782         = b` in the `statement of the CVXOPT conelp program
 783         <http://cvxopt.org/userguide/coneprog.html#linear-cone-programs>`_.
 784
 785         This method is static because the dimensions and entries of
 786         ``b`` are known beforehand, and don't depend on any other
 787         properties of the game.
 788
 789         .. warning::
 790
 791             It is not safe to cache any of the matrices passed to
 792             CVXOPT, because it can clobber them.
 793
 794         Returns
 795         -------
 796
 797         matrix
 798             A ``1``-by-``1`` matrix containing a single entry ``1``.
 799
 800         Examples
 801         --------
 802
 803             >>> from dunshire import *
 804             >>> K = NonnegativeOrthant(3)
 805             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 806             >>> e1 = [1,2,3]
 807             >>> e2 = [1,1,1]
 808             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 809             >>> print(SLG.b())
 810             [1.0000000]
 811             <BLANKLINE>
 812
 813         """
 814         return matrix([1], tc='d')
 815
 816
 817     def player1_start(self):
 818         """
 819         Return a feasible starting point for player one.
 820
 821         This starting point is for the CVXOPT formulation and not for
 822         the original game. The basic premise is that if you scale
 823         :meth:`e2` by the reciprocal of its squared norm, then you get a
 824         point in :meth:`K` that makes a unit inner product with
 825         :meth:`e2`. We then get to choose the primal objective function
 826         value such that the constraint involving :meth:`L` is satisfied.
 827
 828         Returns
 829         -------
 830
 831         dict
 832             A dictionary with two keys, ``'x'`` and ``'s'``, which
 833             contain the vectors of the same name in the CVXOPT primal
 834             problem formulation.
 835
 836             The vector ``x`` consists of the primal objective function
 837             value concatenated with the strategy (for player one) that
 838             achieves it. The vector ``s`` is essentially a dummy
 839             variable, and is computed from the equality constraing in
 840             the CVXOPT primal problem.
 841
 842         """
 843         p = self.e2() / (norm(self.e2()) ** 2)
 844         dist = self.K().ball_radius(self.e1())
 845         nu = - self._L_specnorm()/(dist*norm(self.e2()))
 846         x = matrix([nu, p], (self.dimension() + 1, 1))
 847         s = - self.G()*x
 848
 849         return {'x': x, 's': s}
 850
 851
 852     def player2_start(self):
 853         """
 854         Return a feasible starting point for player two.
 855
 856         This starting point is for the CVXOPT formulation and not for
 857         the original game. The basic premise is that if you scale
 858         :meth:`e1` by the reciprocal of its squared norm, then you get a
 859         point in :meth:`K` that makes a unit inner product with
 860         :meth:`e1`. We then get to choose the dual objective function
 861         value such that the constraint involving :meth:`L` is satisfied.
 862
 863         Returns
 864         -------
 865
 866         dict
 867             A dictionary with two keys, ``'y'`` and ``'z'``, which
 868             contain the vectors of the same name in the CVXOPT dual
 869             problem formulation.
 870
 871             The ``1``-by-``1`` vector ``y`` consists of the dual
 872             objective function value. The last :meth:`dimension` entries
 873             of the vector ``z`` contain the strategy (for player two)
 874             that achieves it. The remaining entries of ``z`` are
 875             essentially dummy variables, computed from the equality
 876             constraint in the CVXOPT dual problem.
 877
 878         """
 879         q = self.e1() / (norm(self.e1()) ** 2)
 880         dist = self.K().ball_radius(self.e2())
 881         omega = self._L_specnorm()/(dist*norm(self.e1()))
 882         y = matrix([omega])
 883         z2 = q
 884         z1 = y*self.e2() - self.L().trans()*z2
 885         z = matrix([z1, z2], (self.dimension()*2, 1))
 886
 887         return {'y': y, 'z': z}
 888
 889
 890     def _L_specnorm(self):
 891         """
 892         Compute the spectral norm of :meth:`L` and cache it.
 893
 894         The spectral norm of the matrix :meth:`L` is used in a few
 895         places. Since it can be expensive to compute, we want to cache
 896         its value. That is not possible in :func:`specnorm`, which lies
 897         outside of a class, so this is the place to do it.
 898
 899         Returns
 900         -------
 901
 902         float
 903             A nonnegative real number; the largest singular value of
 904             the matrix :meth:`L`.
 905
 906         Examples
 907         --------
 908
 909             >>> from dunshire import *
 910             >>> from dunshire.matrices import specnorm
 911             >>> L = [[1,2],[3,4]]
 912             >>> K = NonnegativeOrthant(2)
 913             >>> e1 = [1,1]
 914             >>> e2 = e1
 915             >>> SLG = SymmetricLinearGame(L,K,e1,e2)
 916             >>> specnorm(SLG.L()) == SLG._L_specnorm()
 917             True
 918
 919         """
 920         if self._L_specnorm_value is None:
 921             self._L_specnorm_value = specnorm(self.L())
 922         return self._L_specnorm_value
 923
 924
 925     def tolerance_scale(self, solution):
 926         r"""
 927
 928         Return a scaling factor that should be applied to
 929         :const:`dunshire.options.ABS_TOL` for this game.
 930
 931         When performing certain comparisons, the default tolerance
 932         :const:`dunshire.options.ABS_TOL` may not be appropriate. For
 933         example, if we expect ``x`` and ``y`` to be within
 934         :const:`dunshire.options.ABS_TOL` of each other, than the inner
 935         product of ``L*x`` and ``y`` can be as far apart as the spectral
 936         norm of ``L`` times the sum of the norms of ``x`` and
 937         ``y``. Such a comparison is made in :meth:`solution`, and in
 938         many of our unit tests.
 939
 940         The returned scaling factor found from the inner product
 941         mentioned above is
 942
 943         .. math::
 944
 945             \left\lVert L \right\rVert_{2}
 946             \left( \left\lVert \bar{x} \right\rVert
 947                    + \left\lVert \bar{y} \right\rVert
 948             \right),
 949
 950         where :math:`\bar{x}` and :math:`\bar{y}` are optimal solutions
 951         for players one and two respectively. This scaling factor is not
 952         formally justified, but attempting anything smaller leads to
 953         test failures.
 954
 955         .. warning::
 956
 957             Optimal solutions are not unique, so the scaling factor
 958             obtained from ``solution`` may not work when comparing other
 959             solutions.
 960
 961         Parameters
 962         ----------
 963
 964         solution : Solution
 965             A solution of this game, used to obtain the norms of the
 966             optimal strategies.
 967
 968         Returns
 969         -------
 970
 971         float
 972             A scaling factor to be multiplied by
 973             :const:`dunshire.options.ABS_TOL` when
 974             making comparisons involving solutions of this game.
 975
 976         Examples
 977         --------
 978
 979         The spectral norm of ``L`` in this case is around ``5.464``, and
 980         the optimal strategies both have norm one, so we expect the
 981         tolerance scale to be somewhere around ``2 * 5.464``, or
 982         ``10.929``::
 983
 984             >>> from dunshire import *
 985             >>> L = [[1,2],[3,4]]
 986             >>> K = NonnegativeOrthant(2)
 987             >>> e1 = [1,1]
 988             >>> e2 = e1
 989             >>> SLG = SymmetricLinearGame(L,K,e1,e2)
 990             >>> SLG.tolerance_scale(SLG.solution())
 991             10.929...
 992
 993         """
 994         norm_p1_opt = norm(solution.player1_optimal())
 995         norm_p2_opt = norm(solution.player2_optimal())
 996         scale = self._L_specnorm()*(norm_p1_opt + norm_p2_opt)
 997
 998         # Don't return anything smaller than 1... we can't go below
 999         # out "minimum tolerance."
1000         return max(1, scale)
1001
1002
1003     def solution(self):
1004         """
1005         Solve this linear game and return a :class:`Solution`.
1006
1007         Returns
1008         -------
1009
1010         Solution
1011             A :class:`Solution` object describing the game's value and
1012             the optimal strategies of both players.
1013
1014         Raises
1015         ------
1016         GameUnsolvableException
1017             If the game could not be solved (if an optimal solution to its
1018             associated cone program was not found).
1019
1020         PoorScalingException
1021             If the game could not be solved because CVXOPT crashed while
1022             trying to take the square root of a negative number.
1023
1024         Examples
1025         --------
1026
1027         This example is computed in Gowda and Ravindran in the section
1028         "The value of a Z-transformation"::
1029
1030             >>> from dunshire import *
1031             >>> K = NonnegativeOrthant(3)
1032             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
1033             >>> e1 = [1,1,1]
1034             >>> e2 = [1,1,1]
1035             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1036             >>> print(SLG.solution())
1037             Game value: -6.172...
1038             Player 1 optimal:
1039               [0.551...]
1040               [0.000...]
1041               [0.448...]
1042             Player 2 optimal:
1043               [0.448...]
1044               [0.000...]
1045               [0.551...]
1046
1047         The value of the following game can be computed using the fact
1048         that the identity is invertible::
1049
1050             >>> from dunshire import *
1051             >>> K = NonnegativeOrthant(3)
1052             >>> L = [[1,0,0],[0,1,0],[0,0,1]]
1053             >>> e1 = [1,2,3]
1054             >>> e2 = [4,5,6]
1055             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1056             >>> print(SLG.solution())
1057             Game value: 0.031...
1058             Player 1 optimal:
1059               [0.031...]
1060               [0.062...]
1061               [0.093...]
1062             Player 2 optimal:
1063               [0.125...]
1064               [0.156...]
1065               [0.187...]
1066
1067         This is another Gowda/Ravindran example that is supposed to have
1068         a negative game value::
1069
1070             >>> from dunshire import *
1071             >>> from dunshire.options import ABS_TOL
1072             >>> L = [[1, -2], [-2, 1]]
1073             >>> K = NonnegativeOrthant(2)
1074             >>> e1 = [1, 1]
1075             >>> e2 = e1
1076             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1077             >>> SLG.solution().game_value() < -ABS_TOL
1078             True
1079
1080         The following two games are problematic numerically, but we
1081         should be able to solve them::
1082
1083             >>> from dunshire import *
1084             >>> L = [[-0.95237953890954685221, 1.83474556206462535712],
1085             ...      [ 1.30481749924621448500, 1.65278664543326403447]]
1086             >>> K = NonnegativeOrthant(2)
1087             >>> e1 = [0.95477167524644313001, 0.63270781756540095397]
1088             >>> e2 = [0.39633793037154141370, 0.10239281495640320530]
1089             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1090             >>> print(SLG.solution())
1091             Game value: 18.767...
1092             Player 1 optimal:
1093               [0.000...]
1094               [9.766...]
1095             Player 2 optimal:
1096               [1.047...]
1097               [0.000...]
1098
1099         ::
1100
1101             >>> from dunshire import *
1102             >>> L = [[1.54159395026049472754, 2.21344728574316684799],
1103             ...      [1.33147433507846657541, 1.17913616272988108769]]
1104             >>> K = NonnegativeOrthant(2)
1105             >>> e1 = [0.39903040089404784307, 0.12377403622479113410]
1106             >>> e2 = [0.15695181142215544612, 0.85527381344651265405]
1107             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1108             >>> print(SLG.solution())
1109             Game value: 24.614...
1110             Player 1 optimal:
1111               [6.371...]
1112               [0.000...]
1113             Player 2 optimal:
1114               [2.506...]
1115               [0.000...]
1116
1117         This is another one that was difficult numerically, and caused
1118         trouble even after we fixed the first two::
1119
1120             >>> from dunshire import *
1121             >>> L = [[57.22233908627052301199, 41.70631373437460354126],
1122             ...      [83.04512571985074487202, 57.82581810406928468637]]
1123             >>> K = NonnegativeOrthant(2)
1124             >>> e1 = [7.31887017043399268346, 0.89744171905822367474]
1125             >>> e2 = [0.11099824781179848388, 6.12564670639315345113]
1126             >>> SLG = SymmetricLinearGame(L,K,e1,e2)
1127             >>> print(SLG.solution())
1128             Game value: 70.437...
1129             Player 1 optimal:
1130               [9.009...]
1131               [0.000...]
1132             Player 2 optimal:
1133               [0.136...]
1134               [0.000...]
1135
1136         And finally, here's one that returns an "optimal" solution, but
1137         whose primal/dual objective function values are far apart::
1138
1139             >>> from dunshire import *
1140             >>> L = [[ 6.49260076597376212248, -0.60528030227678542019],
1141             ...      [ 2.59896077096751731972, -0.97685530240286766457]]
1142             >>> K = IceCream(2)
1143             >>> e1 = [1, 0.43749513972645248661]
1144             >>> e2 = [1, 0.46008379832200291260]
1145             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1146             >>> print(SLG.solution())
1147             Game value: 11.596...
1148             Player 1 optimal:
1149               [ 1.852...]
1150               [-1.852...]
1151             Player 2 optimal:
1152               [ 1.777...]
1153               [-1.777...]
1154
1155         """
1156         try:
1157             opts = {'show_progress': False}
1158             soln_dict = solvers.conelp(self.c(),
1159                                        self.G(),
1160                                        self.h(),
1161                                        self.C().cvxopt_dims(),
1162                                        self.A(),
1163                                        self.b(),
1164                                        primalstart=self.player1_start(),
1165                                        dualstart=self.player2_start(),
1166                                        options=opts)
1167         except ValueError as error:
1168             if str(error) == 'math domain error':
1169                 # Oops, CVXOPT tried to take the square root of a
1170                 # negative number. Report some details about the game
1171                 # rather than just the underlying CVXOPT crash.
1172                 printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1173                 raise PoorScalingException(self)
1174             else:
1175                 raise error
1176
1177         # The optimal strategies are named ``p`` and ``q`` in the
1178         # background documentation, and we need to extract them from
1179         # the CVXOPT ``x`` and ``z`` variables. The objective values
1180         # :math:`nu` and :math:`omega` can also be found in the CVXOPT
1181         # ``x`` and ``y`` variables; however, they're stored
1182         # conveniently as separate entries in the solution dictionary.
1183         p1_value = -soln_dict['primal objective']
1184         p2_value = -soln_dict['dual objective']
1185         p1_optimal = soln_dict['x'][1:]
1186         p2_optimal = soln_dict['z'][self.dimension():]
1187
1188         # The "status" field contains "optimal" if everything went
1189         # according to plan. Other possible values are "primal
1190         # infeasible", "dual infeasible", "unknown", all of which mean
1191         # we didn't get a solution.
1192         #
1193         # The "infeasible" ones are the worst, since they indicate
1194         # that CVXOPT is convinced the problem is infeasible (and that
1195         # cannot happen).
1196         if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
1197             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1198             raise GameUnsolvableException(self, soln_dict)
1199
1200         # For the game value, we could use any of:
1201         #
1202         #   * p1_value
1203         #   * p2_value
1204         #   * (p1_value + p2_value)/2
1205         #   * the game payoff
1206         #
1207         # We want the game value to be the payoff, however, so it
1208         # makes the most sense to just use that, even if it means we
1209         # can't test the fact that p1_value/p2_value are close to the
1210         # payoff.
1211         payoff = self.payoff(p1_optimal, p2_optimal)
1212         soln = Solution(payoff, p1_optimal, p2_optimal)
1213
1214         # The "optimal" and "unknown" results, we actually treat the
1215         # same. Even if CVXOPT bails out due to numerical difficulty,
1216         # it will have some candidate points in mind. If those
1217         # candidates are good enough, we take them. We do the same
1218         # check for "optimal" results.
1219         #
1220         # First we check that the primal/dual objective values are
1221         # close enough because otherwise CVXOPT might return "unknown"
1222         # and give us two points in the cone that are nowhere near
1223         # optimal. And in fact, we need to ensure that they're close
1224         # for "optimal" results, too, because we need to know how
1225         # lenient to be in our testing.
1226         #
1227         if abs(p1_value - p2_value) > self.tolerance_scale(soln)*ABS_TOL:
1228             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1229             raise GameUnsolvableException(self, soln_dict)
1230
1231         # And we also check that the points it gave us belong to the
1232         # cone, just in case...
1233         if (p1_optimal not in self._K) or (p2_optimal not in self._K):
1234             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1235             raise GameUnsolvableException(self, soln_dict)
1236
1237         return soln
1238
1239
1240     def condition(self):
1241         r"""
1242         Return the condition number of this game.
1243
1244         In the CVXOPT construction of this game, two matrices ``G`` and
1245         ``A`` appear. When those matrices are nasty, numerical problems
1246         can show up. We define the condition number of this game to be
1247         the average of the condition numbers of ``G`` and ``A`` in the
1248         CVXOPT construction. If the condition number of this game is
1249         high, you can problems like :class:`PoorScalingException`.
1250
1251         Random testing shows that a condition number of around ``125``
1252         is about the best that we can solve reliably. However, the
1253         failures are intermittent, and you may get lucky with an
1254         ill-conditioned game.
1255
1256         Returns
1257         -------
1258
1259         float
1260             A real number greater than or equal to one that measures how
1261             bad this game is numerically.
1262
1263         Examples
1264         --------
1265
1266         >>> from dunshire import *
1267         >>> K = NonnegativeOrthant(1)
1268         >>> L = [[1]]
1269         >>> e1 = [1]
1270         >>> e2 = e1
1271         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1272         >>> SLG.condition()
1273         1.809...
1274
1275         """
1276         return (condition_number(self.G()) + condition_number(self.A()))/2
1277
1278
1279     def dual(self):
1280         r"""
1281         Return the dual game to this game.
1282
1283         If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
1284         then its dual is :math:`G^{*} =
1285         \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
1286         is symmetric, :math:`K^{*} = K`.
1287
1288         Examples
1289         --------
1290
1291             >>> from dunshire import *
1292             >>> K = NonnegativeOrthant(3)
1293             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
1294             >>> e1 = [1,1,1]
1295             >>> e2 = [1,2,3]
1296             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1297             >>> print(SLG.dual())
1298             The linear game (L, K, e1, e2) where
1299               L = [  1  -1 -12]
1300                   [ -5   2 -15]
1301                   [-15  -3   1],
1302               K = Nonnegative orthant in the real 3-space,
1303               e1 = [ 1]
1304                    [ 2]
1305                    [ 3],
1306               e2 = [ 1]
1307                    [ 1]
1308                    [ 1]
1309
1310         """
1311         # We pass ``self.L()`` right back into the constructor, because
1312         # it will be transposed there. And keep in mind that ``self._K``
1313         # is its own dual.
1314         return SymmetricLinearGame(self.L(),
1315                                    self.K(),
1316                                    self.e2(),
1317                                    self.e1())