X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=dunshire%2Fgames.py;h=6480d7d31153afcc419d331e0cd416a64253a188;hb=2a59d8c4a2d7723a52de15e3d1f8ffce76bb948a;hp=46092c380eca141ff993313bd30ec55989a32ed8;hpb=1e02fd12b64e090e0b0ab0d3fecbd9c1b18d0fcf;p=dunshire.git diff --git a/dunshire/games.py b/dunshire/games.py index 46092c3..6480d7d 100644 --- a/dunshire/games.py +++ b/dunshire/games.py @@ -12,8 +12,6 @@ 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: """ @@ -222,7 +220,7 @@ class SymmetricLinearGame: e2 = [ 1] [ 2] [ 3], - Condition((L, K, e1, e2)) = 63.669790. + Condition((L, K, e1, e2)) = 31.834... Lists can (and probably should) be used for every argument:: @@ -241,7 +239,7 @@ class SymmetricLinearGame: [ 1], e2 = [ 1] [ 1], - Condition((L, K, e1, e2)) = 3.414214. + 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 @@ -264,7 +262,7 @@ class SymmetricLinearGame: [ 1], e2 = [ 1] [ 1], - Condition((L, K, e1, e2)) = 3.414214. + 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 @@ -286,7 +284,7 @@ class SymmetricLinearGame: [ 1], e2 = [ 1] [ 1], - Condition((L, K, e1, e2)) = 12.147542. + Condition((L, K, e1, e2)) = 6.073... >>> L = cvxopt.matrix(L) >>> print(L) [ 1 3] @@ -302,7 +300,7 @@ class SymmetricLinearGame: [ 1], e2 = [ 1] [ 1], - Condition((L, K, e1, e2)) = 12.147542. + Condition((L, K, e1, e2)) = 6.073... """ def __init__(self, L, K, e1, e2): @@ -324,8 +322,6 @@ 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): @@ -346,7 +342,7 @@ class SymmetricLinearGame: str(self._K), indented_e1, indented_e2, - self._condition()) + self.condition()) def _zero(self): @@ -354,11 +350,35 @@ class SymmetricLinearGame: Return a column of zeros that fits ``K``. This is used in our CVXOPT construction. + + .. warning:: + + It is not safe to cache any of the matrices passed to + CVXOPT, because it can clobber them. + + Returns + ------- + + matrix + A ``K.dimension()``-by-``1`` column vector of zeros. + + Examples + -------- + + >>> from dunshire import * + >>> K = NonnegativeOrthant(3) + >>> L = identity(3) + >>> e1 = [1,1,1] + >>> e2 = e1 + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG._zero()) + [0.0000000] + [0.0000000] + [0.0000000] + + """ - 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 + return matrix(0, (self._K.dimension(), 1), tc='d') def _A(self): @@ -367,22 +387,136 @@ class SymmetricLinearGame: This matrix ``A`` appears on the right-hand side of ``Ax = b`` in the statement of the CVXOPT conelp program. + + .. warning:: + + It is not safe to cache any of the matrices passed to + CVXOPT, because it can clobber them. + + Returns + ------- + + matrix + A ``1``-by-``(1 + K.dimension())`` row vector. Its first + entry is zero, and the rest are the entries of ``e2``. + + Examples + -------- + + >>> from dunshire import * + >>> K = NonnegativeOrthant(3) + >>> L = [[1,1,1],[1,1,1],[1,1,1]] + >>> e1 = [1,1,1] + >>> e2 = [1,2,3] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG._A()) + [0.0000000 1.0000000 2.0000000 3.0000000] + + """ 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. + + .. warning:: + + It is not safe to cache any of the matrices passed to + CVXOPT, because it can clobber them. + + Returns + ------- + + matrix + A ``2*K.dimension()``-by-``1 + K.dimension()`` matrix. + + Examples + -------- + + >>> from dunshire import * + >>> K = NonnegativeOrthant(3) + >>> L = [[4,5,6],[7,8,9],[10,11,12]] + >>> e1 = [1,2,3] + >>> e2 = [1,1,1] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> 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] + [ 1.0000000 -4.0000000 -5.0000000 -6.0000000] + [ 2.0000000 -7.0000000 -8.0000000 -9.0000000] + [ 3.0000000 -10.0000000 -11.0000000 -12.0000000] + + """ I = identity(self._K.dimension()) return append_row(append_col(self._zero(), -I), append_col(self._e1, -self._L)) + def _try_solution(self, c, h, C, b, tolerance): + # Actually solve the thing and obtain a dictionary describing + # what happened. + try: + solvers.options['show_progress'] = options.VERBOSE + solvers.options['abs_tol'] = tolerance + soln_dict = solvers.conelp(c,self._G(),h,C,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 + # 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():] + + # 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) + + def solution(self): """ Solve this linear game and return a :class:`Solution`. @@ -419,13 +553,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:: @@ -439,13 +573,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 @@ -464,67 +598,36 @@ class SymmetricLinearGame: # value in the statement of the CVXOPT conelp program. c = matrix([-1, self._zero()]) - # Actually solve the thing and obtain a dictionary describing - # what happened. try: - 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 + # First try with a stricter tolerance. Who knows, it might work. + return self._try_solution(c, h, C.cvxopt_dims(), b, + tolerance = options.ABS_TOL / 10) - # 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():] - - # 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). - 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): - raise GameUnsolvableException(self, soln_dict) - - return Solution(p1_value, p1_optimal, p2_optimal) + except (PoorScalingException, GameUnsolvableException): + # Ok, that didn't work. Let's try it with the default. + return self._try_solution(c, h, C.cvxopt_dims(), b, + tolerance = options.ABS_TOL) - def _condition(self): + 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 sum of the condition numbers of ``G`` and ``A`` in the + 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 -------- @@ -534,13 +637,13 @@ class SymmetricLinearGame: >>> e1 = [1] >>> e2 = e1 >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> actual = SLG._condition() - >>> expected = 3.6180339887498953 + >>> actual = SLG.condition() + >>> expected = 1.8090169943749477 >>> abs(actual - expected) < options.ABS_TOL True """ - return condition_number(self._G()) + condition_number(self._A()) + return (condition_number(self._G()) + condition_number(self._A()))/2 def dual(self): @@ -573,7 +676,7 @@ class SymmetricLinearGame: e2 = [ 1] [ 1] [ 1], - Condition((L, K, e1, e2)) = 88.953530. + Condition((L, K, e1, e2)) = 44.476... """ # We pass ``self._L`` right back into the constructor, because