From: Michael Orlitzky Date: Fri, 11 Nov 2016 20:23:37 +0000 (-0500) Subject: Add a ball_radius() method for cones and use it to compute starting points. X-Git-Tag: 0.1.0~37 X-Git-Url: http://gitweb.michael.orlitzky.com/?p=dunshire.git;a=commitdiff_plain;h=2132f293d3ab198630f9fa26151eed52b21512fb Add a ball_radius() method for cones and use it to compute starting points. --- diff --git a/dunshire/cones.py b/dunshire/cones.py index 7fc8bf9..af2415b 100644 --- a/dunshire/cones.py +++ b/dunshire/cones.py @@ -77,6 +77,44 @@ class SymmetricCone: raise NotImplementedError + def ball_radius(self, point): + """ + Return the radius of a ball around ``point`` in this cone. + + Since a radius cannot be negative, the ``point`` must be + contained in this cone; otherwise, an error is raised. + + Parameters + ---------- + + point : matrix + A point contained in this cone. + + Returns + ------- + + float + A radius ``r`` such that the ball of radius ``r`` centered at + ``point`` is contained entirely within this cone. + + Raises + ------ + + NotImplementedError + Always, this method must be implemented in subclasses. + + Examples + -------- + + >>> K = SymmetricCone(5) + >>> K.ball_radius(matrix([1,1,1,1,1])) + Traceback (most recent call last): + ... + NotImplementedError + + """ + raise NotImplementedError + def dimension(self): """ @@ -201,6 +239,64 @@ class NonnegativeOrthant(SymmetricCone): return all([x > -options.ABS_TOL for x in point]) + def ball_radius(self, point): + """ + Return the radius of a ball around ``point`` in this cone. + + Since a radius cannot be negative, the ``point`` must be + contained in this cone; otherwise, an error is raised. + + The minimum distance from ``point`` to the complement of this + cone will always occur at its projection onto that set. It is + not hard to see that the projection is directly along one of the + coordinates, and so the minimum distance from ``point`` to the + boundary of this cone is the smallest coordinate of + ``point``. Therefore any radius less than that will work; we + divide it in half somewhat arbitrarily. + + Parameters + ---------- + + point : matrix + A point contained in this cone. + + Returns + ------- + + float + A radius ``r`` such that the ball of radius ``r`` centered at + ``point`` is contained entirely within this cone. + + Raises + ------ + + TypeError + If ``point`` is not a :class:`cvxopt.base.matrix`. + + TypeError + If ``point`` has the wrong dimensions. + + ValueError + if ``point`` is not contained in this cone. + + Examples + -------- + + >>> K = NonnegativeOrthant(5) + >>> K.ball_radius(matrix([1,2,3,4,5])) + 0.5 + + """ + if not isinstance(point, matrix): + raise TypeError('the given point is not a cvxopt.base.matrix') + if not point.size == (self.dimension(), 1): + raise TypeError('the given point has the wrong dimensions') + if not point in self: + raise ValueError('the given point does not lie in the cone') + + return min(list(point)) / 2.0 + + class IceCream(SymmetricCone): """ @@ -305,6 +401,80 @@ class IceCream(SymmetricCone): return norm(radius) < (height + options.ABS_TOL) + def ball_radius(self, point): + """ + Return the radius of a ball around ``point`` in this cone. + + Since a radius cannot be negative, the ``point`` must be + contained in this cone; otherwise, an error is raised. + + The minimum distance from ``point`` to the complement of this + cone will always occur at its projection onto that set. It is + not hard to see that the projection is at a "down and out" angle + of ``pi/4`` towards the outside of the cone. If one draws a + right triangle involving the ``point`` and that projection, it + becomes clear that the distance between ``point`` and its + projection is a factor of ``1/sqrt(2)`` times the "horizontal" + distance from ``point`` to boundary of this cone. For simplicity + we take ``1/2`` instead. + + Parameters + ---------- + + point : matrix + A point contained in this cone. + + Returns + ------- + + float + A radius ``r`` such that the ball of radius ``r`` centered at + ``point`` is contained entirely within this cone. + + Raises + ------ + + TypeError + If ``point`` is not a :class:`cvxopt.base.matrix`. + + TypeError + If ``point`` has the wrong dimensions. + + ValueError + if ``point`` is not contained in this cone. + + Examples + -------- + + The height of ``x`` is one (its first coordinate), and so the + radius of the circle obtained from a height-one cross section is + also one. Note that the last two coordinates of ``x`` are half + of the way to the boundary of the cone, and in the direction of + a 30-60-90 triangle. If one follows those coordinates, they hit + at ``(1, sqrt(3)/2, 1/2)`` having unit norm. Thus the + "horizontal" distance to the boundary of the cone is ``(1 - + norm(x)``, which simplifies to ``1/2``. And rather than involve + a square root, we divide by two for a final safe radius of + ``1/4``. + + >>> from math import sqrt + >>> K = IceCream(3) + >>> x = matrix([1, sqrt(3)/4.0, 1/4.0]) + >>> K.ball_radius(x) + 0.25 + + """ + if not isinstance(point, matrix): + raise TypeError('the given point is not a cvxopt.base.matrix') + if not point.size == (self.dimension(), 1): + raise TypeError('the given point has the wrong dimensions') + if not point in self: + raise ValueError('the given point does not lie in the cone') + + height = point[0] + radius = norm(point[1:]) + return (height - radius) / 2.0 + class SymmetricPSD(SymmetricCone): r""" diff --git a/dunshire/games.py b/dunshire/games.py index f9877b3..bb808cb 100644 --- a/dunshire/games.py +++ b/dunshire/games.py @@ -5,7 +5,7 @@ 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, IceCream, NonnegativeOrthant +from .cones import CartesianProduct from .errors import GameUnsolvableException, PoorScalingException from .matrices import (append_col, append_row, condition_number, identity, inner_product, norm, specnorm) @@ -820,26 +820,9 @@ class SymmetricLinearGame: :meth:`L` is satisfied. """ p = self.e2() / (norm(self.e2()) ** 2) - - # Compute the distance from p to the outside of K. - if isinstance(self.K(), NonnegativeOrthant): - # How far is it to a wall? - dist = min(list(self.e1())) - elif isinstance(self.K(), IceCream): - # How far is it to the boundary of the ball that defines - # the ice-cream cone at a given height? Now draw a - # 45-45-90 triangle and the shortest distance to the - # outside of the cone should be 1/sqrt(2) of that. - # It works in R^2, so it works everywhere, right? - # We use "2" because it's better numerically than sqrt(2). - height = self.e1()[0] - radius = norm(self.e1()[1:]) - dist = (height - radius) / 2 - else: - raise NotImplementedError - + dist = self.K().ball_radius(self.e1()) nu = - specnorm(self.L())/(dist*norm(self.e2())) - x = matrix([nu,p], (self.dimension() + 1, 1)) + x = matrix([nu, p], (self.dimension() + 1, 1)) s = - self._G()*x return {'x': x, 's': s} @@ -850,29 +833,12 @@ class SymmetricLinearGame: Return a feasible starting point for player two. """ q = self.e1() / (norm(self.e1()) ** 2) - - # Compute the distance from p to the outside of K. - if isinstance(self.K(), NonnegativeOrthant): - # How far is it to a wall? - dist = min(list(self.e2())) - elif isinstance(self.K(), IceCream): - # How far is it to the boundary of the ball that defines - # the ice-cream cone at a given height? Now draw a - # 45-45-90 triangle and the shortest distance to the - # outside of the cone should be 1/sqrt(2) of that. - # It works in R^2, so it works everywhere, right? - # We use "2" because it's better numerically than sqrt(2). - height = self.e2()[0] - radius = norm(self.e2()[1:]) - dist = (height - radius) / 2 - else: - raise NotImplementedError - + dist = self.K().ball_radius(self.e2()) omega = specnorm(self.L())/(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)) + z = matrix([z1, z2], (self.dimension()*2, 1)) return {'y': y, 'z': z}