X-Git-Url: https://gitweb.michael.orlitzky.com/?p=dunshire.git;a=blobdiff_plain;f=dunshire%2Fgames.py;h=ea7a64f6b8e6451a808b464494c11e9be9f0de78;hp=08a6d6f342b59f93d26fbdeae320034f95004321;hb=0274de467062ab29d2a41d2a91ec0b28fcd95c8d;hpb=eb940c4312540b7db2f43a8c4067aff6e94928de diff --git a/dunshire/games.py b/dunshire/games.py index 08a6d6f..ea7a64f 100644 --- a/dunshire/games.py +++ b/dunshire/games.py @@ -4,14 +4,15 @@ Symmetric linear games and their solutions. This module contains the main :class:`SymmetricLinearGame` class that knows how to solve a linear game. """ - from cvxopt import matrix, printing, solvers from .cones import CartesianProduct from .errors import GameUnsolvableException, PoorScalingException -from .matrices import append_col, append_row, condition_number, identity -from . import options +from .matrices import (append_col, append_row, condition_number, identity, + inner_product, norm, specnorm) +from .options import ABS_TOL, FLOAT_FORMAT, DEBUG_FLOAT_FORMAT + +printing.options['dformat'] = FLOAT_FORMAT -printing.options['dformat'] = options.FLOAT_FORMAT class Solution: """ @@ -22,7 +23,7 @@ class Solution: -------- >>> print(Solution(10, matrix([1,2]), matrix([3,4]))) - Game value: 10.0000000 + Game value: 10.000... Player 1 optimal: [ 1] [ 2] @@ -178,11 +179,15 @@ class SymmetricLinearGame: ---------- L : list of list of float - A matrix represented as a list of ROWS. This representation - agrees with (for example) SageMath and NumPy, but not with CVXOPT - (whose matrix constructor accepts a list of columns). - - K : :class:`SymmetricCone` + A matrix represented as a list of **rows**. This representation + agrees with (for example) `SageMath `_ + and `NumPy `_, but not with CVXOPT (whose + matrix constructor accepts a list of columns). In reality, ``L`` + can be any iterable type of the correct length; however, you + should be extremely wary of the way we interpret anything other + than a list of rows. + + K : dunshire.cones.SymmetricCone The symmetric cone instance over which the game is played. e1 : iterable float @@ -219,8 +224,7 @@ class SymmetricLinearGame: [ 1], e2 = [ 1] [ 2] - [ 3], - Condition((L, K, e1, e2)) = 31.834... + [ 3] Lists can (and probably should) be used for every argument:: @@ -238,8 +242,7 @@ class SymmetricLinearGame: e1 = [ 1] [ 1], e2 = [ 1] - [ 1], - Condition((L, K, e1, e2)) = 1.707... + [ 1] The points ``e1`` and ``e2`` can also be passed as some other enumerable type (of the correct length) without much harm, since @@ -261,8 +264,7 @@ class SymmetricLinearGame: e1 = [ 1] [ 1], e2 = [ 1] - [ 1], - Condition((L, K, e1, e2)) = 1.707... + [ 1] However, ``L`` will always be intepreted as a list of rows, even if it is passed as a :class:`cvxopt.base.matrix` which is @@ -283,8 +285,7 @@ class SymmetricLinearGame: e1 = [ 1] [ 1], e2 = [ 1] - [ 1], - Condition((L, K, e1, e2)) = 6.073... + [ 1] >>> L = cvxopt.matrix(L) >>> print(L) [ 1 3] @@ -299,8 +300,7 @@ class SymmetricLinearGame: e1 = [ 1] [ 1], e2 = [ 1] - [ 1], - Condition((L, K, e1, e2)) = 6.073... + [ 1] """ def __init__(self, L, K, e1, e2): @@ -322,6 +322,8 @@ class SymmetricLinearGame: if not self._e2 in K: raise ValueError('the point e2 must lie in the interior of K') + # Initial value of cached method. + self._L_specnorm_value = None def __str__(self): @@ -332,17 +334,213 @@ class SymmetricLinearGame: ' L = {:s},\n' \ ' K = {!s},\n' \ ' e1 = {:s},\n' \ - ' e2 = {:s},\n' \ - ' Condition((L, K, e1, e2)) = {:f}.' - indented_L = '\n '.join(str(self._L).splitlines()) - indented_e1 = '\n '.join(str(self._e1).splitlines()) - indented_e2 = '\n '.join(str(self._e2).splitlines()) + ' e2 = {:s}' + indented_L = '\n '.join(str(self.L()).splitlines()) + indented_e1 = '\n '.join(str(self.e1()).splitlines()) + indented_e2 = '\n '.join(str(self.e2()).splitlines()) return tpl.format(indented_L, - str(self._K), + str(self.K()), indented_e1, - indented_e2, - self.condition()) + indented_e2) + + + def L(self): + """ + Return the matrix ``L`` passed to the constructor. + + Returns + ------- + + matrix + The matrix that defines this game's :meth:`payoff` operator. + + Examples + -------- + + >>> from dunshire import * + >>> K = NonnegativeOrthant(3) + >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]] + >>> e1 = [1,1,1] + >>> e2 = [1,2,3] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG.L()) + [ 1 -5 -15] + [ -1 2 -3] + [-12 -15 1] + + + """ + return self._L + + + def K(self): + """ + Return the cone over which this game is played. + + Returns + ------- + + SymmetricCone + The :class:`SymmetricCone` over which this game is played. + + Examples + -------- + + >>> from dunshire import * + >>> K = NonnegativeOrthant(3) + >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]] + >>> e1 = [1,1,1] + >>> e2 = [1,2,3] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG.K()) + Nonnegative orthant in the real 3-space + + """ + return self._K + + + def e1(self): + """ + Return player one's interior point. + + Returns + ------- + + matrix + The point interior to :meth:`K` affiliated with player one. + + Examples + -------- + + >>> from dunshire import * + >>> K = NonnegativeOrthant(3) + >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]] + >>> e1 = [1,1,1] + >>> e2 = [1,2,3] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG.e1()) + [ 1] + [ 1] + [ 1] + + + """ + return self._e1 + + + def e2(self): + """ + Return player two's interior point. + + Returns + ------- + + matrix + The point interior to :meth:`K` affiliated with player one. + + Examples + -------- + + >>> from dunshire import * + >>> K = NonnegativeOrthant(3) + >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]] + >>> e1 = [1,1,1] + >>> e2 = [1,2,3] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG.e2()) + [ 1] + [ 2] + [ 3] + + + """ + return self._e2 + + + def payoff(self, strategy1, strategy2): + r""" + Return the payoff associated with ``strategy1`` and ``strategy2``. + + The payoff operator takes pairs of strategies to a real + number. For example, if player one's strategy is :math:`x` and + player two's strategy is :math:`y`, then the associated payoff + is :math:`\left\langle L\left(x\right),y \right\rangle \in + \mathbb{R}`. Here, :math:`L` denotes the same linear operator as + :meth:`L`. This method computes the payoff given the two + players' strategies. + + Parameters + ---------- + + strategy1 : matrix + Player one's strategy. + + strategy2 : matrix + Player two's strategy. + + Returns + ------- + + float + The payoff for the game when player one plays ``strategy1`` + and player two plays ``strategy2``. + + Examples + -------- + + The value of the game should be the payoff at the optimal + strategies:: + + >>> from dunshire import * + >>> K = NonnegativeOrthant(3) + >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]] + >>> e1 = [1,1,1] + >>> e2 = [1,1,1] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> soln = SLG.solution() + >>> x_bar = soln.player1_optimal() + >>> y_bar = soln.player2_optimal() + >>> SLG.payoff(x_bar, y_bar) == soln.game_value() + True + + """ + return inner_product(self.L()*strategy1, strategy2) + + + def dimension(self): + """ + Return the dimension of this game. + + The dimension of a game is not needed for the theory, but it is + useful for the implementation. We define the dimension of a game + to be the dimension of its underlying cone. Or what is the same, + the dimension of the space from which the strategies are chosen. + + Returns + ------- + + int + The dimension of the cone :meth:`K`, or of the space where + this game is played. + + Examples + -------- + + The dimension of a game over the nonnegative quadrant in the + plane should be two (the dimension of the plane):: + + >>> from dunshire import * + >>> K = NonnegativeOrthant(2) + >>> L = [[1,-5],[-1,2]] + >>> e1 = [1,1] + >>> e2 = [1,4] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> SLG.dimension() + 2 + + """ + return self.K().dimension() def _zero(self): @@ -360,7 +558,7 @@ class SymmetricLinearGame: ------- matrix - A ``K.dimension()``-by-``1`` column vector of zeros. + A ``self.dimension()``-by-``1`` column vector of zeros. Examples -------- @@ -378,15 +576,16 @@ class SymmetricLinearGame: """ - return matrix(0, (self._K.dimension(), 1), tc='d') + return matrix(0, (self.dimension(), 1), tc='d') - def _A(self): - """ + def A(self): + r""" Return the matrix ``A`` used in our CVXOPT construction. - This matrix ``A`` appears on the right-hand side of ``Ax = b`` - in the statement of the CVXOPT conelp program. + This matrix :math:`A` appears on the right-hand side of + :math:`Ax = b` in the `statement of the CVXOPT conelp program + `_. .. warning:: @@ -397,8 +596,8 @@ class SymmetricLinearGame: ------- matrix - A ``1``-by-``(1 + K.dimension())`` row vector. Its first - entry is zero, and the rest are the entries of ``e2``. + A ``1``-by-``(1 + self.dimension())`` row vector. Its first + entry is zero, and the rest are the entries of :meth:`e2`. Examples -------- @@ -409,21 +608,22 @@ class SymmetricLinearGame: >>> e1 = [1,1,1] >>> e2 = [1,2,3] >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> print(SLG._A()) + >>> print(SLG.A()) [0.0000000 1.0000000 2.0000000 3.0000000] """ - return matrix([0, self._e2], (1, self._K.dimension() + 1), 'd') + return matrix([0, self.e2()], (1, self.dimension() + 1), 'd') - def _G(self): + def G(self): r""" Return the matrix ``G`` used in our CVXOPT construction. - Thus matrix ``G`` appears on the left-hand side of ``Gx + s = h`` - in the statement of the CVXOPT conelp program. + Thus matrix :math:`G` appears on the left-hand side of :math:`Gx + + s = h` in the `statement of the CVXOPT conelp program + `_. .. warning:: @@ -434,7 +634,7 @@ class SymmetricLinearGame: ------- matrix - A ``2*K.dimension()``-by-``1 + K.dimension()`` matrix. + A ``2*self.dimension()``-by-``(1 + self.dimension())`` matrix. Examples -------- @@ -445,7 +645,7 @@ class SymmetricLinearGame: >>> e1 = [1,2,3] >>> e2 = [1,1,1] >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> print(SLG._G()) + >>> print(SLG.G()) [ 0.0000000 -1.0000000 0.0000000 0.0000000] [ 0.0000000 0.0000000 -1.0000000 0.0000000] [ 0.0000000 0.0000000 0.0000000 -1.0000000] @@ -455,17 +655,19 @@ class SymmetricLinearGame: """ - I = identity(self._K.dimension()) - return append_row(append_col(self._zero(), -I), - append_col(self._e1, -self._L)) + identity_matrix = identity(self.dimension()) + return append_row(append_col(self._zero(), -identity_matrix), + append_col(self.e1(), -self.L())) - def _c(self): - """ + def c(self): + r""" Return the vector ``c`` used in our CVXOPT construction. - The column vector ``c`` appears in the objective function - value ```` in the statement of the CVXOPT conelp program. + The column vector :math:`c` appears in the objective function + value :math:`\left\langle c,x \right\rangle` in the `statement + of the CVXOPT conelp program + `_. .. warning:: @@ -476,7 +678,7 @@ class SymmetricLinearGame: ------- matrix - A ``K.dimension()``-by-``1`` column vector. + A :meth:`dimension`-by-``1`` column vector. Examples -------- @@ -487,7 +689,7 @@ class SymmetricLinearGame: >>> e1 = [1,2,3] >>> e2 = [1,1,1] >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> print(SLG._c()) + >>> print(SLG.c()) [-1.0000000] [ 0.0000000] [ 0.0000000] @@ -498,12 +700,13 @@ class SymmetricLinearGame: return matrix([-1, self._zero()]) - def _C(self): + def C(self): """ Return the cone ``C`` used in our CVXOPT construction. - The cone ``C`` is the cone over which the conelp program takes - place. + This is the cone over which the `CVXOPT conelp program + `_ + takes place. Returns ------- @@ -520,7 +723,7 @@ class SymmetricLinearGame: >>> e1 = [1,2,3] >>> e2 = [1,1,1] >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> print(SLG._C()) + >>> print(SLG.C()) Cartesian product of dimension 6 with 2 factors: * Nonnegative orthant in the real 3-space * Nonnegative orthant in the real 3-space @@ -528,12 +731,13 @@ class SymmetricLinearGame: """ return CartesianProduct(self._K, self._K) - def _h(self): - """ + def h(self): + r""" Return the ``h`` vector used in our CVXOPT construction. - The ``h`` vector appears on the right-hand side of :math:`Gx + s - = h` in the statement of the CVXOPT conelp program. + The :math:`h` vector appears on the right-hand side of :math:`Gx + + s = h` in the `statement of the CVXOPT conelp program + `_. .. warning:: @@ -544,7 +748,7 @@ class SymmetricLinearGame: ------- matrix - A ``2*K.dimension()``-by-``1`` column vector of zeros. + A ``2*self.dimension()``-by-``1`` column vector of zeros. Examples -------- @@ -555,7 +759,7 @@ class SymmetricLinearGame: >>> e1 = [1,2,3] >>> e2 = [1,1,1] >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> print(SLG._h()) + >>> print(SLG.h()) [0.0000000] [0.0000000] [0.0000000] @@ -568,12 +772,19 @@ class SymmetricLinearGame: return matrix([self._zero(), self._zero()]) - def _b(self): - """ + + @staticmethod + def b(): + r""" Return the ``b`` vector used in our CVXOPT construction. - The vector ``b`` appears on the right-hand side of :math:`Ax = - b` in the statement of the CVXOPT conelp program. + The vector :math:`b` appears on the right-hand side of :math:`Ax + = b` in the `statement of the CVXOPT conelp program + `_. + + This method is static because the dimensions and entries of + ``b`` are known beforehand, and don't depend on any other + properties of the game. .. warning:: @@ -595,7 +806,7 @@ class SymmetricLinearGame: >>> e1 = [1,2,3] >>> e2 = [1,1,1] >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> print(SLG._b()) + >>> print(SLG.b()) [1.0000000] @@ -603,152 +814,190 @@ class SymmetricLinearGame: return matrix([1], tc='d') - def _try_solution(self, tolerance): + def player1_start(self): """ - Solve this linear game within ``tolerance``, if possible. + Return a feasible starting point for player one. - This private function is the one that does all of the actual - work for :meth:`solution`. This method accepts a ``tolerance``, - and what :meth:`solution` does is call this method twice with - two different tolerances. First it tries a strict tolerance, and - then it tries a looser one. + This starting point is for the CVXOPT formulation and not for + the original game. The basic premise is that if you scale + :meth:`e2` by the reciprocal of its squared norm, then you get a + point in :meth:`K` that makes a unit inner product with + :meth:`e2`. We then get to choose the primal objective function + value such that the constraint involving :meth:`L` is satisfied. - .. warning:: + Returns + ------- - If you try to be smart and precompute the matrices used by - this function (the ones passed to ``conelp``), then you're - going to shoot yourself in the foot. CVXOPT can and will - clobber some (but not all) of its input matrices. This isn't - performance sensitive, so play it safe. + dict + A dictionary with two keys, ``'x'`` and ``'s'``, which + contain the vectors of the same name in the CVXOPT primal + problem formulation. - Parameters - ---------- + The vector ``x`` consists of the primal objective function + value concatenated with the strategy (for player one) that + achieves it. The vector ``s`` is essentially a dummy + variable, and is computed from the equality constraing in + the CVXOPT primal problem. + + """ + p = self.e2() / (norm(self.e2()) ** 2) + dist = self.K().ball_radius(self.e1()) + nu = - self._L_specnorm()/(dist*norm(self.e2())) + x = matrix([nu, p], (self.dimension() + 1, 1)) + s = - self.G()*x - tolerance : float - The absolute tolerance to pass to the CVXOPT solver. + return {'x': x, 's': s} + + + def player2_start(self): + """ + Return a feasible starting point for player two. + + This starting point is for the CVXOPT formulation and not for + the original game. The basic premise is that if you scale + :meth:`e1` by the reciprocal of its squared norm, then you get a + point in :meth:`K` that makes a unit inner product with + :meth:`e1`. We then get to choose the dual objective function + value such that the constraint involving :meth:`L` is satisfied. Returns ------- - :class:`Solution` - A :class:`Solution` object describing the game's value and - the optimal strategies of both players. + dict + A dictionary with two keys, ``'y'`` and ``'z'``, which + contain the vectors of the same name in the CVXOPT dual + problem formulation. - Raises - ------ - GameUnsolvableException - If the game could not be solved (if an optimal solution to its - associated cone program was not found). + The ``1``-by-``1`` vector ``y`` consists of the dual + objective function value. The last :meth:`dimension` entries + of the vector ``z`` contain the strategy (for player two) + that achieves it. The remaining entries of ``z`` are + essentially dummy variables, computed from the equality + constraint in the CVXOPT dual problem. - PoorScalingException - If the game could not be solved because CVXOPT crashed while - trying to take the square root of a negative number. + """ + q = self.e1() / (norm(self.e1()) ** 2) + dist = self.K().ball_radius(self.e2()) + omega = self._L_specnorm()/(dist*norm(self.e1())) + y = matrix([omega]) + z2 = q + z1 = y*self.e2() - self.L().trans()*z2 + z = matrix([z1, z2], (self.dimension()*2, 1)) + + return {'y': y, 'z': z} + + + def _L_specnorm(self): + """ + Compute the spectral norm of :meth:`L` and cache it. + + The spectral norm of the matrix :meth:`L` is used in a few + places. Since it can be expensive to compute, we want to cache + its value. That is not possible in :func:`specnorm`, which lies + outside of a class, so this is the place to do it. + + Returns + ------- + + float + A nonnegative real number; the largest singular value of + the matrix :meth:`L`. Examples -------- - This game can be solved easily, so the first attempt in - :meth:`solution` should succeed:: - >>> from dunshire import * - >>> from dunshire.matrices import norm - >>> from dunshire.options import ABS_TOL - >>> K = NonnegativeOrthant(3) - >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]] - >>> e1 = [1,1,1] - >>> e2 = [1,1,1] - >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> s1 = SLG.solution() - >>> s2 = SLG._try_solution(options.ABS_TOL) - >>> abs(s1.game_value() - s2.game_value()) < ABS_TOL - True - >>> norm(s1.player1_optimal() - s2.player1_optimal()) < ABS_TOL - True - >>> norm(s1.player2_optimal() - s2.player2_optimal()) < ABS_TOL + >>> from dunshire.matrices import specnorm + >>> L = [[1,2],[3,4]] + >>> K = NonnegativeOrthant(2) + >>> e1 = [1,1] + >>> e2 = e1 + >>> SLG = SymmetricLinearGame(L,K,e1,e2) + >>> specnorm(SLG.L()) == SLG._L_specnorm() True - This game cannot be solved with the default tolerance, but it - can be solved with a weaker one:: + """ + if self._L_specnorm_value is None: + self._L_specnorm_value = specnorm(self.L()) + return self._L_specnorm_value + + + def tolerance_scale(self, solution): + r""" + + Return a scaling factor that should be applied to + :const:`dunshire.options.ABS_TOL` for this game. + + When performing certain comparisons, the default tolerance + :const:`dunshire.options.ABS_TOL` may not be appropriate. For + example, if we expect ``x`` and ``y`` to be within + :const:`dunshire.options.ABS_TOL` of each other, than the inner + product of ``L*x`` and ``y`` can be as far apart as the spectral + norm of ``L`` times the sum of the norms of ``x`` and + ``y``. Such a comparison is made in :meth:`solution`, and in + many of our unit tests. + + The returned scaling factor found from the inner product + mentioned above is + + .. math:: + + \left\lVert L \right\rVert_{2} + \left( \left\lVert \bar{x} \right\rVert + + \left\lVert \bar{y} \right\rVert + \right), + + where :math:`\bar{x}` and :math:`\bar{y}` are optimal solutions + for players one and two respectively. This scaling factor is not + formally justified, but attempting anything smaller leads to + test failures. + + .. warning:: + + Optimal solutions are not unique, so the scaling factor + obtained from ``solution`` may not work when comparing other + solutions. + + Parameters + ---------- + + solution : Solution + A solution of this game, used to obtain the norms of the + optimal strategies. + + Returns + ------- + + float + A scaling factor to be multiplied by + :const:`dunshire.options.ABS_TOL` when + making comparisons involving solutions of this game. + + Examples + -------- + + The spectral norm of ``L`` in this case is around ``5.464``, and + the optimal strategies both have norm one, so we expect the + tolerance scale to be somewhere around ``2 * 5.464``, or + ``10.929``:: >>> from dunshire import * - >>> from dunshire.options import ABS_TOL - >>> L = [[ 0.58538005706658102767, 1.53764301129883040886], - ... [-1.34901059721452210027, 1.50121179114155500756]] + >>> L = [[1,2],[3,4]] >>> K = NonnegativeOrthant(2) - >>> e1 = [1.04537193228494995623, 1.39699624965841895374] - >>> e2 = [0.35326554172108337593, 0.11795703527854853321] + >>> e1 = [1,1] + >>> e2 = e1 >>> SLG = SymmetricLinearGame(L,K,e1,e2) - >>> print(SLG._try_solution(ABS_TOL / 10)) - Traceback (most recent call last): - ... - dunshire.errors.GameUnsolvableException: Solution failed... - >>> print(SLG._try_solution(ABS_TOL)) - Game value: 9.1100945 - Player 1 optimal: - [-0.0000000] - [ 8.4776631] - Player 2 optimal: - [0.0000000] - [0.7158216] + >>> SLG.tolerance_scale(SLG.solution()) + 10.929... """ - try: - opts = {'show_progress': options.VERBOSE, 'abstol': tolerance} - soln_dict = solvers.conelp(self._c(), - self._G(), - self._h(), - self._C().cvxopt_dims(), - self._A(), - self._b(), - options=opts) - except ValueError as error: - if str(error) == 'math domain error': - # Oops, CVXOPT tried to take the square root of a - # negative number. Report some details about the game - # rather than just the underlying CVXOPT crash. - raise PoorScalingException(self) - else: - raise error - - # The optimal strategies are named ``p`` and ``q`` in the - # background documentation, and we need to extract them from - # the CVXOPT ``x`` and ``z`` variables. The objective values - # :math:`nu` and :math:`omega` can also be found in the CVXOPT - # ``x`` and ``y`` variables; however, they're stored - # conveniently as separate entries in the solution dictionary. - p1_value = -soln_dict['primal objective'] - p2_value = -soln_dict['dual objective'] - p1_optimal = soln_dict['x'][1:] - p2_optimal = soln_dict['z'][self._K.dimension():] + norm_p1_opt = norm(solution.player1_optimal()) + norm_p2_opt = norm(solution.player2_optimal()) + scale = self._L_specnorm()*(norm_p1_opt + norm_p2_opt) - # The "status" field contains "optimal" if everything went - # according to plan. Other possible values are "primal - # infeasible", "dual infeasible", "unknown", all of which mean - # we didn't get a solution. The "infeasible" ones are the - # worst, since they indicate that CVXOPT is convinced the - # problem is infeasible (and that cannot happen). - if soln_dict['status'] in ['primal infeasible', 'dual infeasible']: - raise GameUnsolvableException(self, soln_dict) - elif soln_dict['status'] == 'unknown': - # When we get a status of "unknown", we may still be able - # to salvage a solution out of the returned - # dictionary. Often this is the result of numerical - # difficulty and we can simply check that the primal/dual - # objectives match (within a tolerance) and that the - # primal/dual optimal solutions are within the cone (to a - # tolerance as well). - # - # The fudge factor of two is basically unjustified, but - # makes intuitive sense when you imagine that the primal - # value could be under the true optimal by ``ABS_TOL`` - # and the dual value could be over by the same amount. - # - if abs(p1_value - p2_value) > tolerance: - raise GameUnsolvableException(self, soln_dict) - if (p1_optimal not in self._K) or (p2_optimal not in self._K): - raise GameUnsolvableException(self, soln_dict) - - return Solution(p1_value, p1_optimal, p2_optimal) + # Don't return anything smaller than 1... we can't go below + # out "minimum tolerance." + return max(1, scale) def solution(self): @@ -758,7 +1007,7 @@ class SymmetricLinearGame: Returns ------- - :class:`Solution` + Solution A :class:`Solution` object describing the game's value and the optimal strategies of both players. @@ -785,11 +1034,11 @@ class SymmetricLinearGame: >>> e2 = [1,1,1] >>> SLG = SymmetricLinearGame(L, K, e1, e2) >>> print(SLG.solution()) - Game value: -6.1724138 + Game value: -6.172... Player 1 optimal: - [ 0.551...] - [-0.000...] - [ 0.448...] + [0.551...] + [0.000...] + [0.448...] Player 2 optimal: [0.448...] [0.000...] @@ -805,7 +1054,7 @@ class SymmetricLinearGame: >>> e2 = [4,5,6] >>> SLG = SymmetricLinearGame(L, K, e1, e2) >>> print(SLG.solution()) - Game value: 0.0312500 + Game value: 0.031... Player 1 optimal: [0.031...] [0.062...] @@ -815,16 +1064,177 @@ class SymmetricLinearGame: [0.156...] [0.187...] + This is another Gowda/Ravindran example that is supposed to have + a negative game value:: + + >>> from dunshire import * + >>> from dunshire.options import ABS_TOL + >>> L = [[1, -2], [-2, 1]] + >>> K = NonnegativeOrthant(2) + >>> e1 = [1, 1] + >>> e2 = e1 + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> SLG.solution().game_value() < -ABS_TOL + True + + The following two games are problematic numerically, but we + should be able to solve them:: + + >>> from dunshire import * + >>> L = [[-0.95237953890954685221, 1.83474556206462535712], + ... [ 1.30481749924621448500, 1.65278664543326403447]] + >>> K = NonnegativeOrthant(2) + >>> e1 = [0.95477167524644313001, 0.63270781756540095397] + >>> e2 = [0.39633793037154141370, 0.10239281495640320530] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG.solution()) + Game value: 18.767... + Player 1 optimal: + [0.000...] + [9.766...] + Player 2 optimal: + [1.047...] + [0.000...] + + :: + + >>> from dunshire import * + >>> L = [[1.54159395026049472754, 2.21344728574316684799], + ... [1.33147433507846657541, 1.17913616272988108769]] + >>> K = NonnegativeOrthant(2) + >>> e1 = [0.39903040089404784307, 0.12377403622479113410] + >>> e2 = [0.15695181142215544612, 0.85527381344651265405] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG.solution()) + Game value: 24.614... + Player 1 optimal: + [6.371...] + [0.000...] + Player 2 optimal: + [2.506...] + [0.000...] + + This is another one that was difficult numerically, and caused + trouble even after we fixed the first two:: + + >>> from dunshire import * + >>> L = [[57.22233908627052301199, 41.70631373437460354126], + ... [83.04512571985074487202, 57.82581810406928468637]] + >>> K = NonnegativeOrthant(2) + >>> e1 = [7.31887017043399268346, 0.89744171905822367474] + >>> e2 = [0.11099824781179848388, 6.12564670639315345113] + >>> SLG = SymmetricLinearGame(L,K,e1,e2) + >>> print(SLG.solution()) + Game value: 70.437... + Player 1 optimal: + [9.009...] + [0.000...] + Player 2 optimal: + [0.136...] + [0.000...] + + And finally, here's one that returns an "optimal" solution, but + whose primal/dual objective function values are far apart:: + + >>> from dunshire import * + >>> L = [[ 6.49260076597376212248, -0.60528030227678542019], + ... [ 2.59896077096751731972, -0.97685530240286766457]] + >>> K = IceCream(2) + >>> e1 = [1, 0.43749513972645248661] + >>> e2 = [1, 0.46008379832200291260] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG.solution()) + Game value: 11.596... + Player 1 optimal: + [ 1.852...] + [-1.852...] + Player 2 optimal: + [ 1.777...] + [-1.777...] + """ try: - # First try with a stricter tolerance. Who knows, it might - # work. If it does, we prefer that solution. - return self._try_solution(options.ABS_TOL / 10) + opts = {'show_progress': False} + soln_dict = solvers.conelp(self.c(), + self.G(), + self.h(), + self.C().cvxopt_dims(), + self.A(), + self.b(), + primalstart=self.player1_start(), + dualstart=self.player2_start(), + options=opts) + except ValueError as error: + if str(error) == 'math domain error': + # Oops, CVXOPT tried to take the square root of a + # negative number. Report some details about the game + # rather than just the underlying CVXOPT crash. + printing.options['dformat'] = DEBUG_FLOAT_FORMAT + raise PoorScalingException(self) + else: + raise error - except (PoorScalingException, GameUnsolvableException): - # Ok, that didn't work. Let's try it with the default - # tolerance, and whatever happens, happens. - return self._try_solution(options.ABS_TOL) + # The optimal strategies are named ``p`` and ``q`` in the + # background documentation, and we need to extract them from + # the CVXOPT ``x`` and ``z`` variables. The objective values + # :math:`nu` and :math:`omega` can also be found in the CVXOPT + # ``x`` and ``y`` variables; however, they're stored + # conveniently as separate entries in the solution dictionary. + p1_value = -soln_dict['primal objective'] + p2_value = -soln_dict['dual objective'] + p1_optimal = soln_dict['x'][1:] + p2_optimal = soln_dict['z'][self.dimension():] + + # The "status" field contains "optimal" if everything went + # according to plan. Other possible values are "primal + # infeasible", "dual infeasible", "unknown", all of which mean + # we didn't get a solution. + # + # The "infeasible" ones are the worst, since they indicate + # that CVXOPT is convinced the problem is infeasible (and that + # cannot happen). + if soln_dict['status'] in ['primal infeasible', 'dual infeasible']: + printing.options['dformat'] = DEBUG_FLOAT_FORMAT + raise GameUnsolvableException(self, soln_dict) + + # For the game value, we could use any of: + # + # * p1_value + # * p2_value + # * (p1_value + p2_value)/2 + # * the game payoff + # + # We want the game value to be the payoff, however, so it + # makes the most sense to just use that, even if it means we + # can't test the fact that p1_value/p2_value are close to the + # payoff. + payoff = self.payoff(p1_optimal, p2_optimal) + soln = Solution(payoff, p1_optimal, p2_optimal) + + # The "optimal" and "unknown" results, we actually treat the + # same. Even if CVXOPT bails out due to numerical difficulty, + # it will have some candidate points in mind. If those + # candidates are good enough, we take them. We do the same + # check for "optimal" results. + # + # First we check that the primal/dual objective values are + # close enough because otherwise CVXOPT might return "unknown" + # and give us two points in the cone that are nowhere near + # optimal. And in fact, we need to ensure that they're close + # for "optimal" results, too, because we need to know how + # lenient to be in our testing. + # + if abs(p1_value - p2_value) > self.tolerance_scale(soln)*ABS_TOL: + printing.options['dformat'] = DEBUG_FLOAT_FORMAT + raise GameUnsolvableException(self, soln_dict) + + # And we also check that the points it gave us belong to the + # cone, just in case... + if (p1_optimal not in self._K) or (p2_optimal not in self._K): + printing.options['dformat'] = DEBUG_FLOAT_FORMAT + raise GameUnsolvableException(self, soln_dict) + + return soln def condition(self): @@ -836,8 +1246,12 @@ class SymmetricLinearGame: can show up. We define the condition number of this game to be the average of the condition numbers of ``G`` and ``A`` in the CVXOPT construction. If the condition number of this game is - high, then you can expect numerical difficulty (such as - :class:`PoorScalingException`). + high, you can problems like :class:`PoorScalingException`. + + Random testing shows that a condition number of around ``125`` + is about the best that we can solve reliably. However, the + failures are intermittent, and you may get lucky with an + ill-conditioned game. Returns ------- @@ -855,13 +1269,11 @@ class SymmetricLinearGame: >>> e1 = [1] >>> e2 = e1 >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> actual = SLG.condition() - >>> expected = 1.8090169943749477 - >>> abs(actual - expected) < options.ABS_TOL - True + >>> SLG.condition() + 1.809... """ - return (condition_number(self._G()) + condition_number(self._A()))/2 + return (condition_number(self.G()) + condition_number(self.A()))/2 def dual(self): @@ -893,14 +1305,13 @@ class SymmetricLinearGame: [ 3], e2 = [ 1] [ 1] - [ 1], - Condition((L, K, e1, e2)) = 44.476... + [ 1] """ - # We pass ``self._L`` right back into the constructor, because + # We pass ``self.L()`` right back into the constructor, because # it will be transposed there. And keep in mind that ``self._K`` # is its own dual. - return SymmetricLinearGame(self._L, - self._K, - self._e2, - self._e1) + return SymmetricLinearGame(self.L(), + self.K(), + self.e2(), + self.e1())