From 2a59d8c4a2d7723a52de15e3d1f8ffce76bb948a Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Thu, 3 Nov 2016 15:56:22 -0400 Subject: [PATCH] Add more docs for a few private methods and implement the do-over. Some examples and return type documentation were added for the _zero(), _A(), and _G() methods of SymmetricLinearGame. The _zero() method was also uncached, since attempting to cache _A() and _G() revealed that CVXOPT clobbers them. The zero vector appears safe, but let's play it safe. This also implements a two-stage retry for the solution of a game. Our default tolerance is now 1e-6, but when we're trying to solve a game, we try it with 1e7 first. If we can solve it with the stricter tolerance, then great. If not, we use the default. --- dunshire/games.py | 210 +++++++++++++++++++++++++++++++++------------- 1 file changed, 150 insertions(+), 60 deletions(-) diff --git a/dunshire/games.py b/dunshire/games.py index c2ca68e..6480d7d 100644 --- a/dunshire/games.py +++ b/dunshire/games.py @@ -322,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): @@ -352,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): @@ -365,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`. @@ -462,61 +598,15 @@ 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: - 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 - # 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) > 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) + # 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) + + 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): -- 2.44.2