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         """
 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         """
 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 ``L`` and cache it.
 852         """
 853         if self._L_specnorm_value is None:
 854             self._L_specnorm_value = specnorm(self.L())
 855         return self._L_specnorm_value
 856
 857     def tolerance_scale(self, solution):
 858         # Don't return anything smaller than 1... we can't go below
 859         # out "minimum tolerance."
 860         norm_p1_opt = norm(solution.player1_optimal())
 861         norm_p2_opt = norm(solution.player2_optimal())
 862         scale = self._L_specnorm()*(norm_p1_opt + norm_p2_opt)
 863         return max(1, scale/2.0)
 864
 865
 866     def solution(self):
 867         """
 868         Solve this linear game and return a :class:`Solution`.
 869
 870         Returns
 871         -------
 872
 873         :class:`Solution`
 874             A :class:`Solution` object describing the game's value and
 875             the optimal strategies of both players.
 876
 877         Raises
 878         ------
 879         GameUnsolvableException
 880             If the game could not be solved (if an optimal solution to its
 881             associated cone program was not found).
 882
 883         PoorScalingException
 884             If the game could not be solved because CVXOPT crashed while
 885             trying to take the square root of a negative number.
 886
 887         Examples
 888         --------
 889
 890         This example is computed in Gowda and Ravindran in the section
 891         "The value of a Z-transformation"::
 892
 893             >>> from dunshire import *
 894             >>> K = NonnegativeOrthant(3)
 895             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 896             >>> e1 = [1,1,1]
 897             >>> e2 = [1,1,1]
 898             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 899             >>> print(SLG.solution())
 900             Game value: -6.172...
 901             Player 1 optimal:
 902               [0.551...]
 903               [0.000...]
 904               [0.448...]
 905             Player 2 optimal:
 906               [0.448...]
 907               [0.000...]
 908               [0.551...]
 909
 910         The value of the following game can be computed using the fact
 911         that the identity is invertible::
 912
 913             >>> from dunshire import *
 914             >>> K = NonnegativeOrthant(3)
 915             >>> L = [[1,0,0],[0,1,0],[0,0,1]]
 916             >>> e1 = [1,2,3]
 917             >>> e2 = [4,5,6]
 918             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 919             >>> print(SLG.solution())
 920             Game value: 0.031...
 921             Player 1 optimal:
 922               [0.031...]
 923               [0.062...]
 924               [0.093...]
 925             Player 2 optimal:
 926               [0.125...]
 927               [0.156...]
 928               [0.187...]
 929
 930         This is another Gowda/Ravindran example that is supposed to have
 931         a negative game value::
 932
 933             >>> from dunshire import *
 934             >>> from dunshire.options import ABS_TOL
 935             >>> L = [[1, -2], [-2, 1]]
 936             >>> K = NonnegativeOrthant(2)
 937             >>> e1 = [1, 1]
 938             >>> e2 = e1
 939             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 940             >>> SLG.solution().game_value() < -ABS_TOL
 941             True
 942
 943         The following two games are problematic numerically, but we
 944         should be able to solve them::
 945
 946             >>> from dunshire import *
 947             >>> L = [[-0.95237953890954685221, 1.83474556206462535712],
 948             ...      [ 1.30481749924621448500, 1.65278664543326403447]]
 949             >>> K = NonnegativeOrthant(2)
 950             >>> e1 = [0.95477167524644313001, 0.63270781756540095397]
 951             >>> e2 = [0.39633793037154141370, 0.10239281495640320530]
 952             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 953             >>> print(SLG.solution())
 954             Game value: 18.767...
 955             Player 1 optimal:
 956               [0.000...]
 957               [9.766...]
 958             Player 2 optimal:
 959               [1.047...]
 960               [0.000...]
 961
 962         ::
 963
 964             >>> from dunshire import *
 965             >>> L = [[1.54159395026049472754, 2.21344728574316684799],
 966             ...      [1.33147433507846657541, 1.17913616272988108769]]
 967             >>> K = NonnegativeOrthant(2)
 968             >>> e1 = [0.39903040089404784307, 0.12377403622479113410]
 969             >>> e2 = [0.15695181142215544612, 0.85527381344651265405]
 970             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 971             >>> print(SLG.solution())
 972             Game value: 24.614...
 973             Player 1 optimal:
 974               [6.371...]
 975               [0.000...]
 976             Player 2 optimal:
 977               [2.506...]
 978               [0.000...]
 979
 980         This is another one that was difficult numerically, and caused
 981         trouble even after we fixed the first two::
 982
 983             >>> from dunshire import *
 984             >>> L = [[57.22233908627052301199, 41.70631373437460354126],
 985             ...      [83.04512571985074487202, 57.82581810406928468637]]
 986             >>> K = NonnegativeOrthant(2)
 987             >>> e1 = [7.31887017043399268346, 0.89744171905822367474]
 988             >>> e2 = [0.11099824781179848388, 6.12564670639315345113]
 989             >>> SLG = SymmetricLinearGame(L,K,e1,e2)
 990             >>> print(SLG.solution())
 991             Game value: 70.437...
 992             Player 1 optimal:
 993               [9.009...]
 994               [0.000...]
 995             Player 2 optimal:
 996               [0.136...]
 997               [0.000...]
 998
 999         And finally, here's one that returns an "optimal" solution, but
1000         whose primal/dual objective function values are far apart::
1001
1002             >>> from dunshire import *
1003             >>> L = [[ 6.49260076597376212248, -0.60528030227678542019],
1004             ...      [ 2.59896077096751731972, -0.97685530240286766457]]
1005             >>> K = IceCream(2)
1006             >>> e1 = [1, 0.43749513972645248661]
1007             >>> e2 = [1, 0.46008379832200291260]
1008             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1009             >>> print(SLG.solution())
1010             Game value: 11.596...
1011             Player 1 optimal:
1012               [ 1.852...]
1013               [-1.852...]
1014             Player 2 optimal:
1015               [ 1.777...]
1016               [-1.777...]
1017
1018         """
1019         try:
1020             opts = {'show_progress': False}
1021             soln_dict = solvers.conelp(self._c(),
1022                                        self._G(),
1023                                        self._h(),
1024                                        self.C().cvxopt_dims(),
1025                                        self.A(),
1026                                        self.b(),
1027                                        primalstart=self.player1_start(),
1028                                        dualstart=self.player2_start(),
1029                                        options=opts)
1030         except ValueError as error:
1031             if str(error) == 'math domain error':
1032                 # Oops, CVXOPT tried to take the square root of a
1033                 # negative number. Report some details about the game
1034                 # rather than just the underlying CVXOPT crash.
1035                 printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1036                 raise PoorScalingException(self)
1037             else:
1038                 raise error
1039
1040         # The optimal strategies are named ``p`` and ``q`` in the
1041         # background documentation, and we need to extract them from
1042         # the CVXOPT ``x`` and ``z`` variables. The objective values
1043         # :math:`nu` and :math:`omega` can also be found in the CVXOPT
1044         # ``x`` and ``y`` variables; however, they're stored
1045         # conveniently as separate entries in the solution dictionary.
1046         p1_value = -soln_dict['primal objective']
1047         p2_value = -soln_dict['dual objective']
1048         p1_optimal = soln_dict['x'][1:]
1049         p2_optimal = soln_dict['z'][self.dimension():]
1050
1051         # The "status" field contains "optimal" if everything went
1052         # according to plan. Other possible values are "primal
1053         # infeasible", "dual infeasible", "unknown", all of which mean
1054         # we didn't get a solution.
1055         #
1056         # The "infeasible" ones are the worst, since they indicate
1057         # that CVXOPT is convinced the problem is infeasible (and that
1058         # cannot happen).
1059         if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
1060             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1061             raise GameUnsolvableException(self, soln_dict)
1062
1063         # For the game value, we could use any of:
1064         #
1065         #   * p1_value
1066         #   * p2_value
1067         #   * (p1_value + p2_value)/2
1068         #   * the game payoff
1069         #
1070         # We want the game value to be the payoff, however, so it
1071         # makes the most sense to just use that, even if it means we
1072         # can't test the fact that p1_value/p2_value are close to the
1073         # payoff.
1074         payoff = self.payoff(p1_optimal, p2_optimal)
1075         soln = Solution(payoff, p1_optimal, p2_optimal)
1076
1077         # The "optimal" and "unknown" results, we actually treat the
1078         # same. Even if CVXOPT bails out due to numerical difficulty,
1079         # it will have some candidate points in mind. If those
1080         # candidates are good enough, we take them. We do the same
1081         # check (perhaps pointlessly so) for "optimal" results.
1082         #
1083         # First we check that the primal/dual objective values are
1084         # close enough (one could be low by ABS_TOL, the other high by
1085         # it) because otherwise CVXOPT might return "unknown" and give
1086         # us two points in the cone that are nowhere near optimal.
1087         #
1088         if abs(p1_value - p2_value) > self.tolerance_scale(soln)*ABS_TOL:
1089             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1090             raise GameUnsolvableException(self, soln_dict)
1091
1092         # And we also check that the points it gave us belong to the
1093         # cone, just in case...
1094         if (p1_optimal not in self._K) or (p2_optimal not in self._K):
1095             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1096             raise GameUnsolvableException(self, soln_dict)
1097
1098         return soln
1099
1100
1101     def condition(self):
1102         r"""
1103         Return the condition number of this game.
1104
1105         In the CVXOPT construction of this game, two matrices ``G`` and
1106         ``A`` appear. When those matrices are nasty, numerical problems
1107         can show up. We define the condition number of this game to be
1108         the average of the condition numbers of ``G`` and ``A`` in the
1109         CVXOPT construction. If the condition number of this game is
1110         high, then you can expect numerical difficulty (such as
1111         :class:`PoorScalingException`).
1112
1113         Returns
1114         -------
1115
1116         float
1117             A real number greater than or equal to one that measures how
1118             bad this game is numerically.
1119
1120         Examples
1121         --------
1122
1123         >>> from dunshire import *
1124         >>> K = NonnegativeOrthant(1)
1125         >>> L = [[1]]
1126         >>> e1 = [1]
1127         >>> e2 = e1
1128         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1129         >>> SLG.condition()
1130         1.809...
1131
1132         """
1133         return (condition_number(self._G()) + condition_number(self.A()))/2
1134
1135
1136     def dual(self):
1137         r"""
1138         Return the dual game to this game.
1139
1140         If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
1141         then its dual is :math:`G^{*} =
1142         \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
1143         is symmetric, :math:`K^{*} = K`.
1144
1145         Examples
1146         --------
1147
1148             >>> from dunshire import *
1149             >>> K = NonnegativeOrthant(3)
1150             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
1151             >>> e1 = [1,1,1]
1152             >>> e2 = [1,2,3]
1153             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1154             >>> print(SLG.dual())
1155             The linear game (L, K, e1, e2) where
1156               L = [  1  -1 -12]
1157                   [ -5   2 -15]
1158                   [-15  -3   1],
1159               K = Nonnegative orthant in the real 3-space,
1160               e1 = [ 1]
1161                    [ 2]
1162                    [ 3],
1163               e2 = [ 1]
1164                    [ 1]
1165                    [ 1],
1166               Condition((L, K, e1, e2)) = 44.476...
1167
1168         """
1169         # We pass ``self.L()`` right back into the constructor, because
1170         # it will be transposed there. And keep in mind that ``self._K``
1171         # is its own dual.
1172         return SymmetricLinearGame(self.L(),
1173                                    self.K(),
1174                                    self.e2(),
1175                                    self.e1())