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         """
 581         Return the matrix ``A`` used in our CVXOPT construction.
 582
 583         This matrix ``A``  appears on the right-hand side of ``Ax = b``
 584         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 ``G`` appears on the left-hand side of ``Gx + s = h``
 621         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         """
 660         Return the vector ``c`` used in our CVXOPT construction.
 661
 662         The column vector ``c``  appears in the objective function
 663         value ``<c,x>`` in the statement of the CVXOPT conelp program.
 664
 665         .. warning::
 666
 667             It is not safe to cache any of the matrices passed to
 668             CVXOPT, because it can clobber them.
 669
 670         Returns
 671         -------
 672
 673         matrix
 674             A ``self.dimension()``-by-``1`` column vector.
 675
 676         Examples
 677         --------
 678
 679             >>> from dunshire import *
 680             >>> K = NonnegativeOrthant(3)
 681             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 682             >>> e1 = [1,2,3]
 683             >>> e2 = [1,1,1]
 684             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 685             >>> print(SLG.c())
 686             [-1.0000000]
 687             [ 0.0000000]
 688             [ 0.0000000]
 689             [ 0.0000000]
 690             <BLANKLINE>
 691
 692         """
 693         return matrix([-1, self._zero()])
 694
 695
 696     def C(self):
 697         """
 698         Return the cone ``C`` used in our CVXOPT construction.
 699
 700         The cone ``C`` is the cone over which the conelp program takes
 701         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 ``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 normalize
 813         :meth:`e2`, then you get a point in :meth:`K` that makes a unit
 814         inner product with :meth:`e2`. We then get to choose the primal
 815         objective function value such that the constraint involving
 816         :meth:`L` is satisfied.
 817         """
 818         p = self.e2() / (norm(self.e2()) ** 2)
 819         dist = self.K().ball_radius(self.e1())
 820         nu = - self._L_specnorm()/(dist*norm(self.e2()))
 821         x = matrix([nu, p], (self.dimension() + 1, 1))
 822         s = - self.G()*x
 823
 824         return {'x': x, 's': s}
 825
 826
 827     def player2_start(self):
 828         """
 829         Return a feasible starting point for player two.
 830         """
 831         q = self.e1() / (norm(self.e1()) ** 2)
 832         dist = self.K().ball_radius(self.e2())
 833         omega = self._L_specnorm()/(dist*norm(self.e1()))
 834         y = matrix([omega])
 835         z2 = q
 836         z1 = y*self.e2() - self.L().trans()*z2
 837         z = matrix([z1, z2], (self.dimension()*2, 1))
 838
 839         return {'y': y, 'z': z}
 840
 841
 842     def _L_specnorm(self):
 843         """
 844         Compute the spectral norm of :meth:`L` and cache it.
 845
 846         The spectral norm of the matrix :meth:`L` is used in a few
 847         places. Since it can be expensive to compute, we want to cache
 848         its value. That is not possible in :func:`specnorm`, which lies
 849         outside of a class, so this is the place to do it.
 850
 851         Returns
 852         -------
 853
 854         float
 855             A nonnegative real number; the largest singular value of
 856             the matrix :meth:`L`.
 857
 858         """
 859         if self._L_specnorm_value is None:
 860             self._L_specnorm_value = specnorm(self.L())
 861         return self._L_specnorm_value
 862
 863
 864     def tolerance_scale(self, solution):
 865         r"""
 866         Return a scaling factor that should be applied to ``ABS_TOL``
 867         for this game.
 868
 869         When performing certain comparisons, the default tolernace
 870         ``ABS_TOL`` may not be appropriate. For example, if we expect
 871         ``x`` and ``y`` to be within ``ABS_TOL`` of each other, than the
 872         inner product of ``L*x`` and ``y`` can be as far apart as the
 873         spectral norm of ``L`` times the sum of the norms of ``x`` and
 874         ``y``. Such a comparison is made in :meth:`solution`, and in
 875         many of our unit tests.
 876
 877         The returned scaling factor found from the inner product mentioned
 878         above is
 879
 880         .. math::
 881
 882             \left\lVert L \right\rVert_{2}
 883             \left( \left\lVert \bar{x} \right\rVert
 884                    + \left\lVert \bar{y} \right\rVert
 885             \right),
 886
 887         where :math:`\bar{x}` and :math:`\bar{y}` are optimal solutions
 888         for players one and two respectively. This scaling factor is not
 889         formally justified, but attempting anything smaller leads to
 890         test failures.
 891
 892         .. warning::
 893
 894             Optimal solutions are not unique, so the scaling factor
 895             obtained from ``solution`` may not work when comparing other
 896             solutions.
 897
 898         Parameters
 899         ----------
 900
 901         solution : Solution
 902             A solution of this game, used to obtain the norms of the
 903             optimal strategies.
 904
 905         Returns
 906         -------
 907
 908         float
 909             A scaling factor to be multiplied by ``ABS_TOL`` when
 910             making comparisons involving solutions of this game.
 911
 912         """
 913         norm_p1_opt = norm(solution.player1_optimal())
 914         norm_p2_opt = norm(solution.player2_optimal())
 915         scale = self._L_specnorm()*(norm_p1_opt + norm_p2_opt)
 916
 917         # Don't return anything smaller than 1... we can't go below
 918         # out "minimum tolerance."
 919         return max(1, scale)
 920
 921
 922     def solution(self):
 923         """
 924         Solve this linear game and return a :class:`Solution`.
 925
 926         Returns
 927         -------
 928
 929         :class:`Solution`
 930             A :class:`Solution` object describing the game's value and
 931             the optimal strategies of both players.
 932
 933         Raises
 934         ------
 935         GameUnsolvableException
 936             If the game could not be solved (if an optimal solution to its
 937             associated cone program was not found).
 938
 939         PoorScalingException
 940             If the game could not be solved because CVXOPT crashed while
 941             trying to take the square root of a negative number.
 942
 943         Examples
 944         --------
 945
 946         This example is computed in Gowda and Ravindran in the section
 947         "The value of a Z-transformation"::
 948
 949             >>> from dunshire import *
 950             >>> K = NonnegativeOrthant(3)
 951             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 952             >>> e1 = [1,1,1]
 953             >>> e2 = [1,1,1]
 954             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 955             >>> print(SLG.solution())
 956             Game value: -6.172...
 957             Player 1 optimal:
 958               [0.551...]
 959               [0.000...]
 960               [0.448...]
 961             Player 2 optimal:
 962               [0.448...]
 963               [0.000...]
 964               [0.551...]
 965
 966         The value of the following game can be computed using the fact
 967         that the identity is invertible::
 968
 969             >>> from dunshire import *
 970             >>> K = NonnegativeOrthant(3)
 971             >>> L = [[1,0,0],[0,1,0],[0,0,1]]
 972             >>> e1 = [1,2,3]
 973             >>> e2 = [4,5,6]
 974             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 975             >>> print(SLG.solution())
 976             Game value: 0.031...
 977             Player 1 optimal:
 978               [0.031...]
 979               [0.062...]
 980               [0.093...]
 981             Player 2 optimal:
 982               [0.125...]
 983               [0.156...]
 984               [0.187...]
 985
 986         This is another Gowda/Ravindran example that is supposed to have
 987         a negative game value::
 988
 989             >>> from dunshire import *
 990             >>> from dunshire.options import ABS_TOL
 991             >>> L = [[1, -2], [-2, 1]]
 992             >>> K = NonnegativeOrthant(2)
 993             >>> e1 = [1, 1]
 994             >>> e2 = e1
 995             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 996             >>> SLG.solution().game_value() < -ABS_TOL
 997             True
 998
 999         The following two games are problematic numerically, but we
1000         should be able to solve them::
1001
1002             >>> from dunshire import *
1003             >>> L = [[-0.95237953890954685221, 1.83474556206462535712],
1004             ...      [ 1.30481749924621448500, 1.65278664543326403447]]
1005             >>> K = NonnegativeOrthant(2)
1006             >>> e1 = [0.95477167524644313001, 0.63270781756540095397]
1007             >>> e2 = [0.39633793037154141370, 0.10239281495640320530]
1008             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1009             >>> print(SLG.solution())
1010             Game value: 18.767...
1011             Player 1 optimal:
1012               [0.000...]
1013               [9.766...]
1014             Player 2 optimal:
1015               [1.047...]
1016               [0.000...]
1017
1018         ::
1019
1020             >>> from dunshire import *
1021             >>> L = [[1.54159395026049472754, 2.21344728574316684799],
1022             ...      [1.33147433507846657541, 1.17913616272988108769]]
1023             >>> K = NonnegativeOrthant(2)
1024             >>> e1 = [0.39903040089404784307, 0.12377403622479113410]
1025             >>> e2 = [0.15695181142215544612, 0.85527381344651265405]
1026             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1027             >>> print(SLG.solution())
1028             Game value: 24.614...
1029             Player 1 optimal:
1030               [6.371...]
1031               [0.000...]
1032             Player 2 optimal:
1033               [2.506...]
1034               [0.000...]
1035
1036         This is another one that was difficult numerically, and caused
1037         trouble even after we fixed the first two::
1038
1039             >>> from dunshire import *
1040             >>> L = [[57.22233908627052301199, 41.70631373437460354126],
1041             ...      [83.04512571985074487202, 57.82581810406928468637]]
1042             >>> K = NonnegativeOrthant(2)
1043             >>> e1 = [7.31887017043399268346, 0.89744171905822367474]
1044             >>> e2 = [0.11099824781179848388, 6.12564670639315345113]
1045             >>> SLG = SymmetricLinearGame(L,K,e1,e2)
1046             >>> print(SLG.solution())
1047             Game value: 70.437...
1048             Player 1 optimal:
1049               [9.009...]
1050               [0.000...]
1051             Player 2 optimal:
1052               [0.136...]
1053               [0.000...]
1054
1055         And finally, here's one that returns an "optimal" solution, but
1056         whose primal/dual objective function values are far apart::
1057
1058             >>> from dunshire import *
1059             >>> L = [[ 6.49260076597376212248, -0.60528030227678542019],
1060             ...      [ 2.59896077096751731972, -0.97685530240286766457]]
1061             >>> K = IceCream(2)
1062             >>> e1 = [1, 0.43749513972645248661]
1063             >>> e2 = [1, 0.46008379832200291260]
1064             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1065             >>> print(SLG.solution())
1066             Game value: 11.596...
1067             Player 1 optimal:
1068               [ 1.852...]
1069               [-1.852...]
1070             Player 2 optimal:
1071               [ 1.777...]
1072               [-1.777...]
1073
1074         """
1075         try:
1076             opts = {'show_progress': False}
1077             soln_dict = solvers.conelp(self.c(),
1078                                        self.G(),
1079                                        self.h(),
1080                                        self.C().cvxopt_dims(),
1081                                        self.A(),
1082                                        self.b(),
1083                                        primalstart=self.player1_start(),
1084                                        dualstart=self.player2_start(),
1085                                        options=opts)
1086         except ValueError as error:
1087             if str(error) == 'math domain error':
1088                 # Oops, CVXOPT tried to take the square root of a
1089                 # negative number. Report some details about the game
1090                 # rather than just the underlying CVXOPT crash.
1091                 printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1092                 raise PoorScalingException(self)
1093             else:
1094                 raise error
1095
1096         # The optimal strategies are named ``p`` and ``q`` in the
1097         # background documentation, and we need to extract them from
1098         # the CVXOPT ``x`` and ``z`` variables. The objective values
1099         # :math:`nu` and :math:`omega` can also be found in the CVXOPT
1100         # ``x`` and ``y`` variables; however, they're stored
1101         # conveniently as separate entries in the solution dictionary.
1102         p1_value = -soln_dict['primal objective']
1103         p2_value = -soln_dict['dual objective']
1104         p1_optimal = soln_dict['x'][1:]
1105         p2_optimal = soln_dict['z'][self.dimension():]
1106
1107         # The "status" field contains "optimal" if everything went
1108         # according to plan. Other possible values are "primal
1109         # infeasible", "dual infeasible", "unknown", all of which mean
1110         # we didn't get a solution.
1111         #
1112         # The "infeasible" ones are the worst, since they indicate
1113         # that CVXOPT is convinced the problem is infeasible (and that
1114         # cannot happen).
1115         if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
1116             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1117             raise GameUnsolvableException(self, soln_dict)
1118
1119         # For the game value, we could use any of:
1120         #
1121         #   * p1_value
1122         #   * p2_value
1123         #   * (p1_value + p2_value)/2
1124         #   * the game payoff
1125         #
1126         # We want the game value to be the payoff, however, so it
1127         # makes the most sense to just use that, even if it means we
1128         # can't test the fact that p1_value/p2_value are close to the
1129         # payoff.
1130         payoff = self.payoff(p1_optimal, p2_optimal)
1131         soln = Solution(payoff, p1_optimal, p2_optimal)
1132
1133         # The "optimal" and "unknown" results, we actually treat the
1134         # same. Even if CVXOPT bails out due to numerical difficulty,
1135         # it will have some candidate points in mind. If those
1136         # candidates are good enough, we take them. We do the same
1137         # check (perhaps pointlessly so) for "optimal" results.
1138         #
1139         # First we check that the primal/dual objective values are
1140         # close enough (one could be low by ABS_TOL, the other high by
1141         # it) because otherwise CVXOPT might return "unknown" and give
1142         # us two points in the cone that are nowhere near optimal.
1143         #
1144         if abs(p1_value - p2_value) > self.tolerance_scale(soln)*ABS_TOL:
1145             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1146             raise GameUnsolvableException(self, soln_dict)
1147
1148         # And we also check that the points it gave us belong to the
1149         # cone, just in case...
1150         if (p1_optimal not in self._K) or (p2_optimal not in self._K):
1151             printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1152             raise GameUnsolvableException(self, soln_dict)
1153
1154         return soln
1155
1156
1157     def condition(self):
1158         r"""
1159         Return the condition number of this game.
1160
1161         In the CVXOPT construction of this game, two matrices ``G`` and
1162         ``A`` appear. When those matrices are nasty, numerical problems
1163         can show up. We define the condition number of this game to be
1164         the average of the condition numbers of ``G`` and ``A`` in the
1165         CVXOPT construction. If the condition number of this game is
1166         high, then you can expect numerical difficulty (such as
1167         :class:`PoorScalingException`).
1168
1169         Returns
1170         -------
1171
1172         float
1173             A real number greater than or equal to one that measures how
1174             bad this game is numerically.
1175
1176         Examples
1177         --------
1178
1179         >>> from dunshire import *
1180         >>> K = NonnegativeOrthant(1)
1181         >>> L = [[1]]
1182         >>> e1 = [1]
1183         >>> e2 = e1
1184         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1185         >>> SLG.condition()
1186         1.809...
1187
1188         """
1189         return (condition_number(self.G()) + condition_number(self.A()))/2
1190
1191
1192     def dual(self):
1193         r"""
1194         Return the dual game to this game.
1195
1196         If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
1197         then its dual is :math:`G^{*} =
1198         \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
1199         is symmetric, :math:`K^{*} = K`.
1200
1201         Examples
1202         --------
1203
1204             >>> from dunshire import *
1205             >>> K = NonnegativeOrthant(3)
1206             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
1207             >>> e1 = [1,1,1]
1208             >>> e2 = [1,2,3]
1209             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1210             >>> print(SLG.dual())
1211             The linear game (L, K, e1, e2) where
1212               L = [  1  -1 -12]
1213                   [ -5   2 -15]
1214                   [-15  -3   1],
1215               K = Nonnegative orthant in the real 3-space,
1216               e1 = [ 1]
1217                    [ 2]
1218                    [ 3],
1219               e2 = [ 1]
1220                    [ 1]
1221                    [ 1]
1222
1223         """
1224         # We pass ``self.L()`` right back into the constructor, because
1225         # it will be transposed there. And keep in mind that ``self._K``
1226         # is its own dual.
1227         return SymmetricLinearGame(self.L(),
1228                                    self.K(),
1229                                    self.e2(),
1230                                    self.e1())