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         Return a scaling factor that should be applied to :const:`ABS_TOL`
 928         for this game.
 929
 930         When performing certain comparisons, the default tolerance
 931         :const:`ABS_TOL` may not be appropriate. For example, if we expect
 932         ``x`` and ``y`` to be within :const:`ABS_TOL` of each other,
 933         than the inner product of ``L*x`` and ``y`` can be as far apart
 934         as the spectral norm of ``L`` times the sum of the norms of
 935         ``x`` and ``y``. Such a comparison is made in :meth:`solution`,
 936         and in many of our unit tests.
 937
 938         The returned scaling factor found from the inner product
 939         mentioned above is
 940
 941         .. math::
 942
 943             \left\lVert L \right\rVert_{2}
 944             \left( \left\lVert \bar{x} \right\rVert
 945                    + \left\lVert \bar{y} \right\rVert
 946             \right),
 947
 948         where :math:`\bar{x}` and :math:`\bar{y}` are optimal solutions
 949         for players one and two respectively. This scaling factor is not
 950         formally justified, but attempting anything smaller leads to
 951         test failures.
 952
 953         .. warning::
 954
 955             Optimal solutions are not unique, so the scaling factor
 956             obtained from ``solution`` may not work when comparing other
 957             solutions.
 958
 959         Parameters
 960         ----------
 961
 962         solution : Solution
 963             A solution of this game, used to obtain the norms of the
 964             optimal strategies.
 965
 966         Returns
 967         -------
 968
 969         float
 970             A scaling factor to be multiplied by :const:`ABS_TOL` when
 971             making comparisons involving solutions of this game.
 972
 973         Examples
 974         --------
 975
 976         The spectral norm of ``L`` in this case is around ``5.464``, and
 977         the optimal strategies both have norm one, so we expect the
 978         tolerance scale to be somewhere around ``2 * 5.464``, or
 979         ``10.929``::
 980
 981             >>> from dunshire import *
 982             >>> L = [[1,2],[3,4]]
 983             >>> K = NonnegativeOrthant(2)
 984             >>> e1 = [1,1]
 985             >>> e2 = e1
 986             >>> SLG = SymmetricLinearGame(L,K,e1,e2)
 987             >>> SLG.tolerance_scale(SLG.solution())
 988             10.929...
 989
 990         """
 991         norm_p1_opt = norm(solution.player1_optimal())
 992         norm_p2_opt = norm(solution.player2_optimal())
 993         scale = self._L_specnorm()*(norm_p1_opt + norm_p2_opt)
 994
 995         # Don't return anything smaller than 1... we can't go below
 996         # out "minimum tolerance."
 997         return max(1, scale)
 998
 999
1000     def solution(self):
1001         """
1002         Solve this linear game and return a :class:`Solution`.
1003
1004         Returns
1005         -------
1006
1007         :class:`Solution`
1008             A :class:`Solution` object describing the game's value and
1009             the optimal strategies of both players.
1010
1011         Raises
1012         ------
1013         GameUnsolvableException
1014             If the game could not be solved (if an optimal solution to its
1015             associated cone program was not found).
1016
1017         PoorScalingException
1018             If the game could not be solved because CVXOPT crashed while
1019             trying to take the square root of a negative number.
1020
1021         Examples
1022         --------
1023
1024         This example is computed in Gowda and Ravindran in the section
1025         "The value of a Z-transformation"::
1026
1027             >>> from dunshire import *
1028             >>> K = NonnegativeOrthant(3)
1029             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
1030             >>> e1 = [1,1,1]
1031             >>> e2 = [1,1,1]
1032             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1033             >>> print(SLG.solution())
1034             Game value: -6.172...
1035             Player 1 optimal:
1036               [0.551...]
1037               [0.000...]
1038               [0.448...]
1039             Player 2 optimal:
1040               [0.448...]
1041               [0.000...]
1042               [0.551...]
1043
1044         The value of the following game can be computed using the fact
1045         that the identity is invertible::
1046
1047             >>> from dunshire import *
1048             >>> K = NonnegativeOrthant(3)
1049             >>> L = [[1,0,0],[0,1,0],[0,0,1]]
1050             >>> e1 = [1,2,3]
1051             >>> e2 = [4,5,6]
1052             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1053             >>> print(SLG.solution())
1054             Game value: 0.031...
1055             Player 1 optimal:
1056               [0.031...]
1057               [0.062...]
1058               [0.093...]
1059             Player 2 optimal:
1060               [0.125...]
1061               [0.156...]
1062               [0.187...]
1063
1064         This is another Gowda/Ravindran example that is supposed to have
1065         a negative game value::
1066
1067             >>> from dunshire import *
1068             >>> from dunshire.options import ABS_TOL
1069             >>> L = [[1, -2], [-2, 1]]
1070             >>> K = NonnegativeOrthant(2)
1071             >>> e1 = [1, 1]
1072             >>> e2 = e1
1073             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1074             >>> SLG.solution().game_value() < -ABS_TOL
1075             True
1076
1077         The following two games are problematic numerically, but we
1078         should be able to solve them::
1079
1080             >>> from dunshire import *
1081             >>> L = [[-0.95237953890954685221, 1.83474556206462535712],
1082             ...      [ 1.30481749924621448500, 1.65278664543326403447]]
1083             >>> K = NonnegativeOrthant(2)
1084             >>> e1 = [0.95477167524644313001, 0.63270781756540095397]
1085             >>> e2 = [0.39633793037154141370, 0.10239281495640320530]
1086             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1087             >>> print(SLG.solution())
1088             Game value: 18.767...
1089             Player 1 optimal:
1090               [0.000...]
1091               [9.766...]
1092             Player 2 optimal:
1093               [1.047...]
1094               [0.000...]
1095
1096         ::
1097
1098             >>> from dunshire import *
1099             >>> L = [[1.54159395026049472754, 2.21344728574316684799],
1100             ...      [1.33147433507846657541, 1.17913616272988108769]]
1101             >>> K = NonnegativeOrthant(2)
1102             >>> e1 = [0.39903040089404784307, 0.12377403622479113410]
1103             >>> e2 = [0.15695181142215544612, 0.85527381344651265405]
1104             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1105             >>> print(SLG.solution())
1106             Game value: 24.614...
1107             Player 1 optimal:
1108               [6.371...]
1109               [0.000...]
1110             Player 2 optimal:
1111               [2.506...]
1112               [0.000...]
1113
1114         This is another one that was difficult numerically, and caused
1115         trouble even after we fixed the first two::
1116
1117             >>> from dunshire import *
1118             >>> L = [[57.22233908627052301199, 41.70631373437460354126],
1119             ...      [83.04512571985074487202, 57.82581810406928468637]]
1120             >>> K = NonnegativeOrthant(2)
1121             >>> e1 = [7.31887017043399268346, 0.89744171905822367474]
1122             >>> e2 = [0.11099824781179848388, 6.12564670639315345113]
1123             >>> SLG = SymmetricLinearGame(L,K,e1,e2)
1124             >>> print(SLG.solution())
1125             Game value: 70.437...
1126             Player 1 optimal:
1127               [9.009...]
1128               [0.000...]
1129             Player 2 optimal:
1130               [0.136...]
1131               [0.000...]
1132
1133         And finally, here's one that returns an "optimal" solution, but
1134         whose primal/dual objective function values are far apart::
1135
1136             >>> from dunshire import *
1137             >>> L = [[ 6.49260076597376212248, -0.60528030227678542019],
1138             ...      [ 2.59896077096751731972, -0.97685530240286766457]]
1139             >>> K = IceCream(2)
1140             >>> e1 = [1, 0.43749513972645248661]
1141             >>> e2 = [1, 0.46008379832200291260]
1142             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1143             >>> print(SLG.solution())
1144             Game value: 11.596...
1145             Player 1 optimal:
1146               [ 1.852...]
1147               [-1.852...]
1148             Player 2 optimal:
1149               [ 1.777...]
1150               [-1.777...]
1151
1152         """
1153         try:
1154             opts = {'show_progress': False}
1155             soln_dict = solvers.conelp(self.c(),
1156                                        self.G(),
1157                                        self.h(),
1158                                        self.C().cvxopt_dims(),
1159                                        self.A(),
1160                                        self.b(),
1161                                        primalstart=self.player1_start(),
1162                                        dualstart=self.player2_start(),
1163                                        options=opts)
1164         except ValueError as error:
1165             if str(error) == 'math domain error':
1166                 # Oops, CVXOPT tried to take the square root of a
1167                 # negative number. Report some details about the game
1168                 # rather than just the underlying CVXOPT crash.
1169                 printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1170                 raise PoorScalingException(self)
1171             else:
1172                 raise error
1173
1174         # The optimal strategies are named ``p`` and ``q`` in the
1175         # background documentation, and we need to extract them from
1176         # the CVXOPT ``x`` and ``z`` variables. The objective values
1177         # :math:`nu` and :math:`omega` can also be found in the CVXOPT
1178         # ``x`` and ``y`` variables; however, they're stored
1179         # conveniently as separate entries in the solution dictionary.
1180         p1_value = -soln_dict['primal objective']
1181         p2_value = -soln_dict['dual objective']
1182         p1_optimal = soln_dict['x'][1:]
1183         p2_optimal = soln_dict['z'][self.dimension():]
1184
1185         # The "status" field contains "optimal" if everything went
1186         # according to plan. Other possible values are "primal
1187         # infeasible", "dual infeasible", "unknown", all of which mean
1188         # we didn't get a solution.
1189         #
1190         # The "infeasible" ones are the worst, since they indicate
1191         # that CVXOPT is convinced the problem is infeasible (and that
1192         # cannot happen).
1193         if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
1194             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1195             raise GameUnsolvableException(self, soln_dict)
1196
1197         # For the game value, we could use any of:
1198         #
1199         #   * p1_value
1200         #   * p2_value
1201         #   * (p1_value + p2_value)/2
1202         #   * the game payoff
1203         #
1204         # We want the game value to be the payoff, however, so it
1205         # makes the most sense to just use that, even if it means we
1206         # can't test the fact that p1_value/p2_value are close to the
1207         # payoff.
1208         payoff = self.payoff(p1_optimal, p2_optimal)
1209         soln = Solution(payoff, p1_optimal, p2_optimal)
1210
1211         # The "optimal" and "unknown" results, we actually treat the
1212         # same. Even if CVXOPT bails out due to numerical difficulty,
1213         # it will have some candidate points in mind. If those
1214         # candidates are good enough, we take them. We do the same
1215         # check for "optimal" results.
1216         #
1217         # First we check that the primal/dual objective values are
1218         # close enough because otherwise CVXOPT might return "unknown"
1219         # and give us two points in the cone that are nowhere near
1220         # optimal. And in fact, we need to ensure that they're close
1221         # for "optimal" results, too, because we need to know how
1222         # lenient to be in our testing.
1223         #
1224         if abs(p1_value - p2_value) > self.tolerance_scale(soln)*ABS_TOL:
1225             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1226             raise GameUnsolvableException(self, soln_dict)
1227
1228         # And we also check that the points it gave us belong to the
1229         # cone, just in case...
1230         if (p1_optimal not in self._K) or (p2_optimal not in self._K):
1231             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1232             raise GameUnsolvableException(self, soln_dict)
1233
1234         return soln
1235
1236
1237     def condition(self):
1238         r"""
1239         Return the condition number of this game.
1240
1241         In the CVXOPT construction of this game, two matrices ``G`` and
1242         ``A`` appear. When those matrices are nasty, numerical problems
1243         can show up. We define the condition number of this game to be
1244         the average of the condition numbers of ``G`` and ``A`` in the
1245         CVXOPT construction. If the condition number of this game is
1246         high, you can problems like :class:`PoorScalingException`.
1247
1248         Random testing shows that a condition number of around ``125``
1249         is about the best that we can solve reliably. However, the
1250         failures are intermittent, and you may get lucky with an
1251         ill-conditioned game.
1252
1253         Returns
1254         -------
1255
1256         float
1257             A real number greater than or equal to one that measures how
1258             bad this game is numerically.
1259
1260         Examples
1261         --------
1262
1263         >>> from dunshire import *
1264         >>> K = NonnegativeOrthant(1)
1265         >>> L = [[1]]
1266         >>> e1 = [1]
1267         >>> e2 = e1
1268         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1269         >>> SLG.condition()
1270         1.809...
1271
1272         """
1273         return (condition_number(self.G()) + condition_number(self.A()))/2
1274
1275
1276     def dual(self):
1277         r"""
1278         Return the dual game to this game.
1279
1280         If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
1281         then its dual is :math:`G^{*} =
1282         \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
1283         is symmetric, :math:`K^{*} = K`.
1284
1285         Examples
1286         --------
1287
1288             >>> from dunshire import *
1289             >>> K = NonnegativeOrthant(3)
1290             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
1291             >>> e1 = [1,1,1]
1292             >>> e2 = [1,2,3]
1293             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1294             >>> print(SLG.dual())
1295             The linear game (L, K, e1, e2) where
1296               L = [  1  -1 -12]
1297                   [ -5   2 -15]
1298                   [-15  -3   1],
1299               K = Nonnegative orthant in the real 3-space,
1300               e1 = [ 1]
1301                    [ 2]
1302                    [ 3],
1303               e2 = [ 1]
1304                    [ 1]
1305                    [ 1]
1306
1307         """
1308         # We pass ``self.L()`` right back into the constructor, because
1309         # it will be transposed there. And keep in mind that ``self._K``
1310         # is its own dual.
1311         return SymmetricLinearGame(self.L(),
1312                                    self.K(),
1313                                    self.e2(),
1314                                    self.e1())