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 math import sqrt
   8
   9 from cvxopt import matrix, printing, solvers
  10 from .cones import CartesianProduct, IceCream, NonnegativeOrthant
  11 from .errors import GameUnsolvableException, PoorScalingException
  12 from .matrices import (append_col, append_row, condition_number, identity,
  13                        inner_product, norm, specnorm)
  14 from . import options
  15
  16 printing.options['dformat'] = options.FLOAT_FORMAT
  17
  18 class Solution:
  19     """
  20     A representation of the solution of a linear game. It should contain
  21     the value of the game, and both players' strategies.
  22
  23     Examples
  24     --------
  25
  26         >>> print(Solution(10, matrix([1,2]), matrix([3,4])))
  27         Game value: 10.000...
  28         Player 1 optimal:
  29           [ 1]
  30           [ 2]
  31         Player 2 optimal:
  32           [ 3]
  33           [ 4]
  34
  35     """
  36     def __init__(self, game_value, p1_optimal, p2_optimal):
  37         """
  38         Create a new Solution object from a game value and two optimal
  39         strategies for the players.
  40         """
  41         self._game_value = game_value
  42         self._player1_optimal = p1_optimal
  43         self._player2_optimal = p2_optimal
  44
  45     def __str__(self):
  46         """
  47         Return a string describing the solution of a linear game.
  48
  49         The three data that are described are,
  50
  51           * The value of the game.
  52           * The optimal strategy of player one.
  53           * The optimal strategy of player two.
  54
  55         The two optimal strategy vectors are indented by two spaces.
  56         """
  57         tpl = 'Game value: {:.7f}\n' \
  58               'Player 1 optimal:{:s}\n' \
  59               'Player 2 optimal:{:s}'
  60
  61         p1_str = '\n{!s}'.format(self.player1_optimal())
  62         p1_str = '\n  '.join(p1_str.splitlines())
  63         p2_str = '\n{!s}'.format(self.player2_optimal())
  64         p2_str = '\n  '.join(p2_str.splitlines())
  65
  66         return tpl.format(self.game_value(), p1_str, p2_str)
  67
  68
  69     def game_value(self):
  70         """
  71         Return the game value for this solution.
  72
  73         Examples
  74         --------
  75
  76             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  77             >>> s.game_value()
  78             10
  79
  80         """
  81         return self._game_value
  82
  83
  84     def player1_optimal(self):
  85         """
  86         Return player one's optimal strategy in this solution.
  87
  88         Examples
  89         --------
  90
  91             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  92             >>> print(s.player1_optimal())
  93             [ 1]
  94             [ 2]
  95             <BLANKLINE>
  96
  97         """
  98         return self._player1_optimal
  99
 100
 101     def player2_optimal(self):
 102         """
 103         Return player two's optimal strategy in this solution.
 104
 105         Examples
 106         --------
 107
 108             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
 109             >>> print(s.player2_optimal())
 110             [ 3]
 111             [ 4]
 112             <BLANKLINE>
 113
 114         """
 115         return self._player2_optimal
 116
 117
 118 class SymmetricLinearGame:
 119     r"""
 120     A representation of a symmetric linear game.
 121
 122     The data for a symmetric linear game are,
 123
 124       * A "payoff" operator ``L``.
 125       * A symmetric cone ``K``.
 126       * Two points ``e1`` and ``e2`` in the interior of ``K``.
 127
 128     The ambient space is assumed to be the span of ``K``.
 129
 130     With those data understood, the game is played as follows. Players
 131     one and two choose points :math:`x` and :math:`y` respectively, from
 132     their respective strategy sets,
 133
 134     .. math::
 135         \begin{aligned}
 136             \Delta_{1}
 137             &=
 138             \left\{
 139               x \in K \ \middle|\ \left\langle x, e_{2} \right\rangle = 1
 140             \right\}\\
 141             \Delta_{2}
 142             &=
 143             \left\{
 144               y \in K \ \middle|\ \left\langle y, e_{1} \right\rangle = 1
 145             \right\}.
 146         \end{aligned}
 147
 148     Afterwards, a "payout" is computed as :math:`\left\langle
 149     L\left(x\right), y \right\rangle` and is paid to player one out of
 150     player two's pocket. The game is therefore zero sum, and we suppose
 151     that player one would like to guarantee himself the largest minimum
 152     payout possible. That is, player one wishes to,
 153
 154     .. math::
 155         \begin{aligned}
 156             \text{maximize }
 157               &\underset{y \in \Delta_{2}}{\min}\left(
 158                 \left\langle L\left(x\right), y \right\rangle
 159               \right)\\
 160             \text{subject to } & x \in \Delta_{1}.
 161         \end{aligned}
 162
 163     Player two has the simultaneous goal to,
 164
 165     .. math::
 166         \begin{aligned}
 167             \text{minimize }
 168               &\underset{x \in \Delta_{1}}{\max}\left(
 169                 \left\langle L\left(x\right), y \right\rangle
 170               \right)\\
 171             \text{subject to } & y \in \Delta_{2}.
 172         \end{aligned}
 173
 174     These goals obviously conflict (the game is zero sum), but an
 175     existence theorem guarantees at least one optimal min-max solution
 176     from which neither player would like to deviate. This class is
 177     able to find such a solution.
 178
 179     Parameters
 180     ----------
 181
 182     L : list of list of float
 183         A matrix represented as a list of ROWS. This representation
 184         agrees with (for example) SageMath and NumPy, but not with CVXOPT
 185         (whose matrix constructor accepts a list of columns).
 186
 187     K : :class:`SymmetricCone`
 188         The symmetric cone instance over which the game is played.
 189
 190     e1 : iterable float
 191         The interior point of ``K`` belonging to player one; it
 192         can be of any iterable type having the correct length.
 193
 194     e2 : iterable float
 195         The interior point of ``K`` belonging to player two; it
 196         can be of any enumerable type having the correct length.
 197
 198     Raises
 199     ------
 200
 201     ValueError
 202         If either ``e1`` or ``e2`` lie outside of the cone ``K``.
 203
 204     Examples
 205     --------
 206
 207         >>> from dunshire import *
 208         >>> K = NonnegativeOrthant(3)
 209         >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 210         >>> e1 = [1,1,1]
 211         >>> e2 = [1,2,3]
 212         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 213         >>> print(SLG)
 214         The linear game (L, K, e1, e2) where
 215           L = [  1  -5 -15]
 216               [ -1   2  -3]
 217               [-12 -15   1],
 218           K = Nonnegative orthant in the real 3-space,
 219           e1 = [ 1]
 220                [ 1]
 221                [ 1],
 222           e2 = [ 1]
 223                [ 2]
 224                [ 3],
 225           Condition((L, K, e1, e2)) = 31.834...
 226
 227     Lists can (and probably should) be used for every argument::
 228
 229         >>> from dunshire import *
 230         >>> K = NonnegativeOrthant(2)
 231         >>> L = [[1,0],[0,1]]
 232         >>> e1 = [1,1]
 233         >>> e2 = [1,1]
 234         >>> G = SymmetricLinearGame(L, K, e1, e2)
 235         >>> print(G)
 236         The linear game (L, K, e1, e2) where
 237           L = [ 1  0]
 238               [ 0  1],
 239           K = Nonnegative orthant in the real 2-space,
 240           e1 = [ 1]
 241                [ 1],
 242           e2 = [ 1]
 243                [ 1],
 244           Condition((L, K, e1, e2)) = 1.707...
 245
 246     The points ``e1`` and ``e2`` can also be passed as some other
 247     enumerable type (of the correct length) without much harm, since
 248     there is no row/column ambiguity::
 249
 250         >>> import cvxopt
 251         >>> import numpy
 252         >>> from dunshire import *
 253         >>> K = NonnegativeOrthant(2)
 254         >>> L = [[1,0],[0,1]]
 255         >>> e1 = cvxopt.matrix([1,1])
 256         >>> e2 = numpy.matrix([1,1])
 257         >>> G = SymmetricLinearGame(L, K, e1, e2)
 258         >>> print(G)
 259         The linear game (L, K, e1, e2) where
 260           L = [ 1  0]
 261               [ 0  1],
 262           K = Nonnegative orthant in the real 2-space,
 263           e1 = [ 1]
 264                [ 1],
 265           e2 = [ 1]
 266                [ 1],
 267           Condition((L, K, e1, e2)) = 1.707...
 268
 269     However, ``L`` will always be intepreted as a list of rows, even
 270     if it is passed as a :class:`cvxopt.base.matrix` which is
 271     otherwise indexed by columns::
 272
 273         >>> import cvxopt
 274         >>> from dunshire import *
 275         >>> K = NonnegativeOrthant(2)
 276         >>> L = [[1,2],[3,4]]
 277         >>> e1 = [1,1]
 278         >>> e2 = e1
 279         >>> G = SymmetricLinearGame(L, K, e1, e2)
 280         >>> print(G)
 281         The linear game (L, K, e1, e2) where
 282           L = [ 1  2]
 283               [ 3  4],
 284           K = Nonnegative orthant in the real 2-space,
 285           e1 = [ 1]
 286                [ 1],
 287           e2 = [ 1]
 288                [ 1],
 289           Condition((L, K, e1, e2)) = 6.073...
 290         >>> L = cvxopt.matrix(L)
 291         >>> print(L)
 292         [ 1  3]
 293         [ 2  4]
 294         <BLANKLINE>
 295         >>> G = SymmetricLinearGame(L, K, e1, e2)
 296         >>> print(G)
 297         The linear game (L, K, e1, e2) where
 298           L = [ 1  2]
 299               [ 3  4],
 300           K = Nonnegative orthant in the real 2-space,
 301           e1 = [ 1]
 302                [ 1],
 303           e2 = [ 1]
 304                [ 1],
 305           Condition((L, K, e1, e2)) = 6.073...
 306
 307     """
 308     def __init__(self, L, K, e1, e2):
 309         """
 310         Create a new SymmetricLinearGame object.
 311         """
 312         self._K = K
 313         self._e1 = matrix(e1, (K.dimension(), 1))
 314         self._e2 = matrix(e2, (K.dimension(), 1))
 315
 316         # Our input ``L`` is indexed by rows but CVXOPT matrices are
 317         # indexed by columns, so we need to transpose the input before
 318         # feeding it to CVXOPT.
 319         self._L = matrix(L, (K.dimension(), K.dimension())).trans()
 320
 321         if not self._e1 in K:
 322             raise ValueError('the point e1 must lie in the interior of K')
 323
 324         if not self._e2 in K:
 325             raise ValueError('the point e2 must lie in the interior of K')
 326
 327
 328
 329     def __str__(self):
 330         """
 331         Return a string representation of this game.
 332         """
 333         tpl = 'The linear game (L, K, e1, e2) where\n' \
 334               '  L = {:s},\n' \
 335               '  K = {!s},\n' \
 336               '  e1 = {:s},\n' \
 337               '  e2 = {:s},\n' \
 338               '  Condition((L, K, e1, e2)) = {:f}.'
 339         indented_L = '\n      '.join(str(self.L()).splitlines())
 340         indented_e1 = '\n       '.join(str(self.e1()).splitlines())
 341         indented_e2 = '\n       '.join(str(self.e2()).splitlines())
 342
 343         return tpl.format(indented_L,
 344                           str(self.K()),
 345                           indented_e1,
 346                           indented_e2,
 347                           self.condition())
 348
 349
 350     def L(self):
 351         """
 352         Return the matrix ``L`` passed to the constructor.
 353
 354         Returns
 355         -------
 356
 357         matrix
 358             The matrix that defines this game's :meth:`payoff` operator.
 359
 360         Examples
 361         --------
 362
 363             >>> from dunshire import *
 364             >>> K = NonnegativeOrthant(3)
 365             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 366             >>> e1 = [1,1,1]
 367             >>> e2 = [1,2,3]
 368             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 369             >>> print(SLG.L())
 370             [  1  -5 -15]
 371             [ -1   2  -3]
 372             [-12 -15   1]
 373             <BLANKLINE>
 374
 375         """
 376         return self._L
 377
 378
 379     def K(self):
 380         """
 381         Return the cone over which this game is played.
 382
 383         Returns
 384         -------
 385
 386         SymmetricCone
 387             The :class:`SymmetricCone` over which this game is played.
 388
 389         Examples
 390         --------
 391
 392             >>> from dunshire import *
 393             >>> K = NonnegativeOrthant(3)
 394             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 395             >>> e1 = [1,1,1]
 396             >>> e2 = [1,2,3]
 397             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 398             >>> print(SLG.K())
 399             Nonnegative orthant in the real 3-space
 400
 401         """
 402         return self._K
 403
 404
 405     def e1(self):
 406         """
 407         Return player one's interior point.
 408
 409         Returns
 410         -------
 411
 412         matrix
 413             The point interior to :meth:`K` affiliated with player one.
 414
 415         Examples
 416         --------
 417
 418             >>> from dunshire import *
 419             >>> K = NonnegativeOrthant(3)
 420             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 421             >>> e1 = [1,1,1]
 422             >>> e2 = [1,2,3]
 423             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 424             >>> print(SLG.e1())
 425             [ 1]
 426             [ 1]
 427             [ 1]
 428             <BLANKLINE>
 429
 430         """
 431         return self._e1
 432
 433
 434     def e2(self):
 435         """
 436         Return player two's interior point.
 437
 438         Returns
 439         -------
 440
 441         matrix
 442             The point interior to :meth:`K` affiliated with player one.
 443
 444         Examples
 445         --------
 446
 447             >>> from dunshire import *
 448             >>> K = NonnegativeOrthant(3)
 449             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 450             >>> e1 = [1,1,1]
 451             >>> e2 = [1,2,3]
 452             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 453             >>> print(SLG.e2())
 454             [ 1]
 455             [ 2]
 456             [ 3]
 457             <BLANKLINE>
 458
 459         """
 460         return self._e2
 461
 462
 463     def payoff(self, strategy1, strategy2):
 464         r"""
 465         Return the payoff associated with ``strategy1`` and ``strategy2``.
 466
 467         The payoff operator takes pairs of strategies to a real
 468         number. For example, if player one's strategy is :math:`x` and
 469         player two's strategy is :math:`y`, then the associated payoff
 470         is :math:`\left\langle L\left(x\right),y \right\rangle` \in
 471         \mathbb{R}. Here, :math:`L` denotes the same linear operator as
 472         :meth:`L`. This method computes the payoff given the two
 473         players' strategies.
 474
 475         Parameters
 476         ----------
 477
 478         strategy1 : matrix
 479             Player one's strategy.
 480
 481         strategy2 : matrix
 482             Player two's strategy.
 483
 484         Returns
 485         -------
 486
 487         float
 488             The payoff for the game when player one plays ``strategy1``
 489             and player two plays ``strategy2``.
 490
 491         Examples
 492         --------
 493
 494         The value of the game should be the payoff at the optimal
 495         strategies::
 496
 497             >>> from dunshire import *
 498             >>> from dunshire.options import ABS_TOL
 499             >>> K = NonnegativeOrthant(3)
 500             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 501             >>> e1 = [1,1,1]
 502             >>> e2 = [1,1,1]
 503             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 504             >>> soln = SLG.solution()
 505             >>> x_bar = soln.player1_optimal()
 506             >>> y_bar = soln.player2_optimal()
 507             >>> abs(SLG.payoff(x_bar, y_bar) - soln.game_value()) < ABS_TOL
 508             True
 509
 510         """
 511         return inner_product(self.L()*strategy1, strategy2)
 512
 513
 514     def dimension(self):
 515         """
 516         Return the dimension of this game.
 517
 518         The dimension of a game is not needed for the theory, but it is
 519         useful for the implementation. We define the dimension of a game
 520         to be the dimension of its underlying cone. Or what is the same,
 521         the dimension of the space from which the strategies are chosen.
 522
 523         Returns
 524         -------
 525
 526         int
 527             The dimension of the cone :meth:`K`, or of the space where
 528             this game is played.
 529
 530         Examples
 531         --------
 532
 533         The dimension of a game over the nonnegative quadrant in the
 534         plane should be two (the dimension of the plane)::
 535
 536             >>> from dunshire import *
 537             >>> K = NonnegativeOrthant(2)
 538             >>> L = [[1,-5],[-1,2]]
 539             >>> e1 = [1,1]
 540             >>> e2 = [1,4]
 541             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 542             >>> SLG.dimension()
 543             2
 544
 545         """
 546         return self.K().dimension()
 547
 548
 549     def _zero(self):
 550         """
 551         Return a column of zeros that fits ``K``.
 552
 553         This is used in our CVXOPT construction.
 554
 555         .. warning::
 556
 557             It is not safe to cache any of the matrices passed to
 558             CVXOPT, because it can clobber them.
 559
 560         Returns
 561         -------
 562
 563         matrix
 564             A ``self.dimension()``-by-``1`` column vector of zeros.
 565
 566         Examples
 567         --------
 568
 569             >>> from dunshire import *
 570             >>> K = NonnegativeOrthant(3)
 571             >>> L = identity(3)
 572             >>> e1 = [1,1,1]
 573             >>> e2 = e1
 574             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 575             >>> print(SLG._zero())
 576             [0.0000000]
 577             [0.0000000]
 578             [0.0000000]
 579             <BLANKLINE>
 580
 581         """
 582         return matrix(0, (self.dimension(), 1), tc='d')
 583
 584
 585     def A(self):
 586         """
 587         Return the matrix ``A`` used in our CVXOPT construction.
 588
 589         This matrix ``A``  appears on the right-hand side of ``Ax = b``
 590         in the statement of the CVXOPT conelp program.
 591
 592         .. warning::
 593
 594             It is not safe to cache any of the matrices passed to
 595             CVXOPT, because it can clobber them.
 596
 597         Returns
 598         -------
 599
 600         matrix
 601             A ``1``-by-``(1 + self.dimension())`` row vector. Its first
 602             entry is zero, and the rest are the entries of ``e2``.
 603
 604         Examples
 605         --------
 606
 607             >>> from dunshire import *
 608             >>> K = NonnegativeOrthant(3)
 609             >>> L = [[1,1,1],[1,1,1],[1,1,1]]
 610             >>> e1 = [1,1,1]
 611             >>> e2 = [1,2,3]
 612             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 613             >>> print(SLG.A())
 614             [0.0000000 1.0000000 2.0000000 3.0000000]
 615             <BLANKLINE>
 616
 617         """
 618         return matrix([0, self.e2()], (1, self.dimension() + 1), 'd')
 619
 620
 621
 622     def _G(self):
 623         r"""
 624         Return the matrix ``G`` used in our CVXOPT construction.
 625
 626         Thus matrix ``G`` appears on the left-hand side of ``Gx + s = h``
 627         in the statement of the CVXOPT conelp program.
 628
 629         .. warning::
 630
 631             It is not safe to cache any of the matrices passed to
 632             CVXOPT, because it can clobber them.
 633
 634         Returns
 635         -------
 636
 637         matrix
 638             A ``2*self.dimension()``-by-``(1 + self.dimension())`` matrix.
 639
 640         Examples
 641         --------
 642
 643             >>> from dunshire import *
 644             >>> K = NonnegativeOrthant(3)
 645             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 646             >>> e1 = [1,2,3]
 647             >>> e2 = [1,1,1]
 648             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 649             >>> print(SLG._G())
 650             [  0.0000000  -1.0000000   0.0000000   0.0000000]
 651             [  0.0000000   0.0000000  -1.0000000   0.0000000]
 652             [  0.0000000   0.0000000   0.0000000  -1.0000000]
 653             [  1.0000000  -4.0000000  -5.0000000  -6.0000000]
 654             [  2.0000000  -7.0000000  -8.0000000  -9.0000000]
 655             [  3.0000000 -10.0000000 -11.0000000 -12.0000000]
 656             <BLANKLINE>
 657
 658         """
 659         identity_matrix = identity(self.dimension())
 660         return append_row(append_col(self._zero(), -identity_matrix),
 661                           append_col(self.e1(), -self.L()))
 662
 663
 664     def _c(self):
 665         """
 666         Return the vector ``c`` used in our CVXOPT construction.
 667
 668         The column vector ``c``  appears in the objective function
 669         value ``<c,x>`` in the statement of the CVXOPT conelp program.
 670
 671         .. warning::
 672
 673             It is not safe to cache any of the matrices passed to
 674             CVXOPT, because it can clobber them.
 675
 676         Returns
 677         -------
 678
 679         matrix
 680             A ``self.dimension()``-by-``1`` column vector.
 681
 682         Examples
 683         --------
 684
 685             >>> from dunshire import *
 686             >>> K = NonnegativeOrthant(3)
 687             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 688             >>> e1 = [1,2,3]
 689             >>> e2 = [1,1,1]
 690             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 691             >>> print(SLG._c())
 692             [-1.0000000]
 693             [ 0.0000000]
 694             [ 0.0000000]
 695             [ 0.0000000]
 696             <BLANKLINE>
 697
 698         """
 699         return matrix([-1, self._zero()])
 700
 701
 702     def C(self):
 703         """
 704         Return the cone ``C`` used in our CVXOPT construction.
 705
 706         The cone ``C`` is the cone over which the conelp program takes
 707         place.
 708
 709         Returns
 710         -------
 711
 712         CartesianProduct
 713             The cartesian product of ``K`` with itself.
 714
 715         Examples
 716         --------
 717
 718             >>> from dunshire import *
 719             >>> K = NonnegativeOrthant(3)
 720             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 721             >>> e1 = [1,2,3]
 722             >>> e2 = [1,1,1]
 723             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 724             >>> print(SLG.C())
 725             Cartesian product of dimension 6 with 2 factors:
 726               * Nonnegative orthant in the real 3-space
 727               * Nonnegative orthant in the real 3-space
 728
 729         """
 730         return CartesianProduct(self._K, self._K)
 731
 732     def _h(self):
 733         """
 734         Return the ``h`` vector used in our CVXOPT construction.
 735
 736         The ``h`` vector appears on the right-hand side of :math:`Gx + s
 737         = h` in the statement of the CVXOPT conelp program.
 738
 739         .. warning::
 740
 741             It is not safe to cache any of the matrices passed to
 742             CVXOPT, because it can clobber them.
 743
 744         Returns
 745         -------
 746
 747         matrix
 748             A ``2*self.dimension()``-by-``1`` column vector of zeros.
 749
 750         Examples
 751         --------
 752
 753             >>> from dunshire import *
 754             >>> K = NonnegativeOrthant(3)
 755             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 756             >>> e1 = [1,2,3]
 757             >>> e2 = [1,1,1]
 758             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 759             >>> print(SLG._h())
 760             [0.0000000]
 761             [0.0000000]
 762             [0.0000000]
 763             [0.0000000]
 764             [0.0000000]
 765             [0.0000000]
 766             <BLANKLINE>
 767
 768         """
 769
 770         return matrix([self._zero(), self._zero()])
 771
 772
 773     @staticmethod
 774     def b():
 775         """
 776         Return the ``b`` vector used in our CVXOPT construction.
 777
 778         The vector ``b`` appears on the right-hand side of :math:`Ax =
 779         b` in the statement of the CVXOPT conelp program.
 780
 781         This method is static because the dimensions and entries of
 782         ``b`` are known beforehand, and don't depend on any other
 783         properties of the game.
 784
 785         .. warning::
 786
 787             It is not safe to cache any of the matrices passed to
 788             CVXOPT, because it can clobber them.
 789
 790         Returns
 791         -------
 792
 793         matrix
 794             A ``1``-by-``1`` matrix containing a single entry ``1``.
 795
 796         Examples
 797         --------
 798
 799             >>> from dunshire import *
 800             >>> K = NonnegativeOrthant(3)
 801             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 802             >>> e1 = [1,2,3]
 803             >>> e2 = [1,1,1]
 804             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 805             >>> print(SLG.b())
 806             [1.0000000]
 807             <BLANKLINE>
 808
 809         """
 810         return matrix([1], tc='d')
 811
 812
 813     def player1_start(self):
 814         """
 815         Return a feasible starting point for player one.
 816
 817         This starting point is for the CVXOPT formulation and not for
 818         the original game. The basic premise is that if you normalize
 819         :meth:`e2`, then you get a point in :meth:`K` that makes a unit
 820         inner product with :meth:`e2`. We then get to choose the primal
 821         objective function value such that the constraint involving
 822         :meth:`L` is satisfied.
 823         """
 824         p = self.e2() / (norm(self.e2()) ** 2)
 825
 826         # Compute the distance from p to the outside of K.
 827         if isinstance(self.K(), NonnegativeOrthant):
 828             # How far is it to a wall?
 829             dist = min(list(self.e1()))
 830         elif isinstance(self.K(), IceCream):
 831             # How far is it to the boundary of the ball that defines
 832             # the ice-cream cone at a given height? Now draw a
 833             # 45-45-90 triangle and the shortest distance to the
 834             # outside of the cone should be 1/sqrt(2) of that.
 835             # It works in R^2, so it works everywhere, right?
 836             height = self.e1()[0]
 837             radius = norm(self.e1()[1:])
 838             dist = (height - radius) / sqrt(2)
 839         else:
 840             raise NotImplementedError
 841
 842         nu = - specnorm(self.L())/(dist*norm(self.e2()))
 843         x = matrix([nu,p], (self.dimension() + 1, 1))
 844         s = - self._G()*x
 845
 846         return {'x': x, 's': s}
 847
 848
 849     def player2_start(self):
 850         """
 851         Return a feasible starting point for player two.
 852         """
 853         q = self.e1() / (norm(self.e1()) ** 2)
 854
 855         # Compute the distance from p to the outside of K.
 856         if isinstance(self.K(), NonnegativeOrthant):
 857             # How far is it to a wall?
 858             dist = min(list(self.e2()))
 859         elif isinstance(self.K(), IceCream):
 860             # How far is it to the boundary of the ball that defines
 861             # the ice-cream cone at a given height? Now draw a
 862             # 45-45-90 triangle and the shortest distance to the
 863             # outside of the cone should be 1/sqrt(2) of that.
 864             # It works in R^2, so it works everywhere, right?
 865             height = self.e2()[0]
 866             radius = norm(self.e2()[1:])
 867             dist = (height - radius) / sqrt(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())