X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=dunshire%2Fgames.py;h=c2ca68ea7faee4f751410a1dbf27c31517400d9f;hb=aea70b13ddca62b0f6294083fd8130495d9a331e;hp=80f5f8ae2a46d09f15e34546c511edeca286f4fa;hpb=bdb596b84a06d0c97e39d42586a51fc36ba44186;p=dunshire.git diff --git a/dunshire/games.py b/dunshire/games.py index 80f5f8a..c2ca68e 100644 --- a/dunshire/games.py +++ b/dunshire/games.py @@ -7,13 +7,11 @@ knows how to solve a linear game. from cvxopt import matrix, printing, solvers from .cones import CartesianProduct -from .errors import GameUnsolvableException -from .matrices import append_col, append_row, identity +from .errors import GameUnsolvableException, PoorScalingException +from .matrices import append_col, append_row, condition_number, identity from . import options printing.options['dformat'] = options.FLOAT_FORMAT -solvers.options['show_progress'] = options.VERBOSE - class Solution: """ @@ -221,7 +219,8 @@ class SymmetricLinearGame: [ 1], e2 = [ 1] [ 2] - [ 3]. + [ 3], + Condition((L, K, e1, e2)) = 31.834... Lists can (and probably should) be used for every argument:: @@ -239,7 +238,8 @@ class SymmetricLinearGame: e1 = [ 1] [ 1], e2 = [ 1] - [ 1]. + [ 1], + Condition((L, K, e1, e2)) = 1.707... The points ``e1`` and ``e2`` can also be passed as some other enumerable type (of the correct length) without much harm, since @@ -261,7 +261,8 @@ class SymmetricLinearGame: e1 = [ 1] [ 1], e2 = [ 1] - [ 1]. + [ 1], + Condition((L, K, e1, e2)) = 1.707... However, ``L`` will always be intepreted as a list of rows, even if it is passed as a :class:`cvxopt.base.matrix` which is @@ -282,7 +283,8 @@ class SymmetricLinearGame: e1 = [ 1] [ 1], e2 = [ 1] - [ 1]. + [ 1], + Condition((L, K, e1, e2)) = 6.073... >>> L = cvxopt.matrix(L) >>> print(L) [ 1 3] @@ -297,7 +299,8 @@ class SymmetricLinearGame: e1 = [ 1] [ 1], e2 = [ 1] - [ 1]. + [ 1], + Condition((L, K, e1, e2)) = 6.073... """ def __init__(self, L, K, e1, e2): @@ -319,6 +322,10 @@ class SymmetricLinearGame: if not self._e2 in K: raise ValueError('the point e2 must lie in the interior of K') + # Cached result of the self._zero() method. + self._zero_col = None + + def __str__(self): """ Return a string representation of this game. @@ -327,11 +334,51 @@ class SymmetricLinearGame: ' L = {:s},\n' \ ' K = {!s},\n' \ ' e1 = {:s},\n' \ - ' e2 = {:s}.' + ' 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()) - return tpl.format(indented_L, str(self._K), indented_e1, indented_e2) + + return tpl.format(indented_L, + str(self._K), + indented_e1, + indented_e2, + self.condition()) + + + def _zero(self): + """ + Return a column of zeros that fits ``K``. + + This is used in our CVXOPT construction. + """ + if self._zero_col is None: + # Cache it, it's constant. + self._zero_col = matrix(0, (self._K.dimension(), 1), tc='d') + return self._zero_col + + + def _A(self): + """ + 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. + """ + return matrix([0, self._e2], (1, self._K.dimension() + 1), 'd') + + + def _G(self): + r""" + Return the matrix ``G`` used in our CVXOPT construction. + + Thus matrix ``G``that appears on the left-hand side of ``Gx + s = h`` + in the statement of the CVXOPT conelp program. + """ + I = identity(self._K.dimension()) + return append_row(append_col(self._zero(), -I), + append_col(self._e1, -self._L)) def solution(self): @@ -351,6 +398,10 @@ class SymmetricLinearGame: If the game could not be solved (if an optimal solution to its associated cone program was not found). + PoorScalingException + If the game could not be solved because CVXOPT crashed while + trying to take the square root of a negative number. + Examples -------- @@ -366,13 +417,13 @@ class SymmetricLinearGame: >>> print(SLG.solution()) Game value: -6.1724138 Player 1 optimal: - [ 0.5517241] - [-0.0000000] - [ 0.4482759] + [ 0.551...] + [-0.000...] + [ 0.448...] Player 2 optimal: - [0.4482759] - [0.0000000] - [0.5517241] + [0.448...] + [0.000...] + [0.551...] The value of the following game can be computed using the fact that the identity is invertible:: @@ -386,13 +437,13 @@ class SymmetricLinearGame: >>> print(SLG.solution()) Game value: 0.0312500 Player 1 optimal: - [0.0312500] - [0.0625000] - [0.0937500] + [0.031...] + [0.062...] + [0.093...] Player 2 optimal: - [0.1250000] - [0.1562500] - [0.1875000] + [0.125...] + [0.156...] + [0.187...] """ # The cone "C" that appears in the statement of the CVXOPT @@ -403,29 +454,29 @@ class SymmetricLinearGame: # Ax = b in the statement of the CVXOPT conelp program. b = matrix([1], tc='d') - # A column of zeros that fits K. - zero = matrix(0, (self._K.dimension(), 1), tc='d') - # The column vector "h" that appears on the right-hand side of # Gx + s = h in the statement of the CVXOPT conelp program. - h = matrix([zero, zero]) + h = matrix([self._zero(), self._zero()]) # The column vector "c" that appears in the objective function # value in the statement of the CVXOPT conelp program. - c = matrix([-1, zero]) - - # The matrix "G" that appears on the left-hand side of Gx + s = h - # in the statement of the CVXOPT conelp program. - G = append_row(append_col(zero, -identity(self._K.dimension())), - append_col(self._e1, -self._L)) - - # The matrix "A" that appears on the right-hand side of Ax = b - # in the statement of the CVXOPT conelp program. - A = matrix([0, self._e2], (1, self._K.dimension() + 1), 'd') + c = matrix([-1, self._zero()]) # Actually solve the thing and obtain a dictionary describing # what happened. - soln_dict = solvers.conelp(c, G, h, C.cvxopt_dims(), A, b) + try: + solvers.options['show_progress'] = options.VERBOSE + solvers.options['abs_tol'] = options.ABS_TOL + soln_dict = solvers.conelp(c, self._G(), h, + C.cvxopt_dims(), self._A(), b) + except ValueError as e: + if str(e) == '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 e # The optimal strategies are named ``p`` and ``q`` in the # background documentation, and we need to extract them from @@ -454,6 +505,12 @@ class SymmetricLinearGame: # 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) > options.ABS_TOL: raise GameUnsolvableException(self, soln_dict) if (p1_optimal not in self._K) or (p2_optimal not in self._K): @@ -462,6 +519,43 @@ class SymmetricLinearGame: return Solution(p1_value, p1_optimal, p2_optimal) + def condition(self): + r""" + Return the condition number of this game. + + In the CVXOPT construction of this game, two matrices ``G`` and + ``A`` appear. When those matrices are nasty, numerical problems + 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`). + + Returns + ------- + + float + A real number greater than or equal to one that measures how + bad this game is numerically. + + Examples + -------- + + >>> from dunshire import * + >>> K = NonnegativeOrthant(1) + >>> L = [[1]] + >>> e1 = [1] + >>> e2 = e1 + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> actual = SLG.condition() + >>> expected = 1.8090169943749477 + >>> abs(actual - expected) < options.ABS_TOL + True + + """ + return (condition_number(self._G()) + condition_number(self._A()))/2 + + def dual(self): r""" Return the dual game to this game. @@ -491,7 +585,8 @@ class SymmetricLinearGame: [ 3], e2 = [ 1] [ 1] - [ 1]. + [ 1], + Condition((L, K, e1, e2)) = 44.476... """ # We pass ``self._L`` right back into the constructor, because