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, IceCream, NonnegativeOrthant
   9 from .errors import GameUnsolvableException, PoorScalingException
  10 from .matrices import (append_col, append_row, condition_number, identity,
  11                        inner_product, norm, specnorm)
  12 from . import options
  13
  14 printing.options['dformat'] = options.FLOAT_FORMAT
  15
  16 class Solution:
  17     """
  18     A representation of the solution of a linear game. It should contain
  19     the value of the game, and both players' strategies.
  20
  21     Examples
  22     --------
  23
  24         >>> print(Solution(10, matrix([1,2]), matrix([3,4])))
  25         Game value: 10.000...
  26         Player 1 optimal:
  27           [ 1]
  28           [ 2]
  29         Player 2 optimal:
  30           [ 3]
  31           [ 4]
  32
  33     """
  34     def __init__(self, game_value, p1_optimal, p2_optimal):
  35         """
  36         Create a new Solution object from a game value and two optimal
  37         strategies for the players.
  38         """
  39         self._game_value = game_value
  40         self._player1_optimal = p1_optimal
  41         self._player2_optimal = p2_optimal
  42
  43     def __str__(self):
  44         """
  45         Return a string describing the solution of a linear game.
  46
  47         The three data that are described are,
  48
  49           * The value of the game.
  50           * The optimal strategy of player one.
  51           * The optimal strategy of player two.
  52
  53         The two optimal strategy vectors are indented by two spaces.
  54         """
  55         tpl = 'Game value: {:.7f}\n' \
  56               'Player 1 optimal:{:s}\n' \
  57               'Player 2 optimal:{:s}'
  58
  59         p1_str = '\n{!s}'.format(self.player1_optimal())
  60         p1_str = '\n  '.join(p1_str.splitlines())
  61         p2_str = '\n{!s}'.format(self.player2_optimal())
  62         p2_str = '\n  '.join(p2_str.splitlines())
  63
  64         return tpl.format(self.game_value(), p1_str, p2_str)
  65
  66
  67     def game_value(self):
  68         """
  69         Return the game value for this solution.
  70
  71         Examples
  72         --------
  73
  74             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  75             >>> s.game_value()
  76             10
  77
  78         """
  79         return self._game_value
  80
  81
  82     def player1_optimal(self):
  83         """
  84         Return player one's optimal strategy in this solution.
  85
  86         Examples
  87         --------
  88
  89             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  90             >>> print(s.player1_optimal())
  91             [ 1]
  92             [ 2]
  93             <BLANKLINE>
  94
  95         """
  96         return self._player1_optimal
  97
  98
  99     def player2_optimal(self):
 100         """
 101         Return player two's optimal strategy in this solution.
 102
 103         Examples
 104         --------
 105
 106             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
 107             >>> print(s.player2_optimal())
 108             [ 3]
 109             [ 4]
 110             <BLANKLINE>
 111
 112         """
 113         return self._player2_optimal
 114
 115
 116 class SymmetricLinearGame:
 117     r"""
 118     A representation of a symmetric linear game.
 119
 120     The data for a symmetric linear game are,
 121
 122       * A "payoff" operator ``L``.
 123       * A symmetric cone ``K``.
 124       * Two points ``e1`` and ``e2`` in the interior of ``K``.
 125
 126     The ambient space is assumed to be the span of ``K``.
 127
 128     With those data understood, the game is played as follows. Players
 129     one and two choose points :math:`x` and :math:`y` respectively, from
 130     their respective strategy sets,
 131
 132     .. math::
 133         \begin{aligned}
 134             \Delta_{1}
 135             &=
 136             \left\{
 137               x \in K \ \middle|\ \left\langle x, e_{2} \right\rangle = 1
 138             \right\}\\
 139             \Delta_{2}
 140             &=
 141             \left\{
 142               y \in K \ \middle|\ \left\langle y, e_{1} \right\rangle = 1
 143             \right\}.
 144         \end{aligned}
 145
 146     Afterwards, a "payout" is computed as :math:`\left\langle
 147     L\left(x\right), y \right\rangle` and is paid to player one out of
 148     player two's pocket. The game is therefore zero sum, and we suppose
 149     that player one would like to guarantee himself the largest minimum
 150     payout possible. That is, player one wishes to,
 151
 152     .. math::
 153         \begin{aligned}
 154             \text{maximize }
 155               &\underset{y \in \Delta_{2}}{\min}\left(
 156                 \left\langle L\left(x\right), y \right\rangle
 157               \right)\\
 158             \text{subject to } & x \in \Delta_{1}.
 159         \end{aligned}
 160
 161     Player two has the simultaneous goal to,
 162
 163     .. math::
 164         \begin{aligned}
 165             \text{minimize }
 166               &\underset{x \in \Delta_{1}}{\max}\left(
 167                 \left\langle L\left(x\right), y \right\rangle
 168               \right)\\
 169             \text{subject to } & y \in \Delta_{2}.
 170         \end{aligned}
 171
 172     These goals obviously conflict (the game is zero sum), but an
 173     existence theorem guarantees at least one optimal min-max solution
 174     from which neither player would like to deviate. This class is
 175     able to find such a solution.
 176
 177     Parameters
 178     ----------
 179
 180     L : list of list of float
 181         A matrix represented as a list of ROWS. This representation
 182         agrees with (for example) SageMath and NumPy, but not with CVXOPT
 183         (whose matrix constructor accepts a list of columns).
 184
 185     K : :class:`SymmetricCone`
 186         The symmetric cone instance over which the game is played.
 187
 188     e1 : iterable float
 189         The interior point of ``K`` belonging to player one; it
 190         can be of any iterable type having the correct length.
 191
 192     e2 : iterable float
 193         The interior point of ``K`` belonging to player two; it
 194         can be of any enumerable type having the correct length.
 195
 196     Raises
 197     ------
 198
 199     ValueError
 200         If either ``e1`` or ``e2`` lie outside of the cone ``K``.
 201
 202     Examples
 203     --------
 204
 205         >>> from dunshire import *
 206         >>> K = NonnegativeOrthant(3)
 207         >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 208         >>> e1 = [1,1,1]
 209         >>> e2 = [1,2,3]
 210         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 211         >>> print(SLG)
 212         The linear game (L, K, e1, e2) where
 213           L = [  1  -5 -15]
 214               [ -1   2  -3]
 215               [-12 -15   1],
 216           K = Nonnegative orthant in the real 3-space,
 217           e1 = [ 1]
 218                [ 1]
 219                [ 1],
 220           e2 = [ 1]
 221                [ 2]
 222                [ 3],
 223           Condition((L, K, e1, e2)) = 31.834...
 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           Condition((L, K, e1, e2)) = 1.707...
 243
 244     The points ``e1`` and ``e2`` can also be passed as some other
 245     enumerable type (of the correct length) without much harm, since
 246     there is no row/column ambiguity::
 247
 248         >>> import cvxopt
 249         >>> import numpy
 250         >>> from dunshire import *
 251         >>> K = NonnegativeOrthant(2)
 252         >>> L = [[1,0],[0,1]]
 253         >>> e1 = cvxopt.matrix([1,1])
 254         >>> e2 = numpy.matrix([1,1])
 255         >>> G = SymmetricLinearGame(L, K, e1, e2)
 256         >>> print(G)
 257         The linear game (L, K, e1, e2) where
 258           L = [ 1  0]
 259               [ 0  1],
 260           K = Nonnegative orthant in the real 2-space,
 261           e1 = [ 1]
 262                [ 1],
 263           e2 = [ 1]
 264                [ 1],
 265           Condition((L, K, e1, e2)) = 1.707...
 266
 267     However, ``L`` will always be intepreted as a list of rows, even
 268     if it is passed as a :class:`cvxopt.base.matrix` which is
 269     otherwise indexed by columns::
 270
 271         >>> import cvxopt
 272         >>> from dunshire import *
 273         >>> K = NonnegativeOrthant(2)
 274         >>> L = [[1,2],[3,4]]
 275         >>> e1 = [1,1]
 276         >>> e2 = e1
 277         >>> G = SymmetricLinearGame(L, K, e1, e2)
 278         >>> print(G)
 279         The linear game (L, K, e1, e2) where
 280           L = [ 1  2]
 281               [ 3  4],
 282           K = Nonnegative orthant in the real 2-space,
 283           e1 = [ 1]
 284                [ 1],
 285           e2 = [ 1]
 286                [ 1],
 287           Condition((L, K, e1, e2)) = 6.073...
 288         >>> L = cvxopt.matrix(L)
 289         >>> print(L)
 290         [ 1  3]
 291         [ 2  4]
 292         <BLANKLINE>
 293         >>> G = SymmetricLinearGame(L, K, e1, e2)
 294         >>> print(G)
 295         The linear game (L, K, e1, e2) where
 296           L = [ 1  2]
 297               [ 3  4],
 298           K = Nonnegative orthant in the real 2-space,
 299           e1 = [ 1]
 300                [ 1],
 301           e2 = [ 1]
 302                [ 1],
 303           Condition((L, K, e1, e2)) = 6.073...
 304
 305     """
 306     def __init__(self, L, K, e1, e2):
 307         """
 308         Create a new SymmetricLinearGame object.
 309         """
 310         self._K = K
 311         self._e1 = matrix(e1, (K.dimension(), 1))
 312         self._e2 = matrix(e2, (K.dimension(), 1))
 313
 314         # Our input ``L`` is indexed by rows but CVXOPT matrices are
 315         # indexed by columns, so we need to transpose the input before
 316         # feeding it to CVXOPT.
 317         self._L = matrix(L, (K.dimension(), K.dimension())).trans()
 318
 319         if not self._e1 in K:
 320             raise ValueError('the point e1 must lie in the interior of K')
 321
 322         if not self._e2 in K:
 323             raise ValueError('the point e2 must lie in the interior of K')
 324
 325
 326
 327     def __str__(self):
 328         """
 329         Return a string representation of this game.
 330         """
 331         tpl = 'The linear game (L, K, e1, e2) where\n' \
 332               '  L = {:s},\n' \
 333               '  K = {!s},\n' \
 334               '  e1 = {:s},\n' \
 335               '  e2 = {:s},\n' \
 336               '  Condition((L, K, e1, e2)) = {:f}.'
 337         indented_L = '\n      '.join(str(self.L()).splitlines())
 338         indented_e1 = '\n       '.join(str(self.e1()).splitlines())
 339         indented_e2 = '\n       '.join(str(self.e2()).splitlines())
 340
 341         return tpl.format(indented_L,
 342                           str(self.K()),
 343                           indented_e1,
 344                           indented_e2,
 345                           self.condition())
 346
 347
 348     def L(self):
 349         """
 350         Return the matrix ``L`` passed to the constructor.
 351
 352         Returns
 353         -------
 354
 355         matrix
 356             The matrix that defines this game's :meth:`payoff` operator.
 357
 358         Examples
 359         --------
 360
 361             >>> from dunshire import *
 362             >>> K = NonnegativeOrthant(3)
 363             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 364             >>> e1 = [1,1,1]
 365             >>> e2 = [1,2,3]
 366             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 367             >>> print(SLG.L())
 368             [  1  -5 -15]
 369             [ -1   2  -3]
 370             [-12 -15   1]
 371             <BLANKLINE>
 372
 373         """
 374         return self._L
 375
 376
 377     def K(self):
 378         """
 379         Return the cone over which this game is played.
 380
 381         Returns
 382         -------
 383
 384         SymmetricCone
 385             The :class:`SymmetricCone` over which this game is played.
 386
 387         Examples
 388         --------
 389
 390             >>> from dunshire import *
 391             >>> K = NonnegativeOrthant(3)
 392             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 393             >>> e1 = [1,1,1]
 394             >>> e2 = [1,2,3]
 395             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 396             >>> print(SLG.K())
 397             Nonnegative orthant in the real 3-space
 398
 399         """
 400         return self._K
 401
 402
 403     def e1(self):
 404         """
 405         Return player one's interior point.
 406
 407         Returns
 408         -------
 409
 410         matrix
 411             The point interior to :meth:`K` affiliated with player one.
 412
 413         Examples
 414         --------
 415
 416             >>> from dunshire import *
 417             >>> K = NonnegativeOrthant(3)
 418             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 419             >>> e1 = [1,1,1]
 420             >>> e2 = [1,2,3]
 421             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 422             >>> print(SLG.e1())
 423             [ 1]
 424             [ 1]
 425             [ 1]
 426             <BLANKLINE>
 427
 428         """
 429         return self._e1
 430
 431
 432     def e2(self):
 433         """
 434         Return player two's interior point.
 435
 436         Returns
 437         -------
 438
 439         matrix
 440             The point interior to :meth:`K` affiliated with player one.
 441
 442         Examples
 443         --------
 444
 445             >>> from dunshire import *
 446             >>> K = NonnegativeOrthant(3)
 447             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 448             >>> e1 = [1,1,1]
 449             >>> e2 = [1,2,3]
 450             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 451             >>> print(SLG.e2())
 452             [ 1]
 453             [ 2]
 454             [ 3]
 455             <BLANKLINE>
 456
 457         """
 458         return self._e2
 459
 460
 461     def payoff(self, strategy1, strategy2):
 462         r"""
 463         Return the payoff associated with ``strategy1`` and ``strategy2``.
 464
 465         The payoff operator takes pairs of strategies to a real
 466         number. For example, if player one's strategy is :math:`x` and
 467         player two's strategy is :math:`y`, then the associated payoff
 468         is :math:`\left\langle L\left(x\right),y \right\rangle` \in
 469         \mathbb{R}. Here, :math:`L` denotes the same linear operator as
 470         :meth:`L`. This method computes the payoff given the two
 471         players' strategies.
 472
 473         Parameters
 474         ----------
 475
 476         strategy1 : matrix
 477             Player one's strategy.
 478
 479         strategy2 : matrix
 480             Player two's strategy.
 481
 482         Returns
 483         -------
 484
 485         float
 486             The payoff for the game when player one plays ``strategy1``
 487             and player two plays ``strategy2``.
 488
 489         Examples
 490         --------
 491
 492         The value of the game should be the payoff at the optimal
 493         strategies::
 494
 495             >>> from dunshire import *
 496             >>> from dunshire.options import ABS_TOL
 497             >>> K = NonnegativeOrthant(3)
 498             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 499             >>> e1 = [1,1,1]
 500             >>> e2 = [1,1,1]
 501             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 502             >>> soln = SLG.solution()
 503             >>> x_bar = soln.player1_optimal()
 504             >>> y_bar = soln.player2_optimal()
 505             >>> abs(SLG.payoff(x_bar, y_bar) - soln.game_value()) < ABS_TOL
 506             True
 507
 508         """
 509         return inner_product(self.L()*strategy1, strategy2)
 510
 511
 512     def dimension(self):
 513         """
 514         Return the dimension of this game.
 515
 516         The dimension of a game is not needed for the theory, but it is
 517         useful for the implementation. We define the dimension of a game
 518         to be the dimension of its underlying cone. Or what is the same,
 519         the dimension of the space from which the strategies are chosen.
 520
 521         Returns
 522         -------
 523
 524         int
 525             The dimension of the cone :meth:`K`, or of the space where
 526             this game is played.
 527
 528         Examples
 529         --------
 530
 531         The dimension of a game over the nonnegative quadrant in the
 532         plane should be two (the dimension of the plane)::
 533
 534             >>> from dunshire import *
 535             >>> K = NonnegativeOrthant(2)
 536             >>> L = [[1,-5],[-1,2]]
 537             >>> e1 = [1,1]
 538             >>> e2 = [1,4]
 539             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 540             >>> SLG.dimension()
 541             2
 542
 543         """
 544         return self.K().dimension()
 545
 546
 547     def _zero(self):
 548         """
 549         Return a column of zeros that fits ``K``.
 550
 551         This is used in our CVXOPT construction.
 552
 553         .. warning::
 554
 555             It is not safe to cache any of the matrices passed to
 556             CVXOPT, because it can clobber them.
 557
 558         Returns
 559         -------
 560
 561         matrix
 562             A ``self.dimension()``-by-``1`` column vector of zeros.
 563
 564         Examples
 565         --------
 566
 567             >>> from dunshire import *
 568             >>> K = NonnegativeOrthant(3)
 569             >>> L = identity(3)
 570             >>> e1 = [1,1,1]
 571             >>> e2 = e1
 572             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 573             >>> print(SLG._zero())
 574             [0.0000000]
 575             [0.0000000]
 576             [0.0000000]
 577             <BLANKLINE>
 578
 579         """
 580         return matrix(0, (self.dimension(), 1), tc='d')
 581
 582
 583     def A(self):
 584         """
 585         Return the matrix ``A`` used in our CVXOPT construction.
 586
 587         This matrix ``A``  appears on the right-hand side of ``Ax = b``
 588         in the statement of the CVXOPT conelp program.
 589
 590         .. warning::
 591
 592             It is not safe to cache any of the matrices passed to
 593             CVXOPT, because it can clobber them.
 594
 595         Returns
 596         -------
 597
 598         matrix
 599             A ``1``-by-``(1 + self.dimension())`` row vector. Its first
 600             entry is zero, and the rest are the entries of ``e2``.
 601
 602         Examples
 603         --------
 604
 605             >>> from dunshire import *
 606             >>> K = NonnegativeOrthant(3)
 607             >>> L = [[1,1,1],[1,1,1],[1,1,1]]
 608             >>> e1 = [1,1,1]
 609             >>> e2 = [1,2,3]
 610             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 611             >>> print(SLG.A())
 612             [0.0000000 1.0000000 2.0000000 3.0000000]
 613             <BLANKLINE>
 614
 615         """
 616         return matrix([0, self.e2()], (1, self.dimension() + 1), 'd')
 617
 618
 619
 620     def _G(self):
 621         r"""
 622         Return the matrix ``G`` used in our CVXOPT construction.
 623
 624         Thus matrix ``G`` appears on the left-hand side of ``Gx + s = h``
 625         in the statement of the CVXOPT conelp program.
 626
 627         .. warning::
 628
 629             It is not safe to cache any of the matrices passed to
 630             CVXOPT, because it can clobber them.
 631
 632         Returns
 633         -------
 634
 635         matrix
 636             A ``2*self.dimension()``-by-``(1 + self.dimension())`` matrix.
 637
 638         Examples
 639         --------
 640
 641             >>> from dunshire import *
 642             >>> K = NonnegativeOrthant(3)
 643             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 644             >>> e1 = [1,2,3]
 645             >>> e2 = [1,1,1]
 646             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 647             >>> print(SLG._G())
 648             [  0.0000000  -1.0000000   0.0000000   0.0000000]
 649             [  0.0000000   0.0000000  -1.0000000   0.0000000]
 650             [  0.0000000   0.0000000   0.0000000  -1.0000000]
 651             [  1.0000000  -4.0000000  -5.0000000  -6.0000000]
 652             [  2.0000000  -7.0000000  -8.0000000  -9.0000000]
 653             [  3.0000000 -10.0000000 -11.0000000 -12.0000000]
 654             <BLANKLINE>
 655
 656         """
 657         identity_matrix = identity(self.dimension())
 658         return append_row(append_col(self._zero(), -identity_matrix),
 659                           append_col(self.e1(), -self.L()))
 660
 661
 662     def _c(self):
 663         """
 664         Return the vector ``c`` used in our CVXOPT construction.
 665
 666         The column vector ``c``  appears in the objective function
 667         value ``<c,x>`` in the statement of the CVXOPT conelp program.
 668
 669         .. warning::
 670
 671             It is not safe to cache any of the matrices passed to
 672             CVXOPT, because it can clobber them.
 673
 674         Returns
 675         -------
 676
 677         matrix
 678             A ``self.dimension()``-by-``1`` column vector.
 679
 680         Examples
 681         --------
 682
 683             >>> from dunshire import *
 684             >>> K = NonnegativeOrthant(3)
 685             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 686             >>> e1 = [1,2,3]
 687             >>> e2 = [1,1,1]
 688             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 689             >>> print(SLG._c())
 690             [-1.0000000]
 691             [ 0.0000000]
 692             [ 0.0000000]
 693             [ 0.0000000]
 694             <BLANKLINE>
 695
 696         """
 697         return matrix([-1, self._zero()])
 698
 699
 700     def C(self):
 701         """
 702         Return the cone ``C`` used in our CVXOPT construction.
 703
 704         The cone ``C`` is the cone over which the conelp program takes
 705         place.
 706
 707         Returns
 708         -------
 709
 710         CartesianProduct
 711             The cartesian product of ``K`` with itself.
 712
 713         Examples
 714         --------
 715
 716             >>> from dunshire import *
 717             >>> K = NonnegativeOrthant(3)
 718             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 719             >>> e1 = [1,2,3]
 720             >>> e2 = [1,1,1]
 721             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 722             >>> print(SLG.C())
 723             Cartesian product of dimension 6 with 2 factors:
 724               * Nonnegative orthant in the real 3-space
 725               * Nonnegative orthant in the real 3-space
 726
 727         """
 728         return CartesianProduct(self._K, self._K)
 729
 730     def _h(self):
 731         """
 732         Return the ``h`` vector used in our CVXOPT construction.
 733
 734         The ``h`` vector appears on the right-hand side of :math:`Gx + s
 735         = h` in the statement of the CVXOPT conelp program.
 736
 737         .. warning::
 738
 739             It is not safe to cache any of the matrices passed to
 740             CVXOPT, because it can clobber them.
 741
 742         Returns
 743         -------
 744
 745         matrix
 746             A ``2*self.dimension()``-by-``1`` column vector of zeros.
 747
 748         Examples
 749         --------
 750
 751             >>> from dunshire import *
 752             >>> K = NonnegativeOrthant(3)
 753             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 754             >>> e1 = [1,2,3]
 755             >>> e2 = [1,1,1]
 756             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 757             >>> print(SLG._h())
 758             [0.0000000]
 759             [0.0000000]
 760             [0.0000000]
 761             [0.0000000]
 762             [0.0000000]
 763             [0.0000000]
 764             <BLANKLINE>
 765
 766         """
 767
 768         return matrix([self._zero(), self._zero()])
 769
 770
 771     @staticmethod
 772     def b():
 773         """
 774         Return the ``b`` vector used in our CVXOPT construction.
 775
 776         The vector ``b`` appears on the right-hand side of :math:`Ax =
 777         b` in the statement of the CVXOPT conelp program.
 778
 779         This method is static because the dimensions and entries of
 780         ``b`` are known beforehand, and don't depend on any other
 781         properties of the game.
 782
 783         .. warning::
 784
 785             It is not safe to cache any of the matrices passed to
 786             CVXOPT, because it can clobber them.
 787
 788         Returns
 789         -------
 790
 791         matrix
 792             A ``1``-by-``1`` matrix containing a single entry ``1``.
 793
 794         Examples
 795         --------
 796
 797             >>> from dunshire import *
 798             >>> K = NonnegativeOrthant(3)
 799             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 800             >>> e1 = [1,2,3]
 801             >>> e2 = [1,1,1]
 802             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 803             >>> print(SLG.b())
 804             [1.0000000]
 805             <BLANKLINE>
 806
 807         """
 808         return matrix([1], tc='d')
 809
 810
 811     def player1_start(self):
 812         """
 813         Return a feasible starting point for player one.
 814
 815         This starting point is for the CVXOPT formulation and not for
 816         the original game. The basic premise is that if you normalize
 817         :meth:`e2`, then you get a point in :meth:`K` that makes a unit
 818         inner product with :meth:`e2`. We then get to choose the primal
 819         objective function value such that the constraint involving
 820         :meth:`L` is satisfied.
 821         """
 822         p = self.e2() / (norm(self.e2()) ** 2)
 823
 824         # Compute the distance from p to the outside of K.
 825         if isinstance(self.K(), NonnegativeOrthant):
 826             # How far is it to a wall?
 827             dist = min(list(self.e1()))
 828         elif isinstance(self.K(), IceCream):
 829             # How far is it to the boundary of the ball that defines
 830             # the ice-cream cone at a given height? Now draw a
 831             # 45-45-90 triangle and the shortest distance to the
 832             # outside of the cone should be 1/sqrt(2) of that.
 833             # It works in R^2, so it works everywhere, right?
 834             # We use "2" because it's better numerically than sqrt(2).
 835             height = self.e1()[0]
 836             radius = norm(self.e1()[1:])
 837             dist = (height - radius) / 2
 838         else:
 839             raise NotImplementedError
 840
 841         nu = - specnorm(self.L())/(dist*norm(self.e2()))
 842         x = matrix([nu,p], (self.dimension() + 1, 1))
 843         s = - self._G()*x
 844
 845         return {'x': x, 's': s}
 846
 847
 848     def player2_start(self):
 849         """
 850         Return a feasible starting point for player two.
 851         """
 852         q = self.e1() / (norm(self.e1()) ** 2)
 853
 854         # Compute the distance from p to the outside of K.
 855         if isinstance(self.K(), NonnegativeOrthant):
 856             # How far is it to a wall?
 857             dist = min(list(self.e2()))
 858         elif isinstance(self.K(), IceCream):
 859             # How far is it to the boundary of the ball that defines
 860             # the ice-cream cone at a given height? Now draw a
 861             # 45-45-90 triangle and the shortest distance to the
 862             # outside of the cone should be 1/sqrt(2) of that.
 863             # It works in R^2, so it works everywhere, right?
 864             # We use "2" because it's better numerically than sqrt(2).
 865             height = self.e2()[0]
 866             radius = norm(self.e2()[1:])
 867             dist = (height - radius) / 2
 868         else:
 869             raise NotImplementedError
 870
 871         omega = specnorm(self.L())/(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 solution(self):
 881         """
 882         Solve this linear game and return a :class:`Solution`.
 883
 884         Returns
 885         -------
 886
 887         :class:`Solution`
 888             A :class:`Solution` object describing the game's value and
 889             the optimal strategies of both players.
 890
 891         Raises
 892         ------
 893         GameUnsolvableException
 894             If the game could not be solved (if an optimal solution to its
 895             associated cone program was not found).
 896
 897         PoorScalingException
 898             If the game could not be solved because CVXOPT crashed while
 899             trying to take the square root of a negative number.
 900
 901         Examples
 902         --------
 903
 904         This example is computed in Gowda and Ravindran in the section
 905         "The value of a Z-transformation"::
 906
 907             >>> from dunshire import *
 908             >>> K = NonnegativeOrthant(3)
 909             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 910             >>> e1 = [1,1,1]
 911             >>> e2 = [1,1,1]
 912             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 913             >>> print(SLG.solution())
 914             Game value: -6.172...
 915             Player 1 optimal:
 916               [0.551...]
 917               [0.000...]
 918               [0.448...]
 919             Player 2 optimal:
 920               [0.448...]
 921               [0.000...]
 922               [0.551...]
 923
 924         The value of the following game can be computed using the fact
 925         that the identity is invertible::
 926
 927             >>> from dunshire import *
 928             >>> K = NonnegativeOrthant(3)
 929             >>> L = [[1,0,0],[0,1,0],[0,0,1]]
 930             >>> e1 = [1,2,3]
 931             >>> e2 = [4,5,6]
 932             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 933             >>> print(SLG.solution())
 934             Game value: 0.031...
 935             Player 1 optimal:
 936               [0.031...]
 937               [0.062...]
 938               [0.093...]
 939             Player 2 optimal:
 940               [0.125...]
 941               [0.156...]
 942               [0.187...]
 943
 944         This is another Gowda/Ravindran example that is supposed to have
 945         a negative game value::
 946
 947             >>> from dunshire import *
 948             >>> from dunshire.options import ABS_TOL
 949             >>> L = [[1, -2], [-2, 1]]
 950             >>> K = NonnegativeOrthant(2)
 951             >>> e1 = [1, 1]
 952             >>> e2 = e1
 953             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 954             >>> SLG.solution().game_value() < -ABS_TOL
 955             True
 956
 957         The following two games are problematic numerically, but we
 958         should be able to solve them::
 959
 960             >>> from dunshire import *
 961             >>> L = [[-0.95237953890954685221, 1.83474556206462535712],
 962             ...      [ 1.30481749924621448500, 1.65278664543326403447]]
 963             >>> K = NonnegativeOrthant(2)
 964             >>> e1 = [0.95477167524644313001, 0.63270781756540095397]
 965             >>> e2 = [0.39633793037154141370, 0.10239281495640320530]
 966             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 967             >>> print(SLG.solution())
 968             Game value: 18.767...
 969             Player 1 optimal:
 970               [0.000...]
 971               [9.766...]
 972             Player 2 optimal:
 973               [1.047...]
 974               [0.000...]
 975
 976         ::
 977
 978             >>> from dunshire import *
 979             >>> L = [[1.54159395026049472754, 2.21344728574316684799],
 980             ...      [1.33147433507846657541, 1.17913616272988108769]]
 981             >>> K = NonnegativeOrthant(2)
 982             >>> e1 = [0.39903040089404784307, 0.12377403622479113410]
 983             >>> e2 = [0.15695181142215544612, 0.85527381344651265405]
 984             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 985             >>> print(SLG.solution())
 986             Game value: 24.614...
 987             Player 1 optimal:
 988               [6.371...]
 989               [0.000...]
 990             Player 2 optimal:
 991               [2.506...]
 992               [0.000...]
 993
 994         """
 995         try:
 996             opts = {'show_progress': False}
 997             soln_dict = solvers.conelp(self._c(),
 998                                        self._G(),
 999                                        self._h(),
1000                                        self.C().cvxopt_dims(),
1001                                        self.A(),
1002                                        self.b(),
1003                                        primalstart=self.player1_start(),
1004                                        options=opts)
1005         except ValueError as error:
1006             if str(error) == 'math domain error':
1007                 # Oops, CVXOPT tried to take the square root of a
1008                 # negative number. Report some details about the game
1009                 # rather than just the underlying CVXOPT crash.
1010                 printing.options['dformat'] = options.DEBUG_FLOAT_FORMAT
1011                 raise PoorScalingException(self)
1012             else:
1013                 raise error
1014
1015         # The optimal strategies are named ``p`` and ``q`` in the
1016         # background documentation, and we need to extract them from
1017         # the CVXOPT ``x`` and ``z`` variables. The objective values
1018         # :math:`nu` and :math:`omega` can also be found in the CVXOPT
1019         # ``x`` and ``y`` variables; however, they're stored
1020         # conveniently as separate entries in the solution dictionary.
1021         p1_value = -soln_dict['primal objective']
1022         p2_value = -soln_dict['dual objective']
1023         p1_optimal = soln_dict['x'][1:]
1024         p2_optimal = soln_dict['z'][self.dimension():]
1025
1026         # The "status" field contains "optimal" if everything went
1027         # according to plan. Other possible values are "primal
1028         # infeasible", "dual infeasible", "unknown", all of which mean
1029         # we didn't get a solution.
1030         #
1031         # The "infeasible" ones are the worst, since they indicate
1032         # that CVXOPT is convinced the problem is infeasible (and that
1033         # cannot happen).
1034         if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
1035             printing.options['dformat'] = options.DEBUG_FLOAT_FORMAT
1036             raise GameUnsolvableException(self, soln_dict)
1037
1038         # The "optimal" and "unknown" results, we actually treat the
1039         # same. Even if CVXOPT bails out due to numerical difficulty,
1040         # it will have some candidate points in mind. If those
1041         # candidates are good enough, we take them. We do the same
1042         # check (perhaps pointlessly so) for "optimal" results.
1043         #
1044         # First we check that the primal/dual objective values are
1045         # close enough (one could be low by ABS_TOL, the other high by
1046         # it) because otherwise CVXOPT might return "unknown" and give
1047         # us two points in the cone that are nowhere near optimal.
1048         if abs(p1_value - p2_value) > 2*options.ABS_TOL:
1049             printing.options['dformat'] = options.DEBUG_FLOAT_FORMAT
1050             raise GameUnsolvableException(self, soln_dict)
1051
1052         # And we also check that the points it gave us belong to the
1053         # cone, just in case...
1054         if (p1_optimal not in self._K) or (p2_optimal not in self._K):
1055             printing.options['dformat'] = options.DEBUG_FLOAT_FORMAT
1056             raise GameUnsolvableException(self, soln_dict)
1057
1058         # For the game value, we could use any of:
1059         #
1060         #   * p1_value
1061         #   * p2_value
1062         #   * (p1_value + p2_value)/2
1063         #   * the game payoff
1064         #
1065         # We want the game value to be the payoff, however, so it
1066         # makes the most sense to just use that, even if it means we
1067         # can't test the fact that p1_value/p2_value are close to the
1068         # payoff.
1069         payoff = self.payoff(p1_optimal, p2_optimal)
1070         return Solution(payoff, p1_optimal, p2_optimal)
1071
1072
1073     def condition(self):
1074         r"""
1075         Return the condition number of this game.
1076
1077         In the CVXOPT construction of this game, two matrices ``G`` and
1078         ``A`` appear. When those matrices are nasty, numerical problems
1079         can show up. We define the condition number of this game to be
1080         the average of the condition numbers of ``G`` and ``A`` in the
1081         CVXOPT construction. If the condition number of this game is
1082         high, then you can expect numerical difficulty (such as
1083         :class:`PoorScalingException`).
1084
1085         Returns
1086         -------
1087
1088         float
1089             A real number greater than or equal to one that measures how
1090             bad this game is numerically.
1091
1092         Examples
1093         --------
1094
1095         >>> from dunshire import *
1096         >>> K = NonnegativeOrthant(1)
1097         >>> L = [[1]]
1098         >>> e1 = [1]
1099         >>> e2 = e1
1100         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1101         >>> actual = SLG.condition()
1102         >>> expected = 1.8090169943749477
1103         >>> abs(actual - expected) < options.ABS_TOL
1104         True
1105
1106         """
1107         return (condition_number(self._G()) + condition_number(self.A()))/2
1108
1109
1110     def dual(self):
1111         r"""
1112         Return the dual game to this game.
1113
1114         If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
1115         then its dual is :math:`G^{*} =
1116         \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
1117         is symmetric, :math:`K^{*} = K`.
1118
1119         Examples
1120         --------
1121
1122             >>> from dunshire import *
1123             >>> K = NonnegativeOrthant(3)
1124             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
1125             >>> e1 = [1,1,1]
1126             >>> e2 = [1,2,3]
1127             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1128             >>> print(SLG.dual())
1129             The linear game (L, K, e1, e2) where
1130               L = [  1  -1 -12]
1131                   [ -5   2 -15]
1132                   [-15  -3   1],
1133               K = Nonnegative orthant in the real 3-space,
1134               e1 = [ 1]
1135                    [ 2]
1136                    [ 3],
1137               e2 = [ 1]
1138                    [ 1]
1139                    [ 1],
1140               Condition((L, K, e1, e2)) = 44.476...
1141
1142         """
1143         # We pass ``self.L()`` right back into the constructor, because
1144         # it will be transposed there. And keep in mind that ``self._K``
1145         # is its own dual.
1146         return SymmetricLinearGame(self.L(),
1147                                    self.K(),
1148                                    self.e2(),
1149                                    self.e1())