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