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