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