Created
June 28, 2021 15:33
-
-
Save robertcampion/ac5f8fdddb7af631ada53c57849f7b60 to your computer and use it in GitHub Desktop.
Code to solve https://puzzling.stackexchange.com/q/110592
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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