In [2]:
import numpy as np
import cdd
import itertools as it
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.patches import Polygon, Ellipse
from matplotlib.collections import PatchCollection
In [3]:
def h2v(A, b):
"""Converts the halfspace representation Ax <= b of a polyhedron P
to a vertex-representation
:param A: ndarray of shape (m, d)
:param b: ndarray of shape (d)
:returns: ndarray of shape (n, d)
"""
assert len(A) == len(b)
mat = cdd.Matrix([[b_i] + list(-A_i) for (b_i, A_i) in zip(b, A)],
number_type='float')
mat.rep_type = cdd.RepType.INEQUALITY
vertices = cdd.Polyhedron(mat).get_generators()
result = []
for v in vertices:
assert v[0] == 1.0
result.append(v[1:])
return np.array(result)
In [11]:
def proj(ket):
return ket[:, None] * ket.conj()[None, :]
In [12]:
def plot_poly(A, b, ax=None):
ax = ax if ax is not None else pl.gca()
vertices = h2v(A, b)
poly = Polygon(vertices, closed=True)
p = PatchCollection([poly])
ax.add_collection(p)
def plot_hyperplane(c, x, ax=None):
ax = ax if ax is not None else pl.gca()
t = np.linspace(-1, 1, 20)
xs = x[0] + t * c[1]
ys = x[1] - t * c[0]
ax.plot(xs, ys)
ax.arrow(x[0], x[1], c[0], c[1], fc='k', ec='k')
def plot_ellipse(x, A, ax=None):
ax = ax if ax is not None else pl.gca()
vals, vecs = np.linalg.eigh(A)
phi = np.arccos(vecs[0][0])
ell = Ellipse(x, width=vals[0], height=vals[1], angle=phi / np.pi * 180,
fill=False, linestyle='solid', linewidth=1, edgecolor='k')
ax.add_artist(ell)
In [13]:
def normalize_hrep(A, b):
norm = np.linalg.norm(A, axis=-1)
return A / norm[:, None], b / norm
In [14]:
SQUARE2D = (np.array([[1,0], [-1,0], [0,1], [0,-1]], dtype=np.float64),
np.ones(4, dtype=np.float64))
SQUARE2D = normalize_hrep(*SQUARE2D)
In [15]:
def contained(x, poly, delta=1e-6):
A, b = poly
violated = A.dot(x) > b + delta
if not np.any(violated):
return True
else:
return A[violated][0]
In [16]:
POINT = np.array([ 1.26403501, 1.26403501])
plot_poly(*SQUARE2D)
pl.scatter([POINT[0]], [POINT[1]])
res = contained(POINT, SQUARE2D)
if not (res is True):
plot_hyperplane(res, POINT)
In [17]:
def ellipsoid_iter(c, oracle, ball):
n = len(c)
x, R = ball
A = R**2 * np.eye(n)
while True:
sep = oracle(x)
yield x, A, sep
a = c if sep is True else -sep
x += 1. / (n+1) * A.dot(a) / np.sqrt(a.dot(A.dot(a)))
A = n**2 / (n**2 - 1) * (A - 2. / (n+1) * A.dot(proj(a).dot(A)) / a.dot(A.dot(a)))
In [18]:
SCORE = np.array([-2., 1.])
BALL = (np.array([0., 0.]), 2.)
solution = ellipsoid_iter(SCORE, lambda x: contained(x, SQUARE2D), BALL)
for x, A, c in it.islice(solution, 0, 20, 2):
plot_poly(*SQUARE2D)
pl.scatter([x[0]], [x[1]], c='r')
plot_ellipse(x, A)
if not (c is True):
plot_hyperplane(c, x)
pl.show()
print(x, np.linalg.det(A))
In [ ]:
Ellipsis()