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