Code to solve https://puzzling.stackexchange.com/q/110592

disks.py

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()