Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
from scipy.spatial import Voronoi, Delaunay
from scipy import optimize
import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import norm
from itertools import islice
from time import time
from multiprocessing import Process, Queue
from queue import Empty
N = 23
R = 4
rng = np.random.default_rng()
def normalize(v):
return v / norm(v)
def in_circle(p):
return p.dot(p) <= R**2
def perp_2d(x, axis=None):
# return np.array((-x[1], x[0]))
return np.concatenate((
-np.take(x, (1,), axis=axis),
np.take(x, (0,), axis=axis)),
axis=axis)
def squared_norm(x, axis=None):
return np.sum(np.square(x), axis=axis)
def line_circle_intersection(p1, p2, finite=True):
'''
|(1-t)p1 + t p2|^2 == R^2
|dp t + p1|^2 - R^2 == 0 (dp = p2-p1)
dp.dp t^2 + 2 dp.p1 t + p1.p1 - R^2
at^2 - 2bt + c == 0 (a = dp.dp, b = -dp.p1, c = p1.p1 - R^2)
delta = b^2-ac = (dp.p1)^2 - (dp.dp)(p1.p1) + (dp.dp)R^2
t = (b +/- sqrt(delta)) / a
'''
dp = p2 - p1
b = -p1.dot(dp)
a = dp.dot(dp)
delta = R**2 * a - np.cross(p1, p2)**2
if delta < 0:
return []
delta_sqrt = np.sqrt(delta)
pts = []
for s in (-1, +1):
t = (b + s*delta_sqrt)/a
if t >= 0 and (not finite or t <= 1):
pts.append(p1 + t*dp)
return pts
def circle_circle_intersection(p1, p2, r):
dp = (p2 - p1)/2
d_sq = np.dot(dp, dp)
if d_sq > r**2:
return None
m = (p1 + p2) / 2
perp = perp_2d(dp)
perp = perp * np.sqrt(r**2/d_sq - 1)
return (m + perp, m - perp)
def min_cover_radius(points):
points = np.reshape(points, (-1, 2))
# make sure all points are inside the circle
if not all(in_circle(p) for p in points):
return max(norm(p) for p in points)
# filter out near duplicate points that may cause an issue with
# computing the voronoi diagram
points_dedup = []
for point in points:
if all(norm(point - existing) > 1e-3*R for existing in points_dedup):
points_dedup.append(point)
points = points_dedup
# if there are too few points, we cannot compute the voronoi diagram;
# but the result would be bad anyway so we just reject those cases by
# returning a high value
if len(points) < 4:
return R
vor = Voronoi(points)
center = np.mean(vor.vertices, axis=0)
radius = min(norm(p) for p in points)
# for every edge of the voronoi diagram that crosses the circle, compute
# its intersection and compare that points distance to the centers of
# the adjacent cells
for (i1, i2), vert in zip(vor.ridge_points, vor.ridge_vertices):
p1 = vor.points[i1]
p2 = vor.points[i2]
assert len(vert) == 2
# if a vertex index is less than zero, the edge goes to infinity
v1 = vor.vertices[vert[0]] if vert[0] >= 0 else None
v2 = vor.vertices[vert[1]] if vert[1] >= 0 else None
assert v1 is not None or v2 is not None
in1 = v1 is not None and in_circle(v1)
in2 = v2 is not None and in_circle(v2)
# test if the edge crosses the circle
if not in1 or not in2:
finite = True
# make v1 finite and v2 finite or infinite
if v1 is None:
(v1, v2) = (v2, v1)
# find effective endpoint of half-infinite line
if v2 is None:
# the edge is perpendicular to the line between the two points
n = p2 - p1
t = np.array((-n[1], n[0]))
# m = (p1 + p2) / 2
# if the edge direction points towards the center, reverse it
if (v1-center).dot(t) < 0:
t = -t
# compute "effective" endpoint
v2 = v1 + t
finite = False
for v in line_circle_intersection(v1, v2, finite):
radius = max(radius, norm(v - p1), norm(v - p2))
# for every vertex of the voronoi diagram inside the circle, compare its
# distance to the centers of the adjacent cells
for p, ri in zip(vor.points, vor.point_region):
for vi in vor.regions[ri]:
if vi >= 0:
v = vor.vertices[vi]
if in_circle(v):
radius = max(radius, norm(v - p))
return radius
class plot_thread:
def __init__(self, queue, interval):
self.queue = queue
self.interval = interval
def __call__(self):
try:
self.fig, self.ax = plt.subplots()
timer = self.fig.canvas.new_timer(interval=(1000*self.interval))
timer.add_callback(self.callback)
timer.start()
plt.show()
except KeyboardInterrupt:
pass
self.terminate()
def callback(self):
while True:
try:
xr = self.queue.get_nowait()
if xr is None:
self.terminate()
return False
x, r = xr
x = np.reshape(x, (-1, 2))
self.ax.clear()
self.ax.scatter(x[:, 0], x[:, 1])
self.ax.add_patch(plt.Circle((0, 0), R, fill=False))
for i, p in enumerate(x):
self.ax.add_patch(plt.Circle(p, r, fill=False))
self.ax.annotate(str(i), p)
# plt.show(block=False)
self.ax.set_xlim((-1.25*R, 1.25*R))
self.ax.set_ylim((-1.25*R, 1.25*R))
self.ax.set_title(f'r = {r:.6f}, R = {R/r:.6f}')
self.ax.set_aspect('equal')
except Empty:
break
self.fig.canvas.draw()
return True
def terminate(self):
plt.close(self.fig)
class plotting_wrapper:
def __init__(self, fun):
self.fun = fun
self.count = 0
self.best_value = None
self.best_points = None
self.new_best = False
self.last_plot_time = None
self.plot_interval = 1 # seconds
self.plot_queue = Queue()
self.plotter = plot_thread(self.plot_queue, self.plot_interval)
self.plot_thread = Process(target=self.plotter)
self.plot_thread.start()
def __call__(self, x):
self.count += 1
r = self.fun(x)
self.store_solution(x, r)
return r
def store_solution(self, x, r):
# store best for printing
if self.best_value is None or r < self.best_value:
self.best_value = r
self.best_points = np.reshape(x, (-1, 2))
self.new_best = True
# print best periodically
now = time()
if self.last_plot_time is None or self.last_plot_time + 1 < now:
self.last_plot_time = now
# print(f'{self.count:,d} {self.best:.6f}')
if self.new_best:
self.print_best()
self.new_best = False
def print_best(self):
r = self.best_value
x = self.best_points
print(f'current best (r = {r:.6f}, R = {R/r:.6f}):')
for i, p in enumerate(x):
print(f'{i:2d}, {p[0]:8.5f}, {p[1]:8.5f}')
self.plot_queue.put((x, r))
def terminate(self):
self.plot_queue.put(None)
self.plot_thread.join()
def __del__(self):
self.terminate()
class triple_point_constr:
'''
r >= circumradius of triangle formed by circle centers
circumradius = (|p1-p2| |p2-p3| |p3-p1|) / (4 area)
area = 1/2 | x1y2 - x1y3 + x2y3 - x2y1 + x3y1 - x3y2 |
so:
(2r(x1y2 - x2y1 + x2y3 - x3y2 + x3y1 - x1y3))^2 - (x1^2+y1^2)(x2^2+y2^2)(x3^2+y3^2) >= 0
'''
def __init__(self, indexes):
self.indexes = np.array(indexes)
def lb(self):
# (m,)
return np.zeros((self.indexes.shape[0],))
def ub(self):
# (m,)
return np.full((self.indexes.shape[0],), np.inf)
def fun(self, x):
# (m,)
p = np.reshape(x[:-1], (-1, 2))
r = x[-1]
p = p[self.indexes]
q = np.roll(p, 1, axis=1)
d = p-q
area = np.sum(np.cross(p, q, axisa=2, axisb=2), axis=1)
sqnorm = squared_norm(d, axis=2)
return np.square(2*r*area) - np.prod(sqnorm, axis=1)
def jacobian(self, x):
# (m,n); d(fun)/dx
m = len(self.indexes)
n = len(x)
k = (n-1)//2
p = np.reshape(x[:-1], (k, 2))
r = x[-1]
q = p[self.indexes]
d = q - np.roll(q, -1, axis=1)
s = squared_norm(d, axis=2)[..., np.newaxis]
area = np.sum(q[:, :, 0]*np.roll(d[:, :, 1], -1, axis=1), axis=1)
perp = perp_2d(d, axis=2)
jac_elem = -8*(r**2)*np.roll(perp, -1, axis=1)*area[..., np.newaxis, np.newaxis] - \
2*np.roll(s, -1, axis=1)*(d*np.roll(s, 1, axis=1) -
np.roll(d, 1, axis=1)*s)
jac = np.zeros((m, k, 2))
i = np.repeat(np.arange(m)[..., np.newaxis], 3, axis=1)
j = self.indexes
jac[i, j] = jac_elem
jac = np.reshape(jac, (m, n-1))
return np.concatenate((jac, 8*r*np.square(area)[..., np.newaxis]), axis=1)
def hessian(self, x, v):
# (n,n)
pass
class edge_point_constr:
'''
|intersection of circles|^2 >= R^2
dp = (p2 - p1)/2
d_sq = np.dot(dp, dp)
m = (p1 + p2) / 2
perp = np.array((-dp[1], dp[0]))
intersection = m +/- perp * np.sqrt(r**2/d_sq - 1)
|intersection|^2 = m.m + 2*abs(m.perp sqrt(r**2/d_sq - 1)) + dp.dp (r**2/d_sq - 1) >= R^2
2*abs(m.perp sqrt(r**2/d_sq - 1)) >= R^2 - m.m - (r**2 - dp.dp)
(x1y2-x2y1)^2*(r**2/d_sq - 1) >= (R^2 - r**2 - p1.p2)^2
(x1y2-x2y1)^2*(r**2 - d_sq) >= (R^2 - r**2 - p1.p2)^2*d_sq
'''
def __init__(self, indexes):
self.indexes = np.array(indexes)
def lb(self):
# (m,)
return np.zeros((self.indexes.shape[0],))
def ub(self):
# (m,)
return np.full((self.indexes.shape[0],), np.inf)
def fun(self, x):
# (m,)
p = np.reshape(x[:-1], (-1, 2))
r = x[-1]
p = p[self.indexes]
dp = (p[:, 1] - p[:, 0])/2
det = np.cross(p[:, 0], p[:, 1], axisa=1, axisb=1)
dp_sq = squared_norm(dp, axis=1)
dot = np.sum(p[:, 0]*p[:, 1], axis=1)
return np.square(det)*(r**2 - dp_sq) - np.square(R**2 - r**2 - dot)*dp_sq
def jacobian(self, x):
'''
1/2 (-("cross"^2 + ("dot" + r^2 - R^2)^2) {x1 - x2, y1 - y2} -
"dsq" ("dot" + r^2 - R^2) {x2, y2} -
"cross" ("dsq" - 4 r^2) {-y2, x2})
1/2 (-("cross"^2 + ("dot" + r^2 - R^2)^2) {x2 - x1, y2 - y1} -
"dsq" ("dot" + r^2 - R^2) {x1, y1} + <===
"cross" ("dsq" - 4 r^2) {-y1, x1})
'''
# (m,n); d(fun)/dx
m = len(self.indexes)
n = len(x)
k = (n-1)//2
p = np.reshape(x[:-1], (k, 2))
j_p = np.reshape(np.eye(2*k, n), (k, 2, n))
r = x[-1]
j_r = np.zeros((n,))
j_r[-1] = 1
p = p[self.indexes]
j_p = j_p[self.indexes]
dp = (p[:, 1] - p[:, 0])/2
j_dp = (j_p[:, 1]-j_p[:, 0])/2
det = np.cross(p[:, 0], p[:, 1], axisa=1, axisb=1)
j_det = \
p[:, 1, 1, np.newaxis]*j_p[:, 0, 0] + \
p[:, 0, 0, np.newaxis]*j_p[:, 1, 1] - \
p[:, 1, 0, np.newaxis]*j_p[:, 0, 1] - \
p[:, 0, 1, np.newaxis]*j_p[:, 1, 0]
dp_sq = squared_norm(dp, axis=1)
j_sq = np.sum(2 * dp[..., np.newaxis] * j_dp, axis=1)
dot = np.sum(p[:, 0]*p[:, 1], axis=1)
j_dot = np.sum(p[:, 1, :, np.newaxis]*j_p[:, 0] +
p[:, 0, :, np.newaxis]*j_p[:, 1], axis=1)
return 2*det[..., np.newaxis]*j_det*(r**2 - dp_sq)[..., np.newaxis] + \
np.square(det)[..., np.newaxis]*(2*r[..., np.newaxis]*j_r - j_sq) - \
(2*(R**2 - r**2 - dot)[..., np.newaxis]*(-2*r[..., np.newaxis]*j_r-j_dot)*dp_sq[..., np.newaxis] +
np.square(R**2 - r**2 - dot)[..., np.newaxis]*j_sq)
def hessian(self, x, v):
# (n,n)
pass
def local_minimize(fun, x0, *args, **kwargs):
p0 = np.reshape(x0, (-1, 2))
r0 = fun(x0)
print('starting local optimization...')
triang = Delaunay(p0)
assert len(triang.coplanar) == 0
triple_points = []
for idxs in triang.simplices:
p = p0[idxs]
q = np.roll(p, 1, axis=0)
if all(norm(p-q, axis=1) <= 2*r0):
triple_points.append(idxs)
triple_points = np.array(triple_points)
edge_points = []
for i1, p1 in enumerate(p0):
for i2, p2 in islice(enumerate(p0), i1 + 1, None):
qs = circle_circle_intersection(p1, p2, r0)
if qs is None:
continue
for q in qs:
if norm(q) >= R:
edge_points.append((i1, i2))
# ax.add_patch(plt.Circle(q, 0.1, color='blue'))
constr = []
if len(triple_points) > 0:
constr.append(triple_point_constr(triple_points))
if len(edge_points) > 0:
constr.append(edge_point_constr(edge_points))
x = np.concatenate((x0.flatten(), (r0,)))
n = len(x)
j = np.zeros((n,))
j[-1] = 1
h = np.zeros((n, n))
result = optimize.minimize(lambda x: x[-1], x, method='SLSQP',
jac=lambda _: j, # hess=lambda _: h,
constraints=[{'type': 'ineq', 'fun': c.fun,
'jac': c.jacobian} for c in constr]
)
# print(result)
x = result.x[:-1]
r = result.x[-1]
r_actual = fun(x)
# print(r_actual, r, r_actual-r)
# if r < R and r_actual < R:
# assert(abs(r_actual-r) < 1e-6)
print(
f'improved from {r0:.6f} to {r:.6f} / {r_actual:.6f}: {result.message}')
result.fun = r_actual
result.x = x
return result
if __name__ == '__main__':
fun = plotting_wrapper(min_cover_radius)
try:
bounds = [(-R, +R) for _ in range(2*N)]
res = optimize.dual_annealing(fun, bounds,
maxiter=10_000, initial_temp=7000,
local_search_options={'method': local_minimize})
fun.print_best()
print('done!')
except KeyboardInterrupt:
pass
fun.terminate()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment